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

Technical Summary: An End-to-End Reproducible RNA-Seq Workflow

Problem Statement
While RNA sequencing (RNA-seq) is a principal technology for transcriptome-wide analysis, the transition from raw sequencing reads to biologically interpretable results often lacks transparency and reproducibility. Although mature software exists for individual analytical steps, the integration of these tools into a fully reproducible, end-to-end workflow remains challenging. Issues such as undocumented manual processing, inconsistent software versions, and fragmented file management can hinder independent verification. Furthermore, moving from statistically significant gene lists to biologically meaningful interpretation often requires complementary analytical approaches (e.g., combining threshold-based differential expression with rank-based gene-set enrichment) that are not always integrated into standard pipelines.

Methodology
The authors developed a modular, version-controlled bulk RNA-seq workflow designed to span the entire analytical spectrum from data acquisition to biological interpretation. The workflow was implemented in a Linux environment using Conda/Miniforge for dependency management and Git for version control.

• Data Acquisition and Quality Control: Raw sequencing data (six single-end libraries from the GSE206185 dataset) were retrieved from the NCBI SRA using SRA Toolkit. Quality assessment was performed using FastQC, with metrics aggregated via MultiQC.
• Quantification and Summarization: A decoy-aware reference index was constructed using GENCODE v25 and the GRCh38 genome assembly. Transcript-level quantification was performed using Salmon (v2.4.1) with selective alignment. Transcript abundances were summarized to the gene level using tximport (v1.34.0), generating length-scaled count estimates.
• Differential Expression (DE): Statistical modeling was conducted using DESeq2 (v1.46.0). The workflow included low-expression filtering (retaining genes with ≥10 counts in ≥2 samples), normalization via size-factor estimation, and Wald testing. Log2 fold changes were shrunk using apeglm to improve effect-size estimation. Significance was defined by a Benjamini–Hochberg adjusted P value <0.05 and an absolute shrunken log2 fold change ≥1.
• Functional and Pathway Analysis: The workflow integrated multiple enrichment strategies:
  • Over-representation Analysis (ORA): Gene Ontology Biological Process (GO-BP) and KEGG pathway analyses were performed using clusterProfiler on significantly up- and downregulated genes.
  • Gene Set Enrichment Analysis (GSEA): Preranked GSEA was conducted using fgsea with MSigDB Hallmark gene sets, utilizing genome-wide ranked lists based on Wald statistics rather than thresholded gene lists.
• Reciprocal and Leading-Edge Analysis: A specific analytical module was developed to identify "reciprocal" pathways—those significantly enriched in opposite directions between the mutation state (Y537S vs. WT) and the treatment state (OTX015 vs. Y537S). Leading-edge genes contributing to these reciprocal signals were extracted and visualized.
• Reproducibility: All scripts, environment specifications (YAML), reference checksums, and source data matrices were preserved under Git version control and archived in Zenodo.

Key Contributions

1. A Reproducible Framework: The paper provides a complete, documented, and version-controlled computational workflow that transforms raw RNA-seq data into publication-quality figures and source data. It addresses the "black box" nature of many bioinformatics analyses by preserving the relationship between raw data, intermediate outputs, and final interpretations.
2. Integrated Analytical Strategy: The workflow demonstrates the utility of combining threshold-based DE analysis with rank-based GSEA and reciprocal pathway analysis to reveal coordinated transcriptional remodeling that might be missed by gene-list inspection alone.
3. Application to ESR1-Mutant Breast Cancer: The framework was applied to a specific biological context: ESR1 Y537S-mutant breast cancer cells treated with the BET inhibitor OTX015. This served as a worked example to validate the workflow's ability to detect mutation-associated signatures and treatment-induced remodeling.

Results

• Data Quality and Quantification: All six libraries passed quality control with high mapping rates (92.80%–94.13%) using the decoy-aware Salmon index. After filtering, 16,503 genes were retained for statistical analysis.
• Global Transcriptomic Structure: Principal Component Analysis (PCA) showed strong separation of the three experimental conditions (WT, Y537S, Y537S+OTX015), with PC1 and PC2 explaining 75.84% of the variance.
• Differential Expression:
  • Y537S vs. WT: 846 DEGs (401 up, 445 down).
  • Y537S+OTX015 vs. Y537S: 2,359 DEGs (959 up, 1,400 down), indicating extensive remodeling.
  • Y537S+OTX015 vs. WT: 2,539 DEGs, confirming that OTX015 treatment did not simply restore the wild-type transcriptome.
• Pathway Enrichment and Reciprocity:
  • Y537S State: Characterized by strong positive enrichment of estrogen-response (Early and Late) and MYC-target pathways, and negative enrichment of p53, hypoxia, and TNFα/NF-κB pathways.
  • OTX015 Treatment: Induced a reciprocal shift. Estrogen-response and MYC-target pathways became negatively enriched, while p53, hypoxia, and TNFα/NF-κB pathways became positively enriched.
  • Leading-Edge Genes: Specific genes (e.g., PDZK1, EGR3, AREG for estrogen response; CDKN2B, SAT1 for p53) exhibited strong reciprocal expression changes, providing gene-level support for the pathway-level findings.

Significance and Claims
The authors claim that this study establishes a transparent and reproducible framework for transforming raw bulk RNA-seq data into differential-expression results and pathway-level interpretations. The significance of the work lies in two areas:

1. Methodological: It provides a practical, version-controlled pipeline that integrates diverse tools (Salmon, DESeq2, clusterProfiler, fgsea) to ensure that analyses can be independently reproduced, audited, and extended.
2. Biological: The application to the ESR1 Y537S dataset demonstrates how integrated gene-level and enrichment analyses can reveal coordinated, reciprocal transcriptional changes. Specifically, the study shows that OTX015 treatment opposes the direction of major Y537S-associated transcriptional programs (estrogen response and MYC targets) and activates opposing programs (p53, hypoxia), without completely restoring the wild-type state.

The authors explicitly state that while the workflow generates experimentally testable hypotheses regarding the transcriptional response of ESR1-mutant cells to BET inhibition, the computational reanalysis alone does not prove direct mechanistic causality (e.g., direct inhibition of mutant ESR1 by OTX015) or cellular phenotype reversal. They emphasize that the findings represent selective reversal of major transcriptional programs rather than complete normalization, and that independent biological validation is required to confirm these hypotheses.