How to Align and Quantify Your Own RNA-seq Data

Oak
rna-seq alignment salmon star quantification
A studio-lit iridescent salmon-like specimen in a glass column of fluid, its tail ending in a glowing five-pointed radial fin, against black.

By the end of this you will have two complementary outputs from the same FASTQ files: a transcript-level abundance table from Salmon (estimated counts, effective lengths, and TPM for roughly 250,000 GENCODE transcripts), and a coordinate-sorted genome BAM from STAR with splice junction counts and per-gene read counts. The first is what you want for expression questions, the second is what you want for anything positional: splice junction usage, allele-specific expression against your whole-genome variants, coverage inspection of a single locus. You need paired-end FASTQs (we assume 2x100 bp or 2x150 bp polyA-selected or rRNA-depleted libraries), a machine with at least 64 GB of RAM and 500 GB of free disk if you intend to run STAR, and a working conda or container setup. Everything below is command-line, single-sample, and reproducible on a laptop only if you skip the genome alignment.

A word on what “alignment” means here, since the term covers two quite different operations. Classical sequence alignment finds the arrangement of matches, mismatches, and gaps between two sequences that maximizes a scoring function, which Smith-Waterman solves exactly by dynamic programming in time proportional to the product of the two lengths. Nothing that reads a whole sequencing run does that exhaustively. Practical read mappers use seed-and-extend: find short exact or near-exact matches with an index, then extend only those candidate locations with a banded dynamic programming step, which is the same heuristic BLAST introduced and which remains the target of hardware acceleration work decades later 1. RNA-seq adds one twist, which is that a read from a mature mRNA can span an exon-exon boundary and therefore does not align contiguously to the genome. Either you align to the genome and allow large deletions at splice junctions (STAR, HISAT2), or you align to transcript sequences where the splicing has already been applied (Salmon, kallisto, Bowtie2 against a transcriptome FASTA).

1. Pick the reference and build the transcript FASTA

Use GENCODE, not RefSeq, and not Ensembl unless you have a reason. GENCODE and Ensembl annotations are the same underlying gene set, but GENCODE ships the exact file combination you want in one place: a primary assembly genome FASTA, a comprehensive GTF on that assembly, and transcript sequences. Avoid the full analysis set with alt contigs for quantification, because alt haplotypes create spurious multi-mapping. Download the primary assembly FASTA and the comprehensive annotation GTF for a recent release and keep the release number in every filename you produce. Changing annotation release between timepoints in a longitudinal series will shift TPMs for affected genes by more than real biology does.

REL=47
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/GRCh38.primary_assembly.genome.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/gencode.v${REL}.annotation.gtf.gz
gunzip GRCh38.primary_assembly.genome.fa.gz gencode.v${REL}.annotation.gtf.gz

# Derive transcript sequences from the GTF so they are guaranteed consistent with it
gffread -w transcripts.v${REL}.fa -g GRCh38.primary_assembly.genome.fa gencode.v${REL}.annotation.gtf

Deriving transcripts with gffread rather than downloading the prebuilt transcript FASTA guarantees that every transcript ID in your quantification exists in the GTF you use downstream, which removes an entire class of tximport failures. Keep the comprehensive annotation rather than the basic set. The extra retained-intron and processed-transcript entries soak up reads that would otherwise be misassigned to the protein-coding isoform, and if you only care about gene-level output you will sum them away anyway.

2. Trim and QC the reads

Run fastp once per sample. It does adapter detection, quality trimming, polyG trimming for NovaSeq two-color chemistry, and produces a JSON report you can parse programmatically.

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o trim_R1.fastq.gz  -O trim_R2.fastq.gz \
  --detect_adapter_for_pe \
  --trim_poly_g \
  --cut_tail --cut_tail_mean_quality 20 \
  --length_required 36 \
  --thread 8 \
  --json fastp.json --html fastp.html

Do not quality-trim aggressively. Soft-clipping in the aligner handles low-quality read ends better than truncation does, and over-trimming shortens reads until they multi-map. The numbers to read out of fastp.json are the duplication rate, the insert size peak, and read1_after_filtering.total_reads. For a personal profile we want at least 30 million read pairs passing filter for gene-level work and closer to 60-100 million if you intend to look at isoform ratios, because isoform-level estimates for genes with many similar transcripts are variance-limited, not bias-limited.

3. Build a decoy-aware Salmon index

This is the step people skip and then regret. Salmon quantifies against transcript sequences, so a read originating from an intron or an unannotated intergenic region has nowhere correct to go and will be assigned to whichever transcript it resembles most. Adding the whole genome as “decoy” sequence gives those reads a home, and Salmon’s selective alignment scores the transcriptomic mapping against the genomic one before accepting it. The cost is index build time and memory.

grep "^>" GRCh38.primary_assembly.genome.fa | cut -d " " -f1 | sed 's/>//g' > decoys.txt
cat transcripts.v47.fa GRCh38.primary_assembly.genome.fa > gentrome.fa

salmon index \
  -t gentrome.fa \
  -d decoys.txt \
  -i salmon_idx_v47 \
  -k 31 \
  --gencode \
  -p 16

--gencode strips the pipe-delimited GENCODE FASTA headers down to the bare ENST identifier, which is what you want for joining to a tx2gene table. -k 31 is the default and correct choice for 100 bp or longer reads. Drop to -k 23 only for reads shorter than 50 bp, and expect mapping rate to rise while specificity falls. Budget around 30 GB of RAM and 40 minutes on 16 cores for the index build, and keep the index directory, because rebuilding it per sample is a waste.

4. Quantify with Salmon

salmon quant \
  -i salmon_idx_v47 \
  -l A \
  -1 trim_R1.fastq.gz -2 trim_R2.fastq.gz \
  --seqBias --gcBias \
  --numBootstraps 30 \
  --validateMappings \
  -p 16 \
  -o quant_sample01

-l A lets Salmon infer library type, and you should check what it inferred in quant_sample01/lib_format_counts.json. A standard dUTP stranded protocol reports ISR. If it reports IU or a near-even split of strand-compatible and strand-incompatible fragments, your library is unstranded, which is fine but costs you the ability to resolve overlapping antisense genes. --seqBias and --gcBias fit models for random-hexamer priming bias at fragment ends and for GC-dependent fragment abundance, and both are cheap enough to always enable. --numBootstraps 30 produces the inferential replicates you need later to decide whether an isoform-level difference between two of your own timepoints exceeds the technical uncertainty of the estimate itself, which for tightly related isoforms is often large.

The two outputs that matter are quant.sf (TransactionID, Length, EffectiveLength, TPM, NumReads) and logs/salmon_quant.log, where the line to find is the mapping rate. For a clean polyA human library against a decoy-aware index, expect 75-90 percent. Below about 65 percent something is wrong and step 8 tells you where to look.

5. Align to the genome with STAR

Salmon never produces alignments, so anything positional requires a real spliced aligner. We use STAR for personal data because it emits a transcriptome-coordinate BAM, per-gene counts, and a splice junction table in a single pass, and because its two-pass mode re-indexes on junctions discovered in the first pass, which improves sensitivity for novel junctions in your own sample.

STAR --runMode genomeGenerate \
  --genomeDir star_idx_v47 \
  --genomeFastaFiles GRCh38.primary_assembly.genome.fa \
  --sjdbGTFfile gencode.v47.annotation.gtf \
  --sjdbOverhang 100 \
  --runThreadN 16

STAR --genomeDir star_idx_v47 \
  --readFilesIn trim_R1.fastq.gz trim_R2.fastq.gz \
  --readFilesCommand zcat \
  --twopassMode Basic \
  --outSAMtype BAM SortedByCoordinate \
  --quantMode TranscriptomeSAM GeneCounts \
  --outFilterMultimapNmax 20 \
  --outFilterMismatchNoverReadLmax 0.04 \
  --alignSJoverhangMin 8 \
  --alignSJDBoverhangMin 1 \
  --alignIntronMin 20 --alignIntronMax 1000000 \
  --outSAMattributes NH HI AS nM NM MD \
  --runThreadN 16 \
  --outFileNamePrefix sample01.

Set --sjdbOverhang to read length minus one (99 for 2x100, 149 for 2x150). The index needs roughly 32 GB of RAM to build and about 30 GB resident to align, which is the practical reason to run this on a rented machine rather than a laptop. Include MD and NM in the SAM attributes if you plan to do variant-aware work. The ReadsPerGene.out.tab file gives you unstranded, forward, and reverse counts in columns 2 through 4; pick the column matching the library type Salmon inferred, and if you pick wrong you will see roughly 10-20 percent of the expected counts, which is a satisfyingly obvious failure.

For gene-level differential work across your own timepoints we prefer the Salmon estimates over the STAR gene counts, because Salmon distributes multi-mapping reads probabilistically across isoforms rather than discarding ambiguous assignments, and because effective length correction matters when library insert sizes drift between batches. Use the STAR BAM for inspection and positional questions.

6. Collapse to genes with tximport, and understand what TPM is

Transcript-level estimates are noisy and gene-level estimates are stable, so most interpretation happens at the gene level. tximport does the aggregation correctly, accounting for the fact that transcripts within a gene have different effective lengths.

library(tximport); library(GenomicFeatures); library(readr)

txdb  <- makeTxDbFromGFF("gencode.v47.annotation.gtf")
k     <- keys(txdb, keytype = "TXNAME")
tx2gene <- AnnotationDbi::select(txdb, k, "GENEID", "TXNAME")

files <- c(sample01 = "quant_sample01/quant.sf")
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

countsFromAbundance = "lengthScaledTPM" gives you gene-level counts that are already corrected for changes in average transcript length between samples, which is the form to hand to limma-voom. If you plan to use DESeq2 instead, pass countsFromAbundance = "no" and let DESeq2 consume txi$counts with txi$length as a normalization offset via DESeqDataSetFromTximport.

On the recurring question of TPM versus RPKM: both divide read counts by transcript length and by sequencing depth, but in opposite orders. RPKM divides by total mapped reads first and then by kilobases of transcript, so the per-sample sum of RPKM values depends on the length distribution of what happened to be expressed, and the sum is different in every sample. TPM divides by length first, then rescales so that the values sum to exactly one million, which makes each TPM a fraction of the transcript molecules in the sample. That fixed denominator is what makes TPM comparable across samples. Neither metric handles composition: if hemoglobin transcripts occupy 60 percent of your blood library, every other gene’s TPM is suppressed by roughly that factor even at constant absolute abundance. For comparing your own timepoints, quantify with TPM for intuition and run the statistics on counts with median-of-ratios or TMM normalization.

7. QC the alignment before you interpret anything

Run these three and record the numbers in a table you keep across timepoints, because the absolute values matter less than their drift.

samtools flagstat sample01.Aligned.sortedByCoord.out.bam

picard CollectRnaSeqMetrics \
  I=sample01.Aligned.sortedByCoord.out.bam \
  O=sample01.rnaseq_metrics.txt \
  REF_FLAT=refFlat.v47.txt \
  RIBOSOMAL_INTERVALS=rRNA.v47.interval_list \
  STRAND_SPECIFICITY=SECOND_READ_TRANSCRIPTION_STRAND

The thresholds we treat as pass/fail for human blood or tissue polyA libraries: uniquely mapped reads above 85 percent of input, PCT_MRNA_BASES above 0.70, PCT_RIBOSOMAL_BASES below 0.05, PCT_INTRONIC_BASES below 0.25, and MEDIAN_5PRIME_TO_3PRIME_BIAS between 0.7 and 1.3. High intronic fraction with normal mapping rate usually means DNA contamination or incomplete splicing from a degraded sample. A 3’ bias above 1.5 means RNA degradation, which compresses long transcripts specifically and will produce a false pattern of “long genes down” that has fooled many people.

The check we consider mandatory for personal data is identity concordance with your own genome. Call genotypes at common exonic SNPs from the RNA BAM and compare against your WGS VCF. If the concordance is not above 95 percent at covered sites, you are looking at someone else’s transcriptome.

bcftools mpileup -f GRCh38.primary_assembly.genome.fa \
  -R common_exonic_snps.bed -Ou sample01.Aligned.sortedByCoord.out.bam \
| bcftools call -m -Oz -o rna_calls.vcf.gz
bcftools index rna_calls.vcf.gz
bcftools gtcheck -g wgs_genotypes.vcf.gz rna_calls.vcf.gz

As a second, cheaper sanity check, look at XIST and RPS4Y1 TPMs and confirm they match the expected pattern for your karyotype. Both checks take minutes and catch the failure mode that no downstream statistics can rescue.

8. Interpret in the context of your other layers, carefully

The reason to align your transcriptome rather than only your genome is that expression is the layer where a variant either does something or does not. A rare splice-region variant in your VCF is a hypothesis, and the STAR SJ.out.tab junction counts are the observation that either supports it or does not. Similarly, allele-specific expression at a heterozygous site is measurable only when you have both the phased genome and the RNA BAM in the same coordinate space, which is the practical argument for running STAR alongside Salmon. Expression state is partially recoverable even from other modalities: nucleosome footprints at transcription start sites in cell-free plasma DNA carry enough signal to infer which genes are expressed in the cells contributing that DNA 2. Integrating transcript abundance with proteomics and metabolite data is an active and unsettled area, and current approaches range from straightforward correlation across timepoints to learned joint embeddings 3 4.

What RNA expression from a blood draw cannot tell you is what is happening in a tissue you did not sample. Whole blood transcriptomes are dominated by the cell composition of that draw, so a shift in neutrophil fraction will move thousands of genes in a way that looks like a coordinated biological program and is not. Deconvolve, or at minimum record a differential count alongside the sample. And any expression pattern that you are tempted to read as evidence of disease is a question for a clinician working from validated clinical assays, because none of the above is a diagnostic test.

Common problems

Low Salmon mapping rate (under 65 percent) has four usual causes, in order of frequency. Adapter contamination that fastp missed, which shows up as a spike in the fastp overrepresented sequences; check for it before anything else. An index built from a different genome build than your annotation, which produces transcripts with wrong sequences; rebuild both from the same release. rRNA-heavy libraries, where most reads have no polyA transcript to map to, visible as high PCT_RIBOSOMAL_BASES. And non-human contamination, which you can test by mapping a subsample of 100,000 reads with --writeUnmappedNames and running the unmapped reads through a classifier.

Transcript-level estimates that swing wildly between replicates or timepoints are usually real uncertainty rather than a bug. Genes with many highly similar isoforms distribute reads almost arbitrarily among them, and the bootstrap replicates you generated in step 4 will show it. Aggregate to gene level, or restrict isoform claims to genes where the inferential variance is small.

A common annotation mismatch arises from GENCODE version suffixes. quant.sf contains ENST00000456328.2 while many downstream tables carry ENST00000456328. Strip suffixes consistently at one point in your pipeline, never at two, and never in a way that could collapse two distinct transcripts.

Finally, STAR out-of-memory kills during two-pass mode are common on machines with 32 GB. The genome plus junction database does not fit comfortably alongside sorting buffers. Either run on 64 GB, or set --outSAMtype BAM Unsorted and sort separately with samtools sort -m 2G -@ 8, which keeps peak memory predictable.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Kaiyong Zhao, Xiaowen Chu. G-BLASTN: accelerating nucleotide alignment by graphics processors. Bioinformatics, 2014. https://doi.org/10.1093/bioinformatics/btu047 ↩

  2. Peter Ulz, Gerhard G Thallinger, Martina Auer, et al. Inferring expressed genes by whole-genome sequencing of plasma DNA. Nature Genetics, 2016. https://doi.org/10.1038/ng.3648 ↩

  3. Heng-Rui Liu. AI-driven integration of multi-omics and multimodal data for precision medicine. Medical Data Mining, 2026. https://doi.org/10.53388/mdm202609001 ↩

  4. Amir Ebrahimi, Alireza Fotuhi Siahpirani, Hesam Montazeri. sCIN: a contrastive learning framework for single-cell multi-omics data integration. Briefings in Bioinformatics, 2025. https://doi.org/10.1093/bib/bbaf411 ↩