Tree Building Algorithm that treats gaps as deletions

Question

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.

Answer 1

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.