A transmembrane protein family gene signature for overall survival prediction in osteosarcoma

The transmembrane (TMEM) protein family is constituted by a large number of proteins that span the lipid bilayer. Dysregulation of TMEM protein genes widely occurs and is associated with clinical outcomes of patients with multiple tumors. Nonetheless, the significance of TMEM genes in the prognosis prediction of patients with osteosarcoma remains largely unclear. Here, we comprehensively analyzed TMEM protein family genes in osteosarcoma using public resources and bioinformatics methods. Prognosis-related TMEM protein family genes were identified by the univariate Cox regression analysis and were utilized to construct a signature based on six TMEM protein family genes (TMEM120B, TMEM147, TMEM9B, TMEM8A, TMEM59, and TMEM39B) in osteosarcoma. The prognostic signature stratified patients into high- and low-risk groups, and validation in the internal and external cohorts confirmed the risk stratification ability of the signature. Functional enrichment analyses of differentially expressed genes between high- and low-risk groups connected immunity with the prognostic signature. Moreover, we found that M2 and M0 macrophages were the most abundant infiltrated immune cell types in the immune microenvironment, and samples of the high-risk group showed a decreased proportion of M2 macrophages. Single-sample gene set enrichment analysis revealed that the scores of neutrophils and Treg were markedly lower in the high-risk group than these in the low-risk group in The Cancer Genome Atlas and GSE16091 cohorts. As for the related immune functions, APC co-inhibition and cytolytic activity exhibited fewer active levels in the high-risk group than that in the low-risk group in both cohorts. Of the six TMEM genes, the expression of TMEM9B was lower in the high-risk group than in the low-risk group and was positively associated with the overall survival of osteosarcoma patients. In conclusion, our TMEM protein family gene-based signature is a novel and clinically useful prognostic biomarker for osteosarcoma patients, and TMEM9B might be a potential therapeutic target in osteosarcoma.

Introduction

Osteosarcoma is the most common primary bone cancer that mainly affects individuals at the age of 15–19 years, with a second peak in them at 75–79 years (Rickel et al., 2017). It is a very rare disease that has an incidence of about one to three individuals annually per million (Kansara et al., 2014). Osteosarcoma is a highly malignant cancer and has a propensity for local invasion and early metastasis (Sayles et al., 2019; Sheng et al., 2021). The most common site of distant metastasis is the lung, and lung-only metastases account for approximately 82% of metastatic cases (Gill and Gorlick, 2021). The standard therapy for osteosarcoma patients includes surgery resection and chemotherapy, which has led to a dramatic increase in the overall survival of patients with a localized disease, and the five-year survival rate increased from less than 20 to 70% over the past 30 years (Belayneh et al., 2021; Chen et al., 2021). However, unfortunately, 20–30% of osteosarcoma patients are metastatic or recurrent cases, with the five-year survival rate being less than 20% and remaining stagnant (Saraf et al., 2018). The traditional clinicopathological characteristics including gender, age, and TNM stage are the basis for designing management schedules and prognosis prediction for osteosarcoma patients. However, significant variations in the clinical outcomes of patients have been discovered even for those cases that receive the standard management and harbor similar clinicopathological characteristics (Gill and Gorlick, 2021). Thus, biomarkers that could accurately predict the prognosis in osteosarcoma are urgently needed, which might be beneficial for personalized management and clinical decision.

The transmembrane (TMEM) protein family is constituted by a large number of proteins that span the lipid bilayer and these proteins are components of the biological membranes including the membranes of endoplasmic reticulum, lysosome, mitochondria, and the Golgi apparatus (Beasley et al., 2021; Koteluk et al., 2021). TMEMs are wildly expressed in various types of tissues and are supposed to function as channels to permit the transport of different substances across them (Zhang et al., 2022). Though the function of TMEMs remains largely unknown, emerging evidence reveals the vital roles of TMEMs in tumor occurrence and progression (Schmit and Michiels, 2018; Marx et al., 2020). Recent research studies have reported that the aberrantly expressed TMEM genes in tumors could serve as tumor suppressors or oncogenes and that TMEMs were involved in the regulation of cell proliferation, invasion, metastasis, and chemoresistance (Xu et al., 2019; Zhang et al., 2021). Moreover, several TMEMs were correlated with the overall survival of patients and could act as prognostic biomarkers in multiple tumors. For example, Shiraishi et al. (2021) reported that higher expression of TMEM180, a colorectal cancer-specific molecule, predicted worse overall survival in patients with stage III colorectal cancer. A further example includes TMEM106C, which functions as an oncogene, and the upregulation of TMEM106C was associated with poor prognosis in hepatocellular carcinoma (Duan et al., 2021). A better understanding of TMEM protein family genes thus opens perspectives for the identification of prognostic markers in osteosarcoma.

In the present study, we conducted a comprehensive analysis of TMEM protein family genes in osteosarcoma using public resources and bioinformatics methods. The prognosis-related TMEM protein family genes were identified and utilized for the construction of a prognostic signature in osteosarcoma. The risk stratification ability of the TMEM protein family gene-based signature was validated in both internal and external cohorts. Moreover, we analyzed the association of the prognostic signature with tumor immune cell infiltration and the immune microenvironment.

Materials and methods

Data sources

The expression profiles including the RNA-sequencing data and clinical data of patients with osteosarcoma were collected from the cancer genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The GSE16091 dataset containing RNA-seq data and corresponding clinicopathological features was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and was employed for external validation. Patients without complete clinical information were not considered for the study. A total of 86 and 34 osteosarcoma patients were selected from the TCGA and GEO databases for further analysis, respectively.

Signature construction and validation

In the TCGA osteosarcoma cohort (hereafter referred to as the entire cohort), TMEM protein family genes were subjected to the univariate Cox regression analysis to identify prognosis-related genes using the survival package in R upon the threshold of P-value less than 0.05. Then, the entire cohort was randomly separated into a training and a testing cohort at a ratio of approximately 1:1. Then, prognostic TMEM protein family genes were subjected to least absolute shrinkage and selection operator (LASSO) regression analysis to avoid overfitting. Multivariate Cox regression analysis was performed to further screen the candidate genes and calculate the corresponding regression coefficients. A prognostic signature was ultimately constructed based on the linear combination of gene expression levels and regression coefficients. The risk score of each case in the training, testing, entire, and GSE16091 cohorts was calculated using the following formula: $\mathrm{risk}\mathrm{score=}{\displaystyle \sum }_{i=1}^n\left(Coefi\ast Expi\right)$. Here, Coefi is the regression coefficient of the selected TMEM protein family gene and Expi is the expression level of the selected TMEM protein family gene. Patients in each cohort (training, testing, entire, and GSE16091 cohorts) were stratified into high- and low-risk groups according to the median risk score value in the training cohort. Kaplan–Meier survival analysis was performed to compare the overall survival between high- and low-risk groups. The sensitivity and specificity of the TMEM protein family gene-based signature were determined through a time-dependent receiver operating characteristic (ROC) curve analysis.

Nomogram construction

A prognostic nomogram was developed using the rms package to quantitatively analyze the overall survival of osteosarcoma patients. The nomogram integrated clinical factors including gender and age, and the signature-derived the risk score. Calibration curves were plotted to assess the prediction performance of the nomogram by analyzing the consistency of the nomogram-predicted survival with the actual survival.

Functional enrichment analysis

Osteosarcoma patients of the TCGA entire cohort were stratified into high- and low-risk groups based on the prognostic signature. Differentially expressed genes (DEGs) were screened out according to the following criteria: P-value < 0.05 and |log2(fold change)| >1. The clusterProfiler package in R was utilized to perform gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analyses of these DEGs.

Gene set enrichment analysis

For gene set enrichment analysis (GSEA), patients of the TCGA and GSE16091 cohorts were classified into high- and low-risk groups based on the prognostic signature. Enriched pathways in the high- or low-risk groups were identified using the GSEA software (version GSEA 4.0.2). Pathways with NOM P-value < 0.05 and |NES| > 1 were regarded as significantly enriched.

Immune cell infiltration and immune microenvironment

The CIBERSORT algorithm was utilized to quantify the proportions of 22 immune cell subtypes infiltrated in each osteosarcoma sample of the TCGA and GSE16091 cohorts according to the gene expression profiles. Next, differences in the abundance of infiltrated immune cell types were compared between high- and low-risk groups in both cohorts. Moreover, single sample GSEA (ssGSEA) was performed using the gsva package in R to compare the enrichment score of immune cells and related immune functions between different subgroups in both TCGA and GSE16091 cohorts.

Statistical analysis

All the statistical analyses and visualization of the results in the present study were conducted using R software (version 4.1.0) and corresponding packages. The Kaplan–Meier method with a log-rank test was utilized to compare the overall survival in different subgroups. In all instances, differences were considered statistically significant when the P-value was less than 0.05.

Results

Identification of the prognostic transmembrane protein family genes in osteosarcoma

A total of 249 well-defined TMEM protein family genes were enrolled in the present study. At first, we conducted a univariate Cox regression analysis to explore the association of TMEM protein family genes with the overall survival of osteosarcoma patients using the dataset from the TCGA database. As shown in Figure 1A, twenty-six prognostic TMEM protein family genes were identified. Of these 26 genes, 15 genes (TMEM114, TMEM239, TMEM210, TMEM61, TMEM125, TMEM198, TMEM65, TMEM59, TMEM200C, TMEM120B, TMEM229A, TMEM136, TMEM74B, TMEM147, and TMEM119) were regarded as risk factors (hazard ratio >1), while the other 11 genes (TMEM43, TMEM127, TMEM8A, TMEM51, TMEM39B, TMEM251, TMEM9B, TMEM216, TMEM131L, TMEM150B, and, TMEM53) were identified as protective factors (hazard ratio <1). Figure 1B shows the expression profiles of these prognostic TMEM protein family genes and their correlation based on the expression profile is shown in Figure 1C.

FIGURE 1

FIGURE 1. Identification of prognosis-related transmembrane (TMEM) protein family genes in osteosarcoma. (A) Prognosis-related TMEM protein family genes were identified by performing a univariate Cox regression analysis. (B) Expression profiles of the prognosis-related TMEM protein family genes. (C) Correlation of the TMEM protein family genes based on the gene expressions.

Construction of a prognostic signature based on transmembrane protein family genes

In the TCGA training cohort, the prognostic TMEM protein family genes were subjected to LASSO-Cox regression analysis (Figures 2A and B), which led to the construction of a prognostic signature based on six TMEM protein family genes. The regression coefficients of the six genes are shown in Figure 2C. The risk score of each individual was determined by a linear combination of the expression levels of the six genes and their regression coefficients. The formula was as follows: risk score = TMEM120B × 0.523 + TMEM147 × 0.054 + TMEM9B × (−0.422) + TMEM8A × (−0.300) + TMEM59 × 0.077 + TMEM39B × (−0.257). Next, the risk score of each case was calculated and it allowed patients to be assigned to the high- and low-risk groups according to the median risk score value (Figure 2D). Figure 2E shows the survival status and survival time of osteosarcoma patients, and it suggests that the overall survival of patients in the high-risk group seemed to be worse than that in the low-risk group. Figure 2F shows the expression profiles of the six TMEM protein family genes in the high- and low-risk groups. Kaplan–Meier survival analysis demonstrated that patients in the high-risk group had markedly short overall survival compared with patients in the low-risk group (Figure 2G). Then, we performed a time-dependent ROC curve analysis to evaluate the predictive reliability of the prognostic signature. The area under the curve (AUC) values of 1-, 2-, and 3-year overall survival were 0.903, 0.948, and 0.931, respectively (Figure 2H).

FIGURE 2

FIGURE 2. Construction of a TMEM protein family gene-based signature in osteosarcoma. (A,B) LASSO regression and multivariate Cox regression analyses. (C) The regression coefficients of the six TMEM protein family genes. (D,E) The distribution of risk scores and survival status of patients in The Cancer Genome Atlas (TCGA) training cohort. (F) Theheatmap shows the expression profiles of the six genes in the high- and low-risk groups. (G) Kaplan–Meier survival curve for overall survival in the high- and low-risk groups. (H) Time-dependent receiver operating characteristic (ROC) curve for 1-, 2-, and 3-year overall survival.

Validation of the transmembrane protein family gene-based prognostic signature in internal and external cohorts

To evaluate the accuracy of the six TMEM protein family gene-based signature, the testing and entire cohorts were utilized for internal validation. Using the same formula as for the training cohort, the risk score of each patient in the testing and entire cohorts was calculated. Then, patients were divided into high- and low-risk groups using the median risk score in the training cohort as the cutoff value (Figures 3A and B). The survival status of each patient in the testing cohort and the entire cohort is described in Figures 3C and D, and it suggests that the mortality rate was markedly increased in the high-risk group than that in the low-risk group. The expression patterns of the six TMEM protein family genes in the testing cohort and the entire cohort are shown in Figures 3E and F. Kaplan–Meier survival analyses indicated that the overall survival of the high-risk group was significantly shorter than that of the low-risk group (Figures 3G and H). The AUC values of the ROC curves at 1, 2, and 3 years were 0.905, 0.884, and 0.840 in the testing cohort (Figure 3I) and 0.902, 0.908, and 0.891 in the entire cohort (Figure 3J). Furthermore, patients in the entire cohort were classified into different subgroups according to their clinical characteristics (gender and age). Kaplan–Meier survival analysis revealed that patients in the high-risk group had worse overall survival compared with those in the low-risk group, which was consistent in all the subgroups (Figures 4A–D).

FIGURE 3

FIGURE 3. Validation of the TMEM protein family gene-based signature in internal cohorts. (A,B) Distribution of risk scores of patients in the testing and entire cohorts. (C,D) distribution of survival time and survival status of patients in the testing and entire cohorts. (E,F) Expression profiles of the six genes in the testing and entire cohorts. (G,H) Kaplan–Meier survival analysis for the comparison of the overall survival between the high- and the low-risk group. (I,J) Time-dependent ROC curve for 1-, 2-, and 3-year overall survival in the testing and entire cohorts.

FIGURE 4

FIGURE 4. Kaplan–Meier survival curve analysis in different subgroups. (A,B) Comparison of the overall survival between high- and low-risk groups in subgroups classified by gender in the TCGA cohort. (C,D) Comparison of the overall survival between high- and low-risk groups in subgroups classified by age in the TCGA cohort.

We then employed the GSE16091 cohort for external validation. Using the same formula aforementioned, the risk score of each individual in the GSE16091 cohort was determined and these patients were further stratified into high- and low-risk groups according to the median risk score value obtained from the training cohort (Figure 5A). The survival time and survival status of patients in the GSE16091 cohort are shown in Figure 5B, and it suggests that patients in the high-risk group had a higher mortality rate than those in the low-risk group. Figure 5C shows a comparison of the expression levels of the six TMEM protein family genes between high- and low-risk groups. Kaplan–Meier survival analysis demonstrated that the overall survival was worse in the high-risk group than that in the low-risk group (Figure 5D). The AUC values were 0.771, 0.750, and 0.736 at 1, 2, and 3 years, respectively (Figure 5E). Taken together, these analyses revealed the prognostic robustness of the six TMEM protein family gene-based signature.

FIGURE 5

FIGURE 5. Validation of the TMEM protein family gene-based signature in external cohorts. (A) Dot plots exhibit risk score distribution in the GSE16091 cohort. (B) Dot plots show the comparison of the survival time and survival status of osteosarcoma patients in the GSE16091 cohort. (C) Heatmap for comparison of the gene expression levels in the high- and low-risk groups. (D) Kaplan–Meier survival analysis for the comparison of the overall survival between the high- and low-risk group. (E) Time-dependent ROC curve for 1-, 2-, and 3-year overall survival in the GSE16091 cohort.

Independent analysis of the prognostic signature and establishment of a nomogram

To determine whether the TMEM protein family gene-based signature could be an independent prognostic factor for osteosarcoma patients, we conducted univariate and multivariate Cox regression analyses in the TCGA dataset. In both univariate and multivariate Cox regression analyses, the risk score, derived from the prognostic signature, was the only indicator of overall survival in patients with osteosarcoma (Table 1). Moreover, we established a nomogram using the relative risk score, gender, and age as variables, and a higher point was related to worse prognosis on the nomogram (Figure 6A). The calibration diagram suggested that the predictive overall survival of the nomogram showed satisfactory consistency with the actual overall survival (Figure 6B). We also plotted decision curves to assess the clinical utility of the nomogram and found that it harbored comparable net benefit for predicting the 1-, 3-, and 5-year survival rates with the signature-derived risk score (Figure 6C), which further indicated that our signature performed well in predicting the prognosis of osteosarcoma patients.

TABLE 1

TABLE 1. Variables associated with the overall survival in osteosarcoma: univariate and multivariate analyses.

FIGURE 6

FIGURE 6. Construction and validation of a prognostic nomogram in osteosarcoma. (A) A nomogram was developed by integrating age, gender, and relative risk score. (B) Calibration curves exhibit the prediction performance of the nomogram in the TCGA cohort. (C) Decision curves for 1, 2, and 3 years.

Functional analyses based on the prognostic signature

To further investigate the potential function of the TMEM protein family gene-based signature, we utilized R package edgeR to identify DEGs between groups stratified by the risk score. A total of 561 DEGs were screened out between the high- and low-risk groups in the TCGA cohort with the criteria of P-value < 0.05 and |log2FC| > 1. Of the 561 DEGs, 211 were upregulated and 350 were downregulated in the high-risk group compared to those in the low-risk group (Figure 7A). The expression profiles of the DEGs in high- and low-risk groups are shown in Figure 7B. Then, these DEGs were subjected to GO and KEGG enrichment analyses. In the biological process category, DEGs were mainly enriched in positive regulation of cell activation, T-cell activation, positive regulation of leukocyte activation, positive regulation of lymphocyte activation, and humoral immune response. In terms of cellular components, the external side of the plasma membrane, collagen-containing extracellular matrix, and endocytic vesicle were markedly enriched. As for the molecular function, signaling receptor activator activity, receptor ligand activity, and G protein-coupled receptor binding were significantly related to the prognostic signature (Figure 7C). KEGG enrichment analysis revealed that pathways including cytokine–cytokine receptor interaction, cell adhesion molecules, rheumatoid arthritis, hematopoietic cell lineage, and viral protein interaction with cytokine and cytokine receptors were significantly enriched (Figure 7D). Meanwhile, we also performed GSEA to reveal signal pathways associated with the prognostic signature, and the results suggested that immune-related pathways including antigen processing and presentation, B-cell receptor signaling pathway, complement and coagulation cascades, cytokine–cytokine receptor interaction, JAK/STAT signaling pathway, and natural killer cell-mediated cytotoxicity were significantly enriched in both TCGA and GSE16091 cohorts (Figures 8A and B). These analyses illustrated that our TMEM protein family gene-based signature was significantly related to immunity.

FIGURE 7

FIGURE 7. Identification of risk-related differentially expressed genes and functional enrichment analysis. (A) The volcano plot exhibits differentially expressed genes (DEGs) between the high- and low-risk groups. (B) Heatmap showing the expression profiles of the DEGs. (C–D) GO and Kyoto encyclopedia of genes and genomes enrichment analyses of the DEGs.

FIGURE 8

FIGURE 8. Gene set enrichment analysis between high- and low-risk groups in the TCGA cohort (A) and the GSE16091 cohort (B).

Prognostic signature was associated with the immune microenvironment and immune cell infiltration

Based on the functional analyses, we further explore the association of the prognostic signature with the tumor immune microenvironment and immune cell infiltration in osteosarcoma. First, we employed the CIBERSORT algorithm to quantify the infiltrated immune cells in high- and low-risk groups. The proportions of 22 immune cell types in osteosarcoma samples of the TCGA and GSE16091 cohorts are shown in Figures 9A and B, and the results revealed that M2 macrophages and M0 macrophages were the most abundant immune cell types in the immune microenvironment. Moreover, we found that high-risk patients with osteosarcoma had a decreased M2 macrophage proportion in both TCGA and GSE16091 cohorts (Figures 9C and D). Furthermore, we also conducted ssGSEA to compare the enrichment scores of immune cells and related immune functions in high- and low-risk groups. In the TCGA cohort, the scores of most immune cell types were significantly different between the high- and low-risk groups, especially the DCs, macrophages, neutrophils, T-helper cells, and TIL (Figure 9E). In the GSE16091 cohort, the scores of neutrophils and Treg were markedly lower in the high-risk group than those in the low-risk group (Figure 9F), which were consistent in both cohorts. As for the related immune functions, APC co-inhibition and cytolytic activity exhibited fewer active levels in the high-risk group than that in the low-risk group in both cohorts (Figures 9G and H).

FIGURE 9

FIGURE 9. The TMEM protein family gene-based signature was correlated with the tumor immune microenvironment. (A,B) Distribution of 22 immune cell types infiltrated in osteosarcoma samples of the TCGA and GSE16091 cohort. (C,D) Comparison of the infiltrated immune cells in the high- and low-risk groups of the TCGA and the GSE16091 cohort. (E,F) Enrichment scores for the immune cell were compared between the high- and low-risk groups in the TCGA and GSE16091 cohorts. (G,H) Enrichment scores for the related-immune function were compared between the high- and low-risk groups in the TCGA and GSE16091 cohorts.

Expression and Kaplan–Meier survival analyses of the six transmembrane protein family genes

At last, we conducted expression and Kaplan–Meier survival analyses of the six TMEM protein family genes using datasets from the TCGA cohort and the GSE16091 cohort. In both cohorts, the expression levels of TMEM8A and TMEM9B were lower in the high-risk group than those in the low-risk group (Figures 10A and B), while the expression levels of TMEM39B and TMEM147 exhibited no significant difference between the two groups (Figures 10C and F). Besides, the expression levels of TMEM59 or TMEM120B were higher in the high-risk group in the TCGA cohort or the GSE16091 cohort (Figures 10D and E). Moreover, Kaplan–Meier survival analysis revealed that TMEM9B was positively associated with the overall survival of osteosarcoma patients in both TCGA and GSE16091 cohorts (Figure 11B). Higher expression levels of TMEM120B and TMEM147 predicted worse overall survival in patients of the TCGA cohort, while they were not associated with the prognosis of osteosarcoma patients in the GSE16091 cohort (Figures 11E and F). In both cohorts, the expression levels of TMEM8A, TMEM39B, and TMEM59 were not related to the overall survival of osteosarcoma patients in the Kaplan–Meier survival analyses (Figures 11A, C, and D).

FIGURE 10

FIGURE 10. Expression of TMEM8A (A), TMEM9B (B), TMEM39B (C), TMEM59 (D), TMEM120B (E), and TMEM147 (F) in high and low-risk groups.

FIGURE 11

FIGURE 11. Kaplan–Meier survival curve analysis of TMEM8A (A), TMEM9B (B), TMEM39B (C), TMEM59 (D), TMEM120B (E), and TMEM147 (F) in the TCGA and GSE16091 cohorts.

Discussion

Osteosarcoma is a kind of malignant tumor with high heterogeneity and it is still a challenge to accurately predict the clinical outcome of patients with osteosarcoma even if they are under standard management (Schiavone et al., 2019; Zhou et al., 2020). In general, the traditional prognostic stratification method is based on the clinicopathological features of tumors, such as tumor size, site, local invasion, and metastasis. However, the accuracy and sensitivity of these features in predicting the prognosis of osteosarcoma patients are limited. With the popularization of high-throughput sequencing technology in cancer research, a growing number of genomic data were uploaded to open-sourced online platforms (Reuter et al., 2015; Hong et al., 2020). Re-analyzing of datasets in the public databases is being a popular and effective tool for the identification of potential therapeutic targets and prognostic biomarkers. In recent times, an expanding list of biomarkers, based on a predetermined gene expression signature, are identified in multiple cancers (Deng et al., 2021; Xu et al., 2021). For example, Wang et al. (2022) developed a TNF family-based signature in diffuse gliomas with regard to prognosis and further explored the association of the signature with the tumor immune microenvironment.

It is known that TMEM proteins participated in oncogenesis and progression. The upregulated or downregulated TMEMs in cancers contributed to cell proliferation, migration, invasion, and metastasis, and these TMEMs were regarded as prognostic markers in multiple types of cancers and therapeutic targets for cancer treatments (Liu et al., 2019; Marx et al., 2020). However, the expression profiles of TMEM genes as well as the corresponding clinical significance in osteosarcoma remain to be elucidated. We here systematically analyzed TMEM protein family genes in osteosarcoma using a dataset from the TCGA database. A total of twenty-six prognosis-related TMEMs were identified by performing the univariate Cox regression analysis and were further utilized to construct a six-gene signature for prognosis prediction in osteosarcoma. Validation in the internal cohort and external GSE16091 cohort demonstrated that the TMEM gene-based signature could discriminate osteosarcoma prognosis with high accuracy. In addition, the signature-derived risk score was the only independent prognostic factor in osteosarcoma as revealed by the univariate and multivariate Cox regression analyses. We also developed a nomogram by integrating the relative risk score and clinical features including gender and age, which could be utilized to conduct personalized survival prediction for each case with osteosarcoma and might be helpful for designing management schedules and decision-making. These analyses suggested that our TMEM protein family gene-based signature is a novel and clinically useful prognostic marker for osteosarcoma patients.

To reveal the potential biological mechanism of our TMEM protein family gene-based signature, we identified risk-related DEGs and performed a functional enrichment analysis. In the GO analysis, immune-related biological processes such as positive regulation of cell activation, T-cell activation, positive regulation of leukocyte activation, positive regulation of lymphocyte activation, and humoral immune response were significantly enriched. KEGG enrichment analysis suggested that risk-related DEGs were mainly enriched in cytokine–cytokine receptor interaction, cell adhesion molecules, rheumatoid arthritis, hematopoietic cell lineage, and viral protein interaction with cytokine and cytokine receptors. Further annotation of the TMEM protein family gene-based signature via GSEA indicated that immune-related pathways including antigen processing and presentation, B-cell receptor signaling pathway, complement and coagulation cascades, cytokine–cytokine receptor interaction, JAK/STAT signaling pathway, and natural killer cell-mediated cytotoxicity were significantly enriched in both cohorts. Thus, our analyses connected our signature with the immunity in osteosarcoma. Then, we explored the association of the TMEM protein family gene-based signature with tumor immune cell infiltration and immune microenvironment. As revealed by the CIBERSORT algorithm, M2 macrophages and M0 macrophages were the most abundant infiltrated immune cell types in osteosarcoma, and patients at high-risk had a decreased M2 macrophage proportion in both cohorts. M2 macrophages are generally considered to promote tumor growth. However, a recent study also suggested that the presence of CD163-positive M2-polarized macrophages is essential for the inhibition of osteosarcoma progression (Gomez-Brouchet et al., 2017), which is in contrast to what is observed in other solid tumors. Thus, the heterogeneous role of infiltrated macrophages in various types of tumors needed to be further explored. A previous study also verified that higher infiltrated M2 macrophages were associated with improved outcomes in patients with osteosarcoma (Zhang et al., 2020). Therefore, infiltrated M2 macrophages in the tumor immune microenvironment were a predictor of prognosis in osteosarcoma. Moreover, we also performed ssGSEA to compare the enrichment scores of immune cells and related immune functions in high- and low-risk groups. In both cohorts, the scores of neutrophils and Treg were markedly lower in the high-risk group than those in the low-risk group. As for the related immune functions, APC co-inhibition and cytolytic activity exhibited fewer active levels in the high-risk group than those in the low-risk group in both cohorts. Taken together, these analyses suggested the immune-suppressive status of the high-risk group, and our TMEM protein family gene-based signature might be utilized to precite the effect of immune therapy.

Of the six TMEM protein family genes comprised in our signature, the expression level of TMEM9B was lower in the high-risk group than that in the low-risk group and a lower expression of TMEM9B predicted worse overall survival in patients with osteosarcoma, suggesting that TMEM9B might act as a protective factor in osteosarcoma. Besides, we found that the expression of TMEM9B was decreased in osteosarcoma cell lines compared to that in mesenchymal stem cells by analyzing the GSE70414 dataset in the GEO database (data not shown). TMEM9B, mainly located in the lysosome, was a key component of the TNF signaling cascade and was required for the production of proinflammatory cytokines such as IL-6 and IL-8 (Dodeller et al., 2008). In addition, TMEM9B was reported to be the downstream effector of the p53-p21 and p16-pRB tumor suppressor pathways, and cells, silencing of TMEM9B, would bypass senescence (Rovillain et al., 2011). Despite these findings, the role of TMEM9B in malignant tumors had not been reported. Our analyses revealed that TMEM9B might exert a tumor-suppressive role in osteosarcoma. Further experiments are needed to explore the effect of TMEM9B knockdown or overexpression on the malignant behaviors of osteosarcoma cells and cell senescence. TMEM120B and its paralog TMEM120A were located in the nuclear membrane (Ke et al., 2021). TMEM120B and TMEM120A were reported to be highly expressed in fat and participated in the regulation of adipocyte differentiation (Batrakou et al., 2015). Besides, TMEM120B protein was significantly altered in the aqueous humor of patients with primary open-angle glaucoma (Sharma et al., 2018). However, the exact function of TMEM120B in malignant tumors required further exploration. TMEM147 is a highly conserved membrane protein and is wildly expressed in mammalian tissues and cells (Christodoulou et al., 2020). The cellular sublocalization of TMEM147 varies in different kinds of cells (Maimaris et al., 2021). Though the function of TMEM147 remains largely unknown, it is suggested to be implicated in the regulation of cell proliferation, cell apoptosis, and transcription of target genes (Li et al., 2016). In colon cancer, the TMEM147 expression was significantly increased and might represent a biomarker (Feng et al., 2019). Here, we found that the expression of TMEM147 negatively correlated with the prognosis of osteosarcoma patients in the TCGA cohort. The role of TMEM147 deserved further investigation and it might be a novel therapeutic target in osteosarcoma. TMEM59 is a ubiquitously expressed TMEM protein in human tissues and cells (Liu et al., 2020). In function, TMEM59 could interact with FZD, promoting the formation of multimeric WNT-FZD assemblies and positively regulating the activity of WNT signaling (Gerlach et al., 2018). Besides, TMEM59 mediated autophagy through interacting with ATG16L1 (Boada-Romero et al., 2013). Until now, the role of TMEM59 in tumors is relatively understudied and further experiments need to be performed to investigate the function of TMEM59 in the malignant behaviors of tumor cells.

This study has several limitations. First, not all the TMEM protein family genes were included in our study due to their low expression or lack of them in datasets from the TCGA and GEO public databases. Second, all the osteosarcoma patients used in our study were retrospective cases. Validation of the TMEM protein family gene-based signature in a real-world cohort is necessary. Besides, the effect of knockdown or overexpression of the six TMEM genes, especially TMEM9B, on the malignant behaviors of osteosarcoma cells should be further explored and we will conduct it in the future.

In all, we here for the first time investigate the clinical significance of TMEM protein family genes in osteosarcoma. We developed a prognostic signature based on six TMEM protein family genes, which exhibited satisfactory predictive performance in osteosarcoma. Meanwhile, our TMEM protein family gene-based signature was associated with immune cell infiltration and the immune microenvironment. Of the six TMEM protein family genes, lower expression of TMEM9B predicted worse overall survival in patients with osteosarcoma. TMEM9B might be a potential therapeutic target in osteosarcoma.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, and further inquiries can be directed to the corresponding author.

Author contributions

WX designed the study. YD, XZ, and WY conducted bioinformatics analysis, wrote the manuscript, and were responsible for language revisions. All authors reviewed the manuscript.

Acknowledgments

We sincerely thank the TCGA project for use of their data.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Batrakou, D. G., De Las Heras, J. I., Czapiewski, R., Mouras, R., and Schirmer, E. C. (2015). TMEM120A and B: Nuclear envelope transmembrane proteins important for adipocyte differentiation. PLoS One 10, e0127712. doi:10.1371/journal.pone.0127712

PubMed Abstract | CrossRef Full Text | Google Scholar

Beasley, H. K., Rodman, T. A., Collins, G. V., Hinton, A., and Exil, V. (2021). TMEM135 is a novel regulator of mitochondrial dynamics and physiology with implications for human health conditions. Cells 10, 1750. doi:10.3390/cells10071750

PubMed Abstract | CrossRef Full Text | Google Scholar

Belayneh, R., Fourman, M. S., Bhogal, S., and Weiss, K. R. (2021). Update on osteosarcoma. Curr. Oncol. Rep. 23, 71. doi:10.1007/s11912-021-01053-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Boada-Romero, E., Letek, M., Fleischer, A., Pallauf, K., Ramón-Barros, C., and Pimentel-Muiños, F. X. (2013). TMEM59 defines a novel ATG16L1-binding motif that promotes local activation of LC3. Embo J. 32, 566–582. doi:10.1038/emboj.2013.8

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Liu, R., Wang, W., Wang, C., Zhang, N., Shao, X., et al. (2021). Advances in targeted therapy for osteosarcoma based on molecular classification. Pharmacol. Res. 169, 105684. doi:10.1016/j.phrs.2021.105684

PubMed Abstract | CrossRef Full Text | Google Scholar

Christodoulou, A., Maimaris, G., Makrigiorgi, A., Charidemou, E., Lüchtenborg, C., Ververis, A., et al. (2020). TMEM147 interacts with lamin B receptor, regulates its localization and levels, and affects cholesterol homeostasis. J. Cell. Sci. 133, jcs245357. doi:10.1242/jcs.245357

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, Y., Yuan, W., Ren, E., Wu, Z., Zhang, G., and Xie, Q. (2021). A four-methylated LncRNA signature predicts survival of osteosarcoma patients based on machine learning. Genomics 113, 785–794. doi:10.1016/j.ygeno.2020.10.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Dodeller, F., Gottar, M., Huesken, D., Iourgenko, V., and Cenni, B. (2008). The lysosomal transmembrane protein 9B regulates the activity of inflammatory signaling pathways. J. Biol. Chem. 283, 21487–21494. doi:10.1074/jbc.M801908200

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, J., Qian, Y., Fu, X., Chen, M., Liu, K., Liu, H., et al. (2021). TMEM106C contributes to the malignant characteristics and poor prognosis of hepatocellular carcinoma. Aging 13, 5585–5606. doi:10.18632/aging.202487

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, Y., Li, Y., Li, L., Wang, X., and Chen, Z. (2019). Identification of specific modules and significant genes associated with colon cancer by weighted gene co‑expression network analysis. Mol. Med. Rep. 20, 693–700. doi:10.3892/mmr.2019.10295

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerlach, J. P., Jordens, I., Tauriello, D. V. F., van ‘t Land-Kuper, I., Bugter, J. M., Noordstra, I., et al. (2018). TMEM59 potentiates Wnt signaling by promoting signalosome formation. Proc. Natl. Acad. Sci. U.S.A. 115, E3996–e4005. doi:10.1073/pnas.1721321115

PubMed Abstract | CrossRef Full Text | Google Scholar

Gill, J., and Gorlick, R. (2021). Advancing therapy for osteosarcoma. Nat. Rev. Clin. Oncol. 18, 609–624. doi:10.1038/s41571-021-00519-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Gomez-Brouchet, A., Illac, C., Gilhodes, J., Bouvier, C., Aubert, S., Guinebretiere, J.-M., et al. (2017). CD163-positive tumor-associated macrophages and CD8-positive cytotoxic lymphocytes are powerful diagnostic markers for the therapeutic stratification of osteosarcoma patients: An immunohistochemical analysis of the biopsies from the French OS2006 phase 3 trial. Oncoimmunology 6, e1331193. doi:10.1080/2162402x.2017.1331193

PubMed Abstract | CrossRef Full Text | Google Scholar

Hong, M., Tao, S., Zhang, L., Diao, L.-T., Huang, X., Huang, S., et al. (2020). RNA sequencing: New technologies and applications in cancer research. J. Hematol. Oncol. 13, 166. doi:10.1186/s13045-020-01005-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kansara, M., Teng, M. W., Smyth, M. J., and Thomas, D. M. (2014). Translational biology of osteosarcoma. Nat. Rev. Cancer 14, 722–735. doi:10.1038/nrc3838

PubMed Abstract | CrossRef Full Text | Google Scholar

Ke, M., Yu, Y., Zhao, C., Lai, S., Su, Q., Yuan, W., et al. (2021). Cryo-EM structures of human TMEM120A and TMEM120B. Cell. Discov. 7, 77. doi:10.1038/s41421-021-00319-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Koteluk, O., Bielicka, A., Lemańska, Ż., Jóźwiak, K., Klawiter, W., Mackiewicz, A., et al. (2021). The landscape of transmembrane protein family members in head and neck cancers: Their biological role and diagnostic utility. Cancers 13, 4737. doi:10.3390/cancers13194737

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Yuan, C., Wang, L., Lu, M., Wang, Y., Wen, Y., et al. (2016). Transmembrane protein 147 (TMEM147): Another partner protein of Haemonchus contortus galectin on the goat peripheral blood mononuclear cells (PBMC). Parasites Vectors 9, 355. doi:10.1186/s13071-016-1640-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, W. B., Han, F., Huang, Y. S., Chen, H. Q., Chen, J. P., Wang, D. D., et al. (2019). TMEM196 hypermethylation as a novel diagnostic and prognostic biomarker for lung cancer. Mol. Carcinog. 58, 474–487. doi:10.1002/mc.22942

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Ning, J., Zheng, X., Meng, J., Han, L., Zheng, H., et al. (2020). TMEM59 interacts with TREM2 and modulates TREM2-dependent microglial activities. Cell. Death Dis. 11, 678. doi:10.1038/s41419-020-02874-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Maimaris, G., Christodoulou, A., Santama, N., and Lederer, C. W. (2021). Regulation of ER composition and extent, and putative action in protein networks by ER/NE protein TMEM147. Ijms 22, 10231. doi:10.3390/ijms221910231

PubMed Abstract | CrossRef Full Text | Google Scholar

Marx, S., Dal Maso, T., Chen, J.-W., Bury, M., Wouters, J., Michiels, C., et al. (2020). Transmembrane (TMEM) protein family members: Poorly characterized even if essential for the metastatic process. Seminars Cancer Biol. 60, 96–106. doi:10.1016/j.semcancer.2019.08.018

CrossRef Full Text | Google Scholar

Reuter, J. A., Spacek, D. V., and Snyder, M. P. (2015). High-throughput sequencing technologies. Mol. Cell. 58, 586–597. doi:10.1016/j.molcel.2015.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Rickel, K., Fang, F., and Tao, J. (2017). Molecular genetics of osteosarcoma. Bone 102, 69–79. doi:10.1016/j.bone.2016.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Rovillain, E., Mansfield, L., Lord, C. J., Ashworth, A., and Jat, P. S. (2011). An RNA interference screen for identifying downstream effectors of the p53 and pRB tumour suppressor pathways involved in senescence. BMC Genomics 12, 355. doi:10.1186/1471-2164-12-355

PubMed Abstract | CrossRef Full Text | Google Scholar

Saraf, A. J., Fenger, J. M., and Roberts, R. D. (2018). Osteosarcoma: Accelerating progress makes for a hopeful future. Front. Oncol. 8, 4. doi:10.3389/fonc.2018.00004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sayles, L. C., Breese, M. R., Koehne, A. L., Leung, S. G., Lee, A. G., Liu, H.-Y., et al. (2019). Genome-informed targeted therapy for osteosarcoma. Cancer Discov. 9, 46–63. doi:10.1158/2159-8290.Cd-17-1152

PubMed Abstract | CrossRef Full Text | Google Scholar

Schiavone, K., Garnier, D., Heymann, M.-F., and Heymann, D. (2019). The heterogeneity of osteosarcoma: The role played by cancer stem cells. Adv. Exp. Med. Biol. 1139, 187–200. doi:10.1007/978-3-030-14366-4_11

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmit, K., and Michiels, C. (2018). TMEM proteins in cancer: A review. Front. Pharmacol. 9, 1345. doi:10.3389/fphar.2018.01345

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, S., Bollinger, K. E., Kodeboyina, S. K., Zhi, W., Patton, J., Bai, S., et al. (2018). Proteomic alterations in aqueous humor from patients with primary open angle glaucoma. Ophthalmol. Vis. Sci. 59, 2635–2643. doi:10.1167/iovs.17-23434

CrossRef Full Text | Google Scholar

Sheng, G., Gao, Y., Yang, Y., and Wu, H. (2021). Osteosarcoma and metastasis. Front. Oncol. 11, 780264. doi:10.3389/fonc.2021.780264

PubMed Abstract | CrossRef Full Text | Google Scholar

Shiraishi, T., Ikeda, K., Tsukada, Y., Nishizawa, Y., Sasaki, T., Ito, M., et al. (2021). High expression of TMEM180, a novel tumour marker, is associated with poor survival in stage III colorectal cancer. BMC Cancer 21, 302. doi:10.1186/s12885-021-08046-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q.-W., Lin, W.-W., and Zhu, Y.-J. (2022). Comprehensive analysis of a TNF family based-signature in diffuse gliomas with regard to prognosis and immune significance. Cell. Commun. Signal 20, 6. doi:10.1186/s12964-021-00814-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, F., Huang, X., Li, Y., Chen, Y., and Lin, L. (2021). m6A-related lncRNAs are potential biomarkers for predicting prognoses and immune responses in patients with LUAD. Mol. Ther. - Nucleic Acids 24, 780–791. doi:10.1016/j.omtn.2021.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, J., Su, Z., Ding, Q., Shen, L., Nie, X., Pan, X., et al. (2019). Inhibition of proliferation by knockdown of transmembrane (TMEM) 168 in glioblastoma cells via suppression of wnt/β-catenin pathway. Oncol. Res. 27, 819–826. doi:10.3727/096504018x15478559215014

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Zheng, J.-H., Lin, Z.-H., Lv, H.-Y., Ye, Z.-M., Chen, Y.-P., et al. (2020). Profiles of immune cell infiltration and immune-related genes in the tumor microenvironment of osteosarcoma. Aging 12, 3486–3501. doi:10.18632/aging.102824

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, N., Pan, H., Liang, X., Xie, J., and Han, W. (2022). The roles of transmembrane family proteins in the regulation of store-operated Ca2+ entry. Cell. Mol. Life Sci. 79, 118. doi:10.1007/s00018-021-04034-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, S., Dai, H., Li, W., Wang, R., Wu, H., Shen, M., et al. (2021). TMEM116 is required for lung cancer cell motility and metastasis through PDK1 signaling pathway. Cell. Death Dis. 12, 1086. doi:10.1038/s41419-021-04369-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Yang, D., Yang, Q., Lv, X., Huang, W., Zhou, Z., et al. (2020). Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma. Nat. Commun. 11, 6322. doi:10.1038/s41467-020-20059-6

PubMed Abstract | CrossRef Full Text | Google Scholar