# Tree Building Algorithm that treats gaps as deletions

I'm part of a nanopore sequencing experiment that will sequence several generations of viruses. The intent is to perform directed evolution by putting selective pressure on these viruses and tracking the various mutations that occur. See here for the paper this methodology is based on.

I've been working on a software package that can compile consensus sequences from longread nanopore sequencing (see here). However, there's a wrinkle with generating phylogenetic trees from this type of data. It appears that most tree building software ignores indels. This makes sense, as gaps could easily indicate a lack of coverage from illumina sequencing methods or similar issues rather than betray the existence of an insertion/deletion. However, the nature of our sequencing method is such that a gaps really should be evaluated as mutational changes that we want to keep track of in a phylogenetic tree.

We've explored using NextStrain to generate our trees and visuals, but the tree building algorithms they employ too suffer from this gap assumption (see the issue I raised here). These include IQTree, RAxML and PhyML. I have not had much luck finding an algorithm or software that treats gaps as indels when building a phylogenetic tree. Does anyone know of such an algorithm or package I could employ? I figure I could take a stab at writing my own algorithm, but consider that a Plan B resort. No sense re-inventing the wheel and all that.

• Various distance measures can do this given an MSA. See for example academic.oup.com/bioinformatics/article/28/4/495/212883. As long as you have a transition matrix that deals with gaps it should work fine. You can build trees with e.g. neighbor joining from distance matrices. References on ML methods for dealing with indels can be found in this paper. Apr 16 at 0:41
• What the OP is saying is correct wft mainstream phylogenetics: the reason is they model point mutations - not alignment reliability. There are numerous methods (> 2 contrasting methods) to do what the OP seeks, its quite easy. There's no need to write their own algorithm and would be discourage (account for lots of parameters). Generic NJ is okish (its pairwise), its not quite what the OP wants. The virus cited in the Cell paper is Sindbis virus (alphavirus) ... very interesting. I think I'm probably going to sit out on this question, but I definitely look forward to seeing the paper. Upvoted.
– M__
Apr 16 at 3:26
• Speaking personally it will take me a good half day to understand the VEGAS system (Cell paper), which - whilst I get the gist: using viruses as replication vectors to accelerate adaptation, in particularly I think it leverages copy number thresholds of capsid proteins; I don't understand the nuance of the wet-lab system.
– M__
Apr 18 at 3:04

I need to do a more thorough testing to determine just how effective this answer may be, but I think I've put together a workflow that is close enough to an appropriate answer that I'd post it here.

In my research I came across gappy, (GitHub here), a piece of software that identifies potentially relevant gaps in a phylogenetic analysis and returns a PHYLIP-formatted file that can be run through most tree-building programs. Gappy accounts for indels, but not substitutions, so the output would need to be combined with another analysis to generate the tree I want.

Fortunately, IQTree allows for a partitioned analysis of mixed data (see here), so by feeding the original alignment AND the gappy output through IQTree you should be able to generate a tree that accounts for both substitutions and indels.

Again, more thorough testing is required, but I tested this possibility using an MSA example taken from this website. Saving the FASTA-formatted example as alignment.fasta, and running that file through gappy2, you get an output file called GAPPY2_alignment.fasta_2-inf.phy. Thus if you create a NEXUS-formatted file like the following:

#nexus
begin sets;
charset part1 = alignment.fasta: *;
charset part2 = GAPPY2_alignment.fasta_2-inf.phy: *;
end;


and run it through IQTree using the command iqtree -p example.nex, you get the following tree:


+---ENA|CAA23748|CAA23748.1
|
+-------------ENA|CAA24095|CAA24095.1
|
|              +-------------------------------------------ENA|BAA20512|BAA20512.1
+--------------|
+-----------------------ENA|CAA28435|CAA28435.1


Now, disclaimer: this is the exact same tree you get when you run alignment.fasta through IQTree, hence why more thorough testing is required. The sample size for this example is quite small, and it's therefore not surprising that accounting for indels does not change the outcome in this scenario. However, I think this serves as a solid baseline workflow, where future updates will be more tweaks than complete overhauls, if that makes sense.

Let me know if you see any glaring flaws in my logic, but I think I'm going to close this question. I may provide edits to this answer in the future if significant tweaks/overhauls wind up being done.