How to Align Your Own RNA-Seq Reads to the Human Genome

By the end of this you will have, for each RNA-seq library, a coordinate-sorted BAM file with splice-aware alignments to GRCh38, a gene-level count table, a transcript-level abundance table, and a QC report that tells you whether the numbers mean anything. You need the raw paired-end FASTQ files from your sequencing provider and a Linux machine with at least 64 GB of RAM and 500 GB of free disk. STAR’s human index alone needs roughly 32 GB of RAM to load, and a 40-million-read library produces a BAM of a few gigabytes. You also need 16 or more cores if you want this to finish in an hour rather than a day, plus conda or a container runtime. Everything below is command-line work on your own data. Nothing here is a clinical interpretation, and expression changes you find should go to a physician before they influence any decision about your health.

1. Know which kind of alignment you are doing, and pick the reference

The word “alignment” covers two different problems, and search results mix them constantly. Multiple sequence alignment (MSA) takes a handful of homologous protein or gene sequences from different species and arranges them column by column to expose conserved residues and evolutionary relationships. That is the problem Clustal Omega, MAFFT, and MUSCLE solve, and it remains genuinely hard: benchmark studies show that even the best MSA programs degrade badly on large, divergent, or domain-shuffled sequence sets, with no single method dominating across test categories.¹ Extending those methods from small protein families to complete genomic sequences was recognized two decades ago as a distinct computational problem with its own requirements.²

RNA-seq alignment is the other problem. You have tens of millions of short reads sampled from a transcriptome, typically 100 or 150 bp and paired. You want to place each one on a 3.1-billion-base reference genome, allowing a read to be split across an intron because it crossed an exon-exon junction. This is short-read mapping with spliced alignment. The tools are STAR and HISAT2, plus Salmon and kallisto for quantification without a genomic BAM.

For the reference, we use the GENCODE primary assembly rather than the full analysis set with alternate haplotypes. Alt contigs create multimapping reads that most RNA-seq quantifiers handle poorly. Download the primary assembly FASTA and the matching comprehensive GTF from the GENCODE FTP site, and keep the release number in every directory name you create:

# from the GENCODE human FTP, release 44 (or newer)
gunzip -c GRCh38.primary_assembly.genome.fa.gz > GRCh38.primary.fa
gunzip -c gencode.v44.primary_assembly.annotation.gtf.gz > gencode.v44.gtf

Pin the release. If you are building a longitudinal profile and you re-align a later blood draw against GENCODE v47 while the first draw was aligned against v44, some of the differences you see between timepoints will be annotation changes rather than biology.

2. Inspect and trim the FASTQs

Run FastQC first, then trim only if the report justifies it. Modern library preps on patterned flow cells often need nothing beyond adapter removal, and aggressive quality trimming loses real signal.

fastqc -t 8 sample1_R1.fastq.gz sample1_R2.fastq.gz
fastp -i sample1_R1.fastq.gz -I sample1_R2.fastq.gz \
  -o sample1_R1.trim.fq.gz -O sample1_R2.trim.fq.gz \
  --detect_adapter_for_pe --trim_poly_g --length_required 36 \
  --thread 8 --json sample1.fastp.json --html sample1.fastp.html

--trim_poly_g matters on two-colour chemistry (NovaSeq, NextSeq), where a dark cycle at the end of a read is called as G and produces long poly-G tails that will not align. --length_required 36 drops reads too short to place uniquely. Note the read length reported by FastQC, because you need it for the index in step 4.

3. Determine strandedness before you count anything

Getting strandedness wrong is the most common way to produce a count matrix that looks plausible and is wrong, typically losing 50 to 90 percent of your reads or assigning them to the antisense gene. Do not trust the kit name alone. Align a small subset first, then ask RSeQC:

# convert the GTF to BED12 once
gtfToGenePred -genePredExt gencode.v44.gtf gencode.v44.genePred
genePredToBed gencode.v44.genePred gencode.v44.bed

infer_experiment.py -r gencode.v44.bed -i sample1.subset.bam -s 400000

If the fraction explained by 1+-,1-+,2++,2-- is above 0.85, the library is reverse-stranded (the dUTP protocol used by Illumina TruSeq Stranded and most modern kits). That corresponds to column 4 of STAR’s ReadsPerGene.out.tab, featureCounts -s 2, and Salmon library type ISR. If the other pattern dominates, the library is forward-stranded. That means column 3, -s 1, and ISF. If both are near 0.5, it is unstranded. That means column 2, -s 0, and IU. Write the answer down in your sample sheet.

4. Build the STAR index

Build once per combination of genome, annotation and read length. The --sjdbOverhang value should be read length minus one, so 100 for 101 bp reads and 149 for 150 bp reads. The default of 100 works acceptably for most read lengths, but if you have 150 bp reads it is worth building the matching index.

STAR --runMode genomeGenerate --runThreadN 16 \
  --genomeDir star_idx_gencode44_sjdb149 \
  --genomeFastaFiles GRCh38.primary.fa \
  --sjdbGTFfile gencode.v44.gtf \
  --sjdbOverhang 149

Expect roughly 45 minutes on 16 cores and a peak of about 32 GB of RAM. On a machine with less memory, add --genomeSAsparseD 2, which halves the suffix array density and cuts memory at some cost in speed. The alternative is HISAT2, whose index fits in roughly 8 GB. In the aligner comparisons we find credible, STAR and HISAT2 both map the large majority of reads correctly and downstream gene-level conclusions are more sensitive to the quantification step than to which of the two spliced aligners produced the BAM.³ We use STAR because its two-pass mode handles novel junctions well and because --quantMode TranscriptomeSAM gives you a transcriptome BAM for free.

5. Align with STAR in two-pass mode

Two-pass alignment runs the whole library once to discover splice junctions, inserts those junctions into the index, and re-aligns. It measurably improves the placement of reads spanning novel or lowly-covered junctions, which is exactly where single-sample analysis is weakest.

STAR --genomeDir star_idx_gencode44_sjdb149 \
  --readFilesIn sample1_R1.trim.fq.gz sample1_R2.trim.fq.gz \
  --readFilesCommand zcat --runThreadN 16 \
  --twopassMode Basic \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMattributes NH HI AS nM NM MD \
  --outSAMstrandField intronMotif \
  --quantMode GeneCounts TranscriptomeSAM \
  --outFilterType BySJout \
  --outFilterMultimapNmax 20 \
  --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 \
  --outFilterMismatchNmax 999 --outFilterMismatchNoverReadLmax 0.04 \
  --alignIntronMin 20 --alignIntronMax 1000000 \
  --alignMatesGapMax 1000000 \
  --limitBAMsortRAM 30000000000 \
  --outFileNamePrefix sample1.

Those filter settings are the ENCODE long-RNA parameters and they are a reasonable default. --outFilterType BySJout re-filters reads so that only junctions surviving the final filter are kept, which reduces spurious novel junctions. --outFilterMismatchNoverReadLmax 0.04 allows four percent mismatches, generous enough for personal variation and RNA editing. On 16 cores, a 40-million-read paired library takes roughly 25 to 40 minutes.

Then index the BAM:

samtools index -@ 8 sample1.Aligned.sortedByCoord.out.bam

sample1.ReadsPerGene.out.tab now holds gene-level counts in four columns. They are gene ID, unstranded, forward and reverse. Take the column your step 3 result dictates.

6. Quantify transcripts with Salmon as a second, independent estimate

Gene-level counts from STAR are fine for most purposes, but they discard isoform structure and they assign multimapping reads crudely. Salmon’s selective alignment against a decoy-aware index gives transcript-level abundances in a few minutes. It uses an explicit model for GC bias, positional bias, and fragment length. Choice of quantification tool changes which genes come out as differentially expressed more than choice of aligner does, so running two methods and checking that they agree is cheap insurance.³

grep "^>" GRCh38.primary.fa | cut -d " " -f1 | sed 's/>//' > decoys.txt
cat gencode.v44.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz
salmon index -t gentrome.fa.gz -d decoys.txt -p 16 -i salmon_idx_v44 --gencode

salmon quant -i salmon_idx_v44 -l ISR \
  -1 sample1_R1.trim.fq.gz -2 sample1_R2.trim.fq.gz \
  -p 16 --seqBias --gcBias --posBias \
  -o sample1_salmon

The decoy set (the whole genome) is what stops intronic and intergenic reads from being force-fitted onto transcripts, which otherwise inflates the apparent expression of genes with similar sequence. Index construction needs about 16 GB of RAM and half an hour. Collapse to gene level in R with tximport:

library(tximport); library(readr)
tx2gene <- read_tsv("tx2gene_gencode_v44.tsv")  # transcript_id, gene_id
txi <- tximport(c(s1="sample1_salmon/quant.sf"), type="salmon",
                tx2gene=tx2gene, countsFromAbundance="lengthScaledTPM")

Then correlate txi$counts against STAR’s column for genes with more than ten counts. Spearman correlation below about 0.95 on log counts means something is off, usually strandedness or a mismatch between the annotation used for the index and the one used for tx2gene.

7. Run post-alignment QC and read the numbers

An alignment that completes is not an alignment that is usable. Collect metrics with Picard and RSeQC, then aggregate with MultiQC.

picard CollectRnaSeqMetrics \
  I=sample1.Aligned.sortedByCoord.out.bam \
  O=sample1.rnaseq_metrics.txt \
  REF_FLAT=gencode.v44.refFlat.txt \
  STRAND=SECOND_READ_TRANSCRIPTION_STRAND \
  RIBOSOMAL_INTERVALS=rRNA.interval_list

read_distribution.py -i sample1.Aligned.sortedByCoord.out.bam -r gencode.v44.bed
multiqc .

For a polyA-selected library from whole blood or PBMCs, we look first at uniquely mapped reads above 80 percent of input. We also want PCT_CODING_BASES plus PCT_UTR_BASES above roughly 0.7, PCT_INTRONIC_BASES below about 0.2, and PCT_RIBOSOMAL_BASES below 0.05. Finally, MEDIAN_5PRIME_TO_3PRIME_BIAS should be near 1.0. A high intronic fraction with an otherwise clean run usually means nuclear or pre-mRNA contamination or a total-RNA rather than polyA prep. A 3’ bias well below 1.0 indicates RNA degradation, and degraded samples systematically under-report long transcripts, which you will misread as biology if you compare them against an intact sample from another timepoint.

There are two blood-specific checks. Globin transcripts such as HBB, HBA1, and HBA2 can consume 50 to 70 percent of reads in whole blood without globin depletion. Mitochondrial transcripts above roughly 20 percent suggest cell stress during collection or processing. Both leave you with far less usable depth for everything else, and both are properties of the sample. Neither is a property of the alignment.

8. Wire it together so you can reproduce it in a year

Every step above is a tool with its own version, index, and set of flags. The assembled workflow is where irreproducibility enters. This problem is well documented across omics, and the practical answer is a workflow manager that pins versions and records provenance rather than a folder of shell scripts.⁴ We run nf-core/rnaseq under Nextflow with a container profile:

nextflow run nf-core/rnaseq -r 3.14.0 \
  --input samplesheet.csv --outdir results \
  --fasta GRCh38.primary.fa --gtf gencode.v44.gtf \
  --aligner star_salmon --profile docker \
  -resume

That pipeline does everything in steps 2 through 7 with the same tools, plus MarkDuplicates, RSeQC, Qualimap, and a MultiQC report. We recommend running the manual commands once so you understand what each stage does, then switching to the pipeline for every subsequent sample so that the twelfth blood draw is processed identically to the first.

Common problems

Very low unique mapping rate (under 60 percent) with a high “too short” percentage in Log.final.out almost always means adapter contamination or rRNA. The genome is usually not at fault. Check the fastp adapter report, and blast a few hundred unmapped reads against nt to identify what they are before you change alignment parameters.

STAR exiting with a segmentation fault or “EXITING because of fatal error: not enough memory” during sorting is a --limitBAMsortRAM problem. Set it explicitly to a value below your available RAM, and remember that the genome index itself occupies about 32 GB before sorting starts.

Counts that look reasonable but sit an order of magnitude below expectation for well-known genes: strandedness. Re-read step 3 and take the correct column.

A large disagreement between STAR gene counts and Salmon gene counts concentrated in gene families (HLA, immunoglobulins, ribosomal proteins) is expected rather than an error. These are multimapping-heavy regions where the two tools resolve ambiguity differently, and neither is exactly right. For HLA specifically, use a dedicated genotype-aware quantifier rather than either general tool.

Deduplication: for standard bulk RNA-seq without unique molecular identifiers, do not deduplicate. High-expression transcripts legitimately produce identical fragments, and marking them as duplicates biases against exactly the genes you can measure best. Mark duplicates for the QC metric, then quantify from the un-deduplicated BAM.

Finally, comparing timepoints. Within one person, the biological signal you are after is often a 1.3- to 2-fold change. That is the same order as batch effects from library prep, sequencing run, and time of day of the blood draw. Keep the protocol fixed, record collection time, and consider spiking in ERCC controls so you have an external scale. Where a change looks meaningful, a physician needs to be the one who decides whether it means anything about your health.

Footnotes

1. Julie D. Thompson, Benjamin Linard, Odile Lecompte, et al. A Comprehensive Benchmark Study of Multiple Sequence Alignment Methods: Current Challenges and Future Perspectives. PLoS ONE, 2011. https://doi.org/10.1371/journal.pone.0018093 ↩

2. Odile Lecompte, Julie D Thompson, Frédéric Plewniak, et al. Multiple alignment of complete sequences (MACS) in the post-genomic era. Gene, 2001. https://doi.org/10.1016/s0378-1119(01)00461-9 ↩

3. Isaac D. Raplee, Alexei V. Evsikov, Caralina Marín de Evsikova. Aligning the Aligners: Comparison of RNA Sequencing Data Alignment and Gene Expression Quantification Tools for Clinical Breast Cancer Research. Journal of Personalized Medicine, 2019. https://doi.org/10.3390/jpm9020018 ↩ ↩²

4. Tianle Ma, Aidong Zhang. Omics Informatics: From Scattered Individual Software Tools to Integrated Workflow Management Systems. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 2017. https://doi.org/10.1109/tcbb.2016.2535251 ↩