How to Run Your Own Bulk RNA-Seq Analysis Well

By the end of this guide you will have a gene-level expression matrix for your own samples, a quantification QC report that tells you whether to trust it, a cell-composition estimate for whole blood, and a per-gene longitudinal view of how your own transcriptome moves across draws. You need the raw FASTQ files and the library preparation metadata, a machine with 32 GB of RAM and roughly 100 GB of free disk per dozen samples, and a package manager (we use micromamba). Everything below assumes bulk RNA-seq of whole blood or PBMCs from a single person sampled repeatedly, which is the common case for a personal profile and which has failure modes that generic tutorials written for case-control experiments do not cover.

1. Read the library metadata before you run anything

The single largest determinant of what your data can answer is a set of choices made at the bench, and you cannot recover from most of them computationally. Find out five things.

• whether RNA was captured in a stabilizing tube (PAXgene or Tempus) or from a delayed-processing EDTA tube
• whether the library is poly-A selected or ribosomal-depleted
• whether globin mRNA was depleted
• whether the library is stranded (almost always dUTP reverse-stranded on Illumina today)
• the read configuration and depth

Whole blood is the hard case. Reticulocytes carry enormous amounts of hemoglobin transcript, and without globin depletion HBB, HBA1, HBA2 and HBD can consume the majority of your reads, leaving the rest of the transcriptome badly undersampled. Poly-A selection also requires intact RNA: a RIN (RNA Integrity Number) below about 7 produces 3’ coverage bias that inflates short transcripts and depresses long ones. If you have a choice at collection time, PAXgene tubes plus ribo-depletion plus globin depletion gives you the widest usable dynamic range, at the cost of more reads spent on intronic and non-coding RNA.

Depth targets follow from the question. Around 30 million paired-end reads per sample is adequate for gene-level quantification of moderately expressed genes. Alternative splicing analysis needs considerably more. Its accuracy depends jointly on depth, read length and the quantifier you choose, which is exactly what large reference-material benchmarks have been built to measure ¹. Record all of this in a samples.tsv now, because it becomes the covariate table later.

2. Build the environment and a decoy-aware index

Reproducibility in transcriptomics starts with pinning the reference and the tool versions, since a gene-level count is not comparable across annotation releases. Aging and longitudinal omics work in particular has shown how much apparent biological signal is really version and batch drift. Version pinning is therefore part of the measurement and should not be treated as bookkeeping ².

micromamba create -n rna -c conda-forge -c bioconda \
  salmon=1.10.3 fastp=0.23.4 fastqc=0.12.1 multiqc=1.21 \
  star=2.7.11b samtools=1.19 -y
micromamba activate rna

Use the GRCh38 no-alt analysis set for the genome and a single GENCODE release for the annotation (we use GENCODE v46 and do not mix). Build a decoy-aware Salmon index so that reads originating from unannotated genomic sequence are absorbed by the genome rather than forced onto a transcript:

# transcriptome + genome as decoys
grep "^>" GRCh38.primary_assembly.genome.fa | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v46.transcripts.fa GRCh38.primary_assembly.genome.fa > gentrome.fa

salmon index -t gentrome.fa -d decoys.txt -i salmon_idx_v46 \
  -k 31 --gencode -p 16

Keep -k 31 unless your reads are shorter than 50 bp, in which case drop to 23 or 25. The --gencode flag strips the pipe-delimited GENCODE FASTA headers down to the transcript ID, which saves you a downstream string-parsing bug. Save salmon_idx_v46 and the decoys.txt checksum alongside your results.

3. Trim, QC, and look at the numbers rather than the traffic lights

Read-level QC exists to catch adapter carryover, low-quality tails, and contamination before they propagate into quantification. We run fastp because it does adapter detection, quality trimming and reporting in one pass, and then aggregate with MultiQC.

fastp -i ${S}_R1.fastq.gz -I ${S}_R2.fastq.gz \
  -o trim/${S}_R1.fq.gz -O trim/${S}_R2.fq.gz \
  --detect_adapter_for_pe --trim_poly_g --cut_tail --cut_tail_mean_quality 20 \
  --length_required 36 --thread 8 \
  --json qc/${S}.fastp.json --html qc/${S}.fastp.html

--trim_poly_g matters on NovaSeq and NextSeq two-color chemistry, where the absence of signal is read as G and dark cycles produce long poly-G tails. Do not trim aggressively beyond this. Selective alignment in Salmon handles soft-clippable ends well, and over-trimming shortens reads into the range where multi-mapping between paralogs rises.

The FastQC modules that warn most loudly are usually the least informative for RNA-seq. Duplicate-rate warnings are expected for highly expressed transcripts. Per-base sequence content skew in the first 12 bases reflects random-hexamer priming bias and does not indicate a defect. The modules to take seriously are overrepresented sequences (adapter dimer or rRNA), per-sequence GC content with a second mode (contamination), and a sharp drop in per-tile quality.

4. Quantify with Salmon, and know when to align instead

For gene-level expression, selective alignment against a decoy-aware index is our default. It is fast enough to rerun the whole cohort when the annotation changes, it models fragment-level GC and positional bias, and it gives you a bootstrap estimate of inferential uncertainty for genes with ambiguous multi-mapping.

salmon quant -i salmon_idx_v46 -l A \
  -1 trim/${S}_R1.fq.gz -2 trim/${S}_R2.fq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  --numGibbsSamples 20 --threads 16 \
  -o quant/${S}

Use -l A to let Salmon infer library type, then check quant/${S}/lib_format_counts.json and confirm that the inferred type is consistent across all your samples. A sample that comes back IU when the rest are ISR is a sample-swap or a prep failure. It is not a quirk. Turn on --gcBias always for multi-batch data: fragment GC bias is one of the most common sources of apparent differential expression between library batches.

Genome alignment with STAR is still required for three things: splicing analysis, allele-specific expression, and any visual inspection of coverage in IGV. If you plan to do those, run STAR in addition rather than instead, with --outSAMattributes NH HI AS nM NM MD and --quantMode GeneCounts. Budget about 32 GB of RAM for a human index.

5. Import to gene level and filter once

Transcript-level estimates are the natural output of Salmon, but gene-level counts are more stable and are what most downstream tooling expects. tximport aggregates transcripts to genes while carrying an average transcript-length offset, which corrects for the fact that isoform switching changes the effective length of a gene between samples.

library(tximport); library(DESeq2)
tx2gene <- read.delim("gencode.v46.tx2gene.tsv")  # tx_id, gene_id
files <- file.path("quant", samples$id, "quant.sf"); names(files) <- samples$id
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ batch + timepoint)
keep <- rowSums(counts(dds) >= 10) >= ceiling(0.25 * ncol(dds))
dds <- dds[keep, ]
vsd <- vst(dds, blind = FALSE)

Use lengthScaledTPM when you want counts that are already length-corrected and can be fed to limma-voom, and the default "no" with the length offset when you stay inside DESeq2. Filter once, before normalization, and use a presence threshold rather than a mean threshold so that a gene expressed in a third of your timepoints is not deleted by a long tail of zeros. In whole blood, expect roughly 13,000 to 16,000 genes to survive a filter like this.

6. QC the quantification and do not stop at the reads

This is the step most tutorials skip, and it is where personal longitudinal data goes wrong. Five checks, in order.

First, the mapping rate reported in quant/${S}/logs/salmon_quant.log. For poly-A whole blood on a decoy-aware index, expect roughly 70 to 90 percent. A drop of fifteen points in one sample means degraded RNA or contamination. It does not point to a modelling issue. Second, globin fraction. Sum the TPM of the following genes and record it per sample.

• HBB
• HBA1
• HBA2
• HBD
• HBG1
• HBG2

If it swings from 5 percent to 40 percent across your timepoints, your effective depth is swinging with it and every downstream comparison inherits that. Third, rRNA and mitochondrial fraction, which flag failed depletion and cell stress respectively. Fourth, a sex check using XIST and RPS4Y1, which catches sample swaps immediately. Fifth, genotype concordance: call common SNPs from the STAR BAM with bcftools mpileup restricted to a panel of high-coverage exonic sites and confirm that every sample matches your own genome. For a single-person profile, a swapped tube is otherwise invisible.

Plot a PCA of the VST matrix and colour by batch, by globin fraction, and by draw date. If PC1 tracks globin fraction rather than time, fix the normalization before interpreting anything. Systematic reporting of these metrics is what separates a reusable dataset from a one-off, and it is the same discipline that omics-based risk assessment frameworks have converged on ³.

7. Decide what “high” means for one person

With a single individual there is no case-control contrast, so the comparisons available to you are internal (this draw versus your own baseline) and external (you versus a population reference). Both are useful and both mislead in different ways.

For the internal comparison, compute a per-gene z-score against your own prior draws on the VST scale, requiring at least four or five baseline timepoints before the variance estimate means much. This directly answers the question you care about: which transcripts moved, relative to how much they normally move in you. Pair it with a moderated model in limma, using duplicateCorrelation to account for repeated measures if you ever add a second individual:

library(limma)
v <- voom(counts(dds), model.matrix(~ batch + timepoint, colData(dds)), plot = TRUE)
fit <- eBayes(lmFit(v, model.matrix(~ batch + timepoint, colData(dds))))
topTable(fit, coef = "timepoint", number = 50)

For the external comparison, GTEx whole blood is the usual reference, and it is a rough one. GTEx blood samples were collected postmortem under PAXgene, quantified against an older annotation, and prepared with a different protocol from yours. Compare ranks rather than absolute TPM, restrict to genes with median TPM above 1 in the reference, and treat anything within a factor of two as unresolved. Interpretation of an individual transcriptome against a population distribution is a research activity. It is not a clinical one. If a finding would change what you do about your health, that conversation belongs with a physician or a clinical geneticist working from a validated assay.

8. Estimate cell composition before interpreting any blood gene

Bulk blood RNA is a weighted average over neutrophils, lymphocytes, monocytes and platelets, and the weights move day to day. A two-fold change in a monocyte-restricted transcript usually means you had more monocytes that morning. Run a deconvolution and carry the estimated fractions as covariates rather than trying to interpret raw expression.

We use CIBERSORTx with the LM22 signature matrix on TPM (not VST, not log), or immunedeconv in R which wraps several methods behind one interface. Compare the estimated neutrophil-to-lymphocyte ratio against the differential from a CBC drawn at the same time. If they disagree badly, trust the CBC and treat the deconvolution as unusable for that sample. Adding cell-fraction covariates to the design typically removes the majority of apparent longitudinal variance in blood, which is the point: the residual is what is worth looking at, and it becomes far more interpretable once other layers of data are aligned to the same timepoints ⁴.

9. Treat splicing and variant calls as a separate, harder project

Alternative splicing analysis from RNA-seq needs genome alignment, high depth, and a quantifier chosen with its biases in mind, because event-level estimates diverge substantially between methods on the same reads. Adding a flag to the pipeline is not enough. Benchmarking against reference materials with known ground truth shows that accuracy depends on library protocol, depth and the detection tool together, and that filtering thresholds materially change which events survive ¹. If you go ahead, run rMATS-turbo on STAR BAMs with --novelSS off initially. Require at least 10 junction reads supporting each event, and validate anything interesting by looking at the sashimi plot yourself.

Allele-specific expression is more tractable and more informative for a single person: with your whole genome already sequenced, you have phased heterozygous sites, and phASER or GATK ASEReadCounter will tell you where one allele is systematically silenced. Apply WASP filtering to remove reference-mapping bias before believing any imbalance.

10. Make the whole thing rerunnable

Wrap the pipeline in Snakemake or Nextflow, or use nf-core/rnaseq directly if you want a maintained implementation with the QC steps already wired in. Record the tool versions, the GENCODE release, the index checksum and the sample metadata in a single manifest that ships with the counts matrix. When you add a thirteenth timepoint next year, you will reprocess all thirteen from raw FASTQ with the same versions rather than merging a new matrix into an old one. Documented, versioned workflows are the practical difference between a longitudinal dataset and a pile of quantifications ⁵, and they are what makes the transcriptome usable alongside proteomic and genomic layers later ⁶.

Common problems

Mapping rate drops in one sample. Check the fastp JSON for adapter content and read length first, then the GC histogram for a contaminant mode. Degraded RNA shows up as a shorter insert size distribution and a rising 3’ coverage bias in the STAR output. It does not usually show up as a mapping failure.

PC1 of your VST matrix separates two batches perfectly. Do not regress it out and move on. Confirm the batches differ in library prep date or kit lot, include batch in the design, and if the batches are confounded with the biology you care about, the comparison is not recoverable. Batch structure is the dominant reproducibility problem in longitudinal omics ².

A gene shows a ten-fold jump in one draw. Before believing it, check whether it is a single-isoform gene with low counts (look at the Gibbs samples in Salmon’s output for inferential variance). Check whether it is cell-type restricted (check your deconvolution), and whether the sample’s globin fraction is an outlier. Most dramatic single-timepoint jumps in blood fail at least one of these.

Transcript-level results that disappear at gene level. This is expected when two isoforms of the same gene move in opposite directions. It is real biology and also the regime where quantifier disagreement is largest, so require both a large effect and consistent direction across bootstrap replicates.

TPM compared across studies. Do not. TPM is normalized within a sample and is not comparable across protocols, annotations or depths. Compare ranks, or reprocess the external data through your own pipeline.

Interpreting a pathway enrichment as a finding. Enrichment of an immune signature in whole blood is nearly always a composition effect. Rerun the enrichment on residuals after adjusting for cell fractions, and read any remaining signal as a hypothesis for the next draw rather than a conclusion. Integrating transcript-level signal with protein and clinical measurements taken at the same time is what turns a hypothesis into something with support ⁷.

Footnotes

1. Duo Wang, Jiaxin Zhao, Qingwang Chen, et al. Benchmarking RNA-seq with the Quartet and MAQC reference materials to establish best practices for accurate alternative splicing analysis. Nature Communications, 2026. https://doi.org/10.1038/s41467-026-76380-z ↩ ↩²

2. Param Priya Singh, Bérénice A. Benayoun. Considerations for reproducible omics in aging research. Nature Aging, 2023. https://doi.org/10.1038/s43587-023-00448-4 ↩ ↩²

3. Roland Buesen, Brian N. Chorley, Beatriz da Silva Lima, et al. Applying ‘omics technologies in chemicals risk assessment: Report of an ECETOC workshop. Regulatory Toxicology and Pharmacology, 2017. https://doi.org/10.1016/j.yrtph.2017.09.002 ↩

4. Yehudit Hasin, Marcus Seldin, Aldons Lusis. Multi-omics approaches to disease. Genome Biology, 2017. https://doi.org/10.1186/s13059-017-1215-1 ↩

5. Michele Araújo Pereira, Eddie Luidy Imada, Rafael Lucas Muniz Guedes. RNA‐seq: Applications and Best Practices. Applications of RNA-Seq and Omics Strategies - From Microorganisms to Human Health, 2017. https://doi.org/10.5772/intechopen.69250 ↩

6. Samuel M. Lancaster, Akshay Sanghi, Si Wu, et al. A Customizable Analysis Flow in Integrative Multi-Omics. Biomolecules, 2020. https://doi.org/10.3390/biom10121606 ↩

7. Tonci Ivanisevic, Raj N. Sewduth. Multi-Omics Integration for the Design of Novel Therapies and the Identification of Novel Biomarkers. Proteomes, 2023. https://doi.org/10.3390/proteomes11040034 ↩