RNA-seq data and metadata were downloaded from the TCGA, PanCancer Atlas cohort (https://www.cbioportal.org/), which contains 52 B/AA participants. A case-matched design was performed where the 52 B/AA participants were matched with 52 white participants by choosing the same sex, tumor stage, smoking status, and similar ages. Most outcomes in this study were evaluated using the TCGA cohort. To further assess the association between gene expression levels and OS of lung cancer patients, a second cohort was used (GEO number: GSE101929) (15). This cohort included tumor samples from 16 pairs of B/AA and white non-small cell lung cancer patients, with mRNA expression analyzed through microarray sequencing.

RNA-seq RSEM values were normalized by the VST method and were classified into two levels, i.e., high and low. The Kaplan-Meier survival curve analysis was performed using the ‘survfit’ function in R. P-values were adjusted by the Benjamini–Hochberg procedure. The association between multiple variables in the metadata and OS was tested by the Cox proportional hazards regression using the ‘coxph’ function in R. The association between each variable in the metadata and OS was tested by the Kaplan-Meier survival curve analysis.

Pathway enrichment for the OS-associated genes and genes in each module in co-expression network analysis was measured by the ‘compareCluster’ function in R using the KEGG pathway (16). Pathway enrichment for genes correlated with four potential biomarker genes, i.e., C9orf64, CUBN, REEP2, and MRAS, was tested using the DAVID online tool (17) and enriched biological processes classified by the Gene Ontology database (18) are shown. P-values were adjusted by the Benjamini–Hochberg procedure.

RNA-seq RSEM values were normalized by the average RSEM value of seven housekeeping genes, i.e., ACTB, GAPDH, HPRT1, B2M, TBP, RPL13A, and RPS18. The expression profiles of OS-associated genes were clustered by PCA analysis using the ‘prcomp’ function in R. Alternatively, Bray-Curtis distance among patients was calculated based on the profiles of the 1589 OS-associated genes. Patients were clustered by the NMDS algorithm using the ‘metaMDS’ function in R. The difference of the expression profiles associated with survival status was measured by a PERMANOVA test using the ‘adnois’ function in R.

RNA-seq RSEM values were normalized by the average RSEM value of seven housekeeping genes, i.e., ACTB, GAPDH, HPRT1, B2M, TBP, RPL13A, and RPS18. Unless otherwise specified, to reduce the risk of overfitting in the modeling process, a repeated 2-fold cross-validation with 100 repeats resampling strategy was applied in all predictions. To compare the performance of different algorithms, genes were selected using Spearman’s correlation, and only genes significantly correlated with survival status were used for the prediction. Five algorithms, i.e., random forest, k-Nearest Neighbors, support vector machines with linear kernel, eXtreme gradient boosting, and boosted generalized linear model, were tested using the ‘train’ function in R with the ‘method’ parameter set as ‘rf’, ‘knn’, ‘svmLinear’, ‘xgbLinear’, and ‘glmboost’, respectively. Two feature selection strategies, i.e., Spearman’s correlation and recursive feature elimination (RFE), were performed using ‘cor.test’ and ‘rfe’ functions, respectively, in R.

Spearman’s correlation was employed to measure gene co-expression among VST-normalized RSEM values of all genes. The co-expression network was constructed and visualized using Gephi software (19), applying the ForceAtlas2 algorithm. Modularity analysis was performed using the ‘modularity’ function in Gephi with default parameters to categorize genes into different modules. An overall representation analysis was conducted to assess the enrichment of KEGG pathways in each module using the ‘phyper’ function in R. P-values were adjusted by the Benjamini–Hochberg procedure.

The RSEM values of the RNA-seq data were converted to integers and were input into the DESeq2 software for the differential gene expression analysis (20).

To focus our study on B/AA LUAD patients, we selected the TCGA, PanCancer Atlas cohort (https://www.cbioportal.org/), which contains 52 B/AA participants. A case-matched design was performed where the 52 B/AA participants were matched with 52 white patients by choosing the same sex, tumor stage, smoking status, and similar ages ( Supplementary Table S1 ).

The influence of race, age, sex, tumor stage, and smoking on OS probability was systematically tested by the Cox proportional hazards regression ( Supplementary Figure S1A ). The results showed only a significant association between tumor stage I and OS. Alternatively, the Kaplan-Meier survival curve analysis was performed to test the influence of each individual variable on OS, and similarly, only the tumor stage had a significant association with OS ( Supplementary Figures S1B–F ). Although race (21), age (22), and smoking (23) have been reported to impact OS of lung cancer patients, since the cohort size is small, there could be potential sampling bias and randomness that led to different observations compared with results in large-scale epidemiological studies. Diverse experimental designs could also lead to varied outcomes.

The RNA sequencing data with RSEM values downloaded from the TCGA database was normalized by variance stabilizing transformation (VST). The VST values were classified into two levels, i.e., high and low, and were utilized to determine the association between gene expression and OS. The Kaplan-Meier survival curve analysis discovered four genes, i.e., C9orf64, CUBN, REEP2, and MRAS, with expression significantly (false discovery rate (FDR) ≤ 0.05) associated with OS in the B/AA patients, where high-level expression of C9orf64 shortened the OS but that of CUBN, REEP2, and MRAS led to longer survival ( Figure 1A ). Furthermore, the four identified potential biomarkers were significantly (P- but not FDR-values ≤ 0.05) associated with OS in both B/AA males and females ( Supplementary Figure S2 ) but not in white patients ( Figure 1A and Supplementary Data Sheet 1 ). Thus, the four genes seem to be potential race- but not sex-specific biomarkers. To further assess the reliability of the four identified biomarkers, their associations with OS was evaluated using the same methodology in an independent cohort case-matched between 16 pairs of B/AA and white LUAD patients [GEO number: GSE101929 (15)]. The results in this second cohort were consistent with those observed in the TCGA cohort ( Supplementary Figure S3 ). However, due to the small sample size, the P-values in the second cohort were less significant compared to those in the TCGA cohort. Additionally, since this second cohort was sequenced using microarray technology, C9orf64 was not included in the dataset.

Limited by the sample size, most FDRs of the association between gene expression and OS in the TCGA cohort were larger than 0.05 ( Supplementary Data Sheet 1 ). To explore biological functions associated with OS, the gene set with significant P- but not FDR-values was selected for functional enrichment analysis. Although false-positive results could exist in the analysis outcomes, true-positive results were also included. We identify shared characteristics in the enriched pathways. There were 1589 genes associated with OS (P ≤ 0.05) in the B/AA patients ( Figure 1B ). Surprisingly, 459 genes in the white patients were detected to be associated with OS (P ≤ 0.05), but only 44 OS-associated genes overlapped between the B/AA and white participants ( Figure 1C ).

A functional enrichment analysis on the 1589 B/AA-specific OS-associated genes identified 17 enriched KEGG pathways (P ≤ 0.01). The P-values for the pathway enrichment were insignificant after multiple testing corrections and false-positive outcomes could exist. However, most of the 17 KEGG pathways were related to the immune system, as the top three significantly enriched pathways were B cell receptor signaling, chemokine signaling, and T cell receptor signaling pathways ( Figure 1D ). Other pathways, e.g., Yersinia, Salmonella, and Human immunodeficiency virus 1 infections, might be enriched because genes involved in these pathways overlapped with immunity-related genes in LUAD. Since a similar trend is observed in these enriched pathways, it seems that the immune system plays an essential role in the OS of the B/AA patients. Furthermore, most OS-associated genes in the pathways had high-level expression in the B/AA participants who survived longer ( Figure 1D ). In other words, a more activated immune system in the B/AAs could prolong the time of survival. The only outlier pathway was the PD-L1 expression and PD-1 checkpoint pathway in cancer, where high-level expression of more than half of the genes were associated with shorter survival. When PD-1 on T-cells binds to PD-L1 or PD-L2 on cancer cells, it inhibits the activity of the T-cells (24). Therefore, the result of the PD-1 pathway is consistent with other observations.

The same functional enrichment analysis was performed on the 459 OS-associated genes in the white patients, which only enriched the protein digestion and absorption pathway (P ≤ 0.01) ( Figure 1D ). This pathway seems to serve as a fundamental process that impacts multiple systems, including the immune system. Previous studies have shown a critical influence of the immune system on LUAD in cohorts of predominantly white patients (13, 14). Thus, these results implied that the immune system in the white patients was statistically less essential than in the B/AA LUAD patients in the case-matched TCGA cohort.

Because similar characteristics associated with OS were observed in the pathway enrichment analysis in the B/AA participants ( Figure 1D ), we classified the B/AA patients by the expression of the 1589 OS-associated genes. The RSEM values were normalized relative to a set of housekeeping genes and subjected to PCA analysis. Most B/AA patients who died after 20 months were separated from the survival patients on the PCA plot ( Figure 2A ). Alternatively, the Bray-Curtis distance among the B/AA patients was calculated based on the normalized RSEM values of the 1589 OS-associated genes, and the similarity of the patients was compared by the NMDS algorithm. Again, the B/AA patients were stratified into clearly distinct groups based on their 20-month survival status ( Figure 2B ). Furthermore, an Adonis test demonstrated a significant difference in the Bray-Curtis distance between the expression profiles of the 1589 OS-associated genes among B/AA patients who survived versus those who were deceased after 20 months.

Various machine learning algorithms were employed to predict the 20-month survival status (see details in the Methods). To mitigate overfitting, a resampling strategy involving repeated 2-fold cross-validation with 100 repeats was utilized. In other words, half of the samples were randomly picked up to train a random forest model and the other half were used to evaluate performance. The process was repeated 100 times. To speed up the training process, only genes significantly associated with the 20-month survival status were selected for the modeling. The performance of the prediction was determined by the accuracy. The results suggested that the random forest algorithm had the best performance ( Figure 2C ). Therefore, the random forest algorithm and the same resampling strategy was applied in the modeling process described below.

Considering potential clinical applications, we aimed to streamline the model using a feature selection procedure to identify the most predictive genes. Initially, when all the genes were utilized for prediction, the median area under the ROC curve (auROC) value across 100 repeats was approximately 0.8 ( Figure 2D ). The first feature selection strategy involved a filter approach employing Spearman’s correlation to retain only genes significantly associated with OS. This approach maintained an average of 700 genes in each prediction ( Figure 2E ). Alternatively, a wrapper approach called recursive feature elimination (RFE) was employed, wherein the model was iteratively trained and less important genes were systematically removed, resulting in an average of fewer than 150 genes retained for modeling ( Figure 2F ). The auROC values were slightly but not significantly reduced by the feature selections ( Figure 2D ). Notably, a linear relationship was observed between the likelihood of a particular gene being selected in the recursive feature elimination and the P-value of the gene’s association with OS in the Kaplan-Meier survival curve analysis ( Figure 2G ). These results indicated that genes more significantly associated with OS were more likely to be selected for predicting 20-month survival.

The frequency of genes selected during the RFE procedure for predicting OS in B/AA LUAD patients is shown in Supplementary Data Sheet 3 . To further assess the potential clinical utility of the RFE strategy, only 90 genes, selected at least 20 times across 100 individual RFE iterations ( Supplementary Data Sheet 3 ), were used to predict OS in these B/AA LUAD patients, employing a Leave-One-Out strategy as previously described (25). The resulting auROC of 0.936 demonstrates strong predictive performance with just 90 genes as input features ( Supplementary Figure S4 ). This gene set size is practical for clinical applications, as a custom microarray can efficiently analyze the expression of 90 genes simultaneously.

Overall, the results suggest that the 20-month survival of the 52 B/AA LUAD patients are predictable and genes more significantly associated with OS are more important in the prediction.

Using the same algorithm and resampling strategy, along with all genes, we attempted to predict the survival of the 52 white patients. However, the auROC values obtained were significantly lower compared to the models for B/AA participants ( Figure 2H ). This observation aligns with the findings indicating fewer enriched pathways associated with OS in the white patients, suggesting less shared characteristics in OS-associated gene expression patterns in white than B/AA patients involved in this case-matched cohort.

To investigate the roles of the potential OS-associated biomarker genes in the B/AA patients, a gene co-expression network was constructed using Spearman’s correlation among the VST-normalized RSEM values of all genes. In this network, each node represents a gene, and each edge signifies a significant correlation (FDR ≤ 0.05) between two genes. The network encompassed 10,575 genes, each of which exhibited at least one significant correlation with another gene.

The network was analyzed and visualized by the Gephi software (19). Specifically, the network was generated utilizing the ForceAtlas2 algorithm, and modularity analysis was performed to categorize genes into approximately 80 modules ( Figure 3A ). The top 8 primary modules each comprises at least 279 genes, and the proportions of OS-associated genes within these modules are depicted ( Figure 3B ). Conversely, the remaining modules contain a limited number of genes (≤ 16 genes). To assess the enrichment of OS-associated genes within the 8 modules, an overall representation analysis (ORA) was conducted. This analysis revealed that modules 41 and 45 exhibited significantly higher proportions of OS-associated genes (FDR ≤ 0.05). Notably, modules 41 and 45 were also closely situated in the network ( Figure 3A ), suggesting a more robust correlation between these two modules.

To elucidate the biological functions of the modules, KEGG functional enrichment analysis was conducted. The results revealed that module 41 was significantly enriched (FDR ≤ 0.05) in pathways associated with the immune system, including the chemokine signaling pathway, intestinal immune network for IgA production, cytokine-cytokine receptor interaction, and Th17 cell differentiation ( Figure 4 ). Moreover, most pathways enriched in the OS-associated genes (as shown in Figure 1D ), e.g., B cell receptor signaling, chemokine signaling, T cell receptor signaling pathways, and natural killer cell-mediated cytotoxicity, overlapped with the pathways enriched in module 41. These findings illustrate that the OS-associated genes not only exhibit functional enrichment in the immune system but also demonstrate a relatively higher degree of correlation with other genes involved in the immune system compared to genes in other modules or with other functions. Thus, the results underscore the significance of the immune system in influencing OS in the B/AA patients.

Module 45 contains 15 enriched pathways ( Figure 4 ). The overlapped biological meaning in the 15 pathways revolves around the regulation of cell-ECM interactions, intracellular signaling pathways, and cytoskeletal dynamics, e.g., ECM−receptor interaction pathway, PI3K-Akt and calcium signaling pathway, and regulation of actin cytoskeleton pathway. These pathways have implication in various physiological and pathological processes, including cell adhesion, migration, proliferation, and collectively contribute to the maintenance of cellular homeostasis and the response to extracellular cues in lung cancer (26–28).

C9orf64 has been reported as an epigenetic biomarker that distinguishes patients with NSCLC from those with nonmalignant lung disease at the transcriptional level (29) and has a higher level of methylation in ovarian cancer (30). However, the biological function of C9orf64 has not been characterized in previous studies. In this study, C9orf64 was found to be significantly (FDR ≤ 0.05) correlated with 212 genes. Moreover, these 212 genes exhibited significant (FDR ≤ 0.05) functional enrichment in innate immune response ( Figure 5A ). Therefore, these findings suggest that the expression of C9orf64 is associated with the expression of immunity-related genes. This implies that C9orf64 could potentially serve as a biomarker or play a certain role in the immune system in B/AA LUAD patients.

Similar tests were conducted on CUBN, REEP2, and MRAS. Previous studies have indicated that MRAS can regulate mitogen-activated protein kinase (MAPK)/extracellular signal-regulated kinase and phosphatidylinositol 3-kinase signaling pathways, thereby modulating various cellular processes related to cell growth, survival, and metabolism (31). Consistent with these findings, enrichment analysis revealed that the 113 genes correlated with MRAS in the network were significantly associated with the positive regulation of the protein kinase B and phosphatidylinositol 3-kinase signaling pathways and cell migration ( Figure 5B ). However, the role of MRAS in cancer remains obscure (32). Additionally, there was no pathway significantly enriched in genes correlated with CUBN or REEP2.

A comparable gene co-expression network was constructed for the matched white LUAD patients using the same methods. The majority of genes present in the network for B/AA patients also featured in the network for white participants ( Figure 6A ). However, there was only an 8% overlap between correlations in the networks for the B/AA and white patients ( Figure 6B ). For instance, while 23 C9orf64-associated genes were identified in the white patient network, only one of these genes correlated with C9orf64 in the network for the B/AA participants ( Figure 6C ). The immune response pathway exhibited significant enrichment among the 212 C9orf64-correlated genes in the B/AA patients ( Figure 5A ). Conversely, there was no significantly enriched pathway among the 23 C9orf64-associated genes in the white patient (data not shown). These findings underscore a pronounced racial disparity in gene expression interaction among B/AA and white LUAD patients.

Differentially expressed genes between the B/AA and white LUAD patients were detected by the DESeq2 software (20). The results, including the 20 most significantly differentially expressed genes and the four identified potential B/AA-specific biomarkers in OS, were illustrated in Figure 6D . Interestingly, most genes associated with OS did not show differential expression between the two racial groups ( Figure 6D and Supplementary Data Sheet 1, 2 ). Furthermore, pathway enrichment analysis did not identify any biological functions ( Figure 6E ) that overlapped with OS-associated pathways ( Figure 1D ). These results suggest that the racial disparities highlighted earlier may not stem from differences in gene expression levels but rather from divergent interactions among gene expressions.