How to Analyze Your Own Bulk RNA-seq Data

By the end of this you will have a gene-level count matrix and TPM matrix for every sample you own. You will have a QC report that tells you whether each library is usable, and an estimate of blood cell composition per sample. You will also have either a differential expression table (if you have a real contrast) or a per-gene trajectory across your own timepoints. You need the raw FASTQs (paired-end, gzipped, typically 4–12 GB per sample at 30–50M read pairs). You also need a machine with 16 GB RAM and ~150 GB free disk, and a conda/mamba environment. Everything below runs on a laptop except the STAR path, which needs ~32 GB. The part most tutorials skip, and the part that matters most for personal data, is step 7 onward. Whole-blood expression is dominated by cell composition, and if you do not model that, you will mostly be measuring your neutrophil fraction.

Set up once:

mamba create -n rnaseq -c conda-forge -c bioconda \
  salmon=1.10.3 fastp falco multiqc seqkit samtools star subread
mamba activate rnaseq

1. Verify the FASTQs before you compute anything

Check the checksums your sequencing provider gave you, then look at the actual read structure. Two minutes here saves a day.

md5sum -c checksums.md5
seqkit stats -a -j 4 *.fastq.gz

seqkit stats -a gives read count, min/avg/max length, Q20/Q30 percentages, and N50. Things to confirm: R1 and R2 have identical read counts, read length matches what you paid for (2×100 or 2×150), and Q30 is above ~85%. If R1 and R2 counts differ, the files are truncated or were re-paired badly, and every downstream step will silently produce garbage.

Run per-file QC with falco (a drop-in FastQC replacement that is roughly an order of magnitude faster) and aggregate:

falco -o qc/ *.fastq.gz
multiqc qc/ -o qc/multiqc

The three plots worth reading: per-base sequence quality (a drop below Q20 in the last 10–20 bases of R2 is normal and salmon handles it), overrepresented sequences (adapter read-through, or a single gene like HBB eating your library), and duplication level (high duplication at low complexity means the library was over-amplified from too little input RNA).

Trim adapters and polyG with fastp. Do not aggressively quality-trim. Modern pseudo-aligners and splice-aware aligners handle low-quality tails better than trimming does, and hard trimming biases coverage.

fastp -i S1_R1.fastq.gz -I S1_R2.fastq.gz \
      -o trim/S1_R1.fq.gz -O trim/S1_R2.fq.gz \
      --detect_adapter_for_pe --trim_poly_g \
      --length_required 36 --thread 8 \
      --html qc/S1.fastp.html --json qc/S1.fastp.json

--trim_poly_g matters on NovaSeq/NextSeq two-color chemistry, where the absence of signal is read as G. Skipping it inflates your unmapped fraction.

2. Choose quantification: salmon, not STAR, unless you need the BAM

We use salmon in selective-alignment mode with a decoy-aware index. It gives gene-level abundances that agree closely with STAR + featureCounts, runs in about 5–10 minutes per sample on 8 threads, needs under 16 GB of RAM, and produces bootstrap replicates you can use to quantify technical uncertainty per gene.

Use STAR instead when you need any of: splice junction counts, variant calling from RNA (allele-specific expression), fusion detection, or coverage tracks to eyeball in IGV. STAR’s human index build needs ~32 GB RAM and the alignments are 3–8 GB per sample. For a longitudinal personal profile where the question is “how does gene X move over time,” the BAM is not worth the storage.

Do not use a transcriptome-only index without decoys. Reads from intronic and intergenic sequence will be forced onto transcripts they do not come from, and the inflation is not uniform across genes.

3. Build the decoy-aware index

Use GENCODE. Its annotation and the FASTA are versioned together, which makes the tx2gene mapping trivial.

V=46
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${V}/gencode.v${V}.transcripts.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${V}/GRCh38.primary_assembly.genome.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${V}/gencode.v${V}.primary_assembly.annotation.gtf.gz

# decoys = all genome sequence names
grep "^>" <(gunzip -c GRCh38.primary_assembly.genome.fa.gz) \
  | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v${V}.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_gencode_v${V} --gencode -k 31

--gencode strips the pipe-delimited GENCODE FASTA headers down to the transcript ID. Forget it and your tx2gene join fails. -k 31 is the default and correct for reads ≥75 bp. Drop to -k 23 only for reads under 50 bp. Indexing takes 30–60 minutes and about 20 GB of RAM.

4. Quantify

for s in S1 S2 S3; do
  salmon quant -i salmon_idx_gencode_v46 -l A \
    -1 trim/${s}_R1.fq.gz -2 trim/${s}_R2.fq.gz \
    -p 8 --validateMappings \
    --seqBias --gcBias --posBias \
    --numBootstraps 30 \
    -o quant/${s}
done

-l A lets salmon infer strandedness. Read quant/S1/lib_format_counts.json afterward and confirm it converged on one library type across all samples (ISR for most stranded dUTP kits, IU for unstranded). If one sample infers a different type, that sample was prepped differently or is contaminated.

--gcBias and --seqBias correct fragment-level biases that differ between library preps. On a longitudinal personal dataset where samples were prepped weeks or months apart, these corrections are doing real work and not cosmetic work. The --posBias flag helps when RNA integrity varies between draws, which it will.

--numBootstraps 30 writes inferential replicates. They cost little and let you tell “this gene moved” from “this gene has three reads and the estimate is noise.”

Record the mapping rate from quant/S1/logs/salmon_quant.log. For clean poly-A selected whole blood, expect roughly 75–90%. Below 60%, stop and diagnose (step “Common problems”).

5. Import to R and build the gene matrix

library(tximport); library(GenomicFeatures); library(DESeq2)

txdb <- makeTxDbFromGFF("gencode.v46.primary_assembly.annotation.gtf.gz")
k    <- keys(txdb, keytype = "TXNAME")
tx2gene <- AnnotationDbi::select(txdb, k, "GENEID", "TXNAME")

files <- file.path("quant", samples$id, "quant.sf")
names(files) <- samples$id
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ 1)

tximport with type="salmon" (default countsFromAbundance = "no") passes DESeq2 the estimated counts plus an average transcript-length offset per gene per sample. That offset is the correct way to handle the fact that isoform usage shifts change effective gene length. If you need a matrix for something outside DESeq2, use countsFromAbundance = "lengthScaledTPM" and treat it as counts.

Filter before anything else:

keep <- rowSums(counts(dds) >= 10) >= ceiling(ncol(dds) / 3)
dds  <- dds[keep, ]

In whole blood this typically leaves 13,000–16,000 genes out of ~62,000 GENCODE entries. Most of what you dropped is lncRNA and pseudogene with zero support.

6. QC the matrix as well as the reads

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

Four numbers per sample decide whether a library is usable:

• Library size after filtering. Anything under ~10M assigned reads has poor power for low-expressed genes.
• Globin fraction: HBB, HBA1, HBA2, and HBD. In PAXgene whole blood without globin depletion this routinely runs 30–70% of reads. Every percent of globin is a percent you did not spend on the transcriptome. If you can choose, use globin-depleted library prep.
• Mitochondrial fraction (^MT-). Elevated MT with degraded RIN usually means a slow or warm sample handoff.
• Correlation to the sample-set median (Spearman on vst values). An outlier below ~0.95 when the rest sit above 0.98 is a prep failure and not biology.

On PCA, if PC1 separates your samples by draw date rather than by anything you care about, you are looking at batch and not signal. Write that down before interpreting anything.

7. Estimate cell composition, because in blood it explains most of the variance

Whole-blood bulk RNA-seq is a weighted average of neutrophils, lymphocytes, monocytes, platelets, and residual reticulocytes. A modest shift in neutrophil fraction (which moves with sleep, exercise, a mild infection, and time of day) produces hundreds of “differentially expressed” genes that have nothing to do with regulation inside any cell type.

Run deconvolution with a validated signature matrix. Options: CIBERSORTx with LM22 (22 leukocyte subsets, well characterized, web or Docker), or granulator in R which wraps several deconvolution methods over a common interface so you can check that they agree.

library(granulator)
tpm <- txi$abundance[keep, ]
dec <- deconvolute(m = tpm, sigMatrix = sigMatrix_ABIS_S0)
plot_proportions(dec)

Two things to do with the output. First, look at it directly. A 2× change in monocyte fraction between two timepoints is itself a measurement worth having. Second, feed the estimated proportions back into the model as covariates. Include the two or three proportions that vary most, and always drop one to avoid perfect collinearity (proportions sum to one).

If you have paired CBC differentials from the same draws, use the measured percentages instead of the deconvolved ones. They are direct.

8. Model the contrast you have, and say so when there isn’t one

Case A, you have real groups with replication (for example, six draws in a fasted state and six fed):

design(dds) <- ~ neutrophil_frac + lymphocyte_frac + state
dds <- DESeq(dds)
res <- lfcShrink(dds, coef = "state_fed_vs_fasted", type = "apeglm")
summary(res, alpha = 0.05)

Use apeglm shrinkage before ranking by fold change. Unshrunk log2 fold changes for low-count genes are unbounded and will fill the top of your list with noise.

Case B, repeated draws on you over time with a continuous covariate (weeks on a protocol, a lab value, or a season):

design(dds) <- ~ neutrophil_frac + weeks
dds <- DESeq(dds, test = "LRT", reduced = ~ neutrophil_frac)

The LRT asks which genes change monotonically with weeks beyond composition. For non-monotone shapes, add a spline basis with splines::ns(weeks, df = 3) in the full design and keep the reduced design without it.

Case C, two samples, one condition each. There is no test here. DESeq2 will produce a table and it will be meaningless because there is no estimate of within-condition variance. What you can do: rank genes by log2 fold change among genes with adequate counts, use the salmon bootstrap replicates to see which of those fold changes exceed technical noise, and treat the result as a hypothesis to test with more draws. Say “this is unreplicated” out loud when you look at the list.

The general rule for personal data: a single timepoint is a point estimate with unknown intra-individual variance. Three or more draws under the same conditions give you that variance, and the variance is what turns a number into a signal.

9. Move to pathways, where the signal is more stable

Individual gene measurements in blood are noisy. Coordinated sets are much less so.

library(fgsea); library(msigdbr)

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

h <- msigdbr(species = "Homo sapiens", collection = "H")
pathways <- split(h$gene_symbol, h$gs_name)   # map Ensembl IDs to symbols first

fg <- fgsea(pathways, stats, minSize = 15, maxSize = 500, nPermSimple = 10000)
fg[order(padj)][1:20]

Rank on the Wald or LRT statistic, not on fold change alone, because the statistic already accounts for dispersion. Start with the Hallmark collection (50 sets, low redundancy). Go to Reactome or GO only after Hallmark points somewhere.

For per-sample scores across your timeline, use GSVA::gsva() with ssgseaParam on the vst matrix. That gives you one number per pathway per draw, which plots cleanly over time and is far more reproducible run-to-run than any single gene.

10. Interpretation limits

Blood transcriptome measurements are sensitive to time of day, recent exercise, recent meals, sleep, and any low-grade infection you did not notice. Standardize draw conditions (same time of day, same fasting state, same posture, same tube type) or you are measuring those variables instead of what you intended. An interferon signature or an inflammatory pathway score is a measurement. Expression is not diagnosis, and a score is not a finding about your health. Anything that looks clinically relevant belongs in front of a physician who can order the appropriate confirmatory tests and interpret them in context.

Common problems

Mapping rate under 60%. Run a quick contamination check: salmon quant against a bacterial or mouse index, or blast a few hundred unmapped reads. Common causes are adapter read-through you did not trim, rRNA carryover (check RNA45SN* and RNA18SN* counts), a species mix-up, or an index built without --gencode so nothing joins.

All samples cluster by date on PC1. Batch. If batches are confounded with your condition, the experiment cannot answer the question and no correction fixes it. If they are crossed, add batch to the design. Use limma::removeBatchEffect only for visualization, never on the matrix you feed to DESeq2.

Thousands of significant genes and they are all immune-related. Almost always cell composition. Recheck step 7 and add the proportions as covariates.

tximport throws “transcripts missing from tx2gene”. Your transcript FASTA and GTF are different GENCODE releases, or you forgot --gencode at index time. Rebuild with matched versions.

Library type inferred inconsistently across samples. Check lib_format_counts.json for each. Mixed types means mixed prep protocols. Force a single -l ISR only after you have confirmed the kit chemistry, and note the change.

Enormous dispersion estimates for a handful of genes. Look at the raw counts. Usually one sample carries 90% of the counts for that gene, which is a mapping artifact rather than a biological event. Raise the count filter and re-run.

Bootstrap variance is large for a gene you care about. That gene has close paralogs or heavy isoform ambiguity (HLA genes, immunoglobulins, and gene families are the usual suspects). Report it at the family level, or align with STAR and inspect coverage directly.