OSCR

Coordinated Multicellular Immune Programs and Drug Targets Revealed by Single-Cell Analysis in Driver-Mutated NSCLC.

Code ↔ Paper

40 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 40 matches
  1. [1] § 2. Results › 2.5. Identifying the Driver Mutation Specific Prognostic Signature (DMSP.sig) with Immune Characteristics ↔ code/FigureS5/FigureS5.R, lines 233–261 · score 0.96 · PAX6 ITGB1, USF2 IGF1, PPARG RARRES2, SNAI1 ITGB4, NR1H4, CEBPB IGF1
  2. [2] § 2. Results › 2.5. Identifying the Driver Mutation Specific Prognostic Signature (DMSP.sig) with Immune Characteristics ↔ code/Figure5/Figure5.R, lines 994–1061 · score 0.96 · PAX6 ITGB1, USF2 IGF1, PPARG RARRES2, SNAI1 ITGB4, NR1H4, CEBPB IGF1
  3. [3] § 2. Results › 2.4. Potential Mediators Underlying the Poor-Prognosis CM2 and CM5 ↔ code/Figure4/Figure4.R, lines 319–359 · score 0.81 · Astrocyte_05_HLA, Ligand receptor, ANGPTL, CCL, GALECTIN, MIF
  4. [4] § 4. Materials and Methods › 4.19. Immune Infiltration Analysis ↔ code/Figure5/Figure5.R, lines 813–892 · score 0.81 · get_immu_ratio, corr.test, Immune infiltration, Pearson, GSVA, matched
  5. [5] § 2. Results › 2.5. Identifying the Driver Mutation Specific Prognostic Signature (DMSP.sig) with Immune Characteristics ↔ code/FigureS5/FigureS5.R, lines 233–261 · score 0.80 · SNAI1 ITGB4, RELA IGF1, SP1 ITGA6, risk scores, S5D, ROC
  6. [6] § 4. Materials and Methods › 4.12. Cellular Interaction Analysis ↔ code/FigureS4/FigureS4.R, lines 179–244 · score 0.79 · netAnalysis_signalingRole_heatmap, CellChatDB, netAnalysis_computeCentrality, Network centrality, probabilities, database
  7. [7] § 4. Materials and Methods › 4.12. Cellular Interaction Analysis ↔ code/Figure4/Figure4.R, lines 234–274 · score 0.78 · CellChat, netAnalysis_signalingRole_heatmap, netAnalysis_computeCentrality, Network centrality, weights, interactions
  8. [8] § 2. Results › 2.4. Potential Mediators Underlying the Poor-Prognosis CM2 and CM5 ↔ code/Figure1/Figure1.R, lines 1105–1150 · score 0.75 · Macrophage_06_C3, mast cells, Astrocyte_05_HLA, stromal cells, CHIT1, MS4A7
  9. [9] § 4. Materials and Methods › 4.16. Prognostic Model Construction ↔ code/Figure5/Figure5.R, lines 412–486 · score 0.72 · Cox regression, prognostic model, glmnet, LASSO, split, GSVA
  10. [10] § 4. Materials and Methods › 4.9. Cancer Cell State Identification ↔ code/Figure3/Figure3.R, lines 1–86 · score 0.70 · CopyKat, Euclidean, chromosome, distance, cancer cell, v2
  11. [11] § 2. Results › 2.5. Identifying the Driver Mutation Specific Prognostic Signature (DMSP.sig) with Immune Characteristics ↔ code/Figure5/Figure5.R, lines 994–1061 · score 0.70 · SNAI1 ITGB4, RELA IGF1, SP1 ITGA6, longer, RS, regulons
  12. [12] § 2. Results › 2.1. A Single-Cell Atlas in Non-Small Cell Lung Cancer (NSCLC) Harboring Driver Mutations ↔ code/Figure1/Figure1.R, lines 478–528 · score 0.68 · cell proportions, NK cells, EGFR co mutation, KRAS co mutation, EGFR BM, driver mutation
  13. [13] § 4. Materials and Methods › 4.13. Metabolic Analysis ↔ code/Figure4/Figure4.R, lines 687–703 · score 0.68 · scMetabolism, KEGG metabolic, AUCell, scores
  14. [14] § 4. Materials and Methods › 4.13. Metabolic Analysis ↔ code/FigureS4/FigureS4.R, lines 245–269 · score 0.68 · scMetabolism, KEGG metabolic, AUCell, scores
  15. [15] § 2. Results › 2.1. A Single-Cell Atlas in Non-Small Cell Lung Cancer (NSCLC) Harboring Driver Mutations ↔ code/Figure1/Figure1.R, lines 1105–1150 · score 0.66 · HLA DQB1, stromal cells, NK cells, plasma, DRA, myeloid
  16. [16] § 4. Materials and Methods › 4.1. Sample Characteristics ↔ code/FigureS1/FigureS1.R, lines 30–107 · score 0.64 · EGFR co mutation, KRAS co mutation, MET BM, S1, EGFR BM, tissue
  17. [17] § 4. Materials and Methods › 4.6. Differential Gene Expression and Pathway Enrichment Analysis ↔ code/FigureS2/FigureS2.R, lines 20–57 · score 0.64 · qvalueCutoff, pvalueCutoff, BP, GO, enrichment, enriched
  18. [18] § 2. Results › 2.4. Potential Mediators Underlying the Poor-Prognosis CM2 and CM5 ↔ code/Figure4/Figure4.R, lines 319–359 · score 0.63 · Macrophage_06_C3, Astrocyte_05_HLA, pro, CHIT1, MS4A7, CD68
  19. [19] § 4. Materials and Methods › 4.15. TF Identification ↔ code/Figure5/Figure5.R, lines 176–233 · score 0.63 · Chi square, chisq.test, TFs
  20. [20] § 2. Results › 2.4. Potential Mediators Underlying the Poor-Prognosis CM2 and CM5 ↔ code/Figure4/Figure4.R, lines 181–232 · score 0.62 · communication networks, cell communication, signaling pathways, incoming, outgoing, identity
  21. [21] § 4. Materials and Methods › 4.14. HR Analysis ↔ code/Figure2/Figure2.R, lines 75–162 · score 0.62 · Worse survival, Better survival, age, Cox, GSVA, HR
  22. [22] § 4. Materials and Methods › 4.14. HR Analysis ↔ code/Figure4/Figure4.R, lines 422–503 · score 0.62 · Worse survival, Better survival, age, TCGA, Cox, GSVA
  23. [23] § 4. Materials and Methods › 4.19. Immune Infiltration Analysis ↔ code/Figure2/Figure2.R, lines 1–73 · score 0.62 · corr.test, TCGAplot, Pearson, GSVA, LUSC, LUAD
  24. [24] § 2. Results › 2.2. Cellular Module Analyses Reveal Five Tumor Immune Microenvironment (TIME) Subtypes ↔ code/Figure2/Figure2.R, lines 364–432 · score 0.62 · MHC class, co inhibition, APC, proliferation, antigen, pathways
  25. [25] § 4. Materials and Methods › 4.6. Differential Gene Expression and Pathway Enrichment Analysis ↔ code/Figure3/Figure3.R, lines 403–453 · score 0.61 · qvalueCutoff, pvalueCutoff, BP, GO, enriched, genes
  26. [26] § 2. Results › 2.1. A Single-Cell Atlas in Non-Small Cell Lung Cancer (NSCLC) Harboring Driver Mutations ↔ code/Figure2/Figure2.R, lines 463–546 · score 0.60 · HLA DRA, HLA DQB1, EGFR BM, plasma, NK, astrocytes
  27. [27] § 4. Materials and Methods › 4.6. Differential Gene Expression and Pathway Enrichment Analysis ↔ code/Figure2/Figure2.R, lines 243–293 · score 0.60 · compareCluster, enrichKEGG, fun, enrichment, Gene
  28. [28] § 4. Materials and Methods › 4.17. Spatial Transcriptomics-Based Regulon Co-Localization Analysis ↔ code/FigureS7/FigureS7.R, lines 58–88 · score 0.59 · SpaGene_LR, co localization, tissue
  29. [29] § 4. Materials and Methods › 4.1. Sample Characteristics ↔ code/Figure1/Figure1.R, lines 313–355 · score 0.59 · EGFR co mutation, KRAS co mutation, MET BM, EGFR BM, driver mutation, HER2
  30. [30] § 4. Materials and Methods › 4.16. Prognostic Model Construction ↔ code/Figure5/Figure5.R, lines 412–486 · score 0.58 · Cox regression, prognostic model, GSVA, coefficient, regulons, sig
  31. [31] § 2. Results › 2.5. Identifying the Driver Mutation Specific Prognostic Signature (DMSP.sig) with Immune Characteristics ↔ code/FigureS5/FigureS5.R, lines 188–231 · score 0.58 · regulatory network, gene associated, cancer cells, CNV, HR, risk
  32. [32] § 4. Materials and Methods › 4.6. Differential Gene Expression and Pathway Enrichment Analysis ↔ code/Figure3/Figure3.R, lines 475–511 · score 0.58 · compareCluster, enrichKEGG, fun, Gene
  33. [33] § 2. Results › 2.2. Cellular Module Analyses Reveal Five Tumor Immune Microenvironment (TIME) Subtypes ↔ code/Figure2/Figure2.R, lines 243–293 · score 0.58 · CC chemokines, functional enrichment, antigen, enriched, receptor, Figure 2
  34. [34] § 2. Results › 2.4. Potential Mediators Underlying the Poor-Prognosis CM2 and CM5 ↔ code/FigureS4/FigureS4.R, lines 1–41 · score 0.56 · Astrocyte_05_HLA, S4A, intensity, COL10A1, HIGD1B, MKI67
  35. [35] § 2. Results › 2.3. Cellular Module (CM)2 and CM5 Are Associated with Poor Prognosis ↔ code/Figure3/Figure3.R, lines 567–607 · score 0.55 · KRAS co mutations, GE3, GE6, GE7, GE13, GE2
  36. [36] § 4. Materials and Methods › 4.20. Drug Sensitivity Prediction and Molecular Docking of Core Targets ↔ code/Figure7/Figure7.R, lines 15–62 · score 0.54 · calcPhenotype, GDSC2, Drug, sensitivity, predicted, trained
  37. [37] § 2. Results › 2.6. Identification of Potential Therapeutic Agents Targeting Prognostic Biomarkers ↔ code/Figure7/Figure7.R, lines 15–62 · score 0.53 · FSC231, JW, LY, drug, sensitivity, RELA
  38. [38] § 4. Materials and Methods › 4.5. Cell Module Identification ↔ code/FigureS3/FigureS3.R, lines 1–67 · score 0.51 · ward.D2, v2, corr, v1, matrix, modules
  39. [39] § 2. Results › 2.1. A Single-Cell Atlas in Non-Small Cell Lung Cancer (NSCLC) Harboring Driver Mutations ↔ code/Figure1/Figure1.R, lines 478–528 · score 0.51 · NK cells, EGFR BM, driver mutations, innate, immunity, stromal
  40. [40] § 2. Results › 2.2. Cellular Module Analyses Reveal Five Tumor Immune Microenvironment (TIME) Subtypes ↔ code/Figure2/Figure2.R, lines 748–799 · score 0.51 · PI3K, pathway activity, MAPK, CM1, Figure 2, NSCLC

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 895 lines · 44 KB · MIT · 7 matches

  1. # ------------- Figure 2 --------------
  2. #----Figure 2A----
  3. #co-occurrence analyses
  4. library(Seurat)
  5. library(dplyr)
  6. library(ggplot2)
  7. library(tidyverse)
  8. setwd("./NSCLC/Figure/figure2/")
  9. rm(list = ls())
  10. TIME <- readRDS(file = "./NSCLC/Figure/figure1/TIME.rds")
  11. sce <- TIME
  12. #remove cell populations with low cell counts, low frequencies, or unidentified cell populations
  13. length(table(sce$subcluster))
  14. freq.low <- c("Unidentified")
  15. freq.high <- sce$subcluster[!sce$subcluster %in% freq.low]
  16. sce <- sce[,!sce$subcluster %in% freq.low]
  17. #Calculate the frequency of cell subsets; the row labels are the samples, and the column labels are the cell subsets
  18. write.csv(prop.table(table(sce$orig.ident,sce$subcluster),1),file = "./cell frequency cor.csv")
  19. data.freq <- data.frame(read.csv(file="./cell frequency cor.csv"),row.names = 1)
  20. data.freq <- data.freq[,!colnames(data.freq) %in% freq.low]#Remove cell clusters with low cell counts or frequencies
  21. rowSums(data.freq)#Determine the proportion of the cell population in the sample
  22. #Calculating Correlation
  23. library(corrplot)
  24. library(psych)
  25. library(pheatmap)
  26. #①corr.test(pearson)
  27. data_test <- corr.test(data.freq, method= "pearson")#spearman/kendall
  28. data_cor <- data_test[["r"]]
  29. write.csv(data_cor,file = "./correlation of TIME cells.csv")
  30. #pheatmap
  31. library(circlize)
  32. col_fun <- colorRamp2(
  33. c( -0.5,0, 1),
  34. c("#638fa9","white","#b10928")
  35. #c("#01665e","white", "#8c510a")
  36. )
  37. pheatmap(data_cor,cluster_rows = T,cluster_cols = T,fontsize = 7,main = "correlation of TIME cells",border_color = NA,color = col_fun)
  40. #Perform HR analysis on the various cell subsets in TIME
  41. library(GSVA)
  42. library(TCGAplot)
  43. library(tidyverse)
  44. library(survival)
  45. library(meta)
  46. library(doParallel)
  47. #Set the number of parallel processing cores
  48. registerDoParallel(cores = 15)
  49. subcluster.marker <- read.csv(file="./yang/wang/TIME.markers1_TOP10.csv",row.names = 1)
  50. tpm <- get_all_tpm()
  51. meta <- get_all_meta()
  52. cancers <- c("LUAD","LUSC")
  53. genelist <- split(subcluster.marker$gene, subcluster.marker$cluster)
  54. genelist <- genelist[-c(16,60)]
  55. # Use a `foreach` loop to parallelize the outer loop
  56. final_pan <- foreach(geneset_name = names(genelist),
  57. .combine = bind_rows,
  58. .packages = c("GSVA", "tidyverse", "survival", "meta", "TCGAplot"),
  59. .export = c("tpm", "meta", "cancers")) %dopar% {
  60. current_geneset <- genelist[[geneset_name]]
  61. current_genelist <- list(current_geneset)
  62. names(current_genelist) <- geneset_name
  64. cox_results <- list()
  65. for (cancer in cancers) {
  66. exprSet <- subset(tpm, Group == "Tumor" & Cancer == cancer) %>%
  67. tibble::add_column(ID = stringr::str_sub(rownames(.), 1, 12), .before = "Cancer") %>%
  68. dplyr::filter(!duplicated(ID)) %>%
  69. tibble::remove_rownames() %>%
  70. tibble::column_to_rownames("ID") %>%
  71. dplyr::filter(rownames(.) %in% rownames(subset(meta, Cancer == cancer)))
  73. exprSet <- exprSet[, -(1:2)] %>% as.matrix() %>% t()
  75. # GSVA
  76. gsvapar <- gsvaParam(exprData = exprSet,
  77. geneSets = current_genelist,
  78. kcdf = "Gaussian")
  79. exprSet_gsva <- gsva(gsvapar)
  81. # Survival Analysis
  82. cl <- meta[colnames(exprSet_gsva), ]
  83. cl$symbol <- exprSet_gsva[geneset_name, ]
  85. if (sum(!is.na(cl$time) & !is.na(cl$event)) > 0) {
  86. m <- tryCatch(
  87. coxph(Surv(time, event) ~ symbol + age, data = cl),
  88. error = function(e) NULL
  89. )
  91. if (!is.null(m)) {
  92. beta <- coef(m)
  93. se <- sqrt(diag(vcov(m)))
  94. tmp <- round(cbind(
  95. HR = exp(beta),
  96. se = se,
  97. lower = exp(beta - 1.96 * se),
  98. upper = exp(beta + 1.96 * se),
  99. p = 1 - pchisq((beta/se)^2, 1)
  100. ), 3)
  101. cox_results[[cancer]] <- tmp["symbol", ]
  102. }
  103. }
  104. }
  106. if (length(cox_results) > 0) {
  107. a <- do.call(rbind, cox_results) %>%
  108. as.data.frame() %>%
  109. rownames_to_column("cancer") %>%
  110. select(cancer, HR, se, lower, upper, p)
  112. # Pan-cancer meta analysis
  113. meta_res <- metagen(log(a$HR), a$se, sm = "HR")
  114. pan_row <- data.frame(
  115. cancer = "Pan cancer",
  116. HR = exp(meta_res$TE.random),
  117. lower = exp(meta_res$lower.random),
  118. upper = exp(meta_res$upper.random),
  119. p = meta_res$pval.random,
  120. stringsAsFactors = FALSE
  121. )
  123. a <- a %>% select(cancer, HR, lower, upper, p)
  125. pan <- bind_rows(a, pan_row) %>%
  126. mutate(
  127. celltype = geneset_name,
  128. significance = case_when(
  129. p <= 0.0001 ~ "****",
  130. p <= 0.001 ~ "***",
  131. p <= 0.01 ~ "**",
  132. p <= 0.05 ~ "*",
  133. TRUE ~ ""
  134. ),
  135. color = case_when(
  136. HR > 1 ~ "Worse survival",
  137. HR < 1 ~ "Better survival",
  138. TRUE ~ "Neutral"
  139. )
  140. )
  141. pan
  142. } else {
  143. NULL
  144. }
  145. }
  147. stopImplicitCluster()
  148. subcluster.hr <- subset(final_pan,cancer %in%c("LUAD","LUSC"))
  149. subcluster.hr <- subcluster.hr %>% mutate(cellmodule = case_when(
  150. celltype %in% c("Bn_06_GZMB","pDC_03_GZMB","Plasma_09_IGHG_RRM2","Plasma_05_IGHG_IGHGP","Plasma_01_IGHA_MZB1","Plasma_02_IGHA_IGLC3","CD8T_02_GZMB",
  151. "CD8T_03_MKI67","CD8T_01_CCL5","Treg_FOXP3","B_02_HLA.DQB1","CD4T_02_GNB2L1","CD4T_06_ACSL5","CD8T_04_HLA-B") ~ "module1",
  152. celltype %in% c("Astrocyte_03_ROM1","Endothelial_02_PTMAP2","Fibroblast_01_Myo_HIGD1B","Fibroblast_02_Myo_COL10A1","Fibroblast_03_Myo_MKI67","B_04_PTMAP2","CD4T_04_SNORD3A","Astrocyte_05_HLA-DRB5","Neurophil_01_GOS2","Astrocyte_04_SCGB3A1","Macrophage_06_C3","Astrocyte_01_CAV1","Astrocyte_02_MALAT1","Astrocyte_06_CCL5") ~ "module2",
  153. celltype %in% c("NK_01_FGFBP2","CD4T_03_SCGB3A1","Endothelial_05_IL1RL1","Mast_01_GNB2L1","Mast_03_PBK","cDC_01_CLEC10A","cDC_02_CD207","Macrophage_02_FCN1","Macrophage_04_SPP1",
  154. "Macrophage_09_MKI67","CD4T_05_SFTPB","Plasma_06_IGHG_SFTPB","Endothelial_09_S100A14","Macrophage_08_SCGB3A2","Fibroblast_04_Alveolar_DPT","Macrophage_07_IGKC") ~ "module3",
  155. celltype %in% c("CD4T_01_IL7R","B_03_CD3E","B_05_FYN","Fibroblast_05_Alveolar_MYH11","Endothelial_07_NOTCH3","B_01_HLA-DRA","Plasma_07_IGHG_GRIFIN","Plasma_08_IGHG_GUSBP11","Plasma_03_IGHA_GUSBP11","Plasma_04_IGHA_VPS13D") ~ "module4",
  156. celltype %in% c("Endothelial_03_VCAM1","Endothelial_08_CCL21","Endothelial_01_FCN3","Endothelial_04_DKK2","Fibroblast_06_Alveolar_MS4A7","Endothelial_06_CD68","Mast_02_CHIT1","Macrophage_03_MARCO","Macrophage_01_LGMN","Macrophage_05_CSTB") ~ "module5",
  157. ))
  159. subcluster.hr$cellmodule <- factor(subcluster.hr$cellmodule,levels = c("module1","module2","module3",'module4','module5'))
  160. subcluster.hr <- subcluster.hr[order(subcluster.hr$cellmodule),]
  162. write.csv(subcluster.hr,file = "/data/yang/wang//subcluster.hr.csv")
  164. #----Figure 2B----
  165. #Analysis of Gene Differential Expression Across Cell Modules
  166. #Remove cell clusters that are less relevant and lie outside the module
  167. rm(list = ls())
  168. TIME <- readRDS(file = "./NSCLC/Figure/figure1/TIME.rds")
  169. subsce.corr <- TIME
  170. #corlow <- c("")
  171. #subsce.corr <- subsce.corr[,!subsce.corr$subcluster %in% freq.low]
  172. [email hidden] <- [email hidden] %>% mutate(cellmodule = case_when(
  173. subcluster %in% c("Bn_06_GZMB","pDC_03_GZMB","Plasma_09_IGHG_RRM2","Plasma_05_IGHG_IGHGP","Plasma_01_IGHA_MZB1","Plasma_02_IGHA_IGLC3","CD8T_02_GZMB",
  174. "CD8T_03_MKI67","CD8T_01_CCL5","Treg_FOXP3","B_02_HLA.DQB1","CD4T_02_GNB2L1","CD4T_06_ACSL5","CD8T_04_HLA-B") ~ "module1",
  175. subcluster %in% c("Astrocyte_03_ROM1","Endothelial_02_PTMAP2","Fibroblast_01_Myo_HIGD1B","Fibroblast_02_Myo_COL10A1","Fibroblast_03_Myo_MKI67","B_04_PTMAP2","CD4T_04_SNORD3A","Astrocyte_05_HLA-DRB5","Neurophil_01_GOS2","Astrocyte_04_SCGB3A1","Macrophage_06_C3","Astrocyte_01_CAV1","Astrocyte_02_MALAT1","Astrocyte_06_CCL5") ~ "module2",
  176. subcluster %in% c("NK_01_FGFBP2","CD4T_03_SCGB3A1","Endothelial_05_IL1RL1","Mast_01_GNB2L1","Mast_03_PBK","cDC_01_CLEC10A","cDC_02_CD207","Macrophage_02_FCN1","Macrophage_04_SPP1",
  177. "Macrophage_09_MKI67","CD4T_05_SFTPB","Plasma_06_IGHG_SFTPB","Endothelial_09_S100A14","Macrophage_08_SCGB3A2","Fibroblast_04_Alveolar_DPT","Macrophage_07_IGKC") ~ "module3",
  178. subcluster %in% c("CD4T_01_IL7R","B_03_CD3E","B_05_FYN","Fibroblast_05_Alveolar_MYH11","Endothelial_07_NOTCH3","B_01_HLA-DRA","Plasma_07_IGHG_GRIFIN","Plasma_08_IGHG_GUSBP11","Plasma_03_IGHA_GUSBP11","Plasma_04_IGHA_VPS13D") ~ "module4",
  179. subcluster %in% c("Endothelial_03_VCAM1","Endothelial_08_CCL21","Endothelial_01_FCN3","Endothelial_04_DKK2","Fibroblast_06_Alveolar_MS4A7","Endothelial_06_CD68","Mast_02_CHIT1","Macrophage_03_MARCO","Macrophage_01_LGMN","Macrophage_05_CSTB") ~ "module5",
  180. ))
  181. table(subsce.corr$cellmodule)
  182. getwd()
  183. saveRDS(subsce.corr,file = "./NSCLC/Figure/figure2/subsce.corr.rds")
  184. subsce.corr <- readRDS("./subsce.corr.rds")
  185. #Identify the differentially upregulated genes in each cell module
  186. DefaultAssay(subsce.corr) <- "RNA"
  187. Idents(subsce.corr) <- "cellmodule"
  188. subsce.corr.markers.up <- FindAllMarkers(subsce.corr,
  189. min.pct = 0.25,
  190. only.pos = T,
  191. logfc.threshold = 0.25)
  192. subsce.corr.markers1.up <- subsce.corr.markers.up [subsce.corr.markers.up $p_val_adj < 0.05, ]
  193. subsce.corr.markers1.up_TOP5 <- subsce.corr.markers1.up %>% group_by(cluster) %>% top_n(n = 5, wt = avg_log2FC)
  194. write.csv(subsce.corr.markers1.up,file = "./subsce.corr.markers1.up.csv")
  195. subsce.corr.markers1.up <- read.csv("./subsce.corr.markers1.up.csv",row.names = 1,header = T)
  196. subsce.corr.markers1.up$cluster <- factor(subsce.corr.markers1.up$cluster,levels = c("module1","module2","module3",'module4','module5'))
  197. subsce.corr.markers1.top30 <- subsce.corr.markers1.up %>% group_by(cluster) %>% top_n(n = 30, wt = avg_log2FC)
  198. subsce.corr.markers1.top30 <- subsce.corr.markers1.top30[order(subsce.corr.markers1.top30$cluster), ]
  199. write.csv(subsce.corr.markers1.top30,file = "./NSCLC/Figure/figure2/subsce.corr.markers1.top30.csv")
  200. #install.packages('devtools')
  201. #devtools::install_github('junjunlab/scRNAtoolVis')
  202. library(scRNAtoolVis)
  203. #DEG
  204. DefaultAssay(subsce.corr) <- "RNA"
  205. subsce.corr.markers <- FindAllMarkers(subsce.corr,
  206. min.pct = 0.25,
  207. only.pos = F,
  208. logfc.threshold = 0.25)
  209. subsce.corr.markers1 <- subsce.corr.markers [subsce.corr.markers $p_val_adj < 0.05, ]
  210. subsce.corr.markers1_TOP5 <- subsce.corr.markers1 %>% group_by(cluster) %>% top_n(n = 5, wt = avg_log2FC)
  211. write.csv(subsce.corr.markers1,file = "./subsce.corr.markers1.csv")
  212. #subsce.corr.markers1 <- read.csv("./subsce.corr.markers1.csv",row.names = 1)
  213. subsce.corr.markers1$cluster <- factor(subsce.corr.markers1$cluster,levels = c("module1","module2","module3",'module4','module5'))
  214. #save(subsce.corr,file = "F:/subsce.corr.RData")
  215. mycolors <- c("#f9b3ad","#80C1D7","#70cdbe","#f1dba6","#CEB6E0")
  216. mycolors <- c("#B1C1E3","#92D3EB","#C9E5D5","#E2EDBC","#FAE7F1")
  217. mycolors <- c("#9B7EBD","#70cdbe","#46aec8","#6fc08d","#457baa")
  218. mycolors <- c("#d72e2d","#ea9740","#316339","#7dbcd1","#313663")
  219. mycolors <- c("#a5c469","#d586b3","#8ea1c9","#dd9f85","#88bdaa")
  220. mycolors <- c("#70cdbe","#80C1D7","#968ab7","#d586b3","#eb7e60")
  221. jjVolcano(diffData = subsce.corr.markers1,
  222. log2FC.cutoff = 0.25,
  223. size = 2.5,
  224. fontface = 'italic',
  225. aesCol = c('#8fb4dc','#eb7e60'),
  226. tile.col = mycolors,
  227. #col.type = "adjustP",
  228. topGeneN = 5,
  229. fontface = 'italic',
  230. polar = T)+ ylim(-7,12)
  232. color = colorRampPalette(c("white","#e6573f", "#d23229"))(100)
  233. jjVolcano(diffData = subsce.corr.markers1,
  234. tile.col = mycolors,
  235. size = 3.5,
  236. pSize = 0.9,
  237. topGeneN = 5,
  238. base_size = 20,
  239. fontface = 'italic',
  240. aesCol = c('#8fb4dc','#d23229'),
  241. polar = T) + ylim(-12,10)#+xlim(-10,10)
  243. #----Figure 2C----
  244. #Functional enrichment analysis
  245. library(clusterProfiler)
  246. library(org.Hs.eg.db)
  247. library(dplyr)
  248. library(enrichplot)
  249. group <- data.frame(
  250. gene = subsce.corr.markers.up$gene,
  251. module = subsce.corr.markers.up$cluster
  252. )
  253. #SYMBOL → ENTREZID
  254. gene_id <- bitr(group$gene,
  255. fromType = "SYMBOL",
  256. toType = "ENTREZID",
  257. OrgDb = org.Hs.eg.db)
  259. mydf <- merge(gene_id, group, by.x = "SYMBOL", by.y = "gene")
  260. #KEGG
  261. kk <- compareCluster(
  262. ENTREZID ~ module,
  263. data = mydf,
  264. fun = "enrichKEGG",
  265. organism = "hsa",
  266. pAdjustMethod = "BH",
  267. pvalueCutoff = 0.05
  268. )
  269. dotplot(kk, showCategory = 10) +
  270. scale_shape_manual(values = 18)
  272. #----Figure 2D----
  273. #Gene set scoring
  274. subsce.corr <- readRDS("./subsce.corr.rds")
  275. DefaultAssay(subsce.corr) <- "RNA"
  276. #Antigen-presenting cell co-inhibitory genes
  277. APC_co_inhibition <- list(c("C10orf54","CD274","LGALS9","PDCD1LG2","PVRL3"))
  278. #Co-stimulatory genes of antigen-presenting cells
  279. APC_co_stimulation <- list(c("CD40","CD58","CD70","ICOSLG","SLAMF1","TNFSF14","TNFSF15","TNFSF18","TNFSF4","TNFSF8","TNFSF9"))
  280. #CC chemokine receptor genes
  281. CCR <- list(c("CCL16","TPO","TGFBR2","CXCL2","CCL14","TGFBR3","IL11RA","CCL11","IL4I1","IL33","CXCL12","CXCL10","BMPER","BMP8A","CXCL11","IL21R","IL17B","TNFRSF9","ILF2","CX3CR1","CCR8","TNFSF12","CSF3","TNFSF4","BMP3","CX3CL1","BMP5","CXCR2","TNFRSF10D","BMP2","CXCL14",
  282. "CCL28","CXCL3","BMP6","CCL21","CXCL9","CCL23","IL6","TNFRSF18","IL17RD","IL17D","IL27","CCL7","IL1R1","CXCR4","CXCR2P1","TGFB1I1","IFNGR1","IL9R","IL1RAPL1","IL11","CSF1","IL20RA","IL25","TNFRSF4","IL18","ILF3","CCL20","TNFRSF12A","IL6ST","CXCL13","IL12B","TNFRSF8",
  283. "IL6R","BMPR2","IFNE","IL1RAPL2","IL3RA","BMP4","CCL24","TNFSF13B","CCR4","IL2RA","IL32","TNFRSF10C","IL22RA1","BMPR1A","CXCR5","CXCR3","IFNA8","IL17REL","IFNB1","IFNAR1","TNFRSF1B","CCL17","IFNL1","IL16","IL1RL1","ILK","CCL25","ILDR2","CXCR1","IL36RN","IL34","TGFB1","IFNG",
  284. "IL19","ILKAP","BMP2K","CCR10","ILDR1","EPO","CCR7","IL17C","IL23A","CCR5","IL7","EPOR","CCL13","IL2RG","IL31RA","TNFAIP6","IFNL2","BMP1","IL12RB1","TNFAIP8","IL4R","TNFRSF6B","TNFAIP8L1","TNFRSF10B","IFNL3","CCL5","CXCL6","CXCL1","CCR3","TNFSF11","CSF1R","IL21","IL1RAP","IL12RB2",
  285. "CCL1","IL17RA","CCR1","IL1RN","TNFRSF11B","TNFRSF14","IL13","IL2RB","BMP8B","CCL2","IL24","IL18RAP","TGFBI","TNFSF10","TNFRSF11A","CXCL5","IL5RA","TNFSF9","IL1RL2","TNFRSF13C","IL36G","IL15RA","TNFRSF21","CXCL8","IL22RA2","TNFAIP8L2","IL18R1","IFNLR1","CXCR6","CCL3L3","TNFRSF1A","IL17RE","IFNGR2",
  286. "IL17RC","TNFAIP8L3","ILVBL","TGFBRAP1","CCL4L1","CSF2RA","CCRN4L","CCL26","TNFAIP1","CCRL2","IFNA10","TNFRSF17","IFNA13","IL20","IL18BP",'CCL3L1',"TNFSF12-TNFSF13","IL5","IL23R","IL26","TNF","TGFA","CSF2","IL1F10","CXCL17","TNFSF13","IFNA4","IL37","IL12A","IL7R","IFNA1","IL1A","IL4","IL2",
  287. "CCL22","CSF3R","IL10","IFNK","TGFB2","IL1R2","IL1B","IL17F","IL27RA","IL15","TNFSF8","IL36B","XCL1","CXCL16","TNFRSF19","IL3","CCL3","IFNA2","BMPR1B","IFNA21","TNFSF18","CCL8","IL17RB","TNFRSF25","IL22","IL10RB","IFNAR2","CCL18","IFNA16","CSF2RB","IL36A","TNFAIP3","IL13RA2","IL13RA1","CCR9","TNFRSF10A",
  288. "IFNA7","IFNW1","XCL2","TNFSF14","CCR2","BMP15","BMP10","CCL15-CCL14","TGFBR1","IFNA5","BMP7","IFNA14","IL20RB","IL10RA","IFNA17","CCR6","TGFB3","CCL15","CCL4","CCL27","TNFRSF13B","TNFAIP2","IL31","IL17A","TNFSF15","CCL19","IFNA6","IL9"))
  289. #Checkpoint genes
  290. Check_point <- list(c("IDO1","LAG3","CTLA4","TNFRSF9","ICOS","CD80","PDCD1LG2","TIGIT","CD70","TNFSF9","ICOSLG","KIR3DL1","CD86","PDCD1","LAIR1","TNFRSF8","TNFSF15","TNFRSF14","IDO2","CD276","CD40","TNFRSF4","TNFSF14","HHLA2","CD244",
  291. "CD274","HAVCR2","CD27","BTLA","LGALS9","TMIGD2","CD28","CD48","TNFRSF25","CD40LG","ADORA2A","VTCN1","CD160","CD44","TNFSF18","TNFRSF18","BTNL2","C10orf54","CD200R1","TNFSF4","CD200","NRP1"))
  292. #Genes encoding cell lysis activity
  293. Cytolytic_activity <- list(c("PRF1","GZMA"))
  295. #Human leukocyte antigen genes
  296. HLA <- list(c("HLA-E","HLA-DPB2","HLA-C","HLA-J","HLA-DQB1","HLA-DQB2","HLA-DQA2","HLA-DQA1","HLA-A","HLA-DMA","HLA-DOB","HLA-DRB1","HLA-H","HLA-B","HLA-DRB5","HLA-DOA","HLA-DPB1","HLA-DRA","HLA-DRB6","HLA-L","HLA-F","HLA-G","HLA-DMB","HLA-DPA1"))
  298. #pro-inflammatory genes
  299. Inflammation_promoting <- list(c("CCL5","CD19","CD8B","CXCL10","CXCL13","CXCL9","GNLY","GZMB","IFNG","IL12A","IL12B","IRF1","PRF1","STAT1","TBX21"))
  301. #Major Histocompatibility Complex Class I Gene
  302. MHC_class_I <- list(c("B2M","HLA-A","TAP1"))
  304. #anti-inflammatory genes
  305. Parainflammation <- list(c("CXCL10","PLAT","CCND1","LGMN","PLAUR","AIM2","MMP7","ICAM1","MX2","CXCL9","ANXA1","TLR2","PLA2G2D","ITGA2","MX1","HMOX1",
  306. "CD276","TIRAP","IL33","PTGES","TNFRSF12A","SCARB1","CD14","BLNK","IFIT3",'RETNLB',"IFIT2","ISG15","OAS2","REL","OAS3","CD44","PPARG","BST2","OAS1","NOX1","PLA2G2A","IFIT1",'IFITM3',"IL1RN"))
  308. #T-cell co-inhibitory genes
  309. T_cell_co_inhibition <- list(c("CD160","CD244","CD274","CTLA4","HAVCR2","LAG3","LAIR1","TIGIT"))
  311. #T-cell costimulation
  312. T_cell_co_stimulation <- list(c("CD2","CD226","CD27","CD28","CD40LG","ICOS","SLAMF1","TNFRSF18","TNFRSF25","TNFRSF4","TNFRSF8","TNFRSF9","TNFSF14"))
  314. #Type I interferon
  315. Type_I_IFN_Reponse <- list(c("DDX4","IFIT1","IFIT2","IFIT3","IRF7","ISG20","MX1","MX2","RSAD2","TNFSF10"))
  317. #Type II interferon
  318. Type_II_IFN_Reponse <- list(c("GPR146","SELP","AHR"))
  320. #T-cell status genes
  321. T_naiveness <- list(c("CCR6","CCR7","TCF7","SELL","LEF1","IL7R"))
  322. T_cytotoxicity <- list(c("GZMA","GZMB","GZMH","GZMK","PRF1","CXCL13","GNLY","IFNG","NKG7"))
  323. T_exhaustion <- list(c("PDCD1","CTLA4","HAVCR2","LAG3","TIGIT","LAYN"))
  324. T_proliferation <- list(c("MKI67","TOP2A","CD3D","CD3E"))
  326. #gene list
  327. library(readxl)
  328. genelist <- "./GeneList1.xlsX"
  329. genes <- read_excel(genelist)
  330. #Antigen Presentation and Processing Gene Set
  331. Antigen_Processing_and_Presentation <- subset(genes,Category=="Antigen_Processing_and_Presentation")
  332. Antigen_Processing_and_Presentation <- list(Antigen_Processing_and_Presentation$Symbol)
  334. #TCR signaling
  335. TCRsignalingPathway <- subset(genes,Category=="TCRsignalingPathway")
  336. TCRsignalingPathway <- list(TCRsignalingPathway$Symbol)
  338. #TNF (Tumor Necrosis Factor) gene set
  339. TNF_Family_Members <- subset(genes,Category==c("TNF_Family_Members"))
  340. TNF_Family_Members_Receptors <- subset(genes,Category==c("TNF_Family_Members_Receptors"))
  341. TNF_Family <- rbind(TNF_Family_Members,TNF_Family_Members_Receptors)
  342. TNF_Family <- list(TNF_Family)
  344. #Natural killer cell cytotoxicity
  345. NaturalKiller_Cell_Cytotoxicity <- subset(genes,Category=="NaturalKiller_Cell_Cytotoxicity")
  346. NaturalKiller_Cell_Cytotoxicity <- list(NaturalKiller_Cell_Cytotoxicity)
  348. #Cytokines
  349. Cytokines <- subset(genes,Category=="Cytokines")
  350. Cytokines <- list(Cytokines)
  352. #BCRSignalingPathway
  353. BCRSignalingPathway <- subset(genes,Category=="BCRSignalingPathway")
  354. BCRSignalingPathway <- list(BCRSignalingPathway)
  356. #Interleukins
  357. Interleukins <- subset(genes,Category=="Interleukins")
  358. Interleukins <- list(Interleukins)
  360. #TGFb_Family_Member
  361. TGFb_Family_Member <- subset(genes,Category=="TGFb_Family_Member")
  362. TGFb_Family_Member <- list(TGFb_Family_Member)
  364. #TIMELASER_marker
  365. TIMELASER_marker <- read.csv(file = "./TIMELASER_marker.csv")
  366. TIME_IA <- subset(TIMELASER_marker, TIMELASER.subtype=="TIME-IA")
  367. IA.gene <- list(TIME_IA$Gene)
  368. TIME_ISM <- subset(TIMELASER_marker, TIMELASER.subtype=="TIME-ISM")
  369. ISM.gene <- list(TIME_ISM$Gene)
  370. TIME_ISS <- subset(TIMELASER_marker, TIMELASER.subtype=="TIME-ISS")
  371. ISS.gene <- list(TIME_ISS$Gene)
  372. TIME_IE <- subset(TIMELASER_marker, TIMELASER.subtype=="TIME-IE")
  373. IE.gene <- list(TIME_IE$Gene)
  374. TIME_IR <- subset(TIMELASER_marker, TIMELASER.subtype=="TIME-IR")
  375. IR.gene <- list(TIME_IR$Gene)
  377. library(Seurat)
  378. library(ggplot2)
  379. library(dplyr)
  380. library(ggpubr)
  382. DefaultAssay(subsce.corr) <- "RNA"
  383. gene_sets <- list(
  384. IA = IA.gene,
  385. ISM = ISM.gene,
  386. ISS = ISS.gene,
  387. IE = IE.gene,
  388. IR = IR.gene,
  389. APC_co_inhibition = APC_co_inhibition,
  390. APC_co_stimulation = APC_co_stimulation,
  391. CCR = CCR,
  392. Check_point = Check_point,
  393. Cytolytic_activity = Cytolytic_activity,
  394. HLA = HLA,
  395. Inflammation_promoting = Inflammation_promoting,
  396. MHC_class_I = MHC_class_I,
  397. Parainflammation = Parainflammation,
  398. T_cell_co_inhibition = T_cell_co_inhibition,
  399. T_cell_co_stimulation = T_cell_co_stimulation,
  400. Type_I_IFN_Reponse = Type_I_IFN_Reponse,
  401. Type_II_IFN_Reponse = Type_II_IFN_Reponse,
  402. T_naiveness = T_naiveness,
  403. T_cytotoxicity = T_cytotoxicity,
  404. T_exhaustion = T_exhaustion,
  405. T_proliferation = T_proliferation,
  406. Antigen_Processing_and_Presentation = Antigen_Processing_and_Presentation,
  407. TCRsignalingPathway = TCRsignalingPathway,
  408. TNF_Family = TNF_Family,
  409. NaturalKiller_Cell_Cytotoxicity = NaturalKiller_Cell_Cytotoxicity,
  410. Cytokines = Cytokines,
  411. BCRSignalingPathway = BCRSignalingPathway,
  412. Interleukins = Interleukins,
  413. TGFb_Family_Member = TGFb_Family_Member
  414. )
  415. #AddModuleScore
  416. for (nm in names(gene_sets)) {
  417. subsce.corr <- AddModuleScore(
  418. subsce.corr,
  419. features = gene_sets[[nm]],
  420. ctrl = 100,
  421. name = "score"
  422. )
  423. }
  424. # score1~scoreN
  425. score_cols <- grep("^score", colnames([email hidden]), value = TRUE)
  426. colnames([email hidden])[match(score_cols, colnames([email hidden]))] <- names(gene_sets)
  427. plot_list <- list()
  428. mycolors <- c("#70cdbe","#80C1D7","#968ab7","#d586b3","#eb7e60")
  429. for (nm in names(gene_sets)) {
  431. df <- [email hidden] %>%
  432. mutate(cellmodule = factor(cellmodule))
  434. # Kruskal-Wallis
  435. kw <- kruskal.test(df[[nm]] ~ df$cellmodule)
  437. p_label <- paste0(
  438. "Kruskal-Wallis: χ²(", kw$parameter, ") = ",
  439. format(kw$statistic, digits = 3),
  440. ifelse(kw$p.value < 0.001,
  441. ", p < 0.001",
  442. paste0(", p = ", format(kw$p.value, digits = 3)))
  443. )
  444. p <- ggplot(df, aes(x = cellmodule, y = .data[[nm]], fill = cellmodule)) +
  445. geom_violin(width = 0.8, alpha = 0.85, scale = "width") +
  446. geom_boxplot(width = 0.15, fill = "white", outlier.size = 1) +
  447. scale_fill_manual(values = mycolors) +
  448. theme_bw() +
  449. annotate(
  450. "text",
  451. x = 3,
  452. y = Inf,
  453. label = p_label,
  454. vjust = 1.5,
  455. size = 4,
  456. fontface = "italic"
  457. ) +
  458. labs(title = nm)
  460. plot_list[[nm]] <- p
  461. }
  463. #----Figure 2E----
  464. #Determine the TIME subtype for each patient
  465. write.csv(prop.table(table(subsce.corr$subcluster,subsce.corr$orig.ident),1),file = "./patient.type.csv")
  466. ALL <- data.frame(read.csv(file="./NSCLC/Figure/figure2/patient.type.csv"))#,row.names = 1))
  467. ALL <- ALL %>% mutate(cellmodule = case_when(
  468. X %in% c("Bn_06_GZMB","pDC_03_GZMB","Plasma_09_IGHG_RRM2","Plasma_05_IGHG_IGHGP","Plasma_01_IGHA_MZB1","Plasma_02_IGHA_IGLC3","CD8T_02_GZMB",
  469. "CD8T_03_MKI67","CD8T_01_CCL5","Treg_FOXP3","B_02_HLA.DQB1","CD4T_02_GNB2L1","CD4T_06_ACSL5","CD8T_04_HLA.B") ~ "module1",
  470. X %in% c("Astrocyte_03_ROM1","Endothelial_02_PTMAP2","Fibroblast_01_Myo_HIGD1B","Fibroblast_02_Myo_COL10A1","Fibroblast_03_Myo_MKI67","B_04_PTMAP2","CD4T_04_SNORD3A",
  471. "Astrocyte_05_HLA.DRB5","Neurophil_01_GOS2","Astrocyte_04_SCGB3A1","Macrophage_06_C3","Astrocyte_01_CAV1","Astrocyte_02_MALAT1","Astrocyte_06_CCL5") ~ "module2",
  472. X %in% c("NK_01_FGFBP2","CD4T_03_SCGB3A1","Endothelial_05_IL1RL1","Mast_01_GNB2L1","Mast_03_PBK","cDC_01_CLEC10A","cDC_02_CD207","cDC_02_CD207","Macrophage_02_FCN1","Macrophage_04_SPP1",
  473. "Macrophage_09_MKI67","CD4T_05_SFTPB","Plasma_06_IGHG_SFTPB","Endothelial_09_S100A14","Macrophage_08_SCGB3A2","Fibroblast_04_Alveolar_DPT","Macrophage_07_IGKC") ~ "module3",
  474. X %in% c("CD4T_01_IL7R","B_03_CD3E","B_05_FYN","Fibroblast_05_Alveolar_MYH11","Endothelial_07_NOTCH3","B_01_HLA.DRA","Plasma_07_IGHG_GRIFIN","Plasma_08_IGHG_GUSBP11","Plasma_03_IGHA_GUSBP11","Plasma_04_IGHA_VPS13D") ~ "module4",
  475. X %in% c("Endothelial_03_VCAM1","Endothelial_08_CCL21","Endothelial_01_FCN3","Endothelial_04_DKK2","Fibroblast_06_Alveolar_MS4A7","Endothelial_06_CD68","Mast_02_CHIT1","Macrophage_03_MARCO","Macrophage_01_LGMN","Macrophage_05_CSTB") ~ "module5",
  476. ))
  478. which.max ()
  479. rownames(ALL) <- ALL$X
  480. ALL$X <-NULL
  482. library(stringr)
  483. module1 <- subset(ALL,cellmodule=="module1")
  484. module1$cellmodule <- NULL
  485. module_freq <- data.frame(colSums(module1))
  486. module2 <- subset(ALL,cellmodule=="module2")
  487. module2$cellmodule <- NULL
  488. module_freq[[2]] <- data.frame(colSums(module2))
  489. module3 <- subset(ALL,cellmodule=="module3")
  490. module3$cellmodule <- NULL
  491. module_freq[[3]] <- data.frame(colSums(module3))
  492. module4 <- subset(ALL,cellmodule=="module4")
  493. module4$cellmodule <- NULL
  494. module_freq[[4]] <- data.frame(colSums(module4))
  495. module5 <- subset(ALL,cellmodule=="module5")
  496. module5$cellmodule <- NULL
  497. module_freq[[5]] <- data.frame(colSums(module5))
  498. module_freq[[6]] <- data.frame(rowSums(module_freq))
  499. colnames(module_freq) <- c("module1","module2","module3","module4","module5","sum_freq")
  500. module_freq$module1 <- module_freq$module1/module_freq$sum_freq
  501. module_freq$module2 <- module_freq$module2/module_freq$sum_freq
  502. module_freq$module3 <- module_freq$module3/module_freq$sum_freq
  503. module_freq$module4 <- module_freq$module4/module_freq$sum_freq
  504. module_freq$module5 <- module_freq$module5/module_freq$sum_freq
  505. module_freq[[6]] <- NULL
  506. module_freq$patient <- rownames(module_freq)
  507. row_anno <- module_freq$patient
  508. row_anno <- data.frame(row_anno)
  509. colnames(row_anno) <- "patient"
  510. #row_anno$patient <- rownames(row_anno)
  511. #row_anno$patient <- NULL
  512. row_anno <- row_anno %>% mutate(mutation = case_when(
  513. patient %in% c("P1","P2","P3","P4","P5","P6","P7","P8","P9","P10","P11") ~ "EGFR",
  514. patient %in% c("P12","P13","P14","P15","P16","P17","P18","P19","P20","P21") ~ "EGFR-BM",
  515. patient %in% c("P22","P23") ~ "EGFR-co-mutation",
  516. patient %in% c("P24","P25","P26","P27","P28","P29","P30","P31","P32","P33","P34","P35") ~ "KRAS",
  517. patient %in% c("P36","P37","P38") ~ "KRAS-co-mutation",
  518. patient %in% c("P39") ~ "ALK",
  519. patient %in% c("P40","P41","P42") ~ "ROS1",
  520. patient %in% c("P43") ~ "TP53",
  521. patient %in% c("P44") ~ "MET-BM",
  522. patient %in% c("P45") ~ "HER2"))
  523. rownames(row_anno) <- row_anno$patient
  524. row_anno$mutation<- factor(row_anno$mutation)
  525. rownames(module_freq) <- module_freq$patient
  526. row_anno$patient <- NULL
  527. module_freq$patient <- NULL
  528. #factor(row_anno$mutation)
  529. row_anno <- data.frame(t(row_anno))
  530. module_freq <- data.frame(t(module_freq))
  531. #row_anno$mutation = list(c("EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR","EGFR-co-mutation","EGFR-co-mutation","EGFR-co-mutation","KRAS","KRAS","KRAS","KRAS","KRAS","KRAS","KRAS","KRAS","KRAS","KRAS","KRAS","KRAS","KRAS","KRAS-co-mutation","KRAS-co-mutation","KRAS-co-mutation","ALK","ROS1","ROS1","ROS1","TP53","MET","BRAF"))
  532. module_freq
  533. colSums(module_freq)
  534. row_anno <- data.frame(t(row_anno))
  535. row_anno$mutation <- factor(row_anno$mutation,levels = c("EGFR","EGFR-BM","EGFR-co-mutation","KRAS","KRAS-co-mutation","ALK","ROS1","TP53","MET-BM","HER2"))
  536. library("pheatmap")
  537. ann_colors=list(class=c(L='#009933',R='#CC33CC',F='#FDDCA9'))
  538. mycol2 = list(mutation = c("EGFR"="#c6b7d4","EGFR-BM"="#d44e26","EGFR-co-mutation"="#e3a264","KRAS"="#6fc2d0","KRAS-co-mutation"="#6f9abf","ALK"="#a5c49b","ROS1"="#7266ac","TP53"="#FF9966","MET-BM"="#d84986","HER2"="#2d588e"))
  539. #names(mycol2) <- unique(row_anno$mutation)
  540. module_freq <- data.frame(t(module_freq))
  541. patientorder <- c("P22","P43","P9","P1","P4","P25","P10","P27",
  542. "P18","P15","P38","P14","P44","P12","P13","P3",
  543. "P16","P42","P2","P23","P8","P41","P26","P20","P17","P45","P21","P36","P7","P11","P19","P39","P34","P5",
  544. "P31","P28","P30","P29","P35","P32",
  545. "P24","P40","P37","P33","P6")
  546. df_sorted <- module_freq[patientorder,]
  549. df_sorted <- data.frame(t(df_sorted))
  550. #color = colorRampPalette(c("white","#e6573f", "#d23229"))(100)
  551. color = colorRampPalette(c("white","#0172be"))(100)
  552. pheatmap(df_sorted,
  553. color = color,
  554. annotation_col = row_anno,
  555. annotation_colors = mycol2,
  556. display_numbers = F,
  557. border="white",
  558. cluster_cols = F,
  559. cluster_rows = F,
  560. treeheight_col = 40,
  561. treeheight_row = 45,
  562. border_color = NA)
  564. #----Figure 2F----
  565. #Plotting the Organizational Preference Map for the TIME Subtype
  566. subsce.corr <- readRDS(file = "./NSCLC/Figure/figure2/subsce.corr.rds")
  567. ROIE <- function(crosstab){
  568. ## Calculate the Ro/e value from the given crosstab
  569. ##
  570. ## Args:
  571. #' @crosstab: the contingency table of given distribution
  572. ##
  573. ## Return:
  574. ## The Ro/e matrix
  575. rowsum.matrix <- matrix(0, nrow = nrow(crosstab), ncol = ncol(crosstab))
  576. rowsum.matrix[,1] <- rowSums(crosstab)
  577. colsum.matrix <- matrix(0, nrow = ncol(crosstab), ncol = ncol(crosstab))
  578. colsum.matrix[1,] <- colSums(crosstab)
  579. allsum <- sum(crosstab)
  580. roie <- divMatrix(crosstab, rowsum.matrix %*% colsum.matrix / allsum)
  581. row.names(roie) <- row.names(crosstab)
  582. colnames(roie) <- colnames(crosstab)
  583. return(roie)
  584. }
  585. divMatrix <- function(m1, m2){
  586. ## Divide each element in turn in two same dimension matrixes
  587. ##
  588. ## Args:
  589. #' @m1: the first matrix
  590. #' @m2: the second matrix
  591. ##
  592. ## Returns:
  593. ## a matrix with the same dimension, row names and column names as m1.
  594. ## result[i,j] = m1[i,j] / m2[i,j]
  595. dim_m1 <- dim(m1)
  596. dim_m2 <- dim(m2)
  597. if( sum(dim_m1 == dim_m2) == 2 ){
  598. div.result <- matrix( rep(0,dim_m1[1] * dim_m1[2]) , nrow = dim_m1[1] )
  599. row.names(div.result) <- row.names(m1)
  600. colnames(div.result) <- colnames(m1)
  601. for(i in 1:dim_m1[1]){
  602. for(j in 1:dim_m1[2]){
  603. div.result[i,j] <- m1[i,j] / m2[i,j]
  604. }
  605. }
  606. return(div.result)
  607. }
  608. else{
  609. warning("The dimensions of m1 and m2 are different")
  610. }
  611. }
  612. metadata <- [email hidden]
  614. meta_filt <- metadata[metadata$mutation %in% c("ALK","EGFR","EGFR-BM","EGFR-co-mutation","HER2","KRAS","KRAS-co-mutation","MET-BM","ROS1","TP53"),]
  615. meta_filt$mutation <- factor(as.vector(meta_filt$mutation),levels=c("ALK","EGFR","EGFR-BM","EGFR-co-mutation","HER2","KRAS","KRAS-co-mutation","MET-BM","ROS1","TP53"))
  616. summary <- table(meta_filt[,c('cellmodule','mutation')])
  618. # ro/e
  619. roe <- as.data.frame(ROIE(summary))
  620. library("pheatmap")
  621. pheatmap(roe, display_numbers = TRUE,number_color = "black",cluster_row = FALSE,cluster_col = FALSE,
  622. color = colorRampPalette(c("#e7eb8e","#cd782d","firebrick3"))(50))
  624. #Create a lollipop chart showing organizational preferences for the TIME subtype
  625. roe0 <- as.data.frame(t(roe))
  626. roe0$group <- rownames(roe0)
  628. library(ggpubr)
  629. mycol2 = c("#c6b7d4","#d44e26","#e3a264","#6fc2d0","#6f9abf","#a5c49b","#7266ac","#FF9966","#d84986","#2d588e")
  630. #module1
  631. ggdotchart(roe0, x = "group", y = "module1",
  632. color = "group", # Color by groups
  633. sorting = "descending", # Sort value in descending order
  634. add = "segments", # Add segments from y = 0 to dots
  635. add.params = list(color = "lightgray", size = 2), # Change segment color and size
  636. dot.size = 11, # Large dot size
  637. label = round(roe0$module1,2), # Add mpg values as dot labels
  638. font.label = list(color = "black", size = 9,
  639. vjust = 0.5), # Adjust label parameters
  640. ggtheme = theme_pubr())+
  641. geom_hline(yintercept = 1, linetype = 2, color = "black")+
  642. theme(legend.position='none')+ylab("Ro/e")+ggtitle("Ro/e-module1")+
  643. theme(axis.text.x = element_text(angle = 45,vjust = 1))+
  644. scale_color_manual(values = c("EGFR"="#c6b7d4","EGFR-BM"="#d44e26","EGFR-co-mutation"="#e3a264","KRAS"="#6fc2d0","KRAS-co-mutation"="#6f9abf","ALK"="#a5c49b","ROS1"="#7266ac","TP53"="#FF9966","MET-BM"="#d84986","HER2"="#2d588e"))
  646. #module2
  647. ggdotchart(roe0, x = "group", y = "module2",
  648. color = "group", # Color by groups
  649. sorting = "descending", # Sort value in descending order
  650. add = "segments", # Add segments from y = 0 to dots
  651. add.params = list(color = "lightgray", size = 2), # Change segment color and size
  652. dot.size = 11, # Large dot size
  653. label = round(roe0$module2,2), # Add mpg values as dot labels
  654. font.label = list(color = "black", size = 9,
  655. vjust = 0.5), # Adjust label parameters
  656. ggtheme = theme_pubr())+
  657. geom_hline(yintercept = 1, linetype = 2, color = "black")+
  658. theme(legend.position='none')+ylab("Ro/e")+ggtitle("Ro/e-module2")+
  659. theme(axis.text.x = element_text(angle = 45,vjust = 1))+
  660. scale_color_manual(values = c("EGFR"="#c6b7d4","EGFR-BM"="#d44e26","EGFR-co-mutation"="#e3a264","KRAS"="#6fc2d0","KRAS-co-mutation"="#6f9abf","ALK"="#a5c49b","ROS1"="#7266ac","TP53"="#FF9966","MET-BM"="#d84986","HER2"="#2d588e"))
  662. #module3
  663. ggdotchart(roe0, x = "group", y = "module3",
  664. color = "group", # Color by groups
  665. sorting = "descending", # Sort value in descending order
  666. add = "segments", # Add segments from y = 0 to dots
  667. add.params = list(color = "lightgray", size = 2), # Change segment color and size
  668. dot.size = 11, # Large dot size
  669. label = round(roe0$module3,2), # Add mpg values as dot labels
  670. font.label = list(color = "black", size = 9,
  671. vjust = 0.5), # Adjust label parameters
  672. ggtheme = theme_pubr())+
  673. geom_hline(yintercept = 1, linetype = 2, color = "black")+
  674. theme(legend.position='none')+ylab("Ro/e")+ggtitle("Ro/e-module3")+
  675. theme(axis.text.x = element_text(angle = 45,vjust = 1))+
  676. scale_color_manual(values = c("EGFR"="#c6b7d4","EGFR-BM"="#d44e26","EGFR-co-mutation"="#e3a264","KRAS"="#6fc2d0","KRAS-co-mutation"="#6f9abf","ALK"="#a5c49b","ROS1"="#7266ac","TP53"="#FF9966","MET-BM"="#d84986","HER2"="#2d588e"))
  678. #module4
  679. ggdotchart(roe0, x = "group", y = "module4",
  680. color = "group", # Color by groups
  681. sorting = "descending", # Sort value in descending order
  682. add = "segments", # Add segments from y = 0 to dots
  683. add.params = list(color = "lightgray", size = 2), # Change segment color and size
  684. dot.size = 11, # Large dot size
  685. label = round(roe0$module4,2), # Add mpg values as dot labels
  686. font.label = list(color = "black", size = 9,
  687. vjust = 0.5), # Adjust label parameters
  688. ggtheme = theme_pubr())+
  689. geom_hline(yintercept = 1, linetype = 2, color = "black")+
  690. theme(legend.position='none')+ylab("Ro/e")+ggtitle("Ro/e-module4")+
  691. theme(axis.text.x = element_text(angle = 45,vjust = 1))+
  692. scale_color_manual(values = c("EGFR"="#c6b7d4","EGFR-BM"="#d44e26","EGFR-co-mutation"="#e3a264","KRAS"="#6fc2d0","KRAS-co-mutation"="#6f9abf","ALK"="#a5c49b","ROS1"="#7266ac","TP53"="#FF9966","MET-BM"="#d84986","HER2"="#2d588e"))
  694. #module5
  695. ggdotchart(roe0, x = "group", y = "module5",
  696. color = "group", # Color by groups
  697. sorting = "descending", # Sort value in descending order
  698. add = "segments", # Add segments from y = 0 to dots
  699. add.params = list(color = "lightgray", size = 2), # Change segment color and size
  700. dot.size = 11, # Large dot size
  701. label = round(roe0$module5,2), # Add mpg values as dot labels
  702. font.label = list(color = "black", size = 9,
  703. vjust = 0.5), # Adjust label parameters
  704. ggtheme = theme_pubr())+
  705. geom_hline(yintercept = 1, linetype = 2, color = "black")+
  706. theme(legend.position='none')+ylab("Ro/e")+ggtitle("Ro/e-module5")+
  707. theme(axis.text.x = element_text(angle = 45,vjust = 1))+
  708. scale_color_manual(values = c("EGFR"="#c6b7d4","EGFR-BM"="#d44e26","EGFR-co-mutation"="#e3a264","KRAS"="#6fc2d0","KRAS-co-mutation"="#6f9abf","ALK"="#a5c49b","ROS1"="#7266ac","TP53"="#FF9966","MET-BM"="#d84986","HER2"="#2d588e"))
  711. #Plotting a scatter plot of gene set scores
  712. library(ggplot2)
  713. library(ggsankey)
  714. subsce.corr <- readRDS("./NSCLC/Figure/figure2/subsce.corr.rds")
  715. DefaultAssay(subsce.corr) <- "RNA"
  716. immune.inhibitor <- list(c("ADORA2A","BTLA","CD160","CD244","CD274","CD96","CSF1R","CTLA4","HAVCR2","IDO1","IL10","IL10RB","KDR","KIR2DL1","KIR2DL3","LAG3","LGALS9","PDCD1","PDCD1LG2","TGFB1","TGFBR1","TIGIT","VTCN1"))
  717. Inscore <- AddModuleScore(subsce.corr,
  718. features = immune.inhibitor,
  719. ctrl = 100,
  720. name = "AUCELL")
  721. colnames([email hidden])
  722. colnames([email hidden])[28] <- 'immune.inhibitor'
  724. immune.stimulator <- list(c("CD27","CD276","CD28","CD40","CD40LG","CD48","CD70","CD80","CD86","CXCL12","CXCR4","ENTPD1","HHLA2","ICOS","ICOSLG","IL2RA","IL6","IL6R","KLRC1","KLRK1","LTA","MICB","NT5E","PVR","RAET1E","TMIGD2","TNFRSF13B","TNFRSF13C","TNFRSF14","TNFRSF17","TNFRSF18","TNFRSF25","TNFRSF4","TNFRSF8","TNFRSF9","TNFSF13","TNFSF13B","TNFSF14","TNFSF15","TNFSF18","TNFSF4","TNFSF9","ULBP1"))
  725. Inscore <- AddModuleScore(Inscore,
  726. features = immune.stimulator,
  727. ctrl = 100,
  728. name = "AUCELL")
  729. colnames([email hidden])
  730. colnames([email hidden])[29] <- 'immune.stimulator'
  732. Chemokine.receptor <- list(c("CCR1","CCR2","CCR3","CCR4","CCR5","CCR6","CCR7","CCR8","CCR9","CCR10", "CXCR1","CXCR2","CXCR3","CXCR4","CXCR5","CXCR6","XCR1","CX3R1"))
  733. Inscore <- AddModuleScore(Inscore,
  734. features = Chemokine.receptor,
  735. ctrl = 100,
  736. name = "AUCELL")
  737. colnames([email hidden])
  738. colnames([email hidden])[30] <- 'Chemokine.receptor'
  740. Chemokine <- list(c("CCL1","CCL2","CCL3","CCL4","CCL5","CCL7","CCL8","CCL11","CCL13","CCL14","CCL15","CCL16","CCL17","CCL18","CCL19","CCL20","CCL21","CCL22","CCL23","CCL24","CCL25","CCL26","CCL28","CX3CL1","CXCL1","CXCL2","CXCL3","CXCL5","CXCL6","CXCL8","CXCL9","CXCL10","CXCL11","CXCL12","CXCL13","CXCL14","CXCL16","CXCL17"))
  741. Inscore <- AddModuleScore(Inscore,
  742. features = Chemokine,
  743. ctrl = 100,
  744. name = "AUCELL")
  745. colnames([email hidden])
  746. colnames([email hidden])[31] <- 'Chemokine'
  748. #----Figure 2G----
  749. #The 14 malicious pathways in the decouple package are scored on CM1-CM5
  750. library(Seurat)
  751. library(decoupleR)
  752. library(dplyr)
  753. library(tibble)
  754. library(tidyverse)
  755. corsce <- readRDS("./yang/wang/NSCLC/subsce.corr.rds")
  756. net <- get_progeny(organism = 'human', top = 500)
  757. mat <- as.matrix(corsce@assays$RNA@layers$data)
  758. rownames(mat) <- rownames(corsce)
  759. colnames(mat) <- colnames(corsce)
  760. acts <- run_mlm(mat=mat, net=net, .source='source', .target='target',
  762. .mor='weight', minsize = 5)
  764. acts
  765. corsce[['pathwaysmlm']] <- acts %>%
  766. pivot_wider(id_cols = 'source', names_from = 'condition', values_from = 'score') %>%
  767. column_to_rownames('source') %>%
  768. Seurat::CreateAssayObject()
  770. DefaultAssay(corsce) <- 'pathwaysmlm'
  772. corsce <- ScaleData(corsce)
  773. corsce@assays$pathwaysmlm@data <- corsce@assays$[email hidden]
  774. #module
  775. Idents(corsce) <- "cellmodule"
  776. df <- t(as.matrix(corsce@assays$pathwaysmlm@data)) %>%
  777. as.data.frame() %>%
  778. mutate(cellmodule = Idents(corsce)) %>%
  779. pivot_longer(cols = -cellmodule, names_to = "pathway", values_to = "score") %>%
  780. group_by(cellmodule, pathway) %>%
  781. summarise(mean_score = mean(score), .groups = "drop")
  783. selected_pathways <- c("Androgen","EGFR","Estrogen","Hypoxia","JAK-STAT","MAPK","NFkB","p53","PI3K","TGFb","TNFa","Trail","VEGF","WNT")
  785. df_filtered <- df %>%
  786. filter(pathway %in% selected_pathways)
  788. top_acts_mat <- df_filtered %>%
  789. pivot_wider(names_from = pathway, values_from = mean_score) %>%
  790. column_to_rownames('cellmodule') %>%
  791. as.matrix()
  792. library(pheatmap)
  793. module.top_acts_mat <- top_acts_mat
  794. pheatmap(module.top_acts_mat,
  795. scale = "row",
  796. cluster_rows = TRUE,
  797. cluster_cols = TRUE,
  798. color = colorRampPalette(c("#A5CCE0", "white", "#d66c61"))(100),
  799. main = "Pathway Activity across Cell Modules")
  802. #mutation
  803. Idents(corsce) <- "mutation"
  804. df <- t(as.matrix(corsce@assays$pathwaysmlm@data)) %>%
  805. as.data.frame() %>%
  806. mutate(mutation = Idents(corsce)) %>%
  807. pivot_longer(cols = -mutation, names_to = "pathway", values_to = "score") %>%
  808. group_by(mutation, pathway) %>%
  809. summarise(mean_score = mean(score), .groups = "drop")
  811. selected_pathways <- c("Androgen","EGFR","Estrogen","Hypoxia","JAK-STAT","MAPK","NFkB","p53","PI3K","TGFb","TNFa","Trail","VEGF","WNT")
  813. df_filtered <- df %>%
  814. filter(pathway %in% selected_pathways)
  816. top_acts_mat <- df_filtered %>%
  817. pivot_wider(names_from = pathway, values_from = mean_score) %>%
  818. column_to_rownames('mutation') %>%
  819. as.matrix()
  820. library(pheatmap)
  822. pheatmap(top_acts_mat,
  823. scale = "row",
  824. cluster_rows = TRUE,
  825. cluster_cols = TRUE,
  826. border="white",
  827. border_color = NA,
  828. color = colorRampPalette(c("#A5CCE0", "white", "#d66c61"))(100),
  829. main = "Pathway Activity across Cell Modules")
  831. df.module <- data.frame(module.top_acts_mat)
  832. df.mutation <- data.frame(top_acts_mat)
  833. write.csv(df.module,file = "/data/yang/wang/NSCLC/df.module.csv")
  834. write.csv(df.mutation,file = "/data/yang/wang/NSCLC/df.mutation.csv")
  836. #----Figure 2H----
  837. #Drawing Sankey
  838. #Classify regulon_group expression values into high and low categories
  839. df <- [email hidden][,c(4,27:31)]
  840. df$score <- ifelse(df$immune.inhibitor >= median(df$immune.inhibitor, na.rm = TRUE), "High", "Low")
  841. df$score <- ifelse(df$immune.stimulator >= median(df$immune.stimulator, na.rm = TRUE), "High", "Low")
  842. df$score <- ifelse(df$Chemokine.receptor >= median(df$Chemokine.receptor, na.rm = TRUE), "High", "Low")
  843. df$score <- ifelse(df$Chemokine >= median(df$Chemokine, na.rm = TRUE), "High", "Low")
  844. write.csv(df,file = "./NSCLC/Figure/figure2/df.csv")
  845. library(ggsankey)
  846. library(tidyverse)
  848. level <- unique(df$score)
  849. mutation <- unique(df$mutation)
  850. module <- unique(df$cellmodule)
  851. level_colors <- c(
  852. "Low" = "#ADD4A5",
  853. "High" = "#C9B8DA"
  854. )
  855. mutation_colors <- c(
  856. "EGFR" = "#c6b7d4",
  857. "EGFR-BM" = "#d44e26",
  858. "EGFR-co-mutation" = "#e3a264",
  859. "KRAS" = "#6fc2d0",
  860. "KRAS-co-mutation" = "#6f9abf",
  861. "ALK" = "#a5c49b",
  862. "ROS1" = "#7266ac",
  863. "TP53" = "#FF9966",
  864. "MET-BM" = "#d84986",
  865. "HER2" = "#2d588e"
  866. )
  868. module_colors <- c(
  869. "module1" = "#70cdbe",
  870. "module2" = "#80C1D7",
  871. "module3" = "#968ab7",
  872. "module4" = "#d586b3",
  873. "module5" = "#eb7e60")
  874. level <- names(level_colors)
  875. mutation <- names(mutation_colors)
  877. module <- names(module_colors)
  879. color_palette <- c(module_colors,level_colors,mutation_colors)
  880. df$cellmodule <- factor(df$cellmodule,levels = c("module1","module2","module3","module4","module5"))
  881. df_sankey <- df %>%
  882. select(cellmodule,score, mutation) %>%
  883. make_long(cellmodule,score, mutation)
  885. ggplot(df_sankey, aes(x = x,
  886. next_x = next_x,
  887. node = node,
  888. next_node = next_node,
  889. fill = node,
  890. label = node)) +
  891. scale_fill_manual(values = color_palette) +
  892. geom_sankey(flow.alpha = 0.6, smooth = 6, width = 0.15) +
  893. geom_sankey_text(size = 3, color = "black") +
  894. theme_void() +
  895. theme(legend.position = "none")

Figure2.R at commit 44ddfe7, under MIT · at the source

Overview

Authors: Kuan Yang1, Kaiyue Yang1, Jiasi Wang1, Hang Zhao1, Wenqi Jiang1, Depeng Mu1, Xiao Peng1, Yiming Yan1, Xing Gao1, Jing Bai1, Congxue Hu1, Yunpeng Zhang1,2, Xia Li1,2
ORCID iDs: Xing Gao, Jing Bai, Yunpeng Zhang, Xia Li
  1. College of Bioinformatics Science and Technology, Harbin Medical University, Harbin 150081, China; (K.Y.); (K.Y.); (J.W.); (H.Z.); (W.J.); (D.M.); (X.P.); (Y.Y.); (X.G.); (J.B.); (C.H.)
  2. School of Intelligent Medicine and Technology, Big Data Research Center, Hainan Medical University, No. 3 Xueyuan Road, Longhua District, Haikou 571199, China
Institutions: Harbin Medical University (China); Hainan Medical University (China)
Journal: International journal of molecular sciences, volume 27, issue 9, article 3997
Dates: received 30 March 2026; accepted 27 April 2026; published online 29 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3390/ijms27093997 · PMID 42123577 · PMCID PMC13163347 · OpenAlex W7160146146
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), other condition (population)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Connectivity
Keywords: driver mutated, immune coordination, malignant regulatory network, prognostic biomarkers, single-cell transcriptomics
MeSH: Carcinoma, Non-Small-Cell Lung*, Lung Neoplasms*, Mutation*, Single-Cell Analysis*, Brain Neoplasms, Gene Expression Regulation, Neoplastic, Humans, Prognosis, Signal Transduction, Tumor Microenvironment (* major topic)
Topic: Cancer Immunotherapy and Biomarkers (Oncology, Medicine), according to OpenAlex
Funding: National Natural Science Foundation of China (62502128, 32570792, U23A20166, 62472131); National Science and Technology Major Program (2024ZD0530500); Heilongjiang Postdoctoral Fund (LBH-Z24210); China Postdoctoral Science Foundation (2024M760709); Longjiang New Era Outstanding Doctoral Dissertation Project Grant (LJYXL2024-069); Key Research and Development Program of Heilongjiang Province (2024ZX12C27); Harbin Medical University (2025-KYYWF-ZR0297)
Citations: not cited yet (Europe PMC); 94 references in the paper

Abstract

Oncogenic driver mutations in non-small cell lung cancer (NSCLC) activate defined signaling pathways that sustain tumor growth and influence the immune landscape. Yet, how coordinated interactions among diverse cell populations within the tumor immune microenvironment (TIME) contribute to this process remains largely unresolved. To address this, we profiled approximately 200,000 single cells from 45 treatment-naïve NSCLC patients representing seven major driver mutations. This analysis uncovered five multicellular modules (CM1–5) with distinct functional properties, each linked to specific malignant regulatory programs. Among them, CM2 and CM5 exhibited pronounced invasive features and were associated with unfavorable clinical outcomes. CM2 was predominantly observed in EGFR- and MET-driven brain metastases and was defined by strong crosstalk between astrocytes and myofibroblasts. Factors such as SPP1, PTN, and PSAP, together with metabolic alterations, contributed to a microenvironment supportive of metastatic colonization in the brain. By contrast, CM5 was enriched in ROS1-, KRAS-, and EGFR-mutant tumors and consisted of diverse myeloid and endothelial subsets characterized by immunosuppressive and pro-angiogenic signaling, including MIF, GALECTIN, and RETN, collectively facilitating immune escape and vascular remodeling. We further constructed and validated a driver mutation-specific prognostic signature (DMSP.sig) model integrating receptor–ligand interactions and core transcription factors, which effectively stratified patient survival. Leveraging this model, we also identified potential therapeutic candidates linked to these prognostic features, highlighting opportunities for clinical intervention. In summary, our study delineates how oncogenic drivers give rise to distinct TIME architectures, providing a framework for prognostic assessment and precision immunotherapy in high-risk NSCLC.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 40 matches between paragraphs and lines of code.

Zhangyunpeng1987/NSCLC-DMSPsig

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 44ddfe7e75bb24418c7eab6299ec053583865f57, 5 May 2026
Languages: R (13)
Size: 28 files, 13 scripts
Software Heritage: not archived
Found in: “Code Availability Statement”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (11 files), tidyverse (11 files), Seurat (9 files), ggpubr (7 files), pheatmap (7 files), survival (6 files), clusterProfiler (4 files), patchwork (4 files), circlize (3 files), ComplexHeatmap (3 files), cowplot (2 files), data.table (2 files), psych (2 files), reshape2 (2 files), rstatix (2 files), DESeq2 (1 file), glmnet (1 file), Harmony (1 file), igraph (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
15 files
At the source: github.com/Zhangyunpeng1987/NSCLC-DMSPsig

Code Availability Statement

The R scripts used for data processing, analysis, and visualization in this study are publicly available at: https://github.com/Zhangyunpeng1987/NSCLC-DMSPsig (accessed on 24 April 2026). The repository includes scripts corresponding to all figures presented in the manuscript and is organized by figure to facilitate reproducibility. A session information file is also provided to ensure transparency of the computational environment.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 13 scripts, each with its path and the digest of its content;
  • 40 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability Statement

All data used in this study are publicly available from GEO and TCGA. The datasets analyzed are existing publicly accessible datasets, and detailed information is provided in Table S1. The scRNA-seq datasets used in this study include GSE171145 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE171145), GSE131907 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE131907), GSE148071 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE148071), GSE202371 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE202371), and GSE136246 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE136246). In addition, we included scRNA-seq data from the study with PMID: 35027529, as well as data from the Code Ocean capsule (CO.0121060.V1). Spatial transcriptomics data were obtained from the publicly available EMBL-EBI database, including dataset E-MTAB-13530.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 13 authors, 5 keywords, 10 MeSH terms, 7 funders, 94 references.

Cite

This paper

Yang, K., Yang, K., Wang, J., Zhao, H., Jiang, W., Mu, D., Peng, X., Yan, Y., Gao, X., Bai, J., Hu, C., Zhang, Y., & Li, X. (2026). Coordinated Multicellular Immune Programs and Drug Targets Revealed by Single-Cell Analysis in Driver-Mutated NSCLC. International journal of molecular sciences, 27(9), 3997. https://doi.org/10.3390/ijms27093997

BibTeX

@article{yang2026coordinated,
author = {Yang, Kuan and Yang, Kaiyue and Wang, Jiasi and Zhao, Hang and Jiang, Wenqi and Mu, Depeng and Peng, Xiao and Yan, Yiming and Gao, Xing and Bai, Jing and Hu, Congxue and Zhang, Yunpeng and Li, Xia},
title = {{Coordinated Multicellular Immune Programs and Drug Targets Revealed by Single-Cell Analysis in Driver-Mutated NSCLC}},
journal = {International journal of molecular sciences},
year = {2026},
month = apr,
volume = {27},
number = {9},
pages = {3997},
publisher = {Multidisciplinary Digital Publishing Institute (MDPI)},
issn = {1422-0067},
doi = {10.3390/ijms27093997},
url = {https://doi.org/10.3390/ijms27093997},
pmid = {42123577},
pmcid = {PMC13163347}
}

RIS

TY - JOUR
AU - Yang, Kuan
AU - Yang, Kaiyue
AU - Wang, Jiasi
AU - Zhao, Hang
AU - Jiang, Wenqi
AU - Mu, Depeng
AU - Peng, Xiao
AU - Yan, Yiming
AU - Gao, Xing
AU - Bai, Jing
AU - Hu, Congxue
AU - Zhang, Yunpeng
AU - Li, Xia
TI - Coordinated Multicellular Immune Programs and Drug Targets Revealed by Single-Cell Analysis in Driver-Mutated NSCLC
T2 - International journal of molecular sciences
J2 - Int J Mol Sci
PY - 2026
DA - 2026/04/29
VL - 27
IS - 9
SP - 3997
SN - 1422-0067
PB - Multidisciplinary Digital Publishing Institute (MDPI)
DO - 10.3390/ijms27093997
UR - https://doi.org/10.3390/ijms27093997
LA - en
ER -

CSL-JSON

{
"id": "10.3390/ijms27093997",
"type": "article-journal",
"title": "Coordinated Multicellular Immune Programs and Drug Targets Revealed by Single-Cell Analysis in Driver-Mutated NSCLC",
"container-title": "International journal of molecular sciences",
"author": [
{
"family": "Yang",
"given": "Kuan"
},
{
"family": "Yang",
"given": "Kaiyue"
},
{
"family": "Wang",
"given": "Jiasi"
},
{
"family": "Zhao",
"given": "Hang"
},
{
"family": "Jiang",
"given": "Wenqi"
},
{
"family": "Mu",
"given": "Depeng"
},
{
"family": "Peng",
"given": "Xiao"
},
{
"family": "Yan",
"given": "Yiming"
},
{
"family": "Gao",
"given": "Xing"
},
{
"family": "Bai",
"given": "Jing"
},
{
"family": "Hu",
"given": "Congxue"
},
{
"family": "Zhang",
"given": "Yunpeng"
},
{
"family": "Li",
"given": "Xia"
}
],
"container-title-short": "Int J Mol Sci",
"volume": "27",
"issue": "9",
"page": "3997",
"DOI": "10.3390/ijms27093997",
"PMID": "42123577",
"PMCID": "PMC13163347",
"ISSN": "1422-0067",
"publisher": "Multidisciplinary Digital Publishing Institute (MDPI)",
"URL": "https://doi.org/10.3390/ijms27093997",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
29
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.3390/ijms27104466 [code]
Uncovering the Key Circuit FOSL2/FOS/EGR3/EGR1, Contributing to the Hyperexcitability of Excitatory Neurons in the Epileptic Temporal Cortex and Hippocampus.
Journal: International journal of molecular sciences
In common: glmnet, igraph, circlize, 11 other tools, genetics / omics, 2 authors
[2] doi:10.1016/j.isci.2026.115657 [code]
Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug.
Journal: iScience
In common: glmnet, survival, Harmony, 15 other tools, other condition
[3] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Harmony, rstatix, DESeq2, 13 other tools, genetics / omics
[4] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: survival, rstatix, DESeq2, 12 other tools, genetics / omics, other condition
[5] doi:10.1126/sciadv.aeg3223 [code]
The extreme diversity of retinal amacrine cells has deep evolutionary roots.
Journal: Science advances
In common: glmnet, psych, rstatix, 12 other tools, genetics / omics
[6] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: glmnet, survival, DESeq2, 12 other tools
[7] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Harmony, DESeq2, igraph, 12 other tools, genetics / omics, other condition
[8] doi:10.1038/s41593-026-02367-0 [code]
A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: Harmony, DESeq2, igraph, 12 other tools
[9] doi:10.1038/s41586-026-10629-x [code]
Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: Harmony, rstatix, DESeq2, 11 other tools, genetics / omics
[10] doi:10.1038/s41593-026-02300-5 [code]
Integrated single-cell and spatial transcriptomic profiling in ALS uncovers peripheral-to-central immune infiltration and reprogramming.
Journal: Nature neuroscience
In common: Harmony, DESeq2, igraph, 11 other tools, genetics / omics, other condition

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.