Building an RNA-seq Variant Calling Pipeline

At the end of this you will have three files for one sample: a filtered VCF of SNPs and short indels restricted to transcribed regions, a per-site allele count table you can use for allele-specific expression, and a concordance report comparing the RNA calls against whole-genome calls from the same person. You need the following inputs and software.

• paired-end RNA-seq FASTQs (we assume poly-A selected, 2×150 bp, 40–80 M read pairs)
• the GRCh38 no-ALT analysis set FASTA
• a GENCODE GTF matching that assembly
• dbSNP and the Mills/1000G indel VCFs from the GATK resource bundle
• STAR ≥ 2.7.10
• GATK ≥ 4.5
• samtools
• bcftools

Budget 64 GB RAM for the STAR index and about 250 GB of disk per sample if you keep intermediates. Runtime on 16 cores is roughly two hours end to end, most of it in HaplotypeCaller.

Before writing any of it, be clear about what RNA-seq variant calling can and cannot tell you. You only see variants in bases that were transcribed in the tissue you sampled, at a depth set by expression level, after the read has survived splice-aware alignment. A variant absent from the RNA VCF is not evidence of a reference genotype. RNA variant calling is a complement to whole-genome sequencing and does not substitute for it. Clinical interpretation of any variant belongs with a clinician and a CLIA-validated assay ¹.

1. Pick the reference, and pick it once

Use the GRCh38 no-ALT analysis set (GCA_000001405.15_GRCh38_no_alt_analysis_set.fna). STAR does not handle ALT contigs sensibly, and leaving them in scatters reads that belong in the primary assembly across near-identical haplotypes, which shows up later as heterozygous calls with absurd allele balance across the MHC. Decoys (hs38d1) are fine to keep and help absorb reads from sequences missing in the primary assembly.

Use one GENCODE release for the index, the quantification, and the annotation, and write the version into your file names. Pipeline differences of exactly this kind are the main driver of non-reproducible call sets. They include reference build, annotation version and tool version. The community effort to define functionally equivalent pipelines exists because two reasonable-looking pipelines on the same reads can disagree on thousands of sites ².

samtools faidx GRCh38_no_alt.fa
gatk CreateSequenceDictionary -R GRCh38_no_alt.fa

2. Build the STAR index

--sjdbOverhang should be read length minus 1. For 2×150 bp reads that is 149.

STAR --runMode genomeGenerate \
  --runThreadN 16 \
  --genomeDir star_idx_gencode44_149 \
  --genomeFastaFiles GRCh38_no_alt.fa \
  --sjdbGTFfile gencode.v44.primary_assembly.annotation.gtf \
  --sjdbOverhang 149

This needs about 32–40 GB of RAM and writes roughly 30 GB. If you are on a memory-constrained box, --genomeSAsparseD 2 halves the suffix array footprint at a modest cost in alignment speed.

3. Align with two-pass mode

Two-pass alignment matters here more than it does for expression quantification. In the first pass STAR collects novel splice junctions, inserts them into the index, and re-aligns. Reads that span a junction not in GENCODE get aligned as spliced reads instead of being soft-clipped or forced into a mismatch-heavy ungapped alignment. Those soft-clipped ends are a large source of false SNP clusters near exon boundaries.

STAR --runThreadN 16 \
  --genomeDir star_idx_gencode44_149 \
  --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
  --readFilesCommand zcat \
  --twopassMode Basic \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMattrRGline ID:S1 SM:S1 PL:ILLUMINA LB:lib1 PU:unit1 \
  --outSAMattributes NH HI AS nM MD NM \
  --outFilterMultimapNmax 20 \
  --outFilterMismatchNmax 999 \
  --outFilterMismatchNoverReadLmax 0.04 \
  --alignIntronMin 20 \
  --alignIntronMax 1000000 \
  --alignMatesGapMax 1000000 \
  --alignSJoverhangMin 8 \
  --alignSJDBoverhangMin 1 \
  --quantMode GeneCounts \
  --outFileNamePrefix S1.

Set the read group in the aligner. GATK refuses to run without one, and fixing it afterwards with AddOrReplaceReadGroups is a wasted pass over the BAM. --quantMode GeneCounts gives you a gene-level count table for free, which you will want in step 8 to interpret coverage.

4. Mark duplicates

For expression quantification, deduplicating is usually wrong: highly expressed transcripts generate genuine duplicate start positions. For variant calling it is the opposite. A PCR jackpot on one allele inflates allele depth and produces a confident call from a single original molecule. Mark them, do not remove them, and let GATK’s read filters ignore them.

gatk MarkDuplicates \
  -I S1.Aligned.sortedByCoord.out.bam \
  -O S1.md.bam \
  -M S1.dup_metrics.txt \
  --CREATE_INDEX true

Expect duplicate rates of 20–60% on poly-A libraries. That number is not a library failure by itself. If you have UMIs, use fgbio GroupReadsByUmi plus UmiAwareMarkDuplicatesWithMateCigar instead, which is strictly better here because it separates PCR duplicates from coincidental fragment-start collisions.

5. SplitNCigarReads

HaplotypeCaller does not understand N operators in the CIGAR string. SplitNCigarReads splits each spliced read at the N into separate alignment records and hard-clips the overhangs into the intron.

gatk SplitNCigarReads \
  -R GRCh38_no_alt.fa \
  -I S1.md.bam \
  -O S1.split.bam

Two behaviors to know. First, GATK4 reassigns MAPQ 255 (STAR’s “uniquely mapped”) to 60 by default, so a separate ReassignOneMappingQuality step is unnecessary. You can turn this off with --skip-mapping-quality-transform if you have already remapped the quality yourself. Second, the tool is single-threaded and slow, typically 30–60 minutes on a full transcriptome BAM. It is the natural place to scatter by contig if you are parallelizing.

6. Base quality score recalibration

BQSR works on RNA-seq, with the caveat that the covariate model is estimated from a read population whose error profile includes reverse-transcriptase errors and A-to-I edited bases. Run it after splitting, not before, and always with known sites supplied.

gatk BaseRecalibrator \
  -R GRCh38_no_alt.fa -I S1.split.bam \
  --known-sites Homo_sapiens_assembly38.dbsnp138.vcf \
  --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  --known-sites Homo_sapiens_assembly38.known_indels.vcf.gz \
  -O S1.recal.table

gatk ApplyBQSR \
  -R GRCh38_no_alt.fa -I S1.split.bam \
  --bqsr-recal-file S1.recal.table \
  -O S1.recal.bam

If your reads come off a two-color instrument with binned quality scores (NovaSeq, four bins), the gain from BQSR is small. We still run it because the downstream QD and FS filters assume calibrated qualities.

7. Call variants

gatk HaplotypeCaller \
  -R GRCh38_no_alt.fa \
  -I S1.recal.bam \
  -O S1.raw.vcf.gz \
  --dont-use-soft-clipped-bases true \
  --standard-min-confidence-threshold-for-calling 20 \
  --dbsnp Homo_sapiens_assembly38.dbsnp138.vcf

--dont-use-soft-clipped-bases is the flag that matters most. Without it, the clipped tails left over from reads crossing unannotated junctions get assembled into the local haplotype graph and generate dense false variant clusters at exon ends. Turning it on costs you a little sensitivity for indels near junctions and buys a large reduction in false positives.

Restricting to exons with -L gencode_exons_pad100.interval_list cuts runtime by roughly half. The tradeoff is that you lose calls in retained introns, 3’ UTR extensions, and unannotated transcripts. We call genome-wide and filter afterwards, because intronic coverage is informative about what the library did.

Joint calling across samples is not supported for RNA-seq in the GATK workflow. The GVCF-then-GenotypeGVCFs path assumes a uniform, genotype-independent depth model that RNA violates. In RNA a site is zero-depth because the gene is off. Zero depth does not mean the genome is absent. Groups have assembled joint workflows for RNA anyway, and they can work for population-scale cohorts where you accept that “no call” and “reference” are being conflated. For a single person, call per sample.

8. Filter

VQSR needs thousands of samples and a depth distribution it can model. Neither holds. Use hard filters, including the cluster filter, which removes any variant in a window of 35 bases containing three or more variants. This targets exactly the exon-boundary artifact pattern described above.

gatk VariantFiltration \
  -R GRCh38_no_alt.fa \
  -V S1.raw.vcf.gz \
  -O S1.filt.vcf.gz \
  --window 35 --cluster 3 \
  --filter-name "FS" --filter-expression "FS > 30.0" \
  --filter-name "QD" --filter-expression "QD < 2.0"

Then apply depth and blacklist filters with bcftools. A depth floor of 10 supporting reads is the practical minimum for calling a heterozygote in RNA. Below that, allelic dropout from sampling alone will make many true hets look homozygous. Clinical pipelines use similar per-site depth floors and require orthogonal confirmation for anything reportable ¹.

bcftools view -f PASS S1.filt.vcf.gz \
  | bcftools filter -i 'FORMAT/DP >= 10' \
  | bcftools view -T ^encode_blacklist_v2.bed \
  | bcftools view -T ^rediportal_editing_sites.bed \
  -Oz -o S1.pass.vcf.gz
bcftools index -t S1.pass.vcf.gz

Subtracting known A-to-I editing sites (REDIportal) is specific to RNA. ADAR converts adenosine to inosine, which reverse transcribes as guanosine. An edited site therefore appears as a het A>G call with allele fraction anywhere from 5% to 90%. The call is present in the transcriptome and absent from the genome. In a typical whole-blood or LCL sample this accounts for a few thousand apparent variants, heavily enriched in Alu elements in 3’ UTRs. If you see a pile of A>G and T>C (the minus-strand equivalent) calls with low allele fraction, that is editing. Those calls do not reflect genotype.

Also expect calls in immunoglobulin and T-cell receptor loci in any sample containing lymphocytes. Somatic V(D)J recombination means the DNA in those cells genuinely differs from the reference, and the aligned reads look like dense variant clusters. Exclude IGH, IGK, IGL, TRA and TRB regions unless that is the thing you are studying.

9. Check against your genome

This is the step most tutorials skip and the one that tells you whether your pipeline works. If you have WGS on the same person, compute concordance restricted to regions with adequate RNA coverage.

samtools depth -a -Q 10 S1.recal.bam \
  | awk '$3>=10 {print $1"\t"$2-1"\t"$2}' \
  | bedtools merge -i - > S1.rna_cov10.bed

gatk Concordance \
  -eval S1.pass.vcf.gz \
  -truth wgs_sample.vcf.gz \
  -L S1.rna_cov10.bed \
  -S S1.concordance.tsv

Healthy numbers for a good poly-A library on a well-expressed gene set: SNP precision above 0.95 and recall above 0.90 within the covered interval list, with indel recall meaningfully worse (0.7–0.85 is common). Indels near homopolymers and in low-complexity 3’ UTRs are the weak spot. If SNP precision is below 0.90, look first at whether your cluster filter is on and whether you excluded editing sites.

The covered interval list itself is the useful artifact. It tells you which genes this sample can speak to at all. In a whole-blood RNA sample, 10,000 to 13,000 genes clear 10× at some position. Everything else is dark.

10. Allele-specific expression

Once you have het sites (ideally from WGS, so the genotype is independent of the RNA), count reads per allele.

gatk ASEReadCounter \
  -R GRCh38_no_alt.fa \
  -I S1.recal.bam \
  -V wgs_hets.vcf.gz \
  -O S1.ase.tsv \
  --min-depth-of-non-filtered-base 10 \
  --min-mapping-quality 10 \
  --min-base-quality 20

Raw counts carry reference bias: a read carrying the alternate allele has one more mismatch to the reference and is slightly less likely to align. Unbiased, the reference fraction across all het sites sits near 0.50. In practice you will see 0.52–0.54. Correct it with WASP, which STAR implements natively. Pass your het VCF with --varVCFfile, add vW to --outSAMattributes, and set --waspOutputMode SAMtag. STAR then swaps each read’s alleles, re-aligns, and tags reads that fail to map back to the same place. Keep only vW:i:1 reads.

Allele-specific expression is the part of this pipeline that WGS cannot give you. A heterozygous variant in a splice region or a promoter is a hypothesis about regulation until you observe the two alleles at unequal ratios in the transcript, and combining sequence with functional readouts like this is how rare regulatory effects get resolved ³. Interpretation still takes care: imprinted genes show monoallelic expression by design, and clonal expansion in blood produces apparent ASE with no cis-regulatory cause at all.

Common problems

Everything in the VCF is A>G. RNA editing. Subtract REDIportal and look at the strand distribution. This is the single most common surprise for people coming from DNA pipelines.

Dense variant clusters at exon boundaries. Either --dont-use-soft-clipped-bases is off, --twopassMode Basic is off, or the cluster filter (--window 35 --cluster 3) is not applied. Check all three.

HaplotypeCaller errors on N in CIGAR. You skipped SplitNCigarReads or pointed the caller at the wrong BAM.

Recall looks terrible genome-wide. You are probably computing it over the whole genome instead of over covered regions. Restrict to the rna_cov10.bed interval list. Recall over regions the assay cannot see is meaningless.

Indel calls near polyA stretches in 3’ UTRs. Internal priming and reverse-transcriptase slippage both put artifacts here. Filter with a low-complexity-region mask, or deprioritize indels in homopolymers longer than 6 bases.

Two runs of the same data disagree. Check tool versions, reference build, and GENCODE release before anything else. Differences at this level dominate, which is why pipeline harmonization is treated as a first-class problem rather than a detail ².

A variant looks important. Stop the analysis and take it to a clinician with a genetics service. RNA-derived calls have not been validated for clinical use, allele fractions in RNA do not mean what they mean in DNA, and confirmation on an orthogonal DNA assay is the baseline standard ¹.

Footnotes

1. Daniel C. Koboldt. Best practices for variant calling in clinical sequencing. Genome Medicine, 2020. https://doi.org/10.1186/s13073-020-00791-w ↩ ↩² ↩³

2. Allison Regier, Yossi Farjoun, David E. Larson, et al. Functional equivalence of genome sequencing analysis pipelines enables harmonized variant calling across human genetics projects. Nature Communications, 2018. https://doi.org/10.1038/s41467-018-06159-4 ↩ ↩²

3. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩