2.1 Microarray datasets from the gene expression Omnibus

We conducted a comprehensive search of breast cancer microarray datasets including lung metastasis information in the GEO database from the National Center for Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov/geo/). Only datasets with a sample size greater than 20 were selected for subsequent analysis. Then, the raw intensity files (CEL) of the datasets meeting our criteria for further analyses were downloaded from the GEO database. The robust multiarray average method of the R package “affy” was used to process raw intensity files and generate the gene expression matrices for each selected dataset [13]. The gene expression data of each sample were matched with the clinical information.

2.2 Univariate survival analysis

Log-rank tests for the high and low expression groups of each gene were performed using the R package “survival”. Lung metastasis of breast cancer was considered the outcome event. Genes with a P-value less than 0.01 were deemed candidate genes associated with lung metastasis in breast cancer for modeling.

2.3 Enrichment analysis of GO functions and KEGG pathways

Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted using the WEB-based GEne SeT AnaLysis Toolkit (http://bioinfo.vanderbilt.edu/webgestalt/login.php) to understand the critical biological significance of the identified genes related to lung metastasis in breast cancer.

2.4 Selection of the best genes for constructing a gene prognostic nomogram

Among the genes related to lung metastasis in breast cancer, a robust likelihood-based survival approach was applied to select the best genes for building a gene prognostic nomogram. The whole selection process was implemented with the R package “rbsurv”. Details of the algorithm are summarized as follows:

All samples were randomly divided into the training set with N*(1 − p) samples and the validation set with N*p samples (p = 1/3). Then, the Cox proportional hazards model was used to fit a gene to the training set of samples to obtain the parameter estimate for this gene. Log-likelihood was evaluated with the parameter estimate and the validation set of samples. This process was implemented for each gene.

The above procedure was repeated 10 times; thus, 10 log-likelihoods were obtained for each gene. Next, the best gene g₍₁₎ with the largest mean log-likelihood was selected. All the best lung metastasis survival-related genes were chosen by the robust likelihood-based method.

Let g₍₁₎ be the chosen best gene in the previous step. Adjusting for g₍₁₎, the second best gene was identified by repeating the above two steps. In other words, g₍₁₎ + g_(j) was evaluated for every j, and an optimal two-gene model, g₍₁₎ + g₍₂₎, was chosen. This forward gene selection procedure was continued until fitting was impossible because of the lack of samples. Thusly, a series of K models were built: M₁ = g₍₁₎, M₂ = g₍₁₎ + g₍₂₎, …, M_K−1 = g₍₁₎ + g₍₂₎ + … + g_(K−1), M_K = g₍₁₎ + g₍₂₎ + … + g_(K).

Akaike information criteria (AICs) for all these models were calculated to avoid overfitting, and the optimal model with the smallest AIC was chosen. The model that is best according to AIC is the one that minimizes prediction error [14, 15].

2.5 Construction of the gene prognostic nomogram

The R package “rms” was applied to build the prognostic nomogram based on the expression level of the best genes that were selected by the last step. In the package, the “cph” function was used to construct the COX model. Based on the model, the “nomogram” function was used to generate the prognostic nomogram. The length of the line corresponding to each gene in the prognostic nomogram reflects the contribution of each gene to one patient’s outcome.

2.6 Internal and external validation of the gene prognostic nomogram

After the nomogram was constructed, the training set and the validation set were used as the internal validation dataset and the external validation dataset respectively. For these cohorts, we calculated the area under curve (AUC) and the C-index to test the effectiveness of the gene prognostic nomogram in discriminating the outcome of patients. In addition, we generated Kaplan-Meier curves for the high-risk group and low-risk group determined by the cut-off point of the ROC curve. Univariate and multivariate cox regression were performed with the nomogram and molecular subtypes in validation set. The expression of the 6 genes in the nomogram between tumor and normal tissues was compared at the mRNA level in the TIMER database (https://cistrome.shinyapps.io/timer/). cBioportal for Cancer Genomics was explored to investigate the genetic alterations of the prognostic genes in the model. In validation set, log-rank tests (Events were defined as bone metastases.) were performed for the genes in nomogram to see whether these genes were associated with bone metastases in breast cancer.

2.7 Statistical analyses

Log-rank tests were performed to preliminarily identify candidate genes for modeling with a P-value less than 0.01 using the R package “survival”. The selection of the best genes for constructing a gene prognostic nomogram was implemented with the R package “rbsurv” [14]. The AUC values were calculated using the R package “timeROC” [16]. These analyses were performed using the R Version 3.5.1(http://www.rproject.org). Univariate and multivariate cox regression were performed using SPSS 25.0 (The alpha level was set as 0.05).

3.1 Selection of microarray datasets for further analyses

There were 935 breast cancer datasets in the GEO database of NCBI. Among them, 2 datasets (GSE5327 [17] and GSE2603 [18]) containing complete information about lung metastasis were selected for further analyses. Their platforms are both GPL96. GSE5327 was considered as the training set and GSE2603 was considered as the validation set. The characteristics of all datasets used in this study are shown in Table 1. The 58 cases of GSE5327 are estrogen receptor (ER) negative, no adjuvant treatment and node negative cancer, without data about progesterone receptor (PR) and HER2 (ERBB2). The 82 cases of GSE2603 include ER, PR and HER2 status and have no data about adjuvant treatment.

Table 1

The characteristics of the datasets used in this study.

Dataset	Sample Size	Tissue	Platform
GSE5327	58	Breast Cancer	GPL96
GSE2603	82	Breast Cancer	GPL96

3.2 Genes associated with lung metastasis in breast cancer

We preliminarily identified 319 candidate genes for modeling with a P-value less than 0.01 using the log-rank test (S1 Table). To understand the critical biological significance of the identified genes associated with lung metastasis in breast cancer, enrichment analyses of GO function and KEGG pathways within the identified genes were conducted. The full lists of GO terms are shown in Fig 1A. Of the GO biological process categories, these genes were closely associated with the “biological regulation” and “metabolic process” terms. Of the GO cellular component categories, these genes were closely associated with the “membrane” and “nucleus” terms. Of the GO molecular function categories, these genes were closely associated with the “protein binding” and “ion binding” terms. In addition, the top 10 enriched KEGG pathway terms of the identified genes are listed in Fig 1B.

Fig 1

The enrichment results of the genes associated with lung metastasis.

(A) Each GO Biological Process, Cellular Component and Molecular Function category is represented by a red, blue and green bar, respectively. The height of the bar represents the number of genes observed in the category. (B) The top 10 enriched KEGG pathway terms of the genes associated with lung metastasis.

3.3 Best genes for constructing prognostic nomogram

Applying the partial likelihood of the Cox proportional hazard regression model, we next selected the best lung metastasis-associated genes in breast cancer. We used a cross-validation technique considering the large data variability. Forward selection was implemented to build a series of gene models, and the optimal model was then determined using the minimal AIC. Finally, 6 genes (GOLGB1, TMEM158, CXCL8, MCM5, HIF1AN, and TSPAN31) were selected that could optimally predict the lung metastasis risk of breast cancer patients (Table 2).

Table 2

The best genes predicting LMFS of breast cancer patients.

Gene Symbol	nloglik	AIC	Selected
GOLGB1	79.62	161.23	*
TMEM158	75.41	154.83	*
CXCL8	70.83	147.66	*
MCM5	69.06	146.12	*
HIF1AN	67.28	144.56	*
TSPAN31	65.32	142.65	*
IFT46	64.89	143.78
SLC9A3R1	64.89	145.77
MAPT	64.73	147.46

3.4 The construction of a prognostic nomogram

The R package “rms” was applied to construct the prognostic nomogram based on the expression level of the 6 genes (GOLGB1, TMEM158, CXCL8, MCM5, HIF1AN, and TSPAN31). As shown in Fig 2, "1" represents a high expression level of each gene, and "0" represents a low expression level of each gene. “Points” is the score corresponding to the expression level of a single gene. “Total points” is the sum of the “Points” of the 6 genes, which corresponds to the accurate lung metastasis-free survival rate of each sample. A greater “Total points” value indicates a higher lung metastasis risk for breast cancer patients.

Fig 2

The 6-gene prognostic nomogram based on the expression level of GOLGB1, TMEM158, CXCL8, MCM5, HIF1AN, and TSPAN31.

The high and low expression level of each gene were represented by “1” and “0” respectively. “Points” is the score corresponding to the expression level of a single gene. “Total points” is the sum of the “Points” that 6 genes get. The greater “Total points” value means the higher lung metastasis risk of breast cancer patients.

3.5 Internal and external validation of the prognostic nomogram

The training set GSE5327 and the validation set GSE2603 were used to evaluate the robustness and effectiveness of the gene prognostic nomogram. The estimated median of LMFS is 3542 days (95% CI: 2821–4264) for GSE5327, not reached for GSE2603. We next generated the time-dependent ROC curve (Fig 3) and calculated the AUC and the C-index for these two datasets. The following values were obtained: for the training set GSE5327, the AUC for 3-, 5- and 7-year were 0.94, 0.87 and 0.90 respectively, and the C-index was 0.862 (P<0.0001); for the validation set GSE2603, the AUC for 3-, 5- and 7-year were 0.87, 0.83 and 0.84 respectively, and the C-index was 0.772 (P<0.0001). These results validated the capability of our gene prognostic nomogram to discriminate the outcome of patients according to the prediction risk. On the other hand, each cohort was divided into a high-risk group and a low-risk group determined by the cut-off point of the ROC curve. It should be noted that the Kaplan-Meier curves showed that lung metastasis was more likely to happen in the high-risk group than in the low-risk group, both in the internal validation set (P<0.0001, Fig 4A) and the external validation set (P<0.0001, Fig 4B). Then, we performed univariate and multivariate cox regression with our nomogram and the molecular subtypes (ER, PR and HER2) in GSE2603. The results showed that both our nomogram and ER status were independent factors for breast cancer lung metastasis (Table 3). Hence, our prognostic nomogram based on 6 genes could effectively predict the lung metastasis risk of patients with breast cancer. What’s more, except for CXCL8 which was not found in the TIMER database, the other five genes were significantly differentially expressed between breast cancer and normal breast tissues (Fig 5). In addition, GOLGB1 possessed the most frequent genetic alterations (9%) in cBioPortal for Cancer Genomics (Fig 6). On the other hand, in the validation set GSE2603, we found by log-rank tests that neither the six genes in the nomogram nor the risk predicted by the nomogram were associated with bone metastasis of breast cancer, preliminarily suggesting that these genes and nomogram were specifically associated with lung metastasis of breast cancer (Fig 7). The whole research process is shown in Fig 8.

Fig 3

Performance of the 6-gene prognostic nomogram in discriminating lung metastasis risk of breast cancer patients from the GSE5327 and GSE2603 cohorts.

(A) For the training set GSE5327, the AUC for 3-, 5- and 7-year were 0.94, 0.87 and 0.90 respectively. (B) For the validation set GSE2603, the AUC for 3-, 5- and 7-year were 0.87, 0.83 and 0.84 respectively.

Fig 4

The survival curves of the high-risk and low-risk groups of the GSE5327 and GSE2603 cohorts.

(A) In the GSE5327 cohort, the high-risk group exhibited a higher lung metastasis risk than the low-risk group (P<0.0001, cutoff = 18.7). (B) In the GSE2603 cohort, the high-risk group exhibited a poorer prognosis than the low-risk group (P<0.0001, cutoff = 12.0).

Fig 5

Comparison of expression levels of the 5 prognostic genes between breast cancer and normal breast tissues using data from the TIMER database.

*** < 0.001, ** < 0.01, * < 0.05.

Fig 6

Genetic alterations of the 6 prognostic genes in the nomogram using data from the cBioportal for Cancer Genomics.

Fig 7

The bone metastasis survival curves of the 6 genes and nomogram in GSE2603 cohorts.

(A) GOLGB1. (B) TMEM158. (C) CXCL8. (D) MCM5. (E) HIF1AN. (F) TSPAN31. (G) nomogram. All P-values were greater than 0.05. The cutoff of each gene is its median expression level in the training set. The cutoff of nomogram is consistent with that in Fig 4B.

Fig 8

The process of developing the 6-gene prognostic nomogram.

First, 319 DEGs associated with lung metastasis in breast cancer patients were identified by univariate survival analysis. Enrichment analysis of GO functions and KEGG pathways was conducted based on these genes. Next, a robust likelihood-based survival modeling approach was applied to identify the best genes for prognosis prediction. Then, the gene prognostic nomogram was constructed based on 6 genes (GOLGB1, TMEM158, CXCL8, MCM5, HIF1AN, and TSPAN31). Finally, the 6-gene prognostic nomogram was validated in the training and validation datasets.

Table 3

Univariate and multivariate cox regression with molecular subtypes.

Variables	Univariate cox regression			Multivariate cox regression
	HR	95% CI	P value	HR	95% CI	P value
ER	0.110	0.025–0.491	0.004	0.193	0.038–0.972	0.046
PR	0.017	0.000–1.236	0.062
HER2	0.508	0.114–2.269	0.375
Nomogram	1.094	1.042–1.149	<0.001	1.061	1.003–1.122	0.039