How to Analyze Your Own RNA-Seq Data

By the end of this guide you will have a gene-by-sample expression matrix from your own RNA, and a set of quality metrics that tell you whether to trust it. You will also have a comparison of your expression against a population reference, and a framework for asking whether a gene changed between two of your own timepoints. You need the raw FASTQ files from your sequencing run (paired-end, ideally 2×100 or 2×150 bp, 30 to 60 million read pairs per sample for bulk tissue), roughly 300 GB of free disk, 32 GB of RAM, and 8 or more cores. Software: fastp, salmon, STAR (only if you want splicing or allele-specific work), samtools, and R 4.3+ with tximport, DESeq2, limma, and biomaRt. Everything below assumes Linux or macOS with conda or micromamba available. Install the aligners through conda install -c bioconda salmon star fastp samtools rather than building from source.

1. Find out what library you have

Before you run anything, read the library prep description on your sequencing report, because three details change every downstream decision. First, selection chemistry. Poly(A) selection captures polyadenylated mRNA and discards most non-coding RNA. Ribosomal depletion (Ribo-Zero, RiboCop) keeps introns, lncRNAs, and unprocessed transcripts. That inflates intronic read fraction and changes what a “low” gene count means. Second, strandedness. Most modern kits are reverse-stranded (dUTP chemistry). That is ISR in Salmon notation and --library-type fr-firststrand elsewhere. Third, for blood: whether globin depletion was performed. Whole blood RNA without globin depletion can devote well over half of all reads to HBB, HBA1, HBA2, and HBD, which means a nominal 40 million read pairs might leave you 15 million informative ones.

If the report is silent, let the tools tell you. Run Salmon with -l A on one sample and read the inferred library type out of lib_format_counts.json. It is reliable when you have more than a few million mapped reads.

2. Build a decoy-aware transcriptome index

We recommend selective alignment against a transcriptome with genomic decoys rather than a full spliced genome alignment, for the simple reason that it is five to ten times faster, uses under 16 GB of RAM, and gives transcript-level estimates that are as good or better for expression quantification. The decoy sequences matter: without them, reads originating from unannotated genomic regions get force-assigned to whatever transcript they resemble most, which produces phantom expression in low-expressed genes.

# GENCODE human release 44, primary assembly
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.transcripts.fa.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/GRCh38.primary_assembly.genome.fa.gz

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

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

Use -k 31 for reads of 75 bp or longer and drop to -k 23 only if your reads are 50 bp. The --gencode flag strips the pipe-delimited GENCODE headers down to the ENST identifier, which saves you a string-parsing step later. Index construction takes about 25 minutes and 20 GB of RAM.

3. Trim and QC the reads

fastp does adapter detection, quality trimming, poly-G trimming (needed for NovaSeq and NextSeq two-color chemistry, where a dark base is read as G), and a QC report in one pass.

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 --poly_g_min_len 10 \
  --qualified_quality_phred 15 --unqualified_percent_limit 40 \
  --length_required 36 \
  --thread 8 \
  --json sample.fastp.json --html sample.fastp.html

Resist aggressive quality trimming. Salmon’s selective alignment handles moderately low-quality tails, and hard trimming at Q30 biases coverage toward read starts. The numbers to check in the JSON: duplication.rate (above 0.6 in a poly(A) library from 30M reads suggests low library complexity from a small RNA input), adapter_cutting.adapter_trimmed_reads fraction (above 30% means short inserts and you should expect fragment length estimates near 150 bp), and the per-base N content.

4. Quantify

salmon quant -i salmon_idx_v44 -l A \
  -1 trim_R1.fastq.gz -2 trim_R2.fastq.gz \
  --seqBias --gcBias --posBias \
  --numBootstraps 30 \
  --validateMappings \
  -p 16 -o quant/sample01

--seqBias and --gcBias fit models for random-hexamer priming bias and fragment GC content, which are the two largest systematic distortions in Illumina RNA-seq and which matter more when you are comparing your sample to an external reference processed by someone else. --numBootstraps 30 gives you inferential replicates, the uncertainty in assigning reads among transcripts that share sequence. You will want them when a “changed” gene turns out to be one paralog of a family that Salmon cannot distinguish.

Check logs/salmon_quant.log for the mapping rate. For human poly(A) libraries against GENCODE, expect 75 to 90 percent. Below 60 percent, stop and investigate before interpreting anything: the usual causes are the wrong species index, heavy rRNA or globin carryover, or genomic DNA contamination.

5. Load into R and build the expression matrix

Salmon writes transcript-level estimates. For most personal analysis you want gene-level counts, and tximport with countsFromAbundance = "lengthScaledTPM" gives you a count matrix that already accounts for changes in average transcript length between samples, so you can hand it straight to DESeq2 or limma.

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

files <- file.path("quant", samples$id, "quant.sf")
names(files) <- samples$id

# tx2gene from the GENCODE GTF, two columns: ENST, ENSG
tx2gene <- readr::read_tsv("tx2gene.v44.tsv")

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM",
                ignoreTxVersion = FALSE)

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

# keep genes with >= 10 counts in at least as many samples as the smallest group
keep <- rowSums(counts(dds) >= 10) >= 2
dds  <- dds[keep, ]

vsd <- vst(dds, blind = TRUE)   # variance-stabilized, log2-like scale

The vst output is what you use for clustering, PCA, and any distance calculation. Never run PCA on raw TPM or on log2(TPM+1) if you can avoid it, because the variance of a log-transformed count is strongly dependent on its mean and the low-expressed genes will dominate your distance metric.

6. Run the sanity checks that catch real mistakes

These four checks take ten minutes and catch most of the errors that would otherwise send you chasing a biological story that does not exist.

Sex concordance: compute the VST values for XIST (ENSG00000229807) and RPS4Y1, DDX3Y, UTY. A sample from a person with a Y chromosome has XIST near the floor and the Y genes clearly expressed, and the reverse otherwise. Anything intermediate means sample mix-up, contamination, or index hopping on a shared flow cell.

Globin fraction: sum the TPM of HBB, HBA1 and HBA2, plus HBD, HBG1 and HBG2. In globin-depleted whole blood this should sit in the low single-digit percent. If it is 40 percent or more, your effective sequencing depth is a fraction of what you paid for and low-expressed genes will be noisy.

Ribosomal and mitochondrial fraction: sum TPM over the MT- genes and over rRNA-biotype genes from the GTF. A high mitochondrial fraction in a tissue sample often indicates degraded or stressed material.

Degradation: if you have a Bioanalyzer or TapeStation RIN from the lab, keep it in your sample sheet. RIN below about 7 produces a 3′ coverage bias that systematically deflates long transcripts. If you do not have RIN, approximate it by plotting mean coverage as a function of relative position along transcripts using qualimap rnaseq on a genome alignment.

7. Compare a single sample against a population reference

Personal RNA-seq usually starts with one sample and no internal control group, which is the hardest situation in transcriptomics. The workable approach is to put your sample on the same scale as a public reference cohort from the matching tissue. GTEx v8 whole blood has 755 samples with gene-level TPM available, and the Human Protein Atlas consensus tissue data gives a coarser but broader alternative.

The comparison must be done carefully. Quantify your sample against the same annotation version the reference used, or subset both to the intersection of gene IDs. Then quantile-normalize your sample’s TPM vector to the reference distribution, take log2(x + 1), and compute a per-gene z-score against the reference mean and standard deviation. Genes beyond |z| of 3 are candidates. They are not findings. Expect dozens of them by chance in 15,000 tested genes, plus a systematic inflation from the batch difference between your library prep and GTEx’s.

A more principled version of this is OUTRIDER, which fits an autoencoder to a reference count matrix, controls for latent confounders, and returns per-gene p-values for expression outliers in a new sample. It was built for rare disease diagnostics and it is the right tool if your question is “is any gene in me expressed at an implausible level”. Extreme results here are hypotheses to bring to a clinical geneticist, who can order a validated assay in a CLIA-certified laboratory. Nothing computed this way is a diagnosis.

8. Deconvolve cell type composition before you interpret blood expression

Most of the variance in whole blood RNA-seq between people, and within one person over time, is composition. Composition is the ratio of neutrophils to lymphocytes to monocytes. It shifts with acute infection, exercise, cortisol rhythm, and time of day. If you do not estimate composition you will attribute a composition shift to a regulatory change.

Run CIBERSORTx (web, LM22 signature matrix) or the R package granulator, or xCell for enrichment scores rather than proportions. Feed them TPM, not counts, and not VST values. Then include the leading one or two estimated proportions as covariates in any downstream model. In practice, adding neutrophil and lymphocyte fractions to a design absorbs a large share of between-timepoint variance in blood.

9. Model change across your own timepoints

With repeated sampling of one person, the design differs from a case-control study: your comparison is within subject, which removes genetic background as a confounder and makes moderate effects detectable with few samples. With two timepoints and no replicates, DESeq2 cannot estimate dispersion from the data, so you either supply a dispersion value from a comparable public dataset or restrict yourself to ranking by fold change with bootstrap-derived uncertainty from Salmon.

With four or more timepoints, this becomes tractable:

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

Setting lfcThreshold inside results() tests against a real effect size rather than against zero, which is the right null when your sample size is small and nearly everything is slightly nonzero. apeglm shrinkage keeps low-count genes from topping your fold-change ranking.

Draw blood at a consistent time of day and a consistent fasting state across timepoints. Circadian and postprandial effects on blood transcript levels are large enough to swamp whatever you were trying to measure.

10. Align to the genome when you need splicing or allele-specific expression

Salmon cannot see novel junctions or count reads per allele. If you have whole-genome sequencing on the same person and want to ask whether a rare variant affects splicing or allelic balance, align with STAR in two-pass mode.

STAR --runMode genomeGenerate --genomeDir star_idx \
  --genomeFastaFiles GRCh38.primary_assembly.genome.fa \
  --sjdbGTFfile gencode.v44.annotation.gtf --sjdbOverhang 149 --runThreadN 16

STAR --genomeDir star_idx --readFilesIn trim_R1.fastq.gz trim_R2.fastq.gz \
  --readFilesCommand zcat --twopassMode Basic \
  --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts \
  --waspOutputMode vW --varVCFfile hets.vcf \
  --outSAMattributes NH HI AS nM vW \
  --runThreadN 16 --outFileNamePrefix star/sample01.

The WASP mode flags reads whose mapping changes when the alternate allele is substituted, which removes the reference-mapping bias that otherwise makes every heterozygous site look skewed toward the reference allele. Filter to vW == 1 before counting with gatk ASEReadCounter. For splicing outliers, FRASER models intron excision ratios against a reference cohort and is the counterpart to OUTRIDER for aberrant splicing. Both need a reference panel of at least 50 to 100 samples from the same tissue to behave, and both produce candidate events that require clinical interpretation rather than self-interpretation.

Common problems

Mapping rate below 60 percent. Check the FastQC overrepresented sequences. If they are rRNA, the depletion step failed. If they are adapter dimer, the library had low insert. If neither, confirm you built the index from the same genome build as your annotation, and that you did not mix GENCODE transcript IDs with an Ensembl tx2gene table stripped of version suffixes.

tximport drops most of your transcripts. This is nearly always a version-suffix mismatch between quant.sf (ENST00000456328.2) and your tx2gene (ENST00000456328). Set ignoreTxVersion = TRUE or fix the table.

Every PCA separates by batch instead of biology. Do not delete the batch effect from the matrix you feed into DESeq2. Put the batch term in the design formula. Use limma::removeBatchEffect on the VST matrix only for visualization, and say so in your notes.

A gene looks dramatically changed but sits at 12 counts. Look at the bootstrap standard deviation from Salmon before you believe it. Genes below roughly 20 counts in every sample carry fold-change estimates dominated by sampling noise even after shrinkage.

You want to confirm one specific gene. RT-qPCR remains the practical orthogonal check, and it answers a question RNA-seq handles poorly: precise quantification of one target across many samples at low cost. Reverse transcription converts RNA to cDNA first, so qPCR does measure RNA when run as RT-qPCR. With careful primer design it resolves sequence variants that short-read counting blurs together, including single-nucleotide differences between closely related small RNA isoforms.¹ Its accuracy depends on the reference genes you normalize to, and it has been shown to give quantitative profiles comparable to sequencing when the normalization is handled properly.² Use it to confirm a specific hypothesis. Do not use it to survey the transcriptome.

The transcriptome alone does not tell you about protein abundance. Transcript and protein levels correlate only moderately across genes, and studies that measure both in the same samples routinely find regulation that appears at the protein layer with no transcript signal, particularly for metabolic enzymes.³ If you have proteomics on the same person, analyze the layers jointly rather than sequentially, using a tool built to align features across omics types on a common sample axis.⁴

Footnotes

1. Sandra Franco, Raquel Pluvinet, Jose Francisco Sanchez-Herrero, et al. Rapid and accurate quantification of isomiRs by RT-qPCR. Scientific Reports, 2022. https://doi.org/10.1038/s41598-022-22298-7 ↩

2. Ching Jian, Panu Luukkonen, Hannele Yki-Järvinen, et al. Quantitative PCR provides a simple and accessible method for quantitative microbiota profiling. PLOS ONE, 2020. https://doi.org/10.1371/journal.pone.0227285 ↩

3. Mika Caplan, Karli J. Wittorf, Kasidy K. Weber, et al. Multi-omics reveals mitochondrial metabolism proteins susceptible for drug discovery in AML. Leukemia, 2022. https://doi.org/10.1038/s41375-022-01518-z ↩

4. David Márquez-Oller, Andrea Pauli, Jörg Fallmann. BRIDGE: an interactive application for multi-omics data analysis, visualization and integration. Bioinformatics, 2026. https://doi.org/10.1093/bioinformatics/btag558 ↩