Survival rates of patients with metastatic castration-resistant prostate cancer (mCRPC) are low due to lack of response or acquired resistance to available therapies, such as abiraterone (Abi). A better understanding of the underlying molecular mechanisms is needed to identify effective targets to overcome resistance. Given the complexity of the transcriptional dynamics in cells, differential gene expression analysis of bulk transcriptomics data cannot provide sufficient detailed insights into resistance mechanisms. Incorporating network structures could overcome this limitation to provide a global and functional perspective of Abi resistance in mCRPC. Here, we developed TraRe, a computational method using sparse Bayesian models to examine phenotypically driven transcriptional mechanistic differences at three distinct levels: transcriptional networks, specific regulons, and individual transcription factors (TF). TraRe was applied to transcriptomic data from 46 patients with mCRPC with Abi-response clinical data and uncovered abrogated immune response transcriptional modules that showed strong differential regulation in Abi-responsive compared with Abi-resistant patients. These modules were replicated in an independent mCRPC study. Furthermore, key rewiring predictions and their associated TFs were experimentally validated in two prostate cancer cell lines with different Abi-resistance features. Among them, ELK3, MXD1, and MYB played a differential role in cell survival in Abi-sensitive and Abi-resistant cells. Moreover, ELK3 regulated cell migration capacity, which could have a direct impact on mCRPC. Collectively, these findings shed light on the underlying transcriptional mechanisms driving Abi response, demonstrating that TraRe is a promising tool for generating novel hypotheses based on identified transcriptional network disruptions.

Significance:

The computational method TraRe built on Bayesian machine learning models for investigating transcriptional network structures shows that disruption of ELK3, MXD1, and MYB signaling cascades impacts abiraterone resistance in prostate cancer.

Prostate cancer is among the most frequently diagnosed male cancers in the world (1). Primary treatments, including surgery or radiotherapy, can achieve a cure for many of the afflicted men (2). However, a substantial portion is diagnosed with advanced-stage disease or experience disease recurrence. Although the mainstay of treatment for these cases has been androgen depletion therapy (ADT), data from clinical trials in patients newly diagnosed with metastatic prostate cancer (mPC; refs. 3, 4) have demonstrated a survival advantage of 37% for those receiving ADT in combination with other drugs (e.g., enzalutamide, abiraterone acetate, apalutamide, and docetaxel) over ADT alone. In addition, even though ADT is initially effective, >95% of patients with mPC on this treatment relapse, indicative of castration-resistant prostate cancer (CRPC) development (5), which has a poor 5-year survival. Because of significantly improved survival rates, the first-line therapy for patients with CRPC then becomes either the CYP17A1 inhibitor abiraterone combined with prednisone (AA/P; ref. 6), or the androgen receptor (AR) inhibitor enzalutamide (7). Although results are encouraging, trials with abiraterone (Abi) or enzalutamide treatment have highlighted two major challenges: (i) pre-existing mechanisms of resistance (primary resistance) preclude responses for nearly half of patients with CRPC, and (ii) resistance can develop rapidly in initial responders (acquired resistance; ref. 8).

The interactions between genes and their corresponding pathways drive various cellular functions that are critical in tumor development and response to therapy (9). These relationships, which are typically encoded as networks or modules of genes, provide a concise representation of the transcriptional landscape of the cell (10, 11). It is well established that functional transcriptional networks and their products change in response to different conditions, such as cellular DNA damage or environmental stress (12–14). Hence, the construction and exploration of the topology of these networks and their constituents are compelling approaches to developing and understanding biological mechanisms. In the CRPC context, the AR transcription factor (TF), together with many other cofactors, such as the steroid receptor coactivator/p160 family or the AR-specific MAGE-A11 coactivator (LxxLL and FxxLF motif-containing coactivators, respectively), regulates a series of downstream genes (e.g., prostate-specific antigen or transmembrane protease serine 2) upon ligand binding, followed by translocation into the nucleus (15). Recently, it was also noted that there are specific cofactors in CRPC tumors not present in normal tissue, which can guide AR to specific gene regions to regulate gene transcription (16, 17). In addition, mutations or tumor-specific alterations in pathways could also significantly affect transcriptional networks and downstream signaling. For example, the Wnt pathway, which is critical in differentiation and cell fate determination, has recently been associated with Abi resistance (18).

Another of the main questions is how cells respond to drugs. Although differential gene expression analysis evaluates changes in gene expression under different conditions, the incorporation of network structures and differential network analysis can provide insights that are mechanistically grounded (19). Specifically, signatures of differentially expressed genes are a downstream effect of global cell de-regulation in different response groups, that is, they are the most evident output of more complex system changes or states (20). As such, cellular phenotypes result from the different but precise combinations of activities of large networks of co-regulated genes in which TFs play a key role. Thus, solely analyzing the differential expression of individual genes (without the network context) would miss important information about the cellular system (21–23). For example, both Hudson and colleagues (20) and Califano and colleagues (22) show that master regulators or causal effectors are rarely differentially expressed. Furthermore, TFs are lowly expressed compared with other types of genes and their activity is not directly linked with their relative abundance (24). In addition, Califano and colleagues (22) and Sinha and colleagues (23) showed that gene regulatory networks (GRN) are a powerful tool to describe different cellular states through the reverse engineering of the cells’ transcriptional programs in cancer (22) and in developmental and behavioral biology (23) settings. In this context, drug treatments can activate different functional pathways in tumors based on differences in their underlying transcriptional networks (25, 26). Thus, identifying significant changes among networks from different response groups, also referred to as network rewirings, can help discover novel molecular diagnostic and prognostic signatures (13). As an example resulted from differential network analysis, RNA levels of the prostate cancer biomarker AMACR were discovered to have a positive correlation with the tumor-suppressor gene PTEN in adjacent normal tissue that was no longer present in prostate tumor samples (27).

To mechanistically understand how the differences in transcriptional networks may contribute to Abi response in metastatic CRPC (mCRPC), we proposed and developed TraRe (28), a computational method already available at Bioconductor (RRID: SCR_006442) as an R package, that enables the broad biomedical community to uncover complex transcriptional dynamics from RNA sequencing (RNA-seq) data. The unique value of TraRe is that it attempts to examine phenotypically driven mechanistic differences at three distinct levels: Transcriptional networks, specific regulons, and individual TFs. TraRe's underlying statistical model outperforms alternative models for transcriptional network inference, as shown by benchmarking scenarios with simulated data under several configurations.

We used TraRe to provide mechanistic insights into differential gene regulation associated with Abi response in mCRPC. Specifically, we applied TraRe to the transcriptomic data of 46 mCRPC samples collected before initiation of AA/P treatment from the PROstate Cancer Medically Optimized Genome Enhanced ThErapy (PROMOTE) study that was conducted in a prospective fashion (NCT: 01953640; ref. 18). First, we uncovered transcriptional modules that captured key cellular processes and whose mechanistic rewiring related to Abi response. Importantly, we were able to match these findings in a separate cohort of mCRPC patients. Furthermore, we uncovered TFs ELK3, MYB, and MXD1 that were significantly associated with response-disrupted networks. Finally, we experimentally validated in two different cell lines the key rewiring predictions identified from the response-specific networks.

TraRe: uncovering transcriptional rewiring across phenotypes from RNA sequencing data

Given a gene expression matrix for m-samples, the phenotype class of each sample, and a list of 1,453 TFs selected from the Human Protein Atlas 2020 (v19.3; HPA; RRID:SCR_006710), TraRe infers transcriptional modules (i.e., set of target genes and the driving TFs) and their phenotype-specific network disruptions (rewiring) by measuring changes in TF-target associations among phenotypes (Supplementary Fig. S1A and S1B in Supplementary File S1). To infer such rewired networks, TraRe first infers the overall transcriptional landscape of the samples and then uncovers those transcriptional modules whose topology is rewired between the phenotypes. In addition, TraRe also identifies TFs and regulons (a TF and its corresponding target genes) that are significantly prevalent in the rewired modules. Finally, to evaluate our methods, we generated a simulated dataset with underlying transcriptional modules, some of which were rewired, and its corresponding simulated gene expression data (See Supplementary Note S1 in Supplementary File S1).

Inference of transcriptional modules

First, TraRe uncovers the modules governing the transcriptomic activity of all samples via a module-based approach, where a module represents one or several related biological functions (29). The module-based approach for inferring transcriptional modules has been shown to be more accurate than inferring a unique global network and then isolating individual communities within it (29, 30).

For this, m-samples of gene expression data are initially divided into a matrix ZRt×m for the expression of t target genes and a matrix XRd×m for the expression of the designated d TFs as regulatory genes (Supplementary Fig. S1A). Z is then iteratively partitioned into K gene submodules, where each submodule's expression pattern is approximated by a sparse combination of TF expression profiles (Supplementary Fig. S1B). Specifically, two steps are taken per iteration. In the first step, for each of the K inferred submodules, the expression profile encoding the largest variance across all the submodule genes is selected as its representative. We denoted this representative “gene” expression profile as the eigengene y, as it is the first eigenvector of the sample covariance matrix of the submodule (29). Then, the set of TFs (and their associated weight) driving each submodule (represented as the sparse vector β) is inferred by approximating the expression profile of the eigengene y with a sparse linear combination of TFs. That is,
formula
This inference is done via Variational Bayes Spike Regression (VBSR; ref. 31). In previous work (29), we analyzed module-based approaches for GRN inference, in which several fitting models were evaluated. VBSR was shown to be the best approach for inferring informative networks. VBSR assumptions match GRNs underlying biological features, as target genes are usually regulated by a sparse set of TFs. In the VBSR model, a Slab and Spike prior is used to guarantee the sparsity of the generated fit, which in turn conditions the number of TFs that will regulate a specific target gene. Within the VBSR, βs are computed via a coordinate ascent algorithm that iteratively minimizes the Kullback–Leibler divergence among the approximate posterior distributions |${\hat{q}}_{{\beta }_j}$| and the complete posterior distribution of the VBSR model p(θ|W):
formula
where p(θ|W) is given by (31):
formula

The inferred approximation of the eigengene |$\hat{y}$| becomes the new submodule representative profile. In the second step of each iteration, each target gene is reassigned to the submodule whose representative is closest to the target gene as measured by the absolute value of the cosine similarity. Cosine similarity provides a relationship metric based on the cosine angle of two vectors of dimension m, hence it only evaluates the orientation patterns among genes’ expressions. It enables the match of a gene with a module representative regardless of their magnitude, as the cosine similarity will consider that a gene and a module representative are potentially related if they have similar gene expression patterns, even though expression abundance across samples is different. Finally, the absolute value of cosine similarity is taken to include those genes showing repressive patterns.

Target and TFs are iteratively reassigned to submodules as aforementioned until no further reassignments occur or when a specified number of iterations is reached.

Note that all target genes get assigned to a particular submodule. As such, there may be outliers that are not well fit by the inferred TF weights β. Thus, as a final refinement step for each inferred submodule, the expression of each target gene is modeled from the TFs in the submodule via a new VBSR model. This process yields a weighted and undirected bipartite graph for each submodule, representing its individual transcriptional network, drawing specific connections between target genes and their inferred TFs. During this step, target genes within submodules can be rejected during the fitting. Therefore, outlier target genes will be dropped out due to lack of good model fit, leading to refined submodules with consistent transcriptional networks (Supplementary Fig. S1B).

At the end of this stage, TraRe has generated a set of refined submodules depicting the transcriptional landscape of the provided input gene expression data. It is worth mentioning that TraRe's package includes two additional regression models to link targets and TFs. These models are Linear Regression Model (LM) and Least Absolute Shrinkage and Selection Operator (LASSO) model, and their assessments are described in the results section and Supplementary Notes S2–S3.

Phenotype-associated transcriptional rewiring

After inferring the refined submodules as described above, TraRe assesses them independently for evidence of network rewiring between phenotypes. Hence, the description below focuses on a single submodule. Given the gene expression matrix associated with each submodule (Supplementary Fig. S1C), the gene–gene covariance matrix associated with each binary-labeled phenotype is computed separately, namely Σ0 and Σ1. Then, the distance between the two covariance matrices is computed as the Frobenius norm of the difference between both of them:
formula

This distance s, which we term the rewiring score, is used as a statistic to measure the dissimilarity between the covariance matrices, and thus, the differential transcriptional mechanisms of the submodules across phenotypes. To assess the statistical significance of the rewiring score s, a permutation test is performed, where the rewiring score is evaluated against a null distribution H0 generated via random permutations of the phenotype labeling (Supplementary Fig. S1D). Thus, a submodule is said to be rewired if its rewiring score s satisfies |${P}_{{H}_0}( s ) < \tau $|where τ is a user-defined threshold, set by default to 0.05. Note that a significant dissimilarity over the covariance matrices may indicate fundamental underlying disruptions in gene coexpression between phenotypes.

Uncovering robust transcriptional modules across multiple runs

Unraveling transcriptional modules from RNA-seq data, also known as reverse-engineering the transcriptome, is a complex problem that generally requires heuristic algorithms. Thus, every run of TraRe (with a unique random seed) generates slightly different results. Hence, toward accounting for the implicit variability of TraRe and increasing its generalization capabilities, the module generation process is repeated B times using different subsets of 80% of the samples. Refined submodules that generalize well across different runs and different samples should be similarly captured in each run. Therefore, similar submodules across runs will cluster together, giving rise to consistent transcriptional modules in which spurious patterns coming from a specific group of samples are dropped out. After inferring the rewired submodules from all runs as mentioned above, a similarity matrix is built so that similar rewired submodules are grouped together through hierarchical clustering.

Specifically, the ij similarity value of this matrix is defined on the basis of the log significance (using the hypergeometric test) of the overlap of the pair of submodules:
formula
where N is the total number of genes in the dataset, K and n are the number of genes in the rewired submodules j and i, respectively, and k is the number of common genes between both rewired submodules. Then, hierarchical clustering is performed to group together very similar rewired submodules that have spawned across different runs, yielding robust rewired transcriptional modules (Supplementary Fig. S1E and S1F).

Uncovering rewiring-specific TFs

We classify TFs as rewiring specific if they appear statistically significantly more times in rewired submodules than in non-rewired ones. We set the re-runs of TraRe to B  =  50, and computed an enrichment P value of the overpresence of a given TF in rewired submodules using a hypergeometric test with

 TF not driving TF driving 
Rewired Rw0 Rw1 
Non-rewired |${{\overline{ {Rw}_0}}}$| |${{\overline{ {Rw}_1}}}$| 
 TF not driving TF driving 
Rewired Rw0 Rw1 
Non-rewired |${{\overline{ {Rw}_0}}}$| |${{\overline{ {Rw}_1}}}$| 

as the associated contingency table. Here, Rw and |$\overline {Rw} $| are the rewired and non-rewired submodules, respectively, and their subindexes point to whether the TF is 1 or is not (0) present in the submodule, respectively.

In addition, the odds ratio (OR) test is used to further assess the rewiring specificity of the TFs genes. Specifically, the OR is computed as:
formula
where the Haldane–Anscombe correction is applied because of the sparsity of the contingency table. Finally, rewiring-specific TFs are filtered by a threshold (set to 0.05 for the hypergeometric test and to 1 for the OR test by default) and sorted according to the associated P value.

Uncovering rewiring-specific regulons

As opposed to the previous analysis, where we evaluated the specificity of TFs within rewired submodules, target genes cannot be evaluated using the aforementioned approach, as they are assigned to a unique submodule in each run. Therefore, we chose to assess target genes via the regulons of rewired submodules that they are associated with.

From the inferred rewired submodules, we first extract all regulons, which are composed by a TF and its associated target genes (Supplementary Fig. S1G). Note that these regulons form simple networks and hence can be tested for their potential rewiring applying the same test as we did for the inferred submodules (see Supplementary Fig. S1C and S1D).

Therefore, each regulon will have an associated P value (adjusted with the Benjamini–Hochberg multiple hypothesis correction to compensate for spurious significant P values), which is informative about the network's disruption significance within that set of targets associated with a single TF. Only significant rewired regulons (adjusted P values below 0.05) are taken henceforward.

Along the 50 runs, several regulons may be led by the same TF and contain repeated associated targets; thus, we merged those regulons under the same TF (Supplementary Fig. S1G). Then, we computed the multiplicity for each target, defined as the number of times they appear within the unified regulons, and its significance, defined as the Fisher's combined probability test over the P values associated with the regulons it appears in. Targets whose P values are below a threshold (set to 0.05 by default) are dropped out and the remaining targets are sorted by the multiplicity. The filtering process can drop out every target within a regulon. If this is the case, the associated regulon is removed.

Visualization of the results in TraRe

TraRe includes a set of functions that allow the visualization of the results obtained in the stages of the method in an html report. Moreover, it outputs different summary files and figures that include the submodule similarity matrix in the form of a clustered heatmap with the corresponding dendrogram and the modeled bipartite graphs of the different robust rewired transcriptional modules (and their associated networks) that the clustering algorithm identifies.

Datasets

We have tested our method using clinical and molecular data from two mCRPC datasets: (i) Mayo Clinic clinical trial (NCT: 01953640), PROMOTE, and (ii) the Abida and colleagues (32) study from the work of the Stand Up To Cancer/Prostate Cancer Foundation-International Prostate Cancer Dream Team (SU2C).

PROMOTE study

The PROMOTE study aimed to identify the genomic alterations associated with resistance to AA/P treatment (18). The primary outcome was determined at 12 weeks of therapy using several criteria for progression such as serum prostate-specific antigen measurement, bone, and CT imaging, and symptom assessments. In the original PROMOTE study (18), whole-transcriptome sequencing (RNA-seq) was performed using the TruSeq poly-A selection library and Illumina HiSeq 2000 instruments. For this article, we obtained the RNA-seq from 64 tumor samples and normalized the gene counts using conditional quantile normalization. These gene expression data came from 46 bone and 18 soft tissue metastatic sites before the AA/P treatment began. We performed our TraRe analysis on the total set of genes (21,302 genes, of which, 1,453 are TFs and the rest are considered as target genes) and only the 46 bone samples, from which 29 and 17 were Abi responders (R) and nonresponders (NR), respectively.

Stand Up To Cancer study

In addition to the PROMOTE study, we ran TraRe on the “Stand Up To Cancer” (SU2C) study data, which also consists of patients with mCRPC disease. We downloaded the clinical and sequencing data from the database of Genotypes and Phenotypes (dbGaP) with accession code phs000915.v2.p2. Because the PROMOTE RNA-seq study data were generated using TruSeq poly-A selection library, in the SU2C study, we filtered out the data generated using a hybridization capture-based library. We obtained the RNA-seq gene expression data (FPKM) from 270 SU2C samples generated using the TruSeq poly-A selection library, and removed genes with no expression in more than 90% of the samples. The FPKM data were then quantile normalized and log2 transformed.

Study harmonization

To harmonize the information of both studies, only those patients treated with Abi in the SU2C study were selected. Specifically, 121 patients whose transcriptomic data were available before Abi treatment was used to infer the transcriptional modules.

Note that the number of genes examined from the two matrices differs. Therefore, only common genes (12,791) have been selected when evaluating both datasets jointly. In all other analyses, the original features (genes) have been maintained. To partition the genes in this analysis between TFs and targets, we defined the driver genes as all genes that match the list of TFs selected from the HPA 2020 (v19.3; RRID:SCR_006710; a total of 1453), and the targets as all of the remaining genes.

Enrichment analysis

We used the CommunityAMARETTO algorithm (33) to cluster together similar inferred submodules across different runs of TraRe, and assign them a biological meaning. Using a hypergeometric test with an FDR < 1e−6, a network of overlapping submodules is created. Then, densely connected submodules (i.e., network communities) are identified using the Girvan–Newman algorithm (34) and clustered into transcriptional modules, which are assessed for enrichment (hypergeometric test) with known functional categories in the curated (C2) and Hallmark (H) gene sets of Molecular Signatures Database (MSigDB) v5.2 (35).

For the Gene Ontology (GO) Biological Processes enrichment analysis of the genes included in the rewiring-specific TF list, regulons and rewired transcriptional modules, the following R packages were used: clusterProfiler (v4.0.5; RRID:SCR_016884), org.Hs.eg.db (v3.13.0), AnnotationDbi (v1.54.1), DOSE (v3.18.3), rrvgo (v1.4.4), and GOSemSim (v2.18.1).

Chromatin immunoprecipitation sequencing enrichment

We downloaded from the Remap 2022 database (36), the latest version of nonredundant peaks from Homo sapiens data (see https://remap.univ-amu.fr/storage/remap2022/hg38/MACS2/remap2022_nr_macs2_hg38_v1_0.bed.gz) and then used the R package ChIPpeakAnno (v.3.30.1; RRID:SCR_012828) to map peaks from chromatin immunoprecipitation sequencing (ChIP-seq) experiments with our inferred and most robust rewired TF-target edges. Criteria to consider that a peak corresponds to the selected target are: (i) if the region is upstream of the gene start, no further than 5kb, (ii) if the region maps into some internal features like introns but is not too distant to the nearest feature, a maximum of 5 kb, (iii) if the region overlaps with the start of the described gene or feature, and (iv) if the region is contained within an exon or described feature of the gene. We then evaluated, for every TF, its TF-target ChIP-seq enrichment via a hypergeometric test on the following contingency table:

 Target in ChIP Target NOT in ChIP 
Target in TraRe True Positive False Positive 
Target NOT in TraRe False Negative True Negative 
 Target in ChIP Target NOT in ChIP 
Target in TraRe True Positive False Positive 
Target NOT in TraRe False Negative True Negative 

Note that for this test we used the total set of PROMOTE genes (21,302) as the universe of genes, thus restricting ChIP-seq peaks to only those.

GRN method comparison

We benchmarked the TraRe's GRN inference method against two known GRN inference algorithms: GRNboost2 (37) and ARACNE-AP (38) using the same input data as we used with TraRe (i.e., the total set of genes 21,302 for the 46 bone samples and the same list of TFs as described above). The assessment was done in two ways: With all the samples to mimic TraRe's inference method, and separately for responder and nonresponder to assess phenotype-network differences. When using all samples to generate the GRNs, the network was evaluated by comparing the enrichment of ChIP-seq peaks. Specifically, we checked whether the edges inferred by these two methods were associated with the TFs leading our rewired regulons (i.e., those inferred with TraRe) and then proceeded to calculate the enrichment in ChIP-seq peaks as described above. Because GRNboost2 generates larger regulons and delivers an "importance" metric, we ranked the targets by their importance and analyzed the top 25%. In the case of ARACNE-AP, because the networks were much smaller we tested the entire regulons.

We then used the networks generated independently for each phenotype to test whether these GRN inference methods were able to infer the relationships that we experimentally validated (see next Materials and Methods sections). For this, we extracted the inferred edges that were associated with the selected TFs for validation (ELK3, MYB, MXD1, ZNF3, and ZNF91). We then generated phenotype-driven distributions based on the edges’ weight metric (“importance” for GRNboost2 and “Mutual information” for ARACNE-AP). Finally, we assessed the statistical significance of the targets that we experimentally validated using a P value of 0.05 on those distributions.

Cell culture

22Rv1 and LNCaP cells purchased from the ATCC were cultured in RPMI-1640 medium (Gibco) supplemented with 10% FBS (Sigma) and 1% Pen–Strep. To develop Abi-resistant cell lines (LNCaP-AR or 22Rv1-AR), cells were maintained in the medium supplemented with 5 μmol/L of abiraterone (Selleck Chemicals) until viability reached over 95%. Abi resistance was validated by proliferation assay. For knockdown experiments, cells were seeded in regular RPMI-1640 in 6-well plates, and were transfected using siRNA against ELK3, MXD1, MYB, ZNF3, ZNF91 or nontargeting siRNA (Horizon Discovery) by RNAiMax (Thermo Fisher Scientific) according to the manufacturer's protocol. siRNA sequences can be found in Supplementary Table S1 in Supplementary File S1.

Proliferation assay

36 hours before Abi treatment, cells were seeded in phenol red-free RPMI1640 (Gibco) supplemented with 10% charcoal stripped FBS (Sigma) in 96-well plates, and 50 nmol/L pregnenolone was added after 16 hours. Cells were then treated with Abi or vehicle and monitored for proliferation at days 0, 2, 4, and 6 after initiation of treatment by Cyquant direct assay (Thermo Fisher Scientific). For proliferation after knockdown, cells were trypsinized 24 hours after transfection and reseeded in 96-well plates. Proliferation was monitored at days 1, 3, 5, and 7 after reseeding.

qRT-PCR

Total RNA for qRT-PCR was extracted from cell lines 48 hours after knockdown using the Quick-RNA MiniPrep Kit (Zymo Research) according to the manufacturer's instructions. qRT-PCR was performed using the Power SYBR Green RNA-to-CTTM 1-Step Kit (Life Technologies) and QuantiTect (QIAGEN) or Prime-Time (IDT, Inc.) predesigned qPCR primers (IDT). Gene expression analyses were performed using the ∆∆Ct method, and GAPDH was used as the internal reference. Two independent experiments were performed. Primer sequences are in Supplementary Table S2.

Wound-healing assay

LNCaP, 22Rv1 or their Abi-resistant cells were seeded in 6-well plates at 6×105 density, and transfected with siRNA-targeting ELK3 or scrambled siRNA control. 24 hours after transfection, a scratch wound was made by pipet tip. Images were then taken immediately and after 24 hours.

Data availability

The datasets analyzed in this study are available in the dbGaP repository with accession numbers phs001141.v1.p1 for PROMOTE and phs000915.v2.p2 for SU2C, upon request. SU2C is also available in the cBioPortal for Cancer Genomics (RRID:SCR_014555) database under the name of Metastatic Prostate Adenocarcinoma (SU2C/PCF Dream Team, PNAS 2019) study.

The TraRe package (28) is available at Bioconductor (RRID: SCR_006442) and the GitHub repository ubioinformat/TraRe. Other data generated in this study and source data for figures, tables, and Supplementary Materials as well as scripts used for the analysis are available at the GitHub repository ubioinformat/TraRe-materials. Also, an online-executable capsule is made available at Code Ocean (https://codeocean.com/capsule/4965686/tree/v1) and the link is also provided in the GitHub repository. Other data generated in this study are available upon request from the corresponding author.

Study and method overview

We applied the developed TraRe framework to the gene expression profiles of 46 mCRPC baseline (pretreatment) samples from the PROMOTE study. Abi response annotation for all patients was measured at 3 months after treatment initiation (Fig. 1: Input Data, see Materials and Methods) and patients were separated into two groups: treatment responders and nonresponders (18). Given the availability of pre-Abi treatment gene expression of the PROMOTE samples and a set of known TFs selected from the HPA 2020 (v19.3; RRID:SCR_006710), TraRe first performed an iterative clustering process, which partitions target genes into different submodules, whose expression can be modeled by the expression of a sparse set of TFs (Supplementary Fig. S1A and S1B, see Materials and Methods). To overcome sensitivity to outliers in the data, TraRe was run several times, each run randomly selecting a different subset of the patient samples (80–20 split), as well as a random initialization of the parameters. We then created our transcriptional modules by clustering the similar inferred submodules across the multiple runs (Fig. 1: TraRe). We reasoned that transcriptional modules sharing common sets of genes across sub-sampled runs would uncover a functional landscape of biological processes that are not a consequence of overfitting to a specific set of patients, and therefore, better generalize to the disease population. To validate this hypothesis, transcriptional modules were independently inferred on the SU2C (32) dataset consisting of 121 patients with mCRPC (Fig. 1: Input Data, see Materials and Methods). Module-oriented community detection and functional annotation methods (33) showed that indeed our method found shared modules between PROMOTE and SU2C (Fig. 1: Validation).

Figure 1.

TraRe's workflow overview. Input Data: Datasets used in the study include bulk RNA-seq data from SU2C and PROMOTE studies of mCRPC cohorts before abiraterone treatment, as well as simulated data. Simulated gene modules were created from linear combinations of sampled PROMOTE regulatory TFs with simulated target gene expression drawn from Multivariate Gaussian distributions and visualized with heatmaps and principal component analysis (PCA). TraRe: The tool presented here follows two steps: Step 1 that builds transcriptional submodules and that is used for functional annotation and Step 2 that identifies rewired submodules, ranks their TFs and regulons for association with the response, and examines their phenotype-specific transcriptional modules. Validation: In this study, we validated our findings with comparative analysis of transcriptional modules from the PROMOTE and SUC2 studies, in vitro TF knockdown assays with naïve and Abi-resistant cell lines to assess proliferation and target mRNA expression, and by evaluating model performance in recovering seeded transcriptional modules from simulated data.

Figure 1.

TraRe's workflow overview. Input Data: Datasets used in the study include bulk RNA-seq data from SU2C and PROMOTE studies of mCRPC cohorts before abiraterone treatment, as well as simulated data. Simulated gene modules were created from linear combinations of sampled PROMOTE regulatory TFs with simulated target gene expression drawn from Multivariate Gaussian distributions and visualized with heatmaps and principal component analysis (PCA). TraRe: The tool presented here follows two steps: Step 1 that builds transcriptional submodules and that is used for functional annotation and Step 2 that identifies rewired submodules, ranks their TFs and regulons for association with the response, and examines their phenotype-specific transcriptional modules. Validation: In this study, we validated our findings with comparative analysis of transcriptional modules from the PROMOTE and SUC2 studies, in vitro TF knockdown assays with naïve and Abi-resistant cell lines to assess proliferation and target mRNA expression, and by evaluating model performance in recovering seeded transcriptional modules from simulated data.

Close modal

We next aimed at unraveling those transcriptional networks from PROMOTE showing differential regulation between R and NR. Specifically, for each submodule generated across runs, we performed a statistical rewiring test (see Materials and Methods and Supplementary Fig. S1C and S1D) to identify those whose TF-target association patterns were significantly different between the two response groups (Fig. 1: TraRe). We considered a transcriptional module to be rewired if at least 40% of its submodules showed significant rewiring by our test.

We then sought more insights into the key TFs potentially driving the response-specific rewiring process. We used Bayesian sparse regression (10) models in each submodule to associate each TF to a specific set of target genes. This set of target genes is known as the TF's regulon (Fig. 1: TraRe, Supplementary Fig. S1G, see Materials and Methods). We, thus, were able to rank TFs and their regulon by their associations with the response-specific rewiring process. Finally, we closely investigated the regulons of key rewired transcriptional networks to uncover specific relationships between TFs and targets that behaved distinctly between response groups.

Therefore, in this study we developed TraRe, a computational method that provides a three-tier analysis to identify key biological processes and genes associated with phenotypically driven mechanistic differences: (i) at the network level, by inferring rewired transcriptional networks; (ii) at the regulon level, by identifying specific TF-target relationships that may be associated with phenotypic differences; and (iii) at the single gene level by identifying TFs consistently associated with rewired networks. We concluded this study by experimentally validating a subset of our findings with TF knockdown in Abi-resistant prostate cancer cell lines. We also shared our findings on the ideal selection of regression models for uncovering phenotype-associated transcriptional networks, as well as the success of the proposed rewiring statistical tests on varied simulated benchmark datasets (Fig. 1: Validation).

TraRe recapitulates transcriptional programs associated with mCRPC

We first focused on uncovering the transcriptional landscape and associated key biological processes of mCRPC. To that end, we ran TraRe on the 46 PROMOTE samples from bone biopsies taken previous to the initiation of AA/P treatment. TraRe was run 10 times, each one sampling without replacement 80% of the samples and setting to infer up to 100 submodules (K = 100). This number of submodules was chosen following previous works (29, 30, 39, 40). Nevertheless, we tested different values for the K parameter and lower numbers yielded submodules with impractical sizes (K = 50), whereas larger numbers [K = (200, 300, 500)] generated submodules of sizes unfit for gene set enrichment analysis (Supplementary Fig. S2). Across all runs, TraRe inferred a total of 835 submodules and identified 83 distinct groups of overlapping submodules (hypergeometric test, FDR < 10−5; see Materials and Methods). Of the 83, 38 groups contained submodules shared among at least half of the 10 runs, and these were identified as robust transcriptional modules (See Supplementary Tables S3–S5). These 38 transcriptional modules were functionally annotated with CommunityAMARETTO software (33) using curated C2 and hallmark gene set enrichment analysis (Fig. 2A; see Materials and Methods; Supplementary File S2). Some of them, such as P-M15 and P-M18, contained genes enriched in gene sets related to cell constitutive and housekeeping functions such as G protein–coupled receptors (GPCR) and peptide-ligand binding receptors (P-M15 and P-M18) or cell-cycle and transcription (P-M15), hence they are not explicitly associated to Abi response (Fig. 2B).

Figure 2.

TraRe recapitulates complex regulatory interactions on metastatic castration-resistant prostate cancer. A, Grouping of submodules by CommunityAMARETTO software after running TraRe on PROMOTE dataset with 9 transcriptional modules and their associated functional categories highlighted. B, Gene Ontology gene sets of biological process terms (rows) in which the rewired transcriptional modules (columns) are enriched and colored by functional category and significance. C and D, Categorization of functional annotation enrichments found with CommunityAMARETTO software in selected PROMOTE (C) and SU2C (D) transcriptional modules. E, Similarity heatmap (Jaccard index) between rewired transcriptional modules in PROMOTE study and overlapping SU2C study transcriptional modules.

Figure 2.

TraRe recapitulates complex regulatory interactions on metastatic castration-resistant prostate cancer. A, Grouping of submodules by CommunityAMARETTO software after running TraRe on PROMOTE dataset with 9 transcriptional modules and their associated functional categories highlighted. B, Gene Ontology gene sets of biological process terms (rows) in which the rewired transcriptional modules (columns) are enriched and colored by functional category and significance. C and D, Categorization of functional annotation enrichments found with CommunityAMARETTO software in selected PROMOTE (C) and SU2C (D) transcriptional modules. E, Similarity heatmap (Jaccard index) between rewired transcriptional modules in PROMOTE study and overlapping SU2C study transcriptional modules.

Close modal

We then sought to test our transcriptional modules for differential regulation between 29 Abi responders and 17 nonresponders in the PROMOTE dataset. To achieve this goal, we applied our rewiring statistical test (see Materials and Methods) to all submodules separately. We found that 4 out of the 38 transcriptional modules (P-M1 to P-M4) were enriched with rewired submodules. From the previous enrichment analysis (Fig. 2B), these rewired transcriptional modules were associated with embryological and hematopoietic-related gene sets (P-M1), mitochondrial gene sets (P-M3), and different gene sets related to cancer and reproductive system organs (P-M2 and P-M4). We additionally examined these rewired transcriptional modules for their enrichment in GO Biological Processes and found terms related to the functions mentioned above: Mitochondrial electron transport in P-M3 and activation of different hematologic cells and immune response in P-M1 (neutrophil, leukocytes, etc.). In the case of the other two transcriptional modules, the most important GO-associated functions were related to RNA processing (P-M2) and cell migration, angiogenesis, and extracellular matrix organization (P-M4; Fig. 2C; Supplementary Table S6 in Supplementary File S3).

Other transcriptional modules containing rewired submodules (P-M5, P-M8, and P-M9) were found to be enriched with genes downregulated in prostate cancer samples (as well as different cancer types; Fig. 2A and B). ETV5, a TF from P-M9, is involved in the AR signaling pathway associated with invasion and it has been identified to participate in rare gene fusion permutations in prostate cancer (41). Interestingly, P-M5, with two rewired submodules, was found to be enriched with liver and xenobiotic metabolism-related pathways such as HNF1A, cytochrome P450 and glucuronidation, affecting drug response and efficacy.

To further validate the transcriptional modules uncovered by TraRe in the PROMOTE dataset, we ran TraRe on the 121 mCRPC samples of the SU2C dataset (see Materials and Methods) with the same run settings. We found 49 robust transcriptional modules in the SU2C dataset. Six of these SU2C modules were highly overlapping with modules identified from the PROMOTE dataset. These six modules were enriched with cell-cycle–related genes that were also found to be differentially regulated in different cancer cells (SU2C S-M5), hematologic (S-M4), embryological and cell cycle (S-M1), breast and prostate cancer gene sets (S-M2, S-M3 and S-M6; Fig. 2D). We were especially interested in the significant overlaps of the six SU2C transcriptional modules with our rewired transcriptional modules from PROMOTE (Fig. 2E). P-M1 from PROMOTE was closely related with S-M4 from SU2C, as expected by their common gene set enrichment analysis (hematologic gene sets, Fig. 2D), and also with S-M2, a module enriched with immune gene sets. Similarly, P-M4 (related with breast and prostate cancer) from PROMOTE was related to the modules S-M2, S-M3, and S-M6 from SU2C, primarily enriched with the same cancer gene sets (Fig. 2D). As a final note, we were not able to compare whether the SU2C modules have similar associations to PROMOTE with response-specific rewiring due to the lack of an analogous clinical response measurement after three months of treatment.

Thus, TraRe was able to generate biologically meaningful transcriptional modules from PROMOTE data that were consistent across multiple runs (Fig. 2A). Specifically, the highlighted PROMOTE transcriptional modules P-M1, P-M2, P-M3, and P-M4, which were enriched with differentially rewired submodules between R and NR, captured novel transcriptional associations that may help elucidate the foundations underpinning abiraterone sensitivity and resistance. Similar modules were also observed in an independent mCRPC study (SU2C).

TraRe identifies TFs playing a key role in abiraterone response

Next, we investigated individual TFs that may be driving significant differences in the transcriptional modules between R and NR. We hypothesized that the rewired transcriptional modules inferred in the previous section are likely to be driven by TFs that themselves might not be rewired between the two groups as the conserved behaviors of these TFs could underpin the module that they drive (as we show in the following section). To increase the statistical power required when working at single-gene resolution, we ran TraRe 50 times on the 46 PROMOTE samples with the default settings (see Materials and Methods), yielding 4,088 submodules, from which, 225 were marked as significantly rewired by TraRe. We then selected those TFs that were identified as driving the core submodule expression and that were highly enriched in rewired submodules (see Materials and Methods). The resulting ranking contained 33 TFs that likely played a key role in the rewiring of the transcriptional modules associated with the different Abi treatment outcomes (Table 1).

Table 1.

TFs specific to rewired modules.

Transcription factorPOdds ratio
CEBPE 1.4e−31 205.94 
GATA1 1.6e−29 82.53 
GLI1 4.6e−14 90.74 
KLF1 2.1e−12 125.89 
PRRX1 2.0e−11 30.21 
NFE2 7.4e−10 48.56 
MXD1 2.8e−09 25.14 
ZNF3 6.1e−09 30.89 
MYB 1.9e−06 32.77 
ZNF91 1.0e−05 13.97 
IRF4 1.8e−05 17.63 
PAX6 2.8e−04 7.16 
NR1H4 4.7e−04 5.47 
ELK3 7.7e−04 4.01 
ZNF181 8.9e−04 14.30 
KLF10 1.4e−03 12.10 
ZFP3 1.5e−03 24.38 
ZFHX2 1.5e−03 24.38 
TAL1 1.5e−03 24.38 
ZSCAN22 2.9e−03 6.64 
EPAS1 3.7e−03 6.21 
ZNF174 4.8e−03 13.53 
TFAP4 4.8e−03 13.53 
RUNX1 0.01 4.14 
BATF3 0.01 28.90 
ZNF526 0.01 28.90 
SOX1 0.02 17.33 
ZNF574 0.03 12.38 
ID2 0.03 5.79 
ZNF286A 0.04 9.62 
RBAK 0.04 9.62 
ZNF519 0.04 3.82 
SNAI2 0.05 2.67 
Transcription factorPOdds ratio
CEBPE 1.4e−31 205.94 
GATA1 1.6e−29 82.53 
GLI1 4.6e−14 90.74 
KLF1 2.1e−12 125.89 
PRRX1 2.0e−11 30.21 
NFE2 7.4e−10 48.56 
MXD1 2.8e−09 25.14 
ZNF3 6.1e−09 30.89 
MYB 1.9e−06 32.77 
ZNF91 1.0e−05 13.97 
IRF4 1.8e−05 17.63 
PAX6 2.8e−04 7.16 
NR1H4 4.7e−04 5.47 
ELK3 7.7e−04 4.01 
ZNF181 8.9e−04 14.30 
KLF10 1.4e−03 12.10 
ZFP3 1.5e−03 24.38 
ZFHX2 1.5e−03 24.38 
TAL1 1.5e−03 24.38 
ZSCAN22 2.9e−03 6.64 
EPAS1 3.7e−03 6.21 
ZNF174 4.8e−03 13.53 
TFAP4 4.8e−03 13.53 
RUNX1 0.01 4.14 
BATF3 0.01 28.90 
ZNF526 0.01 28.90 
SOX1 0.02 17.33 
ZNF574 0.03 12.38 
ID2 0.03 5.79 
ZNF286A 0.04 9.62 
RBAK 0.04 9.62 
ZNF519 0.04 3.82 
SNAI2 0.05 2.67 

Interestingly, the highest-ranked TFs were those also found in the most significantly rewired module (P-M1), such as CEBPE, GATA1, or KLF1. Furthermore, using the STRING database (42), we found that the interaction networks between the top-ranked TFs (Fig. 3A) showed multiple connections with enriched connectivity (STRING PPI enrichment P value 1.2e−15) especially between those TFs identified in P-M1. The top enriched pathways associated with those TFs included regulation of hemopoiesis and different blood cell differentiation pathways, specifically myeloid cells, indicating a strong functional relation among these TFs (Fig. 3B; Supplementary Table S7 in Supplementary File S3). These results showed the robustness of the rewiring analysis of the P-M1 because its TFs were also among the top significant rewired TFs, highlighting the capability of TraRe to infer transcriptional modules with TFs and target genes that are associated with the same pathways and biological functions.

Figure 3.

TraRe unraveled transcription factors playing a key role in abiraterone response of mCRPC patients. A, Protein–protein interaction network of TFs specific to rewired submodules. Relationships are based on experimental data (purple), curated databases (blue), coexpression (black), and text mining (yellow). Seventeen out of 33 nodes contain an edge. Figure generated with STRING Database2 web tool v.11.0 (https://version-11-5.string-db.org/cgi/network?networkId=bxdLVX8dezb1. Permanent Link.). B, Enrichment analysis of top Gene Ontology biological process terms with the rewiring-specific TFs using R package clusterProfiler (v4.0.5) and org.Hs.eg.db (v3.13.0). C, Graphical representation of three example regulons (low multiplicity edges filtered for visualization purposes). D, Enrichment analysis of Gene Ontology biological process broad terms (rows) with the target gene sets for each TF-driven regulon (columns).

Figure 3.

TraRe unraveled transcription factors playing a key role in abiraterone response of mCRPC patients. A, Protein–protein interaction network of TFs specific to rewired submodules. Relationships are based on experimental data (purple), curated databases (blue), coexpression (black), and text mining (yellow). Seventeen out of 33 nodes contain an edge. Figure generated with STRING Database2 web tool v.11.0 (https://version-11-5.string-db.org/cgi/network?networkId=bxdLVX8dezb1. Permanent Link.). B, Enrichment analysis of top Gene Ontology biological process terms with the rewiring-specific TFs using R package clusterProfiler (v4.0.5) and org.Hs.eg.db (v3.13.0). C, Graphical representation of three example regulons (low multiplicity edges filtered for visualization purposes). D, Enrichment analysis of Gene Ontology biological process broad terms (rows) with the target gene sets for each TF-driven regulon (columns).

Close modal

To provide a better understanding of the transcriptional modules that were generated by TraRe, we explored the relationships between TFs and target genes at the regulon level. A regulon is defined as a set of target genes associated with one TF (Supplementary Fig. S1G). Therefore, we identified robust rewired regulons across the rewired submodules (see Materials and Methods, Supplementary Table S8 in Supplementary File S3). A graphic example of the regulons corresponding to the three TFs selected for further biological validation (ELK3, MYB, and MXD1) is shown in Fig. 3C. In total, there were 50 significant rewired regulons, and almost all of the top TFs described previously, such as MYB or MXD1, were masters of these regulons. Functional enrichment analysis with GO terms of the regulons shows that the most frequent terms outputted from the over-representation analysis with biological processes were under broader terms related to regulation of immune response and cell migration, extracellular matrix organization, cell adhesion, and angiogenesis, as well as different receptor signaling pathways (Fig. 3D; Supplementary Table S9 in Supplementary File S3). Consistent with the results obtained in the functional annotation of the rewired transcriptional modules, the regulons of the master TFs presented in P-M1, such as CEBPE, GATA or NFE2, were enriched in gene sets related to neutrophils and innate immune response.

Although our framework does not infer direct TF-target relationships as no TF-binding motifs or additional data other than gene expression are used, we assessed the enrichment of the inferred rewired regulons with known DNA-binding sequencing datasets (ChIP-seq) from the Remap 2022 database (36). Of the 32 regulons with ChIP-seq data available, 23 presented 5 or more targets with associated peaks, of which, 56.5% were significantly enriched in TraRe's inferred edges (see Materials and Methods). These numbers were much lower for the other two GRN inference methods evaluated (see last section of the Results for a description of the used methods, and Supplementary File S4—EXCEL, Supplementary Fig. S3 and Supplementary Note S4 for additional details).

In summary, we showed that TraRe was able to uncover key TFs and regulons related to Abi response and associate them with functional pathways. Furthermore, these analyses at both single-gene (TF) level and regulon level seamlessly complemented the findings at module level.

TraRe uncovers mechanistic changes between response-specific transcriptional networks

To better understand how the interactions between TFs and their target genes differ between the samples from Abi R and NR, we focused on the significantly rewired and highly robust transcriptional module P-M1 (see Supplementary Table S3). P-M1 was composed of 11 submodules from the 10 original runs of TraRe. Note that each run would ideally contribute at least one submodule to the transcriptional module P-M1, indicating that P-M1 is consistently discovered across different runs. Indeed, most runs contributed with one submodule to P-M1, with the first and seventh runs splitting the associated genes into two separate submodules (Supplementary Fig. S4). All 11 of the individual submodules were statistically significantly rewired with respect to Abi response, with a total of 9 driving TFs found (Supplementary Fig. S5); GATA1, CEBPE, FOXN1, MYB, TAL1, GLI1, KLF1, MXD1, and NFE2. We merged the 619 unique genes of the 11 individual submodules and used the VBSR model (see Materials and Methods) to redraw the combined transcriptional network using all 46 bone samples (Fig. 4A, right). The resulting network successfully modeled the expression of 222 of the target genes using the 9 different TFs. Similar results were observed using all patient samples and alternative regression methods, see Supplementary Note S5 and Supplementary Table S10 in Supplementary File S1.

Figure 4.

Regulatory signals of P-M1. A, Sparse VBSR regulatory networks for P-M1 (hematopoietic and immune response) built from expression data from only responders, only nonresponders, and all bone samples. Transcription factors (orange squares) are connected to their target genes (blue circles) by edges colored by their regulatory relationships (green, positive; red, negative). B, For each regulator of the transcriptional module, a violin plot showing the regulator's expression among responders (pink) and nonresponders (blue). C, Changes in the expression (t-statistic) of the global all-bone sample targets of each regulator between responders and nonresponders. D, Correlation relationships between regulators and their all bone sample targets among samples of the different response classes. E, Pairwise correlation among P-M1 genes within responders and nonresponders.

Figure 4.

Regulatory signals of P-M1. A, Sparse VBSR regulatory networks for P-M1 (hematopoietic and immune response) built from expression data from only responders, only nonresponders, and all bone samples. Transcription factors (orange squares) are connected to their target genes (blue circles) by edges colored by their regulatory relationships (green, positive; red, negative). B, For each regulator of the transcriptional module, a violin plot showing the regulator's expression among responders (pink) and nonresponders (blue). C, Changes in the expression (t-statistic) of the global all-bone sample targets of each regulator between responders and nonresponders. D, Correlation relationships between regulators and their all bone sample targets among samples of the different response classes. E, Pairwise correlation among P-M1 genes within responders and nonresponders.

Close modal

Then, we created VBSR-based transcriptional networks for the 29 R and 17 NR patients separately (Fig. 4A left and center) to examine the differences in the underlying transcriptional networks for the two groups. Interestingly, the expression levels of the nine driving TFs were, generally, higher and more variable in R than NR (Fig. 4B). Our method was also able to identify important module TFs that themselves were not differentially expressed (GATA1, KLF1, and MYB) between the two patient groups. TAL1, a BHLH DNA-binding TF essential for maintaining the multipotency and quiescence of hemopoietic stem cells (43) and implicated in the development of hemopoietic malignancies, was significantly differentially expressed (P value 0.036) with decreased expression levels in NR. When considering the target genes of the module, we found TFs such as NFE2 and CEBPE, each associated with more than 20% of the 222 targets in the bone samples, whereas others like MYB and FOXN1 were only associated with a few targets. We found that across all TFs, the targets of our module driving TFs were generally more highly expressed in R (Fig. 4C), with NFE2 and CEBPE associated with the targets that were most strongly differentially expressed.

When we turned our attention to the types of relationships between TFs and their targets in P-M1, we observed that there was generally a positive association representing an activatory role of the TFs. In the patients who respond to abiraterone treatment, the positive association between the TFs and targets was much stronger (Fig. 4D). For the 25 targets of TAL1 and the 33 targets of KLF1, the strong positive associations were mostly preserved in R and NR, whereas for other TFs like MXD1 (32 targets) and GLI1 (15 targets), their activation of the targets in NR seemed to be largely lost. The loss of the module cohesion in NR became even more pronounced in the context of associations between all genes in P-M1 (Fig. 4E). A large percentage of TF-target pairs in this transcriptional module that was tightly associated in the responders showed weak or negative associations in the subset of patients that did not show improvement at 3 months (NR).

In summary, TraRe identified important TFs of P-M1 that played distinct roles in phenotype-specific network rewiring, with target genes in P-M1 appearing to be much more strongly regulated in patients that responded to abiraterone treatment and with looser transcriptional module cohesion found in the nonresponders.

The role of ELK3, MYB, and MXD1 in abiraterone response

We experimentally validated using prostate cancer cell lines, our findings that identified regulons with significant transcriptional network disruptions between Abi R and NR. To examine the potential differential-regulation between Abi-sensitive and Abi-resistant settings, we developed LNCaP and 22Rv1 Abi-resistant cell lines, namely LNCaP-AR and 22Rv1-AR (Fig. 5A and B). We started by selecting five regulons (ELK3, MXD1, MYB, ZNF3, and ZNF91) based on the following criteria: (i) regulon significantly rewired between responders and nonresponders (Fisher method combined P < 0.01); (ii) the TF genes were expressed in both LNCaP and 22Rv1 using Cancer Dependency Map Portal (RRID:SCR_017655) database (logTPM > 1); and (iii) TF associated with a series of highly repeatable targets for at least one of the response groups (Fig. 5C; Supplementary Table S11 in Supplementary File S1) that were also expressed on the basis of DepMap database (logTPM > 1). All TFs except ZNF91 were highly coexpressed with their targets in R only (Fig. 5C; Supplementary Fig. S6). We performed siRNA knockdown of each TF (Fig. 5D), and examined the top-ranked predicted targets based on multiplicity. We found that targets in three regulons, namely ELK3,

Figure 5.

Experimental validation of regulons is significantly different between responders and nonresponders. A and B, Proliferation assay of 22Rv1 (A) and LNCaP (B). Parental and Abi-resistant cell lines in the presence of vehicle or abiraterone presented as mean ± SD of three replicates. C, Multiplicity of target genes identified within five selected regulons. D, qRT-PCR measurement of knockdown efficiency of selected transcription factors, as mean ± SD of two replicates. E and F, Quantification of expression of targets after knockdown of selected transcription factors using 22Rv1 (E) and LNCaP (F) and their respective Abi-resistant cell lines by qRT-PCR, as mean ± SD of two replicates. G–L, Proliferation assay of 22Rv1 and 22Rv1 Abi-resistant cell lines after knockdown of ELK3 (G), MXD1 (I), and MYB (K), presented as mean ± SD of three replicates. Proliferation assay of LNCaP and LNCaP Abi-resistant cell lines after knockdown of ELK3 (H), MXD1 (J) and MYB (L). M, KEGG pathway analysis of ELK3-targeted genes. N, Quantification of wound-healing assay. Data are presented as mean ± SD of six replicates. Statistical difference was tested by a paired two-sided Student t test. O, Representative pictures of wound-healing assay of 22Rv1, LNCaP, and their Abi-resistant cell lines after knockdown of ELK3. Pictures were taken at 0 and 24 hours after scratch. n.s., nonsignificant.

Figure 5.

Experimental validation of regulons is significantly different between responders and nonresponders. A and B, Proliferation assay of 22Rv1 (A) and LNCaP (B). Parental and Abi-resistant cell lines in the presence of vehicle or abiraterone presented as mean ± SD of three replicates. C, Multiplicity of target genes identified within five selected regulons. D, qRT-PCR measurement of knockdown efficiency of selected transcription factors, as mean ± SD of two replicates. E and F, Quantification of expression of targets after knockdown of selected transcription factors using 22Rv1 (E) and LNCaP (F) and their respective Abi-resistant cell lines by qRT-PCR, as mean ± SD of two replicates. G–L, Proliferation assay of 22Rv1 and 22Rv1 Abi-resistant cell lines after knockdown of ELK3 (G), MXD1 (I), and MYB (K), presented as mean ± SD of three replicates. Proliferation assay of LNCaP and LNCaP Abi-resistant cell lines after knockdown of ELK3 (H), MXD1 (J) and MYB (L). M, KEGG pathway analysis of ELK3-targeted genes. N, Quantification of wound-healing assay. Data are presented as mean ± SD of six replicates. Statistical difference was tested by a paired two-sided Student t test. O, Representative pictures of wound-healing assay of 22Rv1, LNCaP, and their Abi-resistant cell lines after knockdown of ELK3. Pictures were taken at 0 and 24 hours after scratch. n.s., nonsignificant.

Close modal

MXD1 and MYB, exhibited significantly differential association with downstream targets between Abi-sensitive parental cells and Abi-resistant cells (Fig. 5E and F). Specifically, among the 30 tested targets in these three regulons, 19 were significantly decreased in at least one parental cell line but not in the Abi-resistant line, 10 were significantly decreased in both parental cell lines but not in either Abi-resistant lines after siRNA knockdown. Therefore, we further tested whether these three regulons might differentially affect the proliferation in parental versus Abi-resistant prostate cancer cell lines. Intriguingly, we found that knockdown of ELK3, MXD1, and MYB appeared to negatively affect the proliferation in 22Rv1 and LNCaP parental cell lines but had minimal or no impact on the Abi-resistant cells (Fig. 5GL). These results indicate that the parental cell lines depend more on these TFs to survive.

We also performed pathway analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to examine downstream functionalities affected by these TFs. We found that the targets of ELK3 were enriched in focal adhesion and the Rap1 signaling pathway, which suggested a potential role of ELK3 in cell migration (Fig. 5M). We thus performed wound-healing assay in ELK3 or nontargeting control siRNA-transfected cells. We found that the migration was significantly inhibited in 22Rv1 and LNCaP parental cell lines, but not in the corresponding Abi-resistant cell lines when ELK3 was suppressed (Fig. 5N and O).

Overall, the in silico predictions of TraRe on associations of key TFs on inferred targets were overwhelmingly supported in two models of Abi-naïve and -resistant cell lines when examined for differential response and cell proliferation.

TraRe performance on simulated data benchmarks and comparison with different GRN algorithms

To assess the performance of TraRe, we generated a synthetic expression matrix from 10 simulated transcriptional modules. These were created from randomly sampled sets of TFs driving PROMOTE data modules, and adding multivariate Gaussian noise to generate expression patterns of target genes (see Supplementary Note S1, Fig. 1: Simulated Data)

Three well-known regression models were evaluated for their ability to discover the underlying simulated modules: (i) VBSR (10, 29, 31), which is the one used by TraRe; (ii) LASSO regression, which has been used extensively to relate regulatory genes to their targets (30, 33, 40); and (iii) LM (11, 44), which is the baseline model. For each of these model types, TraRe was run 5 times using a random 80% of the samples to assess its generalization ability. To visualize the inference capabilities of TraRe for each model, we generated alluvial plots where the leftmost groups are the true simulated transcriptional modules, the middle column shows the inferred submodules across the different runs, and the column on the right shows the final inferred transcriptional modules (Fig. 6AC; see Materials and Methods). VBSR outperformed LM and LASSO, as TraRe was able to correctly infer simulated transcriptional modules (Supplementary Note S2).

Figure 6.

TraRe is able to recapitulate rewired transcriptional modules on simulated data. A–C, Alluvial plots showing the flow of true simulated transcriptional modules to the corresponding inferred submodules and the clustering of these into transcriptional modules by similarity for LM, VBSR, and LASSO models. D–F, Heatmaps showing clustered rewired transcriptional modules from labeled-as-rewired simulated inferred submodules (simulated module 1, 4, and 7). G, Models boxplot and paired Wilcoxon test showing the distribution differences significance from the top 50 Jaccard index scores. ***, P < 0.001.

Figure 6.

TraRe is able to recapitulate rewired transcriptional modules on simulated data. A–C, Alluvial plots showing the flow of true simulated transcriptional modules to the corresponding inferred submodules and the clustering of these into transcriptional modules by similarity for LM, VBSR, and LASSO models. D–F, Heatmaps showing clustered rewired transcriptional modules from labeled-as-rewired simulated inferred submodules (simulated module 1, 4, and 7). G, Models boxplot and paired Wilcoxon test showing the distribution differences significance from the top 50 Jaccard index scores. ***, P < 0.001.

Close modal

Furthermore, to evaluate the capability of TraRe to infer rewired modules under different regression models, three of the transcriptional modules were explicitly generated by applying a whitening process to the gene expression data on half of the samples (see Supplementary Note S1). The rewiring test was then applied to the inferred submodules above using a significance threshold of 0.01. The consistency of detection of rewired transcriptional modules was visualized via a heatmap, where a blockwise-diagonal heatmap is expected when the rewired modules are correctly discovered across runs (Fig. 6DF). From the three tested models, VBSR presented the largest similarity (measured by the Jaccard index) to the simulated transcriptional modules (Fig. 6G).

In addition, we evaluated the behavior of VBSR, LASSO, and LM over different noise levels, TF-assigning thresholds, and number of samples (see Supplementary Note S1). For each scenario with a selected combination of these parameters, we plotted the false-positive rate (FPR) to true-positive rate (TPR) curve. VBSR reported the largest TPR/FPR ratios on each of the tested scenarios (Supplementary Note S3 and Supplementary Figs. S7–S10). These findings reinforce the decision to use VBSR as the default model for TraRe analysis.

Finally, we benchmarked on the PROMOTE data the TraRe's GRN inference method against two GRN inference methods: (i) GRNboost2 (37), a known method for GRN discovery, which has also been previously used in single-cell data (45), and (ii) ARACNE-AP (38), a refined version of one of the most used GRN algorithms (21). The assessment was done in two ways: With all the samples (n = 46), to mimic TraRe's inference method, and separately for responder and nonresponder, to assess phenotype-network differences. We tested whether these GRN inference methods were able to infer the relationships that we experimentally validated. Although both methods showed differences between the responder and nonresponder networks, as well as found some of the validated TF-targets relationships, several of these edges showed weak importance values (red vs. dotted lines in Supplementary Figs. S11–S12). Thus, they would be less likely considered for further experimental validation. In addition, we checked the inferred edges against the ChIP-seq data as previously described and found less enrichment compared with TraRe's inferred edges (see Supplementary Note S4 and Supplementary Fig. S3 for additional details).

In this study, we presented TraRe, a computational method to understand mechanistically altered transcriptional dynamics through differential network analysis. We then applied TraRe to the RNA-seq data from the PROMOTE study to understand how differences in transcriptional networks may contribute to abiraterone response in patients with mCRPC.

Differential network analysis can enhance our understanding of network reconfiguration, shedding light on the molecular relationships driving disease progression or clinical treatments (22, 46). Correlation-based estimators have been typically used to analyze gene–gene dependencies within the networks and uncover network disruptions (47, 48). However, these methods are limited to marginal correlation networks (i.e., two nodes at a time) that are calculated separately among response group observations, therefore, lacking the ability to evaluate conserved relationships across multiple groups. To address these limitations, some approaches divide the group-specific conditional dependencies into global and group-specific components, which have been shown to improve performance over previous existing methods (49). However, these methods were proposed for relatively small datasets (hundreds of genes at a time rather than thousands) and do not scale well to larger datasets (50), and therefore, are not suitable for large studies based on RNA-Seq data such as the PROMOTE study presented in this work. Scaling to larger datasets is crucial for the future applicability of rewiring methods. Furthermore, they do not consider the overall rewiring of individual transcriptional networks and instead focus on pairwise gene–gene differential correlations, overlooking key orchestrated interactions among genes. In addition, previous methods to infer transcriptional networks are often applied to elucidate a unique network encompassing the entire transcriptome, making the discovery of individual sub-networks more challenging, which makes them an unfit method for differential network analysis (22).

TraRe uniquely solves these issues by taking a module-based approach for inferring transcriptional networks, which has been recently shown to outperform methods based on inferring a single network (29), and combines it with an efficient, scalable, and robust rewiring inference model that can be applied to large datasets, evidenced by handling large datasets of more than 20K genes. The unique value of TraRe is that it examines phenotypically driven mechanistic differences at three distinct levels: transcriptional networks, specific regulons, and individual TFs.

TraRe, unlike current state-of-the-art methods, goes beyond inferring transcriptional networks, to elucidate the mechanistic transcriptional changes (network rewirings) induced by phenotypic changes. Specifically, TraRe identifies key TFs via VBSR, then uncovers potential transcriptional rewirings associated with a particular phenotype by performing an efficient permutation test on the Frobenius norm of the difference of the estimated covariance matrices of each phenotype. Importantly, TraRe is readily available as an R Bioconductor package facilitating its usage by the broad biomedical community. Although TraRe achieves impressive results on binary phenotype comparison, it currently lacks the capabilities for multiple conditions or phenotypes to be evaluated at the same time. In addition, TraRe is currently designed to work solely on RNA-Seq data and, thus, cannot leverage other data types while inferring the GRNs. This would allow for potentially more robust GRNs to be identified, at the cost, however, of not being able to use the data for evaluation. For instance, ChIP-seq data could be given as prior information to the framework, establishing a peak-based network over which TraRe could build GRNs. Finally, future improvements of TraRe include additional functionalities to integrate novel GRN inference algorithms into the pipeline before the assessment of the differential network analysis.

When applied to the PROMOTE mCRPC dataset, TraRe highlighted key transcriptional modules associated with essential transcriptional dynamics. These transcriptional modules were also recapitulated when applying TraRe to the independent SU2C mCRPC study. More importantly, TraRe identified a transcriptional module (P-M1) associated with immune response that was highly enriched in rewired regulons between responders and nonresponders of Abi treatment. When examining the TFs driving the global expression patterns of this rewired transcriptional module (Supplementary Table S3), we found that TFs CEBPE, GATA1, KLF1, and MYB are the master regulators of granulocytes differentiation (46), erythroid development (47), and other hematopoietic and immune pathways (48). Another of its TFs, TAL1, was strongly differentially expressed and showed an activating relationship with targets in both patient groups whereas the TF NFE2, although not significantly differentially expressed itself, led to strong differential expression of downstream target genes, demonstrating a mostly activating role. On the other hand, the TFs GLI1 and MXD1 were also not significantly differentially expressed but showed a significant loss of correlations with their targets in the NRs. In addition, TFs driving the other rewired transcriptional modules, have been found to have important implications in cancer such as the SREBF2 (P-M2) that is overexpressed in more aggressive prostate cancer cell lines and metastatic prostate cancer (49); SMAD7 (P-M3) that was shown to also promote migration and invasion in prostate cancer cells (51); SOX8 (P-M3) that is known to be involved in cisplatin-chemoresistance in different cancers (52, 53) or SNAI2 (P-M4) that is involved in epithelial–mesenchymal transitions (EMT) in different cancer types, has anti-apoptotic activity (54) and it is a metastasis-promoting TF in breast cancer progression (55). These findings show the complexity of the transcriptional landscape of metastatic cancer tissue and how TFs may have distinct roles in phenotype-specific network rewiring.

AR signaling is essential for the correct development and function of both normal prostate tissue and cancer cells. As such, the AR pathway is the most important target for the mainstay treatment in these tumors and its reactivation through different mechanisms, such as increased expression of ligand-independent transcription variants, is the principal cause of relapse and disease progression (15). Moreover, the glucocorticoid receptor (GR) and AR present similar characteristics in terms of DNA response elements, sharing a set of target genes, and its activation has been described in enzalutamide-exposed xenograft models (56). However, because of incompatible methods of quantification, AR variants cannot be included in TraRe as independent TFs. We did not find GR to be directly implicated in the rewiring of the transcriptional modules identified in our study. One possible reason is that their inferred targets are largely shared between responders and nonresponders as a baseline common network before Abi treatment. In fact, prostate cancer is a complex disease that presents heterogeneous resistant mechanisms, we had identified a number of different mechanisms from our PROMOTE cohorts in our previous publications, including but not limited to the expression of AR variants, elevated expression of Wnt pathways (18), and overexpression of cell-cycle progression genes (57). In addition, leveraging TraRe's inferred information, we identified and validated three TFs ELK3, MYB, and MXD1, whose expressions were not significantly different between Abi R and NR; however, their regulons were significantly disrupted in NR. ELK1, a key TF of PM-2 rewired transcriptional module is a nuclear target of the MAPK signaling cascade and has a potential role in the activation of AR signaling of growth genes in prostate cancer (58). The regulons driven by these TFs in responders probably played important roles in tumorigenesis or progression in the Abi-sensitive samples. For example, we showed that knockdown of the ETS family gene ELK3 suppressed the proliferation of Abi-sensitive prostate cancer cell lines and, more importantly, inhibited migration by suppressing the expression of genes related to focal adhesion and motility, including RAPGEF1 (59), ACTN4 (60), RAB11B (61), MYO18A (62), and PIK3CA (63).

These genes have been previously linked to RAC1–PAK1-mediated E-cadherin stability, EMT, FAK-induced invasion, and metalloproteinase expression. This finding is consistent with Mao and colleagues (64), who showed that silencing of ELK3 in prostate cancer cell lines induced S–M phase arrest, inhibited cell proliferation and migration. However, these associations were not observed in Abi-resistant settings, probably due to the extensive reprogramming of transcriptional networks in the development of the Abi-resistant phenotype. This may help explain why ETS family and ETS-fusion, despite their high prevalence in prostate cancer (65), are controversial in their prognostic values, especially in metastatic settings (66).

Similarly, we found that knockdown of either MYB and MXD1 suppressed proliferation of Abi-sensitive cell lines. MYB is known to be overexpressed in prostate cancer cells. It interacts with AR and sustains AR activity under androgen-depleted conditions (67). MXD1 usually functions as antagonist of MYC–MAX signaling, but it was also reported to mediate HIF1α–induced PI3KAKT activation and chemoresistance in U2OS and MG-63 cells (68). Interestingly, we found that MYB and MXD1 shared a number of genes in their regulons, such as DMTN and MSRB1, the prior of which regulates actin cytoskeletal organization (69), whereas the latter responds to oxidative stress (70). This crosstalk of regulons indicated potential co-regulations between the two TFs, which was not widely reported previously. Again, these signaling cascades have been disrupted in the development of Abi-resistance, probably due to the reprogramming of the transcriptional networks. All these results emphasized the difficulties of therapeutically targeting TFs due to their dynamic regulatory adaption.

From the prostate cancer research perspective, the study presented here delved into the regulatory molecular mechanisms of the abiraterone response in patients with mCRPC, who unfortunately present severe prognoses. We showed that the clinical data provided by the PROMOTE study holds highly valuable information that our method TraRe was able to exploit. However, RNA-seq data provide a mixed signal from the whole-tissue sample, thus it is difficult to distinguish differential signals for a specific cell type (malignant cells vs. microenvironment) comprising the sample. Moreover, although the in vitro studies were able to shed some light and validate some of the hypotheses uncovered by TraRe, studies with more complex models such as tumor-microenvironment coculture and in vivo patient-derived xenografts models are needed to shed further light on regulatory molecular mechanisms of the abiraterone response.

In summary, we envision TraRe as a promising tool for the generation of novel hypotheses based on identified disrupted transcriptional dynamics, and with further refinement could ultimately contribute to the design of personalized therapies.

No disclosures were reported.

C. Blatti: Conceptualization, software, investigation, visualization, methodology, writing–original draft, writing–review and editing. J. de la Fuente: Software, investigation, visualization, methodology, writing–original draft, writing–review and editing. H. Gao: Validation, investigation, visualization, writing–original draft, writing–review and editing. I. Marín-Goñi: Software, investigation, visualization, validation, writing–original draft, writing–review and editing. Z. Chen: Software, investigation, visualization, writing–review and editing. S.D. Zhao: Methodology, writing–review and editing. W. Tan: Supervision, methodology, writing–review and editing. R. Weinshilboum: Supervision, writing–review and editing. K.R. Kalari: Supervision, writing–review and editing. L. Wang: Conceptualization, supervision, validation, writing–review and editing. M. Hernaez: Conceptualization, supervision, methodology, writing–original draft, writing–review and editing.

This work was supported in part by: The Department of Defense of the US—Congressionally Directed Medical Research Programs (grant number W81XWH-20–1-0262 to L. Wang and M. Hernaez); Mayo Clinic Center for Individualized Medicine (MC1351 to L. Wang); A.T. Suharya and D.H. Ghan (to L. Wang); Gail and Joseph Gassner; Marie Skłodowska—Curie Individual Fellowships (grant number 898356 to M. Hernaez); Ayudas Predoctorales Gobierno de Navarra (grant number 0011–0537–2021–000106 to I. Marín-Goñi); Fulbright Predoctoral Research Program (PS00342367 to J. de la Fuente). Funding for open access charge: Department of Defense of the US—Congressionally Directed Medical Research Programs (grant number W81XWH-20–1-0262 to L. Wang and M. Hernaez).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

1.
Siegel
RL
,
Miller
KD
,
Jemal
A
.
Cancer statistics, 2020
.
CA Cancer J Clin
2020
;
70
:
7
30
.
2.
Litwin
MS
,
Tan
HJ
.
The diagnosis and treatment of prostate cancer: a review
.
JAMA
2017
;
317
:
2532
42
.
3.
James
ND
,
de Bono
JS
,
Spears
MR
,
Clarke
N
,
Mason
MD
,
Dearnaley
DP
, et al
.
Adding abiraterone for men with high-risk prostate cancer (PCa) starting long-term androgen deprivation therapy (ADT): survival results from STAMPEDE (NCT00268476)
.
J Clin Oncol
2017
;
35
.
doi: 10.1200/JCO.2017.35.18_suppl.LBA5003
.
4.
Fizazi
K
,
Tran
N
,
Fein
LE
,
Matsubara
N
,
Rodriguez Antolin
A
,
Alekseev
BY
, et al
.
LATITUDE: a phase III, double-blind, randomized trial of androgen deprivation therapy with abiraterone acetate plus prednisone or placebos in newly diagnosed high-risk metastatic hormone-naive prostate cancer
.
J Clin Oncol
2017
;
35
.
doi: 10.1200/JCO.2017.35.18_suppl.LBA3
.
5.
Petrylak
D
.
Therapeutic options in androgen-independent prostate cancer: building on docetaxel
.
BJU Int
2005
;
96
:
41
6
.
6.
de Bono
JS
,
Logothetis
CJ
,
Molina
A
,
Fizazi
K
,
North
S
,
Chu
L
, et al
.
Abiraterone and increased survival in metastatic prostate cancer
.
N Engl J Med
2011
;
364
:
1995
2005
.
7.
Scher
HI
,
Fizazi
K
,
Saad
F
,
Taplin
ME
,
Sternberg
CN
,
Miller
K
, et al
.
Increased survival with enzalutamide in prostate cancer after chemotherapy
.
N Engl J Med
2012
;
367
:
1187
97
.
8.
Danila
DC
,
Morris
MJ
,
de Bono
JS
,
Ryan
CJ
,
Denmeade
SR
,
Smith
MR
, et al
.
Phase II multicenter study of abiraterone acetate plus prednisone therapy in patients with docetaxel-treated castration-resistant prostate cancer
.
J Clin Oncol
2010
;
28
:
1496
.
9.
Boehm
JS
,
Hahn
WC
.
Towards systematic functional characterization of cancer genomes
.
Nat Rev Genet
2011
;
12
:
487
.
10.
Logsdon
BA
,
Gentles
AJ
,
Miller
CP
,
Blau
CA
,
Becker
PS
,
Lee
SI
.
Sparse expression bases in cancer reveal tumor drivers
.
Nucleic Acids Res
2015
;
43
:
1332
44
.
11.
Akavia
UD
,
Litvin
O
,
Kim
J
,
Sanchez-Garcia
F
,
Kotliar
D
,
Causton
HC
, et al
.
An integrated approach to uncover drivers of cancer
.
Cell
2010
;
143
:
1005
17
.
12.
Bandyopadhyay
S
,
Mehta
M
,
Kuo
D
,
Sung
MK
,
Chuang
R
,
Jaehnig
EJ
, et al
.
Rewiring of genetic networks in response to DNA damage
.
Science
2010
;
330
:
1385
9
.
13.
Califano
A
.
Rewiring makes the difference
.
Mol Syst Biol
2011
;
7
:
463
.
14.
Luscombe
NM
,
Babu
MM
,
Yu
H
,
Snyder
M
,
Teichmann
SA
,
Gerstein
M
.
Genomic analysis of regulatory network dynamics reveals large topological changes
.
Nature
2004
;
431
:
308
.
15.
Feng
Q
,
He
B
.
Androgen receptor signaling in the development of castration-resistant prostate cancer
.
Front Oncol
2019
;
9
:
858
.
16.
Sharma
NL
,
Massie
CE
,
Ramos-Montoya
A
,
Zecchini
V
,
Scott
HE
,
Lamb
AD
, et al
.
The androgen receptor induces a distinct transcriptional program in castration-resistant prostate cancer in man
.
Cancer Cell
2013
;
23
:
35
47
.
17.
Augello
MA
,
Liu
D
,
Deonarine
LD
,
Robinson
BD
,
Huang
D
,
Stelloo
S
, et al
.
CHD1 loss alters AR binding at lineage-specific enhancers and modulates distinct transcriptional programs to drive prostate tumorigenesis
.
Cancer Cell
2019
;
35
:
603
17
.
18.
Wang
L
,
Dehm
SM
,
Hillman
DW
,
Sicotte
H
,
Tan
W
,
Gormley
M
, et al
.
A prospective genome-wide study of prostate cancer metastases reveals association of wnt pathway activation and increased cell cycle proliferation with primary resistance to abiraterone acetate-prednisone
.
Ann Oncol
2018 Feb 1
;
29
:
352
60
.
19.
la Fuente
A
.
From ‘differential expression'to ‘differential networking’–identification of dysfunctional regulatory networks in diseases
.
Trends Genet
2010
;
26
:
326
33
.
20.
Hudson
NJ
,
Dalrymple
BP
,
Reverter
A
.
Beyond differential expression: the quest for causal mutations and effector molecules
.
Bmc Genomics
2012
;
13
:
1
16
.
21.
Margolin
AA
,
Nemenman
I
,
Basso
K
,
Wiggins
C
,
Stolovitzky
G
,
Favera
RD
, et al
.
ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context
.
BMC Bioinf
2006
;
7
.
Available from
: https://doi.org/10.1186/1471-2105-7-S1-S7.
22.
Califano
A
,
Alvarez
MJ
.
The recurrent architecture of tumour initiation, progression and drug sensitivity
.
Nat Rev Cancer
2017
;
17
:
116
.
23.
Sinha
S
,
Jones
BM
,
Traniello
IM
,
Bukhari
SA
,
Halfon
MS
,
Hofmann
HA
, et al
.
Behavior-related gene regulatory networks: a new level of organization in the brain
.
Proc Natl Acad Sci
2020
;
117
:
23270
9
.
24.
Vaquerizas
JM
,
Kummerfeld
SK
,
Teichmann
SA
,
Luscombe
NM
.
A census of human transcription factors: function, expression and evolution
.
Nat Rev Genet
2009
;
10
:
252
63
.
25.
Rhinn
H
,
Fujita
R
,
Qiang
L
,
Cheng
R
,
Lee
JH
,
Abeliovich
A
.
Integrative genomics identifies APOE 4 effectors in Alzheimer's disease
.
Nature
2013
;
500
:
45
.
26.
Taylor
IW
,
Linding
R
,
Warde-Farley
D
,
Liu
Y
,
Pesquita
C
,
Faria
D
, et al
.
Dynamic modularity in protein interaction networks predicts breast cancer outcome
.
Nat Biotechnol
2009
;
27
:
199
.
27.
Lai
Y
,
Wu
B
,
Chen
L
,
Zhao
H
.
A statistical method for identifying differential gene–gene co-expression patterns
.
Bioinformatics
2004
;
20
:
3146
55
.
28.
de La Fuente Cedeño
J
,
Marín
I
,
Hernaez
M
,
Blatti
C
.
TraRe: Transcriptional Rewiring
.
2021
.
Available from
: https://github.com/ubioinformat/TraRe/tree/master.
29.
Hernaez
M
,
Blatti
C
,
Gevaert
O
.
Comparison of single and module-based methods for modeling gene regulatory networks
.
Bioinformatics
2020
;
36
:
558
67
.
30.
Gevaert
O
,
van Vooren
S
,
de Moor
B
.
A framework for elucidating regulatory networks based on prior information and expression data
.
Ann N Y Acad Sci
2007
;
1115
:
240
8
.
31.
Logsdon
BA
,
Carty
CL
,
Reiner
AP
,
Dai
JY
,
Kooperberg
C
.
A novel variational Bayes multiple locus Z-statistic for genome-wide association studies with Bayesian model averaging
.
Bioinformatics
2012
;
28
:
1738
44
.
32.
Abida
W
,
Cyrta
J
,
Heller
G
,
Prandi
D
,
Armenia
J
,
Coleman
I
, et al
.
Genomic correlates of clinical outcome in advanced prostate cancer
.
Proc Natl Acad Sci
USA
2019
;
116
:
11428
36
.
33.
Champion
M
,
Brennan
K
,
Croonenborghs
T
,
Gentles
AJ
,
Pochet
N
,
Gevaert
O
.
Module analysis captures pancancer genetically and epigenetically deregulated cancer driver genes for smoking and antiviral response
.
EBioMedicine
2018
;
27
:
156
66
.
34.
Newman
MEJ
,
Girvan
M
.
Finding and evaluating community structure in networks
.
Phys Rev E
2004
;
69
:
26113
.
35.
Liberzon
A
,
Birger
C
,
Thorvaldsdóttir
H
,
Ghandi
M
,
Mesirov
JP
,
Tamayo
P
.
The molecular signatures database hallmark gene set collection
.
Cell Syst
2015
;
1
:
417
25
.
36.
Hammal
F
,
de Langen
P
,
Bergon
A
,
Lopez
F
,
BB
ReMap
.
2022: a database of Human, Mouse, Drosophila and Arabidopsis regulatory regions from an integrative analysis of DNA-binding sequencing experiments
.
Nucleic Acids Res
2022
;
50
:
D316
25
.
37.
Moerman
T
,
Aibar Santos
S
,
Bravo González-Blas
C
,
Simm
J
,
Moreau
Y
,
Aerts
J
, et al
.
GRNBoost2 and Arboreto: efficient and scalable inference of gene regulatory networks
.
Bioinformatics
2018
;
35
:
2159
61
.
38.
Lachmann
A
,
Giorgi
FM
,
Lopez
G
,
Califano
A
.
ARACNe-AP: gene network reverse engineering through adaptive partitioning inference of mutual information
.
Bioinformatics
2016
;
32
:
2233
5
.
39.
Gevaert
O
,
Villalobos
V
,
Sikic
BI
,
Plevritis
SK
.
Identification of ovarian cancer driver genes by using module network integration of multi-omics data. Interface Focus
.
2013
;
3
:
20130013
.
40.
Manolakos
A
,
Ochoa
I
,
Venkat
K
,
Goldsmith
AJ
,
Gevaert
O
.
CaMoDi: a new method for cancer module discovery
.
Bmc Genomics
2014
;
15
:
S8
.
41.
Helgeson
BE
,
Tomlins
SA
,
Shah
N
,
Laxman
B
,
Cao
Q
,
Prensner
JR
, et al
.
Characterization of TMPRSS2: ETV5 and SLC45A3: ETV5 gene fusions in prostate cancer
.
Cancer Res
2008
;
68
:
73
80
.
42.
Szklarczyk
D
,
Gable
AL
,
Lyon
D
,
Junge
A
,
Wyder
S
,
Huerta-Cepas
J
, et al
.
STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets
.
Nucleic Acids Res
2018
;
47
:
D607
13
.
43.
Vagapova
ER
,
Spirin
PV
,
Lebedev
TD
,
Prassolov
VS
.
The role of TAL1 in hematopoiesis and leukemogenesis
.
Acta Naturae
2018
;
10
:
15
23
.
44.
Langfelder
P
,
Horvath
S
.
WGCNA: an R package for weighted correlation network analysis
.
BMC Bioinf
2008
;
9
:
1
13
.
45.
Aibar
S
,
González-Blas
CB
,
Moerman
T
,
Huynh-Thu
VA
,
Imrichova
H
,
Hulselmans
G
, et al
.
SCENIC: single-cell regulatory network inference and clustering
.
Nat Methods
2017
;
14
:
1083
6
.
46.
Shyamsunder
P
,
Shanmugasundaram
M
,
Mayakonda
A
,
Dakle
P
,
Teoh
WW
,
Han
L
, et al
.
Identification of a novel enhancer of CEBPE essential for granulocytic differentiation
.
Blood
2019
;
133
:
2507
17
.
47.
Gilles
L
,
Arslan
AD
,
Marinaccio
C
,
Wen
QJ
,
Arya
P
,
McNulty
M
, et al
.
Downregulation of GATA1 drives impaired hematopoiesis in primary myelofibrosis
.
J Clin Invest
2017
;
127
:
1316
20
.
48.
Gautam
S
,
Fioravanti
J
,
Zhu
W
,
le Gall
JB
,
Brohawn
P
,
Lacey
NE
, et al
.
The transcription factor c-Myb regulates CD8+ T-cell stemness and antitumor immunity
.
Nat Immunol
2019
;
20
:
337
49
.
49.
Li
X
,
Wu
JB
,
Li
Q
,
Shigemura
K
,
Chung
LWK
,
Huang
WC
.
SREBP-2 promotes stem cell-like properties and metastasis by transcriptional activation of c-Myc in prostate cancer
.
Oncotarget
2016
;
7
:
12869
.
50.
Yu
D
,
Lim
J
,
Wang
X
,
Liang
F
,
Xiao
G
.
Enhanced construction of gene regulatory networks using hub gene information
.
BMC Bioinf
2017
;
18
:
186
.
51.
Thakur
N
,
Hamidi
A
,
Song
J
,
Itoh
S
,
Bergh
A
,
Heldin
CH
, et al
.
Smad7 enhances TGFβ-induced transcription of c-Jun and HDAC6 promoting invasion of prostate cancer cells
.
iScience
2020
;
23
:
101470
.
52.
Sun
H
,
Wang
H
,
Wang
X
,
Aoki
Y
,
Wang
X
,
Yang
Y
, et al
.
Aurora-A/SOX8/FOXK1 signaling axis promotes chemoresistance via suppression of cell senescence and induction of glucose metabolism in ovarian cancer organoids and cells
.
Theranostics
2020
;
10
:
6928
.
53.
Xie
SL
,
Fan
S
,
Zhang
SY
,
Chen
WX
,
Li
QX
,
Pan
GK
, et al
.
SOX8 regulates cancer stem-like properties and cisplatin-induced EMT in tongue squamous cell carcinoma by acting on the Wnt/β-catenin pathway
.
Int J Cancer
2018
;
142
:
1252
65
.
54.
Kim
S
,
Yao
J
,
Suyama
K
,
Qian
X
,
Qian
BZ
,
Bandyopadhyay
S
, et al
.
Slug promotes survival during metastasis through suppression of Puma-mediated apoptosis
.
Cancer Res
2014
;
74
:
3695
706
.
55.
Li
W
,
Shen
M
,
Jiang
YZ
,
Zhang
R
,
Zheng
H
,
Wei
Y
, et al
.
Deubiquitinase USP20 promotes breast cancer metastasis by stabilizing SNAI2
.
Genes Dev
2020
;
34
:
1310
5
.
56.
Arora
VK
,
Schenkein
E
,
Murali
R
,
Subudhi
SK
,
Wongvipat
J
,
Balbas
MD
, et al
.
Glucocorticoid receptor confers resistance to antiandrogens by bypassing androgen receptor blockade
.
Cell
2013
;
155
:
1309
22
.
57.
Qin
S
,
Gao
H
,
Kim
W
,
Zhang
H
,
Gu
Y
,
Kalari
KR
, et al
.
Biomarkers for predicting abiraterone treatment outcome and selecting alternative therapies in castration-resistant prostate cancer
.
Clin Pharmacol Ther
2022
;
111
:
1296
306
.
58.
Rosati
R
,
Patki
M
,
Chari
V
,
Dakshnamurthy
S
,
McFall
T
,
Saxton
J
, et al
.
The amino-terminal domain of the androgen receptor co-opts extracellular signal-regulated kinase (ERK) docking sites in ELK1 protein to induce sustained gene activation that supports prostate cancer cell growth
.
J Biol Chem
2016
;
291
:
25983
98
.
59.
Ohba
Y
,
Ikuta
K
,
Ogura
A
,
Matsuda
J
,
Mochizuki
N
,
Nagashima
K
, et al
.
Requirement for C3G-dependent Rap1 activation for cell adhesion and embryogenesis
.
EMBO J
2001
;
20
:
3333
41
.
60.
Nakatsuji
H
,
Nishimura
N
,
Yamamura
R
,
KH
omi
,
Sasaki
T
.
Involvement of actinin-4 in the recruitment of JRAB/MICAL-L2 to cell–cell junctions and the formation of functional tight junctions
.
Mol Cell Biol
2008
;
28
:
3324
35
.
61.
Erasmus
JC
,
Smolarczyk
K
,
Brezovjakova
H
,
Mohd-Naim
NF
,
Lozano
E
,
Matter
K
, et al
.
Rac1-PAK1 regulation of Rab11 cycling promotes junction destabilization
.
J Cell Biol
2021
;
220
:
e202002114
.
62.
Liu
J
,
Ren
G
,
Li
K
,
Liu
Z
,
Wang
Y
,
Chen
T
, et al
.
The Smad4-MYO18A-PP1A complex regulates β-catenin phosphorylation and pemigatinib resistance by inhibiting PAK1 in cholangiocarcinoma
.
Cell Death Differ
2022
;
29
:
818
31
.
63.
Zeng
ZZ
,
Jia
Y
,
Hahn
NJ
,
Markwart
SM
,
Rockwood
KF
,
Livant
DL
.
Role of focal adhesion kinase and phosphatidylinositol 3′-kinase in integrin fibronectin receptor-mediated, matrix metalloproteinase-1–dependent invasion by metastatic prostate cancer cells
.
Cancer Res
2006
;
66
:
8091
9
.
64.
Mao
Y
,
Li
W
,
Hua
B
,
Gu
X
,
Pan
W
,
Chen
Q
, et al
.
Silencing of ELK3 induces SM phase arrest and apoptosis and upregulates SERPINE1 expression reducing migration in prostate cancer cells
.
Biomed Res Int
2020
;
2020
:
2406159
.
65.
Abeshouse
A
,
Ahn
J
,
Akbani
R
,
Ally
A
,
Amin
S
,
Andry
CD
, et al
.
The molecular taxonomy of primary prostate cancer
.
Cell
2015
;
163
:
1011
25
.
66.
Clark
JP
,
Cooper
CS
.
ETS gene fusions in prostate cancer
.
Nat Rev Urol
2009
;
6
:
429
39
.
67.
Drabsch
Y
,
Ramsay
RG
,
Gonda
TJ
.
MYB suppresses differentiation and apoptosis of human breast cancer cells
.
Breast Cancer Res
2010
;
12
:
1
17
.
68.
Zheng
D
,
Wu
W
,
Dong
N
,
Jiang
X
,
Xu
J
,
Zhan
X
, et al
.
Mxd1 mediates hypoxia-induced cisplatin resistance in osteosarcoma cells by repression of the PTEN tumor-suppressor gene
.
Mol Carcinog
2017
;
56
:
2234
44
.
69.
Ye
YP
,
Jiao
HL
,
Wang
SY
,
Xiao
ZY
,
Zhang
D
,
Qiu
JF
, et al
.
Hypermethylation of DMTN promotes the metastasis of colorectal cancer cells by regulating the actin cytoskeleton through Rac1 signaling activation
.
J Exp Clin Cancer Res
2018
;
37
:
1
14
.
70.
Novoselov
S v
,
Kim
HY
,
Hua
D
,
Lee
BC
,
Astle
CM
,
Harrison
DE
, et al
.
Regulation of selenoproteins and methionine sulfoxide reductases A and B1 by age, calorie restriction, and dietary selenium in mice
.
Antioxid Redox Signal
2010
;
12
:
829
38
.
This open access article is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.