Combine or merge haplotypes

[baileyjc]

Hello,

I have tried to investigate my question of whether or not it is possible to combine or merge haplotypes based on the number of mutations between them and from my current search criteria I have not found any helpful information either for or against this idea.

I am working with low-coverage whole genome sequence data and as a result, the mitochondrial genome is recovered with significant depth. I used GATK HaplotypeCaller to create a VCF file of mitochondrial SNPs only. I extracted the genotypes in R and plotted the network using pegas, but, given it is the whole mitochondrial genome, I have only 2 samples of the 150 I am working with that have the same haplotype. I have ~400 SNPs in total. I am therefore wondering if it is possible to merge or combine haplotypes that only differ by a certain number of mutations i.e. ≤ 3. This way my haplotype network would be more informative regarding the possible patterns present and would highlight the relationships that have more mutations between them given I am working with several island populations. As of now the haplotype network is uninformative regarding any relationships between or within populations. I tried to see if I could do it manually but keep running into roadblocks due to my intermediate coding abilities/understanding.

I apologize if my issue has already been addressed elsewhere or if I am not theoretically or empirically sound in my understanding. If need be, I will remove my issue.

Best,
Bailey

[emmanuelparadis]

Hi @baileyjc ,
I suggest you try to cluster the haplotypes by first computing the distances among them (with dist.dna in ape, if you have the complete sequences, or dist.dna(, method = "N") which gives the Hamming distance), then you can do one of these:

• NJ tree (with nj() also in ape),
• simple hierarchical clustering (hclust()),
• multidimensional scaling (MDS, cmdscale()).

These will give you an idea whether you can reliably define clusters (you may have many pairs of haplotypes differing by a few sequence sites but they may form a "continuous chain").

You may also try rmst(): the new code does an exhaustive search of all MSTs which may help to define clusters.

Best,
Emmanuel

[baileyjc]

Thank you for your reply @emmanuelparadis. I took the route you suggested but unfortunately, I still end up with 154 haplotypes for 154 individuals. Instead of identifying by whole steps like previous, rmst identifies steps that are slightly different and instead it's just distances and does not resolve the issue of still having 154 unique haplotypes. I assume this has to do with the vcf file pulling out SNPs and if sequences in the bam file are of unequal length or lack coverage at loci then ambiguous bases are identified as SNPs in the vcf file. I saw another way to investigate this getting allele counts from ANGSD, creating the individual mtDNA genome for each sample and loading those into a fasta file. With the fasta file, I trimmed all mtDNA genomes to the same size given ambiguous bases in the beginning and end for certain individual genomes and then used pegas to process the fasta file. I was able to get 109 haplotypes from 154 samples. I still wish it would be possible to merge haplotypes with only 1 SNP difference between the nearest haplotype to simplify the network to identify the more unique patterns that I am after. Here is my network now:

haplotype_network.pdf

Thank you for your time, support, and creating such an incredible package.

[emmanuelparadis]

Maybe you can try to use hclust() and then cutree() to create the clusters of similar haplotypes. That could be simpler (although there seems to be some clusters visible with rmst()...)

[sjclare]

I wanted to add onto this since I'm looking for a similar feature to collapse highly similar haplotypes. Is there a way to use protein sequence, my aim would be to have any synonymous mutation removed, and keep non-synonymous mutations. Could this be done directly with AAbin or converting AABin back into DNAbin where amino acids only use specific codons (not completely accurate, but maybe an okay workaround)?

[emmanuelparadis]

@sjclare: you may have a look at the function trans() in ape. You can try with the woodmouse data:

R> library(ape)
R> data(woodmouse)
R> X <- trans(woodmouse, 2)
Message d'avis :
Dans trans(woodmouse, 2) :
  sequence length not a multiple of 3: 2 nucleotides dropped

While the 15 sequences are all different, some AA sequences are the same:

R> table(dist.aa(X))

 0  1  2  3  4  5  6  7  8 
22 30 13 16  9  9  4  1  1

It seems 22 pairs have identical cytochrome b. It's a matter of simple R coding to find them.

HTH

E.