Tamoxifen is an estrogen receptor (ER) antagonist that is most commonly used for the treatment of ER-positive breast cancer. However, tamoxifen resistance remains a major cause of cancer recurrence and progression. Here, we aimed to identify hub genes implicated in the progression and prognosis of ER-positive breast cancer following tamoxifen treatment. Microarray data (GSE9893) for 155 tamoxifen-treated primary ER-positive breast cancer samples were obtained from the Gene Expression Omnibus database. In total, 1706 differentially expressed genes (DEGs), including 859 up-regulated and 847 down-regulated genes, were identified between relapse and relapse-free samples. Weighted correlation network analysis clustered genes into 13 modules, among which the tan and blue modules were the most significantly related to prognosis. From these two modules, we further identified and validated two prognosis-related hub genes (G-rich RNA sequence binding factor 1 (GRSF1) and microtubule-associated protein τ (MAPT)) via survival analysis based on several publicly available datasets. High expression of GRSF1 predicted poor prognosis, whereas MAPT indicated favorable outcomes in ER-positive breast cancer. Using breast cancer cell lines and tissue samples, we confirmed that GRSF1 was significantly up-regulated and MAPT was down-regulated in the tamoxifen-resistant group compared with the tamoxifen-sensitive group. The prognostic value of GRSF1 and MAPT was also verified in 48 tamoxifen-treated ER-positive breast cancer patients in our hospital. Gene set enrichment analysis (GSEA) suggested that GRSF1 was potentially involved in RNA degradation and cell cycle pathways, while MAPT was strongly linked to immune-related signaling pathways. Taken together, our findings established novel prognostic biomarkers to predict tamoxifen sensitivity, which may facilitate individualized management of breast cancer.
Breast cancer is a heterogeneous cancer, displaying a variety of molecular features, prognostic patterns, and therapeutic responses . Up to two-thirds of all cases express estrogen receptor (ER), and can be treated using hormone-based therapy. Tamoxifen is a first-generation selective ER modulator that competes with estradiol to bind to ERs, thereby antagonizing the effects of estrogen and inhibiting the growth and proliferation of tumor cells . The administration of tamoxifen greatly minimizes the risk of recurrence of ER-positive breast cancer, particularly for premenopausal women . Unfortunately, approximately 40% of ER-positive patients are less sensitive to tamoxifen treatment, and will eventually relapse with endocrine-resistant phenotypes [4,5]. To date, the exact mechanisms of tamoxifen insensitivity in breast cancer remain largely unknown, and tamoxifen-resistant cancer is difficult to treat, due to lack of therapeutic targets. Since tamoxifen therapy fails for a large number of patients, there is an urgent need to elucidate the molecular mechanisms of tamoxifen resistance, particularly to identify novel potential genes for monitoring treatment efficacy and predicting prognosis.
Co-expression analysis has recently emerged as a powerful technique for mining gene expression profiles in various cancers. As an effective bioinformatics approach, weighted gene co-expression network analysis (WGCNA) is increasingly applied to explore synergistically altered gene sets, and to identify candidate biomarkers associated with clinical parameters [6–8]. In breast cancer, several studies have utilized WGCNA to identify hub genes closely related to clinicopathological traits (e.g., tumor size, grade, and molecular subtypes) and survival outcomes. For example, Tang et al. found that elevated expression of ASPM, TTK, and CDC20 conferred a poorer prognosis in breast cancer , and Jiang et al. identified six key genes (CA12, MLPH, FOXA1, GATA3, XBP1, and MAGED2) that could serve as biomarkers for the prediction of better chemotherapeutic responses and favorable prognosis in patients with breast cancer .
Accordingly, in the present study, we conducted an integrated analysis based on WGCNA to screen out novel prognostic biomarkers associated with tamoxifen response in breast cancer patients. In addition, the expression levels and the prognostic value of candidate hub genes were determined in vitro using cell lines and clinical tissue samples. Our findings may shed light on the underlying mechanisms of tamoxifen resistance in breast cancer, and may provide new prognostic markers to accurately predict tamoxifen response.
Materials and methods
Data collection and processing
The gene expression profile GSE9893 was obtained from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/), and evaluated using the GPL5049 platform . The dataset GSE9893 comprised 155 tamoxifen-treated primary breast cancer samples, of which 52 cases developed recurrent disease (designated the tamoxifen-resistant group). Robust multiarray average background correction and log2 conversion were performed using the ‘affy’ R package. Probes were mapped on to genes using Affymetrix annotation files. Genes matching with multiple probes were averaged to obtain the expression level of the gene. Probes corresponding to multiple genes were deleted.
Analysis of differentially expressed genes
The ‘limma’ R package with the Empirical Bayes method was employed to identify differentially expressed genes (DEGs) between relapse and relapse-free samples. Statistically significant DEGs were defined as |log2 FC| > 1 and P<0.01. The results were visualized by plotting a volcano plot using the ‘ggplot2’ package in R.
Functional enrichment analysis
After identifying DEGs related to tamoxifen sensitivity in breast cancer, the STRING database (https://string-db.org) was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis to determine the biological functions and pathways of tamoxifen resistance-related genes. The cutoff was set at an adjusted P-value of less than 0.05.
Co-expression network construction by WGCNA
Co-expression networks were established to explore modules involved in tamoxifen sensitivity in breast cancer using the ‘WGCNA’ package in R. First, outlier samples were detected using the sample network method. The soft threshold for WGCNA construction was selected such that the constructed network mainly included genes with strong correlations. We then transformed adjacency to a topological overlap matrix (TOM) to examine the connectivity of the network, followed by hierarchical clustering construction based on the TOM dissimilarity, to categorize genes with similar expression profiles into modules. The minimum module size for the gene dendrogram was 50, and other parameters were set to the default values. Finally, analyses of module eigengene, gene significance, and module–trait relationships were performed to identify clinically significant modules.
Selection of hub genes
Hub genes were identified as highly interconnected genes in a module of WGCNA. Tan and blue modules were considered key modules because they were closely related to the metastasis and recurrence of tamoxifen-resistant breast cancer. Hub genes were then screened out according to the absolute value of the Pearson’s correlation co-efficient. The modules of significance were visualized using Cytoscape (version 3.6.0; https://cytoscape.org/). The Cytoscape plugin ‘molecular complex detection’ (MCODE) was applied to detect the most important subnetworks, with a degree cutoff = 2, node score cutoff = 0.2, k-core = 2, and max depth = 100 set as the criteria .
Validation of prognosis-related hub genes by survival analysis
Survival analysis was performed with hub genes to further identify prognosis-associated genes using The Cancer Genome Atlas (TCGA) breast cancer dataset. All breast cancer patients were classified into two groups according to the expression level of a particular gene (high versus low). Kaplan–Meier survival analysis was then performed to compare the overall survival between these groups using the ‘Survival’ package in R. We further validated the survival results associated with each of the candidate hub genes in three independent ER-positive breast cancer cohorts, comprising the datasets GSE3494 and GSE25066, containing 201 and 296 ER-positive breast cancer patients, respectively, as well as the GSE9893 dataset. Results with P-values <0.05 were considered statistically significant.
The breast cancer cell line MCF-7 was purchased from the American Type Culture Collection (Manassas, VA, U.S.A.) and routinely maintained in Dulbecco’s modified Eagle’s medium (DMEM) with 10% fetal bovine serum. Tamoxifen-resistant cells (MCF-7/TAM) were established by continuously culturing MCF-7 cells in the presence of 4 μM 4-hydroxy-tamoxifen (Sigma–Aldrich, Missouri, U.S.A.) for 6 months. All cells were grown at 37°C in an atmosphere containing 5% CO2.
Reverse transcription quantitative polymerase chain reaction
RNA extraction was performed using TRIzol reagent (TaKaRa, Otsu, Japan). Total RNA was reverse transcribed to cDNA using a PrimeScript RT reagent kit (TaKaRa). Subsequently, qPCR was performed with harvested cDNA using the SYBR Green PCR kit (TaKaRa). The relative mRNA levels were calculated using the 2−ΔΔCt method taking GAPDH as the internal control. The primers used for reverse transcription quantitative polymerase chain reaction (RT-qPCR) were as follows: G-rich RNA sequence binding factor 1 (GRSF1) forward, 5′-ACAGGGAAGAAATTGGTAATCG-3′ and reverse, 5′-ACCATCGTCTACTGCCCTTTC-3′; and microtubule-associated protein τ (MAPT) forward, 5′-AAAGACGGGACTGGAAGCG-3′ and reverse, 5′-GAATCCTGGTGGCGTTGG-3′.
Preparation of ER-positive human breast cancer samples
Twenty-three tamoxifen-resistant and 25 tamoxifen-sensitive paraffin-embedded tumor samples were obtained from ER-positive breast cancer patients who underwent surgery at the Department of Breast Surgery, The Second Affiliated Hospital, Zhejiang University School of Medicine, during June 2012 to September 2019. Informed consent was obtained from all patients. The study was performed with the approval of The Human Research Ethics Committee of The Second Affiliated Hospital, Zhejiang University School of Medicine. The present study conformed to the Declaration of Helsinki.
Paraffin-embedded tissue sections were dewaxed with xylene and rehydrated with ethanol, followed by antigen retrieval with EDTA (pH 9.0). Endogenous peroxidase was removed by adding 3% H2O2. The slides were incubated with goat serum and anti-GRSF1 (Abcam, MA, U.S.A., dilution 1:100) or MAPT (Abcam, dilution 1:400) primary antibodies overnight at 4°C. Detection was performed by incubating with horseradish peroxidase (HRP)-linked anti-Rabbit IgG and 3,3′-diaminobenzidine (DAB). The staining intensity was scored as follows: 0, negative; 1, weak; 2, moderate; 3, strong. The percentage of stained cells was scored into four grades: 0, < 5%; 1, 5–25%; 2, 25–50%; 3, >51%. Intensity and percentage scores were multiplied to obtain the final scores (0, 1, 2, 3, 4, 6, or 9), with a cutoff <3 versus ≥3 to indicate low versus high expression, respectively. Staining was scored by two independent observers in our hospital.
Gene set enrichment analysis of hub genes
To explore the molecular mechanisms of identified hub genes on breast cancer, gene set enrichment analysis (GSEA) was carried out with the ER-positive TCGA dataset . The samples were separated into low and high groups in accordance with degree of hub gene expression, and c2.cp.kegg.v5.2.symbols.gmt was selected as a reference gene set. A false discovery rate < 0.05 was designated as the cut-off criteria.
Tumor immune estimation resource database analysis
Because of the essential role of immune infiltration in cancer initiation and progression, we used the tumor immune estimation resource (TIMER) online database to determine the association between tumor-infiltrating immune cells and each hub gene . The six types of immune cells inferred by TIMER included B cells, CD4 T cells, CD8 T cells, dendritic cells, macrophages, and neutrophils. The levels of hub gene expression were visualized by log2 RSEM.
All in vitro experiments were independently repeated three times. Two-tailed Student’s t tests were used to detect differences between groups using SPSS 17.0 software (IBM, NY, U.S.A.). In survival analysis, we used Kaplan–Meier analysis and a Cox proportional hazards regression model to ascertain whether candidate hub genes had an effect on prognosis. The hazard ratio (HR) and 95% confidence interval (CI) were calculated from the regression coefficients and survival curves were plotted using GraphPad Prism 6.0 (GraphPad Software, CA, USA). P-values <0.05 were considered to be significant.
Identification and functional annotation of DEGs
After applying thresholds of |log2 FC| > 1 and P<0.01, we identified 1706 DEGs, including 859 up-regulated and 847 down-regulated genes, between tamoxifen-sensitive and tamoxifen-resistant breast cancer samples. A volcano plot of the DEGs is shown in Supplementary Figure S1.
According to GO enrichment analysis, up-regulated genes were significantly enriched in various biological processes (BPs) including ‘protein targeting to the ER’. The down-regulated genes were primarily enriched in ‘signal release’ and ‘positive regulation of hormone secretion’ (Supplementary Figure S2). According to KEGG analysis of the up-regulated DEGs, ‘ribosome’ and ‘oxidative phosphorylation’ were the most obviously enriched keywords (Supplementary Figure S3).
Weighted co-expression network construction and key module identification
In the present study, 28 abnormal samples were excluded (Figure 1A). The value of soft-thresholding powers (β) = 6 was selected to achieve a relatively scale-free network, which was closer to the real biological network state (Figure 1B,C). We then identified 14 modules via average linkage hierarchical clustering. The DEGs in gray were not included in any module; therefore, we did not perform any functional analysis of the DEGs in gray (Figure 1D). Of these modules, the tan module showed obvious positive correlations with relapse, distant metastasis, and death. Moreover, a significant negative correlation was found between the blue module and poor prognosis (Figure 2). Hence, the tan and blue modules may play essential roles in the BPs of breast cancer tamoxifen resistance. Thus, these modules, as the most related to disease progression, were chosen for further analysis.
Clustering of samples and determination of soft-thresholding power in the WGCNA
Hub modules selection
Identification of hub genes in the tan and blue modules
Hub genes have high connectivity within clinic-related modules, and tend to play critical roles in the molecular mechanisms of tamoxifen resistance. Therefore, we next used Cytoscape to visualize hub gene networks in the tan and blue modules. As shown in Figure 3, 38 and 50 genes with the highest intramodular connectivity in the tan and blue modules, respectively, were screened.
The visualization of hub genes
Identification of prognosis-related hub genes
To further explore the effects of these candidate key genes on prognosis in breast cancer, we conducted survival analysis of 88 hub genes based on TCGA data. High expression of three hub genes, i.e. GRSF1, cytochrome c oxidase subunit 7B (COX7B), and chaperonin containing TCP1 subunit 8 (CCT8), in the tan module were all significantly associated with poor survival outcomes, whereas MAPT and REC8 meiotic recombination protein (REC8) in the blue module all predicted better prognosis in breast cancer when overexpressed (Figure 4). Thus, these five genes were chosen as candidates for further study.
Survival analysis of prognosis-related hub genes in breast cancer patients from TCGA dataset
Validation of prognosis-related hub genes in three ER-positive breast cancer cohorts
Subsequently, we validated the prognostic relevance of the five selected hub genes (GRSF1, COX7B, CCT8, MAPT, and REC8) based on three independent ER-positive breast cancer cohorts (GSE9893, GSE3494, and GSE25066 datasets) involving 652 patients. Two of the five genes, i.e., GRSF1 and MAPT, remained significantly associated with prognosis in patients with ER-positive breast cancer in these datasets (P<0.05; Table 1).
|Gene .||GSE9893 cohort .||GSE3494 ER+ cohort .||GSE25066 ER+ cohort .|
|.||HR .||95% CI .||P-value .||HR .||95% CI .||P-value .||HR .||95% CI .||P-value .|
|Gene .||GSE9893 cohort .||GSE3494 ER+ cohort .||GSE25066 ER+ cohort .|
|.||HR .||95% CI .||P-value .||HR .||95% CI .||P-value .||HR .||95% CI .||P-value .|
*P<0.05 and **P<0.01.
Validation of hub genes in vitro
Finally, we applied RT-qPCR to further verify the expression of candidate hub genes in breast cancer cell lines. The expression level of GRSF1 was higher in tamoxifen-resistant cells, whereas MAPT expression levels were significantly elevated in MCF-7 cells compared with parental MCF-7/TAM cells (Figure 5A). Similarly, the aberrant expression pattern of GRSF1 and MAPT was also verified by immunohistochemical (IHC) staining in 48 clinical tissue samples (23 tamoxifen-resistant versus 25 tamoxifen-sensitive) collected from breast cancer patients receiving tamoxifen treatment in our hospital (GRSF1: P=0.009; MAPT: P=0.017, Figure 5B). Survival analysis demonstrated that up-regulated GRSF1 expression was significantly associated with poorer disease-free survival (DFS) (HR = 2.348, 95% CI: 1.032–5.346, P=0.044, Figure 5C), and MAPT has the capacity to predict favorable DFS in these patients (HR = 0.435, 95% CI: 0.191–0.988, P=0.048, Figure 5D). The data of these independent experiments therefore verified the hypotheses generated using bioinformatics analysis, indicating that GRSF1 and MAPT might play a vital role in tamoxifen resistance in breast cancer.
Validation of prognosis-related hub genes in vitro
To better understand the underlying function of these hub genes, GSEA was carried out and mapped on to KEGG pathways. As illustrated in Figure 6, GRSF1 was mostly involved in ‘ubiquitin-mediated proteolysis’, ‘oocyte meiosis’, ‘RNA degradation’, ‘cell cycle’, and ‘mismatch repair’. MAPT was related to ‘antigen processing and presentation’, ‘natural killer cell-mediated cytotoxicity’, ‘autoimmune thyroid disease’, ‘cell adhesion molecules’, and the ‘proteasome’.
Association of hub gene expression and immune infiltration level
The distribution of tumor-infiltrating cells is highly relevant to tumor progression. Therefore, we evaluated the association of hub genes with immune infiltration level using the TIMER platform. The level of GRSF1 expression was positively correlated with the abundance of infiltrating immune cells. While MAPT expression displayed a significant negative correlation with infiltration degree by B cells, CD4 T cells, CD8 T cells, neutrophils, and dendritic cells (Figure 7). These findings suggest that GRSF1 and MAPT may be involved in immune infiltration in patients with breast cancer.
The correlation between hub genes and immune cell infiltration levels in breast cancer through TIMER
ER-positive breast cancer exhibits a favorable prognosis owing to the efficacy of anti-estrogen drugs, such as tamoxifen . However, one-third of these patients eventually develop tamoxifen resistance, resulting in cancer progression and death [5,15]. Tamoxifen resistance occurs via a complicated series of events, taking place over multiple genes and various signaling pathways. An in-depth elucidation of the biological mechanisms of tamoxifen insensitivity is beneficial to identify novel prognostic biomarkers, and explore effective therapeutic targets towards overcoming tamoxifen resistance. Due to the establishment of large cancer databases, such as the TCGA and GEO databases, researchers have the capacity to investigate large-scale gene expression profiles . In the present study, we screened for hub genes involved in the development of tamoxifen insensitivity that could be used as potential biomarkers to predict tamoxifen response and prognosis in ER-positive breast cancer patients.
In the present study, we first identified 1706 DEGs associated with tamoxifen resistance, including 859 up-regulated and 847 down-regulated genes. These DEGs were primarily enriched in functions such as protein targeting the ER and pathways such as oxidative phosphorylation. As a hormonal transcription factor, ERs regulate target genes to manipulate cell cycle progression and the endocrine response. The activity of ERs is also regulated by multiple proteins, including the transcription factors Ap-1 and FOXA1, which exert different biological functions in response to endocrine treatment [17–19]. Meanwhile, many studies have shown that oxidative phosphorylation is closely correlated with carcinogenesis. Echeverria et al. reported that an oxidative phosphorylation inhibitor delayed residual tumor regrowth for neoadjuvant chemotherapy-resistant patients with breast cancer . Sansone et al. demonstrated that the activation of oxidative phosphorylation promoted the development of hormone therapy-resistant disease . Overall, these studies imply that the DEGs identified in the present study might be closely connected with tumor progression and endocrine efficiency.
Subsequently, we utilized WGCNA to filter highly reliable and biologically significant modules and hub genes that are responsible for tamoxifen resistance from the list of DEGs [22,23]. The WGCNA clustered genes into 13 modules, of which the tan and blue modules were positively and negatively related to clinical traits, respectively. From this analysis, hub genes in these two modules were selected. Subsequent survival analysis showed that high expression of GRSF1 predicted poor prognosis, whereas MAPT was associated with favorable survival outcomes in TCGA breast cancer patients. Notably, we further verified the prognostic value of candidate hub genes in three independent ER-positive breast cancer cohorts. Compared with previous studies [15,25], we found modules and genes that were relevant to malignant phenotypes and favorable clinical features. More importantly, we validated the hypotheses generated by the available databases using tamoxifen-sensitive and tamoxifen-resistant cell lines as well as clinical tissue specimens. We thus identified GRSF1 and MAPT as the most promising candidate genes related to tamoxifen resistance. Subsequently, we further explored the potential roles of hub genes in ER-positive breast cancer using GSEA software. The results showed that several pathways such as ‘RNA degradation’ and ‘cell cycle’ were dysregulated when GRSF1 was aberrantly expressed. The potential mechanism of MAPT is strongly linked to immune-related signaling pathways such as ‘antigen processing and presentation’, ‘natural killer cell-mediated cytotoxicity’, and ‘cytokine–cytokine receptor interaction’. TIMER analysis also indicated a correlation between MAPT expression and immune infiltration, suggesting that MAPT might have a function in tumor immunity.
GRSF1 was initially identified as an RNA-binding protein with high affinity for G-rich sequences. GRSF1 plays critical roles in maintaining mitochondrial function, including mitochondrial translation, mitochondrial ribosome biosynthesis, and mitochondrial noncoding RNA binding [24,25]. At present, only a small number of studies have focused on the role of GRSF1 in cancer. Sun et al. and Yang et al. demonstrated that GRSF1 regulates miRNAs to facilitate oncogenic behaviors, including autophagy and metastasis, in cervical cancer [26,27]. Wang et al. revealed that GRSF1 can accelerate tumorigeneis and metastasis via PI3K/AKT pathway in gastric cancer . Taken together with the results of the present study, these findings imply that GRSF1 might function as a potential oncogene.
MAPT is a gene encoding τ protein, which is implicated in the pathogenesis of several neurodegenerative disorders such as Alzheimer’s disease, Parkinson’s disease, and progressive supranuclear palsy [29,30]. Recent studies have suggested that elevated expression of MAPT predicts better survival outcomes in pediatric neuroblastoma, breast cancer, renal clear cell cancer, and low-grade glioma, which is consistent with the results of the present study [31–34]. Wang et al. reported that MAPT-hypermethylated tumors are closely associated with poor prognosis in patients with colorectal cancer . Interestingly, MAPT is reported to play an essential role in mediating paclitaxel or taxane resistance in various cancers. Rouzier et al. first identified MAPT as a predictor of the response to paclitaxel in breast cancer . MAPT can also determine paclitaxel chemosensitivity by interacting with several miRNAs in gastric cancer and non-small cell lung cancer . Moreover, clinical and in vitro studies have demonstrated that the expression level of MAPT is positively associated with ER expression, and is influenced by ER signaling [38,39]. Taken together, these studies indicate that MAPT clearly plays a complex and possibly cancer-specific role in different cancers, which warrants more in-depth, well-designed investigations.
However, the present study had some limitations. First, tamoxifen resistance is controlled by a complicated regulatory network comprising mRNAs, miRNAs, and long noncoding RNAs; however, since we were restricted by the available datasets, only protein-coding genes/mRNAs were included in the present analysis. The precise roles of MAPT and GRSF1 may only become clear in the context of miRNAs and long noncoding RNAs. Second, despite the validation of key genes in cancer cell lines and tissue samples, our model has not been verified in a sufficiently large clinical cohort, or prospective individual cohorts. Third, we predicted the possible functions of specific genes using the available network information, but the underlying mechanisms of gene networks involved in tamoxifen response warrant further study.
The present study identified gene networks and potential prognostic biomarkers using a systems biology-based WGCNA approach in patients with primary breast cancer treated with tamoxifen. Through a series of bioinformatics analyses and preliminary biological experiments, we identified and verified two novel biomarkers that may be related to the tamoxifen response in ER-positive breast cancer: GRSF1, a prognostic marker for cancer progression, and MAPT, to predict favorable survival outcomes. GSEA suggested that GRSF1 might be involved in RNA degradation and cell cycle pathways, while MAPT was closely linked to immune-related signaling pathways. However, further studies are needed to elucidate the exact molecular mechanisms and characterize the key genes functionally affecting tamoxifen sensitivity in patients with breast cancer.
Publicly available datasets were analyzed in the present study, these can be obtained from GEO (GSE9893, GSE3494, GSE25066) and TCGA.
The authors declare that there are no competing interests associated with the manuscript.
This work was supported by the Natural Science Foundation of Zhejiang Province [grant numbers LQ19H060002, LQ19H160041].
CRediT Author Contribution
Yanyan Wang: Conceptualization, Software, Formal analysis, Validation, Methodology, Writing—original draft. Xiaonan Gong: Resources, Software, Validation, Visualization. Yujie Zhang: Conceptualization, Supervision, Validation, Writing—review & editing.
chaperonin containing TCP1 subunit 8
cytochrome c oxidase subunit 7B
differentially expressed gene
G-rich RNA sequence binding factor 1
gene set enrichment analysis
Kyoto Encyclopedia of Genes and Genomes
microtubule-associated protein τ
REC8 meiotic recombination protein
reverse transcription quantitative polymerase chain reaction
The Cancer Genome Atlas
tumor immune estimation resource
topological overlap matrix
weighted gene co-expression network analysis