An End-to-End Reproducible RNA-Seq Workflow from Raw Sequencing Reads to Differential Expression, Pathway Enrichment, and Biological Interpretation

<title>Abstract</title> <p> RNA sequencing (RNA-seq) enables comprehensive characterization of transcriptional states, but reproducible analysis from raw sequencing reads to biologically interpretable results requires integration of multiple computational tools, reference resources, statistical procedures, and visualization methods. Although mature software exists for individual analytical steps, implementation of a transparent end-to-end workflow can remain challenging, particularly when analyses must be fully reproducible and accompanied by publication-quality figures and source data. We developed a modular and version-controlled bulk RNA-seq workflow spanning raw sequencing data acquisition, quality assessment, transcript quantification, gene-level summarization, differential-expression analysis, exploratory visualization, functional enrichment, and gene-set interpretation. The workflow integrates SRA Toolkit, FastQC, MultiQC, a GENCODE v25/GRCh38 decoy-aware Salmon reference, Salmon, tximport, DESeq2, apeglm, clusterProfiler, fgsea, MSigDB Hallmark gene sets, and R-based visualization. A publicly available breast-cancer RNA-seq dataset comprising six single-end libraries representing MCF-7 wild-type cells, Y537S-mutant cells, and Y537S-mutant cells treated with the BET inhibitor OTX015 was analyzed as a worked example. Differentially expressed genes were defined using a Benjamini–Hochberg adjusted P value &lt; 0.05 and absolute shrunken log ₂ fold change ≥ 1. Gene Ontology Biological Process (GO-BP) and KEGG over-representation analyses were complemented by preranked Hallmark gene-set enrichment analysis. Salmon successfully quantified all six libraries, with mapping rates ranging from 92.80% to 94.13%. Transcript-level estimates were summarized to 57,627 genes, of which 16,503 passed expression filtering and were retained for statistical analysis. Principal component analysis showed strong separation of the three experimental conditions, with PC1 and PC2 explaining 44.51% and 31.33% of total transcriptomic variance, respectively, accounting jointly for 75.84%. ESR1 Y537S-mutant cells differed from wild-type cells at 846 genes, including 401 upregulated and 445 downregulated genes. OTX015-treated Y537S cells differed from untreated Y537S cells at 2,359 genes, including 959 upregulated and 1,400 downregulated genes, whereas comparison of OTX015-treated Y537S cells with wild-type cells identified 2,539 differentially expressed genes. Functional analyses demonstrated broader pathway-level remodeling following OTX015 exposure. Hallmark GSEA identified 11 significant gene sets for Y537S mutant versus wild type, 26 for OTX015-treated Y537S versus wild type, and 20 for OTX015-treated versus untreated Y537S cells at FDR &lt; 0.05. Seven Hallmark programs showed significant reciprocal enrichment between the mutation and treatment comparisons. Early estrogen response changed from NES 3.36 in Y537S versus wild type to NES − 2.32 after OTX015 treatment; late estrogen response showed a similar reversal. MYC Targets V2 changed from NES 2.32 to -2.34, whereas the p53 pathway changed from NES − 1.51 to 2.06. Leading-edge analysis further demonstrated coordinated reciprocal behavior among genes contributing to these pathways. This workflow provides a transparent and reproducible framework for transforming raw bulk RNA-seq data into differential-expression results, pathway-level interpretation, and publication-quality outputs. Its application to ESR1-mutant breast-cancer data demonstrates how integrated gene-level and enrichment analyses can reveal coordinated treatment-associated transcriptional changes. </p>