Discovery of Oncogenic Mediator Genes in Rectal Cancer Chemotherapy Response using Gene Expression Data from Matched Tumor and Patient-Derived Organoid

Integration of gene expression data from matched RC tumors and PDOs

In this study, we employed a multi-tiered, proof-of-principle integration strategy to unravel the complex molecular landscape of rectal cancer and its response to chemotherapy as illustrated in Fig 1. Our rationale is that 5-fluorauracil-based chemotherapy has been and remains the backbone of chemotherapy and chemoradiation in rectal cancer for the past 60 years – however, no new therapies have been added that specifically address the patients with chemoresistant tumors[3, 24–27]. Our approach began with a comprehensive analysis of transcriptome data from 18 RC tissue samples and matched patient-derived organoids (PDOs), focusing on the identification of drug-related genes based on IC50 outcomes according to FOLFOX (5-fluorouracil, leucovorin, and oxaliplatin) treatment of RC PDOs from our RC biorepository. IC50 in RC PDOs has been shown to correlate well with both oncologic and clinical outcomes [22, 28]. The integration extended to a thorough examination of The Cancer Genome Atlas (TCGA) COAD/READ data, where differential expression analyses (DEAs) were conducted to distinguish between tumor and paired normal tissues, uncovering key differentially expressed genes (DEGs). A subsequent gene regulatory network (GRN) analysis provided a structural and functional framework for understanding the interactions and regulatory relationships among these genes. Central to our strategy was the application of Moonlight, a tool adept at analyzing over a hundred biologically relevant processes, to identify altered oncogenic processes in tumors. This analysis pinpointed specific oncogenic mediator genes acting as regulatory hubs within these processes. By amalgamating the insights gleaned from RC tumors and PDOs, TCGA data, and the process-specific gene evaluation via Moonlight, we established a robust, multidimensional approach. This approach not only enhanced our understanding of the oncogenic pathways in rectal cancer but also illuminated potential oncogenic mediators, offering putative targets for sensitizing resistant tumors to chemotherapy.

• Download figure
• Open in new tab

Figure 1

The workflow of the study. It consists of (A) tissue-organoid analysis for screening genes associated to drug-response, (B) Moonlight-based analysis for discovering oncogenic regulator-regulons groups of colorectal cancer, and (C) the integrated analysis for identifying oncogenic mediators as potential indicators of drug responsiveness.

Gene expression consistency between cancer tissues and matched PDOs

In our analysis, all 50 samples, comprising 18 tissue samples and 32 organoid samples, met the stringent criteria set for quality control and were subsequently included in further analysis (put these criteria in Supplementary Material?). Following the removal of genes with minimal expression (detailed in the Methods section), a substantial pool of 17,418 genes was retained for analysis. Initial evaluations focused on the consistency of gene expression within organoid samples derived from individual patients (Fig. 2A). For each set of three replicated organoid samples per patient, a remarkable consistency was observed, with pairwise Spearman’s correlation coefficients invariably exceeding 0.958 (P-value < 0.001), underscoring the reliability of organoid models in reflecting patient-specific gene expression profiles.

• Download figure
• Open in new tab

Figure 2

Reproducibility of patient-derived organoid and tissue-organoid correlation. (A) Sample-wise correlation based on 21 organoid samples derived from seven patients. (B) Principal component analysis (PCA) of 18 tissue and 32 matched organoid samples based on gene expression data. (C) Histograms of gene-wise based on expression data of all genes, and (D) sample-wise correlation coefficients based on all genes (blue) and top-5000 variable genes selected from organoid (orange) of 18 tissue-organoid matched patients.

Further, Principal Component Analysis (PCA) was employed on the gene expression data from the 50 samples, revealing a predominant segregation into two distinct clusters, one encompassing tissue samples and the other comprising organoid samples (Fig. 2B). Notably, the organoid cluster exhibited reduced variability along the first principal component (PC1), but heightened variability along the second principal component (PC2), suggesting a nuanced representation of genetic diversity within these models. Moreover, replicated organoid samples from the same patients were observed to cluster tightly on the PCA plot, reinforcing the reproducibility of the organoid model.

We then probed the correlation in gene expression between primary tumor tissues and their derivative organoid samples (Fig. 2C). In scenarios involving multiple organoid samples from a single patient, the gene expression for the organoids was represented by the arithmetic mean of the samples. A substantial proportion of genes, 16,918 out of 17,418 (97.1%), demonstrated a positive correlation (Spearman’s rho > 0) between tissue and organoid gene expressions. Within this subset, 6,514 genes manifested a statistically significant correlation, marked by a Spearman’s rho > 0.5 and a P-value < 0.05 (Supplementary Table 1). A patient-specific analysis was further conducted, assessing the gene expression correlation between tissue samples and their corresponding organoid samples, considering the full spectrum of genes as well as a subset of highly variable genes (Fig. 2D). This analysis consistently revealed high correlation coefficients across all patient samples, affirming the organoid model’s capacity to faithfully recapitulate the gene expression landscape of the original tumor tissues.

Gene selection for drug responsiveness

To discern genes correlated with drug response, we executed regression analyses on tumor tissues and their corresponding organoids, paired with IC50 values from FOLFOX treatment (see Methods section). This analysis identified a total of 660 genes from tissues (refer to Supplementary Table 2) and 819 genes from organoids (refer to Supplementary Table 3), each exhibiting a significant correlation with drug response (p-value < 0.05). Notably, there was an overlap of 75 genes between the two datasets (Fig. 3A and itemized in Supplementary Table 4). A subsequent sample-wise correlation assessment demonstrated a high level of consistency in the expression of these drug-related genes between tissue and organoid samples, as depicted in Fig. 3B. Of the 75 drug response genes (DRGs) highlighted, 17 were positively linked to drug resistance, while 58 were positively linked to drug sensitivity, as shown in Fig. 3C. Utilizing the expression profiles of these DRGs, unsupervised hierarchical clustering was performed on a combined dataset of tissue and organoid samples. This clustering effectively segregated the samples into two distinct groups, each exhibiting contrasting patterns of drug responsiveness. Furthermore, recognizing the limited sample size in this analysis and aiming for broader applicability in integrated analyses, we also compiled a supplementary list of 458 genes. These genes exhibit a marginally significant correlation with drug response in both tissue and organoid samples, as determined using a relaxed threshold of p-value 0.15 (see Supplementary Table 5).

• Download figure
• Open in new tab

Figure 3

Identification of drug responsiveness genes. (A) Venn’s diagram of drug-related genes (P < 0.05) identified from organoid (green, n = 819) and tissue (yellow, n = 660) samples. (B) Histograms of sample-wise correlation coefficients based on top-5000 variable genes (orange) and 819 drug-related genes (blue) identified from organoid samples. (C) Heat map of gene expression of 75 drug-related genes (P < 0.05 in both tissue and organoid samples), sorted by descending order based on log-IC50 values of Folfox treatment in sample-wise.

Detection of oncogenic regulator-regulons through Moonlight

In our study, we leveraged an enhanced Moonlight workflow to systematically screen for potential oncogenic mediators from the combined TCGA dataset of rectal (READ) and colon (COAD) cancers, acknowledging the relatively limited sample size in each category (Fig. 4A). A cohort of 223 patients, each with matched primary tumor and adjacent normal tissue samples, and documented chemotherapy history, formed the basis of our analysis. Post exclusion of microRNA data, the remaining 17,082 protein-coding genes underwent a Differential Phenotype Analysis (DPA) using Moonlight, which yielded 5,272 genes demonstrating significant differential expression (log2|FC| > 1 and FDR < 0.01, Supplementary Table 6). Subsequent Gene Regulatory Network (GRN) analysis within Moonlight scrutinized these differentially expressed genes (DEGs) as potential regulators, exploring their associated regulons from the entire gene pool of 17,082 candidates. This comprehensive analysis identified 5,199 genes that govern at least one regulon (Supplementary Table 7), with the most influential gene regulating up to 783 regulons.

• Download figure
• Open in new tab

Figure 4

Identification of potentially oncogenic mediators via (A) Moonlight and (B) integrated analysis. (C) Boxplots of six genes as potential indicators for drug response, and (D) four genes out of the six showing significant association with recurrence-free survival (RFS) among colon cancer patients (n = 1,296) from multiple datasets.

To delve deeper into the biological significance of these findings, we conducted modified functional enrichment analysis (FEA) and upstream regulator analysis (URA) using the unique modules from Moonlight, pinpointing regulator-regulon groups implicated in cancer-relevant biological processes (BPs). Initially, Moonlight offered insights into 101 BPs; however, we honed our focus to 25 BPs intimately associated with cancer progression, such as “apoptosis,” “tumor cell line proliferation,” and “tumor growth” (Supplementary Table 8). These BPs were dissected further using Ingenuity Pathway Analysis (IPA) to isolate a subset of genes with strong associations to each process. An automated PubMed literature review annotated these genes, hinting at the direction of their expression alterations during these processes. In this refined FEA and URA, our criteria for significance were stringent: we considered only regulator-regulon groups with a Moonlight |Z-score| over 1 and a Fisher’s exact test P-value below 0.05. This approach led us to identify 2,118 regulator-regulon groups, each potentially playing a pivotal role in cancer-related BPs (Supplementary Table 9), marking a crucial step in understanding the regulatory underpinnings of oncogenic processes in rectal and colon cancers.

Integration of the discoveries from Moonglight and tissue-organoid analyses

In our pursuit to discern specific associations between the 2,118 cancer-related regulator-regulon groups and drug responsiveness, we embarked on an integrated analysis, as depicted in Fig. 4B. This process commenced with a Gene Set Enrichment Analysis (GSEA), wherein the regulator-regulon groups were examined alongside the insights derived from our tissue-organoid investigation. Employing the 2,118 regulator-regulon groups as gene sets, the GSEA utilized two pre-ranked gene lists, each organized in descending order based on the regression’s t-value, one for tissue and the other for organoid analyses. To enhance the robustness of our findings, we focused on 11,198 genes that consistently indicated drug responsiveness trends across both tissue and organoid samples. This GSEA identified 334 regulator-regulon groups that exhibited significant enrichment across both modalities, leading us to classify these 334 entities as oncogenic mediators (OMs). In our analysis, six genes – BTK, NCKAP1L, LAX1, CDK1, RAD51AP1, and NAT2 – emerged as noteworthy for their association with drug responsiveness. Although the threshold for significance was set at a more inclusive level (P-value < 0.15), these genes consistently showed a positive correlation with sensitivity to FOLFOX treatment in both tissue and organoid samples, as demonstrated in Fig. 4C. This suggests a potential link between these genes and chemotherapy response.

To further determine if these genes were associated with an oncologic outcome, we tested the prognostic potential of these six OMs using the Kaplan-Meier Plotter database, with results illustrated in Fig. 4D. This analysis revealed that RAD51AP1, NAT2, and LAX1 were significantly associated with extended recurrence-free survival (RFS), showcasing HR values of 0.54 (logrank P = 6.5e-09), 0.63 (logrank P = 2.9e-05), and 0.70 (logrank P = 1.8e-03), respectively. These findings align well with the drug responsiveness trends observed in our integrated analysis. Interestingly, NCKAP1L, despite showing a significant association with RFS (HR = 1.36, logrank P = 5.5e-05), exhibited a trend that diverged from the drug response analysis.

To delve into the therapeutic prospects of the identified drug-related genes, we assessed the possibility of targeting specific perturbagens against these genes through the Connectivity Map tool. This exploration led to the identification of 89 perturbagens targeting 24 out of the 334 oncogenic mediators and 107 drugs targeting 23 out of the 458 drug-related genes (p < 0.15), as highlighted in our previous tissue-organoid analysis (Fig. 5A and 5B; Supplementary Table 11 and 12). This revelation underscores the potential for therapeutic intervention and paves the way for further investigation into the efficacy of these perturbagens in modulating drug responsiveness in rectal cancer.

• Download figure
• Open in new tab

Figure 5

Heatmap showing availability of perturbagens from the Connectivity Map in columns and their target drug-related genes in rows from (A) Moonlight-based oncogenic mediators and (B) tissue-organoid analysis with P < 0.15.

Human RC tissue and patient-derived organoid

A total of 18 patients diagnosed with rectal cancer (RC) at Memorial Sloan Kettering Cancer Center were selected in our study. For each patient, one tumor tissue was obtained before treatment, and at least one matched patient-derived organoid (PDO) samples were generated. In specific, for each of the eleven out of 18 patients, one single PDO was generated, while to validate the replicability of PDO model, each of the remaining seven patients had triplicated organoid samples. Endoscopic tumor biopsies obtained with biopsy forceps were obtained from RC patients undergoing standard of care diagnostic biopsy or on-treatment biopsy in the outpatient setting. The specimens are collected per protocol and tissue was prioritized as follows: sample sent to pathology, sample for PDO, and remaining saved for banking. Tissues are chopped into ∼1-mm pieces in cold PBS + antibiotics, and after tissue processing, tumor cells are resuspended in 300 μL of Matrigel and plated onto 24-well culture plates (40 μL droplet per well). After the gel solidifies, 700-800 μL of culture medium is added to wells and incubated at 37°C. PDOs are established and maintained in Matrigel as described in Ganesh et al.[22]. Consequently, 18 tissues and 32 matched organoids were involved in this study and were then RNA-sequenced (see below).

Preparation and RNA-sequencing of tissue and organoid samples

Surgically resected rectal tissue or biopsy tissue was washed with ice-cold PBS-Abs buffer (phosphate-buffered saline with antibiotic-antimycotic (Gibco™), gentamicin (Gibco™) and plasmocin (Invivogen). The tissue was divided into two parts. One part was flash-frozen in liquid nitrogen and then reserved at -80 °C for later purposes. The other part was processed for tumoroid culturing purposes, as reported[22]. The patient-derived tumoroids were cultured and expanded in Matrigel for three passages. To RNA, well-growing tumoroids in 4-6 matrigel domes (40ul matrigel/dome) were dissociated with TrypLE Express (Gibco™) as reported[22]. The obtained pellet was further washed with 5 ml of ice-cold PBS. Then, RNA was extracted from the tumoroid pellet using the AllPrep RNA kit (Qiagen) following the manufacturer’s instructions. If the RNA needs to be extracted separately, then the PureLink™ RNA Mini Kit (Invitrogen™) was used. The tissue RNA was extracted from the flash-frozen specimen using the AllPrep RNA kit.

For RNA sequencing, total RNA was prepped with the Nugen Universal Plus mRNA-Seq (M01442 v2) using 750ng via Qubit and 14 cycles PCR. Libraries were sequenced on Illumina Novaseq 6000 and ∼50M PE150 were generated per sample.

Tissue-organoid data analysis

To process raw FASTQ-formatted RNA sequencing data from organoid and tissue samples, we utilized the nf-core/rnaseq tool[29] (version 3.12.0) in Nextflow pipeline[30] (version 22.10.7) with default settings and parameters. The data processing workflow consisted of pre-processing quality control (QC), STAR-based reads alignment to the reference genome (Genome assembly GRCh38: https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001405.26/), and Salmon-based gene expression data quantification. Detailed information of RNA-seq analysis pipeline of Nextflow could be referred to at https://nf-co.re/rnaseq/3.12.0.

Data analysis for count-based gene expression data were then performed using R software[31] (version 4.2.2) in RStudio[32] (version 2023.09.0). We transferred count data to count-per-million (CPM) using edgeR package[33] (version 3.40.2), and retained only unduplicated genes with CPM values simultaneously greater than 0.1 among at least 15% tissue samples, 15% organoid samples and 15% tissue plus organoid samples. We implemented limma-based standardization using calcNormFactors function from edgeR package, and applied a log2 transformation to the CPM values, resulting in normalized data for the downstream analysis.

To validate the reproducibility of the triplicated organoid samples, we created heatmaps to visually represent the Spearman’s correlation matrix. To evaluate the correlation between gene expression in organoid and tissue samples, we computed patient-wise Spearman’s correlation coefficients. For patients with replicated organoid samples, we averaged their gene expression values to align with tissue values. Histograms were then generated to illustrate the distribution of Spearman’s correlation coefficients for all genes and top-5000 variable genes.

To identify genes associated with drug responsiveness (Fig. 1A), we conducted a preliminary screening through differential expression analysis (DEA) using regression models. As multiple replicated samples were generated for a proportion of organoids, we treated the patient identifier as a mixed-effect variable for regression modeling. Specifically, we fitted the model log₂(CPM + 1) ∼ 1|patient.id + ln(IC₅₀) using the lmer function in lme4 R package for organoid samples. For tissue samples, the model log₂(CPM + 1) ∼ ln(IC₅₀) was fitted using the limma-voom method from the limma package [34]. Here, gene expression was measured by log2-transformed CPM. Genes with unadjusted P-values < 0.05 were selected as candidates associated with chemotherapy response.

Furthermore, we identified the intersection of candidate genes from both organoid and tissue data and used them to generate a heatmap utilizing the ComplexHeatmap[35] package.

A modified Moonlight workflow

To identify a set of oncogenic genes based on TCGA cohort[36], we performed customized analysis based on Moonglight pipeline[37] (Fig. 1B). For all human colon (COAD) and rectal (READ) cancer patients from TCGA, high-throughput RNA-seq data normalized by upper quartile (HTseq – FPKM-UQ) and clinical information were downloaded and processed with TCGABiolinks package[38]. To prioritize the investigation of drug response, only patients who received chemotherapy were included in the analysis. For the patients with duplicated primary solid tumor (TP) samples, only the first sample profiled by date was retained. For consistency with our dataset, we excluded genes that encode microRNAs from subsequent analysis.

The MoonlightR::DPA module was employed to identify differentially expressed genes (DEGs) between tumor and normal tissues, with a significance threshold set at false discovery rate (FDR) < 0.01 and log₂|FC| > 1. To uncover potential target regulons associated with each DEG, the MoonlightR::GRN module was utilized to construct a gene regulatory network (GRN) based on mutual information. This GRN comprised regulator-regulon groups, where each DEG was considered a potential regulator, and its regulons were selected from the entire gene pool using algorithms in the GRN module. To further investigate genes associated with cancer-related biological processes (BPs), a subset of BPs was chosen from the original 101-BP list predefined in the Moonlight R package. Within each BP, there was a set of genes that could potentially exhibit upregulation, downregulation, or unspecified changes in expression levels during that specific BP. A selective subset of cancer-related BPs was retained for subsequent analysis in MoonlightR::URA. In this step, the enrichment significance of each combination of regulator-regulon group and BP was assessed using URA-derived Moonlight Z-scores and P-values. Regulator-regulon groups with at least one significantly enriched BP, meeting the threshold of Mooonlight |Z| > 1 and P < 0.05 in the URA step, were selected as candidates for further analysis.

Gene set enrichment analysis (GSEA) using the fgsea[39] package was performed for each regulator-regulons group to integrate the Moonlight-based analysis and tissue-organoid analysis (Fig. 1C). A set of regulons, representing the potentially target genes for each regulator, was utilized as a gene set. Two gene lists, one from organoid analysis and the other from tissue analysis, were employed in GSEA. These lists were ranked in descending order of log₂FC. Only genes exhibiting the same trends in drug response were retained for GSEA in both pre-ranked lists. The gene set size was constrained, with a maximum of 500 genes and a minimum of 5 genes.

Regulators from the regulator-regulons groups that demonstrated significant enrichment patterns (P-value < 0.05) in both organoid-based and tissue-based lists were ultimately selected as oncogenic mediators (OMs) through this integrated analysis.

Survival analysis

To investigate the prognostic performance of genes, we retrieved an online resource namely Kaplan-Meier Plotter (https://kmplot.com/analysis/index.php?p=service&cancer=colon), which provided survival analysis for 1,296 patients with colon cancer from 16 datasets based on their gene expression and survival information [40]. For each gene, patients were divided into two groups, using an automatic classification with the optimal cutoff in gene expression. This method aimed to achieve the lowest log-rank p-value for recurrence-free survival (RFS) between the groups, and based on which, a Kaplan-Meier (KM) plot was generated for each gene. When multiple probes corresponded to the same gene, the probe with the smallest p-value was selected as the representative.

Connectivity Map analysis

We used Connectivity Map[41] online tools to retrieve potential compounds and perturbations for genes of interest. A text file namely compoundinfo_beta.txt containing a full list of perturbagens and their target genes was downloaded via https://clue.io/data/CMap2020#LINCS2020.