Missing nodes cause error when preparing local %ID using GetGlobalAndLocalPercentID.perl

[XieHongX]

Dear Dr. Hiller,

Thank you so much for providing all the resources. I am sorry I still get some issues when applying the provided script.
I am currently preparing the percent identity files for the ForwardGenomics input. I run the command as:

GetGlobalAndLocalPercentID.perl alignment_for_an_element.BDB -treeFile full_species_tree -allowedAncestralNodes my_nodes -requireNoOutgroup -local -global

I have successfully got the global %ID result file, but failed in getting the local %ID result. The error message is like:

ERROR: do not find node XXX in the alignment.

after checking the input and output, I understand that the error was caused because some species are missing in the alignment, which is very normal in my dataset. Those missing species are assigned "NA" values in the global %ID result, but cause errors in local %ID calculation. I wonder if I should simply trim the input tree, so only species exist in the alignment are included (and probably solve the error)? But in this case each element will probably have different number of nodes in the local %ID result file. I don't know if this will cause problems in the formal ForwardGenomics analysis.

Alternatively, maybe I can go back to previous prank alignment step, and try to fill in all missing species in the alignment by inference (not sure if this makes sense and how)?

Thank you for any input.

Best wishes,
Hongxin

[XieHongX]

Self update.

Found that the real issue is because some species is missing in the maf file (without e line), Maf2SpanningSeq_PRANK.perl do not add these empty species, and the output alignment by default rename the ancestral nodes using existing species, which cause conflict in the "GetGlobalAndLocalPercentID.perl" step when I provide the original full species tree. Now the question is, do I try to add all missing species to the maf file? Or I go ahead with current alignment so use different subtree names for each element. Forwardgenomics seem to work on all elements at the same time, but if the naming of the subtree of each element is different from the full tree, running all elements in parallel will not work. Is it feasible to run Forwardgenomics for the elements one by one manually (so I provide different tree each time). Alternatively, I can also modify Maf2SpanningSeq_PRANK.perl to not rename internal node names. Still figuring out the best way to solve the problem

[XieHongX]

Update.

It turns out that the missing node issue was mainly cause by not running mafaddIrow before mafextract. mafaddIrow will add both I and e rows to the maf file, so Maf2SpanningSeq_PRANK.perl can utilise this information to reconstruct the missing element (but still show aligned chains upstream and downstream)

[MichaelHiller]

Hi Hongxin,

sorry, I am just back from travel.
Yes, running mafAddIRows is essential to get the e and i lines.
Does that solve the problem with the missing XXX species?

Michael

[XieHongX]

Hi Michael,
Thanks for the response. Running mafAddIRows does not fully resolve the problem because there exist some species that are fully lost (most likely due to assembly gap). To solve the problem, I used a different subtree for each element to get the percent identity file, and then run forwardGenomics.R for each element separately (only put one element in the "ID.list" file to do this trick). It seems to be working. Does running forwardGenomics.R for the element one by one have the same result as when running them in batch?

Thank
Hongxin

[MichaelHiller]

In principle, mafAddIRows should add an e line also for the species where nothing aligns. But maybe there are cases where no chain spans the locus and then nothing is added for this species. Yes, in these cases, pls prune the tree, as one does not know whether the species has a real deletion or an assembly gap.

FG will process each element individually, so running several individually or a whole batch should be the same. You would just have to do FDR correction over all.

[XieHongX]

Hi I have successfully run forward genomics for each element one by one. However I encountered another issue. It's sequential so I put it here instead of opening another
issue.

I run forward genomics for both global and branch model, but even when I set "--thresholdConserved=0", I got this warning:

"Large indel occured on a conserved internal branch:
No branches data for that element"

and the branch parameters are "NA" in the result file.

I guess this happened because of the "0" values in the local percent identity file (which indicate complete deletions, and was added by mafAddIRows). I was wondering how to solve this issue. Is there other parameters I can set to skip this issue or I need to trim the "0" branches from the input file?

[MichaelHiller]

If I remember correctly, the following may occur. You have a branch where say 70% of your element is deleted, but the remaining 30% could be well conserved in downstream branches. The downstream branches would have high %id values, but we know that they all share the big del, meaning not much of this element is left in the entire lineage.

This can throw off the branch method, which considers all branches.

Probably best to exclude this element.

You can run it with -verbose. It should then output more info on this element