issues converting parallel_singer output to tskit trees

[h-pawar]

hi, many thanks for developing this cool method!

I ran parallel_singer per chromosome, & the output files look like the ARG generation was successful:

 aw_n82_38_0_1_start_nodes_0.txt
 aw_n82_38_0_1_start_branches_0.txt
 aw_n82_38_0_1_start_recombs_0.txt
 aw_n82_38_0_1_start_muts_0.txt
 aw_n82_38_0_1_nodes_0.txt
 aw_n82_38_0_1_branches_0.txt
 aw_n82_38_0_1_recombs_0.txt
 aw_n82_38_0_1_muts_0.txt

But I am having some issues converting the output format to tskit trees, I tried using convert_to_tskit & convert_long_ARG.py

1. using convert_to_tskit, gives an error I guess because it doesn't expect per window outputs (prefix_0_1_nodes_0.txt instead of prefix_nodes_0.txt)

ID=38
input="aw_n82_${ID}"
 outp="trees/aw_n82_${ID}"
convert_to_tskit \
>  -input ${input} \
>  -output ${outp} \
>  -start 0 \
>  -end 99
Traceback (most recent call last):
  File "convert_to_tskit", line 98, in <module>
    main()
    ~~~~^^
  File "convert_to_tskit", line 95, in main
    write_trees(args.input, args.output, args.start, args.end, args.step)
    ~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "convert_to_tskit", line 64, in write_trees
    ts = read_ARG(node_file, branch_file, mutation_file)
  File "convert_to_tskit", line 52, in read_ARG
    tables = read_ts(node_file, branch_file)
  File "convert_to_tskit", line 10, in read_ts
    node_time = np.loadtxt(node_file)
  File "/home/xcj768/miniforge3/envs/tskit-env/lib/python3.14/site-packages/numpy/lib/_npyio_impl.py", line 1384, in loadtxt
    arr = _read(fname, dtype=dtype, comment=comment, delimiter=delimiter,
                converters=converters, skiplines=skiprows, usecols=usecols,
                unpack=unpack, ndmin=ndmin, encoding=encoding,
                max_rows=max_rows, quote=quotechar)
  File "/home/xcj768/miniforge3/envs/tskit-env/lib/python3.14/site-packages/numpy/lib/_npyio_impl.py", line 1011, in _read
    fh = np.lib._datasource.open(fname, 'rt', encoding=encoding)
  File "/home/xcj768/miniforge3/envs/tskit-env/lib/python3.14/site-packages/numpy/lib/_datasource.py", line 192, in open
    return ds.open(path, mode, encoding=encoding, newline=newline)
           ~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/xcj768/miniforge3/envs/tskit-env/lib/python3.14/site-packages/numpy/lib/_datasource.py", line 529, in open
    raise FileNotFoundError(f"{path} not found.")
FileNotFoundError: aw_n82_38_nodes_0.txt not found.

2. using convert_long_ARG.py following Understanding parallel singer errror #21, gives error - mutation's parent is not consistent with the topology

python convert_long_ARG.py \
> -vcf ../../input/aw_n82_38 \
> -output aw_n82_38 \
> -iteration 0
Processing segment 0
Processing segment 1
Processing segment 2
Processing segment 3
Processing segment 4
Processing segment 5
Processing segment 6
Processing segment 7
Processing segment 8
Processing segment 9
Processing segment 10
Processing segment 11
Processing segment 12
Processing segment 13
Processing segment 14
Processing segment 15
Processing segment 16
Processing segment 17
Processing segment 18
Processing segment 19
Processing segment 20
Processing segment 21
Processing segment 22
Processing segment 23
Traceback (most recent call last):
  File "convert_long_ARG.py", line 118, in <module>
    main()
    ~~~~^^
  File "convert_long_ARG.py", line 114, in main
    ts = read_long_ARG(node_files, branch_files, mutation_files, block_coordinates)
  File "convert_long_ARG.py", line 74, in read_long_ARG
    ts = tables.tree_sequence()
  File "/home/xcj768/miniforge3/envs/tskit-env/lib/python3.14/site-packages/tskit/tables.py", line 3528, in tree_sequence
    return tskit.TreeSequence.load_tables(self)
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^
  File "/home/xcj768/miniforge3/envs/tskit-env/lib/python3.14/site-packages/tskit/trees.py", line 4295, in load_tables
    ts.load_tables(tables._ll_tables, build_indexes=build_indexes)
    ~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
_tskit.LibraryError: A mutation's parent is not consistent with the topology of the tree. Use compute_mutation_parents to set the parents correctly.(TSK_ERR_BAD_MUTATION_PARENT)

Which of the two conversion methods would you recommend for the parallel_singer output convert_to_tskit or convert_long_ARG.py? And your input on how to fix the errors is much appreciated!

[nspope]

This is an issue with tskit v1.0 compatibility -- it requires mutation parents to be assigned (i.e. the sequence of inheritance for multiple mutations at the same site), and this is not currently done in either of the conversion scripts afaict. A quick workaround would be to downgrade to the most recent version of tskit below 1.0, but obviously this isn't ideal as then any downstream analyses will be tied to an older tskit version.

Here's a version of convert_to_tskit.py that does the conversion correctly, but that doesn't help you as you'll need convert_to_long_arg.py. When I have a minute I'll put up code for that (likely in a PR).

"""
Modified version of `convert_to_tskit.py` from the SINGER github repository
(github.com/popgenmethods/SINGER) that computes mutation parents correctly so
as to load without error in tskit>=1.0
"""

import sys
import os
import argparse
import numpy as np
import tskit


def read_ts(node_file, edge_file):
    node_time = np.loadtxt(node_file)
    edge_span = np.loadtxt(edge_file)
    edge_span = edge_span[edge_span[:, 2] >= 0, :]
    length = max(edge_span[:, 1])
    tables = tskit.TableCollection(sequence_length=length)
    node_table = tables.nodes
    edge_table = tables.edges
    prev_time = -1
    for t in node_time:
        if (t == 0):
            node_table.add_row(flags=tskit.NODE_IS_SAMPLE)
        else:
            t = max(prev_time + 1e-4, t)
            node_table.add_row(time = t)
            prev_time = t
    parent_indices = np.array(edge_span[:, 2], dtype = np.int32)
    child_indices = np.array(edge_span[:, 3], dtype = np.int32)
    edge_table.set_columns(left = edge_span[:, 0], right = edge_span[:, 1], parent = parent_indices, child = child_indices)
    return tables


def read_mutation(tables, mutation_file):
    mutations = np.loadtxt(mutation_file)
    n = mutations.shape[0]
    mut_pos = 0
    for i in range(n):
        if mutations[i, 0] != mut_pos:
            tables.sites.add_row(position=mutations[i, 0], ancestral_state='0')
            mut_pos = mutations[i, 0]
        site_id = tables.sites.num_rows - 1
        tables.mutations.add_row(site=site_id, node=int(mutations[i, 1]), derived_state=str(int(mutations[i, 3])))
    # added: rebuild mutations table in reverse time order at each position
    mut_time = tables.nodes.time[tables.mutations.node]
    mut_coord = tables.sites.position[tables.mutations.site]
    mut_order = np.lexsort((-mut_time, mut_coord))
    mut_state = tskit.unpack_strings(
        tables.mutations.derived_state, 
        tables.mutations.derived_state_offset,
    )
    mut_state, mut_state_offset = tskit.pack_strings(np.array(mut_state)[mut_order])
    tables.mutations.set_columns(
        site=tables.mutations.site[mut_order],
        node=tables.mutations.node[mut_order],
        time=np.repeat(tskit.UNKNOWN_TIME, tables.mutations.num_rows),
        derived_state=mut_state,
        derived_state_offset=mut_state_offset,
    )


def read_ARG(node_file, branch_file, mutation_file):
    tables = read_ts(node_file, branch_file)
    read_mutation(tables, mutation_file)
    # added: calculate mutation parents and times
    tables.sort()
    tables.build_index()
    tables.compute_mutation_parents()
    tables.compute_mutation_times()
    ts = tables.tree_sequence()
    return ts


def write_fast_trees(input_prefix, output_prefix, start, end, step):
    for i in range(start, end, step):
        trees_file = f"{output_prefix}_{i}.trees"
        node_file = f"{input_prefix}_fast_nodes_{i}.txt"
        branch_file = f"{input_prefix}_fast_branches_{i}.txt"
        mutation_file = f"{input_prefix}_fast_muts_{i}.txt"
        ts = read_ARG(node_file, branch_file, mutation_file)
        ts.dump(trees_file)


def write_trees(input_prefix, output_prefix, start, end, step):
    for i in range(start, end, step):
        trees_file = f"{output_prefix}_{i}.trees"
        node_file = f"{input_prefix}_nodes_{i}.txt"
        branch_file = f"{input_prefix}_branches_{i}.txt"
        mutation_file = f"{input_prefix}_muts_{i}.txt"
        ts = read_ARG(node_file, branch_file, mutation_file)
        ts.dump(trees_file) 


if __name__ == '__main__':
	
	if "snakemake" in globals():
		writer = write_fast_trees if snakemake.params.fast else write_trees
		writer(
			snakemake.params.input, 
			snakemake.params.output, 
			snakemake.params.start, 
			snakemake.params.end, 
			snakemake.params.step,
		)
	else:
		parser = argparse.ArgumentParser(description='Convert to tskit format')
		parser.add_argument('-input', type=str, required=True, help='Prefix of ARG files.')
		parser.add_argument('-output', type=str, required=True, help='Prefix of output files.')
		parser.add_argument('-start', type=int, required=True, help='Start index of the sample.')
		parser.add_argument('-end', type=int, required=True, help='End index of the sample.')
		parser.add_argument('-step', type=int, default=1, help='Step size of subsampling. Default: 1.') 
		parser.add_argument('-fast', action='store_true', help='Use this flag for fast-SINGER samples.')
		if len(sys.argv) == 1:
			parser.print_help(sys.stderr)
			sys.exit(1)
		args = parser.parse_args()
		if args.fast:
			write_fast_trees(args.input, args.output, args.start, args.end, args.step)
		else:
			write_trees(args.input, args.output, args.start, args.end, args.step)

[hannesbecher]

Not pressure @nspope - that replacement for convert_to_long_arg.py. would come in handy for me too...

[hannesbecher]

A quick fix, using tskit 1.0, is to import the tables and process into a TreeSequence:

tabs = tskit.TableCollection.load("mySingerTeesFile.trees")
tabs.sort()
tabs.build_index()
tabs.compute_mutation_parents()
tabs.compute_mutation_times()
ts=tabs.tree_sequence()
ts

[pratikkatte]

This worked for me

def read_ARG(node_file, branch_file, mutation_file):
    tables = read_ts(node_file, branch_file)
    read_mutation(tables, mutation_file)
    
    # Add proper indexing and parent mutation computation before creating tree sequence
    tables.sort()  
    tables.build_index()
    tables.compute_mutation_parents()
    
    ts = tables.tree_sequence()
    return ts

In convert_long_ARG.py

def read_long_ARG(node_files, branch_files, mutation_files, block_coordinates):
    # ... (existing code) ...
    
    # Add proper indexing and parent mutation computation before creating tree sequence
    tables.sort()
    tables.build_index()
    tables.compute_mutation_parents()
    
    ts = tables.tree_sequence()
    return ts