# Gatk Rna Seq Variant Calling


## Key Takeaways

- GATK HaplotypeCaller, specifically in its RNA-seq mode with `ERC GVCF` and `--dont use soft clipped bases` flags, is the recommended engine for variant calling from transcribed sequences. This approach accounts for the unique characteristics of RNA data, such as spliced alignments.
- Essential preprocessing steps for RNA-seq variant calling include using a spliced-aware aligner (e.g., STAR, HISAT2), running `SplitNCigarReads` to handle exon-intron boundaries, and recalibrating base quality scores with `BaseRecalibrator` and `ApplyBQSR`.
- Deduplication of reads in RNA-seq variant calling is a nuanced decision; while standard deduplication is common in DNA-seq to remove PCR artifacts, it can inadvertently remove true biological signals from highly expressed genes in RNA-seq. Many pipelines opt to skip or leniency apply deduplication for variant discovery.
- Hard filtering is generally recommended over VQSR for RNA-seq variant calling due to typically lower variant counts, with specific filters like `QD < 2.0`, `FS > 30.0`, and `ReadPosRankSum < 8.0` being empirically derived and crucial for reducing false positives.
- Key limitations of RNA-seq variant calling include expression-level bias leading to missed variants in lowly expressed genes, false positives arising from mapping errors around splice junctions, and potential masking of heterozygous variants due to allele-specific expression.
- Quality control metrics for RNA-seq variant calls include evaluating the transition/transversion ratio (expected around 2.0-2.2 for human), ensuring sufficient read depth at variant sites, and comparing novel variant calls against known databases like dbSNP to assess systematic errors.

---

Variant calling from RNA sequencing data using the Genome Analysis Toolkit (GATK) is a specialized process that identifies single nucleotide variants and small insertions or deletions from transcribed sequences. This guide provides a practical, source bounded framework for researchers who want to call variants from RNA seq data using GATK best practices. You should use this guide if you have aligned RNA seq reads and need to detect germline or somatic variants, understand the specialized preprocessing steps required for RNA, and interpret results within the limitations of transcriptome wide data. For authoritative background on variant calling principles, see [NCBI Bookshelf](https://www.ncbi.nlm.nih.gov/books/). Detailed training on RNA seq analysis is available from [EMBL EBI Training](https://www.ebi.ac.uk/training/).

## At a Glance

| Aspect | Key Point |
|--------|-----------|
| Core tool | GATK HaplotypeCaller in RNA seq mode |
| Preprocessing requirement | SplitNTrim, reassign mapping qualities, deduplication |
| Alignment requirement | Spliced aware aligner (STAR, HISAT2) |
| Output | VCF with variant annotations |
| Major limitation | Expression level bias, no intronic variants, false positives from mapping errors |
| Recommended quality filter | FisherStrand > 30, QualByDepth < 2.0, ReadPosRankSum |

## Core Concepts and Decision Points

RNA seq variant calling differs fundamentally from DNA variant calling. The transcriptome only covers expressed regions, meaning coverage is highly variable between genes and samples. Alignment must account for splice junctions, and GATK uses a special mode that treats reads as RNA seq evidence rather than DNA evidence. The core decision points involve your choice of aligner, whether to perform deduplication (more on that later), and how to handle strandedness.

GATK's HaplotypeCaller with the ` ERC GVCF` and ` --dont use soft clipped bases` flags is the recommended engine for RNA seq. A comprehensive review of variant calling tools for RNA seq notes that GATK remains widely used despite challenges posed by non uniform coverage and splicing artifacts [A comprehensive review of variant calling tools for RNA seq: Challenges and advances](https://pubmed.ncbi.nlm.nih.gov/42081940/). You must decide whether to call variants per sample individually or jointly with multiple samples. Joint calling improves sensitivity for low frequency variants but increases computational cost. For somatic variant detection, consider using Mutect2 in RNA seq mode, but be aware that RNA seq somatic calling has higher false positive rates due to expression noise. Recent work using the T2T CHM13 reference improved somatic mutation profiling in small cell lung cancer via RNA seq [Improved somatic mutation profiling in small cell lung cancer using RNA sequencing and the T2T CHM13 genome reference](https://pubmed.ncbi.nlm.nih.gov/41485300/). This suggests that reference genome choice matters.

Another decision is handling of duplicate reads. In DNA seq, PCR duplicates are removed to avoid false variant calls. In RNA seq, duplicate reads often arise from biological expression rather than PCR artifacts, so standard deduplication can remove true variant signals from highly expressed genes. However, for variant discovery, many pipelines skip deduplication or apply lenient duplicate removal. You must weigh the trade off: higher sensitivity versus increased false positives from PCR duplicates.

## Practical Workflow or Implementation Steps

### Step 1: Align RNA seq reads with a spliced aware aligner

Use STAR or HISAT2 to generate coordinate sorted BAM files. STAR is preferred for speed and accuracy. Include the `--outSAMtype BAM SortedByCoordinate` option. Your reference must be compatible with GATK. Index the BAM file with `samtools index`. Detailed tutorials for alignment are provided by the [Galaxy Training Network](https://training.galaxyproject.org/).

### Step 2: Preprocess the BAM file with GATK tools

Run `SplitNCigarReads` to break reads at exon intron boundaries and reassign mapping qualities. This step is mandatory because GATK expects continuous alignments. The command is:

```
gatk SplitNCigarReads R reference.fa I aligned.bam O split.bam
```

Then run `BaseRecalibrator` and `ApplyBQSR` to recalibrate base quality scores. RNA seq base recalibration uses known variant sites from dbSNP. Follow GATK best practice for the known sites argument. After recalibration, mark duplicates if you decide to remove them. Use `MarkDuplicates` with `REMOVE_DUPLICATES=false` to produce a metrics file without deleting reads. For most RNA seq variant calling, keeping duplicates is recommended.

### Step 3: Call variants with HaplotypeCaller in RNA seq mode

Use the following command:

```
gatk HaplotypeCaller R reference.fa I recalibrated.bam O raw variants.g.vcf.gz ERC GVCF --dont use soft clipped bases true --standard min confidence threshold for calling 30
```

The ` --dont use soft clipped bases` flag prevents false variants from alignment artifacts. The `standard min confidence threshold for calling` sets the minimum phred scaled confidence for a call. RNA seq requires higher thresholds than DNA due to higher error rates. You can also output a regular VCF by omitting `ERC GVCF`.

### Step 4: Filter variants

Hard filtering is recommended for RNA seq because VQSR (Variant Quality Score Recalibration) typically fails due to the smaller number of variant sites. Generate filter expressions in GATK's `VariantFiltration`:

```
gatk VariantFiltration R reference.fa V raw variants.vcf O filtered.vcf filter name "RNAseq filters" --filter expression "QD < 2.0" --filter expression "FS > 30.0" --filter expression "SOR > 3.0" --filter expression "MQRankSum < 12.5" --filter expression "ReadPosRankSum < 8.0"
```

These thresholds are empirical and may need adjustment based on your data. A study on variant discovery in human cells exposed to mycotoxins used similar filters and reported good precision [Assessment of single nucleotide variant discovery protocols in RNA seq data from human cells exposed to mycotoxins](https://pubmed.ncbi.nlm.nih.gov/36016515/).

### Step 5: Annotate and export

Use `SnpEff` or `ANNOVAR` to annotate variant effects. For downstream analysis in R, use the `VariantAnnotation` package from [Bioconductor](https://bioconductor.org/). Bioconductor provides extensive documentation on reading VCF files and performing statistical tests on variant calls.

## Quality Checks

After filtering, perform several quality checks. First, verify the number of variants per chromosome. RNA seq variant density should reflect gene density, not copy number. Look for an excess of transitions over transversions (Ti/Tv ratio). For human RNA seq, Ti/Tv around 2.0 to 2.2 is expected. Higher ratios suggest false positives from sequencing errors. Second, check the distribution of read depths at variant sites. Variants called at very low depth (less than 5 reads) are unreliable. Third, compare your variant calls to known SNPs from dbSNP. A high proportion of novel variants may indicate systematic alignment errors.

For germline calling, evaluate heterozygosity rates. RNA seq often shows an excess of homozygous calls due to allele specific expression. This is a biological artifact, not a technical error. Confirm with orthogonal DNA data if available. Tools available through [EMBL EBI Training](https://www.ebi.ac.uk/training/) can help with quality metrics. Also consider using RNA seq specific quality control tools like `RNAseqQC` or `RSeQC`. Finally, validate a subset of variants with Sanger sequencing or targeted DNA sequencing to estimate false positive rates.

## Common Mistakes

One common mistake is using default DNA thresholds for filtering. RNA seq data has higher error rates, so apply more stringent values. Another is failing to running `SplitNCigarReads`. Without it, HaplotypeCaller will produce many false positive splice junction variants. A third mistake is using deduplication without consideration of expression levels. This can remove true variants from highly expressed genes. Conversely, not marking duplicates at all for DNA seq pipelines can lead to inflated confidence.

Many users also forget to use the correct reference. The reference genome must match the one used during alignment. Using a different reference causes errors. Some users attempt VQSR on RNA seq data with low variant counts. VQSR requires at least 10,000 known variant sites to train the model. If your target capture or gene expression panel covers only a few hundred genes, skip VQSR. Finally, do not assume that all variants called are real. RNA seq variant calling has known limits.

## Limits and Uncertainty

RNA seq variant calling has several intrinsic limitations. First, coverage is biased by expression. Variants in lowly expressed genes may be missed entirely. Second, mapping errors around splice junctions create false positives. Even with `SplitNCigarReads`, complex splicing patterns can produce incorrect alignments. Third, allele specific expression can mask heterozygous variants. If one allele is not expressed, the variant appears homozygous. Newer methods using deep learning show promise but are not yet standard [A deep learning based RNA seq germline variant caller](https://pubmed.ncbi.nlm.nih.gov/37416509/). Long read RNA seq introduces different challenges, such as higher error rates that require transformation of alignment files [Transformation of alignment files improves performance of variant callers for long read RNA sequencing data](https://pubmed.ncbi.nlm.nih.gov/37095564/).

Multi omics studies combine RNA seq variant calls with other data types, but integration requires careful normalization [Multi omics data of gastric cancer cell lines](https://pubmed.ncbi.nlm.nih.gov/37081404/). Intersecting RNA seq variants with DNA seq variants improves confidence but does not guarantee accuracy. You should always report the applied filters and quality metrics. The field continues to evolve, and no single pipeline is universally correct.

## Frequently Asked Questions

**1. Can I use GATK Mutect2 for RNA seq somatic variant calling?**
Yes, Mutect2 has an RNA seq mode activated with the ` --rna` flag. However, be aware that expression artifacts increase false positive rates. Validate somatic calls using matched DNA or orthogonal methods.

**2. Should I remove duplicate reads for RNA seq variant calling?**
It depends. If your library has high PCR duplication and low complexity, removal may help. But generally, for variant discovery, keep duplicates to retain biological signal. Use metrics to decide.

**3. What reference genome should I use for RNA seq variant calling?**
Use the same reference genome as your aligner. Recent studies show that the T2T CHM13 reference improves variant detection in repetitive regions, but standard GRCh38 remains widely used.

**4. How do I handle strand specific data?**
GATK can infer strand information from the BAM file. Ensure your aligner produced a proper BAM with the correct strand flag. The `SplitNCigarReads` tool preserves strand information.

## Related Clinical & Scientific Guides

* [Observational vs. Experimental Studies: How to Tell Them Apart](/blog/guides/observational-vs-experimental-studies-how-to-tell-them-apart)
* [Astrocyte Single Cell Rna Seq](/blog/guides/astrocyte-single-cell-rna-seq)
* [Structural Genes](/blog/guides/structural-genes)


## References and Further Reading

- [NCBI Bookshelf](https://www.ncbi.nlm.nih.gov/books/) General background on sequencing and variant analysis.
- [EMBL EBI Training](https://www.ebi.ac.uk/training/) Courses on RNA seq analysis and GATK usage.
- [Galaxy Training Network](https://training.galaxyproject.org/) Hands on workflows for RNA seq variant calling.
- [Bioconductor](https://bioconductor.org/) R packages for VCF manipulation and quality control.
- [A comprehensive review of variant calling tools for RNA seq: Challenges and advances](https://pubmed.ncbi.nlm.nih.gov/42081940/) Overview of tool strengths and weaknesses.
- [Improved somatic mutation profiling in small cell lung cancer using RNA sequencing and the T2T CHM13 genome reference](https://pubmed.ncbi.nlm.nih.gov/41485300/) Reference impact on variant detection.
- [A deep learning based RNA seq germline variant caller](https://pubmed.ncbi.nlm.nih.gov/37416509/) Novel approach for improving accuracy.
- [Transformation of alignment files improves performance of variant callers for long read RNA sequencing data](https://pubmed.ncbi.nlm.nih.gov/37095564/) Adaptation for long reads.
- [Assessment of single nucleotide variant discovery protocols in RNA seq data from human cells exposed to mycotoxins](https://pubmed.ncbi.nlm.nih.gov/36016515/) Practical filtering benchmarks.

## Related Articles

- [Protein Synthesis](/blog/guides/protein-synthesis)
- [Incomplete Dominance Gene](/blog/guides/incomplete-dominance-gene)
- [Cell Membrane Function Biology](/blog/guides/cell-membrane-function-biology)
- [Dna Structure](/blog/guides/dna-structure)
- [Protein Structure](/blog/guides/protein-structure)