How to Analyze Your Own RNA-Seq Data

By the end of this guide you will have taken raw RNA sequencing reads from your own blood and produced three things: a gene-level count matrix, a transcript-level abundance table with effective lengths, and a set of normalized expression values you can join against your genome and your blood biomarkers. You need the FASTQ files, which are paired-end and typically 2 × 150 bp. You also need a machine with at least 32 GB of RAM and 200 GB of free disk, Docker or Conda, and a working knowledge of the shell. Nothing here requires a cluster. A single sample of 40 million read pairs quantifies in about ten minutes on eight cores with a lightweight aligner. The harder part is interpretation, because a transcriptome is a snapshot of a regulatory state rather than a fixed text, and a single snapshot has no internal comparison group. We will spend the second half of the guide on what you can and cannot conclude from one, or a few, timepoints.

Before the steps, a word on what the transcriptome is relative to the genome and the proteome, since this is the question most people arrive with. Your genome is effectively invariant across your cells and your life: the same roughly 3.1 billion base pairs in a neutrophil and a hepatocyte. The transcriptome is the set of RNA molecules present in a particular sample of cells at a particular moment, and it varies by tissue, time of day, infection status, fasting state, and how long the blood tube sat on the bench. The proteome is a further step removed, and correlates with transcript abundance only moderately. Translation rate, protein half-life, and post-translational modification all intervene. That is why transcript and protein measurements are treated as complementary layers rather than substitutes for one another¹. The general framing of these layers as distinct measurement technologies with distinct inferential reach is standard in the omics literature².

1. Inspect what you were given

Start by confirming the library type, because every downstream parameter depends on it. Run fastqc and then multiqc over the raw reads:

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

The numbers to read off: total read pairs (you want 30–60 million for reliable gene-level quantification of blood), per-base quality (Q30 or better across most of the read), duplication rate, and adapter content. For whole-blood RNA, also check the GC distribution for a bimodal shape, which often indicates residual globin transcripts from red blood cells. Globin mRNA can consume 50–70% of reads in unenriched whole blood, which effectively throws away most of your sequencing depth for every other gene. If your provider used a globin depletion kit or PAXgene tubes with GLOBINclear, you should see this reflected in the fraction of reads mapping to HBB, HBA1, and HBA2 after step 3.

Then determine strandedness empirically rather than trusting the documentation. The simplest route is to run the quantifier in automatic detection mode (step 3) and read the inferred library type from its log. A mismatch here typically halves your counts and silently distorts antisense-overlapping genes.

2. Trim only what needs trimming

Aggressive trimming does more harm than good for quantification. Modern pseudoalignment and soft-clipping aligners tolerate adapter read-through reasonably well, and trimming shortens reads in ways that bias effective length estimates. We use fastp with light settings, mostly to remove adapters and poly-G artifacts from patterned flow cells:

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 \
  --length_required 36 \
  --qualified_quality_phred 15 \
  --thread 8 \
  --json fastp.json --html fastp.html

Do not enable quality trimming beyond this. Do not deduplicate RNA-seq reads unless your library used unique molecular identifiers, because in RNA-seq high-abundance transcripts legitimately produce identical fragments and removing them compresses your dynamic range at the top end.

3. Quantify with Salmon against a decoy-aware index

We recommend transcript-level quantification with salmon rather than a genome alignment followed by counting, for two reasons. It is roughly an order of magnitude faster, and it models fragment-level GC and positional bias while resolving reads that map to multiple isoforms of the same gene probabilistically instead of discarding them. The cost is that you get no BAM file, so you cannot inspect coverage over a specific exon or call variants from RNA. If you need that, run STAR in parallel later (step 7).

Build a decoy-aware index once, using GENCODE human transcripts plus the primary genome assembly as decoys, so that reads originating from intronic and intergenic sequence are not forced onto transcripts:

# GENCODE v46, human
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_46/gencode.v46.transcripts.fa.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_46/GRCh38.primary_assembly.genome.fa.gz

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

salmon index \
  -t gentrome.fa.gz \
  -d decoys.txt \
  -i salmon_index_gencode46 \
  -k 31 \
  --gencode \
  -p 8

The -k 31 default is right for 150 bp reads. Drop to -k 23 only if your reads are 50 bp or shorter. The --gencode flag strips the pipe-delimited GENCODE headers down to bare transcript IDs, which saves considerable grief later.

Then quantify:

salmon quant \
  -i salmon_index_gencode46 \
  -l A \
  -1 trim_R1.fastq.gz -2 trim_R2.fastq.gz \
  --validateMappings \
  --seqBias --gcBias \
  --numGibbsSamples 20 \
  -p 8 \
  -o quant_sample01

-l A auto-detects library type: check quant_sample01/lib_format_counts.json and confirm the inferred type (ISR for most stranded TruSeq-style kits, IU for unstranded). The --gcBias flag matters more than people expect for blood libraries, where fragment GC content varies with the depletion chemistry. The --numGibbsSamples 20 flag gives you posterior samples of transcript abundance, which is the only sound way to express uncertainty on isoform-level calls. Drop it if disk space is tight. Check the mapping rate in quant_sample01/logs/salmon_quant.log. Above 80% is normal for polyA-selected blood RNA. Below 60% means contamination, wrong species, heavy rRNA carryover, or degraded input.

4. Summarize to gene level and normalize

Transcript-level TPM is what Salmon reports, but gene-level counts are more stable and are what most downstream tools expect. Use tximport in R, which carries the effective transcript lengths through so that length-aware normalization stays correct:

library(tximport)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(GenomicFeatures)

# transcript -> gene map from the same GENCODE release
txdb <- makeTxDbFromGFF("gencode.v46.annotation.gtf.gz")
tx2gene <- AnnotationDbi::select(txdb, keys = keys(txdb, "TXNAME"),
                                keytype = "TXNAME", columns = "GENEID")

files <- c(s01 = "quant_sample01/quant.sf")
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM",
                ignoreTxVersion = TRUE)

Two normalization choices deserve explanation. TPM (transcripts per million) divides by transcript length and then by the library total, which makes values comparable within a sample across genes. It is the right unit for asking “which genes are high in me right now.” It is the wrong unit for comparing the same gene across samples, because the per-million denominator is sensitive to a handful of very high-abundance transcripts. For across-sample comparison use lengthScaledTPM counts fed into DESeq2’s median-of-ratios size factors, or edgeR’s TMM. With a personal longitudinal series this matters concretely: if one draw has residual globin and another does not, raw TPM for every other gene shifts, and the size-factor approach absorbs most of that.

Filter before any comparison. Keep genes with at least 10 counts in at least a quarter of your samples. With a single sample, keep genes above roughly 1 TPM and treat everything below as present/absent at best.

5. Establish a personal baseline before interpreting anything

This is the step most guides skip and the one that determines whether your data is useful. A single transcriptome tells you almost nothing about yourself in isolation, because you have no reference for what is normal for you. Population references help with orientation: GTEx whole blood medians give you a rough expectation, and the Human Protein Atlas consensus tissue data tells you whether a gene is blood-restricted or ubiquitous. But the between-person variance in blood transcript abundance is large, so a gene at the 10th percentile of GTEx may be perfectly stable for you.

What we recommend is three to five draws over several months, collected under matched conditions: same time of day (transcript abundance in blood follows a diurnal rhythm with substantial amplitude for immune and clock genes), same fasting state, same tube type, same handling time between draw and stabilization. Then compute, per gene, your own mean and coefficient of variation across draws. Genes with low personal CV become your informative markers, because a subsequent deviation in one of them is signal. Genes with high personal CV are noise generators and should be discounted regardless of how interesting their annotation is. This framing of a personal baseline as the reference against which later measurements are judged is what makes longitudinal designs so much more informative than single snapshots, and is the operating logic of most serious multi-omics profiling work³.

library(DESeq2)
dds <- DESeqDataSetFromTximport(txi_multi, colData = meta, design = ~ 1)
dds <- estimateSizeFactors(dds)
vsd <- vst(dds, blind = TRUE)

m <- assay(vsd)
personal_cv <- apply(m, 1, function(x) sd(x) / mean(x))
stable_genes <- names(sort(personal_cv))[1:2000]

Use the variance-stabilizing transform (vst) rather than log2(TPM + 1) for any clustering, correlation, or principal component work. Log-plus-pseudocount inflates the apparent variance of low-count genes, which then dominate your PCs.

6. Move from genes to pathways

Individual gene changes in a personal series are hard to interpret and easy to over-read. Aggregate signal at the pathway or cell-type level is more stable and more meaningful. Two approaches are worth running.

First, gene set scoring per sample. Use GSVA or single-sample GSEA over Hallmark and Reactome gene sets to produce a pathway-by-sample matrix, then look at which pathway scores move across your timepoints:

library(GSVA)
library(msigdbr)

hallmark <- msigdbr(species = "Homo sapiens", collection = "H")
gs <- split(hallmark$gene_symbol, hallmark$gs_name)
par <- gsvaParam(exprData = m, geneSets = gs, kcdf = "Gaussian")
scores <- gsva(par)

Second, cell-type deconvolution. Most variation in a whole-blood transcriptome is compositional: a shift in the neutrophil-to-lymphocyte ratio will move thousands of genes without any per-cell regulatory change. Run CIBERSORTx or immunedeconv against the LM22 signature matrix to estimate cell fractions, then check them against the differential white cell counts from your CBC drawn at the same visit. Agreement between the two is a good validity check on the whole pipeline. Disagreement usually means a handling problem, and neutrophil-specific transcripts are the most sensitive to delayed processing.

The interpretive point generalizes: transcriptomic signals become meaningful mainly when anchored to something measured independently, whether a protein level, a metabolite, or a clinical variable. Integrating layers is the standard route to reducing the false-discovery burden that any single high-dimensional layer carries on its own⁴. Concordance across layers, and the biological logic of why two layers should move together, is the substance of multi-omics interpretation rather than a formality⁵.

7. Add a genome alignment when you need allele-level or structural detail

Salmon cannot tell you which allele a transcript came from, whether an exon was skipped, or whether a variant in your genome disrupts splicing. For those questions, align to the genome:

STAR --runThreadN 8 \
  --genomeDir star_index_GRCh38_g46 \
  --readFilesIn trim_R1.fastq.gz trim_R2.fastq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts \
  --outSAMattributes NH HI AS nM NM MD \
  --outFilterMultimapNmax 20 \
  --twopassMode Basic \
  --outFileNamePrefix star_s01_

Two-pass mode improves detection of novel junctions, which is what you want if you are following up a predicted splice-altering variant from your whole-genome data. With the BAM you can run rMATS or SUPPA2 for splicing quantification, or use samtools mpileup at a specific heterozygous site from your VCF to measure allelic imbalance in expression. Allelic imbalance is one of the few things a single RNA-seq sample can establish well, because the two alleles are internally controlled within the same library.

Be careful about statistical power. A meaningful splicing or allelic call needs adequate read depth over the specific junction or site, which for a moderately expressed gene at 40 million read pairs may be only a few dozen reads. Compute the binomial confidence interval before you believe an imbalance.

Common problems

Pipeline choice changes results, and you should know how much. Different combinations of aligner, quantifier, and differential testing method produce overlapping but non-identical gene lists from the same reads. The magnitude of that disagreement is well documented. Treat any single pipeline’s output as one estimate rather than ground truth, and re-run anything you intend to act on through a second quantifier. Reproducible, version-pinned workflows are the practical answer here, which is why regulatory-grade transcriptomics frameworks specify the full processing chain rather than just the tool names⁶. The easiest way to get this for free is to run nf-core/rnaseq with -r 3.14.0 --aligner star_salmon and archive the pipeline_info directory alongside your results.

Annotation version mismatches are the most common silent failure. If your index is GENCODE v46 and your tx2gene map is v44, tximport will drop transcripts without warning and your gene-level counts will be quietly wrong for the affected loci. Pin one release and use it for the index, the GTF, and every downstream annotation lookup.

Sample handling dominates biology more often than people expect. Blood left at room temperature for hours before stabilization shows large, reproducible shifts in neutrophil and stress-response transcripts. If your timepoints were collected at different clinics or with different transit times, that difference will appear as your largest principal component, and no normalization method will remove it cleanly. Record handling metadata and include it in colData.

RNA integrity limits what is interpretable. Ask your provider for the RIN or DV200 value. Below RIN 7, 3’ bias in coverage distorts effective lengths and long transcripts are systematically underestimated. You can detect this from the transcript integrity numbers in RSeQC or from a coverage-versus-position plot in the MultiQC report.

Finally, resist the urge to read clinical meaning into single-gene changes. Transcript abundance in blood is a composite of cell composition, circadian phase, recent activity, and acute immune state, and the mapping from any one gene to a health outcome is generally weak. The measurement is valuable for tracking your own trajectory over time and for generating specific hypotheses to test with an independent assay. If something in your data appears to bear on a medical question, including a candidate variant from your genome or an unexpected immune signature, bring the data to a physician or a clinical geneticist. Interpretation of that kind requires a clinician who can see your history and order confirmatory testing.

Footnotes

1. Dhirendra Kumar, Gourja Bansal, Ankita Narang, et al. Integrating transcriptome and proteome profiling: Strategies and applications. PROTEOMICS, 2016. https://doi.org/10.1002/pmic.201600140 ↩

2. 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 ↩

3. Mohamad Hesam Shahrajabian, Wenli Sun. Survey on Multi-omics, and Multi-omics Data Analysis, Integration and Application. Current Pharmaceutical Analysis, 2023. https://doi.org/10.2174/1573412919666230406100948 ↩

4. Indhupriya Subramanian, Srikant Verma, Shiva Kumar, et al. Multi-omics Data Integration, Interpretation, and Its Application. Bioinformatics and Biology Insights, 2020. https://doi.org/10.1177/1177932219899051 ↩

5. Biswapriya B Misra, Carl Langefeld, Michael Olivier, et al. Integrated omics: tools, advances and future approaches. Journal of Molecular Endocrinology, 2019. https://doi.org/10.1530/jme-18-0055 ↩

6. M.C. Verheijen, T.W. Gant, W. Tong, et al. R-ODAF: Omics data analysis framework for regulatory application. Toxicology Letters, 2021. https://doi.org/10.1016/s0378-4274(21)00539-7 ↩