How to Analyze Your Own Bulk RNA-Seq Data

Oak
rna-seq transcriptomics bioinformatics salmon deseq2
A tall lab machine of stacked glowing water channels with silver fish-like forms swimming upstream through counting gates.

By the end of this guide you will have taken raw paired-end RNA sequencing reads from your own blood and produced a gene-level expression matrix, a set of quality control reports you can defend, and a within-person comparison of expression across time points, with pathway-level summaries. RNA-seq is short for RNA sequencing: you convert the RNA present in a sample into cDNA, sequence it, and count how many reads map to each transcript, which gives a measurement of relative abundance for roughly 15,000 to 20,000 expressed genes at once. That is the difference from qPCR, which quantifies a handful of transcripts you chose in advance with better precision per target. The approach in its modern short-read form dates to the late 2000s and has since become the standard assay for transcript abundance across research and clinical discovery work 1. You need: paired FASTQ files (typically sample_R1.fastq.gz and sample_R2.fastq.gz), a machine with 16 GB of RAM and about 100 GB of free disk, a POSIX shell, conda or mamba, and R 4.3 or later. A laptop is enough for quantification with Salmon. Genome alignment with STAR is not, unless you have 32 GB or more.

1. Inventory what you received and read the metadata first

Before running anything, confirm the shape of the data. For each sample you should have two gzipped FASTQ files, a checksum file, and a metadata record giving library type, collection date, tissue, and library preparation chemistry. The single most consequential metadata field is whether the library was poly(A)-selected or ribo-depleted, because it determines which RNA species you can quantify and how you interpret intronic reads.

ls -lh fastq/
md5sum -c fastq/checksums.md5
zcat fastq/T01_R1.fastq.gz | head -4
zcat fastq/T01_R1.fastq.gz | echo $(( $(wc -l) / 4 )) reads

For bulk blood RNA-seq intended for differential expression, 20 to 30 million paired reads per sample is the usual working range. Below about 10 million you lose power on low-expressed transcripts. Read length matters less than you would expect: 2x100 is comfortable, and 2x50 is adequate for gene-level counting.

Note the tissue. Blood transcriptomes are dominated by the cell composition of the sample, so a change in neutrophil fraction will move hundreds of genes without any change in per-cell regulation. Keep a complete blood count from the same draw if you have one.

2. Build the environment

Pin versions. Reproducing a result six months later is much harder if your tool versions float.

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 pigz
mamba activate rnaseq
salmon --version

Download the reference. We use GENCODE human release 45, primary assembly, with the matching transcript FASTA. Using a single annotation release for every sample you will ever compare is more important than which release you pick.

wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_45/gencode.v45.transcripts.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_45/GRCh38.primary_assembly.genome.fa.gz

3. Trim adapters and run read-level quality control

fastp handles adapter removal, quality trimming, and reporting in one pass and is fast enough that there is no reason to skip it. Poly-G trimming matters on two-colour Illumina chemistry (NovaSeq, NextSeq), where a dark cycle at the end of a read is called as G.

fastp \
  -i fastq/T01_R1.fastq.gz -I fastq/T01_R2.fastq.gz \
  -o trim/T01_R1.fq.gz -O trim/T01_R2.fq.gz \
  --detect_adapter_for_pe \
  --trim_poly_g \
  --qualified_quality_phred 15 --unqualified_percent_limit 40 \
  --length_required 36 \
  --thread 8 \
  --json qc/T01.fastp.json --html qc/T01.fastp.html

Then aggregate everything with MultiQC and read the report before proceeding.

fastqc -t 8 -o qc/ trim/*.fq.gz
multiqc -o qc/ qc/

What we look at, in order: the fraction of reads surviving trimming (expect above 90 percent), duplication rate, GC distribution (a sharp secondary peak often means rRNA or adapter dimer), and per-sample read counts. High duplication in RNA-seq is not automatically a problem, because highly expressed transcripts genuinely produce identical fragments, but a duplication rate that differs sharply between your samples is a signal that library complexity differs and that your comparison is partly a comparison of library quality.

4. Quantify with Salmon using a decoy-aware index

We recommend Salmon’s selective alignment against a decoy-aware transcriptome over full genome alignment for personal bulk RNA-seq. It runs in minutes on a laptop, corrects for fragment GC and sequence bias, and gives transcript-level estimates that aggregate cleanly to genes. The tradeoff is that you get no BAM file, so you cannot inspect splice junctions, call variants from RNA, or detect novel transcripts. If you need any of those, run STAR in addition, not instead.

Build the index once. The decoy set is the genome, which prevents reads originating from unannotated genomic regions from being forced onto a transcript.

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

Then quantify each sample:

salmon quant -i salmon_idx_v45 -l A \
  -1 trim/T01_R1.fq.gz -2 trim/T01_R2.fq.gz \
  --validateMappings --gcBias --seqBias --posBias \
  --numBootstraps 30 \
  -p 8 -o quant/T01

-l A lets Salmon infer library strandedness; check quant/T01/lib_format_counts.json afterward and confirm the inferred type is consistent across samples (typically ISR for dUTP-based stranded kits). --numBootstraps 30 gives you inferential replicates, which matter if you ever move to transcript-level analysis where quantification uncertainty is large. Record the mapping rate from logs/salmon_quant.log. For poly(A) blood libraries we expect 75 to 90 percent. A mapping rate below 60 percent usually means contamination, wrong species, or heavy rRNA carryover.

5. Import to R and build a gene-level count matrix

tximport converts Salmon’s transcript estimates to gene-level counts while carrying an offset that accounts for average transcript length, which is the correct way to hand these numbers to DESeq2.

library(tximport); library(DESeq2); library(AnnotationDbi)
library(org.Hs.eg.db); library(readr)

tx2gene <- read_tsv("tx2gene.gencode.v45.tsv")  # tx_id, gene_id, gene_name
files <- file.path("quant", meta$sample, "quant.sf")
names(files) <- meta$sample

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                ignoreTxVersion = FALSE, countsFromAbundance = "no")
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ timepoint)

Filter before you model. Genes with almost no reads add noise and inflate multiple testing. A defensible rule is to keep genes with at least 10 counts in at least as many samples as your smallest group.

keep <- rowSums(counts(dds) >= 10) >= min(table(meta$timepoint))
dds <- dds[keep, ]
nrow(dds)   # expect roughly 13,000-16,000 for whole blood

6. Run sample sanity checks before any biology

This step catches the errors that silently ruin analyses. Do it every time.

First, a sex check. Plot normalized counts for XIST and RPS4Y1. They separate cleanly and will expose a sample swap or a mislabelled tube instantly.

vsd <- vst(dds, blind = TRUE)
plotPCA(vsd, intgroup = "timepoint")
mat <- assay(vsd)
cor_mat <- cor(mat, method = "spearman")
pheatmap::pheatmap(cor_mat)

For repeated samples from one person, the Spearman correlation between any two time points on variance-stabilized data is typically above 0.95. A sample sitting well below that is either degraded, differently processed, or not from the same person. Second, check globin. In whole blood collected without globin depletion, HBB, HBA1, and HBA2 can absorb a large share of reads, which compresses the dynamic range for everything else. Compute the globin fraction per sample and treat it as a covariate or a reason to exclude a sample if it varies widely. Third, check the mitochondrial fraction and the rRNA fraction as proxies for RNA integrity.

Finally, estimate cell composition. Bulk blood expression is a weighted average of cell types, and deconvolution methods (CIBERSORTx, or a simpler marker-based score) give you fractions you can include in the model. This is the point where single-cell methods are conceptually relevant: single-cell RNA-seq resolves the per-cell states that bulk data averages over, and the field has built extensive step-by-step workflows and interactive servers for that resolution 23. Bulk profiling with composition covariates recovers much of the interpretive value at a fraction of the cost and with far better longitudinal practicality 4.

7. Model expression across time points

With one person, you are not comparing cases to controls. You are comparing states of the same body: fasted versus fed, before versus after a training block, seasonal sampling, infection versus recovery. Two structures are common.

If you have several replicate draws per condition, a standard DESeq2 model works:

design(dds) <- ~ globin_frac + neutrophil_frac + condition
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "post", "pre"), alpha = 0.05)
res <- lfcShrink(dds, coef = "condition_post_vs_pre", type = "apeglm", res = res)
summary(res)

If you have a dense time series with a single sample per point, formal differential testing is not appropriate, because you have no estimate of within-condition variance. The useful analysis is different: build a reference distribution from your own historical samples, then express each new time point as a deviation from your personal baseline.

ref <- assay(vsd)[, meta$phase == "baseline"]
mu <- rowMeans(ref); sdv <- apply(ref, 1, sd)
z  <- (assay(vsd)[, "T09"] - mu) / pmax(sdv, 0.1)
head(sort(z, decreasing = TRUE), 30)

The pmax(sdv, 0.1) floor prevents genes with near-zero baseline variance from producing enormous z-scores off measurement noise. This personal-baseline framing is what makes longitudinal self-measurement worth doing: the comparison is against your own prior state rather than a population mean, which removes most of the genetic and compositional variation that dominates between-person studies.

8. Summarize at the pathway level

Individual gene results from a small number of samples are fragile. Gene set enrichment on the full ranked list is more stable because it aggregates weak, consistent signal across many genes.

library(fgsea); library(msigdbr)

hallmark <- msigdbr(species = "Homo sapiens", collection = "H")
sets <- split(hallmark$gene_symbol, hallmark$gs_name)

stats <- res$log2FoldChange; names(stats) <- res$gene_name
stats <- sort(stats[!is.na(stats)], decreasing = TRUE)

fg <- fgsea(sets, stats, minSize = 15, maxSize = 500, nPermSimple = 10000)
fg[order(padj)][1:15, .(pathway, NES, padj, size)]

Read the normalized enrichment score and the leading-edge genes together. An interferon or inflammatory response set moving in a blood sample is often driven by a handful of genes and by cell composition, so open the leading edge and check it is not three globin-adjacent artifacts. If you want to go further, transcriptomic signal is most interpretable when layered against other molecular measurements from the same person, and there are structured teaching workflows for integrating expression with epigenetic and other omics layers 51.

None of these outputs are diagnostic. A gene set that moves in your data describes a state of your blood transcriptome at one moment, not a disease, and expression values have no established clinical reference ranges for individuals. If something in your results concerns you, take it to a physician along with conventional clinical labs, which are the measurements clinical decisions are built on.

9. Freeze the analysis so future samples are comparable

Write out the count matrix, the session info, the tool versions, and the exact reference files used. Six months from now you will add another time point, and if the annotation release has changed your gene set will shift underneath you.

saveRDS(dds, "results/dds_v45_2026-09.rds")
write_tsv(as_tibble(counts(dds, normalized = TRUE), rownames = "gene_id"),
          "results/normalized_counts.tsv")
writeLines(capture.output(sessionInfo()), "results/sessionInfo.txt")

Common problems

Low mapping rate with high duplication usually means a low-input library that was over-amplified, and the fix is at the bench, not in software. Do not try to rescue it computationally beyond reporting it.

Batch effects that align with sequencing date are the most common false discovery in personal longitudinal data. If all your “before” samples were run in one batch and all your “after” samples in another, the experiment cannot distinguish biology from batch, no matter what you put in the model. Randomize or interleave when you can, and include the batch term when you cannot.

Strandedness inferred inconsistently across samples means different library chemistries got mixed. Check lib_format_counts.json for every sample and reprocess, rather than letting Salmon quietly apply different assumptions.

Enormous fold changes on lowly expressed genes are an artifact of dividing small numbers. Always use lfcShrink before ranking, and treat any gene with a base mean below about 20 as a hypothesis rather than a result.

Version drift between annotation releases silently changes gene identifiers and merges or splits loci. Store the reference FASTA checksums alongside your results.

Finally, cell composition. In whole blood, the single largest driver of variance across your own time points is usually the proportion of neutrophils and lymphocytes in the tube, which shifts with time of day, recent exercise, and acute illness. If you do not measure or estimate it, you will keep rediscovering it as biology.

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. Audrey Minden. Introduction to Omics Technologies for Research, Diagnosis, and Drug Discovery. Biochemistry, Cell and Molecular Biology for Pharmacists and Health Professionals, 2026. https://doi.org/10.1007/978-3-032-20529-2_27 ↩ ↩2

  2. Shaked Slovin, Annamaria Carissimo, Francesco Panariello, et al. Single-Cell RNA Sequencing Analysis: A Step-by-Step Overview. Methods in Molecular Biology, 2021. https://doi.org/10.1007/978-1-0716-1307-8_19 ↩

  3. Andrew Jiang, Klaus Lehnert, Linya You, et al. ICARUS, an interactive web server for single cell RNA-seq analysis. Nucleic Acids Research, 2022. https://doi.org/10.1093/nar/gkac322 ↩

  4. Sibo Zhu, Tao Qing, Yuanting Zheng, et al. Advances in single-cell RNA sequencing and its applications in cancer research. Oncotarget, 2017. https://doi.org/10.18632/oncotarget.17893 ↩

  5. Nathan A Ruprecht, Joshua D Kennedy, Benu Bansal, et al. Transcriptomics and epigenetic data integration learning module on Google Cloud. Briefings in Bioinformatics, 2024. https://doi.org/10.1093/bib/bbae352 ↩