We downloaded the DNA methylation data measured by Illumina Infinium HumanMethylation450 BeadChip (β values), together with the gene expression data measured by RNA-seq (read counts), the somatic mutation data (MC3 public version), and the overall survival (OS) information of 32 tumor types from the Cancer Genome Atlas [TCGA, downloaded from the GDC data portal (portal.gdc.cancer.gov) in October 2020]. There were 7,131 cases in total. The name, abbreviations, and number of cases of each tumor type are listed in Supplementary Table S1 .

The tumor mutation burden (TMB) was calculated as described before with some small modifications (23). In each case, mutations annotated as “Frame Shift Ins”, “Nonsense Mutation”, “Frame Shift Del”, and “Splice Site” were defined as truncation mutations, while those annotated as “Missense Mutation”, “In Frame Del”, “In Frame Ins”, and “Nonstop Mutation” were defined as non-truncation mutations. TMB value was then calculated as

The distribution of calculated TMBs is shown in Supplementary Figure S1 .

The microsatellite instability (MSI) state for each case was calculated by the MSIpred package as stable or unstable (24).

The cases that simultaneously satisfy the following two criteria are defined as responsive to ICI treatment (positive), and the rest were defined as non-responsive to the treatment (negative). First, the TMB value of the case should be higher than the upper quartile of TMB values in all cases ignoring tumor types. Second, the transforming growth factor β (TGF-β) score [defined in (25) as the weighted average of normalized expression of genes in the TGF-β signaling pathway summarized in (26), taking the regularity direction between genes as weights] should be smaller than the median of the values in all cases ignoring tumor types. It has been shown that both the TMB and the TGF-β score were fine and independent measurements of the responsiveness to ICI treatment (27–31), and the combined measurement has also been successfully used in other pan-cancer studies (32). It is also a practical consideration to define the responsiveness indirectly since there is no large pan-cancer ICI treatment cohort with DNA methylation levels measured for building the models to the best of our knowledge. The number of cases labeled as positive in each tumor type is shown in Supplementary Figure S2 . It was noticeable that this definition led to severe class imbalance, with 9.86% cases labeled as positive. This issue was solved by the random oversampling scheme in the model-building step described below.

The NSCLC validation cohort data were retrieved from literature (13) and (33). The clinical information was downloaded from the Supplementary Tables in corresponding literature. The DNA methylation profiles were downloaded from the associated datasets (GSE119144 and GSE126043) in Gene Expression Omnibus (GEO). Missing values were imputed as described in the raw pieces of literature ( Supplementary Table S2 ).

To perform differential analysis of methylation data, each β value was first transformed into M-value by the transformation, M = log ( β 1 − β ) . Then differential analysis was carried out using limma package as common practices (34). Probes with absolute logFC > 1 and adjusted p-value less than 0.01 were determined as differential methylated probes (DMP).

RNA-seq read counts for each gene were directly inputted into the DESeq2 package for differential expression analysis (35). Genes with log ₂ fold change greater than 1 and adjusted p-value less than 0.01 were determined as differentially expressed.

This study was focused on predicting the responsiveness of immune checkpoint inhibitor (ICI) at the pan-cancer level. To achieve this goal, the 7,131 TCGA cases from 32 tumor types with both DNA methylation (measured by Illumina Infinium HumanMethylation450 BeadChip, about 480K probes measured in total) and gene expression (measured by RNA-seq) data were downloaded (in October 2020). Cases with high tumor mutation burden (TMB, see the distributions of TMBs in Supplementary Figure S1 ) and low TGF-β expression were defined as responsive to ICI (positive) and the others as without responsiveness (negative) as previously suggested (Methods, see Supplementary Figure S2A for the number of cases marked as positive in each tumor type) (27–29, 32) since large pan-cancer tumor immunotherapy cohort with DNA methylation level measured was lacking. A series of well-known biomarkers, such as the concentration of CD8⁺ T cells, the expression level of PD-1 and CTLA4 genes, were significantly separated between positive and negative cases (all with p < 0.01, Supplementary Figures S2B–D ), elucidated the rationality of our classification. It was reported that the clustering based on immune infiltration analysis using DNA methylation data in most tumor types grouped cases into two clusters, and the clusters were shown to be related to the outcome of immunotherapy (22). Based on this background knowledge, we designed an efficient feature selection scheme to select 1,495 methylation probes from the total 480K candidate methylation probes, which distinguished cases with ICI responsiveness and those without (see Methods for detail). In the first step of the feature selection, the analysis was restricted to tumor types with more than 5 cases marked as positive (13 types in total). Tumor cases were clustered into two clusters based on the immune infiltration profiles derived by the DNA methylation profiles. The size of the two clusters for each used tumor type is shown in Supplementary Table S1 . Then we tested differences between the two clusters in each tumor type for three indicators for the responsiveness of ICI treatment, overall survival (OS), tumor mutation burden (TMB), and the PD-L1 gene expression level. After applying the Benjamini-Hochberg FDR correction scheme, we obtained 6, 7, and 7 tumor types with significant differences of these indicators, respectively (FDR<0.01, see Supplementary Figures S3A–C for each indicator Supplementary Figure S3D for the overlap of tumor types selected for the three indicators). For each indicator, we selected those methylation probes which differentially methylated with the same direction in at least half of these tumor types with significant differences in the indicators between groups. We got 890, 862, and 534 feature probes for the indicators OS, TMB, and PD-L1 expression, respectively ( Supplementary Figure S3E ). So far, 11 tumor types and more than 85% of positive cases are covered ( Supplementary Figures S2A and S3D ). We then merged the feature probes selected for each indicator to obtain 1,495 final probes (the selected features are listed in Supplementary Table S3 ). The 1,152 probes in the signature used in the immune infiltration analysis (22) were added and yielded a final feature set with 2,546 probes.

The efficiency of the 2,546 features was revealed by the distinct methylation profiles between cases with ICI responsiveness and those without by visual inspection of the methylation pattern of all cases clustered by the methylation pattern ( Figure 2A ). This pattern was also illustrated by the mutual dependence between the responsiveness of neighboring cases along the dendrogram of the hierarchical clustering (p < 10^-52, χ ² test). Moreover, compared to randomly selected, equal-sized control feature sets, clustering based on methylation pattern of selected probes is much better at separating cases with and without ICI responsiveness. In 100 times of such comparisons, the selected features outperformed 92 and 99 random controls when the separation was measured using F ₁ score and Matthews correlation coefficient (MCC) score, comparing the responsiveness and cluster labels when clustering the cases into two clusters with the methylation patterns ( Figure 2B ). Also, a large portion of selected probes were under differential methylation between cases with and without responsiveness. When performing the differential methylation analysis among all 32 tumor types, there were 8.76% of selected probes with differential methylation, compared with 5.22 to 7.62% in the 100 randomly selected control features under the threshold of adjusted p-value less than 0.01 ( Figure 2C ). At last, the GO and KEGG enrichment analysis of genes in which the selected probes located showed enrichment of terms related to immune activity and immunotherapy outcome, such as “T cell activation”, “adaptive immune response”, and “regulation of lymphocyte activation” in the GO enrichment analysis, and “primary immunodeficiency” and “Th1 and Th2 cell differentiation” in KEGG enrichment analysis ( Figure 2D ).

After the prediction potentials for the selected methylation sites for the ICI treatment responsiveness were elucidated, a prediction model was built and tested from various aspects. Here we tested a series of commonly used machine learning models, including the support vector machine (SVM), logistic regression with L ₁ regularization (LR), the random forest (RF), and k-nearest neighbor classifier (kNN). The hyper-parameters of each model were tuned by 5-fold cross-validation. In each model training, the positive cases were oversampled to the size of negative cases to deal with the severe class imbalance (37). The performances of the optimized models were assessed by randomly splitting the raw data into 80% training and 20% test datasets for 100 times. For comparison, we also added a naive predictor based on the clustering results derived from the immune infiltration analysis based on the DNA methylation profiles. In each tumor type, the naive predictor declared the cluster of cases with higher average TMB as positive and the other cluster as negative.

Among the four assessed models, SVM got the highest performance. The fine-tuned SVM model in the 5-fold cross validation (with λ = 13) outperformed the other three models (p = 0.041, 3.83×10^-28 and 7.33×10^-32 for LR, RF, and kNN, respective, paired t test) when the performance was measured using F ₁ score, and Matthews correlation coefficient (MCC, p = 0.011, 1.30×10^-25 and 8.11×10^-33, Figure 3A ). All the machine learning–based models significantly outperformed the naive predictor in both measurements (all with p < 10^-30, Figure 3A ).

The prediction power of the model cannot be achieved by randomly selected probes. To show this, we randomly selected methylation probe sets of the same size as the selected probe set for 100 times. The performances were measured by repeatedly training models on 80% randomly chosen samples and tested on the other 20% 100 times along all the regularization parameter λ values of the SVM model searched in the cross-validation step (1 to 20). The performances of models based on the selected probes were consistently and significantly higher than those of random controls when the performance was measured by F ₁ score, while the same conclusion held under most λs (17 out of 20) when the performance was measured using MCC score ( Figure 3B ).

The probes selected by the four indicators were all important to the model performance. This was shown by retraining and evaluating the SVM model with probes selected by only one indicator. The performances were all significantly decreased compared to the full model no matter measured by the F₁ score or the MCC score ( Figures 3A, C , all with p < 10^-4, paired t tests).

To exclude the possibility that the model performances were due to overfitting, we tested our prediction model with 100 randomly permuted datasets and expected a sharp shrink of model performances in these permutated datasets. We trained and evaluated SVM models with the same super-parameters in each permutated dataset as described above. The performances (measured by both F₁ and MCC scores) of models in these permutated datasets were only slightly higher than the performance measurements calculated by comparing the responsiveness and its random permutation ( Figure 3D ).

At last, we assessed the separation of a series of well-known biomarkers for ICI responsiveness between groups predicted as positive and negative independent from those used to label the cases. The comparisons were made under the test sets of the 100 randomly trained models in the model evaluation step. First, the PD-1 and CTLA4 gene expression levels were significantly higher in cases predicted as positive than those predicted as negative (both with p < 10^-32, Mann Whitney U test, Figures 3E, F ). Second, there were also outstandingly more CD8⁺ T cells in cases in the positive cases than negative ones (p < 10^-52, t test, Figure 3G ). Last, there were higher proportions of cases with microsatellite instability (MSI) (39) in those positive cases than in those negative ones (p < 10^-52, t test, Figure 3H ).

In conclusion, the fine-tuned SVM model based on the selected probes came up with high performances in prediction of the ICI responsiveness in this pan-cancer cohort. The prediction accuracy is remarkably better than randomly selected probe sets of the same size, and we confirmed the improvement cannot be achieved by overfitting.

It was reported previously that the ICI responsiveness could be predicted by the gene expression levels of genes indicating the immune state of the cases in the pan-cancer level (32). Here we compared the performances of the model based on DNA methylation levels with that of models based on gene expression level (32). The gene expression level–based model was built exactly the same as reported before (32), except for taking the advantage of random oversampling to account for the class imbalance. The accuracy of the native built model was consistently higher than that originally reported [the mean and 95% confident interval MCC score under the 100 random test splits were 0.445 (0.371, 0.501) and 0.296(0.287, 0.306), for the native built one and originally reported, respectively] (32).

The performance of the methylation-based model was competitive with that based on gene expression levels. When assessing the performances using 100 times random splits during the model evaluation step, the methylation-based model got better performances when measured by F ₁ score (with mean F₁ score 0.4971 and 0.4844, p = 2.82×10^-4 Wilcoxon signed rank test, Figure 4A ). The performances were not notably different as measured by MCC score (with mean MCC score 0.4516 and 0.4446, p = 0.06, Figure 4B ), while lower when measured using AUC score (with mean AUC score 0.8914 and 0.8949 p = 0.01, Figure 4C ). The comparability between the performances of two models was also indicated by the closely located ROC curves between the two models ( Figure 4D ).

The genes in which the selected methylation sites located were then compared with those selected in the gene expression–based model. An enormous distinction between the two gene sets was observed. Only 189 out of 2,023 genes selected in the methylation model were found in the genes selected in the expression-based model (2,614 in total). The two gene sets also enriched few common GO and KEGG terms ( Figures 2D and 4E ). These observations indicated that the prediction power of DNA methylation and gene expression profiles could be complementary, and combining the two profiles would further enhance the prediction power.

To validate our assumption of complementarity of DNA methylation and gene expression profile, we further developed an SVM model based on the combination of the methylation levels of the selected methylation sites and the gene expression levels in the gene expression–based models. The best super-parameter λ was selected using 5-fold cross-validation. Performances of the selected model were again assessed in a 100 times random split of the cases into 80% training and 20% testing sets.

The combined model surpassed both the methylation-based and gene expression–based ones under all performance measure [all with p < 10^-10 for measurement F ₁, MCC, and AUC scores ( Figures 4A–C ), Wilcoxon signed rank test, with mean scores 0.5773, 0.5340, and 0.9327, respectively]. The high performance of the combined model was also validated by the ROC curves ( Figure 4D ).

These observations have implicated a scheme to enhance the predictive models of ICI responsiveness by assembling the multi-omics data.

To evaluate the performance of the methylation-based prediction model in each tumor type, we tested the model in tumor types with more than 5%, and 20 cases marked as positive (10 types in total: SKCM, BLCA, UCEC, CESC, COAD, LUAD, LIHC, STAD, LUSC, HNSC, ordered by the proportion of cases marked as position, Figure 5A ).

One important characteristic of responsiveness to ICI treatment was that the effective rate differed from tumor type to tumor type (40). As in the TCGA dataset, the proportion of cases marked as positive varied largely, from 37% in SKCM to zero in five other types ( Figure 5A and Supplementary Table S1 ). So, we first investigated that whether the methylation-based model could predict this variation. We calculated the proportion of cases in each tumor type predicted as positive and compared the numbers with the true proportions in the 100 times random test sets. Although due to the high false positive rate, which was the common characteristic of such kind of models (41, 42), the predicted proportion of cases marked as positive were always higher than the ground truths ( Figure 5A ), the two correlated tightly, with average Spearman’s correlation coefficient 0.74 [with 95% confidence interval (0.46,0.93) among the 100 random test sets, Figure 5A , embedded panel].

Next, we assessed the performances of the model in each tumor type using F ₁ and MCC scores. We compared the performances of the model based on the selected probes in the 100 randomly split train and test datasets with the performances of the models based on randomly selected probes of the same number, under the same hyper-parameter λ = 13, which was optimal for the selected probes-based model ( Figure 5B , upper panel). A paired t test was used to compare the performance between the two types of models. As expected, in 4 out of 10 tumor types, the performance of model based on selected probes was notably higher than that of model based on randomly selected probes (with significant level p < 0.1, Figure 5B , lower panel). These tumor types included SKCM and HNSC, which were commonly admitted as tightly related to the immune checkpoint evasion and may benefit from the ICI treatment ( Figure 5B , lower panel) (43, 44). On the other hand, only 2 and 3 tumor types were tested with significantly lower F ₁ and MCC scores than the random ones (p < 0.1) with the paired t test. This result further illustrated the high performance of the pan-cancer prediction model in tumor type level.

At last, we tested the performance of the methylation-based prediction model in an independent validation cohort. The cohort was taken from two newly published research on non-small-cell lung carcinoma (NSCLC) patients accepting anti-PD-1/PD-L1 treatments with clinical responses measured (13, 33). There were 60 and 18 cases included in the two datasets, with 14 and 6 being identified as responsive to the treatments, respectively. It was worth noticing that this cohort was not included in the TCGA cohort used for model building, and mutation burden was tested as a poor predictor for treatment responsiveness (13). This was the only publically available tumor immunotherapy cohort with DNA methylation levels measured to the best of our knowledge. The methylation levels in this cohort were measured using Illumina Infinium HumanMethylation850 BeadChip. We extracted the probes in the feature set existing in this chip and retrained the model in TCGA data. No significant drop of the performance of the model was observed (with average F ₁ and MCC score 0.4949 and 0.4503, with standard deviation 0.0288 and 0.0308). After applying the retrained model in this cohort, it got F ₁ = 0.4255, MCC = 0.1899, and AUC = 0.6742. The performance was also indicated in the ROC curve ( Figure 5C ). The progression-free survival time (PFS) was also significantly prolonged for cases predicted as positive compared with the negative ones ( Figure 5D ), though the p-value (p=0.06, log rank test) was not so significant due to the limited positive cases (only 20 cases predicted as positives). The performance of the model in this totally independent cohort demonstrated its efficiency at both pan-cancer level and for specific tumor type. It also indicated that the information the model caught was indeed the responsiveness itself other than its indicators such as TMB since the model retained its performance when TMB was not predictable to the responsiveness (13).