We discovered today that we made an error in the documentation article that describes the RNAseq Best Practices workflow. The error is not critical but is likely to cause an increased rate of False Positive calls in your dataset.

The error was made in the description of the "Split & Trim" pre-processing step. We originally wrote that you need to reassign mapping qualities to 60 using the ReassignMappingQuality read filter. However, this causes all MAPQs in the file to be reassigned to 60, whereas what you want to do is reassign MAPQs only for good alignments which STAR identifies with MAPQ 255. This is done with a different read filter, called ReassignOneMappingQuality. The correct command is therefore:

java -jar GenomeAnalysisTK.jar -T SplitNCigarReads -R ref.fasta -I dedupped.bam -o split.bam -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS

In our hands we see a bump in the rate of FP calls from 4% to 8% when the wrong filter is used. We don't see any significant amount of false negatives (lost true positives) with the bad command, although we do see a few more true positives show up in the results of the bad command. So basically the effect is to excessively increase sensitivity, at the expense of specificity, because poorly mapped reads are taken into account with a "good" mapping quality, where they would normally be discarded.

This effect will be stronger in datasets with lower overall quality, so your results may vary. Let us know if you observe any really dramatic effects, but we don't expect that to happen.

To be clear, we do recommend re-processing your data if you can, but if that is not an option, keep in mind how this affects the rate of false positive discovery in your data.

We apologize for this error (which has now been corrected in the documentation) and for the inconvenience it may cause you.


sirian


Thanks for the correction! I was actually wondering a little bit why you changed every score.

Wed 11 Jun 2014

sboyle


Thanks for the correction Geraldine!

Wed 11 Jun 2014

kam


Here's a working link for [ReassignOneMappingQuality](https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_engine_filters_ReassignOneMappingQualityFilter.php). Have you considered using the mapping quality formula proposed by the authors of subread? This mapping quality score works for any read mapper. See page 20 in the [Subread User's Guide](http://bioinf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf).

Wed 11 Jun 2014

Geraldine_VdAuwera


@kam This may be a good recommendation to make to the authors of the mappers.

Wed 11 Jun 2014

justinjj


Dear Geraldine, Could you please clarify me is there any difference or issue if I use the mapq as 50 instead of 60 as suggested to run GATK? The bwa aligner higher quality is 60 but the tophat2 (uses bowtie2) provide higher quality/unique mapq as 50 and GATK runs without any error in both cases also by star aligner "--outSAMmapqUnique 50" The RNAseq mappers is already giving meaningful quality score? https://software.broadinstitute.org/gatk/guide/article?id=3891 (So we use the GATK’s ReassignOneMappingQuality read filter to reassign all good alignments to the default value of 60. This is not ideal, and we hope that in the future RNAseq mappers will emit meaningful quality scores, but in the meantime this is the best we can do.) Thanks.

Wed 11 Jun 2014

Geraldine_VdAuwera


Yes, that's fine. If newer versions of the mappers produce MAPQ scores in that range, there is no need to reassign a different value.

Wed 11 Jun 2014

justinjj


Thanks Geraldine.

Wed 11 Jun 2014




At a glance



Follow us on Twitter

GATK Dev Team

@gatk_dev

RT @broadinstitute: .@NIH funds available for investigators in need of cloud computing & storage resources. Apply: https://t.co/IoiZMUNBM8
26 Apr 17
@DataKimist Enjoy! And let us know if we can help.
19 Apr 17
@mjpchaisson Not meant that way - just depending on what you're doing you may want to cite earlier framework or lat… https://t.co/QpIbwRf0bC
18 Apr 17
@cabioinformatic For more recent versions see https://t.co/QCbos5KBWw
15 Apr 17
@thatdnaguy @notigor @David_McGaughey @brent_p Indel Realign is redundant with assembly-based realign done by HC, w… https://t.co/77Lyil7BJY
13 Apr 17

Our favorite tweets from others

best error output: Please do NOT post this error to the GATK forum unless you have really tried to fix it yourself.
4 Apr 17
From the @gatk_dev page describing .vcf files: "Don't write home-brewed VCF parsing scripts. It never ends well” https://t.co/28KcRoV14j
28 Feb 17
Our 3-day course on GATK https://t.co/mtN60KRTyS finished - 38 participants very happy! Big thanks to @gatk_dev team for excellent lessons.
24 Feb 17
@froggleston @dgmacarthur Sounds like ExAC is reaching Uber stage. ‘Uber but for pizza’. ‘ExAC but for wheat’.
14 Jan 17
#ESRenpeinture grad school - postdoc - after postdoc https://t.co/o3vQMgBDgk
6 Jan 17
See more of our favorite tweets...
Search blog by tag

appistry ashg ashg16 benchmarks best-practices bug bug-fixed cloud cluster cnv collaboration community compute conference conferences cram cromwell depthofcoverage diagnosetargets error forum gatk3 gatk4 genotype-refinement genotypegvcfs google grch38 gvcf haploid haplotypecaller help hg38 holiday hts htsjdk ibm intel java8 job job-offer jobs license meetings mutect mutect2 ngs outreach pairhmm parallelism patch pdf performance picard pipeline plans ploidy polyploid poster presentations printreads profile promote release release-notes rnaseq runtime saas script sequencing service slides snow speed status support talks team terminology topstory troll tutorial unifiedgenotyper vcf-gz version-highlights wdl workflow workshop xhmm