OSCR

Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug.

Code ↔ Paper

27 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 27 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § STAR★Methods › Method details › Unsupervised clustering analysis ↔ workline2.R, lines 606–650 · score 0.94 · ConsensusClusterPlus, clusterAlg, pFeature, pItem, maxK, inner
  2. [2] § Results › Immune infiltration differences in gliomas with different disulfidptosis patterns ↔ workline2.R, lines 1546–1581 · score 0.90 · innate immunity, immune infiltration, tumor cells, NK cells, granulocytes, priming
  3. [3] § STAR★Methods › Method details › Immune microenvironment analysis ↔ R/deconvo_tme.R, lines 1–48 · score 0.87 · MCP counter, quanTIseq, xCell, IOBR package, immune cell, EPIC
  4. [4] § STAR★Methods › Method details › Single-cell RNA-seq and in silico perturbation analysis ↔ R/scTenifoldKnk.R, lines 1–110 · score 0.86 · silico knockout, gene regulatory networks, Quality control, scTenifoldKnk, single cell, mitochondrial
  5. [5] § STAR★Methods › Method details › Data collection and preprocessing ↔ workline2.R, lines 418–502 · score 0.83 · CD2AP, SLC7A11, ACTB, DSTN, TLN1, MYL6
  6. [6] § Results › Transcriptional and genetic characteristics of 15 disulfidptosis genes ↔ workline2.R, lines 2–40 · score 0.81 · Missense mutations, CD2AP, SLC7A11, ACTB, DSTN, TLN1
  7. [7] § STAR★Methods › Quantification and statistical analysis ↔ R/batch_kruskal.R, lines 1–60 · score 0.80 · Kruskal Wallis, Benjamini Hochberg, chi squared, rank sum, variables, validation
  8. [8] § STAR★Methods › Method details › Single-cell RNA-seq and in silico perturbation analysis ↔ Supplementary analysis.R, lines 94–134 · score 0.78 · Cell cycle scoring, variable feature, Seurat, Quality control, Harmony, RNA
  9. [9] § Results › Immune infiltration characteristics associated with DisulfidpScore scores in gliomas ↔ R/data_doc.R, lines 51–116 · score 0.74 · epithelial mesenchymal transition, immune checkpoint, immune microenvironment, stromal, activity, cycle
  10. [10] § Results › Deciphering disulfidptosis and IQGAP1-driven networks in gliomas at single-cell resolution ↔ R/scTenifoldKnk.R, lines 1–110 · score 0.73 · silico knockout, gene regulatory network, perturbed genes, scTenifoldKnk, cell
  11. [11] § STAR★Methods › Method details › Construction and validation of the prognostic model ↔ R/BinomialModel.R, lines 1–115 · score 0.72 · cross validation, Model Construction, 0–1, LASSO, splitting, Ridge
  12. [12] § Results › Deciphering disulfidptosis and IQGAP1-driven networks in gliomas at single-cell resolution ↔ Supplementary analysis.R, lines 136–184 · score 0.70 · cell cycle phases, SingleR, UMAP, atlas, LGG, IQGAP1
  13. [13] § STAR★Methods › Method details › Pathway enrichment analysis ↔ workline2.R, lines 1269–1341 · score 0.65 · c2 cp kegg, v2023, MSigDB, Hs, GSEA, symbols
  14. [14] § STAR★Methods › Method details › Construction and validation of the prognostic model ↔ R/PrognosticModel.R, lines 1–117 · score 0.64 · cross validation, LASSO, splitting, Ridge, CV, fold
  15. [15] § Results › Unsupervised clustering to identify different disulfidptosis patterns in gliomas ↔ workline2.R, lines 606–650 · score 0.64 · CD2AP, remove batch, MYL6, ACTN4, CAPZB, INF2
  16. [16] § Results › Signaling pathway differences in gliomas with different disulfidptosis patterns ↔ R/data_doc.R, lines 51–116 · score 0.64 · epithelial mesenchymal transition, Cell cycle, checkpoint, pathways, tumor, regulation
  17. [17] § STAR★Methods › Method details › Pathway enrichment analysis ↔ R/sig_gsea.R, lines 245–295 · score 0.62 · MSigDB, clusterProfiler, Hs, GSEA, database, symbols
  18. [18] § Results › Transcriptional and genetic characteristics of 15 disulfidptosis genes ↔ workline2.R, lines 2–40 · score 0.61 · Missense mutations, SLC7A11, COAD, SNVs, FLNB, FLNA
  19. [19] § STAR★Methods › Quantification and statistical analysis ↔ R/batch_wilcoxon.R, lines 1–59 · score 0.58 · Benjamini Hochberg, Wilcoxon rank sum, variables, validation
  20. [20] § STAR★Methods › Method details › Immune microenvironment analysis ↔ R/iobr_deconvo_pipeline.R, the whole file · a weak match · score 0.58 · IOBR, MCP, xCELL, EPIC, IPS, QuantiSeq
  21. [21] § Results › Model construction based on the disulfidptosis patterns ↔ R/BinomialModel.R, lines 1–115 · score 0.57 · cross validation, model construction, LASSO, Ridge, fold, training
  22. [22] § Results › Predictive ability of DisulfidpScore for chemotherapy drug sensitivity in gliomas ↔ pRRophetic/R/predict_from_cgp.R, lines 146–181 · score 0.56 · drug response, drug sensitivity, AUY922, bortezomib, elesclomol, IC50
  23. [23] § Results › Immune infiltration characteristics associated with DisulfidpScore scores in gliomas ↔ workline2.R, lines 2196–2272 · score 0.55 · cancer immunity cycle, immune functions, DisulfidpScore score, cell
  24. [24] § Results › Predictive ability of DisulfidpScore ↔ R/PrognosticModel.R, lines 239–292 · score 0.53 · dependent ROC curves, predictive accuracy, AUC, prognostic
  25. [25] § STAR★Methods › Method details › Immunotherapy prediction analysis ↔ R/batch_wilcoxon.R, lines 1–59 · score 0.53 · Benjamini Hochberg, Wilcoxon rank sum, validate
  26. [26] § Results › Signaling pathway differences in gliomas with different disulfidptosis patterns ↔ R/sig_gsea.R, lines 1–113 · score 0.53 · conducted gene, Cell cycle, KRAS, hallmark, pathways, regulation
  27. [27] § STAR★Methods › Method details › Construction and validation of the prognostic model ↔ R/PrognosticModel.R, lines 1–117 · score 0.50 · predictive accuracy, ROC, prognostic, curves, regression, Cox

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 · 3,168 lines · 101 KB · no license · 8 matches

  2. #1.1 15个DRG的SNV分析####
  3. library(maftools)
  4. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.1 genomic")
  5. load("TCGA-LGG_maf.rdata")
  6. barcode=as.data.frame(unique(data$Tumor_Sample_Barcode))
  7. colnames(barcode)="sample"
  8. data1=data
  10. load("TCGA-GBM_maf.rdata")
  11. barcode2=as.data.frame(unique(data$Tumor_Sample_Barcode))
  12. colnames(barcode2)="sample"
  13. data2=data
  15. #合并lgg和gbm
  16. data <- rbind(data1,data2)
  17. maf.coad <- data
  18. class(maf.coad) #data.frame
  19. dim(maf.coad) #87957 141
  21. #匹配结果
  22. setwd("/home/data/t060432/HRT/Disulfidptosis/data/1.1 SNV&CNV/snv")
  23. clin <- read.csv("Grade_SNV_barcode.csv",header=TRUE) #660
  24. clin <- clin[-c(434,514),] #去除2个异常样本 TCGA-06-5416-01A-01D-1486-08,TCGA-DU-6392-01A-11D-1705-08
  26. maf.coad=merge(data,clin,by="Tumor_Sample_Barcode") #59090 142
  27. sample=unique(maf.coad$Tumor_Sample_Barcode)
  29. maf <- read.maf(maf.coad,clinicalData = clin)
  31. #瀑布图
  32. gene=c("ACTB","ACTN4","CAPZB","CD2AP","DSTN","FLNA","FLNB","INF2",
  33. "IQGAP1","MYH10","MYH9","MYL6","PDLIM1","SLC7A11","TLN1")
  35. vc_cols = c("#8DD3C7","#80B1D3","#B2DF8A","#33A02C","#FB9A99","#E31A1C","#FDBF6F","#FF7F00",'#CAB2D6')
  36. names(vc_cols) = c('Multi_Hit','Missense_Mutation','Frame_Shift_Del','Nonsense_Mutation',
  37. 'Frame_Shift_Ins','In_Frame_Ins','Splice_Site','In_Frame_Del','Translation_Start_Site')
  38. oncoplot(maf = maf, genes=gene,
  39. draw_titv = F,fontSize = 0.75 ,
  40. colors = vc_cols,bgCol = "transparent")
  42. #互斥/共现
  43. par(oma = c(3, 4, 5, 1))
  44. somaticInteractions(maf = maf,
  45. genes=gene,
  46. pvalue = c(0.05, 0.5),
  47. colPal = "PiYG")
  51. #1.2 15个DRG的CNV分析####
  52. setwd("/home/data/t060432/HRT/Disulfidptosis/data/1.1 SNV&CNV/cnv")
  53. lgg=fread("TCGA-LGG.gistic.tsv",header = T, sep = '\t',data.table = F)
  54. gbm=fread("TCGA-GBM.gistic.tsv",header = T, sep = '\t',data.table = F)
  55. a=unique(colnames(gbm))
  56. gbm=gbm[,a] #614
  58. #基因注释
  59. LGG <- fread('gencode.v22.annotation.gene.probeMap',data.table = F)%>%
  60. select(id,gene)%>%
  61. inner_join(lgg,by=c('id'='Gene Symbol'))%>%
  62. select(-id)%>%
  63. group_by(gene)%>%
  64. summarise_all(mean)%>%
  65. column_to_rownames('gene')
  67. GBM <- fread('gencode.v22.annotation.gene.probeMap',data.table = F)%>%
  68. select(id,gene)%>%
  69. inner_join(gbm,by=c('id'='Gene Symbol'))%>%
  70. select(-id)%>%
  71. group_by(gene)%>%
  72. summarise_all(mean)%>%
  73. column_to_rownames('gene')
  75. #合并
  76. identical(rownames(LGG),rownames(GBM))
  77. gliomas=cbind(LGG,GBM)
  79. gene=c("ACTB","ACTN4","CAPZB","CD2AP","DSTN","FLNA","FLNB","INF2",
  80. "IQGAP1","MYH10","MYH9","MYL6","PDLIM1","SLC7A11","TLN1")
  82. cnvdat <- gliomas[rownames(gliomas) %in% gene,]
  84. #替换为字符串
  85. cnvdat[cnvdat ==1] <- "gain"
  86. cnvdat[cnvdat ==-1] <- "loss"
  87. cnvdat[cnvdat ==0] <- "neutral"
  89. #统计gain和loss
  90. cnvdat <- as.data.frame(t(cnvdat))
  91. cnv <- as.data.frame(t(apply(cnvdat,2,table)))
  93. #绘图
  94. dat <- cnv
  95. dat$Gene <- rownames(dat)
  96. library(ggalt)
  98. #method1
  99. p <- ggplot(aes(x=loss,xend=gain,y=Gene),data=dat)+
  100. geom_dumbbell(colour_x = "green",colour_xend = "red",size_x = 2,size_xend = 2,size=0.5,color="gray",dot_guide = T)+
  101. theme_light()+theme(panel.grid.minor.x =element_blank(),
  102. panel.grid = element_blank(),
  103. legend.position = c("top")
  104. )+ xlab("CNV.frequency(%)")
  105. p
  109. #1.3 15个DRG的染色体圈图####
  110. library(RCircos)
  112. #导入人类染色体数据
  113. data(UCSC.HG38.Human.CytoBandIdeogram)
  114. head(UCSC.HG38.Human.CytoBandIdeogram)
  116. #构建RCircos的core components
  117. cyto.info <- UCSC.HG38.Human.CytoBandIdeogram
  118. RCircos.Set.Core.Components(cyto.info, chr.exclude=NULL,tracks.inside=6, tracks.outside=0)
  120. #染色体图
  121. RCircos.Set.Plot.Area()
  122. RCircos.Chromosome.Ideogram.Plot()
  124. #染色体位置信息在genemap22
  125. setwd("/home/data/t060432/HRT/Disulfidptosis/data/1.3 circos")
  126. data=read.table("15DRG-cricos.txt", head = T)
  128. #将基因与染色体连接
  129. name.col <- 4
  130. side <- "in"
  131. track.num <- 1
  132. RCircos.Gene.Connector.Plot(data, track.num, side)
  134. track.num <- 2
  135. RCircos.Gene.Name.Plot(data, name.col,track.num, side)
  139. #1.4 15个DRG的共表达网络####
  140. library(tidyverse)
  141. library(ggplot2)
  142. library(igraph)
  143. library(ggraph)
  144. library(RColorBrewer)
  145. library(tidygraph)
  146. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  147. load(file="TCGA.RData")
  148. gene=c("ACTB","ACTN4","CAPZB","CD2AP","DSTN","FLNA","FLNB","INF2",
  149. "IQGAP1","MYH10","MYH9","MYL6","PDLIM1","SLC7A11","TLN1")
  150. gene=TCGA[gene,]
  151. gene=as.data.frame(t(gene)) #691
  153. M = cor(gene)
  155. #计算r和p
  156. data=gene
  157. m1=data[,c(1:15)]
  158. m2=data[,c(1:15)]
  159. cor_2_matrix <- function(m1,m2){
  160. apply(m2 , 2, function(x){
  161. unlist(apply(m1, 2,function(y){
  162. cor(as.numeric(x),
  163. as.numeric(y))
  164. }))
  165. })
  166. }
  167. rdf=cor_2_matrix( m1 , m2 ) %>% as.data.frame()
  168. rdf$gene=rownames(rdf)
  169. rdf <- rdf %>% gather(key = 'soure',value = 'r',-gene)
  171. corP_2_matrix <- function(m1,m2){
  172. apply(m2 , 2, function(x){
  173. unlist(apply(m1, 2,function(y){
  174. cor.test(as.numeric(x),
  175. as.numeric(y))$p.value
  176. }))
  177. })
  178. }
  179. pdf=corP_2_matrix( m1 , m2 ) %>% as.data.frame()
  180. pdf$gene=rownames(pdf)
  181. pdf <- pdf %>% gather(key = 'source',value = 'p',-gene)
  183. Toal <- bind_cols(pdf,rdf) %>%.[-c(4,5)]
  184. Toal2 <- filter(Toal,p<0.0001)
  185. colnames(Toal2) <- c('to','from','pvalue','corr')
  187. #加上相关性信息
  188. Toal2$relation=ifelse(Toal2$corr>0,'Positive correlation with P < 0.0001','Negative correlation with P < 0.0001')
  189. #去除相关性为1的行
  190. Toal3 = filter(Toal2, Toal2$corr !='1')
  191. table(Toal3$relation)
  192. #Negative correlation with P < 0.0001 Positive correlation with P < 0.0001
  193. #22 130
  195. #加载cox回归结果
  196. setwd("/home/data/t060432/HRT/Disulfidptosis/data/1.4 co-network")
  197. Toal4 = read.table(file="15DRG-cor-network.txt",header=T)
  199. m_data=Toal3
  200. #节点数据
  201. nodes <- data.frame(name = unique(union(m_data$from, m_data$to)))
  202. identical(Toal4$Gene,nodes$name)
  203. nodes$survival_impact <- Toal4$Cox_test_pvalue
  204. nodes$cluster <- c(rep("Disulfidptosis",15))
  205. nodes$role_type <- Toal4$type
  207. #边数据
  208. edges <- m_data[c("from","to","corr")]
  209. colnames(edges)[3]="Pearson_R"
  210. edges$class <- ifelse(edges$Pearson_R>0, "Positive correlation with P < 0.0001",
  211. "Negative correlation with P < 0.0001")
  212. g <- tbl_graph(nodes = nodes, edges = edges)
  213. class(g)
  215. #绘制图形
  216. colors <- "white"
  218. ggraph(g,layout='linear',circular = TRUE) +
  219. geom_node_point(aes(size=survival_impact,colour = role_type),
  220. alpha = 0.8) +
  221. geom_node_text(aes(x = x*1.15, y=y*1.15, label=name,color=role_type),
  222. angle=0,hjust=0, fontface="bold",size=2.5,family="Times") +
  223. scale_size_continuous(range = c(16, 8)) +
  224. geom_node_point(size = 4,aes(colour = cluster))+
  225. scale_color_manual(values = c(colors,"#0088FF","#FF0033")) +
  226. geom_edge_arc(mapping = aes(edge_width = abs(Pearson_R),
  227. edge_color = class),
  228. strength = 0.02,alpha = 0.6) +
  229. scale_edge_colour_manual(values = c("#abd9e9", "#fec8c9")) +
  230. scale_edge_width_continuous(range = c(0.5,3)) +
  231. theme_graph()
  235. #1.5 15个DRG的单因素cox回归####
  236. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  237. library(forestplot)
  238. library(grid)
  239. library(magrittr)
  240. library(checkmate)
  241. library(data.table)
  242. library(survival)
  243. library(survminer)
  245. gene=c("ACTB","ACTN4","CAPZB","CD2AP","DSTN","FLNA","FLNB","INF2",
  246. "IQGAP1","MYH10","MYH9","MYL6","PDLIM1","SLC7A11","TLN1")
  248. #TCGA
  249. load(file="TCGA.RData")
  250. dat=TCGA[gene,]
  251. dat=as.data.frame(t(dat))
  253. identical(rownames(dat),rownames(OS))
  254. data=cbind(OS,dat)
  256. univar_out = data.frame(matrix(NA,15,5))
  257. rownames(univar_out) = colnames(data)[-(1:2)]
  258. colnames(univar_out) = c("Coeffcient","HR","lower .95","upper .95","P-value")
  259. cox_data = data[,c(3:ncol(data))]
  260. cox_data = cbind(data[,c(1:2)],cox_data)
  261. str(cox_data)
  262. for(i in colnames(cox_data)[-(1:2)]){
  263. cox = coxph(Surv(OS.time, OS) ~ cox_data[,i], data = cox_data)
  264. cox_summ = summary(cox)
  265. univar_out[i,1] = cox_summ$coefficients[,1]
  266. univar_out[i,2] = cox_summ$coefficients[,2]
  267. univar_out[i,3] = cox_summ$conf.int[,3]
  268. univar_out[i,4] = cox_summ$conf.int[,4]
  269. univar_out[i,5] = cox_summ$coefficients[,5]
  270. }
  271. univar_out_0.05 = univar_out[univar_out[,5] < 0.05,]
  272. univar_out_TCGA=univar_out[,c(2,5)]
  274. #可视化
  275. uni <- univar_out
  276. uni[,1:4] <- round(uni[,1:4],digits = 3)
  277. uni_tabletext <- data.frame(matrix(NA,(nrow(uni)),3))
  278. for(i in 1 : nrow(uni)){
  279. uni_tabletext[i,1] = rownames(uni)[i]
  280. uni_tabletext[i,2] = paste(uni[i,2],"(",uni[i,3],"-",uni[i,4],")",sep = "")
  281. ifelse(uni[i,5]<0.001,
  282. uni_tabletext[i,3]<-"<0.001",
  283. uni_tabletext[i,3]<-round(uni[i,5],digits = 3))}
  284. uni_tabletext_title <- rbind(c("","HR","P-value"),uni_tabletext)
  286. #森林图
  287. forestplot(uni_tabletext_title,
  288. mean = c(NA,univar_out[,2]),
  289. graph.pos=4,
  290. upper = c(NA,univar_out[,4]),
  291. lower = c(NA,univar_out[,3]),
  292. align = "c",
  293. boxsize=0.4,
  294. zero=1,
  295. lineheight = "auto",
  296. colgap=unit(8,"mm"),
  297. ci.vertices=TRUE,
  298. xticks = c(0,1,2,4,8),
  299. title="Univariate analysis",
  300. col = fpColors(box = "#7FBC41",lines = "black",zero = "grey"))
  306. #1.6 15个DRG的表达水平####
  307. #总体
  308. library(ggpubr)
  309. library(ggplot2)
  310. library(ggsignif)
  311. library(ggdist)
  312. #1.加载数并处理
  313. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  314. load("GTEx-f-t.RData")
  315. load("TCGA.RData")
  317. gene=c("ACTB","ACTN4","CAPZB","CD2AP","DSTN","FLNA","FLNB","INF2",
  318. "IQGAP1","MYH10","MYH9","MYL6","PDLIM1","SLC7A11","TLN1")
  319. gene <- as.vector(gene)
  321. #合并
  322. normal=tpm[gene,]
  323. gliomas=TCGA[gene,]
  324. Exp=cbind(normal,gliomas)
  326. Exp_plot <- as.data.frame(t(Exp))
  329. #2.加载样本信息
  330. sample=rownames(Exp_plot)
  331. group=c(rep("normal",290),rep("LGG",524),rep("GBM",167))
  332. info <- as.data.frame(cbind(sample,group))
  333. colnames(info)=c("Sample","Type")
  334. Exp_plot$sam=info$Type
  335. Exp_plot$sam <- factor(Exp_plot$sam, levels = c("normal","LGG","GBM"))
  336. expr_use <- na.omit(Exp_plot)
  337. expr_use <-as.data.frame(expr_use)
  338. expr_use_long <- gather(expr_use, gene, Expression, -sam)
  339. table(expr_use_long$gene)
  340. colnames(expr_use_long) <- c("Group","gene","Expression")
  342. #3.绘图
  343. p=ggboxplot(expr_use_long, x="gene", y="Expression", color = "black", fill="Group",
  344. ylab="Gene expression",
  345. xlab="",
  346. legend.title=NULL,
  347. palette = c("#FF0033","#009934","#0088FF"),
  348. width=0.6, add = "none")
  349. p1=p+stat_compare_means(aes(group=Group),
  350. method="wilcox.test",
  351. symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", " ")),
  352. label = "p.signif")
  353. p1
  356. #疾病分期
  357. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  358. load("TCGA.RData")
  360. table(pheno$grade)
  361. #G2 G3 G4
  362. #214 236 167
  363. a=pheno[,c(1,3)]
  364. a=na.omit(a)
  365. sample=a$sample #617
  367. #1.加载数并处理
  368. Exp <- as.data.frame(t(TCGA[,sample]))
  369. #gene=c("ACTB","ACTN4","CAPZB","CD2AP","DSTN","FLNA","FLNB","INF2",
  370. # "IQGAP1","MYH10","MYH9","MYL6","PDLIM1","SLC7A11","TLN1")
  371. gene="TLN1"
  372. gene <- as.vector(gene)
  373. Exp_plot <- as.data.frame(Exp[,gene])
  374. colnames(Exp_plot) = "TLN1"
  375. rownames(Exp_plot) = rownames(Exp)
  377. #2.加载样本信息
  378. info <- a
  379. colnames(info)=c("Sample","Type")
  380. Exp_plot$sam=info$Type
  381. Exp_plot$sam <- factor(Exp_plot$sam, levels = c("G2","G3","G4"))
  382. expr_use <- na.omit(Exp_plot)
  383. expr_use <-as.data.frame(expr_use)
  384. expr_use_long <- gather(expr_use, gene, Expression, -sam)
  385. table(expr_use_long$gene)
  387. colnames(expr_use_long) <- c("Group","gene","Expression")
  389. #设置比较组
  390. my_comparisons <- list(c("G2", "G3"), c("G3", "G4"), c("G2", "G4"))
  392. Custom.color <- c("#FF0033","#009934","#0088FF")
  394. ggplot(expr_use_long, aes(x = Group, y = Expression, fill=Group)) +
  396. geom_boxplot(position = position_nudge(x = 0.14),width=0.1,outlier.size = 0,outlier.alpha =0)+
  397. stat_halfeye(mapping = aes(fill=Group),width = 0.2, .width = 0, justification = -1.2, point_colour = NA,alpha=0.6) +
  398. scale_fill_manual(values = Custom.color)+
  399. scale_color_manual(values = Custom.color)+
  400. xlab(" ") +
  401. ylab("Expression") +
  402. ggtitle("TLN1")+
  403. theme_classic() +
  404. theme(
  405. legend.position = "none",
  406. axis.title.x = element_text(size = 13),
  407. axis.title.y = element_text(size = 13),
  408. axis.text.x = element_text(size = 12,hjust = 0.3),
  409. axis.text.y = element_text(size = 12),
  410. plot.title = element_text(hjust = 0.5)
  411. )+
  412. geom_signif(comparisons = my_comparisons,step_increase = .1,map_signif_level = TRUE,vjust = 0.5,hjust= 0)
  418. #2.1 DRG的多队列单因素cox回归分析####
  419. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  420. library(forestplot)
  421. library(grid)
  422. library(magrittr)
  423. library(checkmate)
  424. library(data.table)
  425. library(survival)
  426. library(survminer)
  428. gene=c("ACTB","ACTN4","CAPZB","CD2AP","DSTN","FLNA","FLNB","INF2",
  429. "IQGAP1","MYH10","MYH9","MYL6","PDLIM1","SLC7A11","TLN1")
  431. #TCGA(8个胶质瘤列队)
  432. load(file="TCGA.RData")
  433. dat=TCGA[gene,]
  434. dat=as.data.frame(t(dat))
  436. identical(rownames(dat),rownames(OS))
  437. data=cbind(OS,dat)
  439. univar_out = data.frame(matrix(NA,15,5))
  440. rownames(univar_out) = colnames(data)[-(1:2)]
  441. colnames(univar_out) = c("Coeffcient","HR","lower .95","upper .95","P-value")
  442. cox_data = data[,c(3:ncol(data))]
  443. cox_data = cbind(data[,c(1:2)],cox_data)
  444. str(cox_data)
  445. for(i in colnames(cox_data)[-(1:2)]){
  446. cox = coxph(Surv(OS.time, OS) ~ cox_data[,i], data = cox_data)
  447. cox_summ = summary(cox)
  448. univar_out[i,1] = cox_summ$coefficients[,1]
  449. univar_out[i,2] = cox_summ$coefficients[,2]
  450. univar_out[i,3] = cox_summ$conf.int[,3]
  451. univar_out[i,4] = cox_summ$conf.int[,4]
  452. univar_out[i,5] = cox_summ$coefficients[,5]
  453. }
  454. univar_out_0.05 = univar_out[univar_out[,5] < 0.05,]
  455. univar_out_TCGA=univar_out[,c(2,5)]
  458. #汇总
  459. sum=cbind(univar_out_TCGA,univar_out_CGGA325,univar_out_CGGA693,univar_out_GSE16011,
  460. univar_out_GSE108474,univar_out_EMTAB3892,univar_out_GSE4271,univar_out_GSE4412)
  461. colnames(sum)=c("TCGA_HR","TCGA_pvalue","CGGA325_HR","CGGA325_pvalue",
  462. "CGGA693_HR","CGGA693_pvalue","GSE16011_HR","GSE16011_pvalue",
  463. "GSE108474_HR","GSE108474_pvalue","EMTAB3892_HR","EMTAB3892_pvalue",
  464. "GSE4271_HR","GSE4271_pvalue","GSE4412_HR","GSE4412_pvalue")
  465. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.1 cox")
  466. write.table(sum,file="cox_result.txt",quote=F,sep="\t")
  469. #可视化
  470. library(pheatmap)
  471. library(RColorBrewer)
  472. results=read.table(file="heatmap.txt")
  473. datacol<-read.table("anno_col.txt",sep = "\t",header=T,row.names=1)
  474. rownames(datacol)
  475. datarowcolor = list(cohort=c(TCGA = "#F1B6DA", CGGA325 = "#8DD3C7",CGGA693 = "#BC80BD", GSE16011 = "#80B1D3",
  476. GSE108474 = "#FB8072",EMTAB3892 = "#BEBADA",GSE4271 = "#FDB462",GSE4412 = "#B3DE69"))
  477. pheatmap(results,color = colorRampPalette(c("#719dc9", "grey90", "#b595bf"))(25),
  478. border_color="grey30",
  479. cluster_rows = F,cluster_cols = F,shown_colnames=F,
  480. annotation_col = datacol,annotation_colors = datarowcolor,
  481. gaps_col = c(1,2,3,4,5,6,7))
  483. #在超过5个队列中cox.test pvalue < 0.05的基因有9个
  484. #ACTN4、CAPZB、CD2AP、FLNA、INF2、IQGAP1、MYH9、MYL6、PDLIM1
  488. #2.2 8个队列数据合并和去批-meta####
  489. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  490. load(file="TCGA.RData")
  491. load(file="CGGA325.RData")
  492. load(file="CGGA693.RData")
  493. load(file="GSE16011.RData")
  494. GSE16011=as.data.frame(t(GSE16011))
  495. load(file="GSE108474.RData")
  496. GSE108474=as.data.frame(t(GSE108474))
  497. load(file="E-MATE-3892.RData")
  498. E_MATE_3892=as.data.frame(t(E_MATE_3892))
  499. load(file="GSE4271.RData")
  500. GSE4271=as.data.frame(t(GSE4271))
  501. load(file="GSE4412.RData")
  502. GSE4412=as.data.frame(t(GSE4412))
  504. #取symbol交集(10495)
  505. intersects <- function (...) {
  506. Reduce(intersect, list(...))
  507. }
  508. a=intersects(rownames(TCGA),rownames(CGGA325),rownames(CGGA693),rownames(GSE16011),
  509. rownames(GSE108474),rownames(E_MATE_3892),rownames(GSE4271),rownames(GSE4412))
  511. TCGA=TCGA[a,]
  512. CGGA325=CGGA325[a,]
  513. CGGA693=CGGA693[a,]
  514. GSE16011=GSE16011[a,]
  515. GSE108474=GSE108474[a,]
  516. E_MATE_3892=E_MATE_3892[a,]
  517. GSE4271=GSE4271[a,]
  518. GSE4412=GSE4412[a,]
  520. identical(rownames(GSE16011),rownames(E_MATE_3892))
  522. meta=cbind(TCGA,CGGA325,CGGA693,GSE16011,
  523. GSE108474,E_MATE_3892,GSE4271,GSE4412) #2522
  525. meta_pheno=rbind(pheno,CGGA325_pheno,CGGA693_pheno,GSE16011_pheno,
  526. GSE108474_pheno,E_MATE_3892_pheno,GSE4271_pheno,GSE4412_pheno)
  528. group=c(rep("TCGA",691),rep("CGGA325",313),rep("CGGA693",657),rep("GSE16011",264),
  529. rep("GSE108474",284),rep("E_MATE_3892",151),rep("GSE4271",77),rep("GSE4412",85))
  531. save(meta,meta_pheno,group,file="meta.RData")
  535. #PCA(去批前)
  536. load(file="meta.RData")
  537. meta=cbind(meta,group)
  539. pca1 <- prcomp(meta[,-ncol(meta)],center = TRUE,scale. = TRUE)
  541. #提取PC score
  542. df1 <- pca1$x
  543. df1 <- as.data.frame(df1)
  545. #提取主成分的方差贡献率,生成坐标轴标题
  546. summ1 <- summary(pca1)
  547. xlab1 <- paste0("PC1(",round(summ1$importance[2,1]*100,2),"%)")
  548. ylab1 <- paste0("PC2(",round(summ1$importance[2,2]*100,2),"%)")
  550. library(ggplot2)
  551. ggplot(data = df1,aes(x = PC1,y = PC2,color = meta$group))+
  552. geom_point(size = 3)+
  553. labs(x = xlab1,y = ylab1,color = "Cohort",title = "Remove Batch Before")+
  554. guides(fill = "none")+
  555. theme_bw()+
  556. scale_colour_manual(values = c("#8DD3C7","#BC80BD","#BEBADA","#FB8072",
  557. "#80B1D3","#FDB462","#B3DE69","#F1B6DA"))+
  558. theme(plot.title = element_text(hjust = 0.5,size = 15),
  559. axis.text = element_text(size = 11),axis.title = element_text(size = 13),
  560. legend.text = element_text(size = 11),legend.title = element_text(size = 13),
  561. plot.margin = unit(c(0.4,0.4,0.4,0.4),'cm'))
  566. #去除批次效应
  567. library(sva)
  568. load(file="meta.RData")
  569. meta=as.data.frame(t(meta))
  570. Batch=group
  571. combat <- ComBat(dat = meta, batch = Batch)
  572. meta=as.data.frame(combat)
  574. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.2 meta")
  575. identical(rownames(meta_pheno),colnames(meta))
  576. save(meta,meta_pheno,group,file="meta.RData")
  578. table(meta_pheno$Grade)
  579. #G2 G3 G4
  580. #599 816 947
  582. gliomas=meta_pheno %>% filter(!is.na(Grade)) #2362
  583. lgg=gliomas[gliomas$Grade=="G2"|gliomas$Grade=="G3",] #1451
  584. gbm=gliomas[gliomas$Grade=="G4",] #947
  586. LGG=meta[,rownames(lgg)]
  587. GBM=meta[,rownames(gbm)]
  589. save(LGG,lgg,file="LGG.RData")
  590. save(GBM,gbm,file="GBM.RData")
  594. #PCA(去批后)
  595. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.2 meta")
  596. load(file="meta.RData")
  597. meta=as.data.frame(t(meta))
  598. meta=cbind(meta,group)
  600. pca1 <- prcomp(meta[,-ncol(meta)],center = TRUE,scale. = TRUE)
  602. #提取PC score
  603. df1 <- pca1$x
  604. df1 <- as.data.frame(df1)
  606. #提取主成分的方差贡献率,生成坐标轴标题
  607. summ1 <- summary(pca1)
  608. xlab1 <- paste0("PC1(",round(summ1$importance[2,1]*100,2),"%)")
  609. ylab1 <- paste0("PC2(",round(summ1$importance[2,2]*100,2),"%)")
  611. library(ggplot2)
  612. ggplot(data = df1,aes(x = PC1,y = PC2,color = meta$group))+
  613. geom_point(size = 3)+
  614. labs(x = xlab1,y = ylab1,color = "Cohort",title = "Remove Batch After")+
  615. guides(fill = "none")+
  616. theme_bw()+
  617. scale_colour_manual(values = c("#8DD3C7","#BC80BD","#BEBADA","#FB8072",
  618. "#80B1D3","#FDB462","#B3DE69","#F1B6DA"))+
  619. theme(plot.title = element_text(hjust = 0.5,size = 15),
  620. axis.text = element_text(size = 11),axis.title = element_text(size = 13),
  621. legend.text = element_text(size = 11),legend.title = element_text(size = 13),
  622. plot.margin = unit(c(0.4,0.4,0.4,0.4),'cm'))
  626. #2.3 无监督聚类(PCA,K-M,pheatmap)####
  627. library(ConsensusClusterPlus)
  628. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.2 meta")
  629. load(file="meta.RData")
  631. gene=c("ACTN4","CAPZB","CD2AP","FLNA","INF2","IQGAP1","MYH9","MYL6","PDLIM1")
  632. mydata=meta[gene,]
  633. mydata=as.matrix(mydata)
  635. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas")
  636. result_km <- ConsensusClusterPlus(mydata,
  637. maxK = 6,
  638. reps = 1000,
  639. pItem = 0.8,
  640. pFeature = 1,
  641. clusterAlg = "km",
  642. distance="euclidean",
  643. innerLinkage="complete",
  644. finalLinkage = "complete",
  645. title="consensus",
  646. plot="png",
  647. seed = 1234,
  648. writeTable=TRUE)
  649. save(result_km,file="meta-km.RData")
  650. load("meta-km.RData")
  652. #选择K值
  653. maxK = 6
  654. Kvec = 2:maxK
  655. x1 = 0.1; x2 = 0.9
  656. PAC = rep(NA,length(Kvec))
  657. names(PAC) = paste("K=",Kvec,sep="")
  659. for(i in Kvec){
  660. M = result_km[[i]]$consensusMatrix
  661. Fn = ecdf(M[lower.tri(M)])
  662. PAC[i-1] = Fn(x2) - Fn(x1)}
  664. #The optimal K
  665. optK = Kvec[which.min(PAC)]
  666. optK #[2]
  668. #查看分簇
  669. clusterNum=2
  670. cluster=result_km[[clusterNum]][["consensusClass"]]
  671. cluster=as.data.frame(cbind(sample=colnames(meta),group=cluster))
  672. table(cluster$group) #1-1341,2-1118
  673. cluster$group=ifelse(cluster$group=="1","Cluster1","Cluster2")
  674. write.csv(cluster,file="cluster.csv",row.names=F)
  677. #PCA可视化
  678. dat=mydata
  679. dat=as.data.frame(t(dat))
  680. dat=cbind(dat,group=cluster$group)
  682. pca1 <- prcomp(dat[,-ncol(dat)],center = TRUE,scale. = TRUE)
  684. #提取PC score
  685. df1 <- pca1$x
  686. df1 <- as.data.frame(df1)
  688. #提取主成分的方差贡献率,生成坐标轴标题
  689. summ1 <- summary(pca1)
  690. xlab1 <- paste0("PC1(",round(summ1$importance[2,1]*100,2),"%)")
  691. ylab1 <- paste0("PC2(",round(summ1$importance[2,2]*100,2),"%)")
  693. library(ggplot2)
  694. ggplot(data = df1,aes(x = PC1,y = PC2,color = dat$group))+
  695. stat_ellipse(aes(fill = dat$group),type = "norm",geom = "polygon",alpha = 0,color = "grey")+
  696. geom_point(size = 3,alpha=0.8)+
  697. labs(x = xlab1,y = ylab1,color = "Cohort",title = "Meta-cohort")+
  698. guides(fill = "none")+
  699. theme_bw()+
  700. scale_colour_manual(values = c("#F1B6DA","#BC80BD"))+
  701. theme(plot.title = element_text(hjust = 0.5,size = 15),
  702. axis.text = element_text(size = 11),axis.title = element_text(size = 13),
  703. legend.text = element_text(size = 11),legend.title = element_text(size = 13),
  704. plot.margin = unit(c(0.4,0.4,0.4,0.4),'cm'))
  707. #生存分析
  708. library(survival)
  709. library(survminer)
  710. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.2 meta")
  711. load(file="meta.RData")
  713. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas")
  714. cluster=read.csv("cluster.csv",header=T)
  716. identical(rownames(meta_pheno),cluster$sample)
  717. rt=cbind(meta_pheno[,c(1:2)],Type=cluster$group)
  719. diff=survdiff(Surv(OS.time, OS) ~ Type, data=rt)
  720. pValue=1-pchisq(diff$chisq, df=1)
  721. if(pValue<0.001){
  722. pValue="p<0.001"
  723. }else{
  724. pValue=paste0("p=", sprintf("%.03f",pValue))
  725. }
  726. fit <- survfit(Surv(OS.time, OS) ~ Type, data = rt)
  727. surPlot=ggsurvplot(fit,
  728. data=rt,
  729. conf.int=F,
  730. pval=pValue,
  731. pval.size=6,
  732. xlab="Time(years)",
  733. #ylab="Overall survival",
  734. legend.title="Cluster",
  735. break.time.by = 5,
  736. palette=c("#F1B6DA","#BC80BD"))
  740. #临床相关性热图
  741. library(pheatmap)
  742. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.2 meta")
  743. load("meta.RData")
  744. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas")
  745. cluster=read.csv("cluster.csv",header=T)
  746. identical(colnames(meta),cluster$sample)
  748. gene=c("ACTB","ACTN4","CAPZB","CD2AP","DSTN","FLNA","FLNB","INF2",
  749. "IQGAP1","MYH10","MYH9","MYL6","PDLIM1","SLC7A11","TLN1")
  751. data=cbind(meta_pheno[,c(1,3,4,5)],Cohort=group,Cluster=cluster$group,t(meta[gene,]))
  752. data$OS=ifelse(data$OS=="1","Dead","Alive")
  753. data$Age=ifelse(data$Age>47,">47","≤47")
  754. data$Gender=ifelse(data$Gender=="Male","Male","Female")
  758. #表达显著性验证
  759. library(tidyr)
  760. library(ggpubr)
  761. Exp_plot <- as.data.frame(data[,7:21])
  762. info <- as.data.frame(cbind(sample=rownames(data),group=data$Cluster))
  763. colnames(info)=c("Sample","Type")
  764. Exp_plot$sam=info$Type
  765. Exp_plot$sam <- factor(Exp_plot$sam, levels = c("Cluster1","Cluster2"))
  766. expr_use <- na.omit(Exp_plot)
  767. expr_use <-as.data.frame(expr_use)
  768. expr_use_long <- gather(expr_use, gene, Expression, -sam)
  769. table(expr_use_long$gene)
  770. colnames(expr_use_long) <- c("Group","gene","Expression")
  772. #绘图
  773. p=ggboxplot(expr_use_long, x="gene", y="Expression", color = "black", fill="Group",
  774. ylab="Gene expression",
  775. xlab="",
  776. legend.title=NULL,
  777. palette = c("#FF0033","#009934","#0088FF"),
  778. width=0.6, add = "none")
  779. p1=p+stat_compare_means(aes(group=Group),
  780. method="wilcox.test",
  781. symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", " ")),
  782. label = "p.signif")
  783. p1
  784. #表达显著性验证
  787. data=na.omit(data) #2294
  788. data = data[order(data[,"Cluster"]),]
  789. str(data)
  790. table(data$Cluster)
  792. rt=as.data.frame(t(data[,c(7:21)]))
  793. Type=data[,c(1:6)]
  795. anno_col=list(OS=c("Alive"="#A6D854","Dead"="#54b345"),
  796. Age=c("≤47"="#ff7f00",">47"="#FFBE7A"),
  797. Gender=c("Female"="#82B0D2","Male"="#a6cee3"),
  798. Grade=c("G2"="#C7E9B4","G3"="#7FCDBB","G4"="#41B6C4"),
  799. Cohort=c("TCGA"="#8DD3C7","CGGA325"="#FFFFB3","CGGA693"="#BEBADA","GSE16011"="#FB8072",
  800. "GSE108474"="#80B1D3","GSE4412"="#FDB462","GSE4271"="#B3DE69","E_MATE_3892"="#FCCDE5"),
  801. Cluster=c("Cluster1"="#F1B6DA","Cluster2"="#BC80BD"))
  803. pheatmap(rt, annotation=Type,
  804. main="Meta-cohort",
  805. color = colorRampPalette(colors = c("#313695","#4575B4","#74ADD1","#ABD9E9","white",
  806. "#F1B6DA","#DE77AE","#C51B7D","#8E0152"))(50),
  807. cluster_cols =F,
  808. cluster_rows = F,
  809. scale="row",
  810. show_colnames=F,
  811. annotation_colors = anno_col,
  812. fontsize=7.5,
  813. fontsize_row=8,
  814. gaps_col = 1250)
  824. #2.4 分簇的差异分析####
  825. library(limma)
  826. library(edgeR)
  827. library(dplyr)
  829. #RNAseq部分
  830. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  831. load(file="meta.RData")
  832. meta=as.data.frame(t(meta))
  833. meta=meta[,c(1768:2459)] #TCGA
  835. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  836. cluster=read.csv("tcga_cluster.csv")
  838. sample=cluster$sample
  839. meta=meta[,sample]
  840. identical(colnames(meta),cluster$sample)
  842. Group = factor(cluster$group,levels = c("Cluster1","Cluster2"))
  843. table(Group)
  845. #构建 DGEList 对象
  846. design = model.matrix(~0+Group)
  847. colnames(design)=levels(Group)
  848. rownames(design)=colnames(meta)
  849. dge = DGEList(counts=meta)
  851. #计算标准化因子
  852. dge = calcNormFactors(dge)
  854. #将矩阵进行voom转化,即计数数据转换为log2-counts per million (logCPM),生成Elist对象
  855. v = voom(dge,design,normalize="quantile")
  857. #转化后数据用于构建linear model
  858. fit = lmFit(v,design)
  860. #定义比较分组的信息
  861. constrasts = paste(rev(levels(Group)),collapse = "-")
  862. cont.matrix = makeContrasts(contrasts=constrasts,levels = design)
  864. #比较每个基因计算出标准差，使用经验贝叶斯方法进行校正并生成其他统计量
  865. fit3 = contrasts.fit(fit,cont.matrix)
  866. fit3 = eBayes(fit3)
  868. #提取差异分析数据
  869. DEG_limma = topTable(fit3, coef=constrasts, n=Inf)
  870. DEG_limma = na.omit(DEG_limma)
  872. #生成显著上下调gene标签列
  873. DEG_limma$group <- case_when(DEG_limma$logFC > 0.5 & DEG_limma$P.Value < 0.05 ~ "Up",
  874. DEG_limma$logFC < -0.5 & DEG_limma$P.Value < 0.05 ~ "Down",
  875. abs(DEG_limma$logFC) <= 0.5 ~ "None",
  876. DEG_limma$P.Value >= 0.05 ~ "None")
  878. table(DEG_limma$group)
  879. #CGGA325, d,519,u,516
  880. #CGGA693, d,499,u,465
  881. #TCGA, d,527,u,574
  883. #提取差异表达基因
  884. gene=DEG_limma[DEG_limma$group=="Up"|DEG_limma$group=="Down",]
  885. gene=as.data.frame(rownames(gene))
  886. colnames(gene)="tcga"
  887. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.4 DEG")
  888. write.csv(gene,file="tcga.csv",row.names = F)
  892. #microarray部分
  893. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  894. load(file="meta.RData")
  895. meta=as.data.frame(t(meta))
  896. meta=meta[,c(1662:1925)] #GSE16011
  898. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  899. cluster=read.csv("emate3892_cluster.csv")
  901. sample=cluster$sample
  902. meta=meta[,sample]
  903. identical(colnames(meta),cluster$sample)
  905. Group = factor(cluster$group,levels = c("Cluster1","Cluster2"))
  906. table(Group)
  908. #构建design
  909. design = model.matrix(~0+Group)
  910. colnames(design)=levels(Group)
  911. rownames(design)=colnames(meta)
  913. #非线性最小二乘分析
  914. fit = lmFit(meta,design)
  916. #定义比较分组的信息
  917. constrasts = paste(rev(levels(Group)),collapse = "-")
  918. cont.matrix = makeContrasts(contrasts=constrasts,levels = design)
  920. #比较每个基因计算出标准差，使用经验贝叶斯方法进行校正并生成其他统计量
  921. fit3 = contrasts.fit(fit,cont.matrix)
  922. fit3 = eBayes(fit3)
  924. #提取差异分析数据
  925. DEG_limma = topTable(fit3, coef=constrasts, n=Inf)
  926. DEG_limma = na.omit(DEG_limma)
  928. #生成显著上下调gene标签列
  929. DEG_limma$group <- case_when(DEG_limma$logFC > 0.5 & DEG_limma$P.Value < 0.05 ~ "Up",
  930. DEG_limma$logFC < -0.5 & DEG_limma$P.Value < 0.05 ~ "Down",
  931. abs(DEG_limma$logFC) <= 0.5 ~ "None",
  932. DEG_limma$P.Value >= 0.05 ~ "None")
  934. table(DEG_limma$group)
  935. #E-MATE-3892, d,385,u,693
  936. #GSE108474, d,634,u,1102
  937. #GSE16011, d,501,u,1701
  938. #GSE4271, d,404,u,435
  939. #GSE4412, d,571,u,772
  941. #提取差异表达基因
  942. gene=DEG_limma[DEG_limma$group=="Up"|DEG_limma$group=="Down",]
  943. gene=as.data.frame(rownames(gene))
  944. colnames(gene)="emate3892"
  945. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.4 DEG")
  946. write.csv(gene,file="emate3892.csv",row.names = F)
  950. #差异基因的单因素cox分析
  951. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.4 DEG/co-gene")
  952. DEG=read.table("DEG_in_7_cohort.txt")
  953. a=DEG$V1
  955. #8个队列单因素cox回归分析
  956. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  957. library(forestplot)
  958. library(grid)
  959. library(magrittr)
  960. library(checkmate)
  961. library(data.table)
  962. library(survival)
  963. library(survminer)
  965. gene=a
  967. #TCGA
  968. load(file="TCGA.RData")
  969. dat=TCGA[gene,]
  970. dat=as.data.frame(t(dat))
  972. identical(rownames(dat),rownames(OS))
  973. data=cbind(OS,dat)
  975. univar_out = data.frame(matrix(NA,294,5))
  976. rownames(univar_out) = colnames(data)[-(1:2)]
  977. colnames(univar_out) = c("Coeffcient","HR","lower .95","upper .95","P-value")
  978. cox_data = data[,c(3:ncol(data))]
  979. cox_data = cbind(data[,c(1:2)],cox_data)
  980. str(cox_data)
  981. for(i in colnames(cox_data)[-(1:2)]){
  982. cox = coxph(Surv(OS.time, OS) ~ cox_data[,i], data = cox_data)
  983. cox_summ = summary(cox)
  984. univar_out[i,1] = cox_summ$coefficients[,1]
  985. univar_out[i,2] = cox_summ$coefficients[,2]
  986. univar_out[i,3] = cox_summ$conf.int[,3]
  987. univar_out[i,4] = cox_summ$conf.int[,4]
  988. univar_out[i,5] = cox_summ$coefficients[,5]
  989. }
  990. univar_out_0.05 = univar_out[univar_out[,5] < 0.05,]
  991. univar_out_TCGA=univar_out_0.05[,c(2,5)]
  994. #取交集
  995. intersects <- function (...) {
  996. Reduce(intersect, list(...))
  997. }
  998. b=intersects(rownames(univar_out_TCGA),rownames(univar_out_CGGA325),rownames(univar_out_CGGA693),rownames(univar_out_GSE16011),
  999. rownames(univar_out_GSE108474),rownames(univar_out_EMTAB3892),rownames(univar_out_GSE4271),rownames(univar_out_GSE4412))
  1001. univar_out_TCGA=univar_out_TCGA[b,]
  1002. univar_out_CGGA325=univar_out_CGGA325[b,]
  1003. univar_out_CGGA693=univar_out_CGGA693[b,]
  1004. univar_out_GSE16011=univar_out_GSE16011[b,]
  1005. univar_out_GSE108474=univar_out_GSE108474[b,]
  1006. univar_out_EMTAB3892=univar_out_EMTAB3892[b,]
  1007. univar_out_GSE4271=univar_out_GSE4271[b,]
  1008. univar_out_GSE4412=univar_out_GSE4412[b,]
  1010. #汇总
  1011. sum=cbind(univar_out_TCGA,univar_out_CGGA325,univar_out_CGGA693,univar_out_GSE16011,
  1012. univar_out_GSE108474,univar_out_EMTAB3892,univar_out_GSE4271,univar_out_GSE4412)
  1013. colnames(sum)=c("TCGA_HR","TCGA_pvalue","CGGA325_HR","CGGA325_pvalue",
  1014. "CGGA693_HR","CGGA693_pvalue","GSE16011_HR","GSE16011_pvalue",
  1015. "GSE108474_HR","GSE108474_pvalue","EMTAB3892_HR","EMTAB3892_pvalue",
  1016. "GSE4271_HR","GSE4271_pvalue","GSE4412_HR","GSE4412_pvalue")
  1017. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.4 DEG/co-gene")
  1018. write.table(sum,file="cox_result.txt",quote=F,sep="\t")
  1021. #可视化
  1022. library(pheatmap)
  1023. library(RColorBrewer)
  1024. results=read.table(file="heatmap.txt")
  1025. datacol<-read.table("anno_col.txt",sep = "\t",header=T,row.names=1)
  1026. rownames(datacol)
  1027. datarowcolor = list(cohort=c(TCGA = "#A6CEE3", CGGA325 = "#1F78B4",CGGA693 = "#B2DF8A", GSE16011 = "#33A02C",
  1028. GSE108474 = "#FB9A99",EMTAB3892 = "#E31A1C",GSE4271 = "#FDBF6F",GSE4412 = "#FF7F00"))
  1029. pheatmap(results,
  1030. color = colorRampPalette(c("#B3DE69", "grey90", "#FF82A8"))(25),
  1031. border_color="grey30",
  1032. fontsize = 6,
  1033. cluster_rows = F,cluster_cols = F,shown_colnames=F,
  1034. annotation_col = datacol,annotation_colors = datarowcolor,
  1035. gaps_col = c(1,2,3,4,5,6,7))
  1040. #2.5 分簇的通路富集分析####
  1041. library(limma)
  1042. library(edgeR)
  1043. library(dplyr)
  1044. library(GSVA)
  1045. library(clusterProfiler)
  1046. library(tidyverse)
  1047. library(msigdbr)
  1048. library(GSEABase)
  1049. library(pheatmap)
  1050. library(BiocParallel)
  1051. library(ggplot2)
  1052. library(ggridges)
  1053. library(RColorBrewer)
  1054. library(org.Hs.eg.db)
  1055. library(enrichplot)
  1056. library(ggraph)
  1057. library(ggrepel)
  1058. library(qusage)
  1059. library(tibble)
  1060. library(stringr)
  1062. #GSVA-hallmark（error)
  1063. #表达矩阵
  1064. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1065. load(file="TCGA.RData")
  1066. gsva_data=TCGA
  1068. #分组信息
  1069. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  1070. cluster=read.csv("tcga_cluster.csv")
  1071. rownames(cluster)=cluster$sample
  1072. cluster=cluster[colnames(TCGA),]
  1074. #gsva
  1075. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.5 function")
  1076. hallmark_gmt <- read.gmt("h.all.v2023.2.Hs.symbols.gmt")
  1077. gene_set<-split(hallmark_gmt$gene,hallmark_gmt$term)
  1079. gsva_result <- gsva(as.matrix(gsva_data), gene_set, method = "gsva",min.sz=1,
  1080. max.sz=Inf,kcdf="Gaussian",parallel.sz=11)
  1083. #通路的差异分析
  1084. gl
  1085. Cluster1="Cluster1"
  1086. Cluster2="Cluster2"
  1088. group_list = cluster$group
  1089. design <- model.matrix(~0+factor(group_list))
  1090. colnames(design) <- levels(factor(group_list))
  1091. rownames(design) <- colnames(gsva_result)
  1092. contrast.matrix <- makeContrasts(contrasts=paste0(Cluster2,'-',Cluster1),
  1093. levels = design)
  1095. fit1 <- lmFit(gsva_result,design)
  1096. fit2 <- contrasts.fit(fit1, contrast.matrix)
  1097. efit <- eBayes(fit2)
  1099. summary(decideTests(efit,lfc=0.5, p.value=0.05))
  1100. tempOutput <- topTable(efit, coef=paste0(Cluster2,'-',Cluster1), n=Inf)
  1101. degs <- na.omit(tempOutput)
  1103. #火山图
  1104. degs1 <- degs %>%
  1105. rownames_to_column("pathway") %>%
  1106. mutate(Type = if_else(P.Value > 0.05,"ns",if_else(logFC >= 0, "up", "down"))) %>%
  1107. arrange(desc(abs(logFC)))
  1108. deg_path <- degs1
  1110. deg_path$pathway <- str_split_fixed(deg_path$pathway,"_",n=2)[,2]
  1111. deg_path$pathway <- gsub(pattern = "_"," ",deg_path$pathway)
  1114. dat_draw <- deg_path[,c(1,2)]
  1115. dat_draw$group <- ifelse(dat_draw$logFC>0 ,1,-1)
  1116. dat_draw$log <- -log10(deg_path$P.Value)*dat_draw$group
  1117. dat_draw$Type <- deg_path$Type
  1119. p <-ggplot(data = deg_path,
  1120. aes(x = logFC,
  1121. y = -log10(P.Value)
  1122. ))+
  1123. geom_point(alpha=0.7, size=4,
  1124. aes(color=Type))+
  1125. ylab("-log10(adj.p.Val)")+
  1126. scale_color_manual(values=c("#56a902", "grey","#d62a9d"))+
  1127. geom_vline(xintercept=0,lty=2,col="grey50",lwd=0.8)+
  1128. geom_hline(yintercept = -log10(0.05),lty=2,col="grey50",lwd=0.8)+
  1129. theme_bw()
  1130. p
  1132. #添加通路名称
  1133. deg_path$symbol=deg_path$pathway
  1134. deg_path=read.csv("deg_path.csv")
  1136. p1 <- p+geom_text_repel(data = deg_path, aes(x = logFC ,
  1137. y = -log10(P.Value),
  1138. label = label),
  1139. size = 3,box.padding = unit(0.5, "lines"),
  1140. point.padding = unit(0.8, "lines"),
  1141. segment.color = "black",
  1142. show.legend = FALSE)
  1143. p1
  1147. #ssGSEA-10条经典癌症通路
  1148. #表达矩阵
  1149. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1150. load(file="TCGA.RData")
  1151. gsva_data=TCGA
  1153. #分组信息
  1154. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  1155. cluster=read.csv("tcga_cluster.csv")
  1156. rownames(cluster)=cluster$sample
  1157. cluster=cluster[colnames(TCGA),]
  1159. #基因集信息
  1160. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.5 function")
  1161. cancer_pathway=read.csv("10_cancer_pathway_PMID30837276.csv")
  1162. cancer_pathway <- separate(cancer_pathway, Pathway, into = c("Pathway", "Type"), sep = "_")
  1164. activated=cancer_pathway[cancer_pathway$Type=="activated",]
  1165. activated <- split(activated$Symbol, activated$Pathway)
  1167. repressed=cancer_pathway[cancer_pathway$Type=="repressed",]
  1168. repressed <- split(repressed$Symbol, repressed$Pathway)
  1170. #ssgsea
  1171. gsva_result_a <- gsva(as.matrix(gsva_data), activated, method = "ssgsea",min.sz=1,
  1172. max.sz=Inf,kcdf="Gaussian",parallel.sz=11)
  1174. gsva_result_r <- gsva(as.matrix(gsva_data), repressed, method = "ssgsea",min.sz=1,
  1175. max.sz=Inf,kcdf="Gaussian",parallel.sz=11)
  1177. #结果处理
  1178. identical(colnames(gsva_result_a),colnames(gsva_result_r))
  1180. original_row_names <- rownames(gsva_result_a)
  1181. new_row_names <- paste(original_row_names, "_a", sep = "")
  1182. rownames(gsva_result_a) <- new_row_names
  1184. data=rbind(gsva_result_a,gsva_result_r)
  1185. write.csv(data,file="dat.csv")
  1187. #富集分数=激活分数-抑制分数
  1188. result=read.csv(file="dat_result.csv",row.names = 1)
  1189. identical(rownames(result),cluster$sample)
  1190. result$group=cluster$group
  1191. result_long <- gather(result, pathway, scores, -group)
  1193. #山峦图
  1194. ggplot(result_long, aes(x = scores, y = pathway, color = group, point_color = group, fill = group)) +
  1195. geom_density_ridges(jittered_points = TRUE, scale = 0.95, rel_min_height = 0.01,
  1196. point_shape = "|", point_size = 3, size = 0.25, position = position_points_jitter(height = 0)) +
  1197. scale_y_discrete(expand = c(0, 0)) + scale_x_continuous(expand = c(0, 0), name = " ") +
  1198. scale_fill_manual(values = c("#D55E0050", "#0072B250"), labels = c("Cluster1",
  1199. "Cluster2")) +
  1200. scale_color_manual(values = c("#D55E00", "#0072B2"), guide = "none") +
  1201. scale_discrete_manual("point_color", values = c("#D55E00", "#0072B2"), guide = "none") +
  1202. coord_cartesian(clip = "off") + guides(fill = guide_legend(override.aes = list(fill = c("#D55E00A0",
  1203. "#0072B2A0"), color = NA, point_color = NA))) + ggtitle(" ") +
  1204. theme_ridges(center = TRUE)+
  1205. theme(legend.position = "top")
  1209. #GSEA
  1210. #表达矩阵
  1211. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1212. load(file="TCGA.RData")
  1214. #分组信息
  1215. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  1216. cluster=read.csv("tcga_cluster.csv")
  1217. rownames(cluster)=cluster$sample
  1218. cluster=cluster[colnames(TCGA),]
  1220. Group = factor(cluster$group,levels = c("Cluster1","Cluster2"))
  1221. table(Group)
  1223. design = model.matrix(~0+Group)
  1224. colnames(design)=levels(Group)
  1225. rownames(design)=colnames(TCGA)
  1226. dge = DGEList(counts=TCGA)
  1228. dge = calcNormFactors(dge)
  1230. v = voom(dge,design,normalize="quantile")
  1232. fit = lmFit(v,design)
  1234. constrasts = paste(rev(levels(Group)),collapse = "-")
  1235. cont.matrix = makeContrasts(contrasts=constrasts,levels = design)
  1237. fit3 = contrasts.fit(fit,cont.matrix)
  1238. fit3 = eBayes(fit3)
  1240. DEG_limma = topTable(fit3, coef=constrasts, n=Inf)
  1241. DEG_limma = na.omit(DEG_limma)
  1243. DEG_limma$group <- case_when(DEG_limma$logFC > 0.5 & DEG_limma$adj.P.Val < 0.01 ~ "Up",
  1244. DEG_limma$logFC < -0.5 & DEG_limma$adj.P.Val < 0.01 ~ "Down",
  1245. abs(DEG_limma$logFC) <= 0.5 ~ "None",
  1246. DEG_limma$adj.P.Val >= 0.01 ~ "None")
  1248. table(DEG_limma$group)
  1249. #Down None Up
  1250. #1828 54433 2126
  1252. gene=DEG_limma[DEG_limma$group=="Up"|DEG_limma$group=="Down",]
  1253. DEGfilter=rownames(gene)
  1255. eg <- bitr(DEGfilter, fromType = "SYMBOL",
  1256. toType = c("ENTREZID", "ENSEMBL", "SYMBOL"),
  1257. OrgDb = "org.Hs.eg.db")
  1259. gene_info <- gene[DEGfilter, ] %>%
  1260. rownames_to_column(var = "SYMBOL") %>%
  1261. inner_join(., eg[, 1:2], by = "SYMBOL") %>%
  1262. arrange(desc(logFC))
  1264. gene_info <- unique(gene_info)
  1265. geneList <- gene_info$logFC
  1266. names(geneList) <- as.character(gene_info$SYMBOL)
  1267. head(geneList)
  1269. #msigdb C2
  1270. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.5 function")
  1271. kegg_gmt <- read.gmt("c2.cp.kegg_legacy.v2023.2.Hs.symbols.gmt")
  1272. kegg_gmt_df <- data.frame(
  1273. TERM = rep(names(kegg_gmt), sapply(kegg_gmt, length)),
  1274. GENE = unlist(kegg_gmt))
  1276. #gsea
  1277. gsea <- GSEA(geneList, TERM2GENE = kegg_gmt)
  1279. gsea_result=gsea@result
  1280. rownames(gsea_result)=NULL
  1282. #Cluster1
  1283. gseaplot2(gsea, geneSetID = c(1,2,6,9,11,14),pvalue_table = F,
  1284. color = c("#E495A5", "#86B875", "#7DB0DD","#FF7F00","#C51B7D","#FF0033"),
  1285. base_size = 11)
  1287. #Cluster2
  1288. gseaplot2(gsea, geneSetID = c(3,4,7,8,10,12,15),pvalue_table = F,
  1289. color = c("#FFCC33", "#86B875", "#7DB0DD","#FF7F00","#C51B7D","#377EB8",
  1290. "#FF0033"),base_size = 11)
  1294. #2.6 分簇的免疫浸润分析####
  1295. library(tidyverse)
  1296. library(IOBR)
  1297. library(reshape2)
  1298. library(ggplot2)
  1299. library(ggpubr)
  1300. library(ggsci)
  1301. library(pheatmap)
  1302. library(cowplot)
  1303. library(ComplexHeatmap)
  1304. library(circlize)
  1305. library(scales)
  1306. library(RColorBrewer)
  1307. library(ggh4x)
  1308. library(GSVA)
  1309. library(tinyarray)
  1311. #IOBR的6种算法免疫浸润
  1312. #与4.0immune区别是分组不同
  1313. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/IOBR")
  1314. load(file="tme_combine.RData")
  1316. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1317. load("TCGA.RData")
  1318. data = TCGA
  1320. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  1321. group=read.csv("tcga_cluster.csv")
  1323. sample=colnames(TCGA)
  1324. rownames(group)=group$sample
  1325. group=group[sample,]
  1326. identical(group$sample,colnames(TCGA))
  1328. #准备可视化
  1329. #合并数据
  1330. immune=cbind(cibersort[,c(1:26)],quantiseq[,2:12],mcp[,c(2:11)],xcell[,c(2:68)],
  1331. epic[,c(2:9)],ips[,c(2:7)])
  1332. rownames(immune)=immune$ID
  1333. immune=immune[,-1]
  1334. immune=t(immune)
  1336. #构建列分组文件
  1337. library(stringr)
  1338. Methods=c(rep('CIBERSORT',25),rep('Quantiseq',11),rep('MCP-counter',10),
  1339. rep('xCELL',67),rep('EPIC',8),rep('IPS',6))
  1340. type=data.frame(Methods=Methods)
  1341. rownames(type)=rownames(immune)
  1343. #标准化函数
  1344. standarize.fun <- function(indata=NULL, halfwidth=NULL, centerFlag=T, scaleFlag=T) {
  1345. outdata=t(scale(t(indata), center=centerFlag, scale=scaleFlag))
  1346. if (!is.null(halfwidth)) {
  1347. outdata[outdata>halfwidth]=halfwidth
  1348. outdata[outdata<(-halfwidth)]= -halfwidth
  1349. }
  1350. return(outdata)
  1351. }
  1353. #标准化免疫浸润数据，便于绘制热图
  1354. plotdata <- standarize.fun(immune,halfwidth = 2)
  1355. type$Methods=factor(type$Methods,levels = c('CIBERSORT','Quantiseq','MCP-counter',
  1356. 'xCELL','EPIC','IPS'))
  1358. #分组信息
  1359. my2=group[order(group$group),]
  1360. rownames(my2)=my2$sample
  1362. #排序
  1363. order2=rownames(my2)
  1364. table(my2$group)
  1365. Method=c('#66C2A5','#FC8D62','#8DA0CB','#E78AC3','#A6D854',"#FFD92F","#E5C494")
  1366. names(Method)=c('CIBERSORT','Quantiseq','MCP-counter','xCELL','EPIC','IPS')
  1368. #绘图
  1369. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.5 immu")
  1370. pdf('6IOBR_heatmap.pdf',height = 14,width = 8)
  1371. Heatmap(plotdata[,order2],name='Z-score',
  1372. top_annotation = HeatmapAnnotation(foo=anno_block(gp=gpar(fill= c('#F1B6DA','#BC80BD')),
  1373. labels=c('Cluster1','Cluster2'),
  1374. labels_gp = gpar(col='white',fontsize=10,fontface='bold'))),
  1375. cluster_rows = F,
  1376. col=colorRamp2(c(-2,0,2),c('#47b8e0','white','#FDAE61')),
  1377. color_space = "RGB",
  1378. cluster_columns = FALSE,border = T,
  1379. row_order=NULL,
  1380. row_names_side = 'left',
  1381. column_order=NULL,
  1382. show_column_names = FALSE,
  1383. row_names_gp = gpar(fontsize = 9),
  1384. column_split = c(rep(1,428),rep(2,263)),
  1385. row_split = type$Methods,
  1386. left_annotation = rowAnnotation(foo=anno_block(gp=gpar(fill= c('#F46D43','#FDAE61','#FEE090','#ABD9E9','#74ADD1',"#4575B4")),
  1387. labels=c('CIBERSORT','Quantiseq','MCP-counter','xCELL','EPIC','IPS'),
  1388. labels_gp = gpar(col='white',fontsize=8,fontface='bold'))),
  1389. gap = unit(1, "mm"),
  1390. column_title = NULL,
  1391. column_title_gp = gpar(fontsize = 10),
  1392. show_heatmap_legend =T,
  1393. heatmap_legend_param=list(labels_gp = gpar(fontsize = 10), border = T,
  1394. title_gp = gpar(fontsize = 10, fontface = "bold")),
  1395. column_gap = unit(2,'mm'),
  1396. row_title = NULL
  1397. )
  1398. dev.off()
  1402. #estimate
  1403. estimate=column_to_rownames(estimate,"ID")
  1404. estimate=estimate[,1:3]
  1405. a <- estimate
  1406. identical(rownames(a),rownames(group))
  1407. a$group=as.factor(group$group)
  1408. a <- a %>% rownames_to_column("sample")
  1409. b <- gather(a,key=category,value = score,-c(group,sample))
  1410. b1=as.data.frame(lapply(b$score,as.numeric)) %>% t() %>% as.data.frame()
  1411. b$score <- b1$V1
  1413. library(ggpubr)
  1414. library(reshape2)
  1415. p=ggviolin(b, x="category", y="score",
  1416. fill="group",#填充
  1417. palette =c("#F1B6DA","#6A3D9A"),
  1418. alpha=0.3,
  1419. add ="boxplot",
  1420. xlab = F,
  1421. legend = "right",
  1422. size = 0.5)
  1424. p+stat_compare_means(aes(group = group),
  1425. method = "wilcox.test",
  1426. label = "p.signif",
  1427. label.y = 4500,
  1428. symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
  1429. symbols = c("***", "**", "*", "ns")))+
  1430. theme(text = element_text(size=10),axis.text.x = element_text(angle=0, hjust=0.5))
  1435. #CIBERSORT免疫细胞
  1436. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/IOBR")
  1437. load(file="tme_combine.RData")
  1439. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1440. load("TCGA.RData")
  1441. data = TCGA
  1443. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  1444. group=read.csv("tcga_cluster.csv")
  1445. rownames(group)=group$sample
  1446. group=group[colnames(TCGA),]
  1448. cibersort=cibersort[,1:23]
  1449. cibersort=column_to_rownames(cibersort,"ID")
  1450. cibersort=t(cibersort)
  1452. rownames_original <- rownames(cibersort)
  1453. rownames_cleaned <- sub("_CIBERSORT", "", rownames_original)
  1454. rownames(cibersort) <- rownames_cleaned
  1456. draw_boxplot(cibersort,group$group,color = c("#F1B6DA","#BC80BD"),
  1457. xlab = " ",ylab = "Scores")
  1462. #免疫周期
  1463. #表达矩阵
  1464. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1465. load("TCGA.RData")
  1466. exp=as.matrix(TCGA)
  1468. #分组信息
  1469. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  1470. group=read.csv("tcga_cluster.csv")
  1471. rownames(group)=group$sample
  1472. group=group[colnames(TCGA),]
  1474. #immu cycle
  1475. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune")
  1476. geneset=read.csv("IOBR_TIME.csv")
  1477. geneset=geneset[1:7,]
  1478. split_genes <- strsplit(geneset$genes, ",")
  1479. new_geneset <- data.frame(
  1480. genesymbol = unlist(split_genes),
  1481. cell = rep(geneset$cell, lengths(split_genes)))
  1482. geneset = split(new_geneset$genesymbol,new_geneset$cell)
  1484. #ssgsea
  1485. re <- gsva(exp, geneset, method="ssgsea",mx.diff=FALSE, verbose=FALSE,
  1486. kcdf="Gaussian")
  1487. draw_boxplot(re,group$group,color = c("#F1B6DA","#BC80BD"),
  1488. xlab = " ",ylab = "Scores")
  1493. #免疫功能
  1494. #表达矩阵
  1495. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1496. load("TCGA.RData")
  1497. exp=as.matrix(TCGA)
  1499. #分组信息
  1500. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  1501. group=read.csv("tcga_cluster.csv")
  1502. rownames(group)=group$sample
  1503. group=group[colnames(TCGA),]
  1505. #immu function
  1506. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/TISIDB")
  1507. geneset = rio::import("ssGSEA-PMID30594216-Signature-13.xlsx",skip = 1)
  1508. geneset <- apply(geneset,2,function(x){return(as.character(na.omit(as.character(x))))})
  1510. #ssgsea
  1511. re <- gsva(exp, geneset, method="ssgsea",mx.diff=FALSE, verbose=FALSE,kcdf="Gaussian")
  1512. draw_boxplot(re,group$group,color = c("#F1B6DA","#BC80BD"),
  1513. xlab = " ",ylab = "Scores")
  1518. #TIME-Kobayashi/Bagaev
  1519. #表达矩阵
  1520. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1521. load("TCGA.RData")
  1522. exp=as.matrix(TCGA)
  1524. #分组信息
  1525. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  1526. group=read.csv("tcga_cluster.csv")
  1527. rownames(group)=group$sample
  1528. group=group[colnames(TCGA),]
  1530. #免疫细胞基因集
  1531. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/TISIDB")
  1532. dat=read.csv("TIME-Kobayashi.csv") #TIME-Kobayashi.csv/TIME-Bagaev.csv
  1533. geneset = split(dat$gene,dat$type)
  1535. #ssgsea
  1536. re <- gsva(exp, geneset, method="ssgsea",mx.diff=FALSE, verbose=FALSE,kcdf="Gaussian")
  1538. identical(colnames(TCGA),colnames(re))
  1539. re=rbind(re,group=group$group)
  1540. re=as.data.frame(t(re))
  1541. re=re[order(re$group),]
  1543. Cluster1 <- subset(re, group == "Cluster1")
  1544. Cluster2 <- subset(re, group == "Cluster2")
  1546. #Kobayashi
  1547. features <- c("Glycolysis", "IFNG response", "Inhibitory cells MDSCs", "Inhibitory cells Tregs", "Inhibitory molecules",
  1548. "Innate immunity", "Priming & activation", "Proliferation", "Recognition of tumor cells", "T cells")
  1549. #Bagaev
  1550. features <- c("Angiogenesis", "Anti tumor microenvironment", "Antigen presentation", "B cells", "CAF",
  1551. "Checkpoint inhibition", "Cytotoxic T and NK cells", "Granulocytes", "MDSC", "Treg",
  1552. "Tumor features","Tumor promotive immune infiltrate")
  1554. Cluster1_means <- sapply(Cluster1[, features], function(x) mean(as.numeric(x)))
  1555. Cluster2_means <- sapply(Cluster2[, features], function(x) mean(as.numeric(x)))
  1557. group_means <- data.frame(Cluster1 = Cluster1_means, Cluster2 = Cluster2_means)
  1558. rownames(group_means) <- features
  1560. #雷达图
  1561. library(fmsb)
  1562. library(ggradar)
  1563. library(cols4all)
  1564. library(ggplot2)
  1566. data=as.data.frame(t(group_means))
  1567. range(data)
  1568. data=rownames_to_column(data,'group')
  1569. mycol2=c("#F1B6DA","#6A3D9A")
  1571. ggradar(data,
  1572. grid.min = 0,grid.mid = 0.6,grid.max = 1.2,
  1573. gridline.mid.colour = 'grey',
  1574. values.radar = NA,
  1575. axis.label.size = 2.5,
  1576. legend.text.size = 8,
  1577. legend.position = 'bottom',
  1578. group.point.size = 1.5,group.line.width = 1,
  1579. fill = TRUE,fill.alpha = 0.3,
  1580. group.colours = mycol2,
  1581. background.circle.colour = 'white',background.circle.transparency = 0.2)
  1586. #2.7 分簇的突变分析####
  1587. #SNV
  1588. #SNV样本信息
  1589. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.1 genomic")
  1590. load("TCGA-LGG_maf.rdata")
  1591. barcode=as.data.frame(unique(data$Tumor_Sample_Barcode))
  1592. colnames(barcode)="sample"
  1593. data1=data
  1595. load("TCGA-GBM_maf.rdata")
  1596. barcode2=as.data.frame(unique(data$Tumor_Sample_Barcode))
  1597. colnames(barcode2)="sample"
  1598. data2=data
  1600. #合并lgg和gbm
  1601. data <- rbind(data1,data2)
  1602. maf.coad <- data
  1603. class(maf.coad) #data.frame
  1604. dim(maf.coad) #87957 141
  1606. #分组信息
  1607. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  1608. group=read.csv("tcga_cluster.csv")
  1610. #匹配结果
  1611. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.5 mutation/snv")
  1612. clin <- read.csv("SNV_barcode-group.csv",header=TRUE) #987
  1613. del <- which(clin$group=="#N/A")
  1614. clin <- clin[-del,] #732
  1615. colnames(clin)=c("Tumor_Sample_Barcode","sample")
  1617. maf.coad=merge(data,clin,by="Tumor_Sample_Barcode") #63987 142
  1618. sample=unique(maf.coad$Tumor_Sample_Barcode)
  1620. #maf文件
  1621. maf <- read.maf(maf.coad,clinicalData = clin)
  1623. #tmb计算
  1624. tmb_table_wt_log = tmb(maf = maf)
  1626. #合并分组信息
  1627. dat=as.data.frame(tmb_table_wt_log[,c(1,4)])
  1628. rownames(dat)=dat$Tumor_Sample_Barcode
  1629. rownames(clin)=clin$Tumor_Sample_Barcode
  1630. dat2=merge(dat,clin,by=0)
  1632. #更新分组名称
  1633. dat2$sample <- factor(dat2$sample, levels = c("Cluster1","Cluster2"))
  1634. p=ggviolin(dat2, x="sample", y="total_perMB_log", fill="sample",
  1635. palette =c("#1F78B4","#E31A1C"),alpha = 0.5,
  1636. add ="boxplot",size = 0.5,
  1637. xlab = "TMB",
  1638. legend = "none")
  1639. p+stat_compare_means(aes(group = sample),
  1640. method = "wilcox.test",label = "p.signif",
  1641. label.x = 1.5,label.y = 1.5,
  1642. symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
  1643. symbols = c("***", "**", "*", "ns")))
  1648. #CNV
  1649. library(data.table)
  1650. #segment数据准备
  1651. #cnv数据
  1652. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.5 mutation/cnv")
  1653. load("TCGA-LGG_CNV.Rdata")
  1654. lgg=data[2:7]
  1655. load("TCGA-GBM_CNV.Rdata")
  1656. gbm=data[2:7]
  1657. gliomas=rbind(lgg,gbm)
  1658. gliomas = gliomas[,c('Sample','Chromosome','Start','End','Num_Probes','Segment_Mean')]
  1659. write.csv(gliomas,file="gliomas.csv")
  1661. #分组信息
  1662. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.3 K-meas/cohort")
  1663. group=read.csv("tcga_cluster.csv")
  1664. write.csv(group,file="group.csv")
  1666. #匹配分组的数据
  1667. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.5 mutation/cnv")
  1668. gliomas=read.csv(file="gliomas-group.csv")
  1669. del <- which(gliomas$group=="#N/A")
  1670. gliomas <- gliomas[-del,]
  1671. colnames(gliomas)=c('Sample','Chromosome','Start','End','Num_Probes','Segment_Mean','group')
  1673. #cluster1
  1674. cluster1=gliomas[gliomas$group=="Cluster1",]
  1675. cluster1=cluster1[,1:6]
  1676. write.table(cluster1,file="Cluster1 segment file.txt",sep="\t",quote=F,row.names=F)
  1678. #cluster2
  1679. cluster2=gliomas[gliomas$group=="Cluster2",]
  1680. cluster2=cluster2[,1:6]
  1681. write.table(cluster2,file="Cluster2 segment file.txt",sep="\t",quote=F,row.names=F)
  1684. #markfile数据准备
  1685. marker_file = read.delim("snp6.na35.remap.hg38.subset.txt")
  1686. marker_file=marker_file[marker_file$freqcnv=="FALSE",]
  1687. marker_file = marker_file[,c(1,2,3)]
  1688. colnames(marker_file)=c("Marker Name", "Chromosome", "Marker Position")
  1689. write.table(marker_file,"marker file.txt",sep="\t",col.names=TRUE,row.name=FALSE)
  1691. #可视化
  1692. library(BSgenome.Hsapiens.UCSC.hg38)
  1693. library(ggplot2)
  1694. library(ggsci)
  1695. library(ggprism)
  1697. #染色体信息
  1698. df <- data.frame(chromName = seqnames(BSgenome.Hsapiens.UCSC.hg38),
  1699. chromlength = seqlengths(BSgenome.Hsapiens.UCSC.hg38)
  1700. )
  1701. df$chromNum <- 1:length(df$chromName)
  1702. df <- df[1:22,]
  1703. df$chromlengthCumsum <- cumsum(as.numeric(df$chromlength))
  1704. df$chormStartPosFrom0 <- c(0,df$chromlengthCumsum[-nrow(df)])
  1705. tmp_middle <- diff(c(0,df$chromlengthCumsum)) / 2
  1706. df$chromMidelePosFrom0 <- df$chormStartPosFrom0 + tmp_middle
  1708. #读取结果
  1709. scores <- read.table("/home/data/t060432/HRT/Disulfidptosis/data/2.5 mutation/cnv/578326-Cluster1/scores.gistic",
  1710. sep="\t",header=T,stringsAsFactors = F)
  1711. chromID <- scores$Chromosome
  1712. scores$StartPos <- scores$Start + df$chormStartPosFrom0[chromID]
  1713. scores$EndPos <- scores$End + df$chormStartPosFrom0[chromID]
  1714. range(scores$frequency)
  1715. scores[scores$Type == "Del", "frequency"] <- scores[scores$Type == "Del", "frequency"] * -1
  1716. range(scores$frequency)
  1718. df$ypos <- rep(c(0.15,0.18),11)
  1719. df$rect_col=ifelse(df$chromNum %% 2 == 0,"grey","white")
  1721. p1=ggplot(scores, aes(x = StartPos,y = frequency))+
  1722. geom_area(aes(group=Type, fill=factor(Type,levels = c("Del","Amp"))))+
  1723. scale_fill_manual(values = c("#1F78B4","#E31A1C"), guide=guide_legend(reverse = T, position = "top"), name="Type")+
  1724. geom_vline(data = df ,mapping=aes(xintercept=chromlengthCumsum),linetype=2)+
  1725. geom_text(data = df,aes(x=chromMidelePosFrom0,y=ypos,label=chromName))+
  1726. scale_x_continuous(expand = c(0,0),limits = c(0,2.9e9),name = NULL,labels = NULL)+
  1727. scale_y_continuous(expand = c(0.02,0.02),limits = c(-1,1),guide = "prism_offset")+
  1728. theme_prism()+
  1729. theme( axis.ticks.x = element_blank(),
  1730. axis.line.x = element_blank())
  1731. p1
  1733. #Frequency的比较
  1734. cluster1 <- read.table("/home/data/t060432/HRT/Disulfidptosis/data/2.5 mutation/cnv/578326-Cluster1/scores.gistic",
  1735. sep="\t",header=T,stringsAsFactors = F)
  1736. cluster1$group=c(rep("cluster1",17374))
  1738. cluster2 <- read.table("/home/data/t060432/HRT/Disulfidptosis/data/2.5 mutation/cnv/578337-Cluster2/scores.gistic",
  1739. sep="\t",header=T,stringsAsFactors = F)
  1740. cluster2$group=c(rep("cluster2",18047))
  1742. gliomas=rbind(cluster1,cluster2)
  1744. ggviolin(gliomas, x = 'group', y = 'Frequency', fill = 'group',palette = c("#1F78B4","#E31A1C"),
  1745. add = 'boxplot', add.params = list(fill = "white")) +
  1746. stat_compare_means(aes(group=group), label = "p.signif", bracket.size=0.5,
  1747. tip.length = 0.01, method = 'wilcox.test')
  1751. #3.0 模型的构建和验证 DisulfidpScore####
  1752. #见ML.R
  1754. #3.1 与其他模型的比较 Compared####
  1755. #见COM.R
  1757. #3.2 临床特征圈图####
  1758. library(cols4all)
  1759. library(ggplot2)
  1761. #TCGA
  1762. #表达矩阵
  1763. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1764. load("TCGA.RData")
  1766. #分组信息
  1767. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  1768. RS_mat <- read.table("RS_mat.txt",sep = "\t",header = T,stringsAsFactors = F,check.names = F,row.names=1)
  1769. bestmodel="Enet[alpha=0.3]"
  1770. RS_mat2=RS_mat[,bestmodel]
  1771. RS_mat2=RS_mat2[1:691]
  1772. group=ifelse(RS_mat2>median(RS_mat2),"High","Low")
  1774. #临床信息
  1775. data=cbind(pheno[,c(3,5,6,8,9,10)],group)
  1776. colnames(data)[1:5]=c("Grade","Age","Gender","IDH","1P19Q")
  1777. dt=na.omit(data)
  1778. dt$Age=as.numeric(dt$Age)
  1779. dt$Age=ifelse(dt$Age > 47,">47","≤47")
  1781. dt$grade=ifelse(dt$Grade=="G2","2",
  1782. ifelse(dt$Grade=="G2","3","4"))
  1783. dt$gender=ifelse(dt$Gender=="male","1","2")
  1784. dt$age=ifelse(dt$Age==">47","2","1")
  1785. dt$idh=ifelse(dt$IDH=="WT","1","2")
  1786. dt$`1p19q`=ifelse(dt$`1P19Q`=="non-codel","1","2")
  1787. dt$mgmtp=ifelse(dt$MGMTp=="Unmethylated","1","2")
  1789. str(dt)
  1790. dt[,8:13] <- apply(dt[,8:13],2,as.numeric)
  1792. #计算p值
  1793. #分组计算百分比
  1794. mydata=dt
  1795. mydata_summary <- mydata %>%
  1796. group_by(group, MGMTp) %>%
  1797. summarise(count = n()) %>%
  1798. group_by(group) %>%
  1799. mutate(percent = count / sum(count) * 100)
  1801. #卡方检验
  1802. chi_square <- chisq.test(table(mydata$group, mydata$MGMTp))
  1803. p_value <- chi_square$p.value
  1804. p_value
  1806. #圆圈图可视化
  1807. #分成高风险组和低风险组
  1808. high_group <- subset(dt, group == "High")
  1809. low_group <- subset(dt, group == "Low")
  1811. #计算百分比
  1812. high_group$fraction = high_group$mgmtp / sum(high_group$mgmtp)
  1813. low_group$fraction = low_group$`1p19q` / sum(low_group$`1p19q`)
  1815. #颜色
  1816. mycol <- c("#F1B6DA","#E78AC3","#C51B7D") #grade
  1817. mycol <- c("#f28147","#fac074") #age
  1818. mycol <- c("#9EBCDA","#377EB8") #gender
  1819. mycol <- c("#4DAF4A","#B8E186") #idh
  1820. mycol <- c("#D73027","#F1766D") #1p19q
  1821. mycol <- c("#88419D","#DECBE4") #mgmtp
  1823. #绘图
  1824. p5=ggplot(low_group, aes(x = 3,
  1825. y = fraction,
  1826. fill = `1P19Q` )) +
  1827. geom_col(width = 1.5) +
  1828. facet_grid(. ~ group) +
  1829. coord_polar(theta = "y") +
  1830. xlim(c(0.2, 3.8)) +
  1831. scale_fill_manual(values = mycol) +
  1832. theme_void() +
  1833. theme(
  1834. strip.text.x = element_text(size = 14),
  1835. legend.title = element_text(size = 15),
  1836. legend.text = element_text(size = 14))
  1841. #3.3 不同评分组的免疫分析#####
  1842. #IOBR(免疫浸润TCGA)
  1843. library(tidyverse)
  1844. library(IOBR)
  1845. library(reshape2)
  1846. library(ggplot2)
  1847. library(ggpubr)
  1848. library(ggsci)
  1849. library(pheatmap)
  1850. library(cowplot)
  1851. library(ComplexHeatmap)
  1852. library(circlize)
  1853. library(scales)
  1854. library(RColorBrewer)
  1855. library(ggh4x)
  1856. library(tidyr)
  1857. library(dplyr)
  1858. library(tibble)
  1860. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1861. load("TCGA.RData")
  1862. data = TCGA
  1864. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune")
  1865. #CIBERSORT
  1866. cibersort <- deconvo_tme(eset = data, method = "cibersort",
  1867. arrays = FALSE, perm = 1000)
  1869. #quantiseq
  1870. quantiseq <- deconvo_tme(eset = data,tumor = TRUE,arrays = FALSE,
  1871. scale_mrna = TRUE, method = "quantiseq")
  1873. #epic
  1874. epic <- deconvo_tme(eset = data, method = "epic", tumor = T, arrays = FALSE)
  1876. #mcp
  1877. mcp <- deconvo_tme(eset = data, method = "mcpcounter")
  1879. #xcell
  1880. xcell <- deconvo_tme(eset = data, method = "xcell", arrays = FALSE)
  1882. #estimate
  1883. estimate <- deconvo_tme(eset = data, method = "estimate")
  1885. #ips
  1886. ips<-deconvo_tme(eset = data, method = "ips", plot= FALSE)
  1888. #timer
  1889. timer <- deconvo_tme(eset = data, method = "timer", group_list = rep("glioma", dim(data)[2]))
  1891. #合并所有分析结果
  1892. tme_combine <- cibersort %>%
  1893. inner_join(quantiseq, "ID") %>%
  1894. inner_join(mcp, "ID") %>%
  1895. inner_join(xcell, "ID") %>%
  1896. inner_join(epic, "ID") %>%
  1897. inner_join(estimate, "ID") %>%
  1898. inner_join(ips, "ID")
  1900. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/IOBR")
  1901. save(cibersort,quantiseq,epic,mcp,xcell,estimate,ips,
  1902. tme_combine,file = 'tme_combine.RData')
  1903. load(file="tme_combine.RData")
  1905. #提取cibersort前22列数据
  1906. cibersort_data <- as.data.frame(cibersort[,1:23])
  1908. #分组信息
  1909. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  1910. RS_mat <- read.table("RS_mat.txt",sep = "\t",header = T,stringsAsFactors = F,check.names = F,row.names=1)
  1911. bestmodel="Enet[alpha=0.3]"
  1912. RS_mat2=RS_mat[,bestmodel]
  1913. RS_mat2=RS_mat2[1:691]
  1914. group=ifelse(RS_mat2>median(RS_mat2),"High","Low")
  1915. group=as.data.frame(cbind(ID=colnames(TCGA),group=group))
  1917. cibersort <- left_join(cibersort_data,group,by="ID")
  1918. cibersort <- melt(cibersort,id.vars=c("ID","group"))
  1919. colnames(cibersort)=c("sample","status1","variable","value")
  1920. cibersort=cibersort[order(cibersort$status1),]
  1921. cibersort$variable <- gsub("_CIBERSORT", "", cibersort$variable)
  1923. cols = c(brewer.pal(12, "Paired"),brewer.pal(8, "Dark2"),brewer.pal(12, "Set3"))
  1924. cols <- IOBR::palettes(category = "random", palette = 4,show_col = F,show_message = F)
  1925. cols=c("#FD8D3C","#FEB24C","#FED976","#FFED6F","#FFFF99","#EDF8B1","#C7E9B4","#A1D99B",
  1926. "#74C476","#7BCCC4","#4EB3D3","#4292C6","#2B8CBE","#9ECAE1","#C6DBEF","#9EBCDA",
  1927. "#8C96C6","#7570B3","#BEBADA","#DECBE4","#FCC5C0","#F768A1","#E7298A","#DD3497")
  1930. plot_tme = function(meltdf){
  1931. p_stack_df = ggplot(meltdf,aes(x = sample, y = value, fill = variable)) +
  1932. geom_bar(stat = "identity")+ theme_test() +
  1933. scale_fill_manual(values = cols) +
  1934. theme(legend.position = 'bottom')+
  1935. theme(plot.title = element_text(size = rel(2),hjust = 0.5),
  1936. axis.text.x = element_blank())+
  1937. theme(axis.ticks =element_blank() )+
  1938. scale_y_continuous(expand = c(0,0))+
  1939. facet_nested(.~status1,drop=T,scale="free",space="free",switch="x",
  1940. strip =strip_nested(background_x = elem_list_rect(fill =c('#C51B7D','#377EB8')),by_layer_x = F))
  1941. # p_stack_df
  1943. p_box_df = ggplot(meltdf,aes(x = variable, y = value,
  1944. fill = variable)) +
  1945. geom_boxplot(width=0.4,lwd=0.2,color='black',outlier.shape=NA,
  1946. position = position_dodge(width = 0.8))+
  1947. theme_bw()+
  1948. theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = .5))+
  1949. theme(legend.position = 'none')+
  1950. scale_fill_manual(values = cols)
  1951. # p_box_df
  1953. p_box_comp = ggplot(meltdf,aes(x = variable, y = value,
  1954. fill = status1,color = status1)) +
  1956. geom_boxplot(width=0.5,lwd=0.2,color='black',outlier.shape=NA,
  1957. position = position_dodge(width = 0.8))+
  1958. theme_bw()+
  1959. theme(legend.position = 'top')+
  1960. theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = .5))+
  1961. stat_compare_means(aes(group=status1,label = after_stat(p.signif)),hide.ns = FALSE,
  1962. method = "t.test")+
  1963. scale_fill_manual(values = c('#C51B7D','#377EB8'))
  1964. # p_box_comp
  1966. p_heat = ggplot(meltdf,aes(x = sample, y =variable ,fill=value))+
  1967. geom_tile()+
  1968. scale_fill_gradientn(colors = c('white','#C51B7D','#377EB8'))+
  1970. theme(axis.text.x = element_blank(),
  1971. axis.ticks.x = element_blank(),
  1972. panel.border = element_rect(fill = 'NA',color = 'black'))+
  1973. labs(fill="Z score", x = NULL, y = NULL)
  1974. # p_heat
  1976. #plot_comb = (p_stack_df)/(p_box_df|p_box_comp)/p_heat
  1977. plot_comb = (p_box_comp)
  1978. return(plot_comb)
  1979. }
  1981. plot_exist_cib = plot_tme(cibersort)
  1982. print(plot_exist_cib)
  1986. #TME signature(TCGA)
  1987. #表达矩阵
  1988. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  1989. load("TCGA.RData")
  1990. exp=as.matrix(TCGA)
  1992. #分组信息
  1993. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  1994. RS_mat <- read.table("RS_mat.txt",sep = "\t",header = T,stringsAsFactors = F,check.names = F,row.names=1)
  1995. bestmodel="Enet[alpha=0.3]"
  1996. RS_mat2=RS_mat[,bestmodel]
  1997. RS_mat2=RS_mat2[1:691]
  1998. group=ifelse(RS_mat2>median(RS_mat2),"High","Low")
  1999. group=as.data.frame(cbind(sample=colnames(TCGA),group=group))
  2000. group=group[order(group$group),]
  2001. rownames(group)=group$sample
  2003. #IOBR中的TME特征
  2004. library(IOBR)
  2005. library(tidyr)
  2006. data=signature_collection
  2007. data2=signature_tme
  2009. #提取需要的特征
  2010. expr_matrix <- matrix(nrow = 119, ncol = 2)
  2011. for (i in 1:length(data2)) {
  2012. expr_matrix[i, 1] <- names(data2)[i]
  2013. expr_matrix[i, 2] <- paste(data2[[i]], collapse = ",") ）
  2014. }
  2016. expr_df <- as.data.frame(expr_matrix)
  2017. colnames(expr_df) <- c("cell", "genes")
  2018. needata=expr_df[c(1,2,3,4,5,6,8,9,10,13,14,15,16,17,18,19,20,
  2019. 29,30,31,32,33,34,35,36,37,38,39,40,
  2020. 41,42,43,55,56,57,58,59,60,
  2021. 61,62,63,64,65,66,
  2022. 94,97,98,99,100,101,102,103,104,
  2023. 111,112,113,114,115,116,117),]
  2024. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/IOBR_TIME")
  2025. write.csv(needata,file="TIME.csv")
  2027. #获得geneset
  2028. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/IOBR_TIME")
  2029. geneset=read.csv("IOBR_TIME.csv")
  2030. split_genes <- strsplit(geneset$genes, ",")
  2031. new_geneset <- data.frame(
  2032. genesymbol = unlist(split_genes),
  2033. cell = rep(geneset$cell, lengths(split_genes)))
  2034. geneset = split(new_geneset$genesymbol,new_geneset$cell)
  2036. #ssgsea
  2037. re <- gsva(exp, geneset, method="ssgsea",mx.diff=FALSE, verbose=FALSE)
  2038. re=as.data.frame(re)
  2039. type=unique(new_geneset$cell)
  2040. re=re[type,]
  2041. write.csv(re,file="TIMEresult.csv")
  2042. identical(colnames(re),colnames(TCGA))
  2047. #IOBR免疫浸润(meta-LGG-GBM)
  2048. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.2 meta")
  2049. load("GBM.RData")
  2050. data = GBM
  2052. #CIBERSORT
  2053. cibersort <- deconvo_tme(eset = data, method = "cibersort",
  2054. arrays = FALSE, perm = 1000)
  2056. #quantiseq
  2057. quantiseq <- deconvo_tme(eset = data,tumor = TRUE,arrays = FALSE,
  2058. scale_mrna = TRUE, method = "quantiseq")
  2060. #epic
  2061. epic <- deconvo_tme(eset = data, method = "epic", tumor = T, arrays = FALSE)
  2063. #mcp
  2064. mcp <- deconvo_tme(eset = data, method = "mcpcounter")
  2066. #xcell
  2067. xcell <- deconvo_tme(eset = data, method = "xcell", arrays = FALSE)
  2069. #estimate
  2070. estimate <- deconvo_tme(eset = data, method = "estimate")
  2072. #ips
  2073. ips<-deconvo_tme(eset = data, method = "ips", plot= FALSE)
  2075. #合并所有分析结果
  2076. tme_combine <- cibersort %>%
  2077. inner_join(quantiseq, "ID") %>%
  2078. inner_join(mcp, "ID") %>%
  2079. inner_join(xcell, "ID") %>%
  2080. inner_join(epic, "ID") %>%
  2081. inner_join(estimate, "ID") %>%
  2082. inner_join(ips, "ID")
  2084. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/IOBR")
  2085. save(cibersort,quantiseq,epic,mcp,xcell,estimate,ips,
  2086. tme_combine,file = 'tme_combine_GBM.RData')
  2090. #TME signature(meta-LGG-GBM)
  2091. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.2 meta")
  2092. load("GBM.RData")
  2094. #分组信息
  2095. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  2096. sample=read.table("sample.txt")
  2097. RS_mat <- read.table("RS_mat.txt",sep = "\t",header = T,stringsAsFactors = F,check.names = F,row.names=1)
  2098. bestmodel="Enet[alpha=0.3]"
  2099. RS_mat2=RS_mat[,bestmodel]
  2100. group=RS_mat2
  2101. group_list=as.data.frame(cbind(sample=sample$sample,group=group))
  2102. rownames(group_list)=group_list$sample
  2104. sample=colnames(GBM)
  2105. group_list=group_list[sample,]
  2106. group_list$group=as.numeric(group_list$group)
  2107. group_list=na.omit(group_list)
  2108. group_list$group=ifelse(group_list$group>median(group_list$group),"High","Low")
  2109. sample=rownames(group_list)
  2110. GBM=GBM[,sample]
  2111. exp=as.matrix(GBM)
  2113. #IOBR中的TME特征
  2114. library(IOBR)
  2115. library(tidyr)
  2116. library(GSVA)
  2118. #获得geneset
  2119. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune")
  2120. geneset=read.csv("IOBR_TIME.csv")
  2121. split_genes <- strsplit(geneset$genes, ",")
  2122. new_geneset <- data.frame(
  2123. genesymbol = unlist(split_genes),
  2124. cell = rep(geneset$cell, lengths(split_genes)))
  2125. geneset = split(new_geneset$genesymbol,new_geneset$cell)
  2127. #ssgsea
  2128. re <- gsva(exp, geneset, method="ssgsea",mx.diff=FALSE, verbose=FALSE)
  2129. re=as.data.frame(re)
  2130. type=unique(new_geneset$cell)
  2131. re=re[type,]
  2132. write.csv(re,file="TIMEresult_meta.csv")
  2133. identical(colnames(re),colnames(meta))
  2137. #结果整合(meta-LGG-GBM)
  2138. #cibersort&estimate
  2139. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/IOBR")
  2140. load("tme_combine_GBM.RData")
  2141. cibersort=cibersort[,1:23]
  2142. cibersort=column_to_rownames(cibersort,"ID")
  2143. cibersort=cibersort[sample,]
  2144. estimate=estimate[,1:4]
  2145. estimate=column_to_rownames(estimate,"ID")
  2146. estimate=estimate[sample,]
  2148. #TIME
  2149. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune")
  2150. TIME=read.csv("TIMEresult_GBM.csv",check.names=F,row.names = 1)
  2151. TIME=as.data.frame(t(TIME))
  2153. identical(rownames(estimate),rownames(cibersort))
  2154. identical(rownames(TIME),rownames(cibersort))
  2156. #合并
  2157. data=as.data.frame(t(cbind(cibersort,estimate,TIME)))
  2158. data=as.data.frame(t(data))
  2159. data$group=group_list$group
  2160. library(dplyr)
  2162. #高低风险组组的平均值
  2163. grouped_data <- data %>%
  2164. group_by(group) %>%
  2165. summarize(across(everything(), mean, na.rm = TRUE))
  2166. grouped_data=t(grouped_data)
  2167. write.csv(grouped_data,"heatmap_GBM.csv")
  2172. #总体可视化
  2173. #cibersort&estimate
  2174. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/IOBR")
  2175. load("tme_combine.RData")
  2176. cibersort=cibersort[,1:23]
  2177. cibersort=column_to_rownames(cibersort,"ID")
  2178. estimate=estimate[,1:4]
  2179. estimate=column_to_rownames(estimate,"ID")
  2181. #TIME
  2182. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune")
  2183. TIME=read.csv("TIMEresult.csv",check.names=F,row.names = 1)
  2184. TIME=as.data.frame(t(TIME))
  2186. identical(rownames(cibersort),rownames(estimate))
  2187. identical(rownames(TIME),rownames(cibersort))
  2189. data=as.data.frame(t(cbind(cibersort,estimate,TIME)))
  2190. standarize.fun <- function(indata=NULL, halfwidth=NULL, centerFlag=T, scaleFlag=T) {
  2191. outdata=t(scale(t(indata), center=centerFlag, scale=scaleFlag))
  2192. if (!is.null(halfwidth)) {
  2193. outdata[outdata>halfwidth]=halfwidth
  2194. outdata[outdata<(-halfwidth)]= -halfwidth}
  2195. return(outdata)}
  2196. data <- standarize.fun(data,halfwidth = 2)
  2198. types=read.csv("type.csv",check.names = F)
  2199. rownames(types)=types$celltype
  2200. table(types$type)
  2201. types$type=factor(types$type,levels = c('cibersort','estimate','cancer immunity cycle',
  2202. 'immune functions','TME signatures'))
  2204. order2=rownames(group)
  2205. table(group$group)
  2206. type=c('#E64B35FF','#4DBBD5FF','#00A087FF','#3C5488FF','#F39B7FFF',"orange","yellow")
  2207. names(type)=c('cibersort','estimate','cancer immunity cycle','immune functions',
  2208. 'TME signatures')
  2210. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune")
  2211. pdf('TIME_heatmap.pdf',height = 8,width = 6)
  2212. Heatmap(data[,order2],name='Z-score',
  2213. top_annotation = HeatmapAnnotation(foo=anno_block(gp=gpar(fill= c('lightgreen','#FB9A99')),
  2214. labels=c('High','Low'),
  2215. labels_gp = gpar(col='white',fontsize=10,fontface='bold'))),
  2216. cluster_rows = F,
  2217. col=colorRamp2(c(-2,0,2),c('lightblue','white','pink')),
  2218. color_space = "RGB",
  2219. cluster_columns = FALSE,border = T,
  2220. row_order=NULL,
  2221. row_names_side = 'left',
  2222. column_order=NULL,
  2223. show_column_names = FALSE,
  2224. row_names_gp = gpar(fontsize = 9),
  2225. column_split = c(rep(1,345),rep(2,346)),
  2226. row_split = types$type,
  2227. left_annotation = rowAnnotation(foo=anno_block(gp=gpar(fill= c("#FF7F00",'#FDAE61','#FED976','#74ADD1',"#4575B4")),
  2228. labels=c('cibersort','estimate','cancer immunity cycle','immune functions',
  2229. 'TME signatures'),
  2230. labels_gp = gpar(col='white',fontsize=8,fontface='bold'))),
  2231. gap = unit(1, "mm"),
  2232. column_title = NULL,
  2233. column_title_gp = gpar(fontsize = 10),
  2234. show_heatmap_legend =T,
  2235. heatmap_legend_param=list(labels_gp = gpar(fontsize = 10), border = T,
  2236. title_gp = gpar(fontsize = 10, fontface = "bold")),
  2237. column_gap = unit(2,'mm'),
  2238. row_title = NULL)
  2239. dev.off()
  2244. #免疫细胞相关性
  2245. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/IOBR")
  2246. load("tme_combine.RData")
  2247. immu_data=column_to_rownames(cibersort,"ID")
  2248. immu_data=immu_data[,1:22]
  2250. #去除cibersort后缀
  2251. col_names <- colnames(immu_data)
  2252. new_col_names <- gsub("_CIBERSORT", "", col_names)
  2253. colnames(immu_data) <- new_col_names
  2255. #分组信息
  2256. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  2257. RS_mat <- read.table("RS_mat.txt",sep = "\t",header = T,stringsAsFactors = F,check.names = F,row.names=1)
  2258. bestmodel="Enet[alpha=0.3]"
  2259. RS_mat2=RS_mat[,bestmodel]
  2260. y=RS_mat2[1:691]
  2262. correlation <- data.frame()
  2263. for(i in 1:length(colnames(immu_data))){
  2264. print(i)
  2265. dd = cor.test(as.numeric(immu_data[,i]),y,method="spearman")
  2266. correlation[i,1] = colnames(immu_data)[i]
  2267. correlation[i,2] = dd$estimate
  2268. correlation[i,3] = dd$p.value}
  2269. colnames(correlation) <- c("Cell","cor","pvalue")
  2271. #定义圆圈颜色的函数
  2272. p.col = c('pink','gold','orange','LimeGreen','darkgreen')
  2273. fcolor = function(x,p.col){
  2274. color = ifelse(x>0.8,p.col[1],ifelse(x>0.6,p.col[2],ifelse(x>0.4,p.col[3],
  2275. ifelse(x>0.2,p.col[4], p.col[5])
  2276. )))
  2277. return(color)
  2278. }
  2280. #定义设置圆圈大小的函数
  2281. p.cex = seq(2.5, 5.5, length=5)
  2282. fcex = function(x){
  2283. x=abs(x)
  2284. cex = ifelse(x<0.1,p.cex[1],ifelse(x<0.2,p.cex[2],ifelse(x<0.3,p.cex[3],
  2285. ifelse(x<0.4,p.cex[4],p.cex[5]))))
  2286. return(cex)
  2287. }
  2289. data=correlation
  2290. #根据pvalue定义圆圈的颜色
  2291. points.color = fcolor(x=data$pvalue,p.col=p.col)
  2292. data$points.color = points.color
  2293. points.cex = fcex(x=data$cor)
  2294. data$points.cex = points.cex
  2295. data=data[order(data$cor),]
  2297. #可视化
  2298. xlim = ceiling(max(abs(data$cor))*10)/10
  2299. layout(mat=matrix(c(1,1,1,1,1,0,2,0,3,0),nc=2),width=c(8,2.2),heights=c(1,2,1,2,1))
  2300. par(bg="white",las=1,mar=c(5,18,2,4),cex.axis=1.5,cex.lab=2)
  2301. plot(1,type="n",xlim=c(-xlim,xlim),ylim=c(0.5,nrow(data)+0.5),xlab="Correlation Coefficient",ylab="",yaxt="n",yaxs="i",axes=F)
  2302. rect(par('usr')[1],par('usr')[3],par('usr')[2],par('usr')[4],col="#F5F5F5",border="#F5F5F5")
  2303. grid(ny=nrow(data),col="white",lty=1,lwd=2)
  2305. #绘制图形的线段
  2306. segments(x0=data$cor,y0=1:nrow(data),x1=0,y1=1:nrow(data),lwd=4)
  2307. #绘制图形的圆圈
  2308. points(x=data$cor,y = 1:nrow(data),col = data$points.color,pch=16,cex=data$points.cex)
  2309. #展示免疫细胞的名称
  2310. text(par('usr')[1],1:nrow(data),data$Cell,adj=1,xpd=T,cex=1.5)
  2312. #展示pvalue
  2313. pvalue.text=ifelse(data$pvalue<0.001,'<0.001',sprintf("%.03f",data$pvalue))
  2314. redcutoff_cor=0
  2315. redcutoff_pvalue=0.05
  2316. text(par('usr')[2],1:nrow(data),pvalue.text,adj=0,xpd=T,col=ifelse(abs(data$cor)>redcutoff_cor & data$pvalue<redcutoff_pvalue,"red","black"),cex=1.5)
  2317. axis(1,tick=F)
  2319. #绘制圆圈大小的图例
  2320. par(mar=c(0,4,3,4))
  2321. plot(1,type="n",axes=F,xlab="",ylab="")
  2322. legend("left",legend=c(0.1,0.2,0.3,0.4,0.5),col="black",pt.cex=p.cex,pch=16,bty="n",cex=2,title="abs(cor)")
  2324. #绘制圆圈颜色的图例
  2325. par(mar=c(0,6,4,6),cex.axis=1.5,cex.main=2)
  2326. barplot(rep(1,5),horiz=T,space=0,border=NA,col=p.col,xaxt="n",yaxt="n",xlab="",ylab="",main="pvalue")
  2327. axis(4,at=0:5,c(1,0.8,0.6,0.4,0.2,0),tick=F)
  2333. #免疫表型比例（TCGA-meta-LGG-GBM）
  2334. library(ImmuneSubtypeClassifier)
  2335. library(dplyr)
  2337. #表达矩阵
  2338. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  2339. load("TCGA.RData")
  2340. tpm=TCGA
  2342. #分组信息
  2343. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  2344. RS_mat <- read.table("RS_mat.txt",sep = "\t",header = T,stringsAsFactors = F,check.names = F,row.names=1)
  2345. bestmodel="Enet[alpha=0.3]"
  2346. RS_mat2=RS_mat[,bestmodel]
  2347. RS_mat2=RS_mat2[1:691]
  2348. group=ifelse(RS_mat2>median(RS_mat2),"High","Low")
  2349. group_list=as.data.frame(cbind(sample=colnames(TCGA),group=group))
  2350. rownames(group_list)=group_list$sample
  2352. Isubtype <- callEnsemble(X = tpm, geneids = 'symbol')[,1:2]
  2353. Isubtype=column_to_rownames(Isubtype,"SampleIDs")
  2354. Isubtype$BestCall=paste0('C',Isubtype$BestCall)
  2356. #查看为匹配的基因有多少
  2357. geneMatchErrorReport(X=tpm, geneid='symbol')
  2359. #可视化
  2360. df=cbind(group_list[rownames(Isubtype),],Isubtype)%>%
  2361. select(ncol(.)-1,ncol(.))%>%
  2362. group_by(group,BestCall)%>%
  2363. summarise(count=n())
  2364. df$group <- factor(df$group, levels = c("Low","High"))
  2365. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/immuphenotype")
  2366. save(df,file="df.RData")
  2367. load("df.RData")
  2369. library(ggplot2)
  2370. library(paletteer)
  2371. ggplot()+
  2372. geom_bar(data =df, aes(x = group, y = count, fill = BestCall),
  2373. stat = "identity",
  2374. position = "fill")+
  2375. theme_classic()+
  2376. scale_fill_paletteer_d("RColorBrewer::Paired")
  2383. #3.4 不同评分组的突变分析####
  2384. library(broom)
  2385. library(maftools)
  2386. library(ggpubr)
  2387. library(stringr)
  2389. #SNV样本信息
  2390. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.1 genomic")
  2391. load("TCGA-LGG_maf.rdata")
  2392. barcode=as.data.frame(unique(data$Tumor_Sample_Barcode))
  2393. colnames(barcode)="sample"
  2394. data1=data
  2396. load("TCGA-GBM_maf.rdata")
  2397. barcode2=as.data.frame(unique(data$Tumor_Sample_Barcode))
  2398. colnames(barcode2)="sample"
  2399. data2=data
  2401. #合并lgg和gbm
  2402. data <- rbind(data1,data2)
  2403. maf.coad <- data
  2404. class(maf.coad) #data.frame
  2405. dim(maf.coad) #87957 141
  2407. #匹配结果
  2408. clin <- read.csv("TCGA_SNV_barcode.csv",header=TRUE) #732
  2409. clin <- clin[-c(499,586),] #去除2个异常样本 TCGA-06-5416-01A-01D-1486-08,TCGA-DU-6392-01A-11D-1705-08
  2411. maf.coad=merge(data,clin,by="Tumor_Sample_Barcode") #36059 142
  2412. sample=unique(maf.coad$Tumor_Sample_Barcode)
  2414. #提取low risk组
  2415. low=maf.coad[maf.coad$sample=="Low",]
  2416. maf <- read.maf(low,clinicalData = clin)
  2418. #提取high risk组
  2419. high=maf.coad[maf.coad$sample=="High",]
  2420. maf <- read.maf(high,clinicalData = clin)
  2422. #瀑布图
  2423. vc_cols = c("#8DD3C7","#80B1D3","#B2DF8A","#33A02C","#FB9A99","#E31A1C","#FDBF6F","#FF7F00",'#CAB2D6')
  2424. names(vc_cols) = c('Multi_Hit','Missense_Mutation','Frame_Shift_Del','Nonsense_Mutation',
  2425. 'Frame_Shift_Ins','In_Frame_Ins','Splice_Site','In_Frame_Del','Translation_Start_Site')
  2426. col = c("#E31A1C","#1F78B4")
  2427. names(col) = c('High','Low')
  2428. oncoplot(maf = maf, top=30, draw_titv = F,fontSize = 0.75 ,
  2429. clinicalFeatures = 'sample',annotationColor = list(sample=col),
  2430. colors = vc_cols,bgCol = "transparent",annoBorderCol = "white",
  2431. sortByAnnotation = TRUE,borderCol=NULL)
  2435. #合并lgg和gbm
  2436. maf <- read.maf(maf.coad,clinicalData = clin)
  2438. #突变数目与risk的相关性
  2439. mutation_count=[email hidden]
  2441. #TCGA分组信息
  2442. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  2443. load(file="TCGA.RData")
  2445. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  2446. RS_mat <- read.table("RS_mat.txt",sep = "\t",header = T,stringsAsFactors = F,check.names = F,row.names=1)
  2447. bestmodel="Enet[alpha=0.3]"
  2448. RS_mat2=RS_mat[,bestmodel]
  2449. RS_mat2=RS_mat2[1:691]
  2450. group=RS_mat2
  2451. group=as.data.frame(cbind(ID=colnames(TCGA),group=group))
  2453. clin$ID <- substr(clin$Tumor_Sample_Barcode, 1, 16)
  2454. #匹配好risk的clin
  2455. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.1 genomic")
  2456. clin_new=read.csv("clin_new.csv")
  2457. clin_new=clin_new[order(clin_new$Tumor_Sample_Barcode),]
  2458. mutation_count=mutation_count[order(mutation_count$Tumor_Sample_Barcode),]
  2459. identical(mutation_count$Tumor_Sample_Barcode,clin_new$Tumor_Sample_Barcode)
  2461. mutation_count$risk=clin_new$X
  2463. #相关性分析
  2464. ggscatter(mutation_count, x = 'risk', y = 'total', size = 4,add = "reg.line",
  2465. color="lightblue",
  2466. add.params = list(color = "#77C034", fill = "#C5E99B", size = 1),
  2467. conf.int = TRUE)+
  2468. stat_cor(method = "spearman", label.x = 0, label.y = 1300, label.sep = "\n") +
  2469. ggtitle("")+
  2470. xlab("DisulfidpScore") +
  2471. ylab("Total mutation count")+
  2472. theme_classic()+
  2473. theme(legend.position = "none")
  2477. #risk组在TMB中的差异
  2478. tmb_table_wt_log = tmb(maf = maf)
  2479. head(tmb_table_wt_log)
  2481. #合并分组信息
  2482. dat=as.data.frame(tmb_table_wt_log[,c(1,4)])
  2483. rownames(dat)=dat$Tumor_Sample_Barcode
  2484. rownames(clin_new)=clin_new$Tumor_Sample_Barcode
  2485. dat2=merge(dat,clin_new,by=0)
  2486. #更新分组名称
  2487. dat2$sample=ifelse(dat2$X>median(dat2$X),"High","Low")
  2488. dat2$sample <- factor(dat2$sample, levels = c("Low","High"))
  2489. p=ggviolin(dat2, x="sample", y="total_perMB_log", fill="sample",
  2490. palette =c("#1F78B4","#E31A1C"),alpha = 0.5,
  2491. add ="boxplot",size = 0.5,
  2492. xlab = "TMB",
  2493. legend = "none")
  2494. p+stat_compare_means(aes(group = sample),
  2495. method = "wilcox.test",label = "p.signif",
  2496. label.x = 1.5,label.y = 1.5,
  2497. symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
  2498. symbols = c("***", "**", "*", "ns")))
  2502. #TMB的生存分析
  2503. library(data.table)
  2504. clinical=OS
  2505. clinical=rownames_to_column(clinical,"patient")
  2507. surv.df <- tmb_table_wt_log %>%
  2508. mutate(patient = str_sub(as.character(.$Tumor_Sample_Barcode),1,16)) %>%
  2509. distinct(patient, .keep_all = T) %>%
  2510. inner_join(clinical, by = "patient") %>%
  2511. mutate(group = if_else(total_perMB > median(total_perMB), "H-TMB","L-TMB")) %>%
  2512. select(patient,OS.time,OS,group)
  2513. colnames(surv.df)[2:3]=c(" times","status")
  2515. library(survival)
  2516. library(survminer)
  2517. times=surv.df$` times`
  2518. f <- survfit(Surv(times,status) ~ group, data = surv.df)
  2519. col=c("#E31A1C","#1F78B4")
  2520. ggsurvplot(f,pval=T,conf.int=F,xlab="Time(years)",palette=col,
  2521. legend.labs = c("H-TMB", "L-TMB"),legend=c(0.8,0.8))
  2525. #TMB结合risk的生存分析
  2526. clinical=OS
  2527. clinical=rownames_to_column(clinical,"patient")
  2529. surv.df <- tmb_table_wt_log %>%
  2530. mutate(patient = str_sub(as.character(.$Tumor_Sample_Barcode),1,16)) %>%
  2531. distinct(patient, .keep_all = T) %>%
  2532. inner_join(clinical, by = "patient") %>%
  2533. select(patient,OS.time,OS,total_perMB)
  2535. sample=surv.df$patient
  2536. group=group[group$ID %in% sample,]
  2537. surv.df$risk=group$group
  2538. colnames(surv.df)[4]="TMB"
  2540. rt=surv.df
  2541. #根据基因表达,对数据分组
  2542. a=if_else(surv.df$TMB > median(surv.df$TMB), "H-TMB","L-TMB")
  2543. b=ifelse(surv.df$risk >median(surv.df$risk), "High-risk", "Low-risk")
  2544. Type=paste(a,"+",b)
  2545. rt=cbind(rt,Type)
  2546. head(rt)
  2548. #生存差异统计
  2549. length=length(levels(factor(Type)))
  2550. diff=survdiff(Surv(OS.time, OS) ~Type,data = rt)
  2551. pValue=1-pchisq(diff$chisq,df=length-1)
  2552. if(pValue<0.001){
  2553. pValue="p<0.001"
  2554. }else{
  2555. pValue=paste0("p=",sprintf("%.03f",pValue))
  2556. }
  2557. fit <- survfit(Surv(OS.time, OS) ~ Type, data = rt)
  2558. head(fit)
  2560. #绘制生存曲线
  2561. col=c("#E31A1C","#FA9FB5","#1F78B4","lightblue")
  2563. surPlot=ggsurvplot(fit,
  2564. data=rt,
  2565. conf.int=F,
  2566. pval=pValue,pval.size=5,
  2567. palette=col,
  2568. legend = c(0.75,0.75),
  2569. xlab="Time(years)",
  2570. break.time.by = 1,
  2571. risk.table.title="",
  2572. risk.table=F,
  2573. risk.table.height=.25)
  2574. surPlot
  2578. #关键基因的互斥/共现
  2579. #High/Low
  2580. par(oma = c(3, 4, 5, 1))
  2581. somaticInteractions(maf = maf, top = 25,
  2582. genes=c("IDH1","TP53","EGFR","PTEN","TTN","ATRX","CIC","FUBP1"),
  2583. pvalue = c(0.01, 0.05),
  2584. colPal = "PiYG")
  2588. #关键基因的突变数量比较
  2589. gene=c("IDH1","TP53","EGFR","PTEN","TTN","ATRX","CIC","FUBP1")
  2591. #提取low risk组
  2592. low=maf.coad[maf.coad$sample=="Low",]
  2593. low=low[low$Hugo_Symbol %in% gene,]
  2594. low <- read.maf(low,clinicalData = clin)
  2596. #提取high risk组
  2597. high=maf.coad[maf.coad$sample=="High",]
  2598. high=high[high$Hugo_Symbol %in% gene,]
  2599. high <- read.maf(high,clinicalData = clin)
  2601. High.vs.Low <- mafCompare(m1 = high, m2 = low, m1Name = 'High',m2Name = 'Low', minMut = 5)
  2602. a=High.vs.Low[[1]]
  2603. forestPlot(mafCompareRes = High.vs.Low, pVal = 1, color = c('royalblue', 'maroon'), geneFontSize = 0.8)
  2609. #3.5 不同评分组的免疫治疗反应预测####
  2610. library(tidyverse)
  2611. library(ggplot2)
  2612. library(ggsci)
  2613. library(ggpubr)
  2614. library(ggsignif)
  2616. #基因矩阵
  2617. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  2618. load("TCGA.RData")
  2619. identical(colnames(TCGA),rownames(OS))
  2620. dat = TCGA
  2622. #标准化
  2623. exp<-t(apply(dat,1,function(x){x-(mean(x))}))
  2625. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/TIDE")
  2626. write.table(exp,file="TIDE.txt",sep="\t",quote=F)
  2628. #读取TIDE结果文件
  2629. TIDE<-read.csv("export.csv",header=T,check.names=F)
  2630. rownames(TIDE)=TIDE$Patient
  2631. sample=colnames(TCGA)
  2632. TIDE=TIDE[sample,]
  2634. #分组信息
  2635. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  2636. RS_mat <- read.table("RS_mat.txt",sep = "\t",header = T,stringsAsFactors = F,check.names = F,row.names=1)
  2637. bestmodel="Enet[alpha=0.3]"
  2638. RS_mat2=RS_mat[,bestmodel]
  2639. RS_mat2=RS_mat2[1:691]
  2640. group=ifelse(RS_mat2>median(RS_mat2),"High","Low")
  2642. TIDE$group=group
  2643. result=TIDE
  2645. #小提琴图展示结果
  2646. #1.TIDE小提琴图
  2647. my_comparisons <- list( c("Low", "High"))
  2648. p1 <- ggviolin(result, x = 'group', y = 'TIDE', fill = 'group',
  2649. palette = c("#A6D854","#FF1493"),alpha = 0.8,
  2650. add = 'boxplot', add.params = list(fill = "white")) +
  2651. stat_compare_means(comparisons = my_comparisons, label = "p.signif", bracket.size=0.5, tip.length = 0.01, method = 'wilcox.test')
  2652. p1
  2654. #2.Dysfunction小提琴图
  2655. p2 <- ggviolin(result, x = 'group', y = 'Dysfunction', fill = 'group',
  2656. palette = c("#A6D854","#FF1493"),alpha = 0.8,
  2657. add = 'boxplot', add.params = list(fill = "white")) +
  2658. stat_compare_means(comparisons = my_comparisons, label = "p.signif", bracket.size=0.5, tip.length = 0.02, method = 'wilcox.test')
  2659. p2
  2661. #3.Exclusion小提琴图
  2662. p3 <- ggviolin(result, x = 'group', y = 'Exclusion', fill = 'group',
  2663. palette = c("#A6D854","#FF1493"),alpha = 0.8,
  2664. add = 'boxplot', add.params = list(fill = "white")) +
  2665. stat_compare_means(comparisons = my_comparisons, label = "p.signif", bracket.size=0.5, tip.length = 0.02, method = 'wilcox.test')
  2666. p3
  2668. #4.MSI小提琴图
  2669. colnames(result)[6]
  2670. colnames(result)[6] <- c('MSI')
  2671. p4 <- ggviolin(result, x = 'group', y = 'MSI', fill = 'group',
  2672. palette = c("#A6D854","#FF1493"),alpha = 0.8,
  2673. add = 'boxplot', add.params = list(fill = "white")) +
  2674. stat_compare_means(comparisons = my_comparisons, label = "p.signif", bracket.size=0.5, tip.length = 0.02, method = 'wilcox.test')
  2675. p4
  2678. #百分比图
  2679. #分组计算百分比
  2680. mydata_summary <- TIDE %>%
  2681. group_by(group, Responder) %>%
  2682. summarise(count = n()) %>%
  2683. group_by(group) %>%
  2684. mutate(percent = count / sum(count) * 100)
  2686. #卡方检验
  2687. chi_square <- chisq.test(table(TIDE$group, TIDE$Responder))
  2688. p_value <- chi_square$p.value
  2689. p_value
  2691. #绘制百分比柱状堆叠图
  2692. mydata_summary$group <- factor(mydata_summary$group, levels = c("Low","High"))
  2693. ggplot(mydata_summary, aes(x = group, y = percent, fill = Responder)) +
  2694. geom_bar(stat = "identity", position = "stack",color = "#f3f4f4") +
  2695. annotate("text", x = 1.5, y = 105, label=expression(""~italic("P=1.204292e-29")), size = 4)+
  2696. geom_text(data = subset(mydata_summary, group == "High"), aes(label = paste0(round(percent), "%")),
  2697. position = position_stack(vjust = 0.5), color = "black", size = 3) +
  2698. geom_text(data = subset(mydata_summary, group == "Low"), aes(label = paste0(round(percent), "%")),
  2699. position = position_stack(vjust = 0.5), color = "black", size = 3) +
  2700. labs(title = "TCGA",x = "",y = "Percentage(%)") +
  2701. scale_fill_manual(values = c("True" = "#FF1493", "False" = "#A6D854")) +
  2702. theme_bw()+
  2703. theme(panel.grid = element_blank(),
  2704. plot.title = element_text(hjust = 0.5, face = "bold"),
  2705. plot.subtitle = element_text(hjust = 0.5, face = "italic"))+
  2706. guides(fill=guide_legend(reverse=TRUE))
  2709. #生存曲线
  2710. rt=cbind(OS,Responder=TIDE$Responder)
  2711. colnames(rt)[1:2]=c("fustat","futime")
  2712. diff=survdiff(Surv(futime, as.numeric(fustat)) ~ Responder, data=rt)
  2713. pValue=1-pchisq(diff$chisq, df=1)
  2714. if(pValue<0.001){
  2715. pValue="p<0.001"
  2716. }else{
  2717. pValue=paste0("p=", sprintf("%.03f",pValue))
  2718. }
  2719. fit <- survfit(Surv(futime, fustat) ~ Responder, data = rt)
  2720. surPlot=ggsurvplot(fit,
  2721. data=rt,
  2722. conf.int=F,
  2723. pval=pValue,
  2724. pval.size=6,
  2725. xlab="Time(years)",
  2726. legend.title="Cluster",
  2727. break.time.by = 5,
  2728. palette=c("#DECBE4", "#B3CDE3"))
  2731. #风险评分与Response
  2732. result$risk=RS_mat2
  2733. my_comparisons <- list( c("True", "False"))
  2734. #风险评分小提琴图
  2735. p12 <- ggviolin(result, x = 'Responder', y = 'risk', fill = 'Responder',
  2736. palette = c("#FF1493","#A6D854"),alpha = 0.8,
  2737. add = 'boxplot', add.params = list(fill = "white")) +
  2738. stat_compare_means(comparisons = my_comparisons, label = "p.signif", bracket.size=0.5, tip.length = 0.02, method = 'wilcox.test')
  2739. p12
  2743. #TIDE分析(meta-LGG-GBM)
  2744. library(tidyverse)
  2745. library(ggplot2)
  2746. library(ggsci)
  2747. library(ggpubr)
  2748. library(ggsignif)
  2750. #基因矩阵
  2751. setwd("/home/data/t060432/HRT/Disulfidptosis/data/2.2 meta")
  2752. load("meta.RData")
  2753. identical(colnames(meta),rownames(meta_pheno))
  2754. dat = meta
  2756. #标准化
  2757. exp<-t(apply(dat,1,function(x){x-(mean(x))}))
  2759. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.0 immune/TIDE")
  2760. write.table(exp,file="TIDE_LGG.txt",sep="\t",quote=F)
  2762. #读取TIDE结果文件
  2763. TIDE<-read.csv("TIDEresult_meta.csv",header=T,check.names=F)
  2764. rownames(TIDE)=TIDE$Patient
  2765. sample=colnames(meta)
  2766. TIDE=TIDE[sample,]
  2768. #分组信息
  2769. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  2770. sample=read.table("sample.txt")
  2771. RS_mat <- read.table("RS_mat.txt",sep = "\t",header = T,stringsAsFactors = F,check.names = F,row.names=1)
  2772. bestmodel="Enet[alpha=0.3]"
  2773. RS_mat2=RS_mat[,bestmodel]
  2774. group=RS_mat2
  2775. group_list=as.data.frame(cbind(sample=sample$sample,group=group))
  2776. rownames(group_list)=group_list$sample
  2778. sample=colnames(meta)
  2779. group_list=group_list[sample,]
  2780. group_list$group=as.numeric(group_list$group)
  2781. group_list=na.omit(group_list)
  2782. group_list$group=ifelse(group_list$group>median(group_list$group),"High","Low")
  2783. sample=rownames(group_list)
  2785. TIDE=TIDE[sample,]
  2786. TIDE$group=group_list$group
  2787. result=TIDE
  2792. #免疫治疗队列
  2793. library(dplyr)
  2794. library(tidyr)
  2795. library(ggpubr)
  2796. library(ggplot2)
  2797. library(survivalsvm)
  2798. library(survminer)
  2799. library(survival)
  2801. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.2 immutherapy")
  2802. load("GBM-PRJNA482620.RData")
  2803. risk=read.table("GBM-PRJNA482620_risk.txt",header=T)
  2805. identical(rownames(pdata),rownames(risk))
  2807. #有无响应百分堆叠图
  2808. mydata=as.data.frame(cbind(risk=risk,type=pdata$response))
  2809. mydata$group=ifelse(mydata[,1]>median(mydata[,1]), "High risk", "Low risk")
  2811. #分组计算百分比
  2812. mydata_summary <- mydata %>%
  2813. group_by(group, type) %>%
  2814. summarise(count = n()) %>%
  2815. group_by(group) %>%
  2816. mutate(percent = count / sum(count) * 100)
  2818. str(mydata_summary)
  2820. #卡方检验
  2821. chi_square <- chisq.test(table(mydata$group, mydata$type))
  2822. p_value <- chi_square$p.value
  2823. p_value
  2825. #绘制百分比柱状堆叠图
  2826. mydata_summary$group <- factor(mydata_summary$group, levels = c("Low risk","High risk"))
  2827. ggplot(mydata_summary, aes(x = group, y = percent, fill = type)) +
  2828. geom_bar(stat = "identity", position = "stack",color = "#f3f4f4") +
  2829. annotate("text", x = 1.5, y = 105, label=expression(""~italic("P=0.03288")), size = 4)+
  2830. geom_text(data = subset(mydata_summary, group == "High risk"), aes(label = paste0(round(percent), "%")),
  2831. position = position_stack(vjust = 0.5), color = "black", size = 3) +
  2832. geom_text(data = subset(mydata_summary, group == "Low risk"), aes(label = paste0(round(percent), "%")),
  2833. position = position_stack(vjust = 0.5), color = "black", size = 3) +
  2834. labs(title = "GBM-PRJNA482620",x = "",y = "Percentage(%)") +
  2835. scale_fill_manual(values = c("R" = "#FF1493", "N" = "#A6D854")) +
  2836. theme_bw()+
  2837. theme(panel.grid = element_blank(),
  2838. plot.title = element_text(hjust = 0.5, face = "bold"),
  2839. plot.subtitle = element_text(hjust = 0.5, face = "italic"))+
  2840. guides(fill=guide_legend(reverse=TRUE))
  2843. #生存曲线
  2844. rt=cbind(pdata[,c(1:2)],risk)
  2845. colnames(rt)[1:2]=c("fustat","futime")
  2846. Type=ifelse(rt$risk>median(rt$risk), "High", "Low")
  2847. rt=cbind(as.data.frame(rt), Type)
  2848. diff=survdiff(Surv(futime, as.numeric(fustat)) ~ Type, data=rt)
  2849. pValue=1-pchisq(diff$chisq, df=1)
  2850. if(pValue<0.001){
  2851. pValue="p<0.001"
  2852. }else{
  2853. pValue=paste0("p=", sprintf("%.03f",pValue))
  2854. }
  2855. fit <- survfit(Surv(futime, fustat) ~ Type, data = rt)
  2856. surPlot=ggsurvplot(fit,
  2857. data=rt,
  2858. conf.int=F,
  2859. pval=pValue,
  2860. pval.size=6,
  2861. xlab="Time(years)",
  2862. #ylab="Overall survival",
  2863. legend.title="Cluster",
  2864. break.time.by = 5,
  2865. palette=c("#FF1493", "#A6D854"))
  2868. #密度分布图
  2869. #计算中位数
  2870. median_risk_N <- median(mydata$risk[mydata$type == "N"])
  2871. median_risk_R <- median(mydata$risk[mydata$type == "R"])
  2872. type=c("N","R")
  2873. grp.mean=c(0.13,-0.095)
  2874. mu=as.data.frame(cbind(type,grp.mean))
  2876. p=ggplot(data=mydata, aes(x=risk, group=type, fill=type)) +
  2877. geom_density(#adjust=1.5,
  2878. alpha=0.5) +
  2879. ggprism::theme_prism(border = T)
  2880. p+scale_fill_manual(values=c("#A6D854","#FF1493"))
  2887. #3.6 不同评分组的药物敏感性预测####
  2888. #CGGA队列化疗比例
  2889. library(dplyr)
  2890. library(tidyr)
  2892. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  2893. load("CGGA693.RData")
  2894. pheno=CGGA693_pheno[,9:10]
  2895. colnames(pheno)=c("Radio_status","Chemo_status")
  2897. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  2898. RS_mat=read.table("RS_mat.txt",sep = "\t", row.names = 1,header = T,check.names = F)
  2899. best_mod = "Enet[alpha=0.3]"
  2900. RS_mat2=data.frame(RS_mat[,best_mod])
  2901. rownames(RS_mat2)=rownames(RS_mat)
  2902. colnames(RS_mat2)="risk"
  2904. pheno2 <- pheno %>% filter(!is.na(Chemo_status)) #Chemo_status/Radio_status
  2906. sample=rownames(pheno2)
  2908. mydata=as.data.frame(cbind(risk=RS_mat2[sample,],type=pheno2[,2]))
  2909. rownames(mydata)=sample
  2910. mydata$group=ifelse(mydata[,1]>median(mydata[,1]), "High risk", "Low risk")
  2911. mydata$type=ifelse(mydata$type=="0","No","Yes")
  2914. #分组计算百分比
  2915. mydata_summary <- mydata %>%
  2916. group_by(group, type) %>%
  2917. summarise(count = n()) %>%
  2918. group_by(group) %>%
  2919. mutate(percent = count / sum(count) * 100)
  2921. #卡方检验
  2922. chi_square <- chisq.test(table(mydata$group, mydata$type))
  2923. p_value <- chi_square$p.value
  2924. p_value
  2926. #绘制百分比柱状堆叠图
  2927. mydata_summary$group <- factor(mydata_summary$group, levels = c("Low risk","High risk"))
  2928. ggplot(mydata_summary, aes(x = group, y = percent, fill = type)) +
  2929. geom_bar(stat = "identity", position = "stack",color = "#f3f4f4") +
  2930. annotate("text", x = 1.5, y = 105, label=expression(""~italic("P=6.305645e-06")), size = 4)+
  2931. geom_text(data = subset(mydata_summary, group == "High risk"), aes(label = paste0(round(percent), "%")),
  2932. position = position_stack(vjust = 0.5), color = "black", size = 3) +
  2933. geom_text(data = subset(mydata_summary, group == "Low risk"), aes(label = paste0(round(percent), "%")),
  2934. position = position_stack(vjust = 0.5), color = "black", size = 3) +
  2935. labs(title = "CGGA693 Chemo_treated",x = "",y = "Percentage%") +
  2936. scale_fill_manual(values = c("Yes" = "#FF1493", "No" = "#A6D854")) +
  2937. theme_bw()+
  2938. theme(panel.grid = element_blank(),
  2939. plot.title = element_text(hjust = 0.5, face = "bold"),
  2940. plot.subtitle = element_text(hjust = 0.5, face = "italic"))+
  2941. guides(fill=guide_legend(reverse=TRUE))
  2946. library(pRRophetic)
  2947. library(ggplot2)
  2948. library(cowplot)
  2949. library(ggsignif)
  2950. library(ggsci)
  2951. library(tidyr)
  2952. library(dplyr)
  2953. library(ggpubr)
  2954. library(ggsci)
  2955. library(ggforce)
  2956. library(tidyverse)
  2957. library(ggpubr)
  2958. library(ggprism)
  2959. library(paletteer)
  2961. setwd("/home/data/t060432/HRT/Disulfidptosis/data/0. RData")
  2962. load("TCGA.RData")
  2963. dat=TCGA
  2965. data(PANCANCER_IC_Tue_Aug_9_15_28_57_2016)
  2966. GCP.drug <- unique(drugData2016$Drug.name)
  2968. #利用pRRopheticPredict函数进行药物敏感性分析
  2969. drug_list <- GCP.drug
  2970. results <- list()
  2972. for (drug in drug_list) {
  2973. tryCatch({
  2974. predictedPtype <- pRRopheticPredict(
  2975. testMatrix = as.matrix(dat),
  2976. drug = drug,
  2977. tissueType = "nervous_system",
  2978. selection = 1,
  2979. batchCorrect = "eb",
  2980. powerTransformPhenotype = T,
  2981. dataset = "cgp2016")
  2982. results[[drug]] <- predictedPtype
  2983. }, error = function(e) {
  2984. cat("Skipping drug:", drug, "\n")
  2985. })
  2986. }
  2988. #保存
  2989. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.3 drug sensitivity")
  2990. save(results,file="pRRophetic.RData")
  2991. load("pRRophetic.RData")
  2993. #整合结果
  2994. res.df <- do.call('rbind',results)
  2996. #提取风险评分分组
  2997. setwd("/home/data/t060432/HRT/Disulfidptosis/data/3.0 DisulfidpScore/Results")
  2998. RS_mat=read.table("RS_mat.txt",check.names = F)
  2999. best_mod = "Enet[alpha=0.3]"
  3000. RS_mat2 = RS_mat[,best_mod][1:691]
  3001. group=ifelse(RS_mat2>median(RS_mat2), "High", "Low")
  3003. data=rbind(res.df,group)
  3004. data=as.data.frame(t(data))
  3005. data[, 1:237] <- apply(data[, 1:237], 2, as.numeric)
  3007. #可视化
  3008. #Temozolomide\Cisplatin
  3009. ggboxplot(data, x = "group", y = "Vinblastine",
  3010. color = "group", notch = TRUE,
  3011. palette = c("#A6D854", "#FF1493"),alpha = 0.75,
  3012. legend = "none",size = 1,fill = "group")+
  3013. stat_compare_means(aes(group=group),label = "p.signif") +
  3014. labs(title = " ", x = NULL, y = "Vinblastine") +
  3015. theme(plot.title = element_text(hjust = 0.5))
  3019. #循环237个药物
  3020. drug <- colnames(data)[1:237]
  3021. p_value_results <- data.frame(drug = character(), p_value = numeric(),
  3022. group = character(), high_mean = numeric(),
  3023. low_mean = numeric(), stringsAsFactors = FALSE)
  3025. for (i in drug) {
  3026. ic50_high <- data[data$group == "High", i]
  3027. ic50_low <- data[data$group == "Low", i]
  3028. high_mean <- mean(ic50_high)
  3029. low_mean <- mean(ic50_low)
  3030. t_test_result <- t.test(ic50_high, ic50_low)
  3031. p_value <- t_test_result$p.value
  3032. if (high_mean > low_mean) {
  3033. group <- "High"} else {
  3034. group <- "Low"}
  3035. p_value_results <- rbind(p_value_results, data.frame(drug = i, p_value = p_value,
  3036. group = group, high_mean = high_mean,
  3037. low_mean = low_mean))}
  3039. #提取p<0.01的结果
  3040. significant_results <- p_value_results[p_value_results$p_value < 0.01, ] #199/209(交集后结果一致)
  3044. #相关性分析
  3045. data=rbind(res.df,risk=RS_mat2)
  3046. data=as.data.frame(t(data))
  3048. data_use=data
  3049. target_gene <- 'risk'
  3050. target_column <- data_use[,target_gene]
  3052. #单基因
  3053. cor_R <- cor(x = target_column,y = data_use[,1],method = 'pearson')
  3054. cor_P <- cor.test(x = target_column,y = data_use[,1])$p.value
  3055. result_1 <- data.frame(target_gene = target_gene,
  3056. gene_symbol = 'risk', #对应基因名
  3057. cor_R = cor_R,
  3058. cor_P = cor_P)
  3060. #批量计算基因间的相关性
  3061. result <- data.frame("target_gene" = character(),
  3062. "gene_symbol" = character(),
  3063. "cor_R" = numeric(),
  3064. "cor_P" = numeric())
  3066. gene_list <- colnames(data_use)[1:237]
  3067. for (gene in gene_list) {
  3068. print(gene)
  3069. cor_R <- cor(x = target_column,y = data_use[,gene],method = 'pearson')
  3070. cor_P <- cor.test(x = target_column,y = data_use[,gene])$p.value
  3071. temp_result <- data.frame(target_gene = target_gene,
  3072. gene_symbol = gene,
  3073. cor_R = cor_R,
  3074. cor_P = cor_P)
  3075. result <- rbind(result, temp_result)
  3076. }
  3077. resCor <- na.omit(result[(result$cor_R > 0.5|result$cor_R < -0.5 & result$cor_P < 0.01), ]) #95
  3080. drugnames=resCor$gene_symbol
  3081. selected_drugs <- significant_results[significant_results$drug %in% drugnames, ] #95
  3082. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.3 drug sensitivity")
  3083. write.table(selected_drugs,file="ic50.txt",quote=F,sep="\t")
  3086. #相关性散点图
  3087. ggscatter(data, x = 'risk', y = 'Vinblastine', size = 4,add = "reg.line",
  3088. color="lightblue",
  3089. add.params = list(color = "#77C034", fill = "#C5E99B", size = 1),
  3090. conf.int = TRUE)+
  3091. stat_cor(method = "spearman", label.x = 0, label.y = -3.5, label.sep = "\n") +
  3092. ggtitle("")+
  3093. xlab("DisulfidpScore") +
  3094. ylab("IC50(Vinblastine)")+
  3095. theme_classic()+
  3096. theme(legend.position = "none")
  3100. #可视化
  3101. setwd("/home/data/t060432/HRT/Disulfidptosis/data/4.3 drug sensitivity")
  3102. dat=read.table("ic50plot.txt",sep = "\t", header = TRUE)
  3103. dat_sorted <- arrange(dat, group, ic50)
  3105. ggplot(dat_sorted, aes(x = reorder(drug, ic50), y = ic50, fill = group)) +
  3106. geom_col() +
  3107. labs(y="IC50(H)/IC50(L)-1",x="")+
  3108. scale_fill_manual(values = c("High" = "#FF1493", "Low" = "#A6D854"))+
  3109. theme_classic() +
  3110. theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  3113. #筛选IC50最低的10个药物(差异箱式图)
  3114. needata=selected_drugs[order(selected_drugs$low_mean),]
  3115. top10drug=needata$drug[1:10]
  3117. top10data=data[,top10drug]
  3118. group=ifelse(RS_mat2>median(RS_mat2), "High", "Low")
  3119. info=as.data.frame(cbind(sample=rownames(data),group))
  3120. colnames(info)=c("Sample","Type")
  3121. Exp_plot=top10data
  3122. Exp_plot$sam=info$Type
  3123. Exp_plot$sam <- factor(Exp_plot$sam, levels = c("Low","High"))
  3124. expr_use <- na.omit(Exp_plot)
  3125. expr_use <-as.data.frame(expr_use)
  3126. expr_use_long <- gather(expr_use, gene, Expression, -sam)
  3127. table(expr_use_long$gene)
  3128. colnames(expr_use_long) <- c("Group","gene","Expression")
  3130. my_comparisons <- list(c("Low","High"))
  3132. #可视化
  3133. ggboxplot(expr_use_long, x = "gene", y = "Expression",
  3134. color = "Group",
  3135. palette = c("#A6D854", "#FF1493"),alpha = 0.7,notch = TRUE,
  3136. #legend = "top",
  3137. size = 1,fill = NULL)+
  3138. stat_compare_means(aes(group=Group),label = "p.signif") +
  3139. labs(title = " ", x = NULL, y = "IC50") +
  3140. theme(plot.title = element_text(hjust = 0.5)) +
  3141. theme_bw() +
  3142. theme(panel.grid.major = element_blank(),
  3143. panel.grid.minor = element_blank(),
  3144. axis.title.x = element_blank(),
  3145. legend.position = "top")
  3148. #筛选IC50最低的10个药物(相关性)
  3149. cordata=resCor[resCor$gene_symbol %in% top10drug,]
  3150. cordata=cordata[,c(2:4)]
  3151. cordata$cor_P=-log10(cordata$cor_P)
  3153. ggplot(cordata, aes(x = gene_symbol, y = cor_R)) +
  3154. geom_col(width = 0.04, fill = '#FF1493') +
  3155. geom_point(aes(size = abs(cor_P)), color = '#FF1493') +
  3156. theme_bw() +
  3157. theme(panel.grid.major = element_blank(),
  3158. panel.grid.minor = element_blank(),
  3159. axis.title.x = element_blank(),
  3160. legend.position = "top") +
  3161. scale_fill_manual(values = unique(cordata$gene_symbol)) +
  3162. geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  3163. labs(size = "-log10(pvalue)")

workline2.R at commit 64e79df, no license · at the source

Overview

Authors: Ruiting Huang1, Hailin Li1, Yijing Zhong1, Paimin Zhuo1, Yibei Wang1, Aoting Yang1, Yu Zhang2, Jiao Li3, Ruiquan Xu4,5, Quhuan Li1
ORCID iDs: Paimin Zhuo, Quhuan Li
  1. School of Biology and Biological Engineering, South China University of Technology, Guangzhou, Guangdong 510006, China
  2. The First Clinical School of Gannan Medical University, Ganzhou, Jiangxi 341000, China
  3. Gannan Medical University, Ganzhou, Jiangxi province 341000, China
  4. Department of Urology, First Affiliated Hospital of Gannan Medical University, Ganzhou, Jiangxi 341000, China
  5. Longnan First People's Hospital, Longnan, Jiangxi 341700, China
Institutions: South China University of Technology (China); Gannan Medical University (China); First Affiliated Hospital of Gannan Medical University (China)
Journal: iScience, volume 29, issue 5, article 115657
Dates: received 6 September 2025; accepted 6 April 2026; published online 8 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.isci.2026.115657 · PMID 42164521 · PMCID PMC13185775 · OpenAlex W7151849519
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), other condition (population), clinical / translational (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing
Keywords: bioinformatics, cancer, artificial intelligence applications
Topic: Glioma Diagnosis and Treatment (Genetics, Medicine), according to OpenAlex
Funding: Natural Science Foundation of Guangdong Province (2023A1515010829, 2021A1515010040); National Natural Science Foundation of China (32271360, 31870928)
Citations: not cited yet (Europe PMC); 47 references in the paper

Abstract

Glioma prognosis is challenged by tumor heterogeneity and lack of biomarkers. Disulfidptosis, a novel cell death mechanism induced by disulfide stress, remains poorly understood in gliomas. This study analyzed eight glioma cohorts, identifying two disulfidptosis patterns with distinct genomic alterations, immune microenvironments, and clinical outcomes. A prognostic model—DisulfidpScore—was developed using machine learning, demonstrating robust predictive ability for survival. Crucially, single-cell profiling and virtual knockout analysis revealed elevated disulfidptosis in glioblastoma astrocytes and identified IQGAP1 as a key driver that modulates gene networks governing the cell cycle and neuron-glia interactions. High DisulfidpScore scores correlated with immunosuppressive microenvironments and poorer prognosis but increased chemotherapy sensitivity, whereas low scores indicated better survival and immunotherapy response. The model supports prognostic stratification and personalized treatment for glioma patients.

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

Repositories

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

clarozhong/Rdata

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 64e79df796231c4bb8b95430241efb84860dd453, 29 January 2026
Languages: R (4)
Size: 5 files, 4 scripts
Software Heritage: not archived
Found in: “Data and code availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (3 files), clusterProfiler (2 files), cowplot (2 files), ggpubr (2 files), survival (2 files), tidyverse (2 files), broom (1 file), circlize (1 file), ComplexHeatmap (1 file), data.table (1 file), edgeR (1 file), Harmony (1 file), igraph (1 file), limma (1 file), pheatmap (1 file), reshape2 (1 file), Seurat (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
5 files
At the source: github.com/clarozhong/Rdata.git

IOBR/IOBR

License: GPL
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 635effe4d05f37a244c04f2d5cf2a5e73c38a798, 23 July 2026
Languages: R (98)
Size: 310 files, 98 scripts
Software Heritage: not archived
Found in: the resources table
Holds: README, environment (DESCRIPTION, Dockerfile), continuous integration, documentation, 2 notebooks
Not found: license file, CITATION.cff, tests
Tools: tidyverse (42 files), ggplot2 (21 files), survival (9 files), reshape2 (5 files), glmnet (4 files), ComplexHeatmap (3 files), DESeq2 (3 files), ggpubr (3 files), limma (3 files), Seurat (3 files), clusterProfiler (2 files), patchwork (2 files), circlize (1 file), pROC (1 file), psych (1 file), WGCNA (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
99 files
At the source: github.com/IOBR/IOBR

paulgeeleher/pRRophetic

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 0be5ac3c49e4f34ea563e12a79572e7460423c4f, 16 December 2023
Languages: R (8)
Size: 40 files, 8 scripts
Software Heritage: not archived
Found in: the resources table
Holds: README, environment (pRRophetic/DESCRIPTION), documentation
Not found: license file, CITATION.cff, tests, continuous integration
Tools: car (3 files), ggplot2 (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
9 files
At the source: github.com/paulgeeleher/pRRophetic

cailab-tamu/scTenifoldKnk

License: GPL
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: e46221bc988b9db878c84df979fd636b55b59861, 24 September 2026
Languages: R (49), Python (5), Shell (2), MATLAB (1)
Size: 428 files, 57 scripts
Software Heritage: not archived
Found in: the resources table
Holds: README, environment (DESCRIPTION), documentation
Not found: license file, CITATION.cff, tests, continuous integration
Tools: ggplot2 (27 files), igraph (13 files), ComplexHeatmap (9 files), patchwork (8 files), reshape2 (7 files), Seurat (6 files), circlize (4 files), Harmony (3 files), NumPy (3 files), NetworkX (1 file), SciPy (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
58 files
At the source: github.com/cailab-tamu/scTenifoldKnk

The paper's code and data availability statement is in the Data section.

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:

  • 4 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 167 scripts, each with its path and the digest of its content;
  • 27 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 and code availability

This paper analyzes existing, publicly available data, accessible at TCGA: TCGA-Gliomas; GTEx: Normal brain; CGGA: CGGA325, CGGA693; GlioVis: GSE16011 (https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE16011), GSE108474 (https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE108474), GSE4271 (https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE4271), GSE4412 (https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE4412), and E_MATE_3892; and TIGER: GBM-PRJNA482620.

The code supporting the findings of this study is available in GitHub at https://github.com/clarozhong/Rdata.git.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

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, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 10 authors, 3 keywords, 2 funders, 45 references.

Cite

This paper

Huang, R., Li, H., Zhong, Y., Zhuo, P., Wang, Y., Yang, A., Zhang, Y., Li, J., Xu, R., & Li, Q. (2026). Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug. iScience, 29(5), 115657. https://doi.org/10.1016/j.isci.2026.115657

BibTeX

@article{huang2026integration,
author = {Huang, Ruiting and Li, Hailin and Zhong, Yijing and Zhuo, Paimin and Wang, Yibei and Yang, Aoting and Zhang, Yu and Li, Jiao and Xu, Ruiquan and Li, Quhuan},
title = {{Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug}},
journal = {iScience},
year = {2026},
month = apr,
volume = {29},
number = {5},
pages = {115657},
publisher = {Elsevier},
issn = {2589-0042},
doi = {10.1016/j.isci.2026.115657},
url = {https://doi.org/10.1016/j.isci.2026.115657},
pmid = {42164521},
pmcid = {PMC13185775}
}

RIS

TY - JOUR
AU - Huang, Ruiting
AU - Li, Hailin
AU - Zhong, Yijing
AU - Zhuo, Paimin
AU - Wang, Yibei
AU - Yang, Aoting
AU - Zhang, Yu
AU - Li, Jiao
AU - Xu, Ruiquan
AU - Li, Quhuan
TI - Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug
T2 - iScience
J2 - iScience
PY - 2026
DA - 2026/04/08
VL - 29
IS - 5
SP - 115657
SN - 2589-0042
PB - Elsevier
DO - 10.1016/j.isci.2026.115657
UR - https://doi.org/10.1016/j.isci.2026.115657
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.isci.2026.115657",
"type": "article-journal",
"title": "Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug",
"container-title": "iScience",
"author": [
{
"family": "Huang",
"given": "Ruiting"
},
{
"family": "Li",
"given": "Hailin"
},
{
"family": "Zhong",
"given": "Yijing"
},
{
"family": "Zhuo",
"given": "Paimin"
},
{
"family": "Wang",
"given": "Yibei"
},
{
"family": "Yang",
"given": "Aoting"
},
{
"family": "Zhang",
"given": "Yu"
},
{
"family": "Li",
"given": "Jiao"
},
{
"family": "Xu",
"given": "Ruiquan"
},
{
"family": "Li",
"given": "Quhuan"
}
],
"container-title-short": "iScience",
"volume": "29",
"issue": "5",
"page": "115657",
"DOI": "10.1016/j.isci.2026.115657",
"PMID": "42164521",
"PMCID": "PMC13185775",
"ISSN": "2589-0042",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.isci.2026.115657",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
8
]
]
}
}

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.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: pROC, survival, WGCNA, 20 other tools, other condition
[2] 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: pROC, WGCNA, Harmony, 20 other tools, other condition
[3] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: glmnet, survival, WGCNA, 19 other tools
[4] doi:10.1016/j.xcrm.2026.102682 [code]
TET CpG sequence-context-specific DNA demethylation shapes progression of IDH-mutant gliomas.
Journal: Cell reports. Medicine
In common: glmnet, pROC, survival, 17 other tools, other condition
[5] doi:10.3390/ijms27093997 [code]
Coordinated Multicellular Immune Programs and Drug Targets Revealed by Single-Cell Analysis in Driver-Mutated NSCLC.
Journal: International journal of molecular sciences
In common: glmnet, survival, Harmony, 15 other tools, other condition
[6] 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: WGCNA, Harmony, edgeR, 15 other tools
[7] doi:10.1038/s41586-026-10214-2 [code]
Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: Harmony, edgeR, car, 16 other tools, other condition
[8] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Harmony, edgeR, limma, 16 other tools
[9] doi:10.1016/j.isci.2026.115573 [code]
Female cortical cellular mosaicism underlies shared MeCP2 and PCB impacted gene pathways.
Journal: iScience
In common: WGCNA, edgeR, limma, 15 other tools
[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, edgeR, car, 15 other tools, 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.