#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)
