Processing Your Own RNA-Seq Data from FASTQ to Gene Counts

By the end of this guide you will have taken paired-end FASTQ files from a blood RNA sequencing run and produced three artifacts: a transcript-level quantification table with effective lengths and estimated counts, a gene-level count matrix and TPM matrix collapsed from those transcripts, and a quality-control report that tells you whether the numbers are trustworthy. You need the raw FASTQs (typically 8 to 20 GB gzipped per sample for a 40 to 60 million read library), a machine with at least 16 GB of RAM and 200 GB of free disk, a Linux or macOS shell, and either conda/mamba or Docker. If you plan to align to the genome rather than quantify against a transcriptome, budget 32 GB of RAM for a human STAR index and several hours of CPU time. Everything below is measurement and interpretation. Nothing here is a diagnosis, and any result you intend to act on belongs in a conversation with a physician.

A note on terminology before the commands, since search results for “RNA processing” often mean something else. In cell biology, RNA processing refers to what the cell does to a nascent transcript: capping the 5’ end, splicing out introns, and cleaving and polyadenylating the 3’ end. RNA-seq processing, the subject of this guide, refers to what you do to the sequencing reads that come from those mature transcripts. The two are connected in a practical way. Most human library preparations select for polyadenylated RNA, so the cell’s own 3’ end processing determines what you can see. Transcripts that are not polyadenylated are largely absent from a standard dataset. Examples are most long non-coding RNA, histone mRNAs, mature tRNA and circular RNA. Transfer RNA in particular requires dedicated chemistry because its heavy base modifications and tight structure block ordinary reverse transcription.¹

1. Know what library you have before you touch it

Every downstream choice follows from the library chemistry, so establish it first. Ask the sequencing provider for four facts: selection method (poly-A capture or ribosomal RNA depletion), strandedness, read layout and length, and whether unique molecular identifiers (UMIs) were used. For blood, also ask whether globin mRNA was depleted. Without depletion, hemoglobin transcripts can consume 50 to 70 percent of reads in whole blood, gutting your effective depth for everything else.

Poly-A capture gives you clean mRNA quantification at lower depth and is what we prefer for a longitudinal expression baseline. Ribosomal depletion (Ribo-Zero style) keeps non-polyadenylated and partially degraded RNA, which matters if you care about long non-coding RNA or intronic signal, at the cost of needing roughly twice the depth for the same mRNA precision. Either way, remember the ceiling on interpretation: transcript abundance is a proxy for protein abundance, and the correlation between mRNA and protein levels is moderate at best because translation rate and protein turnover vary by orders of magnitude across genes.² This is the main reason we read RNA alongside proteomics rather than in place of it.

If you do not know the strandedness, you can infer it in step 4 rather than guessing.

2. Set up the environment and reference files

Use a pinned conda environment. Version drift between a quantifier and its index is the single most common source of irreproducible counts.

mamba create -n rnaseq -c conda-forge -c bioconda \
  fastp=0.23.4 fastqc=0.12.1 multiqc=1.21 \
  salmon=1.10.3 star=2.7.11b samtools=1.19 \
  subread=2.0.6 r-base=4.3 bioconductor-tximport
conda activate rnaseq

For references, we use GENCODE for human and match the annotation release to the genome build. Do not mix an Ensembl FASTA with a RefSeq GTF. The chromosome naming alone will cost you an afternoon.

mkdir -p ref && cd ref
BASE=https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_45
wget $BASE/gencode.v45.transcripts.fa.gz
wget $BASE/gencode.v45.annotation.gtf.gz
wget $BASE/GRCh38.primary_assembly.genome.fa.gz

Use the primary assembly, not the full analysis set with alternate haplotypes, unless you have a specific reason. Alt contigs create multimapping that most quantifiers handle poorly by default.

3. Trim adapters and low-quality tails

Run fastp rather than a separate FastQC plus Trimmomatic chain. It does adapter detection, quality trimming, poly-G correction for NovaSeq two-color chemistry, and produces its own QC report in one pass at roughly 100 million reads per hour on 8 threads.

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o trim_R1.fastq.gz  -O trim_R2.fastq.gz \
  --detect_adapter_for_pe \
  --trim_poly_g \
  --cut_tail --cut_tail_mean_quality 20 \
  --length_required 36 \
  --thread 8 \
  --json sample.fastp.json --html sample.fastp.html

Two deliberate choices here. We do not use aggressive quality trimming (--cut_right with a high threshold), because modern quantifiers are tolerant of a few low-quality bases and over-trimming biases effective transcript lengths, which shifts TPM values. We keep --length_required 36 because reads shorter than about 30 nucleotides multimap across paralogs and add noise without adding information.

If your library used UMIs, extract them before trimming with umi_tools extract and a pattern matching your kit, because once you trim you lose the positional anchor. UMI deduplication matters most for low-input libraries where PCR duplicates dominate. For a standard 500 ng blood RNA library at 50 million reads, deduplication typically changes gene-level counts by a few percent.

4. Quantify against the transcriptome with Salmon

For gene expression in a healthy adult with a well-annotated genome, we use selective alignment against the transcriptome rather than a full spliced genome alignment. It runs in 10 to 20 minutes per sample on 8 threads instead of several hours, it handles multimapping between isoforms probabilistically rather than discarding it, and the gene-level counts agree closely with STAR plus featureCounts for the great majority of genes. Build the index with decoy sequences from the genome so that reads originating from unannotated or intronic regions are not force-fit onto real transcripts.

cd ref
grep "^>" <(gunzip -c GRCh38.primary_assembly.genome.fa.gz) \
  | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v45.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz

salmon index \
  -t gentrome.fa.gz -d decoys.txt \
  -p 8 -i salmon_idx_v45 --gencode

The --gencode flag strips the pipe-delimited GENCODE headers down to the transcript ID, which saves you a string-parsing step later. Indexing needs roughly 16 GB of RAM and half an hour.

Now quantify. Set the library type to A on your first sample and let Salmon infer strandedness.

salmon quant \
  -i ref/salmon_idx_v45 -l A \
  -1 trim_R1.fastq.gz -2 trim_R2.fastq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  --numGibbsSamples 20 \
  -p 8 -o quant/sample01

Then read quant/sample01/lib_format_counts.json and confirm the inferred type. For a typical stranded dUTP protocol you should see ISR with over 95 percent of fragments consistent. If you see IU (unstranded) on a library you were told was stranded, stop and resolve that with the lab, because it usually means a sample swap or a mislabeled run.

The bias-correction flags are not optional in our view. GC bias correction in particular changes abundance estimates for GC-extreme transcripts substantially. It is the difference between a longitudinal series that tracks biology and one that tracks which day the library was made. The --numGibbsSamples 20 setting gives you posterior uncertainty per transcript, which you will want when an isoform-level result looks interesting. Transcripts whose Gibbs samples span an order of magnitude are ambiguity and should not be read as measurements.

5. Optionally align to the genome as well

Transcriptome quantification cannot see what is not in the annotation. If you care about novel splice junctions, fusion transcripts, allele-specific expression, or RNA editing, you need a coordinate-sorted BAM. We run STAR in two-pass mode for junction discovery.

STAR --runMode genomeGenerate \
  --genomeDir ref/star_idx \
  --genomeFastaFiles ref/GRCh38.primary_assembly.genome.fa \
  --sjdbGTFfile ref/gencode.v45.annotation.gtf \
  --sjdbOverhang 100 --runThreadN 8

STAR --genomeDir ref/star_idx \
  --readFilesIn trim_R1.fastq.gz trim_R2.fastq.gz \
  --readFilesCommand zcat \
  --twopassMode Basic \
  --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts \
  --outFilterMultimapNmax 20 \
  --outFilterMismatchNoverLmax 0.04 \
  --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 \
  --runThreadN 8 --outFileNamePrefix star/sample01.

The index build needs about 32 GB of RAM. --quantMode GeneCounts writes a per-gene count file with columns for unstranded and both stranded conventions, which is a useful independent check on the Salmon strandedness inference. If the column matching your protocol has three to ten times the counts of the other stranded column, your strandedness call is right.

6. Collapse transcripts to genes

Do this in R with tximport, and use the lengthScaledTPM counts rather than raw estimated counts. That option scales counts by average transcript length across samples, which corrects for the fact that differential isoform usage changes a gene’s effective length even when its total output does not.

library(tximport); library(readr)

tx2gene <- read_tsv("ref/tx2gene.v45.tsv",
                    col_names = c("tx_id", "gene_id", "gene_name"))

files <- file.path("quant", list.files("quant"), "quant.sf")
names(files) <- list.files("quant")

txi <- tximport(files, type = "salmon", tx2gene = tx2gene[,1:2],
                countsFromAbundance = "lengthScaledTPM")

dim(txi$counts)
write.csv(txi$counts, "gene_counts_lengthScaledTPM.csv")

Build tx2gene.v45.tsv straight from the GTF:

zcat ref/gencode.v45.annotation.gtf.gz \
  | awk '$3=="transcript"' \
  | sed -n 's/.*transcript_id "\([^"]*\)".*gene_id "\([^"]*\)".*gene_name "\([^"]*\)".*/\1\t\2\t\3/p' \
  > ref/tx2gene.v45.tsv

You should end up with roughly 62,000 gene IDs, of which about 20,000 are protein-coding. In whole blood at 50 million reads, expect 12,000 to 15,000 genes with a count above 10.

7. Check the QC before you interpret anything

Aggregate everything with MultiQC and read it with specific thresholds in mind.

multiqc quant/ star/ . -o qc/

The numbers we check, in order. Salmon mapping rate should be above 80 percent for poly-A blood RNA. A rate below 70 percent usually means rRNA or globin contamination, genomic DNA carryover, or the wrong annotation release. Duplication rate above 60 percent without UMIs suggests a low-complexity library where your effective depth is much lower than your read count. Ribosomal RNA fraction above 5 percent in a poly-A library means the capture underperformed. And from fastp, a mean insert size below 120 nucleotides points to RNA degradation during collection or storage, which biases quantification toward short transcripts and toward 3’ ends.

Degradation is the failure mode most worth catching, because it produces plausible-looking numbers. If your provider reports an RNA integrity number, treat values below 7 as a caution and below 5 as grounds for re-collecting rather than reinterpreting. A degraded sample in a longitudinal series will look like a real biological shift across hundreds of genes.

8. Normalize for the question you are asking

TPM (transcripts per million) is appropriate for comparing relative abundance of different genes within one sample, and for eyeballing whether a gene is expressed at all. It is the wrong unit for comparing one gene across samples or timepoints, because a small number of very high-abundance transcripts (globin, immunoglobulin, mitochondrial) can shift the per-million denominator and move every other gene in the opposite direction.

For across-sample comparison, put raw counts into DESeq2 and use its median-of-ratios size factors, or use edgeR’s TMM. With a single-subject longitudinal series you will not have biological replicates in the usual sense, so treat DESeq2 output as a ranking and an effect-size estimate rather than a hypothesis test. What we find more informative is a per-gene reference interval built from your own prior timepoints: once you have five or more, a gene’s current value expressed as a z-score against your own history is a more meaningful signal than any single-timepoint comparison to a population cohort.

9. Interpret transcriptome data alongside your other layers

RNA sequencing tells you which transcripts were present in the cells you sampled at the moment you sampled them, which is a dynamic and tissue-specific readout, unlike your genome. That makes it most useful when read against other layers rather than alone. A coding variant in your whole-genome data becomes interpretable when RNA-seq shows whether the affected allele is expressed at all, and expression changes become mechanistically legible when protein measurements and metabolic data point the same direction. Frameworks for combining genomic, transcriptomic, and proteomic layers at the systems level are well developed, and the recurring lesson is that concordance across layers is the evidence that matters, while single-layer signals are hypotheses.³ Reproducible integration is enough of an engineering problem that workflow managers like Nextflow are now standard for wiring the pieces together with pinned containers.⁴

Two structural limits are worth holding in mind. First, blood RNA is a mixture of cell types. An apparent expression change often reflects a shift in cell composition (more neutrophils, fewer lymphocytes) rather than a change in transcription within any cell. Deconvolution tools such as CIBERSORTx estimate those fractions from bulk data and are worth running before you attribute a change to regulation. Single-cell and single-cell multi-omics methods resolve this properly by measuring genotype and phenotype in the same cell, which is why they have become the tool of choice for heterogeneous tissue.⁵ Second, standard poly-A RNA-seq is a narrow window on the transcriptome. Non-polyadenylated species, antisense transcripts, and small RNAs require different chemistry, and the catalog of functional RNA outside annotated mRNA is still expanding.⁶

Common problems

Low mapping rate with high-quality reads almost always means a reference mismatch. Check that your transcriptome FASTA and your tx2gene table come from the same GENCODE release, and check that you used --gencode at index time so transcript IDs match between the index and the mapping file. A version mismatch silently drops a few thousand transcripts.

Strandedness set wrong is the error that produces the most confidently wrong results. Quantifying a stranded library as ISF when it is ISR will cut mapping rates and distort antisense-overlapping genes badly. Always run -l A on the first sample of any new batch and read lib_format_counts.json.

Counts that differ between Salmon and featureCounts by more than about 20 percent for a specific gene usually point to a multimapping paralog family (HLA genes, immunoglobulin loci, ribosomal protein genes) where featureCounts discards ambiguous reads and Salmon distributes them. Neither is wrong. For those gene families, look at the Gibbs posterior spread before you believe any single estimate.

Version skew across timepoints will fake biology. If you reprocess timepoint 1 with GENCODE v41 and timepoint 5 with v45, hundreds of genes change annotation boundaries and their counts move. Reprocess the entire series with one pinned pipeline whenever you change anything, and record the container digest alongside the results.

Finally, a low-input or heavily degraded sample can be technically salvageable and biologically misleading at the same time. The pipeline will run, the report will look acceptable, and the expression profile will be shifted toward short and 3’-proximal transcripts. Check insert size and gene body coverage uniformity before trusting a timepoint, and if a result would change what you do about your health, bring it and the raw data to a physician rather than acting on it.

Footnotes

1. Andrea Orioli. tRNA biology in the omics era: Stress signalling dynamics and cancer progression. BioEssays, 2016. https://doi.org/10.1002/bies.201600158 ↩

2. Scott Kuersten, Agnes Radek, Christine Vogel, et al. Translation regulation gets its ‘omics’ moment. WIREs RNA, 2013. https://doi.org/10.1002/wrna.1173 ↩

3. Gun Tae Jung, Kwang-Pyo Kim, Kwoneel Kim. How to interpret and integrate multi-omics data at systems level. Animal Cells and Systems, 2020. https://doi.org/10.1080/19768354.2020.1721321 ↩

4. Bianka Alexandra Pasat, Eleftherios Pilalis, Katarzyna Mnich, et al. MultiOmicsIntegrator: a nextflow pipeline for integrated omics analyses. Bioinformatics Advances, 2024. https://doi.org/10.1093/bioadv/vbae175 ↩

5. Lei Yu, Xinlei Wang, Quanhua Mu, et al. scONE-seq: A single-cell multi-omics method enables simultaneous dissection of phenotype and genotype heterogeneity from frozen tumors. Science Advances, 2023. https://doi.org/10.1126/sciadv.abp8901 ↩

6. Marc Güell, Eva Yus, Maria Lluch-Senar, et al. Bacterial transcriptomics: what is beyond the RNA horiz-ome?. Nature Reviews Microbiology, 2011. https://doi.org/10.1038/nrmicro2620 ↩