Identification DMPs and Construction DNA Methylation Prognosis Classifier
In TCGA BCa data, there are 21 pairs DNA methylated samples(21 BCa samples and 21 corresponding normal adjacent samples) in total. After removing the outlier samples, 18 pairs DNA methylation samples were selected for further analysis(Figure 2A). Firstly, we applied the R package “CHAMP” to normalize the DNA methylated data and perform the differential analysis. With the criteria of deltaBeta>0.45 and p<0.001, 2285 differentially methylated positions(DMP) were identified(Figure 2B, Supplementary Table S1). Among 415 BCa DNA methylated samples, 4 duplicated, 5 low-quality and 61 NA-containing samples were discarded(Figure 1). Therefore, 345 DNA methylated samples were retained for downstream analysis. The univariate Cox proportional hazard regression model was applied to analyze 2285 differentially methylated positions. As a result, 50 DMPs(p<0.05) were taken into multivariate Cox regression analysis(Supplementary Table S2). Finally, 10 DMPs were applied to construct the DNA methylation prognosis classifier (Supplementary Table S3). Kaplan-Meier survival analysis displayed that BCa patients with high risk had a poor prognosis(p<0.0001, Figure 2C).
KEGG Pathway and GO Term Function Enrichment Analysis
In 345 BCa samples, 341 corresponding RNA sequencing data(169 high risk and 172 low risk) were retrieved from the TCGA database. The differential KEGG pathways between high and low risk groups were neuroactive ligand receptor interaction, chemical carcinogenesis receptor activation, melanoma, thyroid hormone synthesis, primary immunodeficiency, complement and coagulation cascades, tryptophan metabolism and carbohydrate digestion and absorption(Figure 3A). The top three differential cellular components(CC) were external side of plasma membrane, blood microparticle and lamellar body(Figure 3B). Humoral immune response, antimicrobial humoral response and negative regulation of response to wounding were the top three differential biological process(Figure 3C). The differential molecular functions(MF) between groups were shown in Figure 3D.
Evaluation of Tumor-infiltrating Immune Cells
In order to understand the immune microenvironment of BCa tissues, we evaluated 22 tumor infiltrating immune cells of 341 tumor samples by CIBERSORT. With the CIBERSORT P<0.05, the proportion of 22 tumor infiltrating immune cells of 173 patients(88 high risk and 85 low risk) were displayed in Figure 4A. In order to further determine the differences between different groups, we compared 22 tumor infiltrating immune cells between high-risk and low-risk patients(Figure 4B). The results showed that naive B cells, plasma cells, CD8+ T cells and T cell regulatory(Tregs) infiltrated more in low-risk patients.
The Correlation Between DNA Methylation Positions and Host Gene Expression
For the purpose of identifying the correlation between DNA methylation position and host gene expression, we calculated the Pearson correlation coefficient of ten DNA methylation positions and corresponding host gene expression. The results showed that the methylation of cg06551997 and cg07568841 negatively correlated with the expression of SLFN12L and ZNRF2, respectively(Figure 5A, 5B). On the contrary, the expression of COL23A1, NR2E1 and TFAP2B were positively correlated with the methylation level of cg08560734, cg06433023 and cg22282405(Figure 5C, 5C, 5D), separately.
Differences in Sensitivity to Immune/Chemotherapy Between High- and Low-risk BCa
Previous studies have confirmed that CTLA-4, PD-1 and PD-L1 inhibitors show significant anti-tumor responses in the treatment of advanced bladder cancer(19). To identify the sensitivity difference to immune checkpoint inhibitors between high and low risk BCa patients, the TIDE algorithm was applied to evaluate the response of 341 patients to immunotherapy. The results showed that low-risk patients (37.21%, 64/ 172) may be more likely to respond to immunotherapy than high-risk patients (23.08%, 39/169) (P = 0.004753, Figure 6A). Besides, subclass mapping was also applied to evaluate the sensitivity of two BCa subtypes to immunotherapy(CTLA-4 and PD-1). As shown in Figure 6B, PD-1 inhibitors may be more effective for low-risk subtype.
As the important method for BCa therapy, chemotherapy plays a key role in improving the prognosis of patients with BCa. Therefore, we applied the cell line data in GDSC and the transcriptome data in TCGA to predict the IC50 of each sample through ridge regression model. The results indicated that low risk patients with BCa were more sensitive to many chemotherapeutic drugs including methotrexate(Figure 7).
For the precision treatment of BCa, we applied CTRP and PRISM- derived drug response data to identify drug candidates with higher drug sensitivity in patients with high or low risk scores. Firstly, we performed differential drug response analysis between the groups with high-risk scores (top decile) and low risk scores (bottom decile) to identify compounds in the high risk group with lower estimated AUC values or compounds in the low risk group with lower estimated AUC values(log2FC>0.10). And then, potential compounds were selected by evaluating the spearman correlation coefficient between the AUC value and the risk score (Spearman’s r <−0.2 or >0.2 for CTRP and PRISM). The results displayed that high-risk patients were more sensitive to poziotinib. However, mitoxantrone and PI−103 may be more effective for patients with low risk scores(Supplementary Figure 1).
Differences in Single Nucleotide Variants and Copy Number Variation Between High- and Low-risk Bladder Cancer
The mutation is involved in the occurrence and progression of tumors, and affect the prognosis of patients. To explore the SNV, we separately described the single nucleotide variant profiles of the high- and low-risk patients(Figure 8). As shown in Figure 9A, the top ten mutation genes were TP53(15%), TTN(14%), KDM6A(10%), KMT2D(9%), MUC16(9%), PIK3CA(9%), ARID1A(7%), RB1(7%), MACF1(6%) and CDKN1A(5%) in high risk patients. In low risk patients(Figure 9B), the top ten mutation genes were TTN(15%), TP53(11%), KMT2D(10%), MUC16(9%), SYNE1(9%), FGFR3(8%), RB1(7%), KDM6A(7%), FAT4(6%) and ARID1A(5%). In addition, we found there was no statistical difference in TMB between high-risk and low-risk BCa(Supplementary Figure 2). As for the copy number variation, variation occurs on most chromosomes except for chromosomes 10, 15 and Y in high-risk BCa patients(Figure 10).