SEQanswers

Go Back   SEQanswers > Applications Forums > Genomic Resequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
Variant calling: reference positions with no reads don't appear in vcf file MattRSmith Bioinformatics 0 05-05-2017 08:06 AM
Calling structural variants (CNVs) with single-end reads agwe Genomic Resequencing 5 01-18-2016 07:30 AM
Mapping Amplicon Reads to Reference svos Bioinformatics 2 04-24-2015 12:12 AM
Mapping reads to reference genome + count reads of genes cumulonimbus RNA Sequencing 12 10-02-2013 08:07 AM
Combined mapping of RNA-Seq reads originating from multiple species schelhorn RNA Sequencing 7 11-05-2010 08:55 AM

Reply
 
Thread Tools
Old 05-18-2017, 03:15 PM   #1
gruberjd
Member
 
Location: Michigan

Join Date: Nov 2011
Posts: 12
Question Calling variants when mapping reads from different species than reference

Hi all,

My SNP calling pipeline seems to be poorly performing when tasked with mapping reads from a sample that is from a related species from the reference assembly. (This sort of thing is common in plant breeding.) Obviously the more variants that differentiate reference and sample, the harder it is to map reads correctly, so it makes sense that this is harder than finding intraspecies SNPs.

Nonetheless, I am 100% confident it can do better than it is doing. I can see from the bam that there are reads mapping to a given region, and even then most of the time my SNP caller (the last free version of GATK) fails to emit genotypes, even of sites that simply match the reference. Interestingly, the few times that it does emit a genotype, it is highly highly biased towards A&T calls.

Is any of this a classic symptom of a parameter I need to adjust, or a sign that I need to switch SNP callers? Can anyone cite literature that compares SNP callers' performance (or influence of different parameters) in this particular challenging task? I haven't had much luck searching.

Thanks!

Jonathan
gruberjd is offline   Reply With Quote
Old 05-18-2017, 05:16 PM   #2
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,553
Default

I don't have any papers, but I'd be interested in seeing if you get better results using BBMap for mapping and BBMap's variant caller for variant-calling. BBMap can align at short reads at quite low identity, so it's good for cross-species alignment (you can increase sensitivity with the "minid" flag). Assuming you have paired interleaved reads in a single file, and starting with the raw reads, the commands would be something like:

Code:
bbduk.sh in=reads.fq out=clean.fq ref=adapters.fa ktrim=r k=23 mink=11 hdist=1 tpe tbo
bbmap.sh in=clean.fq ref=ref.fa out=mapped.sam
callvariants.sh in=mapped.sam ref=ref.fa ploidy=2 out=vars.vcf
For multiple samples, the callvariants command would be:

Code:
callvariants.sh in=sample1.sam,sample2.sam,sample3.sam ref=ref.fa ploidy=2 out=vars.vcf multisample
At JGI we are using CallVariants a lot these days for various strains of things mapped to a divergent reference (largely Aspergillus, E.coli, and some plants from cross-breeding experiments). I don't really do variant-calling in production any more (just for testing and development) but the people who do at JGI switched to CallVariants because, as they told me, it performs much better than their prior pipelines (GATK, FreeBayes, and a couple others), particularly for variants with low allele fractions (using the "rarity" flag). In general I think it's best to use the defaults, though of course you need to specify ploidy correctly.
Brian Bushnell is offline   Reply With Quote
Old 05-19-2017, 11:49 AM   #3
gruberjd
Member
 
Location: Michigan

Join Date: Nov 2011
Posts: 12
Default

Brian-- Thanks. I'll give BBMap a shot and will certainly report success if I find it!
gruberjd is offline   Reply With Quote
Old 05-19-2017, 11:52 AM   #4
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,553
Default

Great - I look forward to your feedback.
Brian Bushnell is offline   Reply With Quote
Reply

Tags
gatk, snp, variant caller

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:37 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO