How to Align RNA-Seq Reads with STAR

By the end of this you will have a STAR genome index for GRCh38 and a coordinate-sorted BAM per sample. You will also have a per-gene count table, a splice junction table, and a small set of QC numbers that tell you whether the library is worth analyzing. You need a Linux machine with at least 64 GB of RAM and 8 or more cores. The practical floor for a human index is 32 GB, and the sort step wants headroom. You also need roughly 200 GB of free disk for reference files plus intermediate BAMs. Software: STAR 2.7.11b, samtools 1.19 or newer, and either fastp or fastqc for read-level QC. Data: paired-end FASTQ files, the GENCODE primary assembly FASTA, and the matching GENCODE GTF annotation.

1. Check that you mean this kind of star alignment

Two unrelated methods share the name. In multiple sequence alignment, the “star” or “center-star” method picks one sequence as a center and aligns every other sequence to it pairwise, then merges those pairwise alignments into one multiple alignment. It is a classic approximation with a provable bound on sum-of-pairs cost, and variants of the idea show up in genomic alignment as a spanning star forest problem ¹. If that is what you came for, use MAFFT rather than a hand-rolled center-star: it handles thousands of sequences with progressive and iterative refinement strategies and is the default recommendation for DNA and protein multiple alignment today ^{2 3}. Benchmarks across accuracy and runtime consistently place MAFFT among the best general-purpose choices ⁴. For aligning several whole genomes to each other, the tooling is different again. It is built on suffix-tree anchoring and chaining ⁵. Older pipelines like MultiPipMaker are still useful for visualizing multi-species conservation ⁶.

The rest of this guide is about the other STAR: Spliced Transcripts Alignment to a Reference, the RNA-seq read aligner from Dobin and colleagues ⁷. It maps short reads to a reference genome while allowing a single read to span an intron, which is the defining problem of RNA-seq mapping. STAR does this by finding Maximal Mappable Prefixes with an uncompressed suffix array. It extends a read as far as it matches the genome, records that anchor, then restarts the search from the first mismatched base. An exon-exon junction therefore falls out of the seed search itself rather than being imposed afterward ⁷.

2. Decide whether STAR is the right tool

If you want gene- or transcript-level abundance and nothing else, use a lightweight quantifier like salmon or kallisto. They map reads against a transcriptome index rather than the genome, take about a tenth of the CPU time and a few gigabytes of RAM, and give comparable gene counts for well-annotated transcriptomes. They do not produce genomic coordinates. You cannot look at coverage over an intron, call variants, or see anything outside the annotation you fed them.

Use STAR when you need the alignment itself, for any of the following.

• novel splice junctions
• intron retention
• allele-specific expression
• RNA variant calling
• fusion detection
• a coverage track to look at in IGV

Compared with HISAT2, STAR is substantially faster at the cost of memory (about 30 GB resident for human versus roughly 5 GB), and in our experience it reports more novel junctions with a well-tuned two-pass run. HISAT2 is the better choice on a laptop or a small cloud instance. We run STAR because the machines are cheap and the extra output is worth having, and because the two-pass mode described below meaningfully improves junction sensitivity on samples with unusual splicing ⁸.

3. Get the reference files

Use GENCODE’s primary assembly. Do not use the “toplevel” Ensembl FASTA. The toplevel file includes haplotype scaffolds and patch sequences that create spurious multi-mapping and inflate your index. The primary assembly is what you want. It includes the chromosomes, the unlocalized and unplaced scaffolds, and the mitochondrial genome.

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.annotation.gtf.gz
gunzip GRCh38.primary_assembly.genome.fa.gz gencode.v47.annotation.gtf.gz

The FASTA and GTF must come from the same source and release. Mixing a UCSC FASTA (chr1) with an Ensembl GTF (1) is the single most common way to get a run that completes cleanly and produces almost no gene counts, because every contig name lookup fails silently at the quantification step.

4. Build the index

Index generation is a one-time cost of 40 to 90 minutes on 16 threads for human. The important parameter is --sjdbOverhang, which sets the length of genomic sequence flanking each annotated junction that STAR splices together into the index. Set it to your read length minus one.

STAR --runMode genomeGenerate \
     --runThreadN 16 \
     --genomeDir ref/star_gencode47_oh100 \
     --genomeFastaFiles ref/GRCh38.primary_assembly.genome.fa \
     --sjdbGTFfile ref/gencode.v47.annotation.gtf \
     --sjdbOverhang 100 \
     --limitGenomeGenerateRAM 64000000000

For 101 bp reads, 100 is exact. For mixed read lengths across a project, build once at 100 and reuse it. Dobin’s own guidance is that the ideal value is read length minus one but that 100 works well for most human data, and the loss from being slightly off is small ⁹. The finished index directory is roughly 30 GB.

If you have less than 32 GB of RAM, add --genomeSAsparseD 2, which halves the suffix array density and cuts memory to about 16 GB. Mapping gets roughly 30 to 40 percent slower and sensitivity is marginally lower. It is a reasonable trade on a memory-constrained box and a bad trade on a rented 64 GB instance.

5. Look at the reads before you map them

Run fastp on each pair. It trims adapters by paired-end overlap detection without your having to name the adapter, drops low-quality tails, and writes an HTML report in one pass.

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
      -o trim/sample_R1.fq.gz -O trim/sample_R2.fq.gz \
      --detect_adapter_for_pe --cut_tail --cut_tail_mean_quality 20 \
      --length_required 36 --thread 8 \
      --html qc/sample.fastp.html --json qc/sample.fastp.json

Two things to note in the report. Duplication rate above roughly 50 percent on a standard polyA library usually means low input RNA and over-amplification. A GC distribution with a sharp secondary peak often means rRNA carryover or adapter dimer. STAR soft-clips adapter remnants reasonably well, so trimming is not strictly required. Untrimmed adapter still shows up later as a high “too short” rate and is easier to remove now than to explain away later.

On the Sanger versus Illumina question that comes up around FASTQ handling: modern Illumina output is Phred+33, the same offset Sanger FASTQ uses, so no conversion is needed. The old Phred+64 encoding disappeared with Illumina 1.8 in 2011. If you inherit archived data and STAR reports implausible quality values, check with head -4 file.fastq and convert with seqtk seq -Q64 -V.

6. Run the alignment

This is the command we use for a standard human polyA or total-RNA library. It runs two-pass mapping, writes a sorted BAM, and produces a gene count table in the same run.

STAR --runThreadN 16 \
     --genomeDir ref/star_gencode47_oh100 \
     --readFilesIn trim/sample_R1.fq.gz trim/sample_R2.fq.gz \
     --readFilesCommand zcat \
     --outFileNamePrefix out/sample_ \
     --outSAMtype BAM SortedByCoordinate \
     --outBAMsortingThreadN 8 \
     --limitBAMsortRAM 30000000000 \
     --outSAMattributes NH HI AS nM NM MD \
     --outSAMunmapped Within \
     --quantMode GeneCounts \
     --twopassMode Basic \
     --outFilterMultimapNmax 20 \
     --alignSJoverhangMin 8 \
     --alignSJDBoverhangMin 1 \
     --outFilterMismatchNmax 999 \
     --outFilterMismatchNoverReadLmax 0.04 \
     --alignIntronMin 20 \
     --alignIntronMax 1000000 \
     --alignMatesGapMax 1000000

The filter block from --outFilterMultimapNmax down is the ENCODE long-RNA parameter set. It replaces the default mismatch cap with a length-proportional one (4 percent of the combined mate length) and lowers the required overhang for junctions already present in the annotation to 1 base, which recovers reads that barely cross a known junction.

Two-pass mode matters more than any single filter. In the first pass STAR maps everything and collects the junctions it finds, including ones absent from the GTF. It then rebuilds the splice junction database on the fly with those junctions inserted and maps again. Reads with only a few bases on the far side of a novel junction, which the first pass clips or discards, align correctly in the second ⁸. The cost is roughly double the mapping time. For multi-sample projects where you want consistent junction treatment, run pass one across all samples. Then concatenate the filtered SJ.out.tab files and feed them to every sample’s final run with --sjdbFileChrStartEnd. For a single sample, --twopassMode Basic is equivalent and simpler.

If you plan to quantify with RSEM or run salmon in alignment mode, add --quantMode TranscriptomeSAM GeneCounts and --quantTranscriptomeBan IndelSoftclipSingleend, which writes a second BAM in transcript coordinates.

Index the BAM when the run finishes:

samtools index -@ 8 out/sample_Aligned.sortedByCoord.out.bam

7. Read Log.final.out before anything else

STAR writes sample_Log.final.out, a short table that is the most information-dense QC artifact in the pipeline. The numbers we check, in order:

Uniquely mapped reads percentage. For human polyA RNA-seq against GRCh38 with GENCODE, expect 80 to 93 percent. Below 70 percent, stop and find out why before interpreting anything downstream.

Percentage of reads mapped to multiple loci. Typically 5 to 15 percent. A value above 30 percent usually means rRNA or other repeat content, which you can confirm by counting reads over RNA45SN1 and the 5.8S/28S loci.

Percentage of reads unmapped: too short. This is the catch-all bucket and the most informative failure signal. It counts reads where no alignment passed --outFilterScoreMinOverLread (default 0.66 of read length). High values mean adapter contamination, a wrong genome build, heavy degradation, or a contaminating organism. Take the unmapped reads (--outSAMunmapped Within, then samtools view -f 4) and BLAST a few hundred of them: the answer is usually obvious within a minute.

Mismatch rate per base. Around 0.2 to 0.6 percent for good Illumina data. Values above 1 percent point to quality problems or to a sample that is not the species you think it is.

Mapping speed, reported in millions of reads per hour, is useful for sizing future runs. On 16 modern cores expect roughly 200 to 400 million reads per hour in single-pass mode.

On alignment scores generally: STAR’s per-read AS tag is a local alignment score, not a probability, and its absolute value depends on read length. It is not comparable across datasets. Use NH (number of reported alignments) to filter to unique reads and MAPQ 255, which is STAR’s convention for uniquely mapped, rather than trying to set a threshold on AS.

8. Turn ReadsPerGene.out.tab into a count matrix

--quantMode GeneCounts produces sample_ReadsPerGene.out.tab with four columns. The columns are gene ID and counts for an unstranded library, followed by counts for reads on the same strand as the gene and counts for reads on the opposite strand. The first four rows are the special categories N_unmapped, N_multimapping, N_noFeature, and N_ambiguous.

Pick the column by determining strandedness empirically rather than trusting the kit documentation:

awk 'NR>4 {u+=$2; f+=$3; r+=$4} END {print "unstr",u," fwd",f," rev",r}' \
    out/sample_ReadsPerGene.out.tab

The common dUTP-based stranded protocols include Illumina TruSeq Stranded and NEBNext Ultra II Directional. For those, the reverse column will hold roughly 10 to 50 times the forward column, and you use column 4. If the two are within about 20 percent of each other, the library is unstranded and you use column 2. A ratio near 2:1 in either direction is a warning that something went wrong in library prep.

Then assemble the matrix in R, keeping only the counted rows:

files <- list.files("out", pattern="ReadsPerGene.out.tab$", full.names=TRUE)
mat <- sapply(files, function(f) {
  x <- read.delim(f, header=FALSE, skip=4)
  setNames(x[[4]], x[[1]])            # column 4 = reverse-stranded
})
colnames(mat) <- sub("_ReadsPerGene.out.tab$", "", basename(files))

N_noFeature above about 40 percent of mapped reads in a stranded library usually means you chose the wrong strand column or the annotation does not match the genome build.

9. Optional outputs worth knowing about

SJ.out.tab lists every detected junction with the count of uniquely mapping reads crossing it, the maximum spliced overhang, and a flag for whether the junction is annotated. Filtering to junctions with at least 5 unique reads and an overhang of at least 10 bases gives a usable novel-junction set for a single sample.

For fusion detection, STAR emits chimeric alignments with --chimSegmentMin 12 --chimJunctionOverhangMin 8 --chimOutJunctionFormat 1 --chimOutType Junctions, which is the input STAR-Fusion consumes. STAR-Fusion filters those candidates against annotation and paralog databases and is fast enough to run on every sample in a project ¹⁰. Treat any fusion call as a hypothesis requiring orthogonal confirmation, and any clinically relevant result as something to take to a physician or clinical laboratory rather than act on.

For variant calling from RNA, run two-pass STAR and then gatk MarkDuplicates. Follow that with gatk SplitNCigarSplitNCigarReads to break reads at the N operators before HaplotypeCaller. RNA variant calls carry allele-specific expression bias, RNA editing sites, and coverage that tracks expression rather than genotype. They are therefore a poor substitute for DNA genotyping and should be read as measurements rather than as a diagnosis.

Common problems

EXITING because of FATAL ERROR: could not create output file: STAR does not create directories. Make the output prefix directory first, and note that the prefix is a literal string prefix, so --outFileNamePrefix out/sample writes out/sampleLog.final.out unless you end it with an underscore or slash.

EXITING because of fatal error: not enough memory for BAM sorting: raise --limitBAMsortRAM to a value below your available RAM, or drop --outBAMsortingThreadN. Each sorting thread holds its own buffer, so more threads means more memory. Adding threads does not reduce it.

Nearly all counts landing in N_noFeature: almost always a chromosome naming mismatch between FASTA and GTF, or the wrong strand column. Compare samtools idxstats output against cut -f1 annotation.gtf | sort -u.

Very high % of reads unmapped: too short with normal FastQC output: check the genome build. Mapping human reads to GRCh37 when you built the index from GRCh38 still yields high unique mapping, but mapping to the wrong species or to a transcriptome FASTA by mistake collapses this metric immediately.

Two-pass mode producing different junction sets across samples in a project: this is expected, since each sample inserts its own novel junctions. If you are comparing junction usage across samples, use the shared --sjdbFileChrStartEnd approach so every sample maps against the same augmented index.

An index built with the wrong --sjdbOverhang: STAR will run, but reads crossing annotated junctions lose sensitivity. Rebuilding is cheaper than explaining an anomalous junction analysis six weeks later.

Footnotes

1. C. Thach Nguyen, Jian Shen, Minmei Hou, et al. Approximating the Spanning Star Forest Problem and Its Application to Genomic Sequence Alignment. SIAM Journal on Computing, 2008. https://doi.org/10.1137/070682150 ↩

2. Kazutaka Katoh, George Asimenos, Hiroyuki Toh. Multiple Alignment of DNA Sequences with MAFFT. Methods in Molecular Biology, 2009. https://doi.org/10.1007/978-1-59745-251-9_3 ↩

3. K. Katoh, H. Toh. Recent developments in the MAFFT multiple sequence alignment program. Briefings in Bioinformatics, 2008. https://doi.org/10.1093/bib/bbn013 ↩

4. Fabiano Sviatopolk-Mirsky Pais, Patrícia de Cássia Ruy, Guilherme Oliveira, et al. Assessing the efficiency of multiple sequence alignment programs. Algorithms for Molecular Biology, 2014. https://doi.org/10.1186/1748-7188-9-4 ↩

5. Michael Höhl, Stefan Kurtz, Enno Ohlebusch. Efficient multiple genome alignment. Bioinformatics, 2002. https://doi.org/10.1093/bioinformatics/18.suppl_1.s312 ↩

6. S. Schwartz. MultiPipMaker and supporting tools: alignments and analysis of multiple genomic DNA sequences. Nucleic Acids Research, 2003. https://doi.org/10.1093/nar/gkg579 ↩

7. Alexander Dobin, Carrie A. Davis, Felix Schlesinger, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 2012. https://doi.org/10.1093/bioinformatics/bts635 ↩ ↩²

8. Alexander Dobin, Thomas R. Gingeras. Mapping RNA‐seq Reads with STAR. Current Protocols in Bioinformatics, 2015. https://doi.org/10.1002/0471250953.bi1114s51 ↩ ↩²

9. Alexander Dobin, Thomas R. Gingeras. Optimizing RNA-Seq Mapping with STAR. Methods in Molecular Biology, 2016. https://doi.org/10.1007/978-1-4939-3572-7_13 ↩

10. Brian J. Haas, Alex Dobin, Nicolas Stransky, et al. STAR-Fusion: Fast and Accurate Fusion Transcript Detection from RNA-Seq. 2017. https://doi.org/10.1101/120295 ↩