How to Generate and Interpret STAR Gene Counts from Your RNA-seq Data

By the end of this guide you will have a gene-by-sample counts matrix derived from your own FASTQ files, a set of alignment quality metrics you can defend, and a TPM matrix suitable for exploratory work. You need the raw reads as paired-end gzipped FASTQ, a machine with at least 40 GB of RAM and a dozen cores, and STAR 2.7.10b or newer. You also need roughly 100 GB of free disk per sample-batch. You also need a matched genome FASTA and GTF annotation from the same source and release. Everything below assumes human data aligned to GRCh38 with GENCODE annotation, which is what we use, and the same logic transfers to any species with a decent GTF.

The reason to run STAR with --quantMode GeneCounts rather than piping a BAM into a separate counting step is that STAR can produce the counts during the alignment pass at essentially no extra cost, and the resulting numbers match what htseq-count would produce in its default union mode. You get a BAM for inspection and a count table from a single command. Uniformly processed public repositories built on this kind of pipeline exist precisely so that expression values from different studies can be compared without re-deriving them each time ¹.

1. Get matched reference files and build the index

The single most common source of silent failure in RNA-seq quantification is a genome and annotation that disagree about contig names or coordinates. Download both from the same GENCODE release and do not mix an Ensembl GTF with a UCSC FASTA, because one uses 1 and the other uses chr1, and STAR will build an index that quietly annotates nothing.

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

Use the primary assembly. Do not use the full analysis set with alternate haplotypes. Alt contigs duplicate sequence and inflate the multimapping rate for genes in the MHC and other polymorphic regions.

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

--sjdbOverhang should be your read length minus one. For 101 bp reads, 100 is exact. If you have mixed read lengths across a longitudinal set of samples, build one index at the maximum length minus one and reuse it, since the penalty for a slightly oversized overhang is small compared with the confusion of maintaining several indices. Index generation takes about an hour on 16 cores and produces roughly 30 GB of files. Human alignment needs about 32 GB of resident memory, so a 16 GB laptop will not work. If you cannot get that much RAM, use --genomeSAsparseD 2, which halves the suffix array memory footprint at the cost of some speed.

2. Align and count in one pass

The command below uses the ENCODE long-RNA parameter set, which is the default we recommend because it is widely used and its behavior on human data is well characterized.

STAR --runThreadN 16 \
     --genomeDir /ref/star_gencode47_oh100 \
     --readFilesIn sampleA_R1.fastq.gz sampleA_R2.fastq.gz \
     --readFilesCommand zcat \
     --outFileNamePrefix out/sampleA. \
     --outSAMtype BAM SortedByCoordinate \
     --quantMode GeneCounts \
     --outFilterType BySJout \
     --outFilterMultimapNmax 20 \
     --alignSJoverhangMin 8 \
     --alignSJDBoverhangMin 1 \
     --outFilterMismatchNmax 999 \
     --outFilterMismatchNoverReadLmax 0.04 \
     --alignIntronMin 20 \
     --alignIntronMax 1000000 \
     --alignMatesGapMax 1000000 \
     --outSAMattributes NH HI AS nM NM MD \
     --limitBAMsortRAM 20000000000

STAR’s speed is the reason this is practical on a workstation. The original algorithm paper reports alignment of roughly 550 million paired 76 nt reads per hour on a 12-core server, and a typical 40 million read pair library finishes in a few minutes once the index is loaded. If you are processing many samples, add --genomeLoad LoadAndKeep to hold the index in shared memory between runs. Drop it when using --twopassMode, which is incompatible.

Two-pass mode (--twopassMode Basic) runs the alignment twice, inserting novel splice junctions discovered in the first pass into the index for the second. It improves detection of reads spanning unannotated junctions and is worth the roughly doubled runtime if you care about splicing. For gene-level counts on well-annotated human samples the difference is modest, usually under one percent of assigned reads.

If you also want transcript-level estimates later, add --quantMode TranscriptomeSAM GeneCounts, which emits Aligned.toTranscriptome.out.bam for RSEM or similar. Transcript-level quantifiers and alignment-based gene counters do not always agree, and systematic comparisons of alignment and quantification tools show that the choice of tool measurably shifts downstream differential expression results, particularly for lowly expressed and multi-isoform genes ². Pick one pipeline and keep it constant across every timepoint you intend to compare.

3. Read Log.final.out before you read the counts

Open out/sampleA.Log.final.out first. It answers the question of how much misalignment is acceptable, which has no universal threshold but does have recognizable failure signatures.

The metric that matters most is “Uniquely mapped reads %”. For a polyA-selected human library on GRCh38, we expect 85 to 93 percent. Values between 70 and 85 percent are usable but deserve an explanation. Below 70 percent, stop and find out why before you interpret any expression value.

Then read the three rejection categories. ”% of reads unmapped: too short” is the catch-all bucket and is usually the informative one: it rises when adapters were not trimmed, when the library contains a foreign organism, or when read quality collapsed at the 3’ end. A few percent is normal, 20 percent is not. ”% of reads mapped to multiple loci” above 15 to 20 percent usually means residual ribosomal RNA, since rRNA repeats map everywhere. ”% of chimeric reads” should be a fraction of a percent unless you deliberately enabled chimeric detection for fusion discovery.

One more number worth recording per sample is the mean mapped length, reported in the same file. If it is far below your read length, soft-clipping is eating your reads, which again points at adapter contamination or a wrong reference.

4. Choose the correct column in ReadsPerGene.out.tab

sampleA.ReadsPerGene.out.tab has four columns and no header. Column 1 is the gene ID. Column 2 is counts for an unstranded library. Column 3 is counts assuming the read maps to the same strand as the gene (equivalent to htseq-count -s yes). Column 4 is counts assuming the reverse (equivalent to htseq-count -s reverse). The first four rows are not genes but summary lines: N_unmapped, N_multimapping, N_noFeature, and N_ambiguous. A read is counted for a gene only if it overlaps that gene and no other, which is why N_ambiguous exists.

You do not have to guess the strandedness of your library. Compare the total assigned counts across the three columns:

awk 'NR>4 {u+=$2; f+=$3; r+=$4} END {print "unstranded",u; print "forward",f; print "reverse",r}' \
    out/sampleA.ReadsPerGene.out.tab

For a dUTP-based stranded protocol, which covers Illumina TruSeq Stranded and most current kits, column 4 will hold nearly all the signal and column 3 will hold roughly 1 to 5 percent of it. For a genuinely unstranded library, columns 3 and 4 will be close to half of column 2 each. Anything in between, such as a 70:30 split, indicates a partially stranded or degraded library. Treat the sample as unstranded and note the anomaly.

Picking the wrong column typically drops 90 percent of your counts into N_noFeature and produces an expression profile dominated by antisense transcription. The error is not subtle. Check it once per library preparation protocol. The project is not the right unit for this check.

5. Build the counts matrix

With strandedness settled, assemble the per-sample files into one matrix. The following Python is what we use, and it keeps the summary rows aside rather than discarding them, because they are part of the QC record.

import pandas as pd, glob, os

STRAND_COL = 3  # 1-based column index: 2 unstranded, 3 forward, 4 reverse

frames, summaries = {}, {}
for path in sorted(glob.glob("out/*.ReadsPerGene.out.tab")):
    name = os.path.basename(path).split(".")[0]
    df = pd.read_csv(path, sep="\t", header=None, index_col=0,
                     usecols=[0, STRAND_COL - 1], names=["gene", name])
    summaries[name] = df.iloc[:4, 0]
    frames[name] = df.iloc[4:, 0]

counts = pd.DataFrame(frames)
qc = pd.DataFrame(summaries).T
counts.to_csv("counts_gencode47.tsv", sep="\t")
qc.to_csv("star_assignment_qc.tsv", sep="\t")

assigned = counts.sum()
print((assigned / (assigned + qc.sum(axis=1))).round(3))

The final line gives the fraction of all reads that ended up assigned to a unique gene. For a clean polyA human library, 60 to 80 percent is typical. Ribosomal-depleted libraries run lower because a larger share of reads falls in introns and intergenic space, and the two protocols give measurably different views of the same RNA population, with rRNA-depleted data capturing non-polyadenylated transcripts that polyA selection misses ³. Do not compare counts from the two protocols directly.

Keep the GENCODE version suffix on gene IDs (ENSG00000141510.17) in the stored matrix and strip it only at join time. Version numbers are the only record of which annotation produced the counts.

6. Convert counts to TPM

Raw counts are proportional to both expression and gene length, so they are comparable across samples for a given gene but not across genes within a sample. Transcripts per million corrects for length first and then for library size, which makes the values sum to one million in every sample.

The length to use with STAR gene counts is the union exon length: merge all exons of a gene across all its transcripts and sum the merged interval widths. Compute it once per annotation.

import pyranges as pr, numpy as np

gtf = pr.read_gtf("gencode.v47.primary_assembly.annotation.gtf")
exons = gtf[gtf.Feature == "exon"]
lengths = (exons.merge(by="gene_id")
                .lengths_by("gene_id"))          # bases per gene, exons merged
L = lengths.reindex(counts.index) / 1000.0        # kilobases

rpk = counts.div(L, axis=0)
tpm = rpk.div(rpk.sum(axis=0), axis=1) * 1e6
tpm.to_csv("tpm_gencode47.tsv", sep="\t")

Two caveats matter. First, this is a fixed length per gene, whereas transcript-level quantifiers compute an effective length that depends on the fragment length distribution and on which isoforms are expressed. When isoform usage shifts between samples, the fixed-length TPM will be slightly wrong, which is one of the reasons alignment-based and pseudoalignment-based quantification disagree ⁴. Second, never use TPM as input to a differential expression test. DESeq2 and edgeR expect raw integer counts and model the count distribution directly. Use TPM for visualization, clustering, and cross-gene comparison. It also suits feeding expression features into downstream integrative models, which is how expression is typically combined with genomic and proteomic layers in multi-omic analyses ⁵.

7. Sanity-check the matrix

Before trusting anything, run three checks. Confirm that mitochondrial genes account for a plausible share of the library, typically 5 to 30 percent in whole blood and tissue and much higher in stressed or dying cells. Confirm that the top ten genes by TPM look like the tissue you sampled, which for whole blood means hemoglobin genes dominating at tens of thousands of TPM. Finally, cluster the samples on the Spearman correlation of log-transformed TPM and confirm that replicates from the same timepoint are each other’s nearest neighbors. Correlations below about 0.9 between technical replicates of the same tissue indicate a processing problem. Biology does not explain them.

Record the STAR version and the index parameters in a manifest alongside the matrix. Note the GENCODE release and the strand column there as well. Reprocessing the same FASTQ files a year later with a different annotation will change gene-level counts in ways that look like biological drift, and the manifest is what lets you distinguish the two. Packaging the whole chain as a single reproducible unit is increasingly how these pipelines are distributed ⁶.

Expression values are measurements. They are not findings. A gene sitting at an unusual TPM relative to a reference cohort is a signal to investigate, and any interpretation with medical consequence belongs with a clinician who can see the full clinical picture.

Common problems

STAR is killed during alignment with no error message. This is almost always the out-of-memory killer. Check dmesg | tail. Human alignment needs roughly 32 GB plus the BAM sort buffer set by --limitBAMsortRAM, so on a 32 GB machine, set the sort buffer to 8 GB or write unsorted BAM and sort separately with samtools sort.

Nearly all reads land in N_noFeature. Either you read the wrong strand column, or the genome and GTF disagree about chromosome naming. Compare the first field of the GTF with samtools idxstats on the BAM.

The uniquely mapped rate is high but N_multimapping is enormous. Look for rRNA. Count the fraction of reads assigned to RNA45SN1 through RNA45SN5 and the 5S repeats, and if depletion failed you will see it immediately. Reads that multimap are never assigned to a gene by --quantMode GeneCounts regardless of --outFilterMultimapNmax, so a high multimapping rate directly reduces your usable depth.

Counts change after you rebuild the index. If --sjdbOverhang changed, junction-spanning read assignment changes slightly. If the GENCODE release changed, gene models changed and some IDs were retired. Both are expected and both are reasons to freeze the reference for the life of a longitudinal series.

A gene you care about shows zero counts in every sample. Confirm it exists in your GTF under the ID you are querying, and check whether it overlaps another gene on the opposite strand, in which case reads may be landing in N_ambiguous. Genes nested inside introns of other genes are the usual offenders.

Batch effects dominate your clustering. Library preparation date and sequencing run explain more variance than biology in most real datasets. Include them as covariates, and where possible, prepare all samples in a longitudinal series together. Multi-omic studies that combine transcriptomic data across collection events consistently have to model this explicitly before any biological signal is visible ⁷.

Footnotes

1. Mark Ziemann, Antony Kaspi, Assam El-Osta. Digital expression explorer 2: a repository of uniformly processed RNA sequencing data. GigaScience, 2019. https://doi.org/10.1093/gigascience/giz022 ↩

2. 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 ↩

3. Li Chen, Ruirui Yang, Tony Kwan, et al. Paired rRNA-depleted and polyA-selected RNA sequencing data and supporting multi-omics data from human T cells. Scientific Data, 2020. https://doi.org/10.1038/s41597-020-00719-4 ↩

4. Jorge A. Tzec‐Interián, Daianna González‐Padilla, Elsa B. Góngora‐Castillo. Bioinformatics perspectives on transcriptomics: A comprehensive review of bulk and single‐cell RNA sequencing analyses. Quantitative Biology, 2025. https://doi.org/10.1002/qub2.78 ↩

5. Chi-Sheng Wu, Hsin-Pai Li, Chia-Hsun Hsieh, et al. Integrated multi-omics analyses of oral squamous cell carcinoma reveal precision patient stratification and personalized treatment strategies. Cancer Letters, 2025. https://doi.org/10.1016/j.canlet.2025.217482 ↩

6. Ling-Hong Hung, Dylan Baker, Bill Flynn, et al. STAR Suite: Transcriptomics processing in a single binary through AI-assisted development. 2026. https://doi.org/10.64898/2026.03.09.710580 ↩

7. Afshin Beheshti, Kaushik Chakravarty, Homer Fogle, et al. Multi-omics analysis of multiple missions to space reveal a theme of lipid dysregulation in mouse liver. Scientific Reports, 2019. https://doi.org/10.1038/s41598-019-55869-2 ↩