[chiayi]

For maaping diploid maize data, I concatenated all the chromosomes for a given species (e.g. 10 chr in maize) and used it as the mapping reference. I then mapped the BBDuk-trimmed reads to this reference (genome, not transcriptome) using the default setting of BBMap.
For the maize data generated from B73 (B73 is the sequeced ref genome), the ambig rate rages from 3% to 15%. I saw the variation occurred between biologica rep, not within a pair. For example, the ambig rates of Read 1 and Read 2 of a replicate are close to each other, but the ambig rates of different biologica replica could be various, (12%, 3%, 3%; 11%, 11%, 3%). The variatoin between replica almost made me feel like it was due to the technical issue of the libraries and/or biological nature of maize. I could try to map it to the transcriptome and see if there's any improvement. Other than that I wasn't sure how to trace down, and/or if this should be a concern.
For a different set of maize data with mixed background, as for now I still used the same B73 base reference genome for mapping. The ambig ranges went up to 6% - 46%. The increase was expected given the SNPs/indels present between different cultivars. I also observed variation between biological replicates similar to described above.
On a related note, I also have diploid Arabidopsis data which I also mapped to its own genome (TAIR9 genome fasta). I added a -maxindel=2000 to accomodate the compact size of the genome. The ambig rate on average was lower (mostly below 4%). I did not see that radical variation between replicates either. This is consistent with my guess about the ambig rates in maize was partly due to the nature of maize.
Please let me know if any furtehr details would be helpful for you to diagnose. Thank you as always for your input.

[Brian Bushnell]

Normally... I see ambig rates of 3% or lower in diploids like human, and 0.1% or lower in haploid bacteria. I have little experience in mapping RNA to plant genomes, aside from feedback from co-workers. And they mainly map to plant transcriptomes rather than genomes.

Can you describe your mapping protocol in more detail? For example, are you concatenating all genomes and mapping to all simultaneously, or are you experiencing a high claimed ambig rate when mapping to a single reference alone?

[chiayi]

BBMap is very tolerant of SNPs and indels so generally you don't need to do any kind of correction, but aligning to a SNP-corrected reference (assuming the SNPs are homozygous, or at least a majority) will be more accurate than aligning to the base reference.

I share your thoughts so I will at least try a few background to see how it goes. On the other hand, I'm not sure how to compare the accuracy between the results generated from base reference vs SNP corrected reference. Do you have any suggestion?

I'm not really sure where the high ambig rate is coming from. Is this WGS, or are you doing some kind of enrichment, RNA-seq, etc? And are you mapping to the genome or transcriptome?

The reads were generated by NEBNext Ultra RNA library prep kit for Illumina (no enrichment) and sequenced on NextSeq. I mapped the reads to the latest genome annotation (B73 RefGen_v4 reference genome). I should also mention that these are mixed background (hybrid of B73, reference genome, and other cultivar). For B73 alone, the ambig rate was lower at about ~10%.
Is there a 'standard' for ambiguous rate when you look at the stats from BBMap?

[Brian Bushnell]

Maize is the kind of organism where de-novo assembly is extremely difficult. I recommend avoiding that; it will make things much more complicated. I don't think there's any reason to do that, either, unless you find that structural variations are causing problems.

BBMap is very tolerant of SNPs and indels so generally you don't need to do any kind of correction, but aligning to a SNP-corrected reference (assuming the SNPs are homozygous, or at least a majority) will be more accurate than aligning to the base reference.

I'm not really sure where the high ambig rate is coming from. Is this WGS, or are you doing some kind of enrichment, RNA-seq, etc? And are you mapping to the genome or transcriptome?

[chiayi]

Hi Brian and other members,

I'm trying to map the Illumina reads generated from 19 inbred and hybrid mazie cultivars under two conditions for downstream differential expression analysis. Considering the SNPs and indels in the mixed genetic background, my original plan was to generate a 'corrected refernece' using tools like Quiver or Pilon. Then I came across to BBMap and was happy to find its capability to handle SNPs and indels. I used the default setting and the mapped rate so far was over 97%. Meanwhile, I did notice a huge variaion of ambiguous rate, ranging from 6% to >50%. Therefore, I was wondering if it is possible to evaluate which method (directly map to a single ref vs map to 19 SNP corrected ref) would be more accurate for the DE analysis.

The third option would be to re-assemble de-novo assembly for each genetic background. However, I wasn's sure if it's possible to come up with a consensus contig conrrespondance (contigX_background1, contigX_background2..contigX_background19) so that I could monitor each contig/gene of interest in multiple backgrounds.

Any input/thoughts you may have on this would be much appreciated. Thank you in advance.

[skbrimer]

Thanks for the input guys! I will start playing with the updated BBmap as well!

[Brian Bushnell]

Thank you Brain!

Is there a variant calling pipeline that works best with BBmap?

Funny you should ask! I just finished developing a variant-caller, because I was unhappy with the speed, usability, and some other aspects of Samtools mpileup and GATK (I've never used Freebayes, though). Also, it has a default ploidy of 1 which better reflects projects at JGI, though the ploidy can be changed to an arbitrary number.

I just uploaded a new version of BBTools (36.71) with the variant caller (older versions do not support vcf output). You can run it like this:

callvariants.sh in=mapped.sam vcf=vars.vcf ref=ref.fasta ploidy=1

CallVariants works on sam or bam, sorted or unsorted, and it's really fast. Also, it does not care what format the cigar strings are in or whether you trimmed sequence names after the first whitespace. I'd be eager to hear how it compares to FreeBayes, if you'd like to try it out. So far, other people at JGI have indicated that it does a much better job of identifying low-depth variants in high-ploidy organisms than Samtools or GATK (if you adjust the ploidy and minallelefraction flags).

That said, no, I don't have a recommended variant-calling workflow for BBMap, because I have not done rigorous benchmarking of all the possibilities yet. Oh, and to modify my previous post, for variant calling I suggest you also add the "mdtag" flag, which some variant-callers like to have. So, for Freebayes:

bbmap.sh in=reads.fq ref=ref.fa out=mapped.sam mdtag trd sam=1.3

BBMap suite now has a variant caller (callvariants.sh) :-)

Indeed But if you use it, I highly recommend downloading BBMap 36.71, as prior versions lack proper VCF output.

[GenoMax]

If you are going to re-do mapping then clean up the fasta ID's (and regenerate BBMap index) in your reference genome, if there are spaces in the words there.

BBMap suite now has a variant caller (callvariants.sh) :-)

[skbrimer]

Thank you Brain!

Is there a variant calling pipeline that works best with BBmap?

[Brian Bushnell]

Unfortunately, it sounds like Freebayes is trimming the names of the sequences after the first whitespace when loading them from the fasta file, but it's not doing that when looking at the sam file. To fix that problem, unfortunately you'd need to rerun the mapping, adding the "trd" flag (trim read description). Additionally, Freebayes is not compliant with the current Sam specification and does not understand cigar strings with "=" and "X" symbols, so you need to output in legacy Sam format with M symbols. Therefore, for Freebayes, you need to add two extra flags when mapping:

trd sam=1.3

[skbrimer]

Hi Brian,

I like to use freebayes for the vcf creation because all of our stuff is haploid, however I do not know how to work with the ref file that BBMap makes during a mapping. Freebayes accepts a fasta for the ref. The ref folder contains a compressed chromosome file which freebayes doesn't like and when I use the fasta that I feed to BBmap it doesn't like the header mismatch.

Can you tell me how to overcome this issue?

[Brian Bushnell]

Can you confirm that pileup doesn't need sam files to be sorted ?

Yes, that's correct.

[Brian Bushnell]

Pileup uses the header lines to determine the names and lengths of the sequences in the reference. So, they must be present. However, it doesn't really care whether there are extra ones randomly in the file; they just get ignored. So, cat works fine.

[Brian Bushnell]

Pileup won't accept multiple input files, but it will accept standard in. So if the files were sam rather than bam, you could do this:

cat lib1.sam lib2.sam | pileup.sh in=stdin.sam out=stats.txt hist=histogram.txt

or for gzipped files:

cat lib1.sam.gz lib2.sam.gz | pileup.sh in=stdin.sam.gz out=stats.txt hist=histogram.txt

However, I don't see a way to do that with bam files. BBMap will produce sam files instead of bam files if you name the output *.sam instead of *.bam.

[Brian Bushnell]

BBMap can output coverage, by the way. Just add the flag "covstats=covstats.txt". You may also want to add the flag "delcov=f" for this purpose.