Creating a tree sequence object within a C API integrated program

[mawassw]

Hi everyone,

My question pertains to creating a valid tree sequence object from a table collection object within a C-based forward in time evolution simulation program intermittently, before the simplification process, to analyze the trees while they contain all genealogical information.

The issue I am facing is in the following sequence of pseudocode:

table collection object → sort table collection object using tsk_table_collection_sort() → rebuild index of table collection using tsk_table_collection_build_index() → initialize tree sequence object with the table collection object.

Is this sequence of pseudocode correct? I understand that I have to sort the table collection object to get a valid tree sequence but sorting drops the index of the table collection object and I can't create a tree sequence object with an unindexed table collection object. Hence, I am building the index and then initializing the tree sequence object.

NOTE: These steps are working without any error but the tree sequence object I obtain has 0 num_samples which does not make sense. Is the rebuilding of the index possibly leading to this issue?

Thanks for the help! I can provide more information if needed.

[benjeffery]

On first glance this seems right.

Are any of your nodes marked as sample nodes? (i.e. they have the TSK_NODE_IS_SAMPLE in their flags). It's perfectly valid to have a tree sequence with zero samples as they are just "special" nodes. See https://tskit.dev/tskit/docs/stable/data-model.html#node-table

[benjeffery]

I thought the samples should be flagged automatically by simplify though... Can you show us a snippet of your code where you're doing this?

You are correct, the default is True.

[mawassw]

I'll share some code snippets to help clarify things (sorry for the long message).

Re simplification
we use the default way of simplifying
returnvaluefortskfunctions = tsk_table_collection_simplify(&treesequencetablecollection, wholepopulationnodesarray, (2*maxPopSize), 0, NULL);

after this simplification, we reset the indexes of the node array as such:
for (k = 0; k < (2*maxPopSize); k++) {
wholepopulationnodesarray[k] = k;
}

Re flagging nodes
First, our model is based on non-WF model. When we add a row to the node table when a new individual is born, we don't flag that they are part of the sample. We add them to the wholepopulationnodesarray which we use to inform the simplification function. The reason for this, is that initially our goal was to assess the final sample of our runs. But now, since I am intermittently (during a single run) creating a tree sequence from a table collection as such every X generations (followed by simplification after the iterative tree analysis is done):

//sort the table
returnvaluefortskfunctions = tsk_table_collection_sort(&treesequencetablecollection, NULL, 0);
check_tsk_error(returnvaluefortskfunctions);
//reindex the table
returnvaluefortskfunctions = tsk_table_collection_build_index(&treesequencetablecollection, 0);
check_tsk_error(returnvaluefortskfunctions);
//create tree sequence object
tsk_treeseq_t treesequence;
returnvaluefortskfunctions = tsk_treeseq_init(&treesequence, &treesequencetablecollection, 0);
check_tsk_error(returnvaluefortskfunctions);

The issue I was getting was that the created tree sequence I was obtaining had 0 samples. So when I am iteratively analyzing its trees, and use the tsk_treeseq_get_samples(tree->tree_sequence); my analysis broke down. The solution I am testing right now is just the use wholepopulationnodesarray as the sample node to find things like the MRCA of the sample and find fixed mutations in the sample, as well as as attempt to get the time of fixation with our own algorithm.

I can provide more details if necessary.

[jeromekelleher]

But now, since I am intermittently (during a single run) creating a tree sequence from a table collection as such every X generations (followed by simplification after the iterative tree analysis is done):

Ah, you're creating the ts before simplification. That all makes sense then - you just need to set the TSK_NODE_IS_SAMPLE flag for the nodes in the current population before you call tsk_treeseq_init and it should all work fine.

[mawassw]

Thanks Jerome, I actually came to this conclusion recently. Having you reach the same thing is reassuring. Just for clarity on how to go about this: I just use tsk_node_table_update_row function to update the is_sample flag of the nodes that are in the population, correct?
and for the simplification process, I don't need to reset them back since the tsk_table_collection_simplify function unsets it when we provide a sample to the function, correct?

[jeromekelleher]

You can actually update the flags in the columns directly if you prefer - it's a bit simpler and more efficient. Something like

for (j = 0; j < population_size; j++) {
     tables.nodes->flags[population[j]] = TSK_NODE_IS_SAMPLE;
}

But update_row is equivalent.

Simplify will reset the sample flags, correct.