We’re excited to introduce our Best Practices recommendations for calling variants on RNAseq data. These recommendations are based on our classic DNA-focused Best Practices, with some key differences in the early data processing steps, as well as in the calling step.

Best Practices workflow for RNAseq

This workflow is intended to be run per-sample; joint calling on RNAseq is not supported yet, though that is on our roadmap.

Please see the new document here for full details about how to run this workflow in practice.

In brief, the key modifications made to the DNAseq Best Practices focus on handling splice junctions correctly, which involves specific mapping and pre-processing procedures, as well as some new functionality in the HaplotypeCaller.

Now, before you try to run this on your data, there are a few important caveats that you need to keep in mind.

Please keep in mind that our DNA-focused Best Practices were developed over several years of thorough experimentation, and are continuously updated as new observations come to light and the analysis methods improve. We have only been working with RNAseq for a few months, so there are many aspects that we still need to examine in more detail before we can be fully confident that we are doing the best possible thing.

For one thing, these recommendations are based on high quality RNA-seq data (30 million 75bp paired-end reads produced on Illumina HiSeq). Other types of data might need slightly different processing. In addition, we have currently worked only on data from one tissue from one individual. Once we’ve had the opportunity to get more experience with different types (and larger amounts) of data, we will update these recommendations to be more comprehensive.

Finally, we know that the current recommended pipeline is producing both false positives (wrong variant calls) and false negatives (missed variants) errors. While some of those errors are inevitable in any pipeline, others are errors that we can and will address in future versions of the pipeline. A few examples of such errors are given in this article as well as our ideas for fixing them in the future.

We will be improving these recommendations progressively as we go, and we hope that the research community will help us by providing feedback of their experiences applying our recommendations to their data. We look forward to hearing your thoughts and observations!

Return to top


Are the authors of GATK BaseRecalibrator concerned about Post-transcriptional Modifications in RNA-seq? Reverse Transcriptases have difficulty correctly reading modified nucleotides. The Reverse Transcriptase may produce an error at a modified nucleotide when making the cDNA. Illumina will then read the resulting cDNA correctly (and give high quality score). Thus, even though the Illumina reads are correctly reporting the base in the cDNA (with high quality scores), it will be "wrong" compared to the reference, and not masked by dbSNP since it is only a Post-transcriptional modification. This will severely reduce the resulting empirical quality scores calculated by BaseRecalibrator. For DNA-seq, BaseRecalibrator masks "--knownSites" of polymorphism when calculating empirical Quality scores. The "--knownSites" is usually a VCF from e.g. dbSNP. In RNA-seq, do the authors of GATK recommend any kind of VCF with known Post-transcriptional modifications in mRNA?

Thu 6 Mar 2014


I'm afraid we don't have any recommendations for this -- in our hands BQSR performed normally on RNAseq data, but we haven't tested for this specifically. The size of potential effect is linked to how random vs. not random the errors might be, and what is the rate at which they occur. The more random and the lower the rate, the less noticeable any potential effect.

Thu 6 Mar 2014

At a glance

Follow us on Twitter

GATK Dev Team


RT @BroadGenomics: Missed #GATK WKSP at #ASHG16? 10am meet Geraldine at booth 329 - Broad’s GATK Guru @gatk_dev #BroadGenomicsExperts https…
20 Oct 16
RT @NJL_NGS: Broad Institute Workbench workshop now at #ashg16. https://t.co/t6L452mp8a
19 Oct 16
RT @konradjk: We're rebranding a bit! gnomAD now adds WGS regions, as well as doubling the data depth and increasing diversity (5K ASJ!) #A…
19 Oct 16
RT @NJL_NGS: https://t.co/mf0zBZrEfg now launched!
19 Oct 16
RT @dgmacarthur: For those who missed it: we just posted variant frequencies from 126,216 exomes and 15,136 genomes at https://t.co/uMAoxqh
19 Oct 16

Our favorite tweets from others

Asked a question, speaker jokes not to edit your BAM header and expect GATK to let you get away with it.
19 Oct 16
My new hobby: finding incomprehensible diagrams on office whiteboards and adding alarming conclusions to them https://t.co/7RAzBF0kYh
30 Sep 16
I've easily written my first custom ReadFilter for GATK. The @gatk_dev 's toolkit is a great example of programming.
21 Sep 16
@gatk_dev "make it so"
8 Sep 16
it's the nightly build owl for GATK :D https://t.co/OwTRrk6KHA https://t.co/rfmAbdIIQp
11 Aug 16
See more of our favorite tweets...
Search blog by tag

ad appistry ashg ashg16 benchmarks best-practices bug bug-fixed cancer cloud cluster cnv commandline commandlinegatk community compute conferences cram cromwell denovo depthofcoverage diagnosetargets error fix forum gatk3 gatk4 genotype genotype-refinement genotypegvcfs google grch38 gvcf haploid haplotypecaller hg38 holiday hts htsjdk ibm java8 job job-offer jobs license meetings mendelianviolations multithreading mutect mutect2 ngs nt outreach pairhmm parallelism patch performance phone-home picard pipeline plans ploidy polyploid poster presentations printreads profile promote release release-notes rnaseq runtime saas script search selectvariants sequencing service slides snow speed status sting support syntax talks team terminology third-party-tools topstory trivia troll tutorial unifiedgenotyper vcf-gz version-highlights versions vqsr wdl webinar workflow workshop