Phylogenetic Separation

Question

I am having a bothersome issue with my current approach in generating a phylogenetic tree of yeast samples, where samples from two groups (in-house and public NCBI data) form separate clades, despite representing the same biological sample.

I have about 50 strains going into this tree that come from two different groups:

1. sequenced data from our in-house strains
2. sequenced data from NCBI (of samples where we directly obtained our in-house strains)

Simply put, sample_1 from group_1 is directly sourced from sample_1 from group_2 so in theory these two samples should sit right next to each other in a phylogenetic tree. However, my final tree shows two different clades (each with their own branching) which consist of group_1 for one clade and group_2 for the other - so they are completely separate in the tree.

I tried to troubleshoot this by:

• making sure all the samples have the same read lengths (and trimming adapters)
• downsample each sample to the same coverage where I can
• normalizing the multi-sample final VCF file (bcftools norm)
• filtering the VCF file with basic filters (bcftools filter -i 'QUAL > 30 && FORMAT/DP > 10 && MQM > 40 && (QUAL / FORMAT/DP) > 2.0')
• filtering the allele frequency (AF) to be > 0.25

Doing this helped put a handful of samples next to their partner from the other group in a new tree I've made but still not all of them are together the way they should be. I assume the biggest issue comes from the difference in sequencing runs / how I calculate the distance so how would you troubleshoot this error further? Has anyone experienced an issue like this?

This post is already so long but in case it helps, this is my current approach:

• pre-process the sequences as mentioned above (standardize read length and coverage, trim adapters)
• align sequences to reference genome using BWA MEM
• mark duplicates with picard, add read groups with picard, mark duplicates again
• filter the final alignment by quality (samtools view -q 50)
• call variants using freebayes and a list of every bam file to call simultaneously
• use gatks VariantsToTable on the vcf file to get the GT info
• use a custom R script to calculate pairwise distances
• use a custom R script to plot the tree

Answer 1

I have some suggestions. I'm assuming of course that you are using the same reference genome across samples.

1. Go find the positions that lead to this phylogenetic placement. You should be able to do this by crosstabulating genotypes in your VCF file by the artefactual groups. Do they have anything in common? For example, if you have a large number of mitochondrial variants but one of these groups of samples has had mitochondrial reads filtered out somehow, you might see this kind of thing.
2. Directly inspect a subset of alignments from the two groups using IGV. Are there any weird artifacts? Visualizing coverage, do you see anything systematic?
3. Are these 50 NCBI genomes from a single study? Can you just go look at methods from a paper and see if they did anything weird?