RNA-Seq Uncovers Association of Endocrine-Disrupting Chemicals with Hub Genes and Transcription Factors in Aggressive Prostate Cancer

Abstract

This study analyzes publicly available RNA-seq data to comprehensively include the complex heterogeneity of prostate cancer (PCa) etiology. It combines prostate and prostate cancer (PCa) cell lines, representing primary PCa cells, Gleason scores, ages, and PCa of different racial origins. Additionally, some cell lines were exposed to endocrine-disrupting chemicals (EDCs). The research aims to identify hub genes and transcription factors (TFs) of the prostate carcinogenesis pathway as molecular targets for clinical investigations to elucidate EDC-induced aggressiveness and to develop potential biomarkers for their exposure risk assessments. PCa cells rely on androgen receptor (AR)-mediated signaling to survive, develop, and function. Fifteen various RNA-seq datasets were normalized for distribution, and the significance (p-value < 0.05) threshold of differentially expressed genes (DEGs) was set based on |log2FC| ≥ 2 change. Through integrated bioinformatics, we applied cBioPortal, UCSC-Xena, TIMER2.0, and TRRUST platforms, among others, to associate hub genes and their TFs based on their biologically meaningful roles in aggressive prostate carcinogenesis. Among all RNA-Seq datasets, we found 75 overlapping DEGs, with BUB1B (32%) and CCNB1 (29%) genes exhibiting the highest degree of mutation, amplification, and deletion. EDC-associated CCNB1, BUB1B, and CCNA2 in PCa cells exposed to EDCs were consistently shown to be associated with high Gleason scores (≥4 + 3) and in the >60 age group of patients. Selected TFs (E2F4, MYC, and YBX1) were also significantly associated with DEGs (NCAPG, MKI67, CCNA2, CCNB1, CDK1, CCNB2, AURKA, UBE2C, BUB1B) and influenced the overall survival (p-value < 0.05) of PCa cases. This is one of the first comprehensive studies combining 15 publicly available RNA-seq datasets to demonstrate the association of EDC-associated hub genes and their TFs aligning with the aggressive carcinogenic pathways in the higher age group (>60 years) of patients. The findings highlight the potential of these hub genes as candidates for further studies to develop molecular biomarkers for assessing the EDC-related PCa risk, diagnosing PCa aggressiveness, and identifying therapeutic targets.

Keywords: RNA-seq; androgen receptor; differentially expressed genes; endocrine-disrupting chemicals; environmental carcinogenesis; environmental health risk assessment; prostate cancer; transcription factors.

Figure 1

Venn diagram analysis of 4 groups of DEGs in PCa cell lines across RNA-seq datasets represented in blue (A), yellow (B), green (C), and red (D) panels. Panel (A) displays 418 DEGs common to cell lines treated with R1881 treatment (AR agonist). Panel (B) shows 116 DEGs shared among cell lines with EDC exposures, which also included AR modulators. Panel (C) illustrates 2235 DEGs shared between Gleason scores and cells overexpressing MYC (TF) for functional analysis. Panel (D) represents 7472 DEGs overlapping between age and race. A total of 75 (0.8%) overlapping genes are among the RNA-seq datasets.

Figure 2

PPI network analysis of 75 DEGs using STRING. (A) features a network of 38 upregulated genes (nodes) with 366 edges, an average node degree of 19.3, and a local clustering coefficient of 0.78. (B) shows a network of 37 downregulated genes with 282 edges, an average node degree of 15.2, and a local clustering coefficient of 0.88. These networks are derived from RNA-seq datasets with PPI enrichment p-value < 1.0⁻¹⁶.

Figure 3

The figure illustrates the GeneMANIA-generated networks of hub genes (orange circles), their associated genes (blue circles), and the weight percentages of various interaction types (physical, co-expression, co-localization, genetic interaction, pathway, and shared protein domain). (A) shows the networks of 20 upregulated hub genes, and (B) shows 20 downregulated hub genes, each ranked and sized according to their positive real-valued link weight score.

Figure 3

The figure illustrates the GeneMANIA-generated networks of hub genes (orange circles), their associated genes (blue circles), and the weight percentages of various interaction types (physical, co-expression, co-localization, genetic interaction, pathway, and shared protein domain). (A) shows the networks of 20 upregulated hub genes, and (B) shows 20 downregulated hub genes, each ranked and sized according to their positive real-valued link weight score.

Figure 4

Genetic alteration analysis of hub genes from CBioPortal. The figure summarizes genetic alterations in upregulated (A) and downregulated (B) hub genes based on data from CBioPortal, encompassing 26 studies with 10,998 samples. Genetic alteration: red is amplification, blue is deep deletion, light red is high mRNA, light blue is low mRNA, and gray is no alteration. Mutation spectrum: adenine (A), cytosine (C), guanine (G), and thymine (T).

Figure 5

Heatmap analysis of 20 hub genes from TCGA, ICGC, and GTEx. The heatmap visualizes the expression of 20 hub genes, with data sourced from TCGA, ICGC, and GTEx via UCSC-Xena. Color saturation levels correspond to maximum and minimum log2-transformed expression values (red for high expression, blue for low). Gleason scores range from 6 (black) to 10 (yellow), and ages range from 41 (dark gray) to 78 (light gray). Panel (A) represents the 10 upregulated and Panel (B) represents the 10 downregulated groups of hub genes.

Figure 6

Network construction of protein-encoding TF genes using TRRUST V.2 database. (A): Functional E2F4, MYC, and YBX1 gene networks in HumanNet, constructed through modified Bayesian integration. (B): Protein–protein interaction (PPI) networks involving E2F4, MYC, and YBX1, mapped using iRefIndex.