How to Map RNA-seq Reads to the Human Genome

By the end of this guide you will have two complementary outputs from the same FASTQ files: a coordinate-sorted BAM aligned to GRCh38 that you can view in IGV, inspect base by base, and use for splicing or allele-specific work, and a transcript-level quantification table with effective lengths that collapses cleanly to gene-level counts for comparison across samples. You need the raw reads (paired-end FASTQ, gzipped) and a machine with at least 64 GB of RAM and 8 or more cores. You also need roughly 200 GB of free disk per sample-set, and a Linux environment with conda or containers. Everything below assumes human poly-A selected or rRNA-depleted bulk RNA-seq at 2×100 bp or 2×150 bp. If you have single-cell data the logic is similar but the tooling differs, and the analytical challenges shift toward sparsity and cell-type deconvolution rather than alignment itself ¹.

1. Decide which of the two mapping problems you have

There are two distinct tasks that both get called “mapping.” The first is spliced alignment to the genome: placing each read at genomic coordinates, allowing it to span introns, so you can look at splice junctions, novel transcripts, variants in expressed sequence, and coverage over any region you care about. The second is quantification: assigning reads probabilistically to annotated transcripts to estimate abundance, without ever committing a read to a single genomic position.

We run both, always. Genome alignment with STAR gives you a BAM that remains useful for years and answers questions you have not thought of yet. Selective alignment with Salmon gives better transcript-level estimates, corrects for fragment GC and positional bias, and finishes in a few minutes per sample. The tradeoff is disk and time: the STAR pass costs roughly 30 to 60 minutes and 30 GB of RAM per sample against Salmon’s 5 minutes and 16 GB. Running both costs you one extra hour and buys you a cross-check. When the two disagree on a gene by more than a factor of two, that gene has a multi-mapping or annotation problem, and you want to know that before you interpret it.

2. Build the reference from a no-alt analysis set

Reference choice determines more of your final numbers than any aligner parameter. Use the GRCh38 no-alt analysis set, which excludes the alternate haplotype contigs. Alt contigs are near-duplicates of primary chromosomes, and an aligner that is not alt-aware will scatter reads across them, deflating expression on the primary assembly for genes in the MHC and other polymorphic regions. Pair it with GENCODE annotation, which is the same annotation Ensembl ships and is more complete than RefSeq for non-coding transcripts.

mkdir -p ref && cd ref
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/GRCh38.primary_assembly.genome.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/gencode.v47.primary_assembly.annotation.gtf.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/gencode.v47.transcripts.fa.gz
gunzip GRCh38.primary_assembly.genome.fa.gz gencode.v47.primary_assembly.annotation.gtf.gz

Pin the release number and record it next to your results. GENCODE adds and retires transcripts at every release, and a gene-level count computed against v39 is not directly comparable to one computed against v47. For a longitudinal profile, where the whole point is comparing your own samples across months or years, re-quantifying old FASTQs against the current annotation is cheaper than reasoning about annotation drift.

3. Trim and inspect the reads before you align anything

Run fastp once per sample. It does adapter detection, quality trimming, poly-G correction for two-color chemistry (NovaSeq and NextSeq call dark cycles as G, which produces spurious high-quality poly-G tails), and emits a JSON report you can parse later.

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o sample_R1.trim.fastq.gz -O sample_R2.trim.fastq.gz \
  --detect_adapter_for_pe \
  --trim_poly_g \
  --length_required 36 \
  --thread 8 \
  --json sample.fastp.json --html sample.fastp.html

Resist heavier trimming. Aggressive quality trimming (say, -q 30) shortens reads and reduces the fraction that span splice junctions unambiguously. It generally lowers mapping rates rather than raising them. Modern spliced aligners soft-clip low-quality ends on their own. What you are looking for in the fastp report is the read length distribution, the adapter fraction (above roughly 20% suggests short inserts and possible degradation), and the duplication estimate.

4. Build the STAR index with the right sjdbOverhang

STAR builds a suffix array over the genome with splice junctions from the GTF baked in. The one parameter that matters is --sjdbOverhang, which should be your read length minus one. For 2×100 bp reads that is 99, for 2×150 bp it is 149. The value sets how much sequence is retained on each side of an annotated junction in the index, and setting it too low costs you junction-spanning reads.

STAR --runMode genomeGenerate \
  --runThreadN 16 \
  --genomeDir star_idx_gencode47_sjdb99 \
  --genomeFastaFiles ref/GRCh38.primary_assembly.genome.fa \
  --sjdbGTFfile ref/gencode.v47.primary_assembly.annotation.gtf \
  --sjdbOverhang 99 \
  --limitGenomeGenerateRAM 60000000000

Expect around 30 GB of peak RAM and 40 to 90 minutes on 16 threads. The resulting index is roughly 30 GB on disk. If you are memory-constrained, --genomeSAsparseD 2 halves the suffix array footprint at the cost of slower alignment, which is a reasonable trade on a laptop-class machine but not on a server. Name the index directory after the annotation release and the overhang value, because six months from now you will not remember which one you used.

5. Align with STAR in two-pass mode

Two-pass alignment runs the whole cohort once, collects the splice junctions discovered across all samples, and re-aligns with those junctions inserted into the index. This meaningfully improves detection of reads spanning novel or lowly-used junctions, which is exactly where single-pass alignment is weakest.

STAR \
  --genomeDir star_idx_gencode47_sjdb99 \
  --readFilesIn sample_R1.trim.fastq.gz sample_R2.trim.fastq.gz \
  --readFilesCommand zcat \
  --runThreadN 16 \
  --twopassMode Basic \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMattributes NH HI AS nM NM MD \
  --quantMode TranscriptomeSAM GeneCounts \
  --outFilterMultimapNmax 20 \
  --alignSJoverhangMin 8 \
  --alignSJDBoverhangMin 1 \
  --outFilterMismatchNoverReadLmax 0.04 \
  --alignIntronMin 20 \
  --alignIntronMax 1000000 \
  --alignMatesGapMax 1000000 \
  --outFileNamePrefix out/sample_

Those filter values are close to the ENCODE long-RNA settings and are a sensible default for human data. Setting --quantMode GeneCounts writes ReadsPerGene.out.tab with four columns. They are gene ID, counts for unstranded, counts for forward-stranded and counts for reverse-stranded. Do not guess which column to use, determine it empirically in step 7. TranscriptomeSAM emits a BAM in transcript coordinates that RSEM or Salmon’s alignment mode can consume, which is useful as a third cross-check.

Index the BAM afterward with samtools index -@ 8 out/sample_Aligned.sortedByCoord.out.bam so you can load it in IGV.

6. Quantify with a decoy-aware Salmon index

Salmon’s selective alignment scores each read against candidate transcripts and against decoy sequence drawn from the genome, so reads originating in introns or unannotated regions are not forced onto the nearest transcript. Without decoys, intronic and intergenic reads inflate the abundance of whichever transcript they happen to resemble, and the effect is worst for genes with retained introns.

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

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

salmon quant -i salmon_idx_v47 -l A \
  -1 sample_R1.trim.fastq.gz -2 sample_R2.trim.fastq.gz \
  -p 16 --validateMappings --gcBias --seqBias --posBias \
  --numBootstraps 30 \
  -o quant/sample

Keep -k 31 unless your reads are shorter than about 75 bp, in which case drop to 23 or 25. Adding --gcBias corrects fragment-level GC composition bias, which varies between library preps and is one of the larger sources of apparent expression difference between batches. Setting --numBootstraps 30 gives you inferential replicates, which matter for transcript-level estimates where two isoforms share most of their sequence and the quantifier genuinely cannot tell them apart. The --gencode flag strips the pipe-delimited GENCODE FASTA headers down to the transcript ID.

7. Determine strandedness and check alignment quality

Never assume strandedness from the kit name in the order confirmation. Determine it from the data with RSeQC, using a BED12 of your annotation.

infer_experiment.py -r gencode.v47.bed12 -i out/sample_Aligned.sortedByCoord.out.bam -s 400000

A typical dUTP protocol (Illumina Stranded mRNA, NEBNext Ultra II Directional) returns roughly 0.90 or higher for 1+-,1-+,2++,2--. That corresponds to reverse-stranded, which is column 4 of the STAR ReadsPerGene file, -s 2 in featureCounts and ISR in Salmon. An unstranded library splits near 0.5/0.5. Salmon’s -l A infers this automatically, and you should check that its inference in lib_format_counts.json agrees with RSeQC. Using the wrong column typically drops your gene counts by an order of magnitude, which is obvious. Alternatively it silently swaps antisense signal into sense counts, which is not obvious.

Then run the rest of the QC and aggregate it:

samtools flagstat -@ 8 out/sample_Aligned.sortedByCoord.out.bam
read_distribution.py -r gencode.v47.bed12 -i out/sample_Aligned.sortedByCoord.out.bam
geneBody_coverage.py -r housekeeping_hg38.bed -i out/sample_Aligned.sortedByCoord.out.bam -o sample
multiqc . -o multiqc_report

For a good human poly-A library, expect a uniquely mapped rate above 85%, %reads mapped to multiple loci under 10%, exonic fraction above 70%, and a gene body coverage curve that is roughly flat from the 10th to the 90th percentile of transcript length. A uniquely mapped rate in the 60s usually means contamination or the wrong reference. A gene body curve that rises steeply toward the 3′ end means degraded input, and the resulting counts are biased against long transcripts in a way no downstream normalization fully repairs.

8. Collapse to gene level and put the numbers on a usable scale

Bring the Salmon output into R with tximport, which handles the effective-length correction that makes transcript estimates summable to genes.

library(tximport); library(DESeq2)
tx2gene <- read.delim("tx2gene_v47.tsv", header = FALSE,
                      col.names = c("TXNAME", "GENEID"))
txi <- tximport(file.path("quant", samples, "quant.sf"),
                type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ timepoint)
dds <- dds[rowSums(counts(dds) >= 10) >= 3, ]
vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = "timepoint")

Use lengthScaledTPM rather than raw TPM when you plan any count-based modeling, because it preserves the count scale that negative binomial models expect while carrying the length correction. For a longitudinal single-person series, the first plot to make is a principal component projection over time. In our experience PC1 in repeated whole-blood sampling is frequently dominated by cell composition (neutrophil fraction in particular) rather than anything regulatory, so pair the expression data with the differential from the same blood draw before you read anything into a shift. Expression is one layer among several, and its interpretation improves considerably when it is anchored to genome, protein, and metabolite measurements from the same individual ²³.

9. Use your own genome to reduce reference bias

If you have whole-genome sequencing on the same person, you can remove a systematic artifact from allele-level analysis. Reads carrying the alternate allele at a heterozygous site accumulate one extra mismatch against the reference and map slightly worse, which fabricates apparent allelic imbalance. STAR implements WASP-style filtering directly:

STAR --genomeDir star_idx_gencode47_sjdb99 \
  --readFilesIn sample_R1.trim.fastq.gz sample_R2.trim.fastq.gz \
  --readFilesCommand zcat --runThreadN 16 \
  --varVCFfile sample.het.vcf \
  --waspOutputMode SAMtag --outSAMattributes NH HI AS nM vW vA vG \
  --outSAMtype BAM SortedByCoordinate \
  --outFileNamePrefix out/wasp_sample_

Keep reads where the vW tag equals 1. Each read is re-mapped with its alleles swapped, and those that move are discarded. This matters when you are asking whether a specific variant you carry changes expression of its own gene, which is one of the few RNA questions that a single individual can address well. Combining sequence, expression, and other molecular layers within one person is the direction most integrative work has moved, and the computational methods for handling that volume of data are now the limiting step rather than the measurement ⁴².

Common problems

Globin dominates whole blood. If your sample came from a PAXgene or Tempus tube without globin depletion, HBB, HBA1, and HBA2 can consume more than half of your reads, leaving far less depth for everything else. Check the fraction early: grep -E "ENSG00000244734|ENSG00000206172|ENSG00000188536" ReadsPerGene.out.tab. You can computationally exclude those genes before normalization, but you cannot recover the sequencing depth they consumed, so specify globin depletion at collection.

Residual rRNA eats depth silently. Poly-A selection usually leaves rRNA under 5% of reads. If RSeQC or a quick alignment against an rRNA FASTA shows 20% or more, the selection failed or the input was degraded, and your effective library size is much smaller than the read count suggests.

Mapping rate collapses from the wrong reference. A 40% unique mapping rate with high unmapped-too-short usually means either adapter read-through that trimming did not catch, or reads from a different organism or a different assembly version. Take 10,000 unmapped reads and BLAST a handful before you tune aligner parameters.

Multi-mapping hides paralogs and pseudogenes. Gene families with high sequence identity, along with processed pseudogenes that mirror their parents, generate reads the aligner cannot place. STAR’s GeneCounts discards multi-mappers entirely, while Salmon distributes them by expectation-maximization. When the two methods disagree sharply on a gene, this is nearly always why, and the Salmon estimate with bootstrap variance tells you how much confidence to place in it.

Batch effects masquerade as biology. Library prep date, flow cell, and input RNA integrity all shift global expression. In a longitudinal series, randomize or at least record these, and include them as covariates. If all your 2024 samples were prepped together and all your 2025 samples separately, no statistical method will cleanly separate the year from the batch.

Finally, an expression change in your own data is a measurement. It is not a finding about your health. Transcript abundance fluctuates with time of day, recent exercise, sleep, infection, and the exact handling of the tube. If something in the data prompts a clinical question, that question belongs with a physician working from validated clinical assays. Research-grade RNA-seq is not the basis for answering it.

Footnotes

1. Alev Baysoy, Zhiliang Bai, Rahul Satija, et al. The technological landscape and applications of single-cell multi-omics. Nature Reviews Molecular Cell Biology, 2023. https://doi.org/10.1038/s41580-023-00615-w ↩

2. C. Nelson Hayes, Hikaru Nakahara, Atsushi Ono, et al. From Omics to Multi-Omics: A Review of Advantages and Tradeoffs. Genes, 2024. https://doi.org/10.3390/genes15121551 ↩ ↩²

3. Biswapriya B Misra, Carl Langefeld, Michael Olivier, et al. Integrated omics: tools, advances and future approaches. Journal of Molecular Endocrinology, 2019. https://doi.org/10.1530/jme-18-0055 ↩

4. Bonnie Berger, Jian Peng, Mona Singh. Computational solutions for omics data. Nature Reviews Genetics, 2013. https://doi.org/10.1038/nrg3433 ↩