How to Analyze Your Own RNA-Seq Data

By the end of this guide you will have a transcript-level and gene-level expression matrix from your own blood RNA-seq, a quality-control report you can defend, and a set of comparisons that are statistically meaningful given the number of samples you have. You need the raw reads as paired-end FASTQ files, typically 30 to 60 million read pairs per sample, gzipped at 3 to 8 GB per file. You also need a machine with at least 16 GB of RAM and 200 GB of free disk, a conda or container runtime, and a reference transcriptome. If you have only one timepoint, you will end up with an expression profile and rank-based interpretation. If you have three or more timepoints, you can do differential expression, which is where transcriptomics becomes genuinely informative about your own biology. This is measurement and interpretation. Nothing here is a diagnosis, and any result that seems clinically meaningful belongs in front of a physician who can order a validated test.

1. Understand what you are measuring before you touch the data

RNA sequencing quantifies how many copies of each RNA molecule were present in your sample at the moment it was stabilized. Transcriptomics is the broader field, the study of the full set of RNA transcripts in a cell or tissue, and RNA-seq is the dominant measurement technology used to do it. They are not synonyms: microarrays and targeted qPCR panels are also transcriptomics, and RNA-seq can be used for purposes beyond expression, including transcript discovery, splice variant detection, allele-specific expression, and fusion detection ¹.

Two properties of the measurement will shape everything downstream. First, RNA-seq is compositional. You measure proportions rather than absolute concentrations, so a large increase in one transcript mechanically depresses the apparent abundance of everything else. Normalization exists to correct for this and it does so imperfectly. Second, blood RNA is dominated by a handful of cell types, and a change in your neutrophil-to-lymphocyte ratio will produce hundreds of apparently differentially expressed genes that reflect cell composition rather than regulation within any cell. Bulk RNA-seq measures a weighted average across cells, which is the main reason single-cell and multimodal single-cell methods exist ².

People ask whether transcriptomics is better than proteomics. They answer different questions. mRNA levels explain a substantial but far from complete fraction of protein variation, transcripts turn over in minutes to hours while many plasma proteins persist for days, and secreted proteins in plasma come from tissues whose RNA you are not sampling. Read RNA-seq as a readout of current transcriptional state in the sampled tissue, and read plasma proteomics as a readout of circulating protein abundance. Neither substitutes for the other.

2. Set up the environment and get the reference

Use conda (or micromamba, which is faster) with the bioconda channel. Pin versions so your results are reproducible six months from now when you add another timepoint.

micromamba create -n rna -c conda-forge -c bioconda \
  fastp=0.23.4 salmon=1.10.3 fastqc=0.12.1 multiqc=1.21 \
  samtools=1.19 star=2.7.11b r-base=4.3 bioconductor-deseq2=1.42.0 \
  bioconductor-tximeta=1.20.1 bioconductor-fgsea=1.28.0
micromamba activate rna

For the reference, we use GENCODE human release 44 or later, primary assembly. Download the transcript FASTA and the matching GTF annotation, and keep them together with a note of the release number. Mixing annotation versions across timepoints is one of the most common sources of phantom differential expression.

wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.transcripts.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.primary_assembly.annotation.gtf.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/GRCh38.primary_assembly.genome.fa.gz

Build a decoy-aware Salmon index. The decoy sequences are the genome itself, which lets Salmon recognize reads that came from intronic or intergenic regions rather than forcing them onto a transcript. Without decoys you get inflated counts for genes that happen to share sequence with abundant pre-mRNA.

# concatenate transcripts + genome, mark genome sequences as decoys
zcat gencode.v44.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz | gzip > gentrome.fa.gz
zgrep '^>' GRCh38.primary_assembly.genome.fa.gz | cut -d ' ' -f1 | sed 's/>//g' > decoys.txt

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

The --gencode flag strips the pipe-delimited GENCODE headers down to the transcript ID, which saves you a string-parsing headache later. -k 31 is the default k-mer length and is right for reads of 75 bp or longer. Drop to -k 23 only if your reads are shorter than 50 bp. Index building takes roughly 25 to 40 minutes and peaks around 16 GB of RAM.

3. Quality-control and trim the reads

Run FastQC first so you know what you have, then trim. We use fastp because it does adapter detection, quality filtering, polyG trimming, and reporting in one pass at high speed.

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o sample_R1.trim.fastq.gz -O sample_R2.trim.fastq.gz \
  --detect_adapter_for_pe \
  --trim_poly_g \
  --qualified_quality_phred 15 --unqualified_percent_limit 40 \
  --length_required 36 \
  --thread 8 \
  --json sample.fastp.json --html sample.fastp.html

Resist the urge to trim aggressively on quality. Modern aligners and selective-alignment quantifiers handle moderate base-quality decay well, and hard trimming shortens reads in ways that bias transcript length estimation. The --trim_poly_g flag matters if your data came from a NovaSeq or other two-color chemistry instrument, where the absence of signal is read as G and produces runs of spurious Gs at read ends.

Four numbers from the report deserve your attention. Duplication rate above 50 percent in a standard poly-A library suggests low input RNA and over-amplification. Adapter content above a few percent means a short insert size and reduced effective coverage. GC distribution with a second peak often indicates bacterial or index-hopping contamination. And if your library is a globin-depleted whole-blood preparation, check the residual fraction of reads mapping to HBB and HBA1/HBA2 later in the quantification step: unsuccessful depletion can leave 30 to 50 percent of reads on three genes, which destroys your effective depth for everything else.

4. Quantify with Salmon

We recommend selective alignment with Salmon over a full genome alignment for personal expression profiling. It runs in roughly ten minutes per sample on eight cores instead of an hour or more, it uses 4 to 6 GB of RAM instead of 30 or more, and it estimates transcript-level abundances with a probabilistic model that handles multi-mapping reads correctly rather than discarding them. Alignment-based workflows remain the right choice when you need the actual coordinates. That is the case for novel splice junction discovery, fusion detection, or variant calling from RNA ³.

salmon quant \
  -i salmon_index_gencode44 -l A \
  -1 sample_R1.trim.fastq.gz -2 sample_R2.trim.fastq.gz \
  -p 8 --validateMappings \
  --seqBias --gcBias --posBias \
  --numGibbsSamples 20 \
  -o quant/sample

The -l A flag lets Salmon infer library type, and you should check the inferred value in lib_format_counts.json against what your provider says the protocol was. The --gcBias and --seqBias flags fit fragment-level bias models and cost a few minutes of runtime, which is worth it when you are comparing samples prepared in different batches. The --numGibbsSamples 20 flag produces posterior samples so you can see which transcript-level estimates are poorly identified, which happens routinely for genes with many similar isoforms.

The mapping rate in logs/salmon_quant.log is your first real quality gate. For a poly-A selected human blood library against a decoy-aware GENCODE index, expect 75 to 90 percent. Below 65 percent, stop and find out why before interpreting anything.

5. Import to R and build the gene-level matrix

Salmon writes transcript-level estimates in quant.sf: transcript name, length, effective length, TPM, and estimated read count. For most personal interpretation you want gene-level counts, because gene-level estimates are far more stable than isoform-level ones. Use tximeta. It reads the index metadata, identifies the annotation automatically, and attaches transcript-to-gene mapping without you hand-building it.

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

coldata <- data.frame(
  files = file.path("quant", samples, "quant.sf"),
  names = samples,
  timepoint = factor(c("t1","t1","t2","t2","t3","t3")),
  stringsAsFactors = FALSE
)
stopifnot(all(file.exists(coldata$files)))

se  <- tximeta(coldata)                 # transcript level
gse <- summarizeToGene(se)              # gene level, length-corrected
dds <- DESeqDataSet(gse, design = ~ timepoint)

# filter low-count genes: at least 10 counts in at least 3 samples
keep <- rowSums(counts(dds) >= 10) >= 3
dds  <- dds[keep, ]

That filter typically takes you from about 62,000 annotated features down to 13,000 to 16,000 expressed genes in whole blood. The discarded features are mostly pseudogenes, lncRNAs below detection, and transcripts of tissues you did not sample. Filtering before testing improves power because it reduces the multiple-testing burden, and it is standard practice in current recommendations ³.

For any visualization, clustering, or distance calculation, use a variance-stabilized transformation rather than raw TPM or log-CPM. Counts have mean-dependent variance, and untransformed log values give low-expression genes wildly inflated apparent variability.

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

If you have samples from multiple draws, look at PC1 and PC2 before doing anything else. In our experience with personal longitudinal blood data, the first principal component usually tracks white blood cell composition. It does not usually track any condition you are studying. Confirm this by correlating PC1 against marker genes for the following cell types.

• neutrophils (DEFA3, FCGR3B, CSF3R)
• T cells (CD3D, CD3E, IL7R)
• platelets (PPBP, PF4)

If PC1 is a composition axis, include an estimated cell-fraction covariate in your design or accept that your differential expression results describe a shift in cell mixture.

6. Test for changes, with realistic limits on what n permits

With two or more samples per condition you can run DESeq2. Its negative binomial model with empirical Bayes shrinkage of dispersion estimates is what makes small-n RNA-seq tractable at all.

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

Use lfcShrink with apeglm for ranking. Unshrunken log2 fold changes for low-count genes are dominated by sampling noise, and a gene going from 3 counts to 12 counts will otherwise outrank a genuine two-fold change in a well-measured gene. Set a fold-change threshold as well as an adjusted p-value threshold, because with adequate depth you can reach statistical significance on biologically meaningless 10 percent differences.

With a single timepoint you cannot do differential expression at all, and you should not try. What you can do is rank your transcriptome and ask which pathways sit at the extremes relative to a reference cohort, or compare your own timepoints once you have more than one. Pathway-level analysis is also more reproducible than gene-level analysis because it aggregates many weak signals, and it is where most of the interpretable biology of an RNA-seq experiment lives ⁴.

library(fgsea)
stats <- res$log2FoldChange / res$lfcSE       # signed test statistic
names(stats) <- rowData(dds)$gene_name[match(rownames(res), rownames(dds))]
stats <- stats[!is.na(names(stats)) & !duplicated(names(stats))]
fg <- fgsea(pathways = gmtPathways("c2.cp.reactome.v2023.2.Hs.symbols.gmt"),
            stats = stats, minSize = 15, maxSize = 500, nPermSimple = 10000)
head(fg[order(fg$padj), c("pathway","padj","NES","size")], 20)

Rank by the signed test statistic rather than fold change alone, because that ranking incorporates how well each gene was measured. Interpret normalized enrichment scores as directional hypotheses about coordinated transcriptional programs, and treat any single gene inside an enriched set as unconfirmed until you see it in an independent sample.

7. Extend beyond expression when the question warrants it

The same FASTQ files support analyses that transcript-level quantification cannot address. Non-coding RNA is a large and regulatorily important fraction of the transcriptome, and studies that pair mRNA with lncRNA and small-RNA profiling consistently find coordinated non-coding responses that expression-only analysis misses ⁵. Splicing is another axis: differential isoform usage can occur with no change in total gene expression, and detecting it requires either a splicing-aware tool on quantifications or a genome alignment with junction counts.

For anything requiring coordinates, align with STAR and keep the sorted BAM.

STAR --genomeDir star_index_gencode44 \
  --readFilesIn sample_R1.trim.fastq.gz sample_R2.trim.fastq.gz \
  --readFilesCommand zcat --runThreadN 8 \
  --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts \
  --outFilterMultimapNmax 20 --alignSJoverhangMin 8 \
  --outFileNamePrefix star/sample.

STAR needs about 32 GB of RAM for a human index, which is the main reason we do not use it as the default. The catalog of downstream methods is large and it is worth surveying what exists before writing custom code, since mature tools already cover fusion detection, allele-specific expression, RNA editing, and transcript assembly ⁶.

Common problems

Low mapping rate with high-quality reads almost always means a reference mismatch. Check that you built the index from the same GENCODE release as your annotation, that you are using a human index and not a leftover mouse one, and that the library was poly-A selected rather than total RNA. Total RNA libraries contain substantial intronic and rRNA-derived sequence, and against a transcriptome index they map at 40 to 60 percent even when the data are fine.

A dominant first principal component that separates samples by draw date is batch effect, cell composition, or both. Include the batch as a covariate in the DESeq2 design if the batches are not confounded with your comparison of interest. If batch and condition are perfectly confounded, the effect is not separable and no software will rescue it. Plan collections so that each timepoint appears in more than one processing batch.

Thousands of significant genes from a comparison you expected to be subtle is a warning. It is not a result. Inspect the top hits for immune and platelet markers, check the estimated size factors for one sample that is wildly different from the rest, and confirm your RNA integrity numbers. Degraded RNA produces a systematic 3’ bias that mimics differential expression in a length-dependent way.

Isoform-level results that shift when you rerun with different parameters are telling you the estimates are not identified. Look at the Gibbs posterior spread from --numGibbsSamples: transcripts with wide posteriors share too much sequence with their siblings to be resolved at your read length. Aggregate them to the gene level or use longer reads.

Finally, remember that the reference annotation is a moving target and that biological interpretation of RNA-seq is still an area where methods disagree. The reasonable response is to fix your annotation version, record every tool version and flag, and rerun the full pipeline when you add samples rather than merging old quantifications with new ones ⁷.

Footnotes

1. Paul A McGettigan. Transcriptomics in the RNA-seq era. Current Opinion in Chemical Biology, 2013. https://doi.org/10.1016/j.cbpa.2012.12.008 ↩

2. Chenxu Zhu, Sebastian Preissl, Bing Ren. Single-cell multimodal omics: the power of many. Nature Methods, 2020. https://doi.org/10.1038/s41592-019-0691-5 ↩

3. Ana Conesa, Pedro Madrigal, Sonia Tarazona, et al. A survey of best practices for RNA-seq data analysis. Genome Biology, 2016. https://doi.org/10.1186/s13059-016-0881-8 ↩ ↩²

4. Jiung-Wen Chen, Lisa Shrestha, George Green, et al. The hitchhikers’ guide to RNA sequencing and functional analysis. Briefings in Bioinformatics, 2023. https://doi.org/10.1093/bib/bbac529 ↩

5. Sang-Yoon Shin, Jin Seo Jeong, Jae Yun Lim, et al. Transcriptomic analyses of rice (Oryza sativa) genes and non-coding RNAs under nitrogen starvation using multiple omics technologies. BMC Genomics, 2018. https://doi.org/10.1186/s12864-018-4897-1 ↩

6. Geng Chen, Charles Wang, TieLiu Shi. Overview of available methods for diverse RNA-Seq data analyses. Science China Life Sciences, 2011. https://doi.org/10.1007/s11427-011-4255-x ↩

7. RNA-Seq Analysis: Methods, Applications and Challenges. Frontiers Research Topics, 2020. https://doi.org/10.3389/978-2-88963-705-8 ↩