We retrieved the raw transcriptome profiling and clinical data of DLBCL patients and cell lines from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) and the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo). The data without survival time were removed from the data. Expression profiling by microarray of DLBCL and normal B-cell lines was obtained from the GSE12195 (73 DLBCL samples versus 20 normal tonsil samples) and GSE56315 (55 DLBCL samples versus 33 normal tonsil samples) datasets with Affymetrix Human Genome U133 Plus 2.0 Platform. Forty-eight DLBCL samples from TCGA DLBC and 305 from GSE10846 (with Affymetrix Human Genome U133 Plus 2.0 Platform) with complete transcriptional and clinical data were preserved for later analysis. Then, we merged the transcriptional expression of TCGA DLBC and GSE10846 into one dataset and identified it as the entire cohort after the removal of the batch effect using the “limma” package (28). Batch effects were corrected by the ComBat algorithm in the sva package (29). The entire cohort was divided into the training and test cohorts equally by the “caret” package (30). Furthermore, GSE87371 including 223 patients with Affymetrix Human Genome U133 Plus 2.0 Platform was attained as an external validation cohort. We also downloaded simple nucleotide variant (SNV) data of DLBCL from TCGA DLBC, and copy number variation (CNV) data from UCSC Xena (https://xenabrowser.net/datapages/).

We searched “necroptosis” in PubMed (https://pubmed.ncbi.nlm.nih.gov/) and filtered out 198 NRGs from several articles (5, 6, 31, 32). All NRGs are listed in Supplementary Table 1 . With the help of the “limma” and “reshape2” packages (33), we compared the expression of DLBCL and normal cell lines and identified 118 NRGs in GSE12195 and GSE56315. Among the 118 NRGs, 103 NRGs were identified as the shared genes in the entire cohort. Then, we performed Cox analysis and found 25 NRGs that were correlated with survival (p < 0.05) using the “limma”, “survival” (34), and “survminer” packages (35). To explore their connections, the gene network and protein–protein interaction (PPI) were generated by the “igraph” package (36) and STRING (https://cn.string-db.org/). Based on the 25 NRGs, two clusters (named NRG clusters) were defined by the “ConsensusClusterPlus” package. We also drew the Kaplan–Meier (K-M) survival analysis curve using the “survival” and “survminer” packages. Principal component analysis (PCA) plotted by the “limma” and “ggplot2” packages (35) could assess the precision of NRGs. Heatmap including clusters and various clinical data was drawn using the “pheatmap” package (37). Furthermore, we ran the “GSVA” (38) and “GSEABase” packages (39) to clarify the relationship of function and immune cells with 25 NRGs.

The univariate analysis of differentially expressed genes (DEGs) associated with prognosis in NRG clusters was performed using the “limma” package. For these DEGs, we conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis using the “clusterProfiler” package (40). Later, to make it more accurate, we obtained 24 DEGs after the Least Absolute Shrinkage and Selection Operator (LASSO) regression and Cox proportional hazard model construction by the “survival”, “survminer”, and “timeROC” packages (41). The risk score of each sample was calculated based on the expression of the DEGs. The risk score = Coef1*Exp1 + Coef2*Exp2 +…+ Coef24*Exp24. The entire cohort was subdivided into low- and high-risk groups by the median value of risk score. We then drew the K-M curve, risk curve, and survival status curve of the entire, training, test, and external validation cohort by “pheatmap”. In order to estimate its predictive value, we drew the receiver operating characteristic (ROC) curve at 1, 3, and 5 years using the “survival”, “survminer”, and “timeROC” packages. Meanwhile, cluster analysis was performed again based on the DEGs using the “ConsensusClusterPlus” package (42) and the K-M curve of two clusters (named DEG clusters) was made using the “survival” and “survminer” packages. The Sankey diagram between NRG clusters, risk score, DEG clusters, and clinical information was drawn using “limma”, “ggpubr” (https://CRAN.R-project.org/package=ggpubr), and “ggalluvial” (https://CRAN.R-project.org/package=ggalluvial).

To improve the accuracy of the current prognosis score system, we performed univariate and multivariate Cox regression analysis of risk score and clinicopathological characteristics using the “survival” package. The factors with p-value < 0.05 were contained to draw the nomogram, calibration curve, and ROC curve using the “survival” package. As there were two regimens in the dataset, we divided them into the CHOP and RCHOP groups and drew the ROC curves. To evaluate the stability of the predictors in the nomogram, we perform collinearity diagnosis by SPSS.

We performed the “maftools” package (43) to show the TMB and used “Rcircos” (44) to identify the CNV and its location on the chromosome for all NRGs in TCGA DLBC. To determine the immune cells related to DEGs and risk score, CIBERSORT algorithm was employed. The TMB was mapped and evaluated by the “maftools” package between the low- and high-risk groups, so as the TMB and risk score. To figure out whether the risk score affected drug selection, we used the “pRRophetic” package (45) to pick out the targeted drugs between the low- and high-risk groups.

All analyses were conducted by using the R software (version 3.4.0) and SPSS software (version 24.0). The references of all packages were all cited. K-M and Cox regression were used to identify the prognosis-related NRGs. We compared statistical differences between the two groups by the Wilcoxon test. We performed survival analysis by the K-M method and analyzed the difference of OS by log-rank test. LASSO and the Cox proportional hazard model were used to identify the most significant genes. We performed univariate and multivariate analysis and built the nomogram by the Cox regression method. The sensitivity and specificity of the nomogram were implemented by the calibration curve and the area under the curve (AUC) of the ROC curve. Collinearity analysis was achieved by the collinearity diagnosis. Statistical significance was defined as p < 0.05 (unless otherwise specified).

The clinical characteristics of the GEO and TCGA datasets are summarized in Table 1 . The expression, genetic mutation landscape, and the CNV location of NRGs in DLBCL are shown in Figure 1 . There were 118 NRGs that were significantly differentially expressed between normal cells and DLBCL (p < 0.05) (see Supplementary Figure 1 ). Among the 118 NRGs, 34 of them were mutated in 59.46% (22/37) of DLBCL samples ( Figure 1A ). STAT3 was the most frequent mutated NRGs (14%) followed by TNFAIP3, FAS, and STAT6 (11%, 8%, and 8%, respectively). It was not hard to find that the missense mutation was the most common mutation variant, and C>T was the highest SNV class. As for the CNV, nearly all of them were revealed as gain or loss after CNV analysis ( Figure 1B ). CDKN2A, TNFAIP3, BACH2, and MAP3K7 were the most deleted NRGs while FASLG, GLUL, PLA2G4A, and USP21 gained the most. Then, we mapped the exact location of the CNV alteration on chromosomes ( Figure 1C ), which mainly enriched in chromosomes 1, 2, and 9. We identified 25 NRGs according to the gene expression, survival time, and follow-up status by Cox analysis. The K-M curve based on 25 genes is shown in Supplementary Figure 2 . Moreover, the gene network is shown in Figure 1D . The size of each circle indicated its significance in DLBCL. The right part of the circle represented its role in DLBCL. Purple meant it was a risk factor while green meant it was favorable. The line between genes indicated that they were correlated with each other (p < 0.0001). Meanwhile, we also explored the relationship between the proteins encoded by the NRGs ( Figure 1E ). CASP8 could interact with other proteins greatly.

According to the 25 NRGs, we established the NRG clusters from the TCGA DLBC and GEO datasets using the “ConsensusClusterPlus” package ( Figure 2A ). To evaluate the prognosis value of the cluster, we drew the K-M curve, PCA, and heatmap ( Figures 2B–D ). Figure 2B shows that cluster A had a worse outcome than B (p < 0.001). PCA could easily distinguish cluster A from cluster B. The link between gene expression, NRG clusters, and clinical characteristics is shown in Figure 2D . Moreover, we exploited the KEGG pathway and immune cells based on the 25 NRGs by gene set variation analysis (GSVA) and single-sample Gene Set Enrichment Analysis (ssGSEA) as well ( Figures 2E,F ). In cluster A, the downregulated genes were enriched in the TGF-β signal pathway and metabolism of glucose, amino acid, and lipid while the overexpressed genes were enriched in the Notch signal pathway, metabolism, and transductions ( Figure 2E ). As for the immune cells in cluster A, activated dendritic cells (DCs), CD56dim natural killer cells (NK cells), myeloid-derived suppressor cells (MDSCs), monocytes, and plasmacytoid DCs were upregulated while activated CD4⁺T cells, γ/δ T cells, immature B cells, mast cells, NK cells, neutrophils, regulatory T cells, type 17 helper T cells, and type 2 helper T cells were downregulated ( Figure 2F ).

To understand the potential significance between the two clusters based on NRGs, we then used the “limma” and “clusterProfiler” packages to clarify the underlying function of DEGs. GO pathway analysis showed that the DEGs were enriched in neutrophil activation, neutrophil-mediated immunity, neutrophil degranulation, T-cell activation, mononuclear cell differentiation, and lymphocyte differentiation in biological process (BP) ( Figure 3A ). The DEGs were mainly involved in the PI3K-Akt signal pathway, cytokine–cytokine receptor interaction, endocytosis, the chemokine signal pathway, the JAK-STAT3 signal pathway, and the NOD-like receptor signal pathway in KEGG analysis ( Figure 3B ). There were 24 DEGs between the NRG cluster ( Figure 3C and Supplementary Table 2 ). We divided the entire cohort into two DEG clusters ( Figure 3D and Supplementary Figure 3 ). K-M showed a significant difference between the DEG clusters (C, D) (p < 0.001) ( Figure 3E ), indicating the potential of DEGs in DLBCL.

The entire cohort was randomly assigned into the training cohort and test cohort to examine the predictive potential of DEGs. After calculating the risk score, we divided the entire cohort and GSE87371 into the low- and high-risk groups. The K-M curve suggested a significant difference in the entire cohort, training cohort, and test cohort ( Figures 4A,D,G ). The risk curve and survival status curve showed that it could separate the high-risk group from the other ( Figures 4B,E,H ). To assess the prognostic value further, we plotted the ROC curve. The AUCs were 0.814, 0.812, and 0.815 at 1, 3, and 5 years, respectively, in the entire cohort ( Figure 4C ). For the training cohort, the AUCs were 0.957, 0.964, and 0.977, respectively ( Figure 4F ). AUCs in the test cohort were 0.685, 0.641, and 0.668 at 1, 3, and 5 years, respectively ( Figure 4I ). In the external validation cohort, the risk score could distinguish between the low- and high-risk groups (p < 0.001) ( Figures 4J,K ) with AUCs of 0.834, 0.834, and 0.777 at 1, 3, and 5 years, respectively ( Figure 4L ). These results elucidated that the DEG-based risk score could perfectly predict the clinical outcome of DLBCL patients.

The univariate and multivariate Cox regression analysis indicated that age, cell of origin (COO), LDH, and risk score served as the independent prognostic factors ( Figure 5A ), and we constructed the nomogram based on them ( Figure 5B ). The nomogram predicted that OS at 1, 3, and 5 years was approximately close to the observed OS in the calibration curve ( Figure 5C ). To assess the factors in the nomogram, we performed collinearity diagnosis and the results showed that the collinearity tolerance of each predictor was great than 0.1 and the variance inflation factor was less than 10. Condition index was less than fifteen, and the maximum value in the horizontal row was not simultaneously greater than 0.5 in the vertical column ( Supplementary Tables 3 , 4 ). These results indicated that there was no collinearity in the independent predictors. Compared with the COO plus International Prognostic Index (IPI, scoring with available data) or IPI alone, the nomogram could remarkably increase the AUC (the AUC is 0.723, 0.712, 0.552, and 0.537, respectively), suggesting its vital role in the prognosis of DLBCL ( Figure 5D ). Different regimens for DLBCL have an effect on the prognosis as well. Our results showed that the risk score could markedly enhance the efficacy in both CHOP and RCHOP groups ( Figures 5E,F ).

By comparing the risk score between NRG clusters, we found that the risk score of cluster A was significantly higher than that of cluster B (p = 2.3e-06) ( Figure 6A ). As seen in Figure 6B , several NRGs showed a significant difference between the low- and high-risk group. Apart from that, we plotted the Sankey diagram of NRG clusters, DEG clusters, risk score, follow-up status, extranodal site, and Ann Arbor stage ( Figure 6C ). The NRG clusters, DEG clusters, risk score, and follow-up status corresponded well while the clinical characteristics varied widely.

The relationship between NRGs and immune cells had been demonstrated before, and we herein investigated the correlation between DEGs and immune cells. In Figure 7A , T regulatory cells (Treg), γ/δ T cells, and CD4 activated memory T cells ranked most with the DEGs, correlated with 18/24 (75%), 18/24 (75%), and 17/24 (70.83%) genes, respectively ( Figure 7A ). Between the low- and high-risk group, γ/δ T cells and Macrophage 1 (M1) were negatively correlated with risk score while NK resting cells and Macrophage 2 were positively correlated ( Figure 7B ).

In addition, the CNV of NRGs in the whole samples had been displayed above, but they varied between the low- and high-risk group. As was shown in Figure 7C , IGLV3-1, KMT2D, GHV2-70, CARD11, HIST1H1E, MYD88, and MCU16 mutated the most (22%) with 6/9 (66.67%) altered samples in the high-risk group. Intriguingly, IGLV3-1 had the highest mutation (62%), followed by CARD11 (50%), GHV2-70 (38%), IGLC2 (38%), B2M (38%), SOCS1 (38%), and HIST1H1C (38%) with 7/8 (87.5%) altered samples in the low-risk group ( Figure 7D ). Missense mutation was the most common variant classification in both groups.

TMB was a pivotal part in tumors as well. NRG-related TMB appeared to be different apparently between the low- and high-risk groups ( Figure 7E ). In DEG clusters, NRG-related TMB was positively related to risk score ( Figure 7F ).

To testify the therapeutic potential of DEGs, we performed the correlation of different risk groups with current anti-tumor drugs. Twenty drugs were discovered between the low- and high-risk group ( Figure 8 and Supplementary Figure 4 ). We found that NVP-BEZ235 and BX.795 were sensitive in the high-risk group ( Figures 8A,B ) while the rest were more sensitive in the low-risk group.