Chromothripsis status of five types of cancer including BRCA, lung adenocarcinoma (LUAD), ovarian cancer (OV), stomach adenocarcinoma (STAD) and UCEC was retrieved from compbio.med.harvard.edu/chromothripsis. Of which, 68 samples indicated with “High confidence” or “Linked to high confidence” in any chromosome was categorized chromothripsis and 80 samples with “No” in all chromosomes was categorized as non-chromothripsis. Clinical information and RNA-seq expression profiles of TCGA Pan-Cancer Atlas Studies were obtained from www.cbioportal.org. A total number of 8912 tumor samples with RNA-seq expression profiles in the format of linear RSEM were finally included in pan-cancer analysis. Raw data of GSE194040, GSE91061, GSE11636 were downloaded from Gene Expression Omnibus (GEO). Raw count of RNA-seq and responsiveness to atezolizumab monotherapy in advanced bladder cancer was extracted from R package “IMvigor210CoreBiologies”. IC₅₀ values in the natural logarithm of 268 anti-cancer drugs was downloaded from www.cancerrxgene.org/downloads/. The CEL files of Cancer Cell Line Encyclopedia(CCLE) cell lines was obtained from www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-3610/. For RSEM and raw count data, variation stabilizing transformation (VST) implemented in DESeq2 was applied prior to calculation of gene signatures and inference of abundance of infiltrated immune cells. For microarray profiled with affymetrix, function “rma” in R package limma was used to obtain expression matrix at probes level with default parameters. For raw data of Agilent, function “background Correct” and “normalize Between Arrays” were used to correct background and normalization. Function “neqc” was used for preprocessing data profiled with Illumina HumanHT-12 microarrays. The detailed sample size of each dataset used in this study were summarized in Supplementary Table 1 .

Given considerable difference in expression profiles among various cancers, the differentially expressed genes between chromothripsis (n=68) and non-chromothripsis (n=80) was identified regrading each cancer type as block by using edgeR. Gene set enrichment analysis (GSEA) was used to infer enriched pathways with gene rank based on log2 fold change and Molecular Signatures Database hallmark c5 gene sets. DEGs were filtrated with criteria: false discovery rate < 0.05 and absolute log2 fold change > 1. We took the advantage of gene set variation analysis (GSVA) by using DEGs listed in the Supplementary Table 3 . The upregulated and downregulated genes generated two gene sets, which were referred to as “CHPsigP” and “CHPsigN”, respectively. The CHP score for each sample was calculated by taking the difference between the enrichment score of CHPsigP and CHPsigN which were both derived from the GSVA. Gene set variation analysis(GSVA) was utilized to VST normalized RNA-seq or microarray expression matrix with these two gene sets. The parameter mx.diff was set to TRUE in all GSVA analysis. The CHP score for each sample was calculated by subtracting enrichment score of CHPsigN from enrichment score of CHPsigP.

Cytotoxic T lymphocyte (CTL score) (16), interferon-gama(IFN-γ)-related T-cell-inflated gene signature (17), activated dendritic cells (DC score) (18), pan-fibroblast TGF-β response signature (FTR score) (19) was evaluated by calculating geometrical mean of expression of marker genes belonging to corresponding signature. The gene signature of Cytotoxic T lymphocyte was the marker genes of cell subcluster “Cytotoxic lymphocytes” originally presented in MCPcounter (16). For tumor immune dysfunction and exclusion (TIDE score), only dysfunction of T lymphocyte cells was considered and measured by calculating Pearson’s correlation coefficients between gene expression of marker genes and mean z statistic values of these marker genes from Cox regression in at least two cancer types as proposed by the original article (20). Abundance of 19 immune cell populations was evaluated with ConsensusTME. Calculation of all gene signatures and ConsensusTME were performed in VST normalized RNA-seq or microarray expression matrix.

TCGA dataset of BRCA patients with available mutation data was used for this analysis. The dataset was composed of 51 patients, which occurred chromothripsis events (N = 32) or non-chromothripsis events (N = 19). Analysis, summarize and visualize of mutation data were performed by maftools (https://github.com/PoisonAlien/maftools). Mutually exclusive or co-occurring set of genes could be detected using SomaticInteractions function, which performs pair-wise Fisher’s Exact test to detect such significant pair of genes. We also analyzed the differences in APOBEC enriched samples proportion and identified differentially altered genes in two cohorts. We took non-negative matrix factorization to decompose the matrix into n signatures, and shown comparison of similarities of detected mutational signatures against validated COSMIC mutational signatures.

Kruskal-Wallis test with multiple comparisons or Wilcoxon test was used to evaluate differences in CHP score and other gene signatures among groups. In pan-cancer analysis, CHP score was categorized into high and low group with median value in each cancer type. Kaplan-Meier curves along with log rank test were used to evaluate difference of progression free survival(PFS) and overall survival(OS) between high and low CHP score groups. Univariate and multivariate Cox regression was used to evaluate hazard ratio and 95% confidence intervals for both continuous and categorized CHP score and other gene signatures for PFS and OS in pan-cancer analysis. In GSE194040, the difference in pathologically complete response (pCR) rate between high and low CHP score groups in each treatment arm was estimated by using Fisher’s exact probability. Univariate and multivariate logistic regression was used to evaluate the association of CHP score and other genes signatures and pCR. Receiver operating characteristic curve (ROC) analysis was conducted to evaluate the efficacy of CHP score to discriminate samples with chromothripsis from samples without chromothripsis or patients with responsiveness to immune checkpoint inhibitors (ICIs) from patients without response. All tests were two-sided. P value less than 0.05 was considered statistically significant.

Flow chart of this study was shown in Supplementary Figure 1 . To identify the possible gene expression signature of tumors with chromothripsis, as defined by the criteria of the Pan-Cancer Analysis of Whole Genomes (7), patient material that underwent WGS and was confirmed for chromothripsis status was matched to their corresponding in-situ sequencing data in the TCGA database. After excluding tumor types with less than five chromothripsis positive samples, the expression data of the five remaining tumor types were collected ( Supplementary Table 1 ): BRCA (n = 57), LUAD (n = 21), OV (n = 17), UCEC (n = 33) and STAD (n = 21). Patients were subsequently divided into chromothripsis (n=68) and non-chromothripsis (n=80) groups, and follow-up GSEA revealed that many pathways related to anti-tumor immunity, e.g., MHC class II complex, cytokine activity, and the IFN-γ response, were significantly underrepresented in the chromothripsis group ( Figure 1A ). Conversely, pathways related to tumor immunosuppression, including the transforming growth factor-β(TGF-β) and fibroblast growth factor receptor(FGFR) signal pathways, were significantly overrepresented in the chromothripsis group. In addition, we observed an overrepresentation of stemness maintenance in the chromothripsis group ( Figure 1B ) ( Supplementary Table 2 ). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis revealed that among the top 20 pathways, 7 pathways were related to tumor immunity and protein-protein interactions ( Supplementary Figure 2 ).

Further data analysis identified 170 significantly DEGs between the two groups ( Supplementary Table 3 ). Gene Ontology (GO) analysis showed that the DEGs, particularly the down-regulated genes, were significantly enriched in the chromothripsis group for biological processes related to chemokine receptor (CXCR) binding, cytokine activity, chemokine activity, chemokine receptor binding, and MHC class II receptor activity (all adjust p < 0.05) ( Figure 1C ). Consistent with this finding, the IMMPORT database classified the differentially expressed immunoregulatory genes into four groups: Chemokines (CK), Antimicrobials (AMB), Antigen processing and presentation (APP), and Cytokine receptors (CKR). Whereas most of the immunoregulatory genes were significantly decreased in the chromothripsis group, the expression of fibroblast growth factor 17 (FGF17) was significantly up-regulated ( Figure 1D ; Supplementary Table 3 ). By analyzing the activation of immune cells, we found that in the chromothripsis group: (i) T-cell receptor (TCR), B-cell receptor (BCR) richness was lower (P>0.05); (ii) fraction altered, HR defects, and intratumor heterogeneity were higher (P>0.05); and (iii) the aneuploid score, fragments number, TGF-β response and TCR evenness were significantly increased (P<0.05) ( Figure 1E ; Supplementary Table 4 ). These results indicate a likely association between chromothripsis and tumor immune evasion.

To further validate the correlation between chromothripsis and tumor immunity, we used GSVA to develop a CHP score based on a chromothripsis related gene signature, which encompasses all DEGs between the chromothripsis and non-chromothripsis groups. The CHP score for each sample was calculated by subtracting the enrichment score of the downregulated gene sets from the upregulated gene sets (see Materials and Methods). Using ROC analysis, our studies show the ability of the CHP score to distinguish chromothripsis with area under the curve (AUC) values for BRCA (0.892,95%CI=0.799-0.985), UCEC (0.892,95%CI=0.799-0.985), LUAD (0.833,95%CI=0.636-1), OV (0.889,95%CI=0.718-1), and STAD (0.893,95%CI=0.736-1) ( Figure 2A ).

We next evaluated the CHP score using all 33 tumor types (n=8912) in the TCGA. First, uveal melanoma (UVM) was found to have the highest CHP score, whereas hepatocellular carcinoma (LIHC) had the lowest CHP score ( Figure 2B ). Second, Thorsson et al. identified six immune subtypes by pan-cancer analysis (21). We observed that type C5 (immunologically quiet) had the highest CHP score, while C2 (IFN-γ dominant) and C6 (TGF-β dominant) have the two lowest CHP scores ( Figure 2C ). Additionally, we found that the CHP score was negatively correlated with the IFN-γ response, Immune score, DC score, and CTL score and is associated with a low abundance of DC cells, NK cells, and monocytes in all 33 tumors. In partial tumors, such as BRCA and THYM, the CHP score was positively correlated with the abundance of plasma cells, fibroblasts and endothelial cells ( Figure 2D ).

We classified tumors into high and low CHP score groups for each cancer type according to the median value of the CHP score. Notably, the expression of immune checkpoint genes, including programmed cell death protein 1 (PD-1), programmed cell death ligand 1 (PD-L1), and Cytotoxic T-lymphocyte antigen number 4 (CTLA-4), was significantly lower in tumors with a high CHP score in nearly all cancer types (i.e., in 25, 27 or 28 of 33, respectively) ( Figure 2E ), as was the CTL score and the IFN-γ signature ( Figure 2F ) ( Supplementary Table 5 ). Altogether, the results indicate that the CHP score reflects a tumor immunosuppression phenotype, which is characterized by a deficiency in immune cell infiltration and cytotoxicity.

The efficacy of immunotherapy can be affected by many factors including tumor immune cell infiltration and/or activation, the expression of immune checkpoints, and tumor mutational burden. Our results revealed a relationship between chromothripsis and tumor immune evasion. We therefore used four GEO datasets, i.e., GSE91061, GSE111636, IMvigor210, and the immunotherapy arm of GSE194040, to validate the association of the CHP score with the immunotherapy response. The CHP score was found to be significantly higher in the non-responder group in comparison to the responder group in the GSE91061 and GSE194040 datasets ( Figure 3A ; Supplementary Table 6 , 7 ). Moreover, the immunotherapy response rate was significantly lower in the high CHP score group in these two datasets ( Figure 3B ). Survival analysis showed that the CHP score can reveal the prognosis of patients in the GSE91061 dataset, where patients with higher scores had lower PFS (log-rank P < 0.001) and OS (log-rank P= 0.0026). However, this relationship between the CHP score and survival did not hold true in the IMvigor210 dataset (log-rank P > 0.05) ( Figure 3C ). ROC analysis, which was used to verify the above observations, supported the conclusion that the CHP score predicts immunotherapy effectiveness. Specifically, the results indicate that the CHP score achieves good performance in both the GSE91061 (AUC=0.831,95%CI=0.705-0.956) and the GSE111636 (AUC=0.767,95%CI=0.443-1) dataset, but not in the IMvigor210 dataset (AUC=0.548,95%CI=0.466-0.625) ( Figure 3D ). These findings support the link between chromothripsis and an impaired response to immune checkpoint blockades and indicate the potential value of the CHP score in predicting the efficiency of immunotherapy in LUAD and BRCA cancers.

Many studies have revealed a close relationship between chromothripsis and therapeutic resistance. To explore the predictive ability of the CHP score for other treatments in the setting of neoadjuvant therapy, we used the GSE194040 dataset, which contains information on 988 breast cancer patients who received 13 arms of neoadjuvants ( Supplementary Table 1 ). First, we analyzed the immune cells infiltration and immune activation signatures(see Materials and Methods) of all patients. Similar to the results of the TCGA database, the CHP score is negatively correlated with the abundance of most immune cells, including cytotoxic cells, dendritic cells, and CD4⁺, CD8⁺ T cells, but positively correlated for the abundance of fibroblasts ( Figure 4A ). Next, logistic regression analysis was performed to determine the correlation between various immune signatures and the pCR. Like some known immune signatures (e.g., IFN-γ, immune score and T cell abundance), the CHP score is significantly associated with the response to Paclitaxel+Pembrolizumab treatment (R=-3.34, P=0.002). Furthermore, the CHP score can predict the response to Paclitaxel (R=-1.62, P=0.03), Paclitaxel+AMG386 (R=-3.25, P<0.001), Paclitaxel+ABT888+Carboplatin (R=-3.71, P=0.001), Paclitaxel+Ganetespib (R=-1.62, P=0.038), and Paclitaxel+Ganitumab (R=-1.75, P=0.05) ( Figure 4B ) ( Supplementary Table 8 ). The CHP score was higher in the non-pCR group of patients that received one of the above six treatments ( Figure 4C ). After correcting for hormone receptor and human epidermal growth factor receptor 2 (HER-2) factors, multivariate logistic regression analysis validated that the CHP score was significantly correlated with a poor response to Paclitaxel+Pembrolizumab (R=-2.84, P=0.012), Paclitaxel+AMG386 (R=-2.56, P=0.014), and Paclitaxel+ABT888+Carboplatin (R=-3.16, P=0.011) ( Supplementary Table 9 , Supplementary Figure 2 ). Meanwhile, the pCR of patients with a high CHP score was significantly reduced in these three treatment groups (all p < 0.05) ( Figure 4D , Supplementary Table 10 ). ROC analysis confirmed the predictive power of the CHP score with regards to therapies involving Paclitaxel+Pembrolizumab (AUC=0.728,95%CI=0.609-0.848), Paclitaxel+AMG386 (AUC=0.731,95%CI=0.605-0.856), or Paclitaxel+ABT888+Carboplatin (AUC=0.724,95%CI=0.624-0.824) ( Figure 4E ). The results collectively indicate that the CHP score not only associates with immunotherapy responses but with tolerance to a range of anti-tumor treatment regimens, particularly anti-angiogenic inhibitors and platinum-based compounds in cases of breast cancer.

The Genomics of Drug Sensitivity in Cancer (GDSC) database was next used to screen for potentially sensitizing drugs for tumors with a high CHP score. After excluding down-regulated genes that are mostly related to tumor immune regulation, the enrichment score of CHPsigP was used to identify chromothripsis status instead of the CHP score for each cell line. The Pearson correlation coefficient was then used to determine sensitivity to different drugs, where a significant negative correlation between the CHPsigP score and the half maximal inhibitory concentration (IC₅₀) indicates a higher drug sensitivity. As shown in Figure 5A , the CHPsigP score is positively correlated with the IC₅₀ of insulin like growth factor 1 receptor (IGF1R), receptor tyrosine kinase(RTK), DNA replication, and PI3K-mTOR pathway inhibitors, indicative of a tolerance against these compounds in CHPsigP high tumors. Notably, there was a negative correlation between the CHPsigP score and the IC₅₀ of EGFR signal pathway inhibitors, suggesting that such drugs might be effective sensitizers against CHPsigP high tumors. Indeed, as revealed in ANOVA analysis, Pearson coefficients of drugs targeting EGFR signaling were significantly lower than those of drugs targeting other pathways ( Supplementary Table 11 ). Furthermore, the IC₅₀ of lapatinib, a HER-2 and EGFR tyrosine kinase inhibitor, is significantly negatively correlated with the CHPsigP score in breast (r=-0.457, P<0.001) and gastric cancer (r=-0.529, P=0.004) ( Figure 5B ). A similar relationship was also seen for osimertinib in breast (r=-0.490, P<0.001) and gastric cancer (r=-0.455, P=0.017) ( Figure 5C ).

Many studies have shown that chromothripsis is associated with poor clinical outcomes (7, 22). Univariate cox regression uncovered that the CHP score was significantly correlated with worse PFS in BRCA (HR=1.89,95%CI=1.07-3.36,P=0.029), UCEC (HR=2.79,95%CI=1.5-5.2,P=0.0012), ACC (HR=35.9,95%CI=7.54-170.86,P<0.001), BLCA (HR=1.84,95%CI=1.02-3.31,P=0.042), KIRC (HR=1.89,95%CI=1-3.57,P=0.049), PRAD(HR=3.1,95%CI=1.07-8.95,P=0.038), and SKCM (HR=1.72,95%CI=1.05-2.8,P=0.031) ( Figure 6A ) and with worse OS in BRCA (HR=1.92,95%CI=1.10-3.34,P=0.021), UCEC (HR=2.33,95%CI=1.12-4.86,P=0.024), ACC (HR=11.15,95%CI=1.78-69.87,P=0.010), KIRC (HR=1.85,95%CI=1.01-3.39,P=0.045), SKCM (HR=2.37,95%CI=1.34-4.20,P=0.0031), and SARC (HR=2.33,95%CI=1.05-5.16,P=0.037) ( Figure 6B ). There was a significant difference in PFS between CHP score high and low tumors in cases of ACC, UCEC, BLCA, KICH, and SKCM, and in OS in ACC, KICH, and SKCM ( Figures 6C, D ). Thus, the CHP score is a prognostic factor for several types of cancers, and a high score is associated with a poor outcome in cases of ACC, UCEC, BLCA, KICH, and SKCM.

Genome mutation profiles were used to analyze the potential mechanisms underlying chromothripsis. Due to the limited numbers of chromothripsis positive samples in most cancer type collections, we were restricted to analyzing genomic mutation features in breast cancer only. As shown by Figure 7A , compared with the non-chromothripsis group (n=19), tumors with chromothripsis (n=32) have a higher frequency of mutations in the VCAN and KIAA1109 genes, which are related to immunotherapeutic efficiency and prognosis in several types of cancer (23, 24). Somatic Interactions analysis showed that the chromothripsis group has more co-mutation patterns than the non-chromothripsis group, possibly related to the hallmark chromosome rearrangements ( Figure 7B ). For instance, BRCA2-MUC16, PI3KCA-MUC17, ATRX-ATP1A4 et al. co-mutations were identified in chromothripsis group.

Signature of single base substitution (SBS) analysis was next performed, and the signatures of SBS1 (Spontaneous or enzymatic deamination of 5-methylcytosine), SBS2 (C > T mutation mediate by APOBEC cytidine deaminase), and SBS3 (Defects in DNA-Double strand break(DSB) repair by HR) were found in the chromothripsis group, whereas non-chromothripsis was enriched in SBS5 (Aetiology: Unknown) and SBS13 (C > G mutation mediate by APOBEC cytidine deaminase) ( Figure 7C ). Moreover, the proportion of APOBEC enriched samples is significantly higher in the chromothripsis group than in the non-chromothripsis group (37.5% vs. 15.8%) ( Figure 7D ), and the TCW motif burden of APOBEC enriched samples is also higher in the former group (0.35 vs. 0.24) ( Figure 7E ). These results reveal a specific mutation signature in breast cancers with chromothripsis and indicate that HR defect and APOBEC may play important roles in the occurrence and process of chromothripsis.