Zubair Khalid

Virologist/Molecular Biologist | Veterinarian | Bioinformatician

Conventional & Molecular Virology • Vaccine Development • Computational Biology

Dr. Zubair Khalid is a veterinarian and virologist specializing in conventional and molecular virology, vaccine development, and computational biology. Dedicated to advancing animal health through innovative research and multi-omics approaches.

Dr. Zubair Khalid - Veterinarian, Virologist, and Vaccine Development Researcher specializing in Computational Biology, Multi-omics, Animal Health, and Infectious Disease Research

Category: Guides

Gatk Rnaseq Variant Calling

GATK RNA seq variant calling is the process of identifying single nucleotide variants (SNVs) and small insertions or deletions (indels) from RNA sequencing data using the Genome Analysis Toolkit (GATK) best practices. This guide explains the core concepts, decision points, and a practical workflow. It is written for bioinformatics researchers, computational biologists, and graduate students who need a source bounded framework for variant discovery from transcriptomes. Use this guide when you have aligned RNA seq reads and need to call variants with confidence while understanding the method's limits.

RNA seq variant calling differs from DNA variant calling because the input is derived from spliced transcripts. Splice aware aligners like STAR or HISAT2 must be used, and the calling pipeline must handle reads spanning exon exon junctions. The GATK Best Practices for RNA seq provide a robust workflow, but careful tuning of filters is required to reduce false positives from sequencing errors and transcriptomic noise. For a complete overview of GATK tools and their rationale, see the Galaxy Training Network which offers hands on tutorials.

Before starting, ensure you have a high quality reference genome (e.g., GRCh38 or T2T CHM13) and RNA seq data in FASTQ format from a reliable source such as the NCBI Sequence Read Archive. This guide assumes you are working with paired end short reads, though many steps apply to long read data after appropriate transformation. For a broader context of variant calling challenges, read EMBL EBI Training materials on RNA seq analysis.

At a Glance

Aspect Detail
Purpose Detect SNVs and small indels from RNA seq data
Key Tools GATK (HaplotypeCaller, VariantFiltration), STAR, Picard
Input Aligned RNA seq reads (BAM files)
Output VCF file with filtered variant calls
Main Challenge High false positive rate from splicing artifacts and RNA editing
Recommended Coverage 20 30x (transcriptome wide)
Time Estimate hours to days depending on data size

Core Concepts and Decision Points

GATK RNA seq variant calling relies on a few core concepts: splice aware alignment, base quality score recalibration (BQSR), and variant filtering. BQSR is often skipped for RNA seq because systematic biases differ from DNA data. Instead, hard filtering is used. The decision to use HaplotypeCaller in RNA seq mode is critical. This mode reassembles reads in active regions including those that span splice junctions, reducing false calls at exon boundaries. See the Bioconductor documentation for additional R based quality assessment tools that complement GATK.

Decision points include:

  • Read Group Information: You must assign read groups (sample, library, platform) in the BAM file. Without it, GATK will fail. Use Picard AddOrReplaceReadGroups.
  • Duplicate Marking: Mark optical and PCR duplicates using MarkDuplicates. Duplicates in RNA seq often come from library amplification and should be flagged but not removed because the same transcript can produce identical reads.
  • SplitNCigarReads: This GATK tool splits reads at exon junctions and hard clips bases that overhang introns. It is essential for accurate variant calling.
  • Base Recalibration: Often omitted because RNA editing and sequencing biases are not easily modeled by BQSR. Many workflows apply hard filters instead.
  • Filtering Thresholds: Common filters include Fisher strand bias (FS > 30.0), QualByDepth (QD < 2.0), and ReadPosRankSum (ReadPosRankSum < -8.0). Adjust based on your data.

For a deeper review of challenges in RNA seq variant calling, read this comprehensive review which covers sources of false positives and strategies to mitigate them.

Practical Workflow Steps

The following workflow assumes you have raw FASTQ files from an RNA seq experiment.

Step 1: Alignment with a Splice Aware Aligner

Use STAR to align reads to the reference genome. Generate a genome index with the splice junction database from the GTF annotation file. Run STAR with two pass mode to improve junction detection. Example parameters:

STAR --genomeDir /path/to/index --readFilesIn sample_R1.fastq sample_R2.fastq --outSAMtype BAM SortedByCoordinate --outSAMunmapped Within --outSAMattributes Standard

Step 2: Preprocessing with Picard and GATK

Add read groups using Picard:

java -jar picard.jar AddOrReplaceReadGroups I=Aligned.sortedByCoord.out.bam O=rg_added.bam RGID=id RGLB=library RGPL=ILLUMINA RGPU=unit RGSM=sample

Sort the BAM, index it, and mark duplicates:

java -jar picard.jar MarkDuplicates I=rg_added.bam O=dedup.bam M=dedup_metrics.txt CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT

Step 3: SplitNCigarReads and Read Processing

Run GATK SplitNCigarReads to split reads at exons and clip overhanging intronic sequences:

gatk SplitNCigarReads -R reference.fasta -I dedup.bam -O split.bam

Recalibrate base qualities only if you have known variant sites for RNA (rarely used). Most workflows skip this step.

Step 4: Variant Calling with HaplotypeCaller

Run HaplotypeCaller in RNA seq mode by adding --dont-use-soft-clipped-bases and --emit-ref-confidence GVCF if you plan to combine samples. For a single sample:

gatk HaplotypeCaller -R reference.fasta -I split.bam -O raw_variants.vcf --dont-use-soft-clipped-bases true --standard-min-confidence-threshold-for-calling 20

Step 5: Hard Filtering

Apply hard filters to reduce false positives from RNA editing and sequencing errors:

gatk VariantFiltration -R reference.fasta -V raw_variants.vcf -O filtered_variants.vcf --filter-name "FS30" --filter-expression "FS > 30.0" --filter-name "QD2" --filter-expression "QD < 2.0" --filter-name "ReadPosRankSum" --filter-expression "ReadPosRankSum < -8.0"

You may also filter variants with low depth (DP < 10) or high strand bias. Adapt thresholds after evaluating your data quality.

Step 6: Annotation and Interpretation

Annotate the filtered VCF using SnpEff or VEP to determine the functional impact of variants. For example, missense variants in coding regions are more reliable than intronic calls. For a study on somatic mutations in cancer using RNA seq, see this paper which highlights the use of the T2T CHM13 reference and careful filtering.

Quality Checks and Common Mistakes

Quality checks after calling:

  • Compare the transition/transversion ratio (Ti/Tv) for known exonic variants. Expected ratios are around 2.0 to 3.0 for human exome data. Lower ratios suggest many false positives.
  • Examine the distribution of variant allele frequency (VAF). RNA seq VAF can be skewed by allelic expression. Expect peaks near 0.5, 1.0, but also alleles with low expression.
  • Verify known variants from dbSNP overlap. A high proportion of novel calls without coverage support indicates false positives.

Common mistakes include:

  • Not marking duplicates. This can inflate variant confidence scores.
  • Using default DNA filters for RNA seq. Strand bias filters are critical because RNA seq has strong strand specific signals.
  • Ignoring splice junction regions. Variants called near exon edges often result from misalignment. Use SplitNCigarReads properly.
  • Over filtering. Too aggressive filtering can remove true positive variants. Validate candidate variants with orthogonal methods like Sanger sequencing.
  • Using a poor quality reference genome. The improved somatic mutation profiling study shows that a complete reference like T2T CHM13 improves calling accuracy.

Limits and Interpretation

RNA seq variant calling has inherent limits compared to DNA calling. Variants in genes with low expression may be missed because of insufficient coverage. RNA editing events (e.g., A to I conversions) can appear as mismatches and lead to false variant calls. You cannot easily distinguish germline variants from somatic mutations or RNA editing without a paired DNA sample. For circulating tumor DNA detection from RNA, see this targeted RNAseq study which uses a panel design to improve sensitivity.

Interpretation must account for allele specific expression and alternative splicing. A variant seen in an exon may be absent in some transcript isoforms. Always check expression levels of the transcript context. Deep learning based callers offer some improvements, as described in this deep learning germline caller, but they still require careful validation.

Long read RNA seq data poses additional challenges. Transforming alignment files can help, as shown in this study on long read data. For multi omics integration, see this gastric cancer cell line analysis. Always remember that RNA seq variant calling is a tool for discovery, not diagnosis. Confirm actionable variants with DNA sequencing.

Frequently Asked Questions

Why does GATK require SplitNCigarReads for RNA seq?

SplitNCigarReads splits reads at known splice junctions and hard clips bases that extend into introns. This prevents the caller from treating intronic bases as mismatches and reduces false positive variants near exon boundaries.

Can I use GATK HaplotypeCaller with default parameters for RNA seq?

No. You must set --dont-use-soft-clipped-bases true and adjust the confidence threshold. Default parameters are designed for DNA and will produce many false calls from RNA seq.

How do I distinguish germline from somatic variants in RNA seq?

Without a matched normal DNA sample you cannot reliably separate germline from somatic variants. Consider comparing variant allele frequencies with known dbSNP entries. Low VAF variants may be somatic but could also reflect technical noise or RNA editing.

What coverage is sufficient for RNA seq variant calling?

A minimum of 20x transcriptome wide coverage is recommended. However, lowly expressed genes may remain underpowered. Targeted RNA seq panels with higher depth improve detection, as shown in this study.

References and Further Reading

Related Articles