WGCNA analysis results of modules and genes related to lung cancer brain metastasis
The Cluster tree displaying the relationship among the biological replicates of REPL samples was as shown in Fig. 1a, a total of 9 samples were clustered in BM sample, they are GSM354029, GSM354030, GSM354025, GSM354027, GSM354032, GSM354028, GSM354033, GSM354026 and GSM354031. The soft-threshold power is a key value used to power the correlation of the genes and affects the mean connectivity and the scale independence of co-expression modules. When the soft-threshold power was β = 3, the scale independence was higher than 0.8 and the mean connectivity was higher (Fig. 1b). ME Diss Thres was set to 0.4 and combined with the dynamic clipping tree algorithm to split the module. After the cluster tree was generated, eight different gene co-expression modules were identified in the REPL, represented by different colors (Fig. 1c). According to the correlation between modules and sample traits (Fig. 1d), 5 modules were correlated with lung cancer BM (P༜0.05), including ME Blue (Core = 0.59, P = 8E-5), ME Black (Core = -0.52, P = 7e-4), ME Brown (Core = 0.8, P = 7e-10), and ME Turquois (Core = 0.89, P = 4E-14) and ME Purple (Core = -0.4, P = 0.01) (Fig. 1e). Among them, ME turquois has the highest correlation with sample traits, showing a positive correlation, with a total of 5724 genes, of which 900 genes satisfy GS> 0.2 and MM> 0.8.
Screening, annotation and expression validation of brain-metastasis-related characteristic genes in lung cancer
After batch processing of sequencing data of GSE14108 and GSE31547, a total of 753 differentially expressed genes were obtained in BM and primary lung cancer samples. Among them, 238 genes were up-regulated in BM samples and 555 genes were down-regulated in BM tissues (Fig. 2a). Furthermore, the intersection of the differentially expressed genes of lung cancer brain metastasis samples and the module genes related to lung cancer brain metastasis resulted in 272 genes, which were used for subsequent analysis (Fig. 2b).As can be seen from the GO enrichment bubble diagram, cell adhesion and immune-related biological processes are the main pathways for the enrichment of key genes in the module, such as positive regulation of cell−cell adhesion, Humoral immune response and positive regulation of cell adhesion (Fig. 2c). The molecular function is mainly related to enzyme activity, including RNA polymerase II activating transcription factor binding, serine-type endopeptidase activity and metallopeptidase activity. As can be seen from the enrichment bar chart of KEGG pathway, the major participating signaling pathways are significantly related to diseases and cell adhesion, such as Pertussis, Rheumatoid arthritis, Malaria, Bladder cancer, Legionellosis, Influenza A, TNF signaling pathway, Focal adhesion, Complement and coagulation cascades and Age-RAGE signaling pathway in diabetic complications (Fig. 2d).The results of GSE126548 differentially expressed genes analysis showed that there were 773 differentially expressed genes between BM / primary lung cancer (Fig. 2e). Jvenn was used to intersect the differentially expressed genes of GSE126548 and the key genes of the module. A total of 14 intersection genes were obtained, indicating that these 14 genes may play an important role in lung cancer brain metastasis (Fig. 2f). The heat map results of further cluster analysis showed that the expressions of the 14 signature genes were significantly different in the transferred and primary samples (Fig. 2g).
LASSO regression analysis and SVM screening of characteristic genes in lung cancer brain metastases
In this study, LASSO regression analysis was used to screen for key genes associated with BM in lung cancer. When lambda. Min was minimal, the best diagnostic model was obtained, which contained 4 characteristic genes (AQP1, GSTA1, ROS1, and TGFB1I1) (Fig. 3a). The receiver operating characteristic (ROC) curve evaluates the accuracy and sensitivity of the diagnostic model, When the area under the ROC curve (AUC) value of the diagnostic model was 1, it indicated that the diagnostic model had high diagnostic value (Fig. 3b). SVM-RFE analysis found that when the number of genes changed from 1 to 14, the error rate of the optimal point for predicting BM lung cancer samples and primary lung cancer samples was 0.00865, and the accuracy rate was 0.991, which corresponding to 9 characteristic genes (ROS1, TGFB1I1, AEBP1, GSTA1, FMO3, AQP1, COMP, OLFM1, MMP1) (Fig. 3c). The intersection of characteristic genes screened by LASSO and SVM was considered as the key genes. The results showed that AQP1, GSTA1, ROS1 and TGFB1I1 may be the diagnostic biomarkers of BMS and primary lung cancer, and are the key genes to study the difference between primary lung cancer and BM (Fig. 3d).
Functional prediction of key genes in lung cancer brain metastases
According to the GSEA function prediction results of AQP1, GO biological processes such as immune response regulation and cell activator chromosome separation were significantly overexpressed in the high-expression group of AQP1. In addition, Cell adhesion, immunity and cytokine related pathways were significantly enriched in the group with high AQP1 expression (Fig. 4a). The GSEA function prediction results of GSTA1 showed that the GO biological processes such as cell division and chromosome separation were significantly low expressed in the high expression group of GSTA1. In addition, Cell cycle, disease-related pathways significantly enriched in high expression group (Fig. 4b). The GSEA function prediction results of ROS1 showed that the immune response was significantly high in the ROS1 high expression group. The prediction results of GSEA function of ROS1 showed that the immune response was significantly overexpressed in the group with high ROS1 expression, Cell division, chromosome separation and other GO biological processes were significantly low expressed in the high expression group of ROS1. In addition, cell cycle, disease infection, phagocytosis, immune-related pathways were significantly enriched in ROS1 high expression group (Fig. 4c). The GSEA function prediction results of TGFB1I1 showed that the biological processes such as immune regulation, lymphocyte activation and lymphocyte-mediated immunity were significantly enriched in the TGFB1I1 high expression group, while the GO biological processes such as cell division and chromosome separation were significantly decreased in the TGFB1I1 high expression group. In addition, pathways related to disease infection and phagocytes were significantly enriched in the TGFB1I1 high expression group (Fig. 4d).
Correlation between key genes of lung cancer brain metastases and immune invasion
The CIBERSORT algorithm and LM22 gene marker were used to analyze the proportion of 22 kinds of immune cells in the tumor samples. According to the calculation results of CIBERSORT, we drawn a heat map of the content of immune cells in each sample with P <0.05 (Fig. 5a). At the same time, Pearson correlation is used to analyze the relationship between immune cells. Fig. 5b shows the specific immune correlation coefficient. For example, the correlation coefficient between Plasma cells and B cells memory is -0.57 (Fig. 5c). Then the immune cell differences based on t-test were shown between BM / primary lung cancer tissues (Fig. 5d). Among them, immune cells B cells naive, Plasma cells, T cells CD8, T cells regulatory, Macrophages M1, Mast cells resting and Mast cells activated have significant differences. We further analyzed the correlation of 4 key genes and 7 immune cells with significant differences between the BM / primary lung cancer groups, and the results showed that the key genes were negatively correlated with T cells regulatory, T cells CD8 and Mast cells activated, and positively correlated with other cells (Fig. 5e).
Prognostic significance of key genes in brain metastases from lung cancer
Analysis of the correlation between the expression of key genes and the overall survival of lung cancer patients shows that AQP1, GSTA1 and TGFB1I1 were significantly correlated with the prognosis of patients with lung cancer, and the higher the gene expression, the worse the prognosis of patients. ROS1 is significantly correlated with the prognosis of lung cancer patients, and the higher the gene expression, the better the prognosis of patients (Fig. 6).