Approach for alignment-free phylogenetic tree generation using kmers

Question

I have a general question about building phylogenetic trees without alignment:

I want to build a tree using an alignment-free approach by counting kmers. I have found a lot of information in the literature about this type of practice but I'm having trouble finding approaches that use Illumina short read sequences (which is what I have).

In my first attempt, I was able to generate a dendrogram by counting kmers and clustering them, but it did not group the sequences I had in the “correct/expected” way (these are known strains that are part of specific subclasses related to S. cerevisiae - in theory, this kmer clustering approach should’ve been able to group the strains based on their subclass pretty easily).

There are 21 samples in total with 5 specific subclasses. Does anyone know a way to either cluster or simply calculate the distance between kmers to generate a tree of 20+ Illumina short-read samples?

EDIT: The main reason for doing this in this way is to eliminate any bias from the reference sequence and maintain all the alleles (most of the samples are tetraploid) so I really want to avoid collapsing the data with any type of preprocessing alignment or assembly.

With that being said I would want to use the whole-genome short-read sequences of these 21 samples to count the kmers across the genome for phylogenetic analysis. If it gets too complicated, I would be okay with assembling the reads first before kmer counting.

I used Kitsune to determine the optimal kmer size for the samples which comes out to 15mer. I was able to count 15mers with all the samples using three methods: pure raw reads, assembling the reads with SPAdes, and downsampling the raw reads to 5x coverage. Using Jellyfish and custom scripts I wrote, which didn't end up being too taxing on my computer/server.

Currently, I have 21 '.tsv' files that list the kmer and the kmer frequency for each individual sample and I have one giant matrix that combines that data together so I have the kmer and the kmer count for each sample in one file:

	sample_1	sample_2	sample_3	...	sample_21
AAAAA	123432	345432	56765	...	55677
AAAAG	1023	0	350655	...	1
AAAAT	462	741	2136	...	121
AAAAC	342	123	3453	...	12345
------------	----------	----------	---------------

I have looked briefly at kWIP. Ideally, I would use their distance calculation function to generate a .dist file for my giant matrix - I'll look into what inputs can be used for that function.

Answer 1

A very simple solution is to use a correlation r or r². This establishes a pairwise matrix which is then resolved under neighbour-joining (don't use UPGMA or similar). MegaX will import a pairwise matrix and produce a tree in neighbour-joining. You could even produce a bootstrapped solution (that would require custom code) - bit tricky.

I definitely recommend neighbour-joining for the tree

MegaX is a good implementation of neighbour-joining (it was originally Nei's idea and that was his package). The reason for this are complex and happy to address this in a separate question.

Algorithm Lets say r or r² with Sample 1 versus Sample 2 ... the correlation value goes in one half of a triangular matrix with zeros running down the centre of the matrix. You then proceed for every pairwise comparison until you get a triangle with Sample 1, Sample 2, Sample 3, etc ... along the top and Sample 1, Sample 2, Sample 3, etc ... along the side. That will import to MegaX (its a long while since I've used it BTW).

Distance Methods If I was performing this I wouldn't use r or r²: I would use either Jaccard, Euclidean or Cosine (very trendy right now) distances. Cosine is trendy because taps into tensors (Tensorflow and all that). Scipy will almost certainly do the distance methods listed here for your data and output a nice triangular matrix.

---

Personally if I was coding this I would absolutely perform a bootstrap, but this is a lot of coding.

Issues The issue with distance metrics without bootstrapping is that they give different solutions and there's no way identify which tree is the best tree. Bootstrapping helps a lot here and values > 80% are considered robust topologies. I've never personally known two distance metrics produce an incongruence for the same data under bootstrapping. The distance matrix issue was why everyone moved (and rigidly stuck with) formal models via maximum likelihood or Bayes - but these require alignments.

Summary Think you can go alignment free? Well yeah you can, but remember there's no such thing as a free lunch ;-)