批量生存分析画图 循环画图 必须用pdf 不可以使用jpeg

批量生存分析画图 循环画图 必须用pdf 不可以使用jpeg

load(file ="G:/r/duqiang_IPF/GSE70866—true—_BAL_IPF_donors_RNA-seq/Rdatafor_freibrug.RData")
dim(expr.freiburg)
dim(expr.freiburg_clean)
head(expr.freiburg_clean)
head(meta.14550)
dim(meta.14550)
dim(meta.freiburg)
head(meta.freiburg)
head(meta.14550)

expr_ipf_freiburg=expr.freiburg_clean[,
        which(colnames(expr.freiburg_clean)=='GSM1820739'):length(colnames(expr.freiburg_clean))]
dim(expr_ipf_freiburg)
head(expr_ipf_freiburg)[,1:5]


meta_ipf_freibrug=meta.14550[which(rownames(meta.14550)=='GSM1820739'):which(rownames(meta.14550)=='GSM1820800'),
                                ]
head(meta_ipf_freibrug)[,1:3]
colnames(meta_ipf_freibrug)
colnames(meta_ipf_freibrug)=c('event','time','Sex',
                              'Diagnosis','Age','GEO_acc')
head(meta_ipf_freibrug)
dim(meta_ipf_freibrug)
str(meta_ipf_freibrug)

library(dplyr)
phe=transform(meta_ipf_freibrug,time=as.numeric(time)) %>% transform(event=as.numeric(event))

head(phe)
#save(meta_ipf_freibrug,expr_ipf_freiburg,file = "G:/r/duqiang_IPF/GSE70866—true—_BAL_IPF_donors_RNA-seq/survival_for_genes-2/ipf_surval.RData")
load("G:/r/duqiang_IPF/GSE70866—true—_BAL_IPF_donors_RNA-seq/survival_for_genes-2/ipf_surval.RData")



#生存分析
library(survival)
library(survminer)
# 利用ggsurvplot快速绘制漂亮的生存曲线图 根据性别
sfit <- survfit(data=phe,Surv(time, event)~Sex)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF"),
           risk.table =TRUE,pval =TRUE,
           conf.int =TRUE,xlab ="Time in months", 
           ggtheme =theme_light(), 
           ncensor.plot = TRUE)

#挑选感兴趣的基因做差异分析
phe$CBX4=ifelse(exprSet['CBX4',]>median(exprSet['CBX4',]),'high','low')
head(phe)
table(phe$CBX4)
ggsurvplot(survfit(Surv(time, event)~CBX4, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~phe[,'CBX4'], data=phe), conf.int=F, pval=TRUE)




exprSet=expr_ipf_freiburg %>% transform(as.numeric()) %>% as.matrix()
head(exprSet)[,1:3]
head(str(exprSet))
#挑选感兴趣的基因做差异分析
gene_interested=readClipboard()
head(gene_interested)
library(stringr)
gene_interested=str_split(pattern = ",",gene_interested)[[1]]
gene_interested=gene_interested[-which(gene_interested=="RAB40A")]
gene_interested

gene_interested %in% rownames(exprSet)
table(gene_interested %in% rownames(exprSet))
gene_interested[!(gene_interested %in% rownames(exprSet))]
getwd()
dir.create("G:/r/duqiang_IPF/GSE70866—true—_BAL_IPF_donors_RNA-seq/survival_for_genes")
setwd("G:/r/duqiang_IPF/GSE70866—true—_BAL_IPF_donors_RNA-seq/survival_for_genes")

table(phe$event)

## 批量生存分析 使用  logrank test 方法
#批量输出  成功！！！推荐这个方法
if(1=1){
  mySurv=with(phe,Surv(time, event))
  for (eachgene in gene_interested) {
    phe$group=ifelse(exprSet[eachgene,]>median(exprSet[eachgene,]),'high','low')
    
    p=ggsurvplot(survfit(mySurv~group, data=phe), conf.int=F, pval=TRUE)
    pdf(paste0(eachgene, "_surv.pdf"),width = 5, height = 5)
    print(p, newpage = FALSE)
    dev.off()
  } 
}

head(phe)
getwd()
#批量输出  成功！ 不推荐这个方法
if(1==2){
  dir.create("G:/r/duqiang_IPF/GSE70866—true—_BAL_IPF_donors_RNA-seq/survival_for_genes-1")
  setwd("G:/r/duqiang_IPF/GSE70866—true—_BAL_IPF_donors_RNA-seq/survival_for_genes-1")
  mySurv=with(phe,Surv(time, event))
  for (eachgene in gene_interested) {
    phe[eachgene]=ifelse(exprSet[eachgene,]>median(exprSet[eachgene,]),'high','low')
    
    p=ggsurvplot(survfit(mySurv~phe[,eachgene], data=phe), conf.int=F, pval=TRUE)
    
    pdf(paste0(eachgene, "_surv.pdf"),width = 5, height = 5)
    print(p, newpage = FALSE)
    dev.off()
  } 
  
}



#批量输出  成功！ 不推荐这个方法 出错了
if(1==2){
  dir.create("G:/r/duqiang_IPF/GSE70866—true—_BAL_IPF_donors_RNA-seq/survival_for_genes-2")
  setwd("G:/r/duqiang_IPF/GSE70866—true—_BAL_IPF_donors_RNA-seq/survival_for_genes-2")
 
   mySurv=with(phe,Surv(time, event))
  for (eachgene in gene_interested) {
    
    phe[eachgene]=as.character(ifelse(exprSet[eachgene,]>median(exprSet[eachgene,]),'high','low'))
    
    p=ggsurvplot(survfit(mySurv~phe[,eachgene], data=phe), conf.int=F, pval=TRUE)
    
    pdf(paste0(eachgene, "_surv.pdf"),width = 5, height = 5)
    print(p, newpage = FALSE)
    dev.off()
  } 
  
}




data.frame("a"=c(1,3))
cbind(c(1,3),c(1,5))
rbind(c(1,3),c(1,5))



## 批量生存分析 使用  logrank test 方法



mySurv=with(phe,Surv(time, event))
log_rank_p <- apply(exprSet , 1 , function(gene){
  #gene=exprSet[1,]
  phe$group=ifelse(gene>median(gene),'high','low')  
  data.survdiff=survdiff(mySurv~group,data=phe)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  return(p.val)
})

log_rank_p=sort(log_rank_p)
head(log_rank_p)
boxplot(log_rank_p) 

phe$H6PD=ifelse(exprSet['H6PD',]>median(exprSet['H6PD',]),'high','low')
table(phe$H6PD)
ggsurvplot(survfit(Surv(time, event)~H6PD, data=phe), conf.int=F, pval=TRUE)

# 前面如果我们使用了WGCNA找到了跟grade相关的基因模块，然后在这里就可以跟生存分析的显著性基因做交集
# 这样就可以得到既能跟grade相关，又有临床预后意义的基因啦。