How to Perform Germline Variant Calling from RNA-Seq Data: Challenges and Best Practices
By Dr. Zubair Khalid, DVM, MS, PhD ·

Key Takeaways
- Germline variant calling from RNA-seq leverages transcriptome reads to identify inherited SNVs and indels, particularly useful when DNA samples are unavailable or for studying allele-specific expression.
- The workflow necessitates splice-aware alignment (e.g., STAR, HISAT2) and read splitting at exon junctions to accurately map reads spanning introns, preventing false variant calls at splice sites.
- Allele-specific expression (ASE) is a critical challenge, as differential expression of alleles can skew variant allele frequencies (VAFs) from the expected 50% for heterozygotes, necessitating cautious interpretation of VAFs.
- Post-alignment processing includes marking duplicates and base quality score recalibration (BQSR) to correct systematic errors, followed by variant calling with tools like GATK HaplotypeCaller in RNA-seq mode or deep-learning models like DeepVariant.
- Stringent filtering is paramount, incorporating expression-aware filters to remove variants in poorly expressed genes and applying quality metrics (e.g., QD, FS, MQ) or machine learning classifiers (e.g., VarRNA) to distinguish true variants from artifacts like RNA editing or library preparation errors.
- Orthogonal validation, such as Sanger sequencing or targeted genotyping, is essential for clinically relevant findings due to the inherent limitations of RNA-seq variant calling, including the inability to detect variants in unexpressed genes.
Germline variant calling from RNA sequencing (RNA-seq) data identifies inherited single nucleotide variants (SNVs) and small insertions or deletions (indels) by analyzing transcriptome reads instead of DNA sequencing reads. This approach serves researchers when DNA samples are unavailable, when the functional impact of variants through allele-specific expression is the research focus, or when existing RNA-seq datasets need repurposing for variant discovery. The core challenge is that RNA-seq data contains splicing events, variable transcript abundance, allele-specific expression, and mapping biases that do not affect DNA-based variant calling. A practical workflow combines splice-aware alignment, base quality score recalibration, variant calling with tools designed for transcriptome data, and stringent filtering to remove artifacts arising from RNA processing and expression dynamics. This article provides a concrete workflow for researchers and laboratory professionals who need to call germline variants from RNA-seq data and interpret the results responsibly.
Scope and Reader Context
This article addresses researchers, biology students, laboratory professionals, and life-science practitioners who have RNA-seq data and want to extract germline variant information from it. The primary use case is the analysis of bulk RNA-seq data from a single individual or a small cohort where germline variants are of interest. The workflow described here follows the Genome Analysis Toolkit (GATK) RNA-seq short variant calling pipeline, with attention to the modifications needed for transcriptome data. The article also covers alternative approaches, including deep-learning-based callers and machine-learning classifiers that distinguish germline variants from somatic mutations and artifacts.
The intended practical outcome is a reproducible variant calling workflow that produces a filtered variant call format (VCF) file with variants that can be interpreted in the context of gene expression. Readers should understand the limitations of RNA-seq variant calling, including the inability to detect variants in genes that are not expressed in the sequenced tissue, the confounding effect of allele-specific expression, and the need for orthogonal validation when clinical decisions depend on the results.
At a Glance
The table below summarizes the key decisions in an RNA-seq germline variant calling workflow, the recommended approach, and the rationale for each choice.
| Workflow Step | Recommended Approach | Rationale |
|---|---|---|
| Input data | Paired-end bulk RNA-seq reads in FASTQ format, preferably with strand-specific library preparation | Paired-end reads improve alignment across splice junctions, strand-specific libraries reduce false variant calls from antisense transcription |
| Alignment | Splice-aware aligner such as STAR or HISAT2, producing BAM files with read groups | Splice-aware alignment maps reads spanning exon-exon junctions correctly, which is essential for transcriptome data |
| Preprocessing | Mark duplicates, split reads at exon junctions, base quality score recalibration | Splitting reads at junctions prevents false variant calls at splice sites, recalibration corrects systematic base quality errors |
| Variant calling | GATK HaplotypeCaller in RNA-seq mode, or a deep-learning caller such as DeepVariant with an RNA-seq model | GATK provides a validated per-sample workflow, DeepVariant learns RNA-specific error patterns and can outperform traditional callers |
| Filtering | Expression-aware filters, variant quality score recalibration (VQSR) or hard filters, removal of variants in poorly expressed genes | Filtering removes artifacts from mapping errors, RNA editing, and low expression, expression filters reduce false positives |
| Validation | Compare against DNA-derived variants when available, use orthogonal methods such as Sanger sequencing or targeted genotyping | RNA-seq variant calling has inherent limitations, validation is required for clinically relevant findings |
Understanding the Biological and Technical Challenges
Splicing and Junction Reads
RNA-seq reads originate from mature messenger RNA, which means they span exon-exon junctions after introns have been removed. A read that crosses a splice junction will not align contiguously to the reference genome because the genomic distance between the two exons includes the intron. Splice-aware aligners handle this by allowing reads to be split across exons, but this introduces complexity in variant calling. If the aligner does not correctly identify the junction, reads may map incorrectly, producing false variant calls near splice sites. The GATK RNA-seq workflow addresses this by splitting reads at exon junctions before variant calling, which prevents the caller from interpreting junction-spanning reads as evidence for variants at the genomic breakpoints.
Allele-Specific Expression
Germline variants are present in the DNA of every cell, but their representation in RNA-seq data depends on the expression level of the gene and the relative expression of the two alleles. If one allele is expressed at a much higher level than the other, the variant allele frequency (VAF) in the RNA-seq data will deviate from the expected 50% heterozygote frequency. This phenomenon, called allele-specific expression, can cause a heterozygous germline variant to appear homozygous or can cause a variant to be missed entirely if the alternate allele is not expressed. A study of cancer transcriptomes demonstrated that variant allele frequencies in RNA-seq data can differ substantially from the corresponding DNA exome data, and this discrepancy is particularly pronounced in cancer-driving genes where the variant allele may be expressed at much higher levels than expected based on the exome sequencing data. Researchers must therefore interpret VAFs from RNA-seq data with caution and should not use VAF alone to distinguish germline from somatic variants.
Mapping Biases and Sequence Context
The transcriptome presents unique mapping challenges. Reads from highly expressed genes dominate the dataset, while reads from lowly expressed genes may be sparse or absent. Paralogous genes with high sequence similarity can cause reads to map to the wrong location, producing false variant calls. Additionally, RNA editing events, particularly A-to-I editing, can create apparent variants that are not present in the DNA. A deep-learning-based variant caller trained on RNA-seq data can learn to account for some of these error sources, but no method can completely eliminate them. The study that extended DeepVariant to RNA-seq data showed that the model outperformed existing approaches such as Platypus and GATK, and the authors examined how the model addresses RNA editing events and how additional thresholding can facilitate use in a production pipeline.
Library Preparation Artifacts
RNA-seq library preparation introduces artifacts that are not present in DNA sequencing. Reverse transcription can introduce errors, PCR amplification can create duplicate reads and chimeric molecules, and the fragmentation step can produce reads that do not represent the original transcript. Indels are particularly problematic because PCR-based library preparation can generate artifacts at homopolymer regions and other repetitive sequences. A study of somatic indel detection from tumor RNA-seq data noted that reliable identification of expressed indels is challenging due to artifacts generated in PCR-based library preparation. While that study focused on somatic variants, the same artifacts affect germline indel calling, and researchers should apply stringent filters to indel calls from RNA-seq data.
Core Principles of RNA-Seq Variant Calling
Variant Calling Is Secondary to Expression Analysis
RNA-seq experiments are typically designed for differential expression analysis, not variant calling. The library preparation, sequencing depth, and read length are optimized for quantifying transcript abundance, not for detecting variants. Researchers should recognize that variant calling from RNA-seq data is a secondary analysis that can provide useful information but cannot replace DNA-based variant calling for clinical or diagnostic purposes. The sensitivity of RNA-seq variant calling is limited by the expression level of each gene, and variants in genes that are not expressed in the sequenced tissue will be missed entirely.
The Workflow Must Be Splice-Aware
Every step of the workflow, from alignment to variant calling, must account for the spliced nature of RNA-seq reads. A standard DNA alignment tool will fail to map junction-spanning reads, producing gaps in coverage at exon boundaries and false variant calls at the edges of exons. The GATK RNA-seq short variant calling pipeline includes specific steps for handling spliced reads, and the pipeline has been validated for per-sample variant calling. A 2022 methods paper described how modern GATK commands from distinct workflows can be combined to call variants on RNA-seq samples, providing a detailed tutorial that starts with raw RNA-seq reads and ends with filtered variants, some of which were associated with bovine paratuberculosis.
Expression Context Is Essential for Interpretation
A variant call from RNA-seq data is only meaningful if the gene is expressed in the sequenced tissue. The depth of coverage at a variant site reflects the expression level of the gene, and low coverage sites produce low-confidence calls. Researchers should record the expression level of each gene and the depth of coverage at each variant site, and they should filter out variants in genes with insufficient expression. The interpretation of a variant should also consider whether the variant allele is expressed, because a germline variant that is not expressed in the tissue has no functional consequence in that tissue.
Recommended Workflow for Germline Variant Calling from RNA-Seq
Step 1: Input Data Preparation and Quality Control
The workflow begins with raw RNA-seq reads in FASTQ format. Paired-end reads are strongly recommended because they improve alignment accuracy across splice junctions and provide better evidence for variant calling. Before alignment, run quality control on the raw reads to check for adapter contamination, low-quality bases, and GC bias. The FastQC tool is commonly used for this purpose, and the Galaxy Training Network offers accessible tutorials for quality control and subsequent analysis steps.
Record the following information for each sample:
- Sequencing platform and read length
- Library preparation method, including whether the library is strand-specific
- Number of raw read pairs
- Mean quality scores per base position
- Adapter contamination levels
- GC content distribution
Step 2: Splice-Aware Alignment
Align the cleaned reads to the reference genome using a splice-aware aligner. STAR and HISAT2 are the most commonly used tools for this purpose. The alignment step produces a BAM file that must include read group information, which is required for downstream GATK tools. Read groups identify the sample, library, and lane, and they are essential for base quality score recalibration and for combining multiple libraries from the same sample.
Alignment parameters should be adjusted for the specific data type. For example, the --outSAMstrandField parameter in STAR can be set to produce the strand information needed for downstream analysis. The --outSAMattributes parameter should include the NH (number of reported alignments) and HI (hit index) attributes, which are used for filtering multi-mapping reads.
After alignment, assess mapping statistics:
- Percentage of reads mapped
- Percentage of reads mapped to multiple locations
- Percentage of reads mapped to intergenic regions
- Insert size distribution
- Coverage uniformity across genes
Step 3: Post-Alignment Processing
The BAM file requires several processing steps before variant calling. These steps are implemented in the GATK RNA-seq short variant calling pipeline and are described in the methods paper on combining GATK workflows for RNA-seq variant calling.
First, mark duplicate reads. Duplicates arise from PCR amplification during library preparation, and they should not be counted as independent evidence for a variant. The MarkDuplicates tool identifies reads that start at the same position and have the same orientation, and it flags them as duplicates.
Second, split reads at exon junctions. The SplitNCigarReads tool splits reads that span exon-exon junctions into separate segments, which prevents the variant caller from interpreting the junction as evidence for a deletion or other variant. This step is critical for RNA-seq data and is not needed for DNA-seq data.
Third, perform base quality score recalibration (BQSR). This step uses known variant sites to model the systematic errors in base quality scores and adjusts the scores accordingly. BQSR requires a set of known variants, which can be obtained from databases such as dbSNP, available through the NCBI Data Resources. The NCBI provides official descriptions of its databases, search systems, and sequence resources, including dbSNP and the Sequence Read Archive.
Step 4: Variant Calling
Call variants using GATK HaplotypeCaller with the RNA-seq mode enabled. The RNA-seq mode adjusts the caller's behavior to account for the spliced nature of the reads and the presence of junction-spanning reads. The command is:
gatk HaplotypeCaller \
-R reference.fasta \
-I sample.bam \
-O sample.g.vcf.gz \
-ERC GVCF \
--dont-use-soft-clipped-bases true
The --dont-use-soft-clipped-bases parameter prevents the caller from using soft-clipped bases, which are often artifacts of alignment across splice junctions.
The output is a genomic VCF (gVCF) file that contains the genotype likelihoods for every position in the genome. The gVCF format is used for joint genotyping across multiple samples, which improves the accuracy of variant calls by leveraging information from all samples in a cohort.
Alternative variant callers include deep-learning-based approaches. The DeepVariant RNA-seq model study demonstrated that a deep-learning-based variant caller extended to RNA-seq data produces highly accurate variant calls and outperforms existing approaches such as Platypus and GATK. The study examined factors that influence accuracy, how the model addresses RNA editing events, and how additional thresholding can facilitate use in a production pipeline. Researchers who have large cohorts or who need maximum accuracy may consider using DeepVariant with the RNA-seq model.
Step 5: Joint Genotyping and Filtering
If multiple samples are available, perform joint genotyping using GenotypeGVCFs. Joint genotyping combines the per-sample gVCF files and produces a multi-sample VCF file. The 2022 methods paper on GATK RNA-seq workflows described how modern GATK commands from distinct workflows can be combined to call variants on RNA-seq samples, including the possibility of performing joint genotyping analysis, which is not part of the fully validated per-sample GATK RNA-seq workflow.
After genotyping, filter the variants to remove artifacts. The GATK RNA-seq workflow recommends hard filtering based on the following criteria:
- Quality by depth (QD) less than 2.0
- Fisher strand bias (FS) greater than 60.0
- Mapping quality (MQ) less than 40.0
- Strand odds ratio (SOR) greater than 3.0
- ReadPosRankSumTest less than -8.0
These thresholds are starting points, and researchers should adjust them based on the characteristics of their data. The Galaxy Training Network provides tutorials on variant filtering that can help researchers understand the effect of different thresholds.
Step 6: Expression-Aware Filtering
The final filtering step is specific to RNA-seq data. For each variant, record the depth of coverage and the expression level of the gene. Filter out variants in genes with low expression, because low coverage produces unreliable genotype calls. A common approach is to require a minimum depth of 10 reads at the variant site and a minimum gene expression level, such as a minimum of 1 transcript per million (TPM).
Also consider allele-specific expression. If a variant is heterozygous in the DNA but appears homozygous in the RNA-seq data, this may indicate allele-specific expression instead of a genotyping error. The 2025 study of cancer transcriptomes showed that some variants classified by the VarRNA method exhibit variant allele frequencies distinct from the corresponding DNA exome data, and this phenomenon is prevalent in cancer-driving genes. Researchers should flag variants with extreme VAFs for further investigation.
Alternative Approaches and Tools
Machine-Learning Classifiers for Variant Classification
The VarRNA method study uses RNA-seq data to classify single nucleotide variants and indels from tumor transcriptomes as germline, somatic, or artifact. The method uses two XGBoost machine-learning models trained and validated on pediatric cancer samples with paired tumor and normal DNA exome sequencing data as ground truth. VarRNA identifies 50% of the variants detected by exome sequencing and detects unique RNA variants absent in paired tumor and normal DNA exome data. This approach is particularly useful for cancer research, where distinguishing germline from somatic variants is essential.
Specialized Indel Callers
Indel calling from RNA-seq data is more challenging than SNV calling due to artifacts from PCR-based library preparation. The RNAIndel tool study predicts somatic, germline, and artifact indels from tumor RNA-seq data using features derived from indel sequence context and biological effect in a machine-learning framework. The study reported that RNAIndel robustly predicts 88-100% of somatic indels in five diverse test datasets of pediatric and adult cancers, outperforming the current best-practice for RNA-seq variant calling which had 57% sensitivity but with 14 times more false positives. Researchers who need to call indels from RNA-seq data should consider using RNAIndel or similar specialized tools.
Workflow Management Systems
Reproducibility is a major concern in bioinformatics, and workflow management systems can help ensure that analyses are reproducible. The nf-core documentation describes the standards for pipeline development and usage. A 2023 study on CWL pipelines assembled automated pipelines for RNA-seq, ChIP-seq, and germline variant calling analyses in Common Workflow Language (CWL), demonstrating that CWL-implemented workflows achieved high accuracy in reproducing previously published results and detecting germline SNP and small indel variants. The workflows are publicly available on GitHub and are suitable for the analysis of short-read data.
The Bioconductor project provides official documentation for R packages and workflows for genomic analysis, including packages for variant annotation and visualization. The EMBL-EBI Training offers bioinformatics learning pathways and practical analysis education that can help researchers develop the skills needed for variant calling. The Carpentries lessons provide foundational computing, data, shell, Git, and programming training that is useful for researchers who need to develop the computational skills required for bioinformatics analysis.
Practical Implementation Steps
Setting Up the Analysis Environment
Before starting the analysis, set up a reproducible computing environment. Use a workflow management system such as Nextflow with nf-core pipelines, or use CWL with containerization to ensure that software versions are consistent across analyses. The nf-core documentation describes the standards for pipeline usage and configuration, and the 2023 CWL study demonstrated that containerization overcomes issues of software incompatibility and laborious configuration requirements.
Install the required software:
- A splice-aware aligner such as STAR or HISAT2
- GATK version 4.x
- Picard tools for BAM processing
- A variant caller such as HaplotypeCaller or DeepVariant
- A variant filtering tool such as GATK VariantFiltration or bcftools
- An annotation tool such as ANNOVAR or SnpEff
Running the Workflow
The following steps outline the practical implementation of the workflow. Each step includes the key commands and the expected outputs.
- Quality control: Run FastQC on the raw FASTQ files and review the reports for adapter contamination, low-quality bases, and GC bias.
- Trimming: If adapter contamination is present, trim the reads using Trimmomatic or cutadapt. Record the trimming parameters and the number of reads removed.
- Alignment: Run STAR or HISAT2 to align the reads to the reference genome. Use the appropriate parameters for the library type and read length.
- Post-alignment processing: Use GATK MarkDuplicates, SplitNCigarReads, and BaseRecalibrator to process the BAM file.
- Variant calling: Run HaplotypeCaller with the RNA-seq mode to produce a gVCF file.
- Joint genotyping: If multiple samples are available, run GenotypeGVCFs to produce a multi-sample VCF file.
- Filtering: Apply hard filters to remove low-quality variants. Record the filtering thresholds and the number of variants removed at each step.
- Annotation: Annotate the filtered variants with gene names, variant effects, and population frequencies using ANNOVAR or SnpEff.
- Expression analysis: For each variant, record the expression level of the gene and the depth of coverage at the variant site.
- Validation: If possible, validate a subset of variants using an orthogonal method such as Sanger sequencing or targeted genotyping.
Recording and Documentation
Maintain detailed records of the analysis for reproducibility. The Carpentries lessons emphasize the importance of good data management practices, including version control with Git and documentation of analysis steps. Record the following information:
- Software versions for all tools
- Reference genome version and source
- Alignment parameters
- Filtering thresholds
- Number of variants at each step
- Any deviations from the standard workflow
Records and Measurements
Key Metrics to Track
The following metrics should be tracked throughout the analysis to assess the quality of the variant calls:
| Metric | Description | Interpretation |
|---|---|---|
| Mapping rate | Percentage of reads that map to the reference genome | Low mapping rates indicate contamination or alignment problems |
| Duplicate rate | Percentage of reads marked as duplicates | High duplicate rates reduce effective coverage |
| Mean coverage | Average depth of coverage across expressed genes | Low coverage reduces variant calling sensitivity |
| Transition/transversion ratio | Ratio of transition to transversion variants | Deviations from expected ratios may indicate artifacts |
| Variant density | Number of variants per kilobase | High variant density may indicate mapping errors |
| Heterozygous/homozygous ratio | Ratio of heterozygous to homozygous variants | Extreme ratios may indicate allele-specific expression or contamination |
| Number of variants in exons | Variants that fall within coding regions | These are the most interpretable variants |
Quality Control Thresholds
The following thresholds are commonly used in RNA-seq variant calling workflows. These are starting points, and researchers should adjust them based on their specific data and research questions.
- Minimum mapping quality: 20
- Minimum base quality: 20
- Minimum depth: 10 reads at the variant site
- Minimum gene expression: 1 TPM
- Maximum FS: 60.0
- Minimum QD: 2.0
Common Failure Patterns and Troubleshooting
Failure Pattern 1: Low Variant Calling Sensitivity
If the variant caller identifies very few variants, the most likely cause is low coverage in expressed genes. RNA-seq coverage is highly variable across genes, and genes with low expression will have insufficient depth for variant calling. Check the coverage distribution and consider whether the sequencing depth is adequate for the research question. If the goal is to identify variants in specific genes, consider targeted RNA-seq or DNA sequencing instead.
Failure Pattern 2: High False Positive Rate
A high false positive rate is often caused by mapping errors, RNA editing, or library preparation artifacts. Check the transition/transversion ratio and the number of variants in homopolymer regions. If the false positive rate is high, apply more stringent filters, particularly for strand bias and mapping quality. Consider using a deep-learning-based caller such as DeepVariant, which can learn to account for RNA-specific error patterns.
Failure Pattern 3: Allele-Specific Expression Confounds Genotype Calls
If a variant appears homozygous in the RNA-seq data but is expected to be heterozygous based on population frequencies, this may indicate allele-specific expression. Check the expression level of the gene and the depth of coverage at the variant site. If the gene is expressed and the depth is adequate, the variant may be subject to allele-specific expression, and the genotype call should be interpreted with caution.
Failure Pattern 4: Indel Calling Artifacts
Indel calls from RNA-seq data are frequently artifacts of PCR-based library preparation. If the number of indel calls is unusually high, check the sequence context of the indels. Indels in homopolymer regions are particularly suspect. Consider using a specialized tool such as RNAIndel, which is designed to distinguish true indels from artifacts.
Failure Pattern 5: Batch Effects Across Samples
If variant calls differ systematically across batches of samples, this may indicate batch effects in library preparation or sequencing. Check the duplicate rate, mapping rate, and coverage distribution across batches. If batch effects are present, include batch information in the analysis and consider using a workflow management system to standardize the analysis across batches.
Limitations and Interpretation Boundaries
RNA-Seq Cannot Detect Variants in Unexpressed Genes
The most fundamental limitation of RNA-seq variant calling is that variants can only be detected in genes that are expressed in the sequenced tissue. A germline variant in a gene that is not expressed in the tissue will be missed entirely. Researchers should not use RNA-seq variant calling to rule out the presence of a variant in a specific gene, particularly if the gene is known to have tissue-specific expression.
Variant Allele Frequencies Are Not Reliable for Genotype Inference
The VAF of a variant in RNA-seq data reflects the relative expression of the two alleles, not the genotype. A heterozygous variant can have a VAF close to 0% or 100% if one allele is preferentially expressed. Researchers should not use VAF thresholds to distinguish heterozygous from homozygous variants in RNA-seq data.
RNA Editing Creates Apparent Variants
RNA editing, particularly A-to-I editing, creates apparent variants that are not present in the DNA. These variants are biologically interesting but should not be interpreted as germline variants. The deep-learning-based DeepVariant RNA-seq model study examined how the model addresses RNA editing events, but no method can completely distinguish RNA editing from true variants without DNA data.
Validation Is Required for Clinically Relevant Findings
RNA-seq variant calling is not appropriate for clinical or diagnostic purposes without orthogonal validation. If a variant has clinical implications, validate it using DNA sequencing, Sanger sequencing, or targeted genotyping. The NCBI Data Resources provide access to databases and tools that can help researchers interpret variants, including ClinVar and dbSNP.
Safety and Regulatory Context
Data Privacy and Consent
RNA-seq data may contain sensitive genetic information. Researchers must ensure that the data were collected with appropriate informed consent and that the analysis complies with applicable privacy regulations. The NCBI Data Resources provide guidance on data sharing and privacy, and researchers should follow the data access policies of their institutions and funding agencies.
Reporting of Incidental Findings
RNA-seq variant calling may reveal variants that have clinical implications for the research participant. Researchers should have a plan for handling incidental findings, including whether and how to report them to participants. This plan should be developed in consultation with the institutional review board or ethics committee.
Professional Escalation Criteria
Researchers should escalate to a clinical genetics professional or genetic counselor in the following situations:
- A variant with known clinical significance is identified in a gene associated with a serious disease
- A variant is identified in a gene with pharmacogenomic implications
- The research participant requests information about their genetic variants
- The analysis reveals a potential sample mix-up or contamination
Decision Framework for Selecting an RNA-Seq Variant Calling Strategy
Choosing the right variant calling approach for RNA-seq data depends on the research question, cohort size, available computing resources, and the tolerance for false positives versus false negatives. This section provides a practical decision framework that researchers can apply before starting the analysis, along with a record system for documenting the rationale behind each choice.
Step 1: Define the Primary Research Objective
The first decision determines whether RNA-seq variant calling is appropriate at all. Record the primary objective in the analysis plan before any computation begins.
| Objective | Recommended Approach | Rationale |
|---|---|---|
| Discover variants in a single sample with no DNA data | GATK HaplotypeCaller RNA-seq mode with hard filtering | Validated per-sample workflow with clear filtering guidance |
| Call variants across a cohort for association studies | GATK with joint genotyping or DeepVariant RNA-seq model | Joint genotyping improves accuracy at low-coverage sites, DeepVariant learns RNA-specific error patterns |
| Distinguish germline from somatic variants in cancer data | VarRNA or RNAIndel for indels | Machine-learning classifiers use paired tumor and normal DNA exome data as ground truth |
| Maximize sensitivity for rare variant discovery | DeepVariant RNA-seq model with additional thresholding | The DeepVariant RNA-seq study demonstrated superior accuracy over Platypus and GATK |
| Prioritize specificity to minimize false positives | GATK with stringent hard filters and expression-aware filtering | Conservative thresholds reduce artifacts from splicing and RNA editing |
Step 2: Assess Cohort Size and Sample Availability
The number of samples and the availability of matched DNA data change the optimal strategy.
For a single sample with no matched DNA, use the per-sample GATK RNA-seq workflow. This is the fully validated approach described in the 2022 methods paper on GATK RNA-seq variant calling. The workflow produces a gVCF file that can be analyzed independently or combined with other samples later.
For a cohort of 10 or more samples, consider joint genotyping. The same 2022 methods paper described how modern GATK commands can be combined to perform joint genotyping on RNA-seq samples, which is not part of the fully validated per-sample workflow. Joint genotyping leverages information across samples to improve genotype calls at sites where individual samples have low coverage due to variable gene expression.
For cancer studies with matched tumor and normal DNA exome data, use the VarRNA approach. The VarRNA study trained two XGBoost machine-learning models using paired tumor and normal DNA exome sequencing data as ground truth. This approach distinguishes germline variants from somatic mutations and artifacts, which is essential for cancer research.
Step 3: Evaluate the Target Variant Types
SNVs and indels require different strategies because indels are more prone to artifacts from PCR-based library preparation.
For SNV-focused studies, GATK HaplotypeCaller or DeepVariant with the RNA-seq model are appropriate. The DeepVariant RNA-seq study demonstrated that the deep-learning model produces highly accurate variant calls and outperforms existing approaches such as Platypus and GATK.
For indel-focused studies, use a specialized tool. The RNAIndel study reported that RNAIndel robustly predicts 88-100% of somatic indels in five diverse test datasets, outperforming the current best-practice for RNA-seq variant calling which had 57% sensitivity but with 14 times more false positives. While RNAIndel was developed for somatic indels, the same artifact challenges apply to germline indel calling from RNA-seq data.
Step 4: Consider Computing Resources and Reproducibility Requirements
The choice between running a manual workflow and using a workflow management system depends on the scale of the analysis and the need for reproducibility.
For small analyses with fewer than 10 samples, a manual workflow with documented commands is acceptable. Record all software versions, parameters, and reference genome versions in a laboratory notebook or version-controlled document.
For larger analyses or when reproducibility across institutions is required, use a workflow management system. The nf-core documentation describes community pipeline standards for usage and configuration. A 2023 study on CWL pipelines demonstrated that CWL-implemented workflows achieved high accuracy in reproducing previously published results and detecting germline SNP and small indel variants. The workflows are publicly available on GitHub and use containerization to overcome software incompatibility issues.
Step 5: Document the Decision Rationale
Record the following information for each analysis to ensure that the chosen strategy can be justified and reproduced:
- Primary research objective and the specific question the variant calls will answer
- Number of samples and whether matched DNA data are available
- Target variant types (SNVs, indels, or both)
- Expected expression levels of the genes of interest in the sequenced tissue
- Computing resources available, including memory, storage, and CPU hours
- Reproducibility requirements, such as institutional or journal mandates
- Validation plan, including which variants will be confirmed by orthogonal methods
Comparison of Variant Calling Strategies
The table below compares the main approaches for germline variant calling from RNA-seq data.
| Strategy | Strengths | Limitations | Best Use Case |
|---|---|---|---|
| GATK HaplotypeCaller RNA-seq mode | Validated per-sample workflow, clear filtering guidance, widely documented | Per-sample workflow does not include joint genotyping by default | Single samples or small cohorts where DNA data are unavailable |
| GATK with joint genotyping | Improved accuracy at low-coverage sites, leverages cohort information | Requires combining commands from distinct workflows, not fully validated | Cohorts of 10 or more samples |
| DeepVariant RNA-seq model | Outperforms Platypus and GATK, learns RNA-specific error patterns | Requires deep-learning infrastructure, model training details may be complex | Large cohorts where maximum accuracy is needed |
| VarRNA | Distinguishes germline, somatic, and artifact variants | Requires paired tumor and normal DNA exome data for training | Cancer studies with matched DNA data |
| RNAIndel | High sensitivity for indels, distinguishes true indels from artifacts | Focused on somatic indels, may require adaptation for germline calling | Studies where indel detection is a priority |
Implementation Checklist
Before starting the analysis, complete the following checklist and record the decisions in the analysis plan:
- Confirm that the genes of interest are expressed in the sequenced tissue. Check publicly available expression databases or preliminary expression analysis results.
- Determine whether the research question requires germline variant calls or whether somatic variants are also of interest.
- Verify that the library preparation method and sequencing platform are compatible with the chosen variant calling approach.
- Check that the reference genome version is consistent across all samples and tools.
- Confirm that read group information is present in the BAM files or will be added during alignment.
- Identify the validation method and select the variants that will be validated before the analysis begins.
- Document the software versions and parameters for every tool in the workflow.
- Establish the filtering thresholds and record the rationale for each threshold.
- Define the escalation criteria for variants with potential clinical significance.
- Store all raw data, intermediate files, and final VCF files in a secure location with appropriate access controls.
Common Decision Errors
The following errors are frequently observed when researchers select a variant calling strategy for RNA-seq data.
Error 1: Using a DNA-seq variant calling workflow without modification. DNA-seq workflows do not account for spliced reads, and they will produce false variant calls at exon-exon junctions. Always use a splice-aware aligner and split reads at exon junctions before variant calling.
Error 2: Assuming that VAF reflects genotype. The 2025 VarRNA study demonstrated that variant allele frequencies in RNA-seq data can differ substantially from the corresponding DNA exome data, particularly in cancer-driving genes. Do not use VAF thresholds to determine zygosity.
Error 3: Ignoring expression level when interpreting variant calls. A variant call in a gene with very low expression is unreliable because the depth of coverage is insufficient. Filter out variants in genes below the expression threshold and record the expression level for every variant.
Error 4: Applying the same filtering thresholds to SNVs and indels. Indels are more prone to artifacts from PCR-based library preparation, as noted in the RNAIndel study. Apply more stringent filters to indel calls or use a specialized tool such as RNAIndel.
Error 5: Failing to document the decision rationale. Without documentation, the analysis cannot be reproduced or defended in peer review. Record the rationale for every decision, including the choice of variant caller, filtering thresholds, and validation strategy.
Professional Escalation Criteria
Escalate to a bioinformatics specialist or clinical genetics professional in the following situations:
- The analysis reveals a variant with known clinical significance in a gene associated with a serious disease
- The variant calls will be used for clinical or diagnostic purposes without orthogonal validation
- The analysis produces an unexpectedly high or low number of variants, suggesting a systematic error
- The research participant requests information about their genetic variants
- The analysis reveals a potential sample mix-up or contamination that cannot be resolved with the available data
The EMBL-EBI Training and Galaxy Training Network offer learning pathways that can help researchers develop the skills needed to make informed decisions about variant calling strategies. The Carpentries lessons provide foundational computing and data management training that supports reproducible analysis practices.
Frequently Asked Questions
What is the difference between germline and somatic variant calling from RNA-seq data?
Germline variant calling identifies variants that are present in the DNA of every cell and are inherited from the parents. Somatic variant calling identifies variants that arise during the lifetime of the individual and are present only in specific cells or tissues. In RNA-seq data, germline variants are expected to be present in all expressed genes, while somatic variants are present only in the affected tissue. The distinction is complicated by allele-specific expression, which can cause germline variants to appear somatic and vice versa. Machine-learning classifiers such as VarRNA can help distinguish germline from somatic variants by using features derived from the variant sequence context and expression patterns.
Why is splice-aware alignment necessary for RNA-seq variant calling?
RNA-seq reads originate from mature mRNA, which has had introns removed. A read that spans an exon-exon junction will not align contiguously to the reference genome because the genomic distance between the two exons includes the intron. A splice-aware aligner allows reads to be split across exons, mapping each segment to the correct exon. Without splice-aware alignment, junction-spanning reads would either fail to map or map incorrectly, producing gaps in coverage and false variant calls near splice sites.
Can RNA-seq variant calling replace DNA sequencing for clinical applications?
RNA-seq variant calling cannot replace DNA sequencing for clinical applications. RNA-seq data only provides information about expressed genes, and the depth of coverage is determined by expression level instead of by experimental design. Variants in unexpressed genes are missed entirely, and allele-specific expression can confound genotype calls. For clinical applications, DNA sequencing is the gold standard, and RNA-seq variant calls should be validated with an orthogonal method before any clinical decision is made.
How does allele-specific expression affect variant calling from RNA-seq data?
Allele-specific expression occurs when one allele of a gene is expressed at a higher level than the other allele. This can cause the variant allele frequency in RNA-seq data to deviate from the expected 50% for a heterozygous variant. If the alternate allele is not expressed, the variant will appear homozygous for the reference allele. If the alternate allele is overexpressed, the variant will appear homozygous for the alternate allele. Researchers should not use VAF alone to determine the genotype of a variant from RNA-seq data.
What are the best tools for calling indels from RNA-seq data?
Indel calling from RNA-seq data is more challenging than SNV calling due to artifacts from PCR-based library preparation. The RNAIndel tool is specifically designed for this purpose and uses a machine-learning framework to predict somatic, germline, and artifact indels from tumor RNA-seq data. The study that introduced RNAIndel reported that it outperformed the current best-practice for RNA-seq variant calling, which had lower sensitivity and many more false positives. Researchers who need to call indels from RNA-seq data should consider using RNAIndel or a similar specialized tool.
How should I validate variants identified from RNA-seq data?
Variants identified from RNA-seq data should be validated using an orthogonal method. The most common validation approaches are Sanger sequencing, targeted genotyping, or DNA sequencing of the same sample. When selecting variants for validation, prioritize variants that have clinical implications, variants in genes of interest, and variants with unusual VAFs that may indicate allele-specific expression. The validation results should be recorded and used to assess the accuracy of the RNA-seq variant calling workflow.
What is the role of joint genotyping in RNA-seq variant calling?
Joint genotyping combines information from multiple samples to improve the accuracy of variant calls. The standard GATK RNA-seq workflow is a per-sample workflow that does not include joint genotyping, but a 2022 methods paper described how modern GATK commands can be combined to perform joint genotyping on RNA-seq samples. Joint genotyping is particularly useful when analyzing a cohort of samples, because it leverages information from all samples to improve genotype calls at sites where individual samples have low coverage.
How do I distinguish RNA editing events from true germline variants?
RNA editing events, particularly A-to-I editing, create apparent variants that are not present in the DNA. These events are biologically interesting but should not be interpreted as germline variants. The deep-learning-based DeepVariant RNA-seq model study examined how the model addresses RNA editing events, but no method can completely distinguish RNA editing from true variants without DNA data. Researchers can also use databases of known RNA editing sites to filter out likely editing events.
Related Bioinformatics Guides
- Alternative Splicing Analysis from RNA-Seq Data
- RNA-Seq Databases: Accessing and Using Public RNA-Seq Data
- RNA-Seq Data Analysis in Galaxy: A User-Friendly Platform
- RNA-Seq Data Analysis Workflow: From Raw Reads to Insights
- RNA Sequencing Data Analysis: From Raw Reads to Differential Expression
Related Clinical & Scientific Guides
- A Practical Guide to Detecting Antimicrobial Resistance Genes in Shotgun Metagenomic Data
- Computational Immunology: Modeling the Immune System
- How to Set Hard Filters for Germline Variant Calling: A Practical Guide to GATK Best Practices
References and Further Reading
- NCBI Data Resources. National Center for Biotechnology Information.
- EMBL-EBI Training. European Bioinformatics Institute.
- Bioconductor. Bioconductor Project.
- Galaxy Training Network. Galaxy Project.
- nf-core Documentation. nf-core.
- The Carpentries Lessons. The Carpentries.
- Variant Calling from RNA-seq Data Using the GATK Joint Genotyping Workflow.. Methods in molecular biology (Clifton, N.J.), 2022.
- Variant calling from RNA-Seq data reveals allele-specific differential expression of pathogenic cancer variants.. Communications medicine, 2025.
- A deep-learning-based RNA-seq germline variant caller.. Bioinformatics advances, 2023.
- Software pipelines for RNA-Seq, ChIP-Seq and germline variant calling analyses in common workflow language (CWL).. Frontiers in bioinformatics, 2023.
- RNAIndel: discovering somatic coding indels from tumor RNA-Seq data.. Bioinformatics (Oxford, England), 2020.
This article is educational and does not replace validated analysis plans, institutional policy, clinical interpretation, or specialist review.