I have taken Illumina reads and aligned them to a reference genome using BWA then obtained the corresponding BAM/SAM files. I have also called SNPs, which are in VCF format, and tried to use this file to predict synonymous and nonsynonymous sites (using snpEff), but this will only give me a N/S ratio and I what I really want is the dN/dS ratio. Is there any way do this from the BWA alignments? I am new to NGS genome assembly, so any tips are much appreciated.
From Wikipedia:
In genetics, the Ka/Ks ratio (or ω, dN/dS), is the ratio of the number of Nonsynonymous substitutions per non-synonymous site (Ka) to the number of synonymous substitutions per synonymous site (Ks).
So you additionally need to know the number of non-synonymous sites and the number of synonymous sites.
Hey,
it seems you are calling rather Pn and Ps (number of synonymous and nonsyn polymorphic sites) denn divergent sites with your method. However, dNdS is highly depended on alignment quality. I wouldn't trust bwa alignments myself for calculating such sensitive measurements as dN and dS.
Best,
axel
Hi Stacy,
I have done similar thing. First of all, you need to know where your genes are (GTF file). You are interested in coding sequence (CDS). Based on you VCF file, you will know whether given variant is homozygous/heterozygous, and when you compare coordinates with GTF, if variant falls into coding sequence and whether it causes amino substitution or not.
Hope this helps,
L.
Thanks everyone for the helpful responses. Unfortunately, there does not appear to be GTF files for the genomes I am looking at, so I do not know if this approach will work for me. However, Axel does make a good point about the importance of the alignments -- since I am looking at very similar genomes (i.e. different strains/isolates), I have been able to reconstruct the gene sequences from the VCF file so that I can then generate good alignments for dN/dS calculation. I am hoping this will work.
If you know how to translate your genome, and if you have a sufficiently large number of mutations, you might do the following: (1) calculate substitution frequency; (2) randomly create the same number of mutations, maintaining the substitution frequencies; (3) calculate N/S ratio again. Repeat 1-3 a thousand times and check how your observed N/S ratio ranks against the randomized mutations.
I think I have an R tool somewhere that can do that.
Hi, Try using the tool Annovar (http://www.openbioinformatics.org/annovar/) which can annotate the SNP as synonymous and non-synomous, from this you can always calculate the dN/dS ratio.
Hi Eithan -- well, I am a newbie with genome assembly, so the farthest I have gotten is the VCF file, but this is not an ideal format for translation of the genome. I am working with Plasmodium genomes (malaria parasite), so there is a lot of low complexity sequence and extremely high AT content. I have a reference genome, so what is the best way to obtain the corresponding genes for the new genome? My solution has been to use reconstruct the gene sequences from the VCF file...
Srinivas, from the defintion of dN/dS, I am not sure you can do this with just the SNP data since in order to calculate dN/dS, you need to know the number of synonymous and nonsynonymous sites -- this requires the multiple sequence alignment of the genes, which I currently don't have.
Hi Rohan, thanks for replying, but unfortunately I am not working with human genomes, so I'm not that confident your suggestion will work. I am looking at Plasmodium falciparum genomes, which are much more difficult to work with because of their extremely high AT content and extensive amounts of low complexity sequence.
- Badges
- Science topic
- Similar topics
- Computer Science and Engineering
- Bioinformatics
- Bioinformatic Tools