#install.packages("survminer")
rm(list = ls())
#设置工作目录
setwd("D:/※raw.data/CZM.exo.m6A/2022.2.22")
library(data.table)
#读入文件
gene_exp<-fread("TCGA-PAAD.htseq_fpkm.All5W.txt",data.table = F)
#转为数据框
gene_exp<-as.data.frame(gene_exp)
row.names(gene_exp)<-gene_exp$Ensembl_ID
gene_exp=gene_exp[,-1]
#得到基因在各个样本中的表达数据
A5=t(gene_exp[which(row.names(gene_exp)=="A5"),])
H1=t(gene_exp[which(row.names(gene_exp)=="H1"),])
library(dplyr)
GENE=cbind(ALKBH5,HES1) %>% as.data.frame()
##下面取获得临床数据
clinical<-read.table(file="TCGA-PAAD.survival.tsv",sep = "\t",header = T)
colnames(clinical)[1]="SampleID"
clinical=clinical[,-3]#删除重复患者ID
clinical$time<- clinical$time/30
colnames(clinical)[3]<-"futime.Months"
colnames(clinical)[2]<-"fustat"
#合并数据，得到的数据既有临床数据，又有某个基因的表达数据
svdata=data.frame(clinical[match(rownames(GENE),clinical$SampleID),],GENE)
library(limma)
#对重复样本名取平均表达量
if(sum(duplicated(svdata$SampleID))>0) #判断是否有重复样本
? svdata<-avereps(svdata,ID=svdata$SampleID) %>%? as.data.frame()
svdata<- na.omit(svdata)
row.names(svdata) <- svdata$SampleID
name<-svdata$SampleID
svdata<-svdata[,-1]
#svdata$SampleID[duplicated(svdata$SampleID)==TRUE]
#转为数值型
svdata<-as.data.frame(lapply(svdata,as.numeric))
row.names(svdata) <- name
##目前的数据后，就可以进行生存分析作图了
library(survival)
library(survminer)#找best separation用的是survminer的函数
##对数据集的基因进行bestSeparation统计
res.cut <- surv_cutpoint(svdata, time = "futime.Months",
? ? ? ? ? ? ? ? ? ? ? ? event = "fustat",
? ? ? ? ? ? ? ? ? ? ? ? variables = names(svdata)[3:ncol(svdata)],
? ? ? ? ? ? ? ? ? ? ? ? minprop = 0.3) #默认组内sample不能低于30%
##按照bestSeparation分高低表达
res.cat <- surv_categorize(res.cut)
##统计作图
my.surv <- Surv(res.cat$futime, res.cat$fustat)
pl<-list()
for (i in colnames(res.cat)[3:ncol(svdata)]) {
? group <- res.cat[,i]
? survival_dat <- data.frame(group = group)
? fit <- survfit(my.surv ~ group)
? ##计算HR以及95%CI
? ##修改分组参照
? group <- factor(group, levels = c("low", "high"))
? data.survdiff <- survdiff(my.surv ~ group)
? p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
? HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
? up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
? low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
? #只画出p value<=0.05的基因，如果不想筛选，就删掉下面这行
? #if (p.val>0.05) next
? HR <- paste("Hazard Ratio = ", round(HR,2), sep = "")
? CI <- paste("95% CI: ", paste(round(low95,2), round(up95,2), sep = " - "), sep = "")
? #按照基因表达量从低到高排序，便于取出分界表达量
? svsort <- svdata[order(svdata[,i]),]
? pl[[i]]<-ggsurvplot(fit, data = survival_dat ,
? ? ? ? ? ? ? ? ? ? ? #ggtheme = theme_bw(), #想要网格就运行这行
? ? ? ? ? ? ? ? ? ? ? conf.int = F, #不画置信区间，想画置信区间就把F改成T
? ? ? ? ? ? ? ? ? ? ? #conf.int.style = "step",#置信区间的类型，还可改为ribbon
? ? ? ? ? ? ? ? ? ? ? censor = F, #不显示观察值所在的位置
? ? ? ? ? ? ? ? ? ? ? #palette = c("#D95F02","#1B9E77"), #线的颜色对应高、低
? ? ? ? ? ? ? ? ? ? ? palette = c("red","#008EA0FF"),#线的颜色
? ? ? ? ? ? ? ? ? ? ? legend.title = i,#基因名写在图例题目的位置
? ? ? ? ? ? ? ? ? ? ? font.legend = 11,#图例的字体大小
? ? ? ? ? ? ? ? ? ? ? #font.title = 12,font.x = 10,font.y = 10,#设置其他字体大小
? ? ? ? ? ? ? ? ? ? ? #在图例上标出高低分界点的表达量，和组内sample数量
? ? ? ? ? ? ? ? ? ? ? legend.labs=c(paste0(">",round(svsort[fit$n[2],i],2),"(",fit$n[1],")"),
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? paste0("<",round(svsort[fit$n[2],i],2),"(",fit$n[2],")")),
? ? ? ? ? ? ? ? ? ? ? xlab="Months",#x轴标为Months
? ? ? ? ? ? ? ? ? ? ? #在左下角标出pvalue、HR、95% CI
? ? ? ? ? ? ? ? ? ? ? #太小的p value标为p < 0.001
? ? ? ? ? ? ? ? ? ? ? pval = paste(pval = ifelse(p.val < 0.001, "p < 0.001",
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? paste("p = ",round(p.val,3), sep = "")),
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? HR, CI, sep = "\n"))
? #如果想要一个图保存为一个pdf文件，就把下面这行前面的“#”删掉
? ggsave(paste0(i,".pdf"),width = 4,height = 4)
}
length(pl)
#########################批量出图
#用survminer包自带的函数组图
res <- arrange_ggsurvplots(pl,
? ? ? ? ? ? ? ? ? ? ? ? ? print = T,
? ? ? ? ? ? ? ? ? ? ? ? ? ncol = 1, nrow = 2)#每页纸画几列几行
#保存到pdf文件
ggsave("bestSurvPlot.pdf",res,width = 12,height = 16)
