How to Analyze Your Own RNA-Seq Data

By the end of this guide you will have taken a set of paired-end FASTQ files from a single bulk RNA sequencing run on your own blood and produced three things: a gene-by-sample matrix of estimated counts, a matching matrix of transcript per million (TPM) values with quality-control metrics attached, and a ranked list of the genes and pathways that separate one of your timepoints from another. You need a Linux or macOS machine with at least 16 GB of RAM, though 32 GB is more comfortable. You also need roughly 100 GB of free disk per sample if you keep intermediates, a conda or mamba environment, and R 4.3 or later. The whole pipeline for one sample runs in well under an hour on a modern laptop if you use a lightweight quantifier, which is what we recommend. This is a measurement and interpretation guide. Nothing here is a diagnosis, and any result that looks clinically consequential belongs in front of a physician who can order confirmatory testing.

1. Understand what the transcriptome measurement is

Transcriptome analysis is the quantification of RNA molecules present in a sample at the moment it was collected. Where your genome is effectively fixed, the transcriptome is a snapshot of which genes were being transcribed in the specific cells you sampled, and in what relative abundance. This is why transcriptomics sits between genomics and proteomics in the standard omics stack: the genome describes capacity, the transcriptome describes current program, and the proteome and metabolome describe the executed state ¹. Each layer answers a different question, and the correlation between them is weaker than people expect, which is the main argument for measuring more than one ².

Two practical consequences follow. First, tissue matters enormously. A blood draw gives you the transcriptome of peripheral leukocytes, which is informative about immune state, inflammation, and interferon signaling, and largely silent about liver or muscle. Second, condition at collection matters. Time of day, recent exercise, a cold coming on, and the time between draw and stabilization all move hundreds of genes. Longitudinal designs are what rescue transcriptomics from being noise: comparing you to you, under known conditions, is far more interpretable than comparing you to a reference cohort ³.

2. Set up the environment and get the reference

Install the tools in one isolated environment. We use mamba because conda’s solver is slow on bioinformatics channels.

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

For quantification we recommend Salmon against a decoy-aware transcriptome index rather than a full spliced alignment with STAR. Salmon’s selective alignment gives accurate transcript-level estimates in a few minutes per sample with about 4 GB of RAM, and the decoy sequences suppress the main failure mode of lightweight mappers, which is spurious assignment of reads that truly originate from unannotated genomic regions. Use STAR only if you need genomic coordinates: novel splice junctions, variant calling from RNA, or fusion detection.

Download a matched transcriptome and genome. Keep the release version in the filename, because mixing annotation versions between timepoints silently changes your gene IDs and breaks comparability.

REL=111
wget https://ftp.ensembl.org/pub/release-${REL}/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz
wget https://ftp.ensembl.org/pub/release-${REL}/fasta/homo_sapiens/ncrna/Homo_sapiens.GRCh38.ncrna.fa.gz
wget https://ftp.ensembl.org/pub/release-${REL}/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
wget https://ftp.ensembl.org/pub/release-${REL}/gtf/homo_sapiens/Homo_sapiens.GRCh38.${REL}.gtf.gz

Build the decoy-aware index. The gentrome is the transcriptome concatenated with the whole genome, and the decoy list tells Salmon which of those sequences are genomic background.

cat Homo_sapiens.GRCh38.cdna.all.fa.gz Homo_sapiens.GRCh38.ncrna.fa.gz \
    Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz > gentrome.fa.gz

grep "^>" <(gunzip -c Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz) \
  | cut -d " " -f 1 | sed 's/>//g' > decoys.txt

salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_idx_r${REL} \
  -k 31 --gencode=false -p 8

Indexing needs roughly 16 GB of RAM and half an hour. Keep -k 31 unless your reads are shorter than 50 bp, in which case drop to -k 23.

3. Quality control and trim the raw reads

Run FastQC on the raw FASTQs first, then trim with fastp. Looking at the untrimmed data is what tells you whether a problem came from the library or from your own trimming choices.

fastqc -t 8 -o qc_raw/ sample_R1.fastq.gz sample_R2.fastq.gz

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 \
  --cut_front --cut_tail --cut_mean_quality 20 \
  --length_required 36 \
  --trim_poly_g \
  --json fastp.json --html fastp.report.html --thread 8

--trim_poly_g matters on NovaSeq and other two-color chemistries, where dark cycles at the end of a read are miscalled as G and inflate mismatch rates. --detect_adapter_for_pe is more reliable than guessing the adapter sequence. Resist aggressive quality trimming beyond this: modern quantifiers handle moderate base-quality decay better than they handle length-truncated reads.

Two FastQC panels are worth reading carefully rather than glancing at. Per-base sequence content that is wild for the first 10 to 12 bases is normal for random-hexamer-primed libraries and is not a reason to trim. Sequence duplication levels above roughly 60 percent in a polyA library, combined with a small number of sequences accounting for a large fraction of reads, usually means low library complexity from a degraded or low-input sample. No amount of computation recovers it.

4. Quantify with Salmon

Run Salmon once per sample. Use automatic library type detection on the first sample, confirm the inferred type in lib_format_counts.json, then hard-code it for the rest so that a mislabeled sample does not quietly get a different setting.

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

A few of these flags carry real weight. --gcBias corrects fragment-level GC content bias introduced by PCR, and it changes fold-change estimates noticeably for GC-extreme genes. --seqBias corrects the random-hexamer priming bias visible in FastQC. --numGibbsSamples 20 gives you posterior samples of the transcript abundances, which is how you find out whether a transcript-level estimate is well determined or is really just one of several indistinguishable isoforms. For most stranded Illumina kits the inferred library type is ISR.

Check the mapping rate in quant/sample01/logs/salmon_quant.log. For a healthy polyA blood library against a decoy-aware human index, expect 75 to 90 percent. Below 60 percent, stop and diagnose before you interpret anything, because the missing reads are rarely missing at random.

5. Import into R and aggregate to gene level

Use tximeta, which reads Salmon’s index checksum and attaches the correct annotation automatically. This removes the most common quiet error in transcriptomics, which is pairing quantifications from one annotation release with gene metadata from another.

library(tximeta); library(SummarizedExperiment); library(DESeq2)

coldata <- data.frame(
  files = file.path("quant", c("sample01","sample02","sample03","sample04"), "quant.sf"),
  names = c("sample01","sample02","sample03","sample04"),
  timepoint = factor(c("t0","t0","t1","t1"), levels = c("t0","t1")),
  stringsAsFactors = FALSE
)
stopifnot(all(file.exists(coldata$files)))

se  <- tximeta(coldata)            # transcript level
gse <- summarizeToGene(se)         # gene level, length-corrected

summarizeToGene does more than sum counts: it computes an average transcript length per gene per sample and stores it, so that a shift in isoform usage does not masquerade as a change in gene expression. DESeq2 uses those lengths as offsets when you build the dataset from gse, which is the main reason to prefer this route over naive summation of TPMs.

Then filter and normalize:

dds <- DESeqDataSet(gse, design = ~ timepoint)
keep <- rowSums(counts(dds) >= 10) >= 2     # >=10 counts in >=2 samples
dds  <- dds[keep, ]
dds  <- DESeq(dds)

Keep raw counts as the input and let DESeq2 handle normalization with its median-of-ratios size factors. Do not feed TPM into a count-based differential expression model. TPM is for comparing genes within a sample, counts plus size factors are for comparing a gene across samples, and swapping them is a routine source of wrong answers.

6. Quality control the quantified matrix

This step is where most real problems surface, and it matters more than any downstream statistic. Start with a variance-stabilizing transform and look at the sample relationships.

vsd <- vst(dds, blind = TRUE)
plotPCA(vsd, intgroup = "timepoint")

# hierarchical clustering on sample distances
d <- dist(t(assay(vsd)))
plot(hclust(d))

# globin and mitochondrial fraction: the two dominant blood artifacts
cpm <- sweep(counts(dds, normalized = TRUE), 2, colSums(counts(dds, normalized = TRUE)), "/")
globin <- c("HBB","HBA1","HBA2","HBD")
colSums(cpm[rownames(cpm) %in% globin, , drop = FALSE]) * 100

In whole blood, globin transcripts can consume 20 to 70 percent of reads unless the library used globin depletion or a PAXgene-style protocol. High globin fraction reduces effective depth and compresses your dynamic range for everything else. Report it per sample, and treat it as a covariate rather than something to ignore. Likewise, a rising mitochondrial read fraction usually tracks handling stress and delayed stabilization.

If PC1 separates your samples by collection date, processing batch, or globin fraction rather than by the condition you care about, you have a batch effect and not biology. Add the covariate to the design (~ globin_bin + timepoint) rather than removing it from the data before testing. Statistical adjustment inside the model preserves the correct standard errors. Visual removal with limma::removeBatchEffect is fine for plots and wrong as an input to testing.

7. Differential expression and pathway interpretation

With a clean matrix, extract shrunken effect sizes. The shrinkage step matters because unshrunken log fold changes for low-count genes are numerically enormous and biologically meaningless.

res <- lfcShrink(dds, coef = "timepoint_t1_vs_t0", type = "apeglm")
res <- res[order(res$padj), ]
sum(res$padj < 0.05 & abs(res$log2FoldChange) > 0.5, na.rm = TRUE)

With n = 2 per group, and especially with n = 1, the p-values are not trustworthy and should be read as a ranking rather than as inference. The central limitation of personal transcriptomics is that you cannot estimate within-condition variance from a single observation. The way out is repetition. Five or six baseline draws under matched conditions, meaning the same time of day and the same fasting state, give you a personal reference distribution. A new sample can then be scored against your own historical variance for each gene, which is a far stronger comparison than a population reference.

For interpretation, prefer a rank-based enrichment method over a cutoff-based one, because it uses the whole ranked list and does not depend on an arbitrary significance threshold.

library(fgsea); library(msigdbr)

stats <- res$stat; names(stats) <- rownames(res)
stats <- stats[!is.na(stats)]

paths <- msigdbr(species = "Homo sapiens", collection = "H")
paths <- split(paths$ensembl_gene, paths$gs_name)

fg <- fgsea(paths, stats, minSize = 15, maxSize = 500, nPermSimple = 10000)
head(fg[order(fg$padj), c("pathway","padj","NES","size")], 15)

The MSigDB hallmark collection is the right starting point because its 50 gene sets are curated to be non-redundant, so you get a readable summary instead of 4,000 overlapping GO terms. In blood, expect the interferon alpha and gamma response, inflammatory response, and TNFA signaling hallmarks to be the most responsive to short-term perturbation. Deconvolution of cell-type proportions (CIBERSORTx, or a marker-based approach) is worth running alongside, since a large share of blood expression change is a shift in cell composition rather than a change in per-cell regulation.

Transcriptome results become considerably more interpretable when read against other layers. Pairing an expression change with a genotype that affects the same gene, or with a protein measurement of its product, distinguishes a regulatory signal from noise. The machine learning and integration strategies for combining these layers are now reasonably well characterized ⁴. The caveat is that transcript abundance predicts protein abundance only moderately, so treat a transcript-only finding as a hypothesis about the protein rather than a measurement of it ². Single-cell methods sharpen this further by resolving which cell population carries a signal, at the cost of more complex informatics ⁵.

Common problems

Low mapping rate, in the 30 to 60 percent range, almost always has one of four causes. Check them in this order. First, wrong species or a contaminated sample, which you can test by running a few thousand reads through a taxonomic classifier. Second, residual ribosomal RNA from a failed depletion, which you can test by mapping against an rRNA reference and counting. Third, untrimmed adapter from a short insert library. Fourth, genomic DNA contamination, visible as reads distributed across introns and requiring a STAR alignment to see. An index built from a transcriptome without decoys inflates the apparent rate while degrading accuracy, so do not “fix” a low rate by removing decoys.

Gene identifier drift causes quiet, hard-to-find breakage. Ensembl IDs with version suffixes (ENSG00000141510.17) will not join against unversioned IDs, and gene symbols are unstable across annotation releases. Store Ensembl gene IDs without versions as your primary key, carry symbols as a display label only, and pin one annotation release for the whole longitudinal series.

High duplication rates get misread constantly. In RNA-seq, duplicate reads are expected because highly expressed transcripts genuinely produce many identical fragments. Do not deduplicate bulk RNA-seq without unique molecular identifiers, because you will preferentially discard signal from your most abundant genes. If the library has UMIs, deduplicate on the UMI. Do not deduplicate on coordinates.

Batch structure that correlates perfectly with your condition of interest is unfixable after the fact. If every t0 sample went through one kit lot and every t1 sample through another, no model separates kit from time. Randomize processing order across conditions when you can, and record lot numbers, draw time, fasting state, recent exercise, and any acute illness for every sample. In a personal longitudinal study this metadata is the difference between a usable series and an uninterpretable one, and the general point holds across all the omics layers ³.

Finally, resist over-reading a single extreme gene. A single-sample outlier at one gene is usually a mapping artifact, a small number of multi-mapped reads in a paralog family, or a genotype in your own genome that disrupts the reference k-mers for that transcript. Require that a finding replicate across timepoints, appear as a coherent pathway-level shift, or corroborate against a second measurement type before you take it seriously. Where a result touches on something with clinical weight, the correct next step is a conversation with a physician and an orthogonal validated assay. It is not a deeper dive into the same FASTQs.

Footnotes

1. Richard P Horgan, Louise C Kenny. ‘Omic’ technologies: genomics, transcriptomics, proteomics and metabolomics. The Obstetrician & Gynaecologist, 2011. https://doi.org/10.1576/toag.13.3.189.27672 ↩

2. Magdalena Zapalska-Sozoniuk, Lukasz Chrobak, Krzysztof Kowalczyk, et al. Is it useful to use several “omics” for obtaining valuable results?. Molecular Biology Reports, 2019. https://doi.org/10.1007/s11033-019-04793-9 ↩ ↩²

3. Michela Gottardi Zamperla, Veronica Barbi, Sara Negri, et al. Omics medicine: what the clinicians should know. European Journal of Internal Medicine, 2026. https://doi.org/10.1016/j.ejim.2026.106759 ↩ ↩²

4. Pedro H. Godoy Sanches, Nicolly Clemente de Melo, Andreia M. Porcari, et al. Integrating Molecular Perspectives: Strategies for Comprehensive Multi-Omics Integrative Data Analysis and Machine Learning Applications in Transcriptomics, Proteomics, and Metabolomics. Biology, 2024. https://doi.org/10.3390/biology13110848 ↩

5. Luonan Chen, Rong Fan, Fuchou Tang. Advanced Single-Cell Omics Technologies and Informatics Tools for Genomics, Proteomics, and Bioinformatics Analysis. Genomics, Proteomics & Bioinformatics, 2021. https://doi.org/10.1016/j.gpb.2021.12.001 ↩