百度360必应搜狗淘宝本站头条
当前位置:网站首页 > 热门文章 > 正文

手把手带你复现NC图表之Figure6 ncl数据处理及绘图

bigegpt 2024-10-30 01:48 8 浏览

复现文章信息:

文章题目:Single-cell analysis reveals prognostic fibroblast subpopulations linked to molecular and immunological subtypes of lung cancer
期刊:Nature Communications
日期:2023年1月31日
DOI: 10.1038/s41467-023-35832-6

复现图——Figure 6

使用多个NSCLC队列进行生存分析

R包载入与数据准备

代码如下:

options(stringsAsFactors = F)
pardefault <- par()
library(Seurat)
library(tidyverse)
library(ggpubr)
library(ggsci)
library(survival)
library(survminer)
library(forestmodel)
data_directory <- "H:\\文献复现\\6\\" 
source(paste0(data_directory, "0_NewFunctions.R"))

load(paste0(data_directory, "BulkData_Zenodo.Rdata"))
load(paste0(data_directory, "IntegratedFibs_Zenodo.Rdata"))

LUAD单变量COX回归

为了检测肌成纤维细胞丰度作为LUAD患者分层预后生物标志物的可能性,使用TCGA-LUAD数据集测试

#为了研究将肌成纤维细胞丰度作为LUAD患者分层的预后生物标志物的潜力,我们使用TCGA-LUAD数据集作为测试队列,以确定将样本分类为肌成纤维细胞高和低的最佳阈值

ys2test <- rev(seq(1,10, by = 1))
cox.zph_OK <- list()
for(DS in levels(Merged.LUADtraits$Dataset.factor)[]){
  surv_data <- Merged.LUADtraits[Merged.LUADtraits$Sample.Subtype == "LUAD" &
                               Merged.LUADtraits$Dataset.factor == DS, ]
  
  cox.zph_p = 0
  for(i in ys2test){
    if (cox.zph_p < 0.05) {
      split_data <- survSplit(Surv(OS_YEARS, OS) ~ 
                                Myo_Fibs.pct,
                              zero = -0.1,
                              data = surv_data[!is.na(surv_data$OS_YEARS), ],
                              cut = i, episode = "tgroup", id = "id")
      model.coxph2 <- coxph(Surv(OS_YEARS, OS) ~ 
                              Myo_Fibs.pct,
                            data = split_data[split_data$tgroup == 1,])
      cox.zph_res = cox.zph(model.coxph2)
      cox.zph_res.sum = c(i, cox.zph_res$table[1,3])
      cox.zph_p <- cox.zph_res$table[1,3]
      names(cox.zph_res.sum) <- c("OS_YEARS", "cox.zph_p")
    }
    cox.zph_OK[[DS]] <- cox.zph_res.sum
  }
}
do.call(rbind, cox.zph_OK)
max_year <- min(do.call(rbind, cox.zph_OK)[,1])
#分别对每个数据集进行Cox回归分析
split_data <- list()
for(DS in levels(Merged.LUADtraits$Dataset.factor)[]){
  surv_data <- Merged.LUADtraits[Merged.LUADtraits$Sample.Subtype == "LUAD" &
                               Merged.LUADtraits$Dataset.factor == DS, ]
  split_data[[DS]] <- survSplit(Surv(OS_YEARS, OS) ~ 
                                  Myo_Fibs.pct + Alveolar_Fibs.pct + Adventitial_Fibs.pct,
                                zero = -0.1,
                                data = surv_data[!is.na(surv_data$OS_YEARS), ],
                                cut = max_year, episode = "tgroup", id = "id")
}
Myo_cox.res <- list()
Myo_cox.res[["TCGA"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                  scale(Myo_Fibs.pct),
                                data =  split_data[["TCGA"]][split_data[["TCGA"]]$tgroup == 1,])
Myo_cox.res[["GSE68465"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                     scale(Myo_Fibs.pct),
                                   data =  split_data[["GSE68465"]][split_data[["GSE68465"]]$tgroup == 1,])
Myo_cox.res[["GSE31210"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                     scale(Myo_Fibs.pct),
                                   data =  split_data[["GSE31210"]][split_data[["GSE31210"]]$tgroup == 1,])
Myo_cox.res[["GSE72094"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                      scale(Myo_Fibs.pct),
                                    data =  split_data[["GSE72094"]][split_data[["GSE72094"]]$tgroup == 1,])
forestmodel::forest_model(model_list = Myo_cox.res, panels = panels)
Myo_LUAD_Zscores <- list()
for(DS in names(Myo_cox.res)){
  Myo_LUAD_Zscores[[DS]] <- summary(Myo_cox.res[[DS]])$coefficients[4:5]
}
Myo_LUAD_pvals <- list()
for(DS in names(Myo_cox.res)){
  Myo_LUAD_pvals[[DS]] <- summary(Myo_cox.res[[DS]])$coefficients[5]
}
Myo_LUAD_pvals
metap::sumz(unlist(Myo_LUAD_pvals))
Alveolar_cox.res <- list()
Alveolar_cox.res[["TCGA"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                  scale(Alveolar_Fibs.pct),
                                data =  split_data[["TCGA"]][split_data[["TCGA"]]$tgroup == 1,])
Alveolar_cox.res[["GSE68465"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                     scale(Alveolar_Fibs.pct),
                                   data =  split_data[["GSE68465"]][split_data[["GSE68465"]]$tgroup == 1,])
Alveolar_cox.res[["GSE31210"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                     scale(Alveolar_Fibs.pct),
                                   data =  split_data[["GSE31210"]][split_data[["GSE31210"]]$tgroup == 1,])
Alveolar_cox.res[["GSE72094"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                      scale(Alveolar_Fibs.pct),
                                    data =  split_data[["GSE72094"]][split_data[["GSE72094"]]$tgroup == 1,])

forestmodel::forest_model(model_list = Alveolar_cox.res, panels = panels)
Alveolar_LUAD_Zscores <- list()
for(DS in names(Alveolar_cox.res)){
  Alveolar_LUAD_Zscores[[DS]] <- summary(Alveolar_cox.res[[DS]])$coefficients[4:5]
}
Alveolar_LUAD_pvals <- list()
for(DS in names(Alveolar_cox.res)){
  Alveolar_LUAD_pvals[[DS]] <- summary(Alveolar_cox.res[[DS]])$coefficients[5]
}
Alveolar_LUAD_pvals
metap::sumz(unlist(Alveolar_LUAD_pvals))
Adventitial_cox.res <- list()
Adventitial_cox.res[["TCGA"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                  scale(Adventitial_Fibs.pct),
                                data =  split_data[["TCGA"]][split_data[["TCGA"]]$tgroup == 1,])
Adventitial_cox.res[["GSE68465"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                     scale(Adventitial_Fibs.pct),
                                   data =  split_data[["GSE68465"]][split_data[["GSE68465"]]$tgroup == 1,])
Adventitial_cox.res[["GSE31210"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                     scale(Adventitial_Fibs.pct),
                                   data =  split_data[["GSE31210"]][split_data[["GSE31210"]]$tgroup == 1,])
Adventitial_cox.res[["GSE72094"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                      scale(Adventitial_Fibs.pct),
                                    data =  split_data[["GSE72094"]][split_data[["GSE72094"]]$tgroup == 1,])
forestmodel::forest_model(model_list = Adventitial_cox.res, panels = panels)
Adventitial_LUAD_Zscores <- list()
for(DS in names(Adventitial_cox.res)){
  Adventitial_LUAD_Zscores[[DS]] <- summary(Adventitial_cox.res[[DS]])$coefficients[4:5]
}

LUAD_zScores.df <- data.frame(Z = do.call(rbind, Myo_LUAD_Zscores)[,1],
                              p = do.call(rbind, Myo_LUAD_Zscores)[,2],
                              SubPop = "Myo",
                              DS = rownames(do.call(rbind, Myo_LUAD_Zscores)))
LUAD_zScores.df <- rbind(LUAD_zScores.df,
                         data.frame(Z = do.call(rbind, Alveolar_LUAD_Zscores)[,1],
                                    p = do.call(rbind, Alveolar_LUAD_Zscores)[,2],
                                    SubPop = "Alveolar",
                                    DS = rownames(do.call(rbind, Alveolar_LUAD_Zscores))))
LUAD_zScores.df <- rbind(LUAD_zScores.df,
                         data.frame(Z = do.call(rbind, Adventitial_LUAD_Zscores)[,1],
                                    p = do.call(rbind, Adventitial_LUAD_Zscores)[,2],
                                    SubPop = "Adventitial",
                                    DS = rownames(do.call(rbind, Adventitial_LUAD_Zscores))))
LUAD_zScores.df$Subtype <- "LUAD"
reshape2::melt(LUAD_zScores.df, id.vars = c("SubPop", "DS"), measure.vars = "Z") %>%
  ggplot(aes(y = value, x = SubPop, fill = SubPop)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(size = 1) + 
  theme_pubr(base_size = 7) +
  scale_fill_manual(values = Fibs_col.palette) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ylab("Good Survival <- CoxPH Zscore -> Poor Survival") +
  theme(axis.title.y = element_blank()) +
  coord_flip()

LUSC单变量COX回归

代码如下:

ys2test <- rev(seq(1,10, by = 1))
cox.zph_OK <- list()
for(DS in unique(Merged.LUSCtraits$Dataset)){
  surv_data <- Merged.LUSCtraits[Merged.LUSCtraits$Sample.Subtype == "LUSC" &
                                   Merged.LUSCtraits$Dataset == DS, ]
  cox.zph_p = 0
  for(i in ys2test){
    if (cox.zph_p < 0.05) {
      split_data <- survSplit(Surv(OS_YEARS, OS) ~ 
                                Myo_Fibs.pct,# + Stage_4cat + Age,
                              zero = -0.1,
                              data = surv_data[!is.na(surv_data$OS_YEARS), ],
                              cut = i, episode = "tgroup", id = "id")
      model.coxph2 <- coxph(Surv(OS_YEARS, OS) ~ 
                              Myo_Fibs.pct,
                            data = split_data[split_data$tgroup == 1,])
      cox.zph_res = cox.zph(model.coxph2)
      cox.zph_res.sum = c(i, cox.zph_res$table[1,3])
      cox.zph_p <- cox.zph_res$table[1,3]
      names(cox.zph_res.sum) <- c("OS_YEARS", "cox.zph_p")
    }
    cox.zph_OK[[DS]] <- cox.zph_res.sum
  }
}
do.call(rbind, cox.zph_OK)
max_year <- 4
#分别对每个数据集进行Cox回归分析
split_data <- list()
for(DS in unique(Merged.LUSCtraits$Dataset)[]){
  surv_data <- Merged.LUSCtraits[Merged.LUSCtraits$Sample.Subtype == "LUSC" &
                                   Merged.LUSCtraits$Dataset == DS, ]
  split_data[[DS]] <- survSplit(Surv(OS_YEARS, OS) ~ 
                                  Myo_Fibs.pct + Alveolar_Fibs.pct + Adventitial_Fibs.pct,
                                zero = -0.1,
                                data = surv_data[!is.na(surv_data$OS_YEARS), ],
                                cut = max_year, episode = "tgroup", id = "id")
}
Myo_cox.res <- list()
Myo_cox.res[["TCGA"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                  scale(Myo_Fibs.pct),
                                data =  split_data[["TCGA"]][split_data[["TCGA"]]$tgroup == 1,])
Myo_cox.res[["GSE157009"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                     scale(Myo_Fibs.pct),
                                   data =  split_data[["GSE157009"]][split_data[["GSE157009"]]$tgroup == 1,])
Myo_cox.res[["GSE157010"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                     scale(Myo_Fibs.pct),
                                   data =  split_data[["GSE157010"]][split_data[["GSE157010"]]$tgroup == 1,])
Myo_cox.res[["GSE4573"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                      scale(Myo_Fibs.pct),
                                    data =  split_data[["GSE4573"]][split_data[["GSE4573"]]$tgroup == 1,])
forestmodel::forest_model(model_list = Myo_cox.res, panels = panels)
Myo_LUSC_Zscores <- list()
for(DS in names(Myo_cox.res)){
  Myo_LUSC_Zscores[[DS]] <- summary(Myo_cox.res[[DS]])$coefficients[4:5]
}
Alveolar_cox.res <- list()
Alveolar_cox.res[["TCGA"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                  scale(Alveolar_Fibs.pct),
                                data =  split_data[["TCGA"]][split_data[["TCGA"]]$tgroup == 1,])
Alveolar_cox.res[["GSE157009"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                       scale(Alveolar_Fibs.pct),
                                     data =  split_data[["GSE157009"]][split_data[["GSE157009"]]$tgroup == 1,])
Alveolar_cox.res[["GSE157010"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                       scale(Alveolar_Fibs.pct),
                                     data =  split_data[["GSE157010"]][split_data[["GSE157010"]]$tgroup == 1,])
Alveolar_cox.res[["GSE4573"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                     scale(Alveolar_Fibs.pct),
                                   data =  split_data[["GSE4573"]][split_data[["GSE4573"]]$tgroup == 1,])
forestmodel::forest_model(model_list = Alveolar_cox.res, panels = panels)
Alveolar_LUSC_Zscores <- list()
for(DS in names(Alveolar_cox.res)){
  Alveolar_LUSC_Zscores[[DS]] <- summary(Alveolar_cox.res[[DS]])$coefficients[4:5]
}
Adventitial_cox.res <- list()
Adventitial_cox.res[["TCGA"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                  scale(Adventitial_Fibs.pct),
                                data =  split_data[["TCGA"]][split_data[["TCGA"]]$tgroup == 1,])
Adventitial_cox.res[["GSE157009"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                       scale(Adventitial_Fibs.pct),
                                     data =  split_data[["GSE157009"]][split_data[["GSE157009"]]$tgroup == 1,])
Adventitial_cox.res[["GSE157010"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                       scale(Adventitial_Fibs.pct),
                                     data =  split_data[["GSE157010"]][split_data[["GSE157010"]]$tgroup == 1,])
Adventitial_cox.res[["GSE4573"]]  <- coxph(Surv(OS_YEARS, OS) ~ 
                                     scale(Adventitial_Fibs.pct),
                                   data =  split_data[["GSE4573"]][split_data[["GSE4573"]]$tgroup == 1,])
forestmodel::forest_model(model_list = Adventitial_cox.res, panels = panels)
Adventitial_LUSC_Zscores <- list()
for(DS in names(Adventitial_cox.res)){
  Adventitial_LUSC_Zscores[[DS]] <- summary(Adventitial_cox.res[[DS]])$coefficients[4:5]
}
LUSC_zScores.df <- data.frame(Z = do.call(rbind, Myo_LUSC_Zscores)[,1],
                              p = do.call(rbind, Myo_LUSC_Zscores)[,2],
                              SubPop = "Myo",
                              DS = rownames(do.call(rbind, Myo_LUSC_Zscores)))
LUSC_zScores.df <- rbind(LUSC_zScores.df,
                         data.frame(Z = do.call(rbind, Alveolar_LUSC_Zscores)[,1],
                                    p = do.call(rbind, Alveolar_LUSC_Zscores)[,2],
                                    SubPop = "Alveolar",
                                    DS = rownames(do.call(rbind, Alveolar_LUSC_Zscores))))
LUSC_zScores.df <- rbind(LUSC_zScores.df,
                         data.frame(Z = do.call(rbind, Adventitial_LUSC_Zscores)[,1],
                                    p = do.call(rbind, Adventitial_LUSC_Zscores)[,2],
                                    SubPop = "Adventitial",
                                    DS = rownames(do.call(rbind, Adventitial_LUSC_Zscores))))
LUSC_zScores.df$Subtype <- "LUSC"

Figure 6A

散点图显示使用不同截点按肌成纤维细胞丰度对TCGA中LUAD25样本进行二分的标准化对数秩生存统计量(与总生存率的相关性)的变化。

Fig_6A <- 
  data.frame(stat = LUAD.Opt_cut$Myo_Fibs.pct$stats,
             x = LUAD.Opt_cut$Myo_Fibs.pct$cuts) %>%
  ggplot(aes(x = x, y = stat, colour = x < LUAD.Opt_cut$cutpoint[1,1])) +
  geom_point(size = 2) +
  theme_pubr(base_size = 15) +
  theme(legend.position = "none") +
  scale_colour_manual(values =  c(pal_d3()(3)[3], "grey70")) +
  xlab("Myo (% of all fibroblasts)") +
  ylab("Standardized Log-Rank statistic") +
  geom_vline(xintercept = LUAD.Opt_cut$cutpoint[1,1], 
             linetype = "dashed") +
  annotate("label", y = 1, size = 5,
           label = paste("Max ranked\nCutpoint =",
                         signif(LUAD.Opt_cut$cutpoint[1,1],3)),
           x =  LUAD.Opt_cut$cutpoint[1,1],
           label.size = NA)

Fig_6A

Figure 6B

显示TCGA-LUAD中25样品中肌成纤维细胞丰度测量值分布的密度图

Fig_6B <- 
  Opt_cut.plot$Myo_Fibs.pct$distribution +
  theme_pubr(base_size = 15) +
  scale_fill_manual(values =  c(pal_d3()(3)[3], "grey70"))+
  scale_colour_manual(values =  c(pal_d3()(3)[3], "grey70")) +
  theme(legend.position = "none", plot.title = element_blank()) +
  xlab("Myo (% of all fibroblasts)")
Fig_6B 

Figure 6C

生存曲线

#生存曲线数据整理
Myo.ggsurvplot_list <- list()
for(i in levels(Merged.LUADtraits$Dataset.factor)){
  cox_res <- summary(coxph(Surv(OS_YEARS, OS) ~ Myo_cat,
                           data = Merged.LUADtraits[Merged.LUADtraits$Sample.Subtype == "LUAD" &
                                                      Merged.LUADtraits$Dataset.factor == i, ]))
  Myo.ggsurvplot_list[[i]] <- ggsurvplot(survfit(Surv(OS_YEARS, OS) ~ Myo_cat,
                                                 data = Merged.LUADtraits[Merged.LUADtraits$Sample.Subtype == "LUAD" &
                                                                            Merged.LUADtraits$Dataset.factor == i, ]),
                                         main = i,
                                         pval = signif(cox_res$logtest[3],3), pval.coord = c(0,0.15), pval.size = 10,
                                         conf.int = F, 
                                         risk.table = TRUE, risk.table.title = element_blank(),
                                         risk.table.fontsize = 5,tables.height = 0.4,
                                         risk.table.pos = "in",
                                         font.tickslab = 17,
                                         censor.shape = 124, censor.size = 4,
                                         legend.labs = c("Low", "High"),
                                         palette = c("grey70",pal_d3()(3)[3]),
                                         ggtheme = theme_pubr(base_size = 17),
                                         risk.table.y.col = T,
                                         legend = c(0.85,0.85), legend.title = element_blank())
  Myo.ggsurvplot_list[[i]]$plot <- Myo.ggsurvplot_list[[i]]$plot + 
    ggtitle(i) + ylab("Overall Survival Rate") +
    theme(axis.text.x = element_blank(), axis.title.x =  element_blank(),
          plot.margin = margin(b= 0)) +
    coord_cartesian(xlim = c(-0.3,NA))
  
  Myo.ggsurvplot_list[[i]]$table <- Myo.ggsurvplot_list[[i]]$table + xlab("Time (Years)")+
    theme(plot.margin = margin(t = 0), plot.background = element_blank()) +
    coord_cartesian(xlim = c(-0.3,NA))
}
#Kaplan-Meier图显示TCGA-LUAD中25队列患者生存率,按肌成纤维细胞丰度使用二分最佳临界点分层
Fig_6C <- 
  ggarrange(Myo.ggsurvplot_list$TCGA$plot + theme(plot.title = element_blank()),
            Myo.ggsurvplot_list$TCGA$table,
            ncol = 1, nrow = 2, heights = c(0.7,0.3), align = "v",
            common.legend = T, legend = "none")


Fig_6C

Figure 6D

将该阈值应用于三个验证队列,证明了一致的显著患者分层

Fig_6D <- 
  ggarrange(Myo.ggsurvplot_list$GSE72094$plot, Myo.ggsurvplot_list$GSE31210$plot,
           Myo.ggsurvplot_list$GSE68465$plot,
          Myo.ggsurvplot_list$GSE72094$table,  Myo.ggsurvplot_list$GSE31210$table,
           Myo.ggsurvplot_list$GSE68465$table,
          ncol = 3, nrow = 2, heights = c(0.7,0.3), align = "v",
          common.legend = T, legend = "none")

Fig_6D

Figure 6E

散点图显示使用不同截点按肺泡成纤维细胞丰度对TCGA中LUAD25样本进行二分的标准化对数秩生存统计量(与总生存率的相关性)的变化。

Fig_6E <- 
  data.frame(stat = LUAD.Opt_cut$Alveolar_Fibs.pct$stats,
             x = LUAD.Opt_cut$Alveolar_Fibs.pct$cuts) %>%
  ggplot(aes(x = x, y = stat, colour = x < LUAD.Opt_cut$cutpoint[2,1])) +
  geom_point(size = 2) +
  theme_pubr(base_size = 15) +
  theme(legend.position = "none") +
  scale_colour_manual(values =  c(pal_d3()(3)[1], "grey70")) +
  xlab("Alveolar (% of all fibroblasts)") +
  ylab("Standardized Log-Rank statistic") +
  geom_vline(xintercept = LUAD.Opt_cut$cutpoint[2,1], 
             linetype = "dashed") +
  annotate("label", y = 1, size = 5,
           label = paste("Max ranked\nCutpoint =",
                         signif(LUAD.Opt_cut$cutpoint[2,1],3)),
           x =  LUAD.Opt_cut$cutpoint[2,1],
           label.size = NA)

Fig_6E

Figure 6F

显示TCGA-LUAD中25样品中肺泡成纤维细胞丰度测量值分布的密度图

Fig_6F <- 
  Opt_cut.plot$Alveolar_Fibs.pct$distribution +
  theme_pubr(base_size = 15) +
  scale_fill_manual(values =  c(pal_d3()(3)[1], "grey70"))+
  scale_colour_manual(values =  c(pal_d3()(3)[1], "grey70")) +
  theme(legend.position = "none", plot.title = element_blank()) +
  xlab("Alveolar (% of all fibroblasts)")

Fig_6F

Figure 6G

代码如下

#生存曲线数据整理
Alv.ggsurvplot_list <- list()
for(i in levels(Merged.LUADtraits$Dataset.factor)){
  cox_res <- summary(coxph(Surv(OS_YEARS, OS) ~ Alv_cat,
                           data = Merged.LUADtraits[Merged.LUADtraits$Sample.Subtype == "LUAD" &
                                                      Merged.LUADtraits$Dataset.factor == i, ]))
  Alv.ggsurvplot_list[[i]] <- ggsurvplot(survfit(Surv(OS_YEARS, OS) ~ Alv_cat,
                                                 data = Merged.LUADtraits[Merged.LUADtraits$Sample.Subtype == "LUAD" &
                                                                            Merged.LUADtraits$Dataset.factor == i, ]),
                                         main = i,
                                         pval = T, pval.coord = c(0,0.15), pval.size = 10,
                                         conf.int = F, 
                                         risk.table = TRUE, risk.table.title = element_blank(),
                                         risk.table.fontsize = 5,tables.height = 0.4,
                                         risk.table.pos = "in",
                                         font.tickslab = 17,
                                         censor.shape = 124, censor.size = 4,
                                         legend.labs = c("Low", "High"),
                                         palette = c("grey70",pal_d3()(3)[1]),
                                         ggtheme = theme_pubr(base_size = 17),
                                         risk.table.y.col = T,
                                         legend = c(0.85,0.85), legend.title = element_blank())
  Alv.ggsurvplot_list[[i]]$plot <- Alv.ggsurvplot_list[[i]]$plot + 
    ggtitle(i) + ylab("Overall Survival Rate") +
    theme(axis.text.x = element_blank(), axis.title.x =  element_blank(),
      plot.margin = margin(b= 0)) +
    coord_cartesian(xlim = c(-0.3,NA))
  
  Alv.ggsurvplot_list[[i]]$table <- Alv.ggsurvplot_list[[i]]$table + xlab("Time (Years)")+
    theme(plot.margin = margin(t = 0), plot.background = element_blank()) +
    coord_cartesian(xlim = c(-0.3,NA))
}
#Kaplan-Meier图显示TCGA-LUAD中25队列患者生存率,按肺泡成纤维细胞丰度使用二分最佳临界点分层
Fig_6G <- 
  ggarrange(Alv.ggsurvplot_list$TCGA$plot + theme(plot.title = element_blank()),
            Alv.ggsurvplot_list$TCGA$table,
            ncol = 1, nrow = 2, heights = c(0.7,0.3), align = "v",
            common.legend = T, legend = "none")

Fig_6G

Figure 6H

将该阈值应用于三个验证队列,证明了一致的显著患者分层

Fig_6H <- 
  ggarrange(Alv.ggsurvplot_list$GSE72094$plot, Alv.ggsurvplot_list$GSE31210$plot,
          Alv.ggsurvplot_list$GSE68465$plot,
          Alv.ggsurvplot_list$GSE72094$table,  Alv.ggsurvplot_list$GSE31210$table,
          Alv.ggsurvplot_list$GSE68465$table,
          ncol = 3, nrow = 2, heights = c(0.7,0.3), align = "v",
          common.legend = T, legend = "none")

Fig_6H

多变量分类分析

森林图显示协变量独立风险比和上述分析的所有LUAD患者队列的4年总生存率多变量COX回归分析的校正p值,使用肌成纤维细胞丰度、疾病分期和患者年龄作为自变量

max_year <- 4
surv_data <- Merged.LUADtraits[Merged.LUADtraits$Dataset.factor %in% c("TCGA", "GSE68465", "GSE72094", "GSE31210"), ]
split_data <- survSplit(Surv(OS_YEARS, OS) ~ 
                          Myo_cat + Stage_4cat + Age,
                        zero = -0.1,
                        data = surv_data[!is.na(surv_data$OS_YEARS), ],
                        cut = max_year, episode = "tgroup", id = "id")
names(split_data)[1:2] <- c("Myo", "Stage")
Myo_cox.res <- coxph(Surv(OS_YEARS, OS) ~ 
                       Myo + Stage + Age,
                     data =  split_data[split_data$tgroup == 1,])
summary(Myo_cox.res)$coefficients

Figure 6I

代码如下

Fig_6I <- forest_model(Myo_cox.res,
                       format_options = forest_model_format_options(
                         banded = T, text_size = 5, point_size = 4)) 

Fig_6I

Figure 6J

森林图显示协变量独立风险比(±95%置信区间)和上述所有LUAD患者队列4年总生存率的多变量考克斯回归分析的校正p值,使用肺泡成纤维细胞丰度、疾病分期和患者年龄作为自变量


surv_data <- Merged.LUADtraits[Merged.LUADtraits$Dataset.factor %in% c("TCGA", "GSE68465", "GSE72094", "GSE31210"), ]
split_data <- survSplit(Surv(OS_YEARS, OS) ~ 
                          Alv_cat + Stage_4cat + Age,
                        zero = -0.1,
                        data = surv_data[!is.na(surv_data$OS_YEARS), ],
                        cut = max_year, episode = "tgroup", id = "id")
names(split_data)[1:2] <- c("Alveolar", "Stage")
Alv_cox.res <- coxph(Surv(OS_YEARS, OS) ~ 
                       Alveolar + Stage + Age,
                     data =  split_data[split_data$tgroup == 1,])

Fig_6J <- forest_model(Alv_cox.res,
                    format_options = forest_model_format_options(
                      banded = T, text_size = 2.5, point_size = 1)) 

Fig_6J

Figure 6

多个LUAD数据集中,肺泡和表皮成纤维细胞丰度与更好的总生存率相关。这种关联在肺泡成纤维细胞中尤其一致,在所有分析的数据集中都是显著的。因此,采用与上述相同的方法来测试将肺泡成纤维细胞丰度作为预后标志物的可能性。同样,这表明将LUAD队列分为肺泡成纤维细胞高或低在分层总生存率方面始终有效,并且这种关联与疾病分期和患者年龄无关。

相关推荐

悠悠万事,吃饭为大(悠悠万事吃饭为大,什么意思)

新媒体编辑:杜岷赵蕾初审:程秀娟审核:汤小俊审签:周星...

高铁扒门事件升级版!婚宴上‘冲喜’老人团:我们抢的是社会资源

凌晨两点改方案时,突然收到婚庆团队发来的视频——胶东某酒店宴会厅,三个穿大红棉袄的中年妇女跟敢死队似的往前冲,眼瞅着就要扑到新娘的高额钻石项链上。要不是门口小伙及时阻拦,这婚礼造型团队熬了三个月的方案...

微服务架构实战:商家管理后台与sso设计,SSO客户端设计

SSO客户端设计下面通过模块merchant-security对SSO客户端安全认证部分的实现进行封装,以便各个接入SSO的客户端应用进行引用。安全认证的项目管理配置SSO客户端安全认证的项目管理使...

还在为 Spring Boot 配置类加载机制困惑?一文为你彻底解惑

在当今微服务架构盛行、项目复杂度不断攀升的开发环境下,SpringBoot作为Java后端开发的主流框架,无疑是我们手中的得力武器。然而,当我们在享受其自动配置带来的便捷时,是否曾被配置类加载...

Seata源码—6.Seata AT模式的数据源代理二

大纲1.Seata的Resource资源接口源码2.Seata数据源连接池代理的实现源码3.Client向Server发起注册RM的源码4.Client向Server注册RM时的交互源码5.数据源连接...

30分钟了解K8S(30分钟了解微积分)

微服务演进方向o面向分布式设计(Distribution):容器、微服务、API驱动的开发;o面向配置设计(Configuration):一个镜像,多个环境配置;o面向韧性设计(Resista...

SpringBoot条件化配置(@Conditional)全面解析与实战指南

一、条件化配置基础概念1.1什么是条件化配置条件化配置是Spring框架提供的一种基于特定条件来决定是否注册Bean或加载配置的机制。在SpringBoot中,这一机制通过@Conditional...

一招解决所有依赖冲突(克服依赖)

背景介绍最近遇到了这样一个问题,我们有一个jar包common-tool,作为基础工具包,被各个项目在引用。突然某一天发现日志很多报错。一看是NoSuchMethodError,意思是Dis...

你读过Mybatis的源码?说说它用到了几种设计模式

学习设计模式时,很多人都有类似的困扰——明明概念背得滚瓜烂熟,一到写代码就完全想不起来怎么用。就像学了一堆游泳技巧,却从没下过水实践,很难真正掌握。其实理解一个知识点,就像看立体模型,单角度观察总...

golang对接阿里云私有Bucket上传图片、授权访问图片

1、为什么要设置私有bucket公共读写:互联网上任何用户都可以对该Bucket内的文件进行访问,并且向该Bucket写入数据。这有可能造成您数据的外泄以及费用激增,若被人恶意写入违法信息还可...

spring中的资源的加载(spring加载原理)

最近在网上看到有人问@ContextConfiguration("classpath:/bean.xml")中除了classpath这种还有其他的写法么,看他的意思是想从本地文件...

Android资源使用(android资源文件)

Android资源管理机制在Android的开发中,需要使用到各式各样的资源,这些资源往往是一些静态资源,比如位图,颜色,布局定义,用户界面使用到的字符串,动画等。这些资源统统放在项目的res/独立子...

如何深度理解mybatis?(如何深度理解康乐服务质量管理的5个维度)

深度自定义mybatis回顾mybatis的操作的核心步骤编写核心类SqlSessionFacotryBuild进行解析配置文件深度分析解析SqlSessionFacotryBuild干的核心工作编写...

@Autowired与@Resource原理知识点详解

springIOCAOP的不多做赘述了,说下IOC:SpringIOC解决的是对象管理和对象依赖的问题,IOC容器可以理解为一个对象工厂,我们都把该对象交给工厂,工厂管理这些对象的创建以及依赖关系...

java的redis连接工具篇(java redis client)

在Java里,有不少用于连接Redis的工具,下面为你介绍一些主流的工具及其特点:JedisJedis是Redis官方推荐的Java连接工具,它提供了全面的Redis命令支持,且...