Molecular signatures identified by integrating gene expression and methylation in non-seminoma and seminoma of testicular germ cell tumours

ABSTRACT

Testicular germ cell tumours (TGCTs) are the most common cancer in young male adults (aged 15 to 40). Unlike most other cancer types, identification of molecular signatures in TGCT has rarely reported. In this study, we developed a novel integrative analysis framework to identify co-methylated and co-expressed genes [mRNAs and microRNAs (miRNAs)] modules in two TGCT subtypes: non-seminoma (NSE) and seminoma (SE). We first integrated DNA methylation and mRNA/miRNA expression data and then used a statistical method, CoMEx (Combined score of DNA Methylation and Expression), to assess differentially expressed and methylated (DEM) genes/miRNAs. Next, we identified co-methylation and co-expression modules by applying WGCNA (Weighted Gene Correlation Network Analysis) tool to these DEM genes/miRNAs. The module with the highest average Pearson’s Correlation Coefficient (PCC) after considering all pair-wise molecules (genes/miRNAs) included 91 molecules. By integrating both transcription factor and miRNA regulations, we constructed subtype-specific regulatory networks for NSE and SE. We identified four hub miRNAs (miR-182-5p, miR-520b, miR-520c-3p, and miR-7-5p), two hub TFs (MYC and SP1), and two genes (RECK and TERT) in the NSE-specific regulatory network, and two hub miRNAs (miR-182-5p and miR-338-3p), five hub TFs (ETS1, HIF1A, HNF1A, MYC, and SP1), and three hub genes (CDH1, CXCR4, and SNAI1) in the SE-specific regulatory network. miRNA (miR-182-5p) and two TFs (MYC and SP1) were common hubs of NSE and SE. We further examined pathways enriched in these subtype-specific networks. Our study provides a comprehensive view of the molecular signatures and co-regulation in two TGCT subtypes.

KEYWORDS:

• Testicular germ cell tumour
• non-seminoma
• seminoma
• microRNA
• methylation
• co-regulation network

Introduction

Testicular germ cell tumours (TGCTs) are the most common cancer in young men (ages 18–45), accounting for about 1% of the male population, and more than 50,000 new cases every year [1–4]. TGCTs can be divided into two subtypes: seminomas (SE) and non-seminomas (NSE), where NSE patients were often diagnosed as mixed subtypes, including SE and subtypes of NSE [1,3,4]. It is of great importance to explore the genomic mechanisms of TGCT to help the decision-making of the treatment of TGCT. Therefore, a number of studies have been performed to analyse the molecular characterization of TGCT. Litchfield et al. discussed the genomic features of TGCT by using a longitudinal approach, from cancer predisposition to long-term survival status, and collected 25 risk loci of TGCT [5]. Later, they identified 19 new risk loci based on the analysis of Genome-Wide Association Studies (GWAS) of a large number of TGCT samples and controls [2]. The epigenetics of TGCTs have also been analysed, including DNA methylation, histone modification, and RNA interference, and furthermore, some DNA methylation-based biomarkers were identified, such as APC, p16, p14, miR-373-3p, and so on [4,6,7]. In our recent work, we performed a systematic investigation of molecular profiles related to telomere length and FFL models by using the molecular expression profiles of TGCT samples from The Cancer Genome Atlas (TCGA) and described it in terms of NSE and SE [8,9]. Recently, Shen et al. studied the molecular characterization of TGCT using multi-omics data, including mRNA, miRNA expression, DNA methylation, copy number variation as well as somatic mutation, and found that DNA methylation and miRNA expression profiles were able to distinguish different subtypes. Three genes (KIT, KRAS, and NRAS) were identified significantly mutated in SE samples [1].

On the other hand, as one of the most studied epigenetic events, high DNA methylation is associated with reduced gene expression, and DNA methylation change in the genes has been considered as a biomarker of complex diseases, including cancer [10]. However, there are problems to be solved, for example, what is the exact relationship between DNA methylation and gene expression? And how does DNA methylation affect disease initialization and development? A number of studies have been performed to identify DNA methylation patterns. Jiao et al. put forwards a method named FEM (Functional Epigenetic Modules) to identify various gene expression modules by integrating DNA methylation and gene expression data [11]. Zhang et al. proposed an algorithm called lmQCM which can find co-methylation modules from pan-cancer data [12]. In our previous work, we proposed a statistical model which can identify significant genes through integrating DNA methylation profile and gene expression profile, and effectively discover tumour markers [13]. We also integrated and analysed DNA methylation as well as gene expression data by using a rule-based approach [14,15]. Ma et al. detected various epigenetic modules from multiple networks through integrating DNA methylation and mRNA expression profiles that classified breast cancer subtypes successfully [16]. Recently, Hao et al. identified DNA methylation markers of breast cancer, colon cancer, liver cancer, and lung cancer via investigating methylation patterns [17]. Furthermore, other related motif/module methods and molecular network analyses have been applied to disease and found interesting results [18,19]. However, there is still no systematic analysis of DNA methylation patterns of TGCT subtypes.

In this work, we investigated and compared methylation patterns of the two main subtypes of TGCT, NSE versus SE, based on our previous method [13]. DNA methylation profiles and mRNA/miRNA expression data were analysed by using a distance score called INTcon [13]. Then, DEM genes/miRNAs were identified with CoMEx score [13]. Since there were no tumour-matched normal samples for TGCT in TCGA project, differential expressed/methylated genes/miRNAs were defined according to their values in NSE samples versus SE samples. We discovered functional co-methylation and co-expression modules of molecules for these DEM genes/miRNAs based on their combined scores using WGCNA (Weighted correlation network analysis) [20]. PCC (Pearson correlation coefficient) was then computed between each pair of genes belonging to the same module. Next, several cluster validity indices were computed to evaluate the quality of clustering in WGCNA.

In addition, we explored various significant molecules considering their transcriptional and post-transcriptional regulations. For this, we evaluated the TF-target gene pairs and miRNA-target gene pairs with specific PCC cut-off and p-value cut-off. After that, we built subtype-specific regulatory networks and identified various hub molecules of each regulatory network. The common hubs between these two networks were identified. Furthermore, we conducted KEGG pathway analysis and extracted significantly enriched KEGG pathways for each subtype. In conclusion, our method is useful for the identification of various potential modules as well as hubs for any kind of disease.

Results

We performed our analysis according to the pipeline in Figure 1. Firstly, we integrated mRNA/miRNA expression profiles and methylation profiles and computed differentially expressed/methylated mRNAs and miRNAs. Then we performed two parallel lines of analysis for these mRNAs and miRNAs. On one hand, we detected molecular modules using WCGNA and evaluated the most significant module. On the other hand, we constructed and analysed subtype-specific gene regulatory networks by integrated expression profiles/methylation profiles and TF-target/miRNA-target pairs.

Figure 1. Flowchart of the analysis of the integration of DNA methylation and expression data for mRNA as well as miRNA

Differentially expressed and methylated gene/miRNA identification and module detection

CoMEx is a statistical test similar to t-test; it uses the fudge factor for integrating the expression and methylation profiles of the transcripts and then generating DEM transcripts (see Materials and Methods). Using the CoMEx analysis, we identified 2739 DEM genes and 28 DEM miRNAs (CoMEx, p-value<0.05) by comparing experimental samples versus control samples. To explore the functional organization of these significant (DEM) molecules, we used WGCNA to detect gene and miRNA modules. Firstly, the integrated expression levels and methylation values of the significant genes and miRNAs were taken as input for the WGCNA. In this study, we compared NSE with SE to identify TGCT subtype-specific features. According to Figure 2(a), we chose the power 3 that was best fit for the scale-free topology. Twenty-three modules were identified (Figure 2(b)), which demonstrated that the molecules could be effectively organized into modules that integrated expression and DNA methylation values. The topmost module has the highest correlation (‘greenyellow’ module) consisting of 41 molecules. We constructed the co-expression network for this module with PCC cut-off 0.6 (Figure 3). Among the 41 molecules, 39 were up-regulated and 2 were down-regulated (Table 1). Silhouette Width plot for NSE versus SE is shown in Figure 4. The x-axis denotes Silhouette Width ([0, 1]), and the y-axis denotes each individual cluster. The average Silhouette Width is obtained as 0.01. Other clustering validity index measurements such as Clustering Co-efficient, Dunn Index, Heterogeneity, Maximum Adjacency Ratio (MAR), average Scaled Connectivity were 0.023, 0.843, 0.369, 0.026, and 0.34, respectively, as shown in Table 2.

Table 1. Up-regulated and down-regulated genes in ‘greenyellow’ module (NSE versus SE, |logFC| >0.5)

Gene	Regulation type	LogFC (in expression data)	AveExpr	p-value (in expression data)
DRAXIN (Clorf187)	Up-regulated & hypo-methylated	2.60	5.89	1.05E-11
PODXL	Up-regulated & hypo-methylated	2.48	11.82	1.34E-14
ITGB5	Up-regulated & hypo-methylated	1.77	10.95	1.98E-14
NOS1	Up-regulated & hypo-methylated	1.77	2.66	3.15E-08
SLC30A2	Up-regulated & hypo-methylated	1.49	5.84	6.40E-05
EPR1	Up-regulated & hypo-methylated	1.48	9.33	1.4E-4
CD9	Up-regulated & hypo-methylated	1.07	13.08	1.93E-05
ISM2	Up-regulated & hypo-methylated	1.07	6.29	4.35E-05
PGM2L1	Up-regulated & hypo-methylated	1.02	7.94	1.26E-10
HOMER1	Up-regulated & hypo-methylated	1.00	7.49	1.35E-07
EDA2 R	Up-regulated & hypo-methylated	0.98	6.95	6.40E-05
AHCYL1	Up-regulated & hypo-methylated	0.94	11.01	1.02E-17
WBP11P1	Up-regulated & hypo-methylated	0.93	4.34	8.75E-06
RBBP4	Up-regulated & hypo-methylated	0.93	12.29	8.36E-11
HSP90B3P	Up-regulated & hypo-methylated	0.89	3.16	3.11E-06
PGAM4	Up-regulated & hypo-methylated	0.87	6.34	1.11E-04
UBE2 H	Up-regulated & hypo-methylated	0.85	9.95	6.80E-12
SYPL1	Up-regulated & hypo-methylated	0.74	10.57	6.63E-08
KPNA6	Up-regulated & hypo-methylated	0.73	10.9	1.54E-10
UBE4B	Up-regulated & hypo-methylated	0.71	10.35	5.47E-11
MAT2A	Up-regulated & hypo-methylated	0.71	11.88	2.49E-09
ECT2	Up-regulated & hypo-methylated	0.69	10.23	3.41E-04
PARD3	Up-regulated & hypo-methylated	0.68	9.98	1.64E-06
RSC1A1	Up-regulated & hypo-methylated	0.67	7.43	2.20E-06
ING3	Up-regulated & hypo-methylated	0.67	8.31	3.67E-09
CACNG5	Up-regulated & hypo-methylated	0.66	1.06	9.25E-04
DHX15	Up-regulated & hypo-methylated	0.64	11.31	1.97E-10
NLN	Up-regulated & hypo-methylated	0.64	8.64	1.23E-04
PPP2R5E	Up-regulated & hypo-methylated	0.63	9.48	1.22E-07
EIF4G1	Up-regulated & hypo-methylated	0.63	13.37	4.47E-06
ARHGEF5	Up-regulated & hypo-methylated	0.63	8.16	1.15E-05
GTF2H3	Up-regulated & hypo-methylated	0.61	7.66	1.79E-04
PCDH1	Up-regulated & hypo-methylated	0.58	9.93	4.02E-04
AXIN2	Up-regulated & hypo-methylated	0.57	10.06	1.14E-04
PTCD2	Up-regulated & hypo-methylated	0.56	6.79	3.51E-07
HSPA8	Up-regulated & hypo-methylated	0.52	14.73	4.63E-03
USP13	Up-regulated & hypo-methylated	0.52	9.20	7.69E-06
RNF121	Up-regulated & hypo-methylated	0.51	8.65	8.56E-08
INTS2	Up-regulated & hypo-methylated	0.50	8.76	6.88E-04
ZNF438	Down-regulated & hyper-methylated	−0.83	7.44	7.49E-07
FAM41 C	Down-regulated & hyper-methylated	−1.02	4.68	3.55E-07

FC: fold change.

Table 2. Different cluster validity index measures during the module detection through average linkage clustering

Cluster validity index measure	Value
Dunn Index	0.843
Average Scaled Connectivity	0.340
Silhouette Distance	0.007
Average Cluster Co-efficient	0.023
Average Maximum Adjacency Ratio	0.026
Density	0.016
Centralization	0.031
Heterogeneity	0.369

Figure 2. Differentially expressed and methylated (DEM) module analysis in two subtypes: non-seminoma (NSE) and seminoma (SE). a. Soft threshold power for NSE versus SE. The first sub-figure (left) presents soft threshold (power) (in x-axis) versus scale-free topology fit, signed R^2 (in y-axis) plot, where the red horizontal line signifies the lower borderline of the standard value of scale-free topology fit, signed R^2 to maintain the scale-free topology. In this study, we obtained power 3 for this purpose. The second sub-figure (right) denotes soft threshold (power) (in x-axis) versus mean connectivity (in y-axis) plot. b. Gene dendrogram and module colours for NSE versus SE obtained by using dynamic tree cut method. The entire tree is cut at default height (0.99). Each colour signifies an individual resultant module

Figure 3. Co-expression module (‘greenyellow’ module, PCC: 0.6). Each node denotes an individual biomolecule where each edge signifies that the connecting node-pair has significant PCC (>0.6 or <-0.6). Insignificant PCC ($\le$0.6 or $\ge$-0.6) between any node-pair signifies that there is no connecting edge between the node-pair

Figure 4. Silhouette Width plot for NSE versus SE. The x-axis denotes Silhouette Width ([0, 1]) and y-axis denotes each individual cluster. The average Silhouette Width is obtained as 0.01

Subtype-specific regulatory network analysis

We evaluated the TF-target gene pairs and miRNA-target gene pairs with PCC cut-off 0.3, and p-value <0.05. Of note, the TFs, miRNAs and target genes in our regulation pairs were significant genes and miRNAs (Figure 5). As TFs may activate or suppress their target genes while miRNAs suppress their target genes, we retained TF-target gene pairs with a |PCC| greater than 0.3, while miRNA-target gene pairs with a PCC less than −0.3. The number of regulations is shown in Table 3. For NSE, there were 72 TF-target gene pairs, including 25 unique TFs and 61 unique target genes and 97 miRNA-target gene pairs, including 7 unique miRNAs and 63 unique target genes, respectively. For SE, there were 64 TF-target gene pairs including 24 unique TFs and 55 unique target genes, and 25 miRNA-target gene pairs including 6 unique miRNAs and 25 unique target genes. We constructed subtype-specific regulatory networks for NSE and SE, respectively.

Table 3. Statistics of regulation pairs in TGCT NSE and SE subtypes (PCC > 0.6, p-value <0.05)

TGCT subtype	Regulation	Relation	# pairs	# regulators	# target genes
NSE	TF-target	Positive	46	19	40
		Negative	26	16	26
		Total	72	25	61
	miRNA-target	Negative	97	7	63
SE	TF-target	Positive	42	19	39
		Negative	22	13	19
		Total	64	24	55
	miRNA-target	Negative	25	6	25

Figure 5. Subtype-specific regulatory networks related to Thyroid hormone signalling pathway. The molecules with red borders were related to the pathway. (a) The mentioned NSE-specific network. (b) The mentioned SE-specific network

For the NSE subtype-specific regulatory network, 142 nodes including 25 TFs, 7 miRNAs, and 120 target genes along with 169 edges were found. Among them, four were hub miRNAs and another four were hub genes. Among these four hub genes, two were TFs and the remaining two were target genes. The four top miRNAs, miR-520b, miR-520 c-3p, 182–5p, and miR-7-5p, were onco-miRNAs for various kinds of cancer including hepatocellular carcinoma, breast cancer, glioblastoma (GBM) and many more [21–23]. They all had higher out-degree, but no in-degree. The out-degree of them was 32, 28, 21, and 11, respectively. The top two TFs, SP1, and MYC, were both oncogenes. They had high out-degree (16 of SP1 and 12 of MYC), and low in-degree (1 of both genes). The top two target genes, TERT and RECK (in-degree 5 for both), were oncogenes.

For the SE-specific regulatory network, 103 nodes including 24 TFs, 6 miRNAs and 73 target genes along with 89 links were identified. Among these molecules, two miRNAs and five TFs were detected as hub. Also three target genes were selected as hubs from 73 target genes. The top two miRNAs, miR-182-5p and miR-338-3p, were onco-miRNAs for different kinds of cancer, including renal cell carcinoma, glioblastoma (GBM), osteosarcoma cells, gastric cancer, and many more. They had higher out-degree (16 and 4, respectively), but no in-degree. Hence, these two miRNAs made more sense than the other molecules. The top three target genes, CDH1, CXCR4, and SNAI1 had the higher in-degree (4, 3, and 3), but no out-degree. Among them, CDH1 and CXCR4 were oncogenes. Three TFs, SP1, HNF1A, and ETS1 had higher out-degree (16, 8, and 5, respectively), but no in-degree. Hence, these three TFs are extremely important.

Pathway analysis

KEGG pathway enrichment analysis was conducted by an online tool Database for Annotation, Visualization and Integrated Discovery (DAVID) [24]. There were 27 and 8 significant KEGG pathways (p < 0.05) enriched in NSE and SE, respectively. Interestingly, we found five enriched pathways that were shared by these two subtypes. They are Pathways in cancer (p = 5.30E-07 for NSE and p = 2.38E-05 for SE), Adherens junction (p = 1.84E-05 for NSE and p = 2.34E-02 for SE), Transcriptional misregulation in cancer (p = 3.60E-03 for NSE and p = 3.20E-03 for SE), Thyroid hormone signalling pathway (p = 4.60E-02 for NSE and p = 1.69E-02 for SE), and Thyroid cancer (p = 4.59E-02 for NSE and p = 2.59E-02 for SE).

Interestingly, with BH correction of p-values, the number of resultant entities for differential analysis or pathway analysis became very low. For differential analysis, we obtained very small number of DEM features (especially, number of DEM miRNAs). Similarly, for pathway analysis, we identified very small number of significant genes/miRNAs (i.e., only 9 BH corrected KEGG pathways for NSE analysis, and only one BH corrected KEGG pathway for SE analysis according to Table 5). Therefore, we considered only p-value without any correction to perform the entire analysis.

Table 5. KEGG pathway enrichment analysis of subtype-specific regulatory network genes using DAVID

Subtype	Pathway name	# Genes	Fold enrichment	P-value	Adjusted p (Bonferroni)	Adjusted p (BH)*
NSE	hsa05200:Pathways in cancer	19	4.0	5.30E-07	7.63E-05	7.63E-05
	hsa05205:Proteoglycans in cancer	13	5.4	3.48E-06	0.0005	0.0003
	hsa04520:Adherens junction	8	9.4	1.84E-05	0.0027	0.0009
	hsa04110:Cell cycle	9	6.0	0.0001	0.0144	0.0036
	hsa05206:MicroRNAs in cancer	13	3.8	0.0001	0.0171	0.0034
	hsa04550:Signalling pathways regulating pluripotency of stem cells	8	4.8	0.0013	0.1699	0.0306
	hsa04510:Focal adhesion	9	3.6	0.0029	0.3442	0.0585
	hsa04350:TGF-beta signalling pathway	6	5.9	0.0031	0.3594	0.0541
	hsa05213:Endometrial cancer	5	8.0	0.0032	0.373	0.0505
	hsa05222:Small cell lung cancer	6	5.9	0.0033	0.3743	0.0458
	hsa05202:Transcriptional misregulation in cancer	8	4.0	0.0036	0.4088	0.0467
	hsa05215:Prostate cancer	6	5.7	0.0038	0.4201	0.0444
	hsa04066:HIF-1 signalling pathway	6	5.1	0.0060	0.5779	0.0642
	hsa05161:Hepatitis B	7	4.0	0.0073	0.6511	0.0725
	hsa04151:PI3 K-Akt signalling pathway	11	2.7	0.0073	0.6528	0.0681
	hsa04390:Hippo signalling pathway	7	3.9	0.0088	0.7209	0.0767
	hsa04015:Rap1 signalling pathway	8	3.2	0.0120	0.8238	0.0971
	hsa04670:Leukocyte transendothelial migration	6	4.2	0.0128	0.8437	0.098
	hsa05130:Pathogenic Escherichia coli infection	4	6.5	0.0223	0.9611	0.1571
	hsa05016:Huntington’s disease	7	3.0	0.0260	0.9776	0.1729
	hsa04010:MAPK signalling pathway	8	2.6	0.0311	0.9894	0.1948
	hsa05203:Viral carcinogenesis	7	2.8	0.0344	0.9935	0.2046
	hsa05210:Colorectal cancer	4	5.4	0.0369	0.9956	0.2098
	hsa05230:Central carbon metabolism in cancer	4	5.2	0.0400	0.9972	0.2172
	hsa05216:Thyroid cancer	3	8.6	0.0459	0.9989	0.2372
	hsa04919:Thyroid hormone signalling pathway	5	3.7	0.0460	0.9989	0.2298
	hsa04022:cGMP-PKG signalling pathway	6	3.0	0.0468	0.999	0.2255
SE	hsa05200:Pathways in cancer	14	4.0	2.38E-05	0.0037	0.0037
	hsa05202:Transcriptional misregulation in cancer	7	4.7	0.0032	0.3888	0.2182
	hsa04919:Thyroid hormone signalling pathway	5	5.0	0.0169	0.9272	0.5824
	hsa04950:Maturity onset diabetes of the young	3	13.1	0.0211	0.9624	0.5597
	hsa04520:Adherens junction	4	6.4	0.0234	0.9741	0.5183
	hsa05216:Thyroid cancer	3	11.7	0.0259	0.9824	0.4901
	hsa04514:Cell adhesion molecules (CAMs)	5	4.0	0.0343	0.9954	0.5365
	hsa04060:Cytokine-cytokine receptor interaction	6	3.0	0.0486	0.9995	0.6165

* p-values were corrected by multiple tests using Benjamini-Hochberg Procedure [25].

Discussion

We constructed subtype-specific regulatory networks for NSE and SE subtypes. There were 142 nodes (25 TFs, 7 miRNAs, and 120 target genes) and 169 links in the NSE-specific regulatory network, while there were 103 nodes (24 TFs, 6 miRNAs, and 73 target genes) and 89 links in the SE-specific regulatory network. We identified four hub miRNAs and four hub genes (two TFs and two target genes) in the NSE-specific regulatory network, and two hub miRNAs and 8 hub genes (five TFs and three target genes) in the SE-specific regulatory network. One miRNA (miR-182-5p) and two genes (SP1 and MYC) were hubs for both of the two subtypes. We explored these hubs and their regulations in our analysis.

In the NSE network, the top four miRNAs (miR-520b, miR-520 c-3p, 182–5p, and miR-7-5p) were onco-miRNAs of many kinds of cancer, including breast cancer, glioblastoma (GBM), hepatocellular carcinoma, and so on [21–23]. They had maximum out-degree (32, 28, 21, and 11, respectively) and no in-degree (Table 4). We sorted the miRNA-target gene pairs according to their absolute PCC value. The top five target genes of miR-520b were RPS6KA1, PON2, ASH1 L, DENND5B, and OPCML. The top five target genes of miR-520 c-3p were RPS6KA1, PON2, DENND5B, ASH1 L, and OPCML. The top five target genes of miR-185-5p were NF1, RECK, CADM1, THBS2, and VCL. The top five target genes of miR-7-5p were BTN2A2, PRDM6, RNF144A, ST3GAL5, and CTSK.

Table 4. Hubs in subtype-specific regulatory networks with respect to in-degree and out-degree centralities

TGCT subtype	Hub type	Gene or miRNA hubs	Out-degree	In-degree
NSE	Out-degree based hubs	miR-520b	32	0
		miR-520 c-3p	28	0
		miR-182-5p	21	0
		SP1	16	1
		MYC	12	1
		miR-7-5p	11	0
	In-degree based hubs	TERT	0	5
		RECK	0	5
SE	Out-degree based hubs	SP1	16	0
		miR-182-5p	16	0
		HNF1A	8	0
		HIF1A	7	2
		MYC	5	2
		ETS1	5	0
		miR-338-3p	4	0
	In-degree based hubs	CDH1	0	4
		CXCR4	0	3
		SNAI1	0	3

In the NSE network, the top two TFs (SP1 and MYC) were oncogenes and had maximum out-degree (16 and 12, respectively) and in-degree of 1 (Table 4). Of the top five regulations for SP1, three target genes (TMPO, ATP2A2, and FGFR1) were activated, and two (OPN1LW and MC4 R) were suppressed. SP1 was activated by another TF ETS1. Of the top five regulations, MYC activated four target genes (TERT, CDC25A, TOP1MT, and PTBP1), while it only regressed one target gene (LDHA). MYC was suppressed by another TF HDAC2. SP1 (Sp1 Transcription Factor) is a zinc finger transcription factor that binds to GC-rich motifs of many promoters. The encoded protein is involved in many cellular processes, including cell differentiation, cell growth, apoptosis, and immune responses, responding to DNA damage, and chromatin remodelling (must provide scientific evidence). SP1 plays an important role during the tumour development, which can activate or suppress the expression of many oncogenes and tumour suppressor genes (TSGs) [26]. MYC (MYC Proto-Oncogene, BHLH Transcription Factor) is a proto-oncogene. Its encoded product plays a role in cell cycle progression, apoptosis, and cellular transformation. MYC is also a telomerase activity-related gene, as was shown in our previous study [8].

In the NSE network, the top two target genes (TERT and RECK) had maximum in-degree of, and both were oncogenes. Four TFs (MYCN, MYC, MAX, and HIF1A) activated TERT, while one TF (ETS1) suppressed TERT (Table 4). It was worth mentioning that MYCN was an important paralog of MYC, and MAX (MYC Associated Factor X) was also a TF that is related to MYC. Therefore, we believed that these three regulation pairs (MYC-TERT, MYCN-TERT, and MAX-TERT) played critical roles for NSE. RECK was activated by SP1, while regressed by three miRNAs (miR-520b, miR-520 c-3p, and miR-182-5p) and a TF (HIF1A). TERT (Telomerase Reverse Transcriptase) is a telomerase activity-related gene, and deregulation of telomerase expression in somatic cells may be involved in oncogenesis. We also found that TERT played important roles for TGCT in our previous study [8]. RECK (Reversion Inducing Cysteine Rich Protein with Kazal Motifs) is suppressed strongly in many tumours and cells transformed by various kinds of oncogenes, and its down-regulation by oncogenic signals may facilitate tumour invasion and metastasis.

In the SE network, the top two miRNAs (miR-182-5p, and miR-338-3p) were onco-miRNAs for many types of cancer, including glioblastoma (GBM), gastric cancer, renal cell carcinoma, osteosarcoma cells, and so on. They had maximum out-degree (16 and 4, respectively) and no in-degree (Table 4). The top five target genes repressed by miR-182-5p were RHOJ, CREB3L1, SLC1A1, OSBPL10, and CADM1. The four target genes repressed by miR-338-3p were SOX4, SLC22A23, HIF1A, and TOM1. We have discussed miR-182-5p, as it was also a hub for the NSE network. miR-338-3p was reported to be related to several types of cancer, such as GBM, cervical cancer, and gastric cancer.

In the SE network, the top five TFs (SP1, HNF1A, HIF1A, MYC, and ETS1) had the maximum out-degree (16, 8, 7, 5, and 5, respectively) and only two TFs (HIF1A and MYC) had in-degree of 2, while the other three TFs had no in-degree (Table 4). Among these five TFs, SP1 and MYC were also hubs in the NSE network, as we discussed above. In the top five regulations, SP1 activated two target genes (GLI1 and LAMC2), while suppressed three target genes (TNFSF13, CYP27A1, and LYZ). Among the top five regulations, HNF1A activated three target genes (CLDN2, S100G, and NEUROG3) as well as suppressed two target genes (SLC22A9 and SNAI1). In the top five regulations, HIF1A activated four target genes (VEGFB, NPPB, PGK1, and CDH1) while suppressed VHL. In the meanwhile, HIF1A was suppressed by a miRNA (miR-338-3p) and a TF (VHL). MYC activated five target genes (EBNA1BP2, TOP1MT, and LMNB2, CXCR4, and HDAC2). MYC was activated by E2F7, while suppressed by HDAC2. ETS1 activated four target genes (TBXAS1, CSF3 R, SPRY2, and THBD) and suppressed SNAI1. Because SP1 and MYC were also the hubs in the NSE network, we discussed the other three TFs. HNF1A (HNF1 Homeobox A) and HIF1A (Hypoxia Inducible Factor 1 Subunit Alpha) were oncogenes. ETS1 (ETS Proto-Oncogene 1, Transcription Factor) was also an oncogene.

In the SE network, the top three target genes (CDH1, CXCR4, and SNAI1) had the maximum in-degree (4, 3, and 3) (Table 4), among which CDH1 and CXCR4 were oncogenes. CDH1 was activated by HIF1A and GATA3, while it was suppressed by KDM5A and SNAI2. CXCR4 was activated by KLF5 and VHL, while it was suppressed by MYC. SNAI1 was suppressed by three TFs (HNF1A, ETS1, and HMGA2). CDH1 (Cadherin 1) and CXCR4 (C-X-C Motif Chemokine Receptor 4) were oncogenes. SNAI1 (Snail Family Transcriptional Repressor 1) was important mesoderm formation in the developing embryo.

We performed KEGG pathway enrichment analysis using the online tool DAVID [18]. We identified 27 and 8 significantly enriched pathways for NSE and SE, respectively (Table 5). Five pathways were enriched in both subtypes, including Pathways in cancer, Transcriptional misregulation in cancer, Adherens junction, Thyroid cancer, and Thyroid hormone signalling pathway. Interestingly, Thyroid hormone signalling pathway was enriched in both NSE and SE networks. It was reported that thyroid hormone was important for testicular development and function [27]. We therefore explored the subtype-specific regulatory networks related to thyroid hormone signalling pathway for NSE and SE, as shown in Figure 5 . Of note, ellipse, triangle, and diamond in the figure refer to TF, miRNA, and target gene, respectively. The molecules with read borders were related to the pathway under consideration. Four TFs (MYC, HIF1A, HDAC2, and CTNNB1) in this pathway regulated target genes, among which MYC and HIF1A regulated the expression levels of many target genes. This result suggested that these two TFs might have played critical roles in both NSE and SE.

In our proposed work, our main aim is to identify hub molecules (miRNAs, genes, and TFs) from each of the two kinds of networks, NSE-specific regulatory network and SE-specific regulatory network, by integrating TCGT expression and methylation data through the consecutive utilization of statistical, clustering, and network analysis. As a result, we identified four hub miRNAs (miR-182-5p, miR-520b, miR-520c-3p, and miR-7-5p), two hub TFs (MYC and SP1), and two genes (RECK and TERT) in the NSE-specific regulatory network, while for the SE-specific regulatory network, we obtained two hub miRNAs (miR-182-5p and miR-338-3p), five hub TFs (ETS1, HIF1A, HNF1A, MYC, and SP1), and three hub genes (CDH1, CXCR4, and SNAI1). To assess the reliability of the results, we focused on the hubs, which are highly connected and unlikely false-positive results, and also the enriched pathways, as well as literature evidence. In pathway enrichment analysis, we focus on the set of genes from our analysis; therefore, the false signal from one molecule would be reduced. We also found the hubs (miRNAs and genes) have been extensively supported by literature evidence. For example, those top four miRNAs (miR-520b, miR-520c-3p, 182–5p, and miR-7-5p) were reported as onco-miRNAs in many cancer types, such as breast cancer, glioblastoma, and hepatocellular carcinoma [21–23]. Furthermore, the top two target genes (TERT and RECK) were oncogenes and TERT was reported to play a significant role in TGCT [8]. In pathway analysis, we first used the normal p-value for each enriched pathway. When we further adjusted the p-value by Benjamini-Hochberg procedure, there were 9 KEGG pathways for NSE and only one KEGG pathway for SE that remained significant. Finally, we can validate these findings by using simulation data and also experiments. For simulation, we can randomly select a similar number of genes, miRNAs, and TFs as we found in NSE and SE analysis, and then perform the network analysis. This procedure will repeat many times (e.g., 100 or 1000) to obtain an empirical p-value. However, such a procedure is complicated as selection of random data should be in caution [28,29]. For experimental validation, we can use gene-editing technology to knock in or knock out of the critical hub molecules in cell lines or animal models and then perform genome-wide gene expression and methylation to examine the altered change.

Conclusions

TGCTs were a well-known malignant disease. A lot of studies had been conducted regarding TCGT. In this study, we performed an extensive analysis to extract the functional co-methylation and co-expression gene/miRNA modules for the TCGT subtypes, NSE and SE, through the integration of the gene/miRNA expression and methylation profiles. First of all, DNA methylation profile and mRNA/miRNA expression data were integrated using a latest score INTcon, and then applied CoMEx score along with the corresponding p-value to identify DEM genes/miRNAs. After that, we generated functional co-methylation and co-expression modules of those molecules for the integrated profile of these DEM genes/miRNAs through WGCNA method. The topmost module with the maximum average PCC comprised of 91 molecules. The cluster validity index measures were also provided to evaluate the quality of the clustering used in WGCNA. In addition, we chose the TF-target gene pairs and miRNA-target gene pairs with p-value cut-off of 0.05 and PCC cut-off of 0.6. We built NSE subtype-based regulatory network as well as SE subtype-based network. We obtained four hub miRNAs and four hub genes for NSE-specific regulatory network. These four hub miRNAs were miR-520b, miR-520c-3p, miR-182-5p, and miR-7-5p that were basically onco-miRNAs for various kinds of cancer including hepatocellular carcinoma, breast cancer, GBM, etc. They all had higher out-degree (32, 28, 21, and 11, respectively), but zero in-degree. Among these four hub genes, two were TFs, SP1, and MYC that were also oncogenes, whereas remaining two, TERT and RECK, were target genes. The out-degrees of these two TFs were 16 and 12, respectively. The in-degree of these two target genes were 5 for both, while out-degree of them was zero for both. These two target genes were also oncogenes. On the other hand, in case of the SE-specific network, we identified two hub miRNAs and eight hub genes. Those miRNAs were miR-182-5p and miR-338-3p that were also onco-miRNAs of different kinds of cancer. They had higher out-degree (16 and 4, respectively), but zero in-degree. Thus, these two miRNAs made more sense than the other molecules. Among those hub genes, five (SP1, HNF1A, HIF1A, MYC, and ETS1) were TFs and three (CDH1, CXCR4, and SNAI1) were target genes. Those three target genes had the higher in-degree (4, 3, and 3, respectively), but zero out-degree. Among them, CDH1 and CXCR4 were oncogenes. Among these five TFs, three (SP1, HNF1A, and ETS1) had higher out-degree (16, 8, and 5, respectively), but without any in-degree. Hence, these three TFs might be interesting. Interestingly, a miRNA (miR-182-5p) and two genes (SP1 and MYC) were common hubs obtained from both these subtypes.

Furthermore, we conducted KEGG pathway analysis, and five pathways (Pathways in cancer, Transcriptional misregulation in cancer, Adherens junction, Thyroid cancer, and Thyroid hormone signalling pathway) were found to be significantly enriched in both these subtypes. In conclusion, this work indicates that our method is highly useful to extract potential gene/miRNA modules as well as hubs for any kind of disease. In the future, we will extend it by integrating more multi-omics profiles to identify hub molecules and evaluate the true positives by using simulation data as well as by experiments.

Materials and methods

Data retrieval and pre-process

Clinical information

We downloaded the TGCT patient clinical pathological information from TCGA (https://xenabrowser.net/datapages/? cohort = TCGA%20Testicular%20Cancer%20(TGCT), accessed date: 20 October 2017). There were a total of 156 samples in the original clinical data file. We filtered the data with the following two conditions as indicated in our previous study [8]: (i) range of the age of the patient between 18 and 45; and (ii) the samples that were histologically defined as NSE or SE based on TCGA clinical data. This resulted in the determination of 48 NSE samples and 55 SE samples, respectively.

Expression profiles

Both of the gene expression and miRNA expression profiles of the patient samples presented above were downloaded from TCGA. For miRNA expression profile, we removed the miRNAs with missing value of more than 10% of the samples. We normalized the gene and miRNA expression levels using standard normalization technique.

Methylation profile

The DNA methylation profile of the patient samples was downloaded from TCGA. DNA methylation profile was measured experimentally using the Illumina Infinium HumanMethylation450 platform, and beta value was calculated for each array probe in each sample via BeadStudio software. We removed the probes with NA values. Furthermore, for the multiple to multiple mapping of probes and genes, we calculated the DNA methylation values for gene level using Principal Component Analysis (PCA). Finally, the DNA methylation values were normalized with z-score.

Transcriptional regulations

We extracted TF-gene pairs from TRANSFAC database [30] (released on 6 April 2016) which provides a comprehensive TF-target relation. There were 12,906 TF-target gene pairs having a total of 901 unique TFs and a total of unique 3,354 unique target genes. Of note, a subset of the target genes was denoted as TFs.

Post-transcriptional regulations

We selected three reliable miRNA-target prediction databases, TargetScan [31] (release 7.1, June 2016), miRanda [32] (released on August 2010), and PITA [33] (released on 9 December 2010). Furthermore, we extracted miRNA-target pairs from miRTarBase [34] (released version 7.0 on 15 September 2017) where the data were curated from low- and high-throughput experimental procedures. Moreover, we retained the pairs if they were present in at least two databases, which resulted in the identification of 170,544 miRNA-target pairs having a total of 697 unique miRNAs and a total of unique 12,507 target genes. Of note, a subset of the target genes was denoted as TFs.

Identification of differentially expressed and methylated (DEM) genes/miRNAs

We extracted differentially expressed and methylated (DEM) genes/miRNAs according to the approach in a previous study [13]. Since there were no tumour-matched normal samples, we defined DEM genes/miRNAs as those with significantly different expression levels and methylation levels between NSE group and SE group. In brief, we firstly normalized the gene/miRNA expression values and DNA methylation values with zero-mean normalization. Then, we integrated these two matrices into a single matrix using INTcon measure [13].

CoMEx score and corresponding p-value [13] were calculated from the integrated singular matrix for each gene/miRNA. Those genes/miRNAs having p-value less than 0.05 were defined as DEM genes/miRNAs. We discovered functional co-methylation and co-expression modules of molecules for these DEM genes/miRNAs based on their combined scores (INTcon values) using WGCNA (Weighted correlation network analysis). The INTcon score and CoMEx score are briefly described below.

Given that up-regulated, down-regulated, hyper-methylated and hypo-methylated genes are termed as DE_up, DE_down, DM_hyper, and DM_hypo, respectively, the combinatorial patterns of these differential conditions in between the expression and methylation scores are (a) (DE_up & DM_hypo), (b) (DE_down & DM_hyper), (c) (DE_up & DM_hyper), and (d) (DE_down & DM_hypo). As high methylation tends to have low gene expression, while low methylation activates gene expression level, we considered only the inversely correlated patterns [(a) and (b)]. The gene expression profile is initially normalized using zero-mean normalization. Similarly, the methylation data is then normalized. Then, INTcon score between two features is computed through the following way: INTcon(i,j) = const*0.5*|NI_exp(i,j)- NI_methyl(i,j)|, where ‘|.|’ denotes the absolute score, ‘NI_exp’ and ‘NI_methyl’ are the normalized expression and normalized methylation scores, respectively. Whenever a gene falls into the pattern (DE_up & DM_hypo) or (DE_down & DM_hyper), ‘const’ is ‘1’ else ‘0.’

CoMEx_i score of a gene i signifies the ratio of the difference in average value of group 1 and group 2 (i.e., avg(INTcon₁)-avg(INTcon₂)) to the sum of the standard error of the means of these two groups (termed as ‘stde_i’). A positive constant value (denoted as ‘ꞓ’) is for continuous dataset, in which each data value of the gene towards a sample signifies the integrated score of its (normalized) methylation and (normalized) gene expression scores estimated depending upon the inverse association between them: CoMEx_i = (avg(INTcon₁)-avg(INTcon₂))/(stde_i + ε). This score is equivalent to t-statistics and inspired by Significant Analysis of Microarrays (SAM) [35,36].

Gene module detection

After extracting the DEM genes and DEM miRNAs, we applied well-known R tool, WGCNA [16] to determine gene/miRNA modules. We also presented the soft threshold power selection plot, and dendrogram plot using dynamic tree cut method. In addition, we evaluated various cluster validity index measures such as Dunn Index, Silhouette Width, Clustering Co-efficient, etc., for the corresponding clustering.

Regulatory property analysis

In addition, we analysed the regulatory properties of these significant genes and miRNAs. We evaluated the transcriptional regulation pairs, i.e., TF-target gene pairs, and post-transcriptional pairs, i.e., miRNA-target gene pairs, for these significant genes and miRNAs with PCC (Pearson Correlation Coefficient). We constructed subtype-specific regulatory networks for NSE and SE, and identified hub nodes based on their in-degree and out-degree, by extending the hub concept introduced by Yu et al. [37]. Finally, the workflow of the proposed pipeline is represented in Figure 1.

Acknowledgments

We thank the lab members of Bioinformatics and Systems Medicine Laboratory (BSML) for their valuable discussion on this project.

Disclosure statement

The authors report no conflict of interests.

Additional information

Funding

This work was partially supported by the National Institutes of Health grant (R01LM012806). We thank the technical support from the Cancer Genomics Core funded by the Cancer Prevention and Research Institute of Texas (CPRIT RP170668 and RP180734). The funder had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.