A Machine Learning Approach to Differentiate Two Specific Breast Cancer Subtypes Using Androgen Receptor Pathway Genes.

Triple-negative breast cancer is a heterogeneous disease with different molecular and histological subtypes. The Androgen receptor is expressed in a portion of triple-negative breast cancer cases and the activation of the androgen receptor pathway is thought to be a molecular subtyping signature as well as a therapeutic target for triple-negative breast cancer. Thus, identification of the androgen receptor pathway status is important for both molecular characterization andclinical management. In this study, we investigate the expression of the androgen receptor pathway in metaplastic breast cancer and luminal androgen receptor subtypes of triple-negative breast cancer and found that the androgen receptor pathway was downregulated in metaplastic breast cancer compared to luminal androgen receptor subtype. Using random forest, we found that the two subtypes of breast cancer can be molecularly classified with the gene expression of the androgen receptor pathway.

Introduction

Breast cancer is a heterogeneous disease with different molecular features and prognoses. Among them, triple-negative breast cancer (TNBC) which is defined by the lack of expression in estrogen receptor (ER), progesterone receptor (PR) and human epidermal receptor 2 (HER2) by immunohistochemical staining has the most limited therapy choice and worst clinical outcome. TNBC can be further classified into subtypes according to histological morphology as well as molecular features. The histological subtypes of TNBC are composed of the commonest invasive ductal carcinoma of no special type (IDC-NST) and other special subtypes including metaplastic breast cancer, adenoid cystic carcinoma, medullary carcinoma and secretory carcinoma. Studies have shown that TNBC of special types as a single group has a worse prognosis than TNBC-NST, indicating the prognostic value of histological subtyping. Metaplastic breast cancer (MBC) was a special subtype of breast cancer accounting for less than 1% of all invasive breast cancer, characterized by the presence of metaplastic components in cancer tissue which is most commonly squamous carcinoma, followed by chondroid and sarcoma components. Most MBC were triple-negative, and study has shown that MBC has a worse prognosis in all clinical stages after treatment compared to other TNBC. Due to the limited cases of MBC, our understanding of their molecular characteristics remains largely unrevealed.

Molecularly, TNBC can also be classified into various subtypes by different algorithms using gene expression data. Though all of the currently applied subtyping algorithms could distinguish a consistent molecular subtype in TNBC which was the luminal androgen receptor (LAR) subtype. LAR accounted for 15%-20% of all TNBC and was characterized by the high expression of the AR gene and enrichment in hormonally regulated pathways. LAR subtype had a relatively low proliferation rate, decreased relapse-free survival and similar distant metastasis-free survival compared with other subtypes and can potentially benefit from anti-AR molecule enzalutamide. Since immunohistochemical stain for AR in TNBC showed that 38%-55% of TNBC has positive AR expression, using AR as a surrogate marker of LAR subtype would reveal low specificity.

Recent studies reported the percentage of AR-positive expression cases in MBC to be 0%, 8.7% and 11% respectively which was significantly lower than that in TNBC-NST, indicating the lack of luminal differentiation in MBC. Genomic mutation characterization of MBC revealed that it harbored a mutation rate of 57% in PI3K/AKT/mTOR pathway which was much higher than the 4% in AR-negative TNBC but closer to 40% in AR-positive TNBC. Thus, whether the low expression of AR in MBC also indicated the downregulation of the AR pathway and the exact molecular difference between the MBC and LAR group remains unknown.

In this study, we analyzed and compared the expression of AR pathway genes in MBC and LAR using data from TCGA. A machine learning approach was used to differentiate MBC and LAR with AR pathway genes.

Results

Clinicopathological Characteristics of the Studied Cohort

A total of 38 cases of LAR and 14 cases of MBC were selected in the TCGA database. The clinicopathological characteristics including age at diagnosis, ethnicity, tumor stage, tumor size and lymph node status were analyzed with no significant difference detected between the two groups (Table 1).

Table 1.

Clinical Features of Selected Patients.

Dependent		LAR	MBC	P value
Age	<50	8 (21.1)	2 (14.3)	0.879
	≥ 50	30 (78.9)	12 (85.7)
Ethnicity	Hispanic or Latino	2 (5.3)	1 (7.1)	0.546
	Not Hispanic or Latino	33 (86.8)	13 (92.9)
	Not reported	3 (7.9)	0 (0.0)
Tumor stage	Stage I	7 (18.4)	2 (21.4)	0.942
	Stage IIa	15 (39.5)	5 (35.7)
	Stage IIb	6 (15.8)	3 (21.4)
	Stage III	1 (2.6)	0 (0.0)
	Stage IIIa	4 (10.5)	2 (14.3)
	Stage IIIb	1 (2.6)	1 (7.1)
	Stage IIIc	3 (7.9)	0 (0.0)
	Not reported	1 (2.6)	0 (0.0)
Tumor	T1	5 (13.2)	2 (14.3)	0.200
	T2	17 (44.7)	2 (14.3)
	T3	1 (2.6)	0 (0.0)
	Tx	1 (2.6)	0 (0.0)
	Not reported	14 (36.8)	10 (71.4)
Lymph node	N0	9 (23.7)	4 (28.6)	0.082
	N1	9 (23.7)	0 (0.0)
	N2	3 (7.9)	0 (0.0)
	N3	3 (7.9)	0 (0.0)
	Not reported	14 (36.8)	10 (71.4)

Androgen Receptor Pathway Genes Were Differentially Expressed in MBC and LAR

A total of 166 genes were identified as the representative genes in the androgen receptor pathway using the Pathway Commons database (Version 12). In addition, recent research has identified another hormonal receptor gene G-protein coupled estrogen receptor (GPER), which was encoded by GPER1. GPER can be activated by hormonal estradiol. Unlike ERalpha and ERbeta which are mostly known to be nuclear receptors, GPER has a seven-transmembrane domain and many studies have confirmed its membrane localization. It was found to be expressed strongly in triple-negative breast cancer and patients younger than 49-years-old. The expression of GPER has reversely correlated with the expression of androgen receptor in TNBC and at the molecular level AR has a repressed regulation on GPER by binding to the promoter of AR genomic region. Thus, GPER1 was also included in our analysis as one of the AR pathway genes. The mRNA expression of genes in the AR pathway was analyzed and compared in the 2 groups. Differentially expressed genes were identified and summarized in Table 2. In total, 32 out of the 167 genes have been found to be differentially expressed between MBC and LAR, including RUNX2, AR and GPER1 (Figure 1). The top 5 genes with the highest significance were RUNX2, SPDEF, FOXA1, DDC and AR. Except for DDC which was a metabolic enzyme, the other 4 genes were all transcription factors that have previously been shown to act intimately with one another. Among them, RUNX2 was the only upregulated gene in MBC and it was reported to inhibit the effect of AR as a transcription factor by promoting the dissociation of AR from the targeted genes. The SPDEF was downstream of AR, whose expression was induced by AR. FOXA1 was the pioneer gene in the AR pathway and acted by loosening the AR-binding DNA region to facilitate the binding of AR.

Table 2.

Differentially Expressed AR Pathway Genes Between MBC and LAR Cancers.a

Name	Ensemble Id	log FC	Ave expr	t	P value	B
RUNX2	ENSG00000124813	1.69	15.63	5.84	3.51E-07	6.48
SPDEF	ENSG00000124664	−4.41	18.47	−5.06	5.67E-06	3.86
FOXA1	ENSG00000129514	−3.81	17.56	−4.35	6.44E-05	1.59
DDC	ENSG00000132437	−5.73	11.54	−4.30	7.47E-05	1.45
AR	ENSG00000169083	−2.79	16.30	−3.87	0.000308289	0.14
FKBP4	ENSG00000004478	−0.81	20.09	−3.61	0.00069105	−0.60
SLC25A4	ENSG00000151729	−0.70	16.84	−3.50	0.000974849	−0.92
ETV5	ENSG00000244405	1.25	16.41	3.26	0.001946062	−1.55
FLNA	ENSG00000196924	0.84	21.03	3.20	0.002348338	−1.72
SMAD3	ENSG00000166949	0.86	16.82	3.15	0.002696079	−1.85
SIRT1	ENSG00000096717	−0.70	17.01	−3.00	0.00413892	−2.23
RCHY1	ENSG00000163743	−0.49	16.04	−2.99	0.004228911	−2.25
TGIF1	ENSG00000177426	−0.51	17.64	−2.96	0.004681079	−2.34
TGFB1I1	ENSG00000140682	0.96	16.89	2.84	0.006513312	−2.64
NCOR2	ENSG00000196498	0.48	18.05	2.81	0.006993992	−2.70
NCOA4	ENSG00000266412	−0.47	19.99	−2.80	0.007114918	−2.72
HSP90AA1	ENSG00000080824	−0.61	22.56	−2.76	0.007991213	−2.82
SVIL	ENSG00000197321	0.69	17.68	2.72	0.008962078	−2.92
SF1	ENSG00000168066	0.23	19.67	2.69	0.009689098	−2.99
PRDX1	ENSG00000117450	−0.56	22.56	−2.62	0.011373878	−3.13
HDAC1	ENSG00000116478	0.39	19.47	2.54	0.014042583	−3.32
GPER1	ENSG00000164850	1.05	13.97	2.50	0.015470474	−3.40
GTF2H2	ENSG00000145736	0.88	12.40	2.46	0.017417516	−3.51
CASP8	ENSG00000064012	−0.56	16.41	−2.39	0.020523257	−3.65
CDC25B	ENSG00000101224	0.63	18.83	2.27	0.027069795	−3.89
KAT5	ENSG00000172977	−0.27	17.29	−2.13	0.038277859	−4.18
AHR	ENSG00000106546	0.67	18.12	2.11	0.039491891	−4.21
CDK1	ENSG00000170312	−0.75	18.05	−2.10	0.040386094	−4.23
CAV1	ENSG00000105974	0.75	18.77	2.07	0.043886593	−4.30
NR0B2	ENSG00000131910	−3.18	5.85	−2.03	0.047145316	−4.36
GTF2F2	ENSG00000188342	0.34	17.31	2.02	0.048112999	−4.37
FHL2	ENSG00000115641	0.75	17.56	2.02	0.048381632	−4.38

a The columns of the table are the gene name, the gene id, the estimated contrast, the expression mean over both groups, contrast t-value, contrast P-value and the estimated log-odds probability ratio (B) that the gene is differentially expressed.

Figure 1.

AR pathway genes were differentially expressed in MBC and LAR. AR was highly expressed in the LAR group while its expression in MBC was low (left panel). The membrane-bound estrogen receptor, GPER1 showed a higher expression in MBC than in LAR (middle panel). As the gene with most significant expression difference, RUNX2 was upregulated in MBC while downregulated in LAR (right panel).

Classification of MBC and LAR Using Random Forest

The above results suggested that MBC and LAR were differently regulated in the AR pathway. Next, we try to directly differentiate the two groups using gene expression data of the AR pathway. Whereas, using the expression data of a single gene was unable to classify the two groups at 100% efficacy as shown in Figure 1. The machine learning approach was reported to be able to achieve good predictive performance for sample classification using gene expression data. Thus, we further tried to look at the effect of androgen receptor pathway genes on classifying the MBC and LAR groups via the random forest algorithm. Random forest is an algorithm for classification developed in 2001 that uses an ensemble of classification trees and it was widely used in the classification using microarray data. In this task, the expression of the 167 AR pathway genes was used as continuous variables to classify the sample as either MBC or LAR (Figure 2A). The prediction accuracy using the random forest algorithm was 100% (Table 3). Genes that contributed to the classification most were listed in Figure 2B and C. The contribution was measured by Mean Decrease Accuracy or Mean Decrease Gini. RUNX2, FKBP4 and UXT were ranked as the top 3 genes by both Mean Decrease Accuracy or Mean Decrease Gini. Interestingly, the UXT gene was not listed in the DEGs between MBC and LAR, Model visualization was performed by displaying decision tree with the most and least nodes (Figure 3). In the simplest decision tree generated by the random forest algorithm which has three nodes, RUNX2 which has the most significant differential expression between MBC and LAR was used as the root node and no other internal node was used.

Figure 2.

Classifying MBC and LAR using random forest algorithm. Clustering of MBC samples (blue) and LAR sample (red) using 167 AR pathway genes (A). Genes that contributed most to the classification were listed using 2 different parameters (B and C).

Table 3.

Classification Accuracy of the Random Forest Model.

Actual classification	Predicted classification
	MBC	LAR
MBC	38	0
LAR	0	14
Prediction accuracy	100%

Figure 3.

Visualization of 2 representative trees with the maximum and minimum nodes generated by random forest. The tree with maximum nodes used SPDEF gene expression value as the root node and the expression of other 9 genes as internal nodes, making the total nodes number to be 21. It was a 2-class split for each root and internal node which was determined by the gene expression value of the specific gene in the node. The cutoff value for the binary split in each node was calculated automatically (A). The tree with minimum nodes used the expression of the RUNX2 gene as the single root and internal node, generating 2 leaf nodes.

In the model construction, a 5-fold cross-validation was also performed for 100 times to avoid overfitting. Average cross-validation error and standard deviation were plotted in Figure 4. It was found that when the number of variables was in the range of 5 to 21, the error of cross-validation reached the minimum value.

Figure 4.

Cross-validation of the random forest algorithm for classification of MBC and LAR. A 5-fold cross-validation was performed for 100 times with the number of variables ranging from 1 to 166. The average value and standard deviation for cross-validation were plotted.

Discussion

AR was expressed in a proportion of TNBC and the activation of AR was thought to be a signature for the LAR subtype of TNBC which can be used as a therapeutic target. Thus, identification of the AR pathway status in TNBC cases was important for both molecular characterization and clinical management. In this study, we showed that the AR pathway was differently regulated in MBC and LAR of TNBC. Moreover, through the random forest, the 2 groups of TNBC can be classified using the expression of AR pathway genes with an accuracy rate of 100%. Although currently, MBC shared the same therapeutic choice with TNBC-NST, The obvious downregulation of the AR pathway in MBC compared to LAR may contribute to its histologic differentiation and aggressive behavior. Also, our research suggests that another hormonal receptor GPER was upregulated in MBC compared with LAR, possibly due to the suppression of the AR pathway. Meanwhile, it also indicated that MBC can possibly be activated by estrogen even though it lacks the expression of ER, PR and AR. Recent studies revealed that MBC has more tumor-infiltrating lymphocytes and showed higher PD-L1 expression in both tumor cells and stromal lymphocytes. Thus, whether MBC has similarity with the immunomodulatory subtype still need to be elucidated. The more sophisticated classification of TNBC would enable us to have a better understanding of its molecular mechanism and promote the development of precision medicine.

This study was limited by the small sample size used due to the rarity of MBC. Moreover, MBC was considered as a single group in our study although the included MBC cases had different metaplastic components.

Materials and Methods

TCGA Data Acquisition and Cohort Selection

TCGA RNA sequencing level 3 normalized data were downloaded from TCGA Data Portal and imported into R (Version 4.0.3) using TCGAbiolinks (Version 2.16.4) functions GDCquery, GDCdownload and GDCprepare for further analysis. Among cases having immunostaining data of ER, PR and HER2, 122 TNBC cases have been selected, among them, there were 14 cases of MBC. Samples that are molecularly classified as LAR was identified in a previous article using Lehmann classifier and were used in this study. In total, there are 38 cases of the LAR subtype of TNBC.

Analysis of Differentially Expressed Genes

The gene list selected in the analysis of the AR pathway was searched in Pathway Commons database. The Fragments Per Kilobase of transcript per Million mapped reads Upper Quartile (FPKM-UQ) RNA-seq data were log2-transformed before further process. The FPKM-UQ was implemented at the GDC on gene-level read counts that were produced by HTSeq and based on a modified version of the FPKM normalization method. The log2-transformed FPKM-UQ data were analyzed using limma (Version 3.44.3) functions lmFit, eBayes and top Table to identify DEGs between MBC samples and LAR samples. Student t-test was utilized to calculate the P values of genes. Genes with P < 0.05 were considered as DEGs.

Random Forest Analysis

The log2-transformed FPKM-UQ data of DEGs in the MBC and LAR samples were imported into the randomForest function of the randomForest package (Version 4.6-14). The randomForest function implements Breiman’s random forest algorithm for classification, the algorithm yields an ensemble that can achieve both low bias and low variance and effectively avoid overfitting. The MDSplot function was implemented for the multi-dimensional scaling plot of the proximity matrix from randomForest. The number of trees (ntree) was set to be 500 by default. Each tree was grown independently, and the final prediction was yielded by the mean value. 70% of the dataset was taken for training and the rest for testing by default.