Bulk RNA-seq analysis: strand-aware quantification from raw FASTQs to DESeq2 contrasts

At a glance

What is bulk RNA sequencing (RNA-seq)?

Bulk RNA-seq aligns or quasi-maps short reads to a reference genome or transcriptome and aggregates signal into a sample-by-gene count matrix. Each library reports average expression across many cells rather than per-cell counts [1]. scRNA-seq instead assigns UMIs to individual barcodes. Compared with microarrays, RNA-seq can capture novel splice junctions and unannotated transcripts when depth permits. The GTEx v8 resource sequenced 17,382 RNA-seq samples across 54 tissue sites from 948 donors and is widely used as a tissue-level expression atlas [4]. Pepkio starts from FASTQs or existing count matrices and returns QC-audited DE results with documented parameters at each stage. Custom inputs, deliverable formats, and analyses outside the standard workflow are agreed at kickoff. See the bulk RNA-seq glossary.

When should you use bulk RNA sequencing (RNA-seq)?

Bulk RNA-seq is a fit when biological variation is at the sample or condition level, such as drug response, genotype, disease stage, or treatment arms, and tissue composition is fairly stable or deconvolution is acceptable.

• Cross-tissue regulation: GTEx v8 profiled 17,382 samples across 54 tissues, linking cis-eQTLs to tissue-specific gene regulation [4].
• Primary vs. metastatic melanoma: Sun et al. [5] analyzed TCGA bulk RNA-seq from 103 primary and 368 metastatic melanoma samples, identifying 246 differentially expressed lncRNAs and 856 differentially expressed mRNAs [5].
• Hypoxia-driven plasticity: Guo et al. [6] used RNA-seq in LNCaP and PC3 prostate cancer lines to show ONECUT2 and hypoxia synergistically activate neuroendocrine marker programs.

How the analysis works — step by step

1. 1. Validate inputs and sample metadata

   Pepkio confirms FASTQ integrity (or count-matrix dimensions), read layout, library type (poly(A)+, rRNA-depleted, strand-specific), and experimental design. Sample IDs, batch, condition, and covariates are recorded in sample_manifest.csv. Designs with fewer than three biological replicates per condition are flagged before DE testing [1].

   Tools and outputs

   Tools used: FastQC 0.12.1 / custom validation scripts

   Output: sample_manifest.csv with library IDs, read counts, strandness flags, and QC notes

2. 2. QC and trim raw reads

   Adapter contamination, low-quality tails, and overrepresented sequences are assessed per library. When trimming is warranted, reads are filtered before alignment. Aggregated metrics are compiled for review [7].

   Tools and outputs

   Tools used: FastQC 0.12.1; fastp 0.24.0; MultiQC 1.25.1

   Output: Per-sample FastQC/fastp reports; multiqc_report.html

3. 3. Align reads to the reference genome

   Paired- or single-end reads are aligned with splice-aware mapping to GRCh38 or GRCm39 with GENCODE annotation. Strand-specific library orientation is set from metadata or inferred from RSeQC when ambiguous [8]. Mapping rates, uniquely mapped reads, and splice junction counts are compared against expected ranges for the library type.

   Tools and outputs

   Tools used: STAR 2.7.11b; samtools 1.21

   Output: Coordinate-sorted, indexed BAM files; star_log_final.txt per sample

4. 4. Quantify gene-level counts

   Gene-level read counts are summarized from aligned BAMs using exon-gene overlap rules consistent with the annotation build. Multi-mapping reads are handled per Subread defaults; parameters are recorded for reproducibility [9].

   Tools and outputs

   Tools used: Subread/featureCounts 2.0.8

   Output: Raw count matrix (counts_raw.csv); featureCounts_summary.txt

5. 5. Audit count matrix quality

   Sample-level library sizes, gene detection rates, and pairwise correlations are computed. PCA and hierarchical clustering are used to separate biological condition from sequencing batch when possible. RSeQC reports read distribution over genomic features and confirms strandness [8]. Outlier samples are flagged before modeling.

   Tools and outputs

   Tools used: DESeq2 1.52.0 (diagnostics); RSeQC 5.0.3

   Output: sample_qc_summary.csv; PCA and correlation heatmaps; RSeQC read-distribution plots

6. 6. Filter low-count genes

   Genes with insufficient counts across samples are removed using a pre-specified rule, typically ≥10 counts in at least n samples, where n equals the smallest replicate group size [3]. Retained and excluded gene counts are logged.

   Tools and outputs

   Tools used: DESeq2 1.52.0

   Output: Filtered count matrix; gene_filter_log.csv

7. 7. Model and test differential expression

   A design matrix encodes condition, batch, and covariates (e.g., sex, RIN, paired blocking). DESeq2 estimates size factors, dispersions, and Wald statistics for each contrast [3]. Paired designs use ~ subject + condition; multi-factor designs include interaction terms when biologically justified.

   Tools and outputs

   Tools used: DESeq2 1.52.0

   Output: DESeqDataSet .rds object; dds_results_<contrast>.csv

8. 8. Apply shrinkage and multiple-testing correction

   Log2 fold-changes can be shrunk with DESeq2 lfcShrink using apeglm or ashr on request. Benjamini-Hochberg adjusted p-values from DESeq2 control the false discovery rate across tested genes [3]. Contrasts with no genes passing a pre-agreed FDR threshold are reported with full ranked tables.

   Tools and outputs

   Tools used: DESeq2 1.52.0; Bioconductor apeglm on request

   Output: deg_results.csv with columns: gene_id, gene_name, baseMean, log2FoldChange, lfcSE, stat, pvalue, padj

9. 9. Run pathway and GO enrichment

   Significant genes (padj threshold agreed at kickoff) are tested for over-representation in GO and KEGG gene sets. Dot plots and bar charts summarize enriched terms with gene counts and adjusted p-values [10].

   Tools and outputs

   Tools used: clusterProfiler 4.20.0

   Output: go_enrichment.csv, kegg_enrichment.csv; enrichment dot plots

10. 10. Package figures, scripts, and Methods draft

    MA plots, volcano plots, and top-gene heatmaps are exported at publication resolution. Commented R scripts reproduce each stage from raw inputs. A journal-formatted Methods section cites software versions and reference builds. Salmon 1.10.3 + tximport 1.34.0 quantification is available when clients provide transcriptome references or prefer quasi-mapping without full BAM alignment [11], scoped at kickoff.

    Tools and outputs

    Tools used: DESeq2 1.52.0; pheatmap; EnhancedVolcano on request

    Output: PDF/SVG figures; R scripts; README; Methods draft; final deliverable bundle

What Pepkio delivers

Technical decisions we make — and why

Alignment + counting: STAR 2.7.11b + featureCounts

Splice-aware genome alignment with explicit exon-gene overlap counting; Salmon + tximport on request for quasi-mapping workflows [9, 11, 12].

DE engine: DESeq2 median-ratio normalization + Wald test

Operates on raw counts with a negative-binomial model; edgeR QL F-test or limma-voom on request [3, 13, 14].

Batch handling: include batch in the design formula

Preferred over post-hoc correction when batch is known and not confounded with condition; ComBat-seq or surrogate variable analysis scoped separately when the design requires it [3].

Gene filtering: ≥10 counts in ≥n samples

Removes genes with near-zero signal before dispersion estimation; DESeq2 independent filtering retained for the test step [3].

Enrichment: clusterProfiler ORA on padj-filtered gene lists

GSEA or ReactomePA extensions on request [10].

Common questions

Related services

• Single-cell RNA-seq — Cell-type resolution when tissue heterogeneity makes bulk averages insufficient.
• Spatial transcriptomics — Retain tissue architecture and microenvironment context that homogenized bulk libraries lose.
• Long-read RNA-seq — Full-length isoform discovery beyond short-read exon counting.
• Experimental design — Replicate planning and power estimation before library prep.
• Custom consulting — Assay and depth selection before sequencing.