Analysis Based on TCGA Data and Single-cell Data, taking TRPM4 as an Example

As we all know, bladder cancer (BLCA) is one of the most aggressive and metastatic malignant tumors in the world, and there are still issues with the current biomarkers for bladder cancer, such as inaccuracy and so on¹. To identify the prognosis of BLCA patients and predict the outcomes of BLCA patients, the search for biomarkers for bladder cancer and the establishment of prognostic models are of great importance. Although people have developed some research methods for biomarkers, most of these methods are currently limited to transcriptomics, which inevitably leads to heterogeneity in samples². Moreover, studying transcriptomics data alone often fails to investigate the roles played by different cell subpopulations, as well as the functions of genes in different pathways at the single-cell level, which makes past research less precise and effective³. Considering that single-cell analysis might be challenging for beginners, an online platform has been introduced for single-cell analysis to help them learn the skills quickly⁴. Last but not least, even if a single gene is found to be a suitable biomarker through transcriptome analysis and single-cell analysis, there is no guarantee that the biomarker will be applicable across all cohorts, so it is necessary to construct prognostic models associated with the biomarker to make the conclusions more universally applicable⁵. The 101 machine learning algorithm refers to the construction of 101 prognosis models using a combination of 10 different machine learning algorithms, with the goal of identifying the optimal prognosis model. The inclusion of algorithms such as random forest, XGBoost, and SVM, which are more capable of handling complex data and exhibit greater stability, has resulted in this combined algorithm demonstrating remarkable stability. In order to make this prognostic model more accurate and relevant to the biomarker, a correlation analysis is first conducted between all the genes and the biomarker, then select around 20-30 genes based on requirements, and subsequently use a 101 machine learning algorithm to build the prognostic model, followed by a series of analyses to finalize the results.

Compared to the biomarkers predicted by simple transcriptomics analysis of the past⁶, the biomarkers of single-cell and transcriptomics analysis allow for an unbiased breakdown of tissues or samples into their basic cellular components. It enables the clear differentiation of different cell types (such as T cells, B cells, and macrophages) and the discovery of new subpopulations (such as exhausted T cells and regulatory T cells), as well as the capture of continuous dynamic processes (such as cell differentiation trajectories)⁴. It is akin to arranging fruits and milk separately, allowing for a clear visualization of each component and its quantity. Compared with single-cohort Cox models⁷, the prognostic model constructed using 101 machine learning is also more accurate and scientifically sound, as it can automatically learn about the complex, non-linear interactions between variables from the data. For instance, the impact of a specific genetic mutation might only be significant in patients of certain ages and tumor sizes. Machine learning models such as random forests and neural networks can automatically capture these high-order interactions without the need for manual specification⁸. The signaling key applicability constraints that the dataset used for analysis must include the vast majority of genes, with the number of genes not being too low, and the number of genes used for constructing prognostic models not being too high, generally maintained around 20-30, being optimal. The quality control standard of the single cell set should be greater than 1000, the UMI count per cell should be greater than 1000, and the gene number per cell should be greater than 500⁹.

Here, a step-by-step approach is provided for the identification of novel biomarkers from public transcriptomic datasets and single-cell datasets, taking the role of TRPM4 in BLCA as an example. Multiple research methods are employed and diverse datasets -- including the Cancer Genome Atlas-Bladder Cancer (TCGA-BLCA) dataset, single-cell dataset GSE145281, and BLCA datasets GSE32894 and GSE31684 -- to enhance the precision and applicability of biomarkers in oncology from multiple perspectives and to advance the study of TRPM4 in BLCA. The TCGA-BLCA dataset was obtained from the University of California, Santa Cruz Xena (UCSC Xena) website. GSE32894 and GSE31684 were downloaded from the Gene Expression Omnibus (GEO), and single-cell data GSE145281 was sourced from Tumor Immune Single-Cell Hub 2 (TISCH2). Additionally, data on BLCA molecular subtypes and treatment responses were extracted from supplementary spreadsheet files of a relevant article. Bioinformatics analysis, single-cell analysis, and machine learning are then conducted, establishing an integrated single-gene research methodology.

NOTE: All the codes used in this article can be found on the website https://github.com/YaoGeng-nmu/Analysis-of-TRPM4-based-on-TCGA-data-and-single-cell-data-in-BLCA/blob/main/code.

1. Preparation of transcriptomics data

1. Preparation of the TCGA-BLCA dataset
   1. Download all the TCGA datasets from the website UCSC Xena (https://xenabrowser.net/datapages/)¹⁰. The datasets downloaded are pre-processed, eliminating the need for additional work such as gene annotation¹¹.
   2. Select the TCGA dataset and then click on the TCGA-BLCA dataset. Then, go to the page of the TCGA-BLCA dataset, click on Gene Expression RNA-seq, and download the clinical data and gene expression profiles from the TCGA-BLCA dataset. In this dataset, ensure that there are 400 BLCA patients' mRNA expression data.
   3. To differentiate tumor samples from adjacent normal tissue samples, separate the samples with the endings 01 and 11 as samples with 01 indicate tumor samples, while the samples with 11 indicate cancer-adjacent tissues (normal samples).
2. Preparation of GSE32894 and GSE31684
   1. Click on the GEO Website (https://www.ncbi.nlm.nih.gov) in order to download BLCA dataset GSE32894 and GSE31684^12,13.
   2. After clicking on the Respective Entries in the upper-right corner, make sure to move to the pages for datasets GSE32894 and GSE31684.
   3. Next, download the respective mRNA data, clinical data, and survival data for GSE32894. Click on the http button of the GSE32894 dataset. Also, download the GPL6947 platform data using the platform button. Click on the Series Matrix File(s) button and download GSE32894 after clicking on the GSE32894_series_matrix.txt.gz button when downloading mRNA data, clinical data, and survival data for GSE32894.
   4. Download the respective mRNA data, clinical data, and survival data for GSE31684. Click on the http button of the GSE31684 dataset. Also, download the GPL570 platform data using the platform button. Click on the Series Matrix File(s) button and download GSE31684 after clicking on the GSE31684_series_matrix.txt.gz button when downloading mRNA data, clinical data, and survival data for GSE31684.
   5. In the filtering step, remove meaningless genes, such as those with an expression level of zero or close to zero. Employ the ComBat method, which is based on the sva package¹⁴, for batch correction.

2. Preparation of single-cell analysis

1. Preparation of GSE145281
   1. To conduct research at the single-cell level, find the BLCA single-cell dataset, considering that transcriptome data alone is insufficient to determine which cell cluster expresses TRPM4 most significantly.
   2. Here, go to TISCH2 (http://tisch.compbio.cn/), an online website for single-cell analysis of tumors. 
   3. According to experimental needs, select the Cancer Type needed, using BLCA as an example in this case.
   4. Click on the BLCA button, choosing the dataset called GSE145281.

3. Preparation of BLCA's molecular subtypes and response to therapeutic choices data

1. Download the data for BLCA's molecular subtypes and response to therapeutic choices from the article referenced here¹⁵.
2. Open the PubMed website (https://pubmed.ncbi.nlm.nih.gov/). Search for the article and click on Supplementary Tables. There is a total of 18 files in this compressed file package.

4. BLCA immune microenvironment analysis of TRPM4

1. Heatmap plot of the expression of 133 immunomodulators in different TRPM4 groups
   1. According to the expression of TRPM4, divide all BLCA samples into the high-TRPM4 group and the low-TRPM4 group.
   2. Copy all the sample names and TRPM4 expression data separately, then sort them in descending order of TRPM4 expression, and categorize the first half of the samples into the High-TRPM4 group and the latter half into the Low-TRPM4 Group in the Group column according to TRPM4 expression. All TCGA-BLCA patients' IDs are listed in Supplementary Table 1.
   3. Install R packages required for the transcriptomic data analysis. Make slight modifications to the provided code, such as changing the file location, and then run the code. Carefully check if the gene IDs and the IDs of the immunomodulators are correctly matched, ensuring that there are no spaces or other inconspicuous factors after the gene IDs, as these are the primary causes of failed heat map generation.
2. Stromal and immunological analysis of tumor samples using the ESTIMATE R package
   1. Make sure that all R packages required for the transcriptomic data analysis are installed. According to the code requirements, name the required data with the corresponding names, and then click the Run button.
      ​NOTE: It is important to remember that this website also offers online analysis capabilities. While it is possible to obtain the desired results by visiting the website (https://bioinformatics.mdanderson.org/estimate/rpackage.html), this method is not suitable for analyzing data that has been self-collected.
   2. Click on the Disease button and choose the Bladder Urothelial Carcinoma option. Then, select the RNA-seq-V2 platform button. Download all samples' stromal score, immune score, and estimate score.
3. Violin plot of the expression between different TRPM4 groups
   1. Install all R packages required for transcriptomic data analysis. All the R packages needed are listed in Supplementary Table 2. According to the code requirements, name the required data with the corresponding names, and then click the Run button.
   2. Make sure to carefully verify that the gene IDs are entered correctly, ensuring there are no spaces or other inconspicuous factors at the end of the gene IDs, as these are the primary causes of failures in the creation of violin plots.
4. Triangle heatmap plot between different TRPM4 groups
   1. Install all R packages required for the transcriptomic data analysis. Make sure to prepare the expression data of several immune checkpoint genes, such as TIGIT and CD80. Open the TCGA-BLCA transcriptome data and then search through each entry individually, given that the number of immune checkpoint genes is not extensive.
   2. According to the code requirements, name the required data with the corresponding names, and then click the Run button.
5. Correlation plot between different genes and TRPM4
   1. Install all R packages required for the transcriptomic data analysis. Fill in TRPM4 as the first gene and add the gene being researched as the second gene, which is CD3E in this case.
   2. After clicking the Run button, make sure to change CD3E to other genes that need to be analyzed, such as CTLA4.

5. Single-cell analysis of TRPM4 in the GSE145281 dataset

1. BLCA GSE145281 dataset analysis using single-cell analysis
   1. Enter the website and click on BLCA. Single-cell analysis allows for a more in-depth exploration of which cell cluster exhibits high expression of TRPM4 compared to transcriptomics analysis. Make sure that the single-cell analysis data is collected from GSE130001 and check against the quality control standard that the cell number per dataset should be greater than 1000, the UMI count per cell should be greater than 1000, and the gene number per cell should be greater than 500. The differentially expressed genes of each cluster compared to all other cells are identified based on the log-transformed fold change (|logFC| >= 0.25), and the clustering resolution should be 0.5.
2. Cell annotation for different cell clusters
   1. Make sure that all cell clusters are annotated. For example, the B cells can be annotated by CD19, MS4A1, CD27, IGHD, IGHM, TCL1A, and FCRL5. All the genes used for cell annotation can be found in Supplementary Table 3.
   2. Then, choose GSE145281 and download all results of TRPM4. Click on the Overall button and download the UMAP plot of GSE145281 of different clusters and cell types.
   3. Click on the Gene button and choose the single gene that is going to be analyzed. Click on the GSEA button and observe the functions of this single gene across various pathways, such as the KEGG pathway and Hallmark pathway^16,17.
   4. Click on the CCI button and find the cell chat plot and cell chat bubble plot of the cluster highly expressing a single gene.
   5. Click on the TF enrichment button to see the enriched Transcript Factor for each cluster.

6. 101 Machine learning to construct a prognostic model relating to TRPM4

1. Install all R packages required for transcriptomic data analysis. Considering that a single machine learning approach might not always be the optimal method for predicting prognosis, a prognosis model can be constructed using 101 machine learning techniques. Conduct a correlation analysis between all the genes in the TCGA-BLCA cohort and GSE32894 and GSE31684 cohorts and TRPM4.
2. Based on the number of genes strongly associated with TRPM4, select all genes with a correlation greater than or equal to 0.6 or less than or equal to -0.6 with TRPM4. Select all genes strongly associated with TRPM4 in all cohorts.
3. Combine all the data from GSE32894 and GSE31684 cohorts to form a validation set and use TCGA-BLCA as the training set.
4. According to the code requirements, name the required data with the corresponding names, and then click the Run button.

It is well-known that the core benefit of using transcriptomics to analyze the immune microenvironment for individual genes is that it allows for a direct correlation between individual genes and the dynamic interactions with immune cells and molecules at the level of gene expression, enabling the direct observation of immune characteristics and the identification of differences in the proportion of immune cells infiltrating between groups with high/low expression of individual genes or changes in the expression of immune checkpoint molecules18. The representative results of bioinformatics analysis are those related to the analysis of the immune microenvironment and immune molecules, such as violin plots, etc. The representative results of single-cell analysis are the target gene expression levels for each cell cluster and the pathway heatmap analysis for each cell cluster. The representative results of 101 machine learning models are the heatmaps of the algorithms used to build prognostic models.

First of all, the heatmap plot shows that 133 immunomodulators are more expressed in the low-TRPM4 group (Figure 1A). According to the calculation of ESTIMATE R package, it's obvious that stromalscore, immunescore, and estimatescore are all lower in the high-TRPM4 group, thus indicating that in samples with high TRPM4 expression, immune infiltration and interstitial infiltration are often insufficient (Figure 1B-D). In addition, the triangular correlation heatmap also shows that TRPM4 is negatively correlated with several immune checkpoints, such as VTCN1, CD274, PDCD1, CTLA4, HAVCR2, TIGIT, IDO1, CD80, CD86, LAIR1, PVR, CD200R1, CD200, LGALS3, CEACAM1, BTLA, ADORA2A, KIR3DL1, KLRC1, CD276, CD47, KLRD1 and LAG3 (Figure 1E). However, we completed the correlation plot between TRPM4 and T cell exhaustion genes, including CTLA4, LAG3, HAVCR2, CD3E, TIGIT, and PDCD1 (Figure 1F).

Although transcriptomics has clearly demonstrated the association between TRPM4 and TME, analysis at this level cannot distinguish which type or subpopulation of cells is responsible for these phenotypes. To precisely pinpoint the cellular origin of the immune phenotype and further clarify whether a specific type of immune cell has undergone a functional change, single-cell transcriptomics analysis is necessary. By analyzing the gene expression characteristics of individual cells, the immune phenotype can be accurately traced back to specific cell subpopulations. Additionally, GSE145281 has been introduced, a single-cell dataset of bladder cancer immune cells, to investigate the role of TRPM4 in the immune cell pathway. All 14462 immune cells are divided into 17 clusters (Figure 2A). All 14462 cells are divided into B cells, CD4 T cells, CD8 T cells, Mono, and NK cells (Figure 2B). Among all 14462 immune cells, 8465 of them are Mono cells, 2677 of them are CD8 T cells, 1771 are CD4 T cells, 566 of them are B cells, and 995 of them are NK cells (Figure 2C). We also describe the distribution of TRPM4 (Figure 2D). The results show that TRPM4 is predominantly expressed in Monocytes. The violin plot shows that TRPM4 is mostly expressed in the C1 and C2 clusters (Figure 2E). The results in Figure 2F also demonstrated the previous result. In addition, the cell chat analysis shows that cluster C1 is associated with cluster C0 and cluster C13, while cluster C2 is associated with cluster C0 and cluster C11 most (Figure 2G-H). The role of TRPM4 is also explored within the Hallmark gene set and depicted using a heatmap. The heatmap result shows that the C1 cluster is mostly enriched in the TNF-α pathway. This suggests that TRPM4 is likely involved in its effects through the TNF-α pathway (Figure 2I). The C2 cluster is mostly enriched in the xenobiotic metabolism pathway.

The single-cell signature analysis also shows that the two channels mentioned above are distributed in relation to TRPM4 (Figure 3A-B). The cell chat bubble plots are also made to describe the cell chat of the C1 and C2 clusters (Figure 3C-D). Lastly, the TF enrichment heatmap plot is generated (Figure 3E). The rank plots indicate that SPI1 could be a transcription factor that regulates gene expression in each cell population (Figure 3F-G). Through single-cell analysis, cell subpopulations have been precisely identified as associated with TRPM4 and clarified the potential mechanisms by which they influence disease progression. These findings provide a solid theoretical foundation for understanding the nature of the disease, but they have not yet been translated into tools that can directly guide clinical decision-making. To bring these insights into practice, further integration of machine learning is required to build prognostic models. This involves collecting data on genes that have a strong correlation with TRPM4 and ultimately developing a model that can reliably predict patient survival. All the genes highly related to TRPM4 are listed in the Supplementary Table 4. Last but not least, to further predict the prognosis of BLCA patients, 101 machine learning models are employed to better describe all patients' prognosis.

In order to avoid the randomness and inaccuracy brought about by relying on a single data source, we utilize multiple BLCA datasets, such as GSE32894 and GSE31684. The results showed that Coxboost+RSF could predict the BLCA patients' prognosis best, and the c-index of this prognostic model can be 0.879 (Figure 4).

In addition to TRPM4, the ESTIMATE R package has also been applied to analyze the stromalscore, immunescore, and ESTIMATEScore within different SDCBP expression groups (Supplementary Figure 1A-C). The triangular correlation heatmap is also made between SDCBP and immune checkpoints to illustrate how to adapt the protocol (Supplementary Figure 1D). The qPCR experiment and the results align with our bioinformatics analysis. Human BLCA cell line T24 and normal bladder cell line SVHUC1 were used. As shown in the result graph, cells with high TRPM4 expression also exhibit high levels of TRPM4 (Supplementary Figure 2).

Figure 1: Bioinformatics analysis of TRPM4 in shaping tumor microenvironment. (A) Heatmap plot of 133 immunomodulators in different TRPM4 groups. (B-D) Violin plots of stromalscore, immunescore and estimatescore in different TRPM4 groups. ****p ≤ 0.0001. (E) Triangular correlation heatmap between TRPM4 and immune checkpoints. (F) Correlation plots between T cell exhaustion genes and TRPM4. Please click here to view a larger version of this figure.

Figure 2: Single-cell analysis of TRPM4 using the GSE145281 dataset. (A) Different clusters of GSE145281. (B) Different immune cells of GSE145281. (C) Pie plot of GSE145281. (D) The distribution of TRPM4. (E-F) Violin plot of TRPM4 in different clusters and immune cells. (G-H) Cell chat analysis of C1 and C2 clusters. (I) Hallmark pathway analysis of different clusters. Please click here to view a larger version of this figure.

Figure 3: Single-cell analysis of TRPM4 using the GSE145281 dataset. (A-B) Single-cell signature analysis using the GSE145281 dataset. (C-D) Bubble plot of C1 and C2 clusters. (E-G) TF enrichment heatmap plots. Please click here to view a larger version of this figure.

Figure 4: 101 machine learning heatmap plot of prognostic model relating to TRPM4. The heatmap plot of different combinations of machine learning algorithms. Please click here to view a larger version of this figure.

Supplementary Figure 1: Bioinformatics analysis of SDCBP in shaping tumor microenvironment. (A-C) Violin plots of stromalscore, immunescore and estimatescore in different SDCBP groups. ****p ≤ 0.0001. (D) Triangular correlation heatmap between SDCBP and immune checkpoints. Please click here to download this File.

Supplementary Figure 2: Gene expression data in cell lines. Please click here to download this File.

Supplementary Table 1: All TCGA-BLCA patients' ID. Please click here to download this File.

Supplementary Table 2: Different versions of these packages used in this article. Please click here to download this File.

Supplementary Table 3: All marker genes used for TISCH2. Please click here to download this File.

Supplementary Table 4: All the genes correlated with TRPM4. Please click here to download this File.

In previous studies, bioinformatics analysis of single genes often encountered issues such as superficiality, inaccuracy, and limited applicability¹⁹. Moreover, previous research on single genes has typically been confined to transcriptome data. Moreover, when focusing solely on transcriptome data, issues such as sample heterogeneity and randomness may arise, making it difficult to identify universal patterns²⁰. Therefore, to address the aforementioned issues, single-cell data have been introduced to better understand the crucial role of individual genes at the pathway level. The key step in this article is that an online website has been introduced called TISCH2 (http://tisch.compbio.cn/) for single-cell analysis⁴. The standard used for quality control is that the cell number per dataset is over 1000, the UMI count per cell is over 1000, and the gene number per cell is over 500. In addition, the 101 machine learning is also the key step used to find the best prognostic model to predict all patients' prognosis and clinical outcomes²¹. Last but not least, do remember to apply several bioinformatics analysis plots to describe the results of a single gene.

Of course, it is inevitable to encounter issues during the process of creating diagrams and performing analyses. If it is discovered that the genes output from the heatmap are few or almost nonexistent, it is essential to carefully check if there are any spaces after the gene IDs, as these spaces can significantly affect the output of the results. If the code frequently throws errors, a significant portion of the issue may stem from personal modifications that are not properly implemented, such as setting the path to a different location on the user's own computer. To resolve this, simply set the path to the location on the user's own machine. Finally, when building a prognostic model using 101 machine learning, patience is crucial. This process can take a considerable amount of time, lasting anywhere from half an hour to an hour, and it is important to ensure that all necessary R packages have been downloaded beforehand.

However, there are still several limitations of this method. Due to time limitations, we only validated this conclusion in BLCA through bioinformatics analysis instead of other cancer types. Moreover, the conclusion may not be applicable to other cancers, necessitating more in-depth and extensive research to provide a comprehensive explanation. The protocols can be richer and more fulfilling by adding more kinds of bioinformatics analysis. In addition, the bioinformatics analysis results should be integrated with experiments; for example, if we find a connection between TRPM4 and TNF-α pathways in the analysis, we can conduct PCR experiments to verify this connection, and the bioinformatics analysis results can serve as a guide for the direction of experimental verification.

In this method article, we also provide convenient download methods for a large amount of data, as well as introductions to single-cell online analysis websites, which greatly reduces the analysis burden for beginners and offers more diverse options for more people. There are also various bioinformatics codes that allow users to experiment with creating their own bioinformatics plots, such as a violin plot, a heatmap plot, a correlation plot, and a Triangular correlation heatmap plot²². Another significant advantage of the research method is that even if this gene is replaced, other single genes can still be studied in a similar manner, making the research findings universally applicable. In the future, we will increase the sample size by collecting data from multiple hospitals and extend the study to additional cancer types. This expansion aims to assess if findings hold true across diverse cancers. Concurrently, a robust prognostic model by integrating TRPM4 with established biomarkers like TMN stage, common in BLCA. This model is designed to enhance patient outcome prediction. Furthermore, machine learning techniques will be employed to build this model, focusing on improving its interpretability and predictive accuracy, which constitutes primary future work.