How to Choose and Build a Reference Genome for RNA-seq

By the end of this guide you will have a directory containing a human genome FASTA, a matching annotation GTF, a STAR index for spliced alignment, and a decoy-aware salmon index for transcript quantification, plus a short set of checks that confirm the three files agree with each other. You need roughly 100 GB of free disk, 32 GB of RAM (64 GB is more comfortable for STAR), a Linux machine or container, and the following tools: samtools, STAR (2.7.x), salmon (1.10 or later), gffread, and wget. If you have your own RNA-seq FASTQs from a sequencing provider, keep them nearby, since the read length determines one index parameter. Everything here is about measurement: getting numbers out of your transcriptome that mean what you think they mean. Nothing here is clinical interpretation, and if a result matters to your health you need a physician and a CLIA-certified laboratory, not a laptop pipeline.

1. Understand what “reference genome” means for RNA-seq

A reference genome is a single assembled consensus sequence for a species, built from a small number of donors and used as a coordinate system that everyone shares. For humans that coordinate system began with the Human Genome Project ¹ and is now maintained by the Genome Reference Consortium as the GRCh series. When you say “chr7:140,753,336,” you are naming a position in that shared system, and the only reason the number is meaningful to anyone else is that they are using the same build.

RNA-seq complicates this because messenger RNA is spliced. Your reads come from mature transcripts in which introns have been removed, so a read spanning an exon-exon junction will not align contiguously to genomic DNA. Two strategies follow from that. You can align to the genome with a splice-aware aligner that allows long gaps and is guided by a junction database, or you can align to a transcriptome FASTA in which splicing is already resolved and treat quantification as an assignment problem. Both strategies need an annotation, a GTF or GFF3 file listing where genes, transcripts, and exons live in genome coordinates, and the quality of that annotation drives the quality of your quantification at least as much as the aligner does ².

So “the reference” is really three coupled artifacts: a genome FASTA, an annotation, and an index built from both. The most common source of silently wrong results is a mismatch among the three.

2. Pick the genome build

Use the GRCh38 no-alt analysis set from NCBI. Concretely:

mkdir -p ref && cd ref
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
gunzip GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
mv GCA_000001405.15_GRCh38_no_alt_analysis_set.fna GRCh38.primary.fa
samtools faidx GRCh38.primary.fa

Three properties make this the right file. It excludes ALT contigs, the alternate haplotype scaffolds that GRCh38 adds for highly polymorphic regions such as the MHC. Aligners that are not ALT-aware will distribute reads across a primary contig and its alternates, turning unique alignments into multi-mappers and depressing expression estimates for exactly the genes you are most likely to care about. It hard-masks the pseudoautosomal regions on chrY, so X and Y homologous sequence is not duplicated. And it uses UCSC-style names (chr1, chrM), which match GENCODE’s annotation files.

The alternative builds are worth naming. Ensembl distributes a primary assembly with bare names (1, MT) that is equivalent in sequence but will not interoperate with a GENCODE GTF without renaming. The full analysis set with decoys and HLA contigs (hs38DH) is the right choice for whole-genome DNA alignment and is unnecessary overhead for RNA. T2T-CHM13v2.0 is a genuinely complete assembly and resolves regions GRCh38 cannot, but it is a different coordinate system with a younger annotation ecosystem, and most public expression references, GTEx quantifications, and clinical variant databases are still on GRCh38. Use CHM13 when you have a specific question about a region GRCh38 gets wrong, and keep GRCh38 as your primary. GRCh37/hg19 is only appropriate when you must match an older dataset exactly.

3. Pick the annotation

Use GENCODE, matching the release to the build:

wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/gencode.v47.primary_assembly.annotation.gtf.gz
gunzip gencode.v47.primary_assembly.annotation.gtf.gz

GENCODE is the annotation used by ENCODE and GTEx, it is comprehensive (on the order of 20,000 protein-coding genes and well over 200,000 transcript isoforms in recent releases), and its chromosome naming matches the NCBI analysis set. RefSeq is more conservative and better curated per transcript, which some clinical pipelines prefer, but it annotates fewer isoforms and will change your gene-level counts. The two are not interchangeable, and comparing counts generated against different annotations is a common and avoidable error.

Choose a release and record it. Annotation grows between releases, and a gene that gains an isoform will get slightly different transcript-level estimates even from identical reads. Pin the version in your pipeline config the way you would pin a dependency.

The comprehensive GTF includes pseudogenes, lncRNAs, and retained-intron transcripts. Keep them. Filtering them out does not make the reads from those loci disappear, it just forces them to be assigned somewhere else, and the presence of processed pseudogenes in the reference is what lets an aligner correctly recognize that a read is ambiguous rather than confidently wrong ³.

4. Verify that the FASTA and GTF agree

Do this before you spend an hour building an index. The check is cheap:

# sequence names in the FASTA
cut -f1 GRCh38.primary.fa.fai | sort > /tmp/fa.names

# sequence names referenced by the GTF
grep -v '^#' gencode.v47.primary_assembly.annotation.gtf | cut -f1 | sort -u > /tmp/gtf.names

# names the GTF uses that the FASTA does not have
comm -13 /tmp/fa.names /tmp/gtf.names

The third command should print nothing. If it prints 1, 2, MT, you have mixed Ensembl and UCSC naming conventions and need to rename one side. If it prints a handful of scaffolds, you downloaded a primary-assembly FASTA against a full annotation, or vice versa.

Also check that the GTF’s coordinates are sane for a few known genes, and confirm the mitochondrial contig is present. Mitochondrial transcripts are often 15 to 30 percent of reads in a polyA library from metabolically active tissue, and a reference missing chrM will quietly push those reads into nuclear mitochondrial pseudogenes.

5. Build the STAR index

STAR is our default spliced aligner. The systematic RGASP evaluation of spliced alignment programs found substantial differences between tools in junction recall and in how they handle indels and short exons, with no single program winning every category, and STAR sits in the group that combines high accuracy with tractable runtime ⁴. Earlier comparative work made the same general point: alignment algorithm choice measurably changes downstream expression and junction calls, so the choice should be deliberate and recorded ⁵.

STAR --runMode genomeGenerate \
     --runThreadN 16 \
     --genomeDir star_idx_GRCh38_v47 \
     --genomeFastaFiles GRCh38.primary.fa \
     --sjdbGTFfile gencode.v47.primary_assembly.annotation.gtf \
     --sjdbOverhang 100 \
     --limitGenomeGenerateRAM 48000000000

Set --sjdbOverhang to read length minus one. For 2x101 bp reads, 100 is exactly right, and 100 is also a reasonable default for mixed read lengths because the penalty for a slightly wrong value is small. Building takes roughly 40 to 90 minutes on 16 cores and produces about 30 GB on disk. If the job is killed, you ran out of RAM: STAR needs about 32 GB for human with annotation, and --limitGenomeGenerateRAM only caps sorting behavior, it does not reduce the requirement.

Passing the GTF at index time is what makes the index junction-aware. The aligner then has a database of known splice sites and can place a read that crosses a junction with only a few bases on one side, which is the case pure genomic alignment handles worst. Approaches that build an explicit junction-augmented reference make the same bet from the other direction, extending the genome with sequences that represent known exon-exon joins so that junction-spanning reads align contiguously ⁶.

For a typical alignment run:

STAR --genomeDir star_idx_GRCh38_v47 \
     --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
     --readFilesCommand zcat \
     --runThreadN 16 \
     --outSAMtype BAM SortedByCoordinate \
     --outSAMattributes NH HI AS nM MD \
     --quantMode GeneCounts TranscriptomeSAM \
     --outFilterMultimapNmax 20 \
     --outFileNamePrefix sample_

--quantMode GeneCounts gives you an HTSeq-equivalent count table for free. TranscriptomeSAM emits a transcriptome-coordinate BAM you can feed to RSEM if you want that route. Raising --outFilterMultimapNmax from the default of 10 to 20 costs little and keeps more information available for tools that model multi-mapping.

6. Build a decoy-aware salmon index

For transcript-level quantification we run salmon in parallel with STAR, because selective alignment against a transcriptome is much faster and handles isoform-level ambiguity with an explicit model rather than a heuristic. The key detail is decoys. A transcriptome-only index has nowhere to put reads originating from introns, unannotated regions, or genomic DNA contamination, so those reads get forced onto whichever transcript looks least bad. Adding the full genome as decoy sequence gives them somewhere honest to go.

gffread gencode.v47.primary_assembly.annotation.gtf \
        -g GRCh38.primary.fa \
        -w gencode.v47.transcripts.fa

grep '^>' GRCh38.primary.fa | cut -d ' ' -f1 | tr -d '>' > decoys.txt
cat gencode.v47.transcripts.fa GRCh38.primary.fa > gentrome.fa

salmon index -t gentrome.fa -d decoys.txt -p 16 -i salmon_idx_GRCh38_v47 -k 31

Keep -k 31 unless your reads are shorter than about 75 bp, in which case drop to 23 or 25. Budget 32 GB of RAM and an hour. If you build the transcript FASTA directly from GENCODE’s own transcripts.fa.gz rather than with gffread, add --gencode to the index command so salmon splits the pipe-delimited FASTA headers and keeps only the ENST identifier.

Then quantify:

salmon quant -i salmon_idx_GRCh38_v47 -l A \
  -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  -p 16 -o sample_salmon

Aggregate transcript estimates to genes with tximport in R, using a transcript-to-gene map derived from the same GTF. Deriving that map from a different GTF version than the index is a real and frequent source of silently dropped genes.

7. Check the mapping before you trust the counts

Read sample_Log.final.out from STAR first. For a good polyA-selected human library against GRCh38, uniquely mapped reads are typically above 85 percent, mismatch rate per base under 0.5 percent, and reads mapped to too many loci under 5 percent. A uniquely-mapped rate in the 60s usually means adapter contamination, rRNA carryover, or the wrong species. A high “too short” percentage usually means untrimmed adapters rather than a reference problem.

Then run Picard CollectRnaSeqMetrics with a refFlat file and rRNA interval list derived from your GTF. The numbers that matter are the fraction of bases in coding, UTR, intronic, and intergenic regions, and the median 5’ to 3’ coverage bias. High intronic fraction points to incomplete splicing or genomic DNA in the library. Strong 3’ bias points to RNA degradation. Both change what your expression estimates mean, and neither is fixed by changing reference.

Depth matters too. Deep sequencing of human B-cell transcriptomes showed that detection of low-abundance transcripts keeps increasing well past the depths typical of routine experiments, so a gene reading zero is often a statement about your library rather than about the cell ⁷.

8. Know when the standard reference is the wrong reference

A single consensus genome is a poor match for any individual at polymorphic sites. When a read carries a non-reference allele, it accumulates mismatches and is slightly less likely to align, which biases allele-specific expression toward the reference allele. Work on calling SNPs from RNA-seq data documents this and the filtering strategies needed to handle it ⁸. If you are measuring allelic imbalance rather than total expression, build a personalized or diploid reference from your own genotypes, or use a tool that models reference bias explicitly.

Clinical-adjacent applications raise the bar further. Reviews of RNA-seq for diagnosis emphasize that the analytical choices, including the annotation, determine whether events like aberrant splicing are detectable at all, and that interpretation requires a clinical framework around the assay ⁹. A splice event you find in your own data is a measurement, and turning it into a conclusion about your health is a job for a clinician working with validated testing.

Direct RNA sequencing on nanopore platforms sidesteps reverse transcription and reads native molecules end to end, which changes the reference question: you are matching full-length isoforms rather than reassembling them from short fragments. It remains higher in per-base error and lower in throughput, so short-read Illumina RNA-seq against GRCh38 stays the practical default for quantification ¹⁰.

Common problems

The index and the GTF disagree. If you rebuild an annotation but reuse an old index, STAR will happily align and produce counts keyed to stale gene identifiers. Name your index directory after both the build and the annotation release, as in star_idx_GRCh38_v47, and never overwrite in place.

Chromosome naming mismatch. Mixing Ensembl FASTA with GENCODE GTF gives zero overlap and an index that builds without error but assigns almost nothing. The comm check in step 4 catches it in seconds.

ALT contigs in the reference. If you used a full GRCh38 FASTA that includes ALT scaffolds with a non-ALT-aware aligner, expression for MHC genes and other polymorphic loci will be depressed and multi-mapping rates inflated. Rebuild against the no-alt analysis set.

Gene-level counts differ between STAR and salmon. Some difference is expected because the two handle multi-mapping reads differently, which is a known and well-characterized source of divergence ³. Large systematic differences, particularly inflated salmon estimates for genes with retained introns, usually mean you built the salmon index without decoys.

Out of memory during index building. STAR needs about 32 GB for human. There is no flag that removes this requirement. Use a larger machine or use a lightweight quantifier.

sjdbOverhang set wrong. A value far from read length minus one degrades junction detection for reads with short overhangs. It is not fatal, but rebuild the index if you switch read lengths substantially, for example from 50 bp to 150 bp.

Downloading a soft-masked genome. Files with sm or rm in the name have repeats in lowercase or replaced by N. Lowercase is harmless for STAR and salmon. Hard-masked files are not, and will lose real reads. Check with grep -c 'N\{100\}' if unsure.

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. P.S Harper. Mapping and sequencing the human genome. Endeavour, 1989. https://doi.org/10.1016/0160-9327(89)90035-5 ↩

2. Manuel Garber, Manfred G Grabherr, Mitchell Guttman, et al. Computational methods for transcriptome annotation and quantification using RNA-seq. Nature Methods, 2011. https://doi.org/10.1038/nmeth.1613 ↩

3. Gabrielle Deschamps-Francoeur, Joël Simoneau, Michelle S. Scott. Handling multi-mapped reads in RNA-seq. Computational and Structural Biotechnology Journal, 2020. https://doi.org/10.1016/j.csbj.2020.06.014 ↩ ↩²

4. The RGASP Consortium, Pär G Engström, Tamara Steijger, et al. Systematic evaluation of spliced alignment programs for RNA-seq data. Nature Methods, 2013. https://doi.org/10.1038/nmeth.2722 ↩

5. Gregory R. Grant, Michael H. Farkas, Angel D. Pizarro, et al. Comparative analysis of RNA-Seq alignment algorithms and the RNA-Seq unified mapper (RUM). Bioinformatics, 2011. https://doi.org/10.1093/bioinformatics/btr427 ↩

6. Yaron S. Butterfield, Maayan Kreitzman, Nina Thiessen, et al. JAGuaR: Junction Alignments to Genome for RNA-Seq Reads. PLoS ONE, 2014. https://doi.org/10.1371/journal.pone.0102398 ↩

7. Jonathan M. Toung, Michael Morley, Mingyao Li, et al. RNA-sequence analysis of human B-cells. Genome Research, 2011. https://doi.org/10.1101/gr.116335.110 ↩

8. Emma M. Quinn, Paul Cormican, Elaine M. Kenny, et al. Development of Strategies for SNP Detection in RNA-Seq Data: Application to Lymphoblastoid Cell Lines and Evaluation Using 1000 Genomes Data. PLoS ONE, 2013. https://doi.org/10.1371/journal.pone.0058815 ↩

9. Guillermo Marco-Puche, Sergio Lois, Javier Benítez, et al. RNA-Seq Perspectives to Improve Clinical Diagnosis. Frontiers in Genetics, 2019. https://doi.org/10.3389/fgene.2019.01152 ↩

10. Marcel C. Van Verk, Richard Hickman, Corné M.J. Pieterse, et al. RNA-Seq: revelation of the messengers. Trends in Plant Science, 2013. https://doi.org/10.1016/j.tplants.2013.02.001 ↩