View Single Post
Old 05-18-2017, 05:16 PM   #2
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
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