Next Article in Journal
Vaccine-Induced T-Cell and Antibody Responses at 12 Months after Full Vaccination Differ with Respect to Omicron Recognition
Next Article in Special Issue
EDA-E7 Activated DCs Induces Cytotoxic T Lymphocyte Immune Responses against HPV Expressing Cervical Cancer in Human Setting
Previous Article in Journal
Global Implications for COVID-19 Vaccine Series Completion: Insights from Real-World Data from the United States
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Vaccines
Volume 10
Issue 9
10.3390/vaccines10091562
vaccines-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Integration of Bulk and Single-Cell RNA-Seq Data to Construct a Prognostic Model of Membrane Tension-Related Genes for Colon Cancer

by
Jiacheng Li
1,2,†,
Yugang Fu
1,2,†,
Kehui Zhang
1 and
Yong Li
1,*
1
Shanghai Municipal Hospital of Traditional Chinese Medicine, Shanghai University of Traditional Chinese Medicine, Shanghai 200071, China
2
Municipal Medical College of Traditional Chinese Medicine, Shanghai University of Traditional Chinese Medicine, Shanghai 200071, China
*
Author to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Vaccines 2022, 10(9), 1562; https://doi.org/10.3390/vaccines10091562
Submission received: 20 August 2022 / Revised: 6 September 2022 / Accepted: 15 September 2022 / Published: 19 September 2022
(This article belongs to the Special Issue Immunology and Immunotherapy in Cancer)
Browse Figures
Versions Notes

Abstract

:
Background: The plasma membrane provides a highly dynamic barrier for cancer cells to interact with their surrounding microenvironment. Membrane tension, a pivotal physical property of the plasma membrane, has attracted widespread attention since it plays a role in the progression of various cancers. This study aimed to identify a prognostic signature in colon cancer from membrane tension-related genes (MTRGs) and explore its implications for the disease. Methods: Bulk RNA-seq data were obtained from The Cancer Genome Atlas (TCGA) database, and then applied to the differentially expressed gene analysis. By implementing a univariate Cox regression and a LASSO-Cox regression, we developed a prognostic model based on four MTRGs. The prognostic efficacy of this model was evaluated in combination with a Kaplan–Meier analysis and receiver operating characteristic (ROC) curve analysis. Moreover, the relationships between the signature and immune cell infiltration, immune status, and somatic mutation were further explored. Lastly, by utilizing single-cell RNA-seq data, cell type annotation, pseudo-time analysis, drug sensitivity, and molecular docking were implemented. Results: We constructed a 4-MTRG signature. The risk score derived from the model was further validated as an independent variable for survival prediction. Two risk groups were divided based on the risk score calculated by the 4-MTRG signature. In addition, we observed a significant difference in immune cell infiltration, such as subsets of CD4 T cells and macrophages, between the high- and low-risk groups. Moreover, in the pseudo-time analysis, TIMP1 was found to be more highly expressed with the progression of time. Finally, three small molecule drugs, elesclomol, shikonin, and bryostatin-1, exhibited a binding potential to TIMP-1. Conclusions: The novel 4-MTRG signature is a promising biomarker in predicting clinical outcomes for colon cancer patients, and TIMP1, a member of the signature, may be a sensitive regulator of the progression of colon cancer.

1. Introduction

Colon cancer remains one of the world’s most common intestinal diseases; it ranks third in terms of incidence, but second in terms of mortality worldwide [1]. In 2020, an estimated 104,610 new cases of colon cancer and 43,340 cases of rectal cancer were predicted to have occurred in the United States [2]. Due to multiple factors, such as diet, lifestyle, and a lack of health-care infrastructure and resources, the incidence and mortality of colon cancer remain on the rise. We are in an era of the rapid development of cancer-screening methods; however, many patients are diagnosed at an advanced stage and suffer a poor survival rate with only a few effective therapeutic targets for colon cancer [3]. Immune checkpoint inhibitor (ICI) drugs and CAR-T cell therapy have provided durable clinical benefits for metastatic colon cancer [4,5]. Therefore, apart from improvements in chemoradiotherapies and surgical instruments, it is also crucial to develop novel diagnostic and prognostic models, as well as possible therapeutic targets, from emerging sequencing data and computational tools.
Membrane tension, defined as the force per unit length acting on a cross-section of the membrane, regulates many vital biological processes [6]. It consists of two parts: one from the in-plane tension in the lipid bilayer, the other from the adhesion between the membrane and the cytoskeleton [7]. Membrane tension organizes a cell’s shape and its motility, regulating the balance between exocytosis and endocytosis [8]. It is known that derailed endocytosis can majorly contribute to several hallmarks of cancer, not only the sustained proliferation of cancer cells but also its enhanced invasiveness and avoidance of apoptosis [9]. A recent study found that reduced membrane tension is a mechanical hallmark of malignant cells. Cells with lower membrane tension levels are more susceptible to tumor invasion and metastasis [10]. Moreover, researchers found that membrane tension could regulate glycolytic rates in cancer cells, shedding light on how a tumor transforms from stiff to soft [11]. However, MTRGs in colon cancer have not yet been identified. Thus, clinical sample-based screenings for MTRGs are necessary for their possible value in the diagnosis and treatment of colon cancer.

2. Materials and Methods

2.1. Data Acquisition

The bulk RNA-seq data were obtained from the TCGA database (https://tcga-data.nci.nih.gov/tcga/, accessed on 18 April 2022) by a thorough examination of COAD-related datasets. Samples from 473 COAD patients were used for further analysis. All samples were normalized by the fragments per kilobase million (FPKM). A total of 427 patients in the cohort provided information on gene mutations and corresponding clinical characteristics, including age, gender, TNM stage, and survival status, which were also extracted from TCGA. A detailed description of the clinical characteristics of these patients is shown in Supplementary Table S1. In addition, the GSE14333 and GSE103479 cohorts were also obtained from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/, accessed on 18 April 2022) database as independent external cohorts for prognostic model validation. The scRNA-seq data were downloaded from GSE161277, including three normal samples (N1, N2, and N3) and three tumor samples (T1, T2, and T3). The tumor samples included 4118, 4383, and 4137 cells, and the normal samples included 3868, 7549, and 4184 cells. A total of 44 MTRGs were retrieved from previously published literature [10,12,13,14,15,16,17,18,19], whose details have been provided in Supplementary Table S2.

2.2. Construction and Validation of the MTRG Prognostic Model

The prognostic risk of each differentially expressed MTRG was firstly assessed using univariate Cox regression, and then features with p-values < 0.05 in the training cohort were defined as prognostic factors. Next, the optimal gene combinations were screened out by least absolute shrinkage and selection operator (LASSO) regression and multivariate Cox analysis to construct the risk score model. The risk score was calculated as the sum of the products of gene expression levels and their coefficients. Subjects were divided into high- and low-risk groups based on the median risk score as a threshold. Then, Kaplan–Meier analysis and ROC curve analysis were applied to assess the prognostic role of the model. We performed univariate and multivariate Cox regression analyses to determine whether risk scores could be an independent prognostic factor for COAD patients using the “survival” R package. As mentioned above, two independent cohorts from the GSE14333 and GSE103479 were used to verify the utility of the prognostic model (Figure 1).

2.3. Visualization of Mutation and CNV Data

Two waterfall diagrams were depicted with the R package “maftools” [20] to explore the somatic mutations between the high- and low-risk groups of COAD patients. Kaplan–Meier analysis was also performed to compare the survival risks between the two groups. Copy Number Variation (CNV) data of the model genes were downloaded from USCS Xena (https://xenabrowser.net/datapages/, accessed on 18 April 2022). Genes in the CNV region were annotated by the Genome Research Consortium Human build 38 (GRCh38) as the reference genome.

2.4. Protein Interaction Network Construction and GSEA Analysis

GeneMANIA (http://genemania.org/, accessed on 18 April 2022) is an online tool for predicting the interactions and functions of genes and gene sets [21]. In the present study, a 4-MTRG-involved protein interaction network was constructed using this web tool, and other potential proteins associated with MTRG were screened and predicted. Then, the “org.Hs.eg.db”, “clusterProfiler”, and “enrichplot” R packages were used to perform Gene Set Enrichment Analysis (GSEA) to identify biologically significant enrichment pathways between high- and low-risk groups.

2.5. Immune Cell Infiltration, Gene Set Enrichment Analysis, and ICI-Related Gene Expression between the High- and Low-Risk Groups

To clarify the immune status between different groups, 22 types of immune cells were identified in each sample by the CIBERSORT algorithm [22]. The infiltration density of these immune cells was calculated groupwise with the R package “limma”. Then, a single-sample gene set enrichment analysis (ssGSEA) was performed using the “GSVA” and “GSEABase” R packages to assess the infiltration fraction of the known 13 immune cell genomes and the known 16 immune activities of immune-related pathways, respectively. The results were calibrated to the range between 0 and 1. The correlation between tumor microenvironment and risk scores was assessed by the stromal score, ESTIMATE score, immune score, and tumor purity obtained from previous ESTIMATE calculations for each sample. In addition, we analyzed the expression levels of immune checkpoint inhibitors (ICIs) and cuproptosis genes between the two groups.

2.6. Dimensionality Reduction, Clustering, and Annotation of scRNA-Seq Data

“Seurat” R package was adopted to convert previously downloaded 10× scRNA-seq data as a Seurat object. All functions below are inherited from Seurat. Quality control was performed on the raw counts by calculating the percentage of mitochondrial and erythrocyte genes and by excluding low-quality cells, followed by homogenization using the “NormalizeData function”. Top 3000 highly variable features were filtered using the “FindVariableFeatures” function, and normalization was completed by the “ScaleData” function. Principal component analysis (PCA), a preliminary linear dimensionality reduction method, was performed on the scaled data with the elimination of the batch effect using Harmony (v1.0) by default. The t-SNE algorithm, a nonlinear dimensionality reduction technique, was performed for cluster identification. Biologically significant cell types were annotated by the “FindAllMarkers” function to find representative genes for each cluster in combination with typical cell markers.

2.7. Subclusters of Major Cell Types and Pseudo-Time Analysis

To validate the expression of the genes at a single cell transcriptome level, we re-clustered the epithelial cell subgroups with the same method. To determine the lineage of non-malignant cells and malignant cells, the “infercnv” R package was used to determine the malignant clusters. The “monocle” R package was adopted for cell trajectory and pseudo-time analysis, with the method “DDRTree” used for dimension reduction. A heatmap was made to visualize the expression level of MTRGs, and a scatter plot was plotted at last to visualize the change in relative gene expression of the 4-MTRG signature over time.

2.8. Validation of Protein Expression by the HPA Database

Immunohistochemistry (IHC) images from the Human Protein Atlas database (HPA, https://www.proteinatlas.org/, accessed on 18 April 2022) were adopted to validate the protein expression of the MTRGs from the prognostic model. The expression level of TIMP1 was compared between tumor tissues and normal tissues.

2.9. Molecular Docking

The “pRRophetic” R package was used to analyze the differences in drug sensitivity between the high- and low-risk groups. It was further used to screen for relevant active ingredients. The structural formulae of the active ingredients were downloaded from the PubChem (https://pubchem.ncbi.nlm.nih.gov/, accessed on 18 April 2022) database. Chem3D software was used to create 3D structures of the active ingredients. The 3D structures of MTRGs were downloaded from the PDB database (http://www.rcsb.org/, accessed on 18 April 2022). PyMOL software (https://pymol.org/2/, accessed on 18 April 2022) was used to perform operations such as dehydration and hydrogenation of proteins. In addition, AutoDockTools (v1.5.7) software was used to search for ligand-binding pockets. Subsequently, the Vina script was applied to calculate the molecular binding energy and display the molecular docking results. Finally, the results were imported into PyMOL software for visualization.

2.10. Statistical Analysis

Wilcoxon sign tests were used to compare the relationship between two groups for continuous variables. Cox and LASSO-Cox regression were used for predictable models. Kaplan–Meier analyses were used to test survival differences between different risk groups. The Spearman test was used to compare qualitative variables if the value of the variable was small. A two-sided p-value < 0.05 was considered significant. All statistical analyses were performed using R software (version 4.1.3).

3. Results

3.1. Identification of 23 Differentially Expressed MTRGs

To determine the expressional changes of MTRGs in the COAD patients, we quantified the MTRG expression from previously downloaded bulk RNA-seq data from TGCA. We used 41 normal samples and randomly chose 41 out of 473 samples from the COAD patients. The differential analysis reported that the expression levels of MTRGs between the tumor and normal samples were distinct (Figure 2A). A total of 23 MTRGs (9 up and 14 down) were finally identified as differentially expressed genes in the COAD group compared with the normal group (Figure 2B; p < 0.05). Additionally, we examined the correlation among the 23 MTRGs, where FLNA, FNBP1, FGFR1, FERMT2, and CAV1 showed a strong correlation between each other (cutoff > 0.85) (Figure 2C).

3.2. Construction and Validation of a 4-MTRG Prognostic Model

First, we applied a univariate Cox regression to evaluate the prognostic effects of the 23 differentially expressed MTRGs; thus, five genes (EZR, CDC42, EGFR, TIMP1, CAV1) were retained (Figure 3A). Then, LASSO–Cox regression analysis was applied to five candidate genes (Figure 3B,C). Ultimately, we constructed an optimal prognostic model with a 4-MTRG signature consisting of CDC42, EGFR, TIMP1, and CAV1 (Supplementary Table S3). A total of 427 COAD patients were assigned to our prognostic model. Then, we calculated the risk score (RS) of each patient according to the following formula:
RS = ( − 0.043 ) × CDC 42 + 0.032 × EFGR + 0.001 × TIMP 1 + 0.011 × CAV 1
Next, all patients were divided into high- and low-risk groups according to the median value of the risk scores. We made a heatmap of the expression level of the 4-MTRG signature in the high- and low-risk groups. The expression of EGFR, TIMP1, and CAV1 was higher in the high-risk group than in the low-risk group, while the expression of CDC42 was the opposite (Figure 3D). The Kaplan–Meier analysis showed that patients in the high-risk group were associated with a worse overall survival (OS) if compared with patients in the low-risk group (Figure 3E). The receiver operating characteristic (ROC) curve showed the potential of the prognostic model in predicting 1-, 3-, and 5-year OS in the entire cohort (areas under the curve are 0.630, 0.625, and 0.615 (Figure 3F)). We also depicted the distribution of the risk score and survival status among the patients (Figure 3G). Moreover, univariate and multivariable Cox regression analyses were utilized to identify whether the model-derived risk score could be an independent predictor of OS. The results of the univariable regression showed that age, AJCC stage, T stage, N stage, and risk score were closely related to OS (Figure 3H; p < 0.05). Similarly, in the multivariable regression, age, AJCC stage, T stage, and risk score maintained their prognostic power (Figure 3I; p < 0.05). In brief, these data demonstrated that the risk score serves as an independent indicator to predict prognoses for patients in the TGCA cohort. Subsequently, we observed similar results in two external cohorts, GSE14333 and GSE103479 (Figure 4A–C). The Kaplan–Meier analysis and ROC curves of each cohort both suggested that the patients from the high-risk group present a lower OS compared with the low-risk group, which validated the robustness and validity of the prognostic model.

3.3. Somatic Mutation Profile of the Tumor Samples

Previously, we divided patients into two groups based on the median risk score calculated by the 4-MTRG prognostic model; however, we have not obtained the genomic profiling of the tumor samples regarding the somatic alterations that drive cancer progression and patient survival. Here, the copy number variation (CNV) profile is represented group-wise using waterfall graphs. (Figure 5A,B). Both groups showed high mutation alteration rates (97.45% for the high-risk group and 94.97% for the low-risk group). Moreover, we reassigned the patients into two groups based on the median of the tumor mutation burden (TMB), an index that indicates the level of somatic mutation in tumor cells. The high TMB group showed a poorer OS than the low TMB group in the Kaplan–Meier analysis (Figure 5C). When the TMB and prognostic risk score were both included in the survival analysis model, the combination of a high TMB and a high risk behaved the worst (Figure 5D). Next, we concentrated on the frequency of the CNVs of the MTRGs. The CNV analysis of the 23 MTRGs suggested that part of the cleavage enzymes had frequent copy number deletions (Figure 5E), and an overall visualization of the genomic position of the MTRG-related CNV is provided in Figure 5F.

3.4. Protein Interaction Network of the 4-MTRG Signature and GSEA Analysis

Here we constructed a protein interaction network for the 4-MTRG signature from the prognostic model (Figure 6A). The regulatory network carried 24 genes, including four MTRGs and an additional twenty genes that have the potential to interact with them. The additional genes were predicted and added automatically through “GeneMANIA”. The correlation between the risk score and the 4-MTRG signature is shown in Figure 6B. Furthermore, we applied GSEA analysis to investigate the relevant biological processes and signaling pathways in the high-risk group (Figure 6C–H). The results show that cancer hallmark-based gene sets, such as epithelial–mesenchymal transition signaling pathway and the JAK–STAT3 signaling pathway, were highly enriched in the high-risk group. Moreover, the high-risk group is involved in biological processes such as the activation of the immune response, adaptive immune response, and inflammatory response. Furthermore, several classical pathways from the KEGG, Reactome, BioCarta, and PID gene sets, including the toll-like receptor pathway, monocyte pathway, TGF-β pathway, PD-1-signaling pathway, and PI3K-AKT pathway, were also highly related to the high-risk group.

3.5. Risk Scores Related to Different Immune Cell Infiltration, Immune Status, and ICIs

The interaction between the risk genes and the immune status was of a certain prognostic significance for colon cancer. We depicted the immune cell infiltration landscape of the two groups with the CIBERSOFT algorithm. Then, we compared the infiltration density of the immune cells between the high- and low-risk groups. The results suggested that activated CD4+ memory T cells and resting CD4+ memory T cells were downregulated in the high-risk group, while regulatory T cells and M0 macrophages were upregulated in the high-risk group (Figure 7A). Moreover, ssGSEA was performed to evaluate immune activity towards 13 immune cell types and 16 relative immune pathways, and the enrichment analysis output was further compared between the two groups. We found that immune cells in the high-risk group had a more active behavior, in which the numbers of B cells, dendritic cells, macrophages, mast cells, neutrophils, plasmacytoid dendritic cells, T helper cells, follicular helper T cells, tumor-infiltrating lymphocytes, and regulatory T cells were significantly higher than the low-risk group (Figure 7B). Moreover, 7 of 13 immune pathways reported significant differences, especially the APC co-stimulation pathway, CCR pathway, HLA pathway, T cell co-stimulation pathway, and type II IFN response pathway (p < 0.001; Figure 7C). To further validate this finding, we applied four tumor microenvironment-related scoring approaches to evaluate the differences of immune status between the two groups. In total, four scores, including the stromal score (substrate cells in the tumor tissue), immune score (immune cell infiltration in the tumor tissue), and estimate score (the summation of stromal and immune scores from individual cases) were adopted. We found that the estimate score, immune score, and stromal score were all significantly higher in the high-risk group, whereas the tumor purity scores were lower in the high-risk group. (p < 0.001; Figure 7D). We also provided a scatter plot to elucidate the correlation between the risk score and the four scores. The estimate score, immune score, and stromal score showed a positive correlation with the risk score, while the tumor purity score showed a negative correlation (Figure 7E). Recently, by taking advantage of cutting-edge tumor immunology, scientists have developed novel drugs targeted at ICIs for treating solid tumors, and a few products targeted at CTLA4 and PD-1 have been put into clinical practice. To further understand the relationship between our prognostic model and ICIs, 25 ICIs were analyzed between the two groups. We found that risk scores had a positive correlation with the expression of ICIs, and all the expression levels of ICIs showed a statistical significance between groups (Supplementary Figure S1). Cuproptosis, a novel, regulated cell death distinct from known mechanisms broadens our knowledge of the homeostatic conditions of the cell [23], and may be involved in the progression of gastrointestinal tumors; hence, we analyzed the expression level of cuproptosis-related genes between the groups. Ten genes were analyzed, eight were upregulated and one downregulated in the high-risk group (Figure 7F), indicating that the cell might have been in a highly active mode and the death mechanism might have been inhibited.

3.6. scRNA-Seq and Cell Type Characteristics of Normal and COAD Samples

We downloaded 10× scRNA-seq data of three tumor samples and three normal samples from the GSE161277 dataset. A total of 23,634 cells were identified after quality control. The top 3000 highly variable characteristics were selected after normalization; the batch effect was calibrated with the “Harmony” package. We then identified the typical markers of 25 distinct clusters after a PCA and t-SNE analysis, and 17 clusters were shown after merging and annotating manually (Figure 8A). Then, we annotated the clusters based on several canonical marker genes for known cell lineages: T cell (CD3D), NK cell (FGFBP2), follicular B cell (MS4A1), native B cell (TCL1A), plasma B cell (MZB1), monocyte (CD14), macrophage (CD68), goblet cell (MUC2), fibroblast (DCN), endothelial cell (VWF), epithelial mix (EPCAM and PROM1), transit-amplifying cell (MKI67), enteroendocrine (CHGA), enterocyte (GUCA2A and GUCA2B), goblet-progenitor cell (SPINK4), absorptive cell (AQP8 and SLC26A3), and mast cell (KIT) (Figure 8B). To obtain a general view of the distribution of the four MTRGs from the prognostic model, we mapped their distribution with a t-SNE plot. The result showed that CDC42, a regulator of both the architecture and motility of the plasma membrane, appeared most broadly; TIMP1 mostly appeared in epithelial mix, macrophages, and monocytes; and CAV1 and EGFR were distributed sparsely (Figure 8C). Notably, epithelial cells were the only cell type that all the 4-MTRGs exhibited expression in. Next, we labeled the epithelial mix either by patient number (P1–P3) or sample type (normal or tumor). Then, we made bar charts (Figure 8D) to present the proportion of the 17 cell types categorized by their labels. We also made t-SNE plots to present the proportion that each label accounts for in the whole area. We could see that a large portion of epithelial mix, transit-amplifying cells, and goblet progenitor cells were from the tumor samples, and the distribution of the different cell types among the patients were heterogeneous (Figure 8E,F).

3.7. Subclusters of the Epithelial Cell and Pseudo-Time analysis

As we know from the descriptive analysis of the previous section, all the four prognostic model-derived MTRGs expressed in the epithelial mix. Therefore, we re-clustered the epithelial mix following the same procedure; the mix was divided into nine clusters (Figure 9A). Then, the “infercnv” R package was used to estimate the proportion of malignant cells based on the clusters. A total of 62.8% of the cells of the epithelial mix were predicted to be malignant (Figure 9B; Supplementary Figure S2). Then, the “monocle” R package was exploited to analyze the cell trajectory and pseudo-time of the epithelial mix. We observed that the epithelial mix was present in all stages and transformed into malignant cells along the timeline (Figure 9C–F). In addition, we created a heatmap to visualize the expression levels of all 23 MTRGs in the epithelial cells; 18 genes reported expressional changes with statistical significance (Figure 9G). Then, we tested the relative expressions of the 4-MTRG signature in the epithelial mix and found that only the expression of TIMP1 had a higher relative increase than at initiation (Figure 9H). With the help of the HPA database, we identified that TIMP1 has a high expression in colorectal tumor tissue (Figure 10). Overall, we found that three genes in the prognostic model (CDC42, CAV1, and TIMP1) had an increasing degree of expression in the epithelial mix from the patients’ samples over time, among which the expression of TIMP1 had a threefold increase at the final stage.

3.8. Drug Sensitivity of the Two Groups and Docking of Drug Candidates to TIMP1

To further explore the difference in the drug-resistance potentials between the high- and low-risk groups, we compared the estimated IC50 levels of small molecule drugs in the two groups with the “pRRophetic” R package. Among these potential drugs, we found that the sensitivity of elesclomol, shikonin, and bryostatin-1 showed significant differences between groups (p < 0.001; Supplementary Figure S3). To further investigate if the potential drugs could bind to TIMP1, the structures of these drugs were downloaded from the PubChem database (Figure 11A–C), and then they were docked to TIMP1. The molecular-docking results showed relatively high affinity scores between the drugs and TIMP1 (Supplementary Table S4). The formation of hydrogen bonds could be observed between the drugs and the predicted pocket of TIMP1 (Figure 11D–F).

4. Discussion

Currently, researchers have developed considerable novel prognostic models for colon cancer based on clinical features, immune cell infiltration, subgroups of the TMN stage, the expression level of mRNA, and non-coding RNA [24,25,26,27]. However, due to the heterogeneity of the disease, prognostic models with high specificity and accuracy are still rare; therefore, more models need to be developed and tested in clinical practice.
MTRGs are a set of membrane tension-related genes that may be involved in carcinogenesis and metastasis. In this study, we retrieved 43 MTRGs from previously published literature. Then, we developed a novel 4-MTRG prognostic model for colon cancer based on differential gene analysis and Cox regression. The signature consists of four genes: EGFR, TIMP1, CAV1, and CDC42. EGFR is a known regulator of colon cancer contributing to tumor carcinogenesis and progression [28]. A recent study found that macrophages marked with phosphorylated EGFR play a crucial role in the development of the inflammation-mediated stages of colon carcinogenesis [29]. Research also revealed that EGFR cooperates with other genes such as HER-2 and YB-1 to promote cell growth and survival [30,31]. As an MTRG, EGFR may regulate membrane tension during its physiological trafficking. Therefore, cancerous cells could take advantage of such a functional link to enhance their oncogenic influence [17]. TIMP1 is a key element in the regulation of ECM remodeling and is secreted through the plasma membrane via exocytosis [32]. Meanwhile, the expression of TIMP1 and CAV1 increases in cirrhosis and hepatocellular carcinoma [12]. This may suggest that these MTRGs may cooperate in regulating the stiffness of tumors. CAV1 is considered a cell surface protein involved in the formation of caveolae, small invaginations of the plasma membrane [33]. An elevated level of CAV1 expression may contribute to colorectal tumor progression by enhancing aerobic glycolysis in colon cancer cells [34]. However, the regulatory mechanism of CAV1 in colon cancer in terms of membrane tension is not well understood. CDC42 has emerged as a key player in cancer metastasis due to its roles in regulating cell division and actin cytoskeletal rearrangements [35]. CDC42 accumulates at the cell’s leading edge dependent on membrane traffic, and controls the dilation of the exocytotic fusion pore by regulating membrane tension [36,37]. Though CDC42 was reported highly expressed in colon cancer and was downregulated by the potential tumor suppressor gene ID4, its role in the metastasis of colon cancer is still unknown [38].
Based on the risk score calculated in the prognostic model, we divided the patients in the TGCA cohort into two risk groups, considered high- and low-risk. The model exhibited accuracy and robustness in predicting survival rates in the TGCA cohort and two validation cohorts. Our 4-MTRG prognostic model is easy to use and might increase the accuracy of survival probability predictions for colon cancer patients. Then, we conducted a series of enrichment analyses and immune-related analyses to further evaluate the mechanism of how this signature regulates colon cancer and interacts with the microenvironment. The result of the GSEA analysis revealed that typical inflammatory pathways such as JAK-STAT3 and novel pathways such as PD-1 rank highly in the high-risk group. Essential inflammatory pathways such as STAT3 are understood to be drivers of inflammation and tumorigenesis, and diving into the crosstalk of these pathways advances our understanding of the complex links between inflammation and colon cancer [39]. In the immune-related analysis, first, we applied the CIBERSOFT algorithm to depict the immune cell infiltration landscape, and we observed less infiltration of CD4 T cells and more infiltration of M0 macrophages and regulatory T cells (Treg) in the patients of the high-risk group. CD4 T cells promote colon cancer stemness via the regulation of stemness genes, which negatively affects patient outcome [40]. The role of Tregs in the progression of colon cancer is controversial: it was shown that Tregs present antitumor immunity through the production of cytokines, such as TGF-β in the early stage, while the infiltration of tumors by Tregs confers growth and metastatic advantages by inhibiting antitumor immunity as the stage progresses [41]. The significance of the high infiltration of M0 macrophages is yet unknown. In the subsequent ssGSEA analysis, we noticed that immune cells in the high-risk group were more active, and many immune-related hallmarks were enriched, such as the T cell co-stimulation pathway and the type II IFN response pathway. This observation corresponds to a recent study showing that different T-cell subsets can express IFN-γ, thereby altering the immune responses in the colorectal cancer microenvironment [42]. Then, four immune-related scores were adopted to evaluate the immune status of the two groups, and the result showed that scores in the high-risk group were significantly high. Previous studies revealed that high immune and stromal scores as well as a high degree of infiltration in macrophages were associated with a poor prognosis of colon cancer [43]. In addition, we evaluated the relationship between the risk score and ICIs, and the relationship between the risk score and the cuproptosis-related genes. The result suggests that patients with high-risk scores might benefit most from drugs targeting ICIs and cuproptosis-related genes.
In the scRNA-seq analysis, we described the landscape of the samples and annotated the clusters with their typical cell markers. We noticed that all the 4-MTRGs exhibited expression in the epithelial mix. Therefore, we re-clustered the epithelial mix and conducted a pseudo-time analysis. The result showed that over half of the epithelial mix was estimated to be malignant by the “infercnv” package. In particular, TIMP1 was the only MTRG that responded to the shift of stage and had a relative threefold increase in expression with the progression of time. TIMP1 is not only considered a prognostic biomarker for various cancers [44] but promotes tumor progression. For example, an increased degree of TIMP1 expression promotes the in vivo growth of both cancer types and stimulates the accumulation of cancer-associated fibroblasts [45]. However, how membrane tension is involved in the dysfunction of TIMP1 is still unknown. Finally, we compared the drug sensitivity levels between the high- and low-risk groups and selected three potential drugs, elesclomol, shikonin, and bryostatin-1. Then, we conducted molecular docking on these drugs to TIMP1. The result showed that all three drugs had an adequate affinity score and the formation of hydrogen bonds was observed. TIMP1 was found to decrease tumor cell sensitivity to multiple anticancer drugs by activating downstream pathways and it also exhibited anti-apoptotic activity [46]. Therefore, discovering and developing novel drugs dependent on alternative pathways of cell death programming is essential. Copper metabolism has a vital role in tumorigenesis and elesclomol, as a copper chelator, could inhibit colon cancer cells by targeting ATP7A and regulating ferroptosis [47]. The second drug, shikonin, was reported to induce apoptosis and autophagy in colon cancer cells by targeting galectin-1 and the JNK-signaling pathway [48]; however, its effect on TIMP1 has not been tested. Overall, we hope these drugs will bind to TIMP1, relieving the progression of colon cancer and the occurrence of drug resistance.
There are still some limitations that must be addressed. Given that our results are based on sequencing datasets and computational simulations, further experimental and clinical research is needed to evaluate the effect of our 4-MTRG prognostic model. Moreover, the underlying mechanisms of these four selected genes in our model should be further explored to help us understand the intrinsic mechanisms involved in the tumorigenesis, progression, and metastasis of colon cancer. Further research on the dynamics of the plasma membrane and its relationship with broad-sense carcinogenesis and metastasis is also necessary.

5. Conclusions

This study demonstrates that MTRGs can be used to classify colon cancer patients based on different clinical and molecular features. A robust 4-MTRG signature model has been presented, which can predict the prognosis of colon cancer patients. Moreover, this study also found that TIMP1, as a membrane tension-related gene, may serve as a therapeutic target in colon cancer patients. The findings of our study provide insights in predicting the prognosis of colon cancer patients, and the drugs we have selected might contribute to treatments in clinical practice.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/vaccines10091562/s1, Table S1. Demographic information of COAD patients in the present study. Table S2. 44 membrane tension-related genes. Table S3. Multivariate COX regression analysis results of model genes. Table S4. Molecular docking scores of the three compounds. Figure S1. The boxplots for the comparison of the immune checkpoints genes between the high-risk and low-risk groups in the colon cancer patients. Figure S2. InferCNV analysis in epithelial cell subsets. Figure S3. Drug sensitivity analysis in the high-risk and low risk groups in colon cancer patients.

Author Contributions

J.L. and Y.L. contributed to the conception and design of the study. J.L. organized the database, performed bioinformatical analysis and visualization. K.Z. checked the statistical analysis. J.L. and Y.F. wrote the first draft of the manuscript. All authors contributed to manuscript revision. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (NO. 81573775; NO. 81873157).

Institutional Review Board Statement

Not applicable.

Data Availability Statement

The original contributions presented in the study are included in the Supplementary Material. Further inquiries can be directed to the corresponding authors.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ahmed, M. Colon Cancer: A Clinician’s Perspective in 2019. Gastroenterol. Res. 2020, 13, 1–10. [Google Scholar] [CrossRef] [PubMed]
  2. Benson, A.B.; Venook, A.P.; Al-Hawary, M.M.; Arain, M.A.; Chen, Y.-J.; Ciombor, K.K.; Cohen, S.; Cooper, H.S.; Deming, D.; Farkas, L.; et al. Colon Cancer, Version 2.2021, NCCN Clinical Practice Guidelines in Oncology. J. Natl. Compr. Cancer Netw. 2021, 19, 329–359. [Google Scholar] [CrossRef] [PubMed]
  3. Bretou, M.; Jouannot, O.; Fanget, I.; Pierobon, P.; Larochette, N.; Gestraud, P.; Guillon, M.; Emiliani, V.; Gasman, S.; Desnos, C.; et al. Cdc42 controls the dilation of the exocytotic fusion pore by regulating membrane tension. Mol. Biol. Cell 2014, 25, 3195–3209. [Google Scholar] [CrossRef]
  4. Chen, B.; Khodadoust, M.S.; Liu, C.L.; Newman, A.M.; Alizadeh, A.A. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. In Cancer Systems Biology, Methods in Molecular Biology; Von Stechow, L., Ed.; Springer: New York, NY, USA, 2018; pp. 243–259. [Google Scholar] [CrossRef]
  5. Cheng, W.-L.; Feng, P.-H.; Lee, K.-Y.; Chen, K.-Y.; Sun, W.-L.; Van Hiep, N.; Luo, C.-S.; Wu, S.-M. The Role of EREG/EGFR Pathway in Tumor Progression. Int. J. Mol. Sci. 2021, 22, 12828. [Google Scholar] [CrossRef]
  6. Chronopoulos, A.; Thorpe, S.D.; Cortes, E.; Lachowski, D.; Rice, A.J.; Mykuliak, V.V.; Róg, T.; Lee, D.A.; Hytönen, V.P.; Del Río Hernández, A.E. Syndecan-4 tunes cell mechanics by activating the kindlin-integrin-RhoA pathway. Nat. Mater 2020, 19, 669–678. [Google Scholar] [CrossRef] [PubMed]
  7. Chugh, P.; Clark, A.G.; Smith, M.B.; Cassani, D.A.D.; Dierkes, K.; Ragab, A.; Roux, P.P.; Charras, G.; Salbreux, G.; Paluch, E. Actin cortex architecture regulates cell surface tension. Nat. Cell Biol. 2017, 19, 689–697. [Google Scholar] [CrossRef]
  8. De Belly, H.; Stubb, A.; Yanagida, A.; Labouesse, C.; Jones, P.H.; Paluch, E.K.; Chalut, K.J. Membrane Tension Gates ERK-Mediated Regulation of Pluripotent Cell Fate. Cell Stem Cell 2021, 28, 273–284.e6. [Google Scholar] [CrossRef]
  9. Deng, X.; Lin, D.; Zhang, X.; Shen, X.; Yang, Z.; Yang, L.; Lu, X.; Yu, L.; Zhang, N.; Lin, J. Profiles of immune-related genes and immune cell infiltration in the tumor microenvironment of diffuse lower-grade gliomas. J. Cell. Physiol. 2020, 235, 7321–7331. [Google Scholar] [CrossRef]
  10. Diz-Muñoz, A.; Fletcher, D.A.; Weiner, O.D. Use the force: Membrane tension as an organizer of cell shape and motility. Trends Cell Biol. 2013, 23, 47–53. [Google Scholar] [CrossRef] [Green Version]
  11. Du, W.; Frankel, T.L.; Green, M.; Zou, W. IFNγ signaling integrity in colorectal cancer immunity and immunotherapy. Cell. Mol. Immunol. 2022, 19, 23–32. [Google Scholar] [CrossRef]
  12. Fu, Z.Y.; Lv, J.H.; Ma, C.Y.; Yang, D.P.; Wang, T. Tissue inhibitor of metalloproteinase-1 decreased chemosensitivity of MDA-435 breast cancer cells to chemotherapeutic drugs through the PI3K/AKT/NF-kB pathway. Biomed. Pharmacother. 2011, 65, 163–167. [Google Scholar] [CrossRef] [PubMed]
  13. Gao, W.; Huang, Z.; Duan, J.; Nice, E.C.; Lin, J.; Huang, C. Elesclomol induces copper-dependent ferroptosis in colorectal cancer cells via degradation of ATP7A. Mol. Oncol. 2021, 15, 3527–3544. [Google Scholar] [CrossRef] [PubMed]
  14. Giannopoulou, E.; Antonacopoulou, A.; Floratou, K.; Papavassiliou, A.G.; Kalofonos, H.P. Dual targeting of EGFR and HER-2 in colon cancer cell lines. Cancer Chemother. Pharmacol. 2009, 63, 973–981. [Google Scholar] [CrossRef] [PubMed]
  15. Del Pulgar, T.G.; Valdes-Mora, F.; Bandrés, E.; Pérez-Palacios, R.; Espina, C.; Cejas, P.; García-Cabezas, M.A.; Nistal, M.; Casado, E.; González-Barón, M.; et al. Cdc42 is highly expressed in colorectal adenocarcinoma and downregulates ID4 through an epigenetic mechanism. Int. J. Oncol. 2008, 33, 185–193. [Google Scholar] [CrossRef]
  16. Gong, Y.; Scott, E.; Lu, R.; Xu, Y.; Oh, W.K.; Yu, Q. TIMP-1 Promotes Accumulation of Cancer Associated Fibroblasts and Cancer Progression. PLoS ONE 2013, 8, e77366. [Google Scholar] [CrossRef] [PubMed]
  17. Ha, T.-K.; Chi, S.-G. CAV1/caveolin 1 enhances aerobic glycolysis in colon cancer cells via activation of SLC2A3/GLUT3 transcription. Autophagy 2012, 8, 1684–1685. [Google Scholar] [CrossRef] [PubMed]
  18. Hardbower, D.M.; Coburn, L.A.; Asim, M.; Singh, K.; Sierra, J.C.; Barry, D.P.; Gobert, A.P.; Piazuelo, M.B.; Washington, M.K.; Wilson, K.T. EGFR-mediated macrophage activation promotes colitis-associated tumorigenesis. Oncogene 2017, 36, 3807–3819. [Google Scholar] [CrossRef]
  19. Issa, I.A.; Noureddine, M. Colorectal cancer screening: An updated review of the available options. World J. Gastroenterol. 2017, 23, 5086–5096. [Google Scholar] [CrossRef]
  20. Keren, K. Membrane tension leads the way. Proc. Natl. Acad. Sci. USA 2011, 108, 14379–14380. [Google Scholar] [CrossRef] [Green Version]
  21. Kryczek, I.; Lin, Y.; Nagarsheth, N.; Peng, D.; Zhao, L.; Zhao, E.; Vatan, L.; Szeliga, W.; Dou, Y.; Owens, S.; et al. IL-22+CD4+ T cells promote colorectal cancer stemness via STAT3 transcription factor activation and induction of the methyltransferase DOT1L. Immunity 2014, 40, 772–784. [Google Scholar] [CrossRef]
  22. Lachowski, D.; Matellan, C.; Gopal, S.; Cortes, E.; Robinson, B.K.; Saiani, A.; Miller, A.F.; Stevens, M.M.; Hernández, A.E.D.R. Substrate Stiffness-Driven Membrane Tension Modulates Vesicular Trafficking via Caveolin-1. ACS Nano 2022, 16, 4322–4337. [Google Scholar] [CrossRef] [PubMed]
  23. Long, A.G.; Lundsmith, E.T.; Hamilton, K.E. Inflammation and Colorectal Cancer. Curr. Colorectal Cancer Rep. 2017, 13, 341–351. [Google Scholar] [CrossRef] [PubMed]
  24. Magee, M.S.; Abraham, T.S.; Baybutt, T.R.; Flickinger, J.C.; Ridge, N.A.; Marszalowicz, G.P.; Prajapati, P.; Hersperger, A.R.; Waldman, S.A.; Snook, A.E. Human GUCY2C-Targeted Chimeric Antigen Receptor (CAR)-Expressing T Cells Eliminate Colorectal Cancer Metastases. Cancer Immunol. Res. 2018, 6, 509–516. [Google Scholar] [CrossRef] [PubMed]
  25. del Maldonado, M.M.; Medina, J.I.; Velazquez, L.; Dharmawardhane, S. Targeting Rac and Cdc42 GEFs in Metastatic Cancer. Front. Cell Dev. Biol. 2020, 8. [Google Scholar] [CrossRef]
  26. Mayakonda, A.; Lin, D.-C.; Assenov, Y.; Plass, C.; Koeffler, H.P. Maftools: Efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018, 28, 1747–1756. [Google Scholar] [CrossRef]
  27. Mellman, I.; Yarden, Y. Endocytosis and Cancer. Cold Spring Harb. Perspect. Biol. 2013, 5, a016949. [Google Scholar] [CrossRef]
  28. MMercier, V.; Larios, J.; Molinard, G.; Goujon, A.; Matile, S.; Gruenberg, J.; Roux, A. Endosomal membrane tension regulates ESCRT-III-dependent intra-lumenal vesicle formation. Nat. Cell Biol. 2020, 22, 947–959. [Google Scholar] [CrossRef]
  29. Nagasu, S.; Sudo, T.; Kinugasa, T.; Yomoda, T.; Fujiyoshi, K.; Shigaki, T.; Akagi, Y. Y-box-binding protein 1 inhibits apoptosis and upregulates EGFR in colon cancer. Oncol. Rep. 2019, 41, 2889–2896. [Google Scholar] [CrossRef]
  30. Osmani, N.; Peglion, F.; Chavrier, P.; Etienne-Manneville, S. Cdc42 localization and cell polarity depend on membrane traffic. J. Cell Biol. 2010, 191, 1261–1269. [Google Scholar] [CrossRef] [Green Version]
  31. Park, J.S.; Burckhardt, C.J.; Lazcano, R.; Solis, L.M.; Isogai, T.; Li, L.; Chen, C.; Gao, B.; Minna, J.D.; Bachoo, R.; et al. Mechanical regulation of glycolysis via cytoskeleton architecture. Nature 2020, 578, 621–626. [Google Scholar] [CrossRef]
  32. Pontes, B.; Monzo, P.; Gauthier, N.C. Membrane tension: A challenging but universal physical parameter in cell biology. Semin. Cell Dev. Biol. 2017, 71, 30–41. [Google Scholar] [CrossRef] [PubMed]
  33. Price, B.; Dennison, C.; Tschesche, H.; Elliott, E. Neutrophil Tissue Inhibitor of Matrix Metalloproteinases-1 Occurs in Novel Vesicles That Do Not Fuse with the Phagosome. J. Biol. Chem. 2000, 275, 28308–28315. [Google Scholar] [CrossRef] [PubMed]
  34. Riggi, M.; Niewola-Staszkowska, K.; Chiaruttini, N.; Colom, A.; Kusmider, B.; Mercier, V.; Soleimanpour, S.; Stahl, M.; Matile, S.; Roux, A.; et al. Decrease in plasma membrane tension triggers PtdIns(4,5)P2 phase separation to inactivate TORC2. Nat. Cell Biol. 2018, 20, 1043–1051. [Google Scholar] [CrossRef] [PubMed]
  35. Shihata, W.A.; Michell, D.L.; Andrews, K.L.; Chin-Dusting, J.P.F. Caveolae: A Role in Endothelial Inflammation and Mechanotransduction? Front. Physiol. 2016, 7, 628. [Google Scholar] [CrossRef]
  36. Sobrero, A.F.; Puccini, A.; Shi, Q.; Grothey, A.; Andrè, T.; Shields, A.F.; Souglakos, I.; Yoshino, T.; Iveson, T.; Ceppi, M.; et al. A new prognostic and predictive tool for shared decision making in stage III colon cancer. Eur. J. Cancer 2020, 138, 182–188. [Google Scholar] [CrossRef]
  37. Song, G.; Xu, S.; Zhang, H.; Wang, Y.; Xiao, C.; Jiang, T.; Wu, L.; Zhang, T.; Sun, X.; Zhong, L.; et al. TIMP1 is a prognostic marker for the progression and metastasis of colon cancer through FAK-PI3K/AKT and MAPK pathway. J. Exp. Clin. Cancer Res. 2016, 35, 148. [Google Scholar] [CrossRef]
  38. Sung, H.; Ferlay, J.; Siegel, R.L.; Laversanne, M.; Soerjomataram, I.; Jemal, A.; Bray, F. Global Cancer Statistics, 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J. Clin. 2021, 71, 209–249. [Google Scholar] [CrossRef]
  39. Thottacherry, J.J.; Kosmalska, A.J.; Kumar, A.; Vishen, A.S.; Elosegui-Artola, A.; Pradhan, S.; Sharma, S.; Singh, P.P.; Guadamillas, M.C.; Chaudhary, N.; et al. Mechanochemical feedback control of dynamin independent endocytosis modulates membrane tension in adherent cells. Nat. Commun. 2018, 9, 4217. [Google Scholar] [CrossRef]
  40. Tsujita, K.; Satow, R.; Asada, S.; Nakamura, Y.; Arnes, L.; Sako, K.; Fujita, Y.; Fukami, K.; Itoh, T. Homeostatic membrane tension constrains cancer cell dissemination by counteracting BAR protein assembly. Nat. Commun. 2021, 12, 5930. [Google Scholar] [CrossRef]
  41. Tsvetkov, P.; Coy, S.; Petrova, B.; Dreishpoon, M.; Verma, A.; Abdusamad, M.; Rossen, J.; Joesch-Cohen, L.; Humeidi, R.; Spangler, R.D.; et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science 2022, 375, 1254–1261. [Google Scholar] [CrossRef]
  42. Warde-Farley, D.; Donaldson, S.L.; Comes, O.; Zuberi, K.; Badrawi, R.; Chao, P.; Franz, M.; Grouios, C.; Kazi, F.; Lopes, C.T.; et al. The GeneMANIA prediction server: Biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. 2010, 38, W214–W220. [Google Scholar] [CrossRef] [PubMed]
  43. Wu, Z.; Lu, Z.; Li, L.; Ma, M.; Long, F.; Wu, R.; Huang, L.; Chou, J.; Yang, K.; Zhang, Y.; et al. Identification and Validation of Ferroptosis-Related LncRNA Signatures as a Novel Prognostic Model for Colon Cancer. Front. Immunol. 2022, 12, 783362. [Google Scholar] [CrossRef] [PubMed]
  44. Yang, X.; Lin, C.; Chen, X.; Li, S.; Li, X.; Xiao, B. Structure deformation and curvature sensing of PIEZO1 in lipid membranes. Nature 2022, 604, 377–383. [Google Scholar] [CrossRef]
  45. Zhang, N.; Peng, F.; Wang, Y.; Yang, L.; Wu, F.; Wang, X.; Ye, C.; Han, B.; He, G. Shikonin induces colorectal carcinoma cells apoptosis and autophagy by targeting galectin-1/JNK signaling axis. Int. J Biol. Sci. 2020, 16, 147–161. [Google Scholar] [CrossRef] [PubMed]
  46. Zhang, X.; Kelaria, S.; Kerstetter, J.; Wang, J. The functional and prognostic implications of regulatory T cells in colorectal carcinoma. J. Gastrointest. Oncol. 2015, 6, 307–313. [Google Scholar] [CrossRef] [PubMed]
  47. Zhou, R.; Zhang, J.; Zeng, D.; Sun, H.; Rong, X.; Shi, M.; Bin, J.; Liao, Y.; Liao, W. Immune cell infiltration as a biomarker for the diagnosis and prognosis of stage I-III colon cancer. Cancer Immunol. Immunother. 2019, 68, 433–442. [Google Scholar] [CrossRef]
  48. Zhu, J.; Kong, W.; Xie, Z. Expression and Prognostic Characteristics of Ferroptosis-Related Genes in Colon Cancer. Int. J. Mol. Sci. 2021, 22, 5652. [Google Scholar] [CrossRef]
Vaccines 10 01562 g001 550
Figure 1. Flowchart of this study. * p < 0.05, ** p < 0.01, *** p < 0.001.
Figure 1. Flowchart of this study. * p < 0.05, ** p < 0.01, *** p < 0.001.
Vaccines 10 01562 g001
Vaccines 10 01562 g002 550
Figure 2. Differentially expressed MTRGs in COAD tissues compared with normal colon tissues and their interactions. (A) Heatmap of the expression of the 23 MTRGs in the tumors and normal tissues of the TCGA dataset. (B) The expression of MTRGs was significantly different in 41 COAD compared with the normal colon pairs. (C) Interaction analysis among the 23 MTRGs. * p < 0.05, ** p < 0.01, *** p < 0.001.
Figure 2. Differentially expressed MTRGs in COAD tissues compared with normal colon tissues and their interactions. (A) Heatmap of the expression of the 23 MTRGs in the tumors and normal tissues of the TCGA dataset. (B) The expression of MTRGs was significantly different in 41 COAD compared with the normal colon pairs. (C) Interaction analysis among the 23 MTRGs. * p < 0.05, ** p < 0.01, *** p < 0.001.
Vaccines 10 01562 g002
Vaccines 10 01562 g003 550
Figure 3. Construction of the prognostic risk model based on the TCGA-COAD cohort: (A) Forest map of 5 MTRGs significantly correlated with OS, identified by univariate cox analysis. (B) Screening of optimal parameters (lambda) in the LASSO regression model based on the TCGA cohort. (C) LASSO coefficient profiles of the 4 MTRGs determined by the optimal lambda. (D) Heatmap of the expression of 4 MTRGs in high- and low-risk groups. (E) Kaplan–Meier curve for the OS of colon cancer patients in the high- and low-risk groups in the TCGA cohort. (F) ROC analysis of the prognostic model for OS and survival status in the TCGA cohort. (G) Scatterplots in the top and bottom illustrate the distribution of the risk score and survival status in the colon cancer patients, respectively. (H,I) Univariate (H) and multivariate (I) Cox regression analyses of the risk score and clinicopathological parameters in the TCGA cohort.
Figure 3. Construction of the prognostic risk model based on the TCGA-COAD cohort: (A) Forest map of 5 MTRGs significantly correlated with OS, identified by univariate cox analysis. (B) Screening of optimal parameters (lambda) in the LASSO regression model based on the TCGA cohort. (C) LASSO coefficient profiles of the 4 MTRGs determined by the optimal lambda. (D) Heatmap of the expression of 4 MTRGs in high- and low-risk groups. (E) Kaplan–Meier curve for the OS of colon cancer patients in the high- and low-risk groups in the TCGA cohort. (F) ROC analysis of the prognostic model for OS and survival status in the TCGA cohort. (G) Scatterplots in the top and bottom illustrate the distribution of the risk score and survival status in the colon cancer patients, respectively. (H,I) Univariate (H) and multivariate (I) Cox regression analyses of the risk score and clinicopathological parameters in the TCGA cohort.
Vaccines 10 01562 g003
Vaccines 10 01562 g004 550
Figure 4. External validation of prognostic risk models: (A) Kaplan–Meier survival analysis of OS between patients with high-risk scores and low-risk scores in the GSE14333 and GSE103479 cohorts. (B) ROC analysis of the prognostic model in the GSE14333 and GSE103479. (C) risk score and survival status in the GSE14333 and GSE103479 cohorts.
Figure 4. External validation of prognostic risk models: (A) Kaplan–Meier survival analysis of OS between patients with high-risk scores and low-risk scores in the GSE14333 and GSE103479 cohorts. (B) ROC analysis of the prognostic model in the GSE14333 and GSE103479. (C) risk score and survival status in the GSE14333 and GSE103479 cohorts.
Vaccines 10 01562 g004
Vaccines 10 01562 g005 550
Figure 5. Somatic mutation and CNV analysis in colon cancer patients: (A,B) Waterfall plots illustrate the gene mutation landscape in high- and low-risk groups. (A) High-risk group; (B) low-risk group. (C) Survival analysis of OS in colon cancer patients between high- and low-TMB groups. (D) Survival analysis of OS in colon cancer patients between high- and low-TMB groups based on risk score. (E) CNV frequency of 23 MTRGs. (F) Genomic position of 23 MTRGs. Bands at the inner circle indicate corresponding expression levels.
Figure 5. Somatic mutation and CNV analysis in colon cancer patients: (A,B) Waterfall plots illustrate the gene mutation landscape in high- and low-risk groups. (A) High-risk group; (B) low-risk group. (C) Survival analysis of OS in colon cancer patients between high- and low-TMB groups. (D) Survival analysis of OS in colon cancer patients between high- and low-TMB groups based on risk score. (E) CNV frequency of 23 MTRGs. (F) Genomic position of 23 MTRGs. Bands at the inner circle indicate corresponding expression levels.
Vaccines 10 01562 g005
Vaccines 10 01562 g006 550
Figure 6. Network regulation and functional enrichment analysis of high-risk groups based on MTRGs prognostic signature. (A) The regulatory network involving 4-MTRGs and twenty potential binding proteins was constructed through GeneMANIA. (B) Correlation analysis of 4 genes. * p < 0.05 (C–H) GSEA showed the significantly enriched Reactome (C), PID (D), Hallmark (E), GOBP (F), KEGG (G), and BioCarta (H) gene sets in the high-risk colon cancer patients.
Figure 6. Network regulation and functional enrichment analysis of high-risk groups based on MTRGs prognostic signature. (A) The regulatory network involving 4-MTRGs and twenty potential binding proteins was constructed through GeneMANIA. (B) Correlation analysis of 4 genes. * p < 0.05 (C–H) GSEA showed the significantly enriched Reactome (C), PID (D), Hallmark (E), GOBP (F), KEGG (G), and BioCarta (H) gene sets in the high-risk colon cancer patients.
Vaccines 10 01562 g006
Vaccines 10 01562 g007 550
Figure 7. The high- and low-risk groups present different immune statuses. (A) The boxplots for the comparison of the 22 immune cells between the high-risk and low-risk groups in colon cancer. (B,C) Immune cells infiltration score (B) and immune-related pathways’ activity (C) in the high- and low-risk groups estimated by ssGSEA. (D) Expression level of the stromal score, ESTIMATE score, immune score, and tumor purity in the high- and low-risk groups. (E) Associations between the risk score and immune cell infiltration levels. (F) Cuproptosis-related genes between the high- and low-risk groups. * p < 0.05, ** p < 0.01, *** p < 0.001.
Figure 7. The high- and low-risk groups present different immune statuses. (A) The boxplots for the comparison of the 22 immune cells between the high-risk and low-risk groups in colon cancer. (B,C) Immune cells infiltration score (B) and immune-related pathways’ activity (C) in the high- and low-risk groups estimated by ssGSEA. (D) Expression level of the stromal score, ESTIMATE score, immune score, and tumor purity in the high- and low-risk groups. (E) Associations between the risk score and immune cell infiltration levels. (F) Cuproptosis-related genes between the high- and low-risk groups. * p < 0.05, ** p < 0.01, *** p < 0.001.
Vaccines 10 01562 g007
Vaccines 10 01562 g008 550
Figure 8. Cell type constitution of the colon normal and adenoma tissues. (A) t-SNE plots of cells from three patients (6 samples). Colors represent cell types. Cells were clustered into 17 sub-clusters based on biological annotation. Each dot represents a single cell. (B) Dot plot of proportion of cells in the respective cluster expressing selected marker genes. Circle size represents the percentage of cells that express the gene, and color represents the average expression value within a cluster. (C) Expression of the 4-MTRGs at the single cell level. (D) Bar plot showing the fraction of cells originating from the 3 normal and 3 tumor samples (left) and the fraction of cells originating from each of the 3 patients (right). (E,F) t-SNE plots of patient origin (E) and sample origin (F).
Figure 8. Cell type constitution of the colon normal and adenoma tissues. (A) t-SNE plots of cells from three patients (6 samples). Colors represent cell types. Cells were clustered into 17 sub-clusters based on biological annotation. Each dot represents a single cell. (B) Dot plot of proportion of cells in the respective cluster expressing selected marker genes. Circle size represents the percentage of cells that express the gene, and color represents the average expression value within a cluster. (C) Expression of the 4-MTRGs at the single cell level. (D) Bar plot showing the fraction of cells originating from the 3 normal and 3 tumor samples (left) and the fraction of cells originating from each of the 3 patients (right). (E,F) t-SNE plots of patient origin (E) and sample origin (F).
Vaccines 10 01562 g008
Vaccines 10 01562 g009 550
Figure 9. Epithelial cell clusters’ and pseudotime analysis. (A,B) t-SNE plot of 3589 epithelial cells, color-coded by their associated cluster (A) or cell classification (B). Note that the cell type was judged to be malignant by inferCNV analysis (Supplementary Figure S2). (C–F) Trajectory reconstruction (C) of all single cells from non-malignant to malignant epithelial cells; states (D), distribution of clusters (E), and cell types (F) on the trajectory are shown. (G) Heatmap showing time-series gene expression (rows) in epithelial cells (columns) for each membrane tension-related gene. (H) Gene expression dynamics of 4-MTRGs are displayed.
Figure 9. Epithelial cell clusters’ and pseudotime analysis. (A,B) t-SNE plot of 3589 epithelial cells, color-coded by their associated cluster (A) or cell classification (B). Note that the cell type was judged to be malignant by inferCNV analysis (Supplementary Figure S2). (C–F) Trajectory reconstruction (C) of all single cells from non-malignant to malignant epithelial cells; states (D), distribution of clusters (E), and cell types (F) on the trajectory are shown. (G) Heatmap showing time-series gene expression (rows) in epithelial cells (columns) for each membrane tension-related gene. (H) Gene expression dynamics of 4-MTRGs are displayed.
Vaccines 10 01562 g009
Vaccines 10 01562 g010 550
Figure 10. Immunohistochemistry of the TIMP1 gene in tumor and normal tissues from the human protein atlas (HPA) database. (A) Protein levels of TIMP1 in tumor tissue. (B) Protein levels of TIMP1 in normal tissue.
Figure 10. Immunohistochemistry of the TIMP1 gene in tumor and normal tissues from the human protein atlas (HPA) database. (A) Protein levels of TIMP1 in tumor tissue. (B) Protein levels of TIMP1 in normal tissue.
Vaccines 10 01562 g010
Vaccines 10 01562 g011 550
Figure 11. Molecular-docking analysis. (A–C) Elesclomol (A), Shikonin (B), and Bryostatin 1 (C) were screened-out by drug sensitivity analysis. Note that drug sensitivity analysis is shown in Supplementary Figure S3. (D–F) 3D structures and binding modes showing the formed hydrogen bonds between the predicted pocket of TIMP1 and Elesclomol (D), Shikonin (E), and Bryostatin_1 (F).
Figure 11. Molecular-docking analysis. (A–C) Elesclomol (A), Shikonin (B), and Bryostatin 1 (C) were screened-out by drug sensitivity analysis. Note that drug sensitivity analysis is shown in Supplementary Figure S3. (D–F) 3D structures and binding modes showing the formed hydrogen bonds between the predicted pocket of TIMP1 and Elesclomol (D), Shikonin (E), and Bryostatin_1 (F).
Vaccines 10 01562 g011
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, J.; Fu, Y.; Zhang, K.; Li, Y. Integration of Bulk and Single-Cell RNA-Seq Data to Construct a Prognostic Model of Membrane Tension-Related Genes for Colon Cancer. Vaccines 2022, 10, 1562. https://doi.org/10.3390/vaccines10091562

AMA Style

Li J, Fu Y, Zhang K, Li Y. Integration of Bulk and Single-Cell RNA-Seq Data to Construct a Prognostic Model of Membrane Tension-Related Genes for Colon Cancer. Vaccines. 2022; 10(9):1562. https://doi.org/10.3390/vaccines10091562

Chicago/Turabian Style

Li, Jiacheng, Yugang Fu, Kehui Zhang, and Yong Li. 2022. "Integration of Bulk and Single-Cell RNA-Seq Data to Construct a Prognostic Model of Membrane Tension-Related Genes for Colon Cancer" Vaccines 10, no. 9: 1562. https://doi.org/10.3390/vaccines10091562

APA Style

Li, J., Fu, Y., Zhang, K., & Li, Y. (2022). Integration of Bulk and Single-Cell RNA-Seq Data to Construct a Prognostic Model of Membrane Tension-Related Genes for Colon Cancer. Vaccines, 10(9), 1562. https://doi.org/10.3390/vaccines10091562

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop