The raw dataset of the CCLE project (DepMap release 23Q4) was downloaded via the DepMap portal (https://depmap.org/portal/). The dataset consisted of six matrices including gene expression, CNV (copy number variation) amplification, CNV loss, DNA mutation of somatic point mutations and indels, pharmacologic sensitivity, and metadata of cell lines. For correlation analysis and survival analysis, the expression profiles of DLBCL patients were downloaded from TCGA (https://www.cancer.gov/tcga) and NCBI Gene Expression Omnibus (GEO) GSE69049.

We prepared the CCLE datasets only for 504 cell lines presented in the pharmacologic sensitivity matrix and built the input matrices of JNMF as follows. The TPMs (transcripts per million) on a log-2 scale quantifying gene expressions were converted ranging from 0 to 1 by a min-max normalization. The DNA mutation profile was converted into a binary matrix where 1 for mutated and 0 for normal. Two binary matrices for CNV gain and loss were constructed by the GISTIC scores in the CNV data; +2 for amplification and −2 for deletion. The pharmacologic sensitivity was also converted into values ranging from 0 (insensitive) to 1 (sensitive) as follows: max X   –   x   /   max X   –   min X , where x is an IC50 (half maximal inhibitory concentration) value and max(X) and min(X) are the maximum and minimum values in the pharmacologic sensitivity profile. A cancer-type matrix in binary format was prepared from the metadata of cell lines.

To generate single-guide RNA (sgRNA) expression vectors targeting COPS5 (COPS5-sgRNA-1 and COPS5-sgRNA-2) or a non-targeting (NT) control, annealed oligonucleotides were cloned into the pLKO5.sgRNA.EFS.tRFP657 vector (Addgene plasmid # 57824; http://n2t.net/addgene:57824), which was a gift from Benjamin Ebert. The Cas9 expression in the cell lines Raji (Burkitt lymphoma cell line) and SLVL (Splenic marginal zone lymphoma cell line) was induced by FUCas9Cherry plasmid, which was a gift from Marco Herold (Addgene plasmid # 70182, http://n2t.net/addgene:70182). Lentiviruses were produced by transient transfection of 293T cells with viral plasmids, along with gag-, pol-, and env-expressing plasmids (pMD2.G and psPAX2) using the calcium-phosphate method (Goyama et al., 2016). pMD2.G (Addgene plasmid #12259; http://n2t.net/addgene:12259) and psPAX2 (Addgene plasmid #12260; http://n2t.net/addgene:12260) were gifts from Didier Trono. The sequences for the sgRNAs are as follows: NT: 5′-cgc​ttc​cgc​ggc​ccg​ttc​aa-3′, COPS5-sgRNA-1: 5′-gtg​atg​cat​gcc​aga​tcg​gg-3′, COPS5-sgRNA-2: 5′-caa​caa​gaa​caa​tat​ccg​ca-3′.

The lymphoma cell lines Raji and SLVL were cultured in RPMI1640 medium supplemented with 10% fetal bovine serum (FBS) and 1% penicillin. 293T cells (CRL-11268, ATCC, Manassas, VA, United States) were cultured in Dulbecco’s modified Eagle’s medium supplemented with 10% FBS and 1% penicillin. These cells were first transduced with the FUCas9Cherry, followed by sorting of mCherry⁺ cells using BD FACSAriaIII (BD Biosciences, San Jose, CA, United States). The Cas9-expressing (mCherry⁺) Raji and SLVL cells were then transduced with the sgRNAs co-expressing tRFP657. The frequency of tRFP657⁺ cells in the cultures was evaluated on Day3, Day6, Day10, Day17 and Day24.

Cells were sorted by FACS Aria (BD Biosciences, San Jose, CA, United States), and the expression of mCherry and tRFP657 was analyzed with FACS CytoFLEX (Beckman Coulter, Brea, California, United States). The cytometry data were analyzed by BD FlowJo software (TREESTAR, Inc., San Carlos, CA. ver.10.8.1).

JNMF, an extension of traditional NMF algorithm, is designed to facilitate the simultaneous decomposition of N datasets (Zhang et al., 2012). The objective function of JNMF is given by: min ∑ i = 1 N X i − W H i F 2 , where X _i is an input matrix with size m×n _i , and F represents the Frobenius norm. The W and H _i represent m×k and k×n _i factorized matrices, respectively. Here, k is the number of clusters to be extracted, namely, common pattern modules (CPMs). To find the optimal W and Hi minimizing the objective function, JNMF updates them based on the traditional multiplication update formulas (Zhang et al., 2012) as follows: W i a = W i a ∑ J = 1 N X J H J T i a W ∑ K = 1 N H K H K T i a H I a μ = H I a μ W T X I a μ W T W H I a μ , I = 1 , ⋯ , N .

To handle missing values in the input matrices, we employed a weighted NMF approach (Zhang et al., 2012): a mask matrix M representing 1 for non-missing and 0 for missing cases is introduced into the JNMF framework. X is accessed by the Hadamard product with M effectively filtering the influence of missing values (Fujita et al., 2018). The objective function and multiplicative update rules of the weighted JNMF are given by: min ∑ i = 1 N M i ∘ X i − W H i F 2 W i a = W i a ∑ J = 1 N M i ∘ X J H J T i a ∑ K = 1 N M i ∘ W H K H K T i a H I a μ = H I a μ W T M i ∘ X I a μ W T M i ∘ W H I a μ , I = 1 , ⋯ , N , where A ∘ B = a i j b i j is Hadamard product.

Given a factorization rank k, JNMF starts with random W and Hs and updates the random matrices toward minimizing the objective function for n iterations step-by-step. To tune the hyperparameters, whether the procedure with k and n stably convergent in repeating t times is monitored. A consensus matrix and its cophenetic correlation coefficient (CCC) (Brunet et al., 2004) evaluate the performance of JNMF under the setting of k and n.

In analyzing the six input matrices of CCLE datasets, our JNMF produces seven matrices: W, H ₁ , H ₂ , H ₃ , H ₄ , H ₅ , and H ₆ . For normalizing the matrices, rather than using the maximum values, z-score normalization was applied to each row and column of W and Hs as follows: z i j = x i j − μ i σ i , where μ i and σ i stand the average and the standard deviation for cell line/drug/mutation/CNV/genes/cancer type feature j, respectively. The feature j is assigned to CPMs if and only if z i j is > +1.96.

To identify the activated pathways within each CPM obtained from JNMF, we employed IPA (Ingenuity Pathway Analysis) (Krämer et al., 2014). IPA is a widely used tool for exploring signaling pathways and biological networks. Specifically, IPA’s upstream analysis and protein-protein interaction analyses were performed to determine general regulators in CPM-activating pathways.

For examining the correlation between candidate biomarkers and known hub genes of DLBCL, we utilized GEPIA2 (Gene Expression Profiling Interactive Analysis 2) (Tang et al., 2019). GEPIA2 is a comprehensive web-based tool that integrates data from the TCGA and GTEx databases. Specifically, we used the “correlation analysis” module within GEPIA2, utilizing the TCGA-DLBCL project dataset, to assess the correlation between the candidate biomarkers and the known hub genes of DLBCL.

The Kaplan-Meier method was used to investigate whether the candidate biomarkers in GSE69049 datasets affected the overall survival (OS) of DLBCL patients treated with chemotherapy. The “survminer” R package was used to explore survival analysis.

We designed a method to find significant sets of multimodal factors (i.e., modules) that have the potency to characterize disease phenotypes. Since multimodality comprises multiple high-dimensional heterogeneous data, we adopted JNMF, a well-proven algorithm for clustering the factors as modules, by modifying it to handle data sparseness efficiently, referred to as weighted JNMF. To interpret the functional importance of the JNMF-detected modules, namely, the common pattern module (CPM), we utilized a pathway analysis by detecting upstream regulators in signaling pathways activated by the modules (Figure 1).

To perform the biomarker discovery using the proposed method, we complied six feature matrices for 504 cancer cell lines by processing the large-scale CCLE datasets, each of which has binary or continuous values: gene expression, CNV amplification, CNV loss, DNA mutation, pharmacologic sensitivity, and cancer type (Supplementary Table S1). Using these matrices as inputs, the weighted JNMF masks sparse elements and generates a factor matrix W given by latent coefficients for the cancer types in the reduced k-dimensional space (i.e., rank). Simultaneously, the six input matrices are reduced to the k dimension generating the matrices H _i (i = 1, … ,6), where each product with W approximates the original input matrix. Given optimal W and Hs, significant feature sets (>+1.96 in z-score) are captured as CPMs by investigating z-score distributions in the factorized matrices.

Subsequently, the gene expression profiles in the JNMF-detected CPMs are analyzed by IPA to identify activated pathways. Meanwhile, the disease subtypes corresponding to the CPM are identified by combining characteristic drugs, DNA mutations, and structural variants in the CPM. Moreover, through IPA upstream analysis and IPA causal network analysis, the upstream regulators of CPMs which are considered candidate biomarkers are investigated. Next, the correlation and survival analyses of candidate biomarkers using TCGA and GEO data examine the clinical significance of the candidate biomarkers. Ultimately, the candidate biomarkers are experimentally validated by the CRISPR-Cas9 system.

To interrogate the ability of our JNMF, we prepared three artificial datasets as used in a previous study (Fujita et al., 2018). We first constructed three matrices with random values imprecating noise: a binary matrix for mutation and continuous matrices for pharmacologic sensitivity and gene expression. Then, we inserted missing values into randomly selected 10% of the entries of each matrix: the missing rates in CCLE datasets were 2.6%–10.5%. Next, we embedded three or four predefined CPMs into the matrices and randomly shuffled the entries in each matrix. Thereby, three artificial input data that are noisy and sparse but include modules were generated (Supplementary Table S2, Figure 2).

Since the embedded CPMs that our JNMF has to detect were three or four, we set the rank k = 4. As shown in Figure 2, our JNMF successfully identified the CPMs by decomposing the input matrices into W, H ₁ , H ₂ , and H ₃ . In addition, the products WH ₁ , WH ₂ , and WH ₃ accurately restored the input matrices. This result demonstrates that our JNMF can uncover hidden relationships within high-dimensional multimodal datasets by reducing the influence of noise and missing values.

To optimize the hyperparameters of JNMF for the six input matrices prepared from CCLE gene expression, CNV amplification, CNV loss, DNA mutation, pharmacologic sensitivity, and cancer cell line type, we investigated the convergence of JNMF during 2000 iterations updating W and Hs (Figure 3A). We noticed that the outputs of JNMF are stable (CCC = 0.72) when k = 40 and manifest substantial consistency across 10 repeats (Figure 3B). Using these hyperparameters, we finally retrieved 40 CPMs corresponding to specific cancer types, such as CPM #1 for hematopoietic and lymphoid malignancies, CPM #7 for breast cancer, CPM #10 for malignant melanoma, and CPM #28 for endometrial cancer (Table1; Figure 3C).

The CPM #1 contained TP53 mutation and the malignant translocations of MYC, BCL2, and BCL6, which are relevant to DLBCL (Chapuy et al., 2018). The CPM #7 included the pharmaceutical agent Lapatinib, HER2 overexpression, and BRCA2 amplification, which is supported by the clinical application of Lapatinib for treating HER2-positive metastatic breast cancer (Xu et al., 2021). The CPM #10 covering almost the entire skin cancer cell lines included the pharmaceutical agent PLX4720 (vemurafenib) and BRAF mutation; the efficacy of vemurafenib has been tested in several clinical trials for treating unresectable or metastatic melanoma with BRAF V600E mutation (Chapman et al., 2011). The CPM #28 confined several known diagnostic and prognostic biomarkers of endometrial cancer, such as the mismatch repair mutation, and the overexpression of MLH1, MSH2, MSH6, and PMS2. In addition, this module included the overexpression and mutation of PTEN and TP53, which significantly contribute to the diagnosis of endometrial cancer (Crosbie et al., 2022).

Collectively, we confirmed that the CPMs are in high concordance with known relationships among variants, medications, and cancers, which suggests the potency of our approach to the discovery of novel biomarkers from multimodal data.

To interpret the functionality of JNMF-detected CPMs, we focused on the CPM #1 that includes DLBCL biomarkers (Table 1). Notably, this module also contained several gene mutations, each of which is known to be involved in pathogenic pathways: MYD88 and CARD11 functioning in the NF-κB pathway, SPEN involved in the NOTCH signaling, and STAT3 which is a pivotal member of the JAK/STAT signaling pathway (Chapman et al., 2011).

Next, we performed a pathway analysis of IPA with the gene expression profile of CPM #1. Consistent with the pathways in which the mutated genes are involved, we found the activation of several signaling pathways of DLBCL (p < 0.05). For example, the activation of the NF-κB pathway causing DLBCL (Odqvist et al., 2014), the deregulation of the JAK-STAT pathway and PI3K-mediated signaling pathway which is the essential contributor to the pathogenesis and poor prognosis of DLBCL (Chapuy et al., 2018) (Supplementary Figure S1). In addition, the sub-networks centered on MYC, TP53, and NF-κB indicated the activation of downstream pathways potentially relevant to DLBCL development (Supplementary Figure S2). Lastly, by performing the IPA upstream analysis and IPA causal network analysis, we investigated upstream regulators (p < 0.05) likely controlling the gene expression of CPM #1 (Tables 2).

Interestingly, COPS5, a subunit of the COP9 signalosome complex, was identified from both analyses and showed regulatory interactions with the overexpressed genes in the CPM #1 including known DLBCL markers (Figure 4). This result supports that COPS5 overexpression affects tumor-negative regulators in diverse cancers (Wang et al., 2016). For example, COPS5 activates MYC by mediating the SKP1-CUL1-F-box protein complex (Lyapina et al., 2001). Also, COPS5 alters cytoplasmic localizations of TP53 and induces the degradation of TP53 (Li et al., 2013). Moreover, COPS5 is co-expressed with STAT3 in cancers (Nishimoto et al., 2013), and MYC mediates this co-regulation associated with poor prognosis (Ok et al., 2014).

Taken together, despite the considerable heterogeneity among DLBCL subtypes (Schmitz et al., 2018), the CPM #1 collectively retained the distinctive characteristics of DLBCL, which emphasizes the importance of understanding the orchestration of multiple oncogenic factors. Moreover, we identified COPS5 using the information of this module as a potential upstream regulator of DLBCL, requiring further validations.

Using the gene expression profiles of 47 DLBCL patients available at TCGA, we sought to confirm the importance of COPS5 in DLBCL patients. Consistent with the results of CCLE data analysis, we observed the positive expression correlation between COPS5 and the marker genes, as well as between TP53 and MYC (Figure 5A). It is noteworthy to mention that all the positive correlations have been reported by previous in vitro studies (Adler et al., 2006; Sitte et al., 2012; Li et al., 2013; Nishimoto et al., 2013; Luo et al., 2022). Interestingly, even BCL6, MYC, and TP53 were grouped in the CPM #1, their expression correlations in the patients were relatively less as shown in Figure 5A, E.g., R = 0.22 between BCL6 and TP53, R = 0.24 between BCL6 and MYC. This result might reflect the high heterogeneity of DLBCL. Therefore, COPS5 may be a key hub gene that efficiently characterizes heterogeneous DLBCL by co-expressing with the marker genes.

Next, we performed the survival analysis to inspect the relationship of COPS5 with the prognosis of DLBCL patients. As shown in Figure 5B, the Kaplan-Meier curve exhibited that the high expression of COPS5 is associated with poor prognosis in the overall survival of the patients treated with chemotherapy alone (p = 0.0168). This result supports that the over expression of COPS5 promotes malignancy.

Finally, Since COPS5 has been shown to be associated with the proliferation of DLBCL in several cell lines derived from subtypes of DLBCL (Pulvino et al., 2015), to further determine the role of COPS5 in B cell lymphoma, we depleted COPS5 in the B cell lymphoma cell lines Raji and SLVL using the CRISPR/Cas9 system. Raji and SLVL cells were transduced with Cas9 together with tRFP657-coexpressing non-targeting (NT) or COPS5-targeting sgRNAs. As shown in Figure 5C, COPS5 depletion showed a strong growth-inhibitory effect in Raji and SLVL cells. These results confirm the key role of COPS5 in the proliferation of malignant B cells (Supplementary Figure S3).