Gatk Rna Seq Variant Calling

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. Detailed training on RNA seq analysis is available from EMBL EBI Training.

At a Glance

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. 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. 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.

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.

Step 5: Annotate and export

Use SnpEff or ANNOVAR to annotate variant effects. For downstream analysis in R, use the VariantAnnotation package from Bioconductor. 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 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. 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.

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. 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.

References and Further Reading

• NCBI Bookshelf General background on sequencing and variant analysis.
• EMBL EBI Training Courses on RNA seq analysis and GATK usage.
• Galaxy Training Network Hands on workflows for RNA seq variant calling.
• Bioconductor R packages for VCF manipulation and quality control.
• A comprehensive review of variant calling tools for RNA seq: Challenges and advances Overview of tool strengths and weaknesses.
• Improved somatic mutation profiling in small cell lung cancer using RNA sequencing and the T2T CHM13 genome reference Reference impact on variant detection.
• A deep learning based RNA seq germline variant caller Novel approach for improving accuracy.
• Transformation of alignment files improves performance of variant callers for long read RNA sequencing data Adaptation for long reads.
• Assessment of single nucleotide variant discovery protocols in RNA seq data from human cells exposed to mycotoxins Practical filtering benchmarks.

Related Articles

• Protein Synthesis
• Incomplete Dominance Gene
• Cell Membrane Function Biology
• Dna Structure
• Protein Structure