step08-gene-survival
Chard Liu
2019年7月25日
载入hubgene数据
Sys.setlocale('LC_ALL','C')## [1] "C"require(tidyverse)## Loading required package: tidyverse## -- Attaching packages ------------------------------------------- tidyverse 1.2.1 --## <U+221A> ggplot2 3.2.0 <U+221A> purrr 0.3.2
## <U+221A> tibble 2.1.3 <U+221A> dplyr 0.8.3
## <U+221A> tidyr 0.8.3 <U+221A> stringr 1.4.0
## <U+221A> readr 1.3.1 <U+221A> forcats 0.4.0## Warning: package 'dplyr' was built under R version 3.6.1## -- Conflicts ---------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()load(file = "F:/Bioinfor_project/Breast/AS_research/AS/result/hubgene.Rdata")
head(data)## # A tibble: 6 x 22
## ID group CCL14 HBA1 CCL16 TUBB3 PAM50 Os_time OS_event RFS_time
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <fct> <int> <int> <int>
## 1 TCGA~ TNBC 3.37 0 0 5.94 Basal 967 1 NA
## 2 TCGA~ TNBC 4.97 3.39 0 6.00 Basal 584 0 NA
## 3 TCGA~ TNBC 4.77 0 0 5.04 Basal 2654 0 NA
## 4 TCGA~ TNBC 4.19 3.74 2.84 5.13 Basal 754 1 NA
## 5 TCGA~ TNBC 0 0 0 5.65 Basal 2048 0 2048
## 6 TCGA~ TNBC 5.34 0 4.07 6.04 Basal 1027 0 1027
## # ... with 12 more variables: RFS_event <int>, age <int>, ER <fct>,
## # PR <fct>, gender <fct>, HER2 <fct>, Margin_status <fct>, Node <int>,
## # M_stage <fct>, N_stage <fct>, T_stage <fct>, `Pathologic stage` <fct>dim(data)## [1] 228 22## 筛选出TNBC
mydata<-data %>%
rename(OS_time=Os_time) %>%
filter(group=="TNBC") 批量绘制生存曲线
library(survival)
library(survminer)## Loading required package: ggpubr## Loading required package: magrittr##
## Attaching package: 'magrittr'## The following object is masked from 'package:purrr':
##
## set_names## The following object is masked from 'package:tidyr':
##
## extractmysurv<-Surv(mydata$OS_time,mydata$OS_event)
p<-list()
## 索引因为这里的基因是3:6，可调整格式固定下来
for(i in colnames(mydata)[3:6]){
## 分组
gene<-as.numeric(as.matrix(mydata[,i]))
grouplist<-factor(ifelse(gene>median(gene),"high","low"))
fit <- survfit(mysurv ~ grouplist, data = mydata)
## 分组比较
data.survdiff <- survdiff(mysurv ~ grouplist)
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]))
## 计算HR CI
HR <- paste("Hazard Ratio = ", round(HR,2), sep = "")
CI <- paste("95% CI: ", paste(round(low95,2), round(up95,2), sep = " - "), sep = "")
### 批量绘图
p[[i]]<-ggsurvplot(fit, data = mydata ,
#ggtheme = theme_bw(), #想要网格就运行这行
conf.int = F, #不画置信区间，想画置信区间就把F改成T
#conf.int.style = "step",#置信区间的类型，还可改为ribbon
censor = F, #不显示观察值所在的位置
palette = c("#00468BFF","#ED0000FF"), #线的颜色对应高、低
#pal_lancet("lanonc")(9)
#[1] "#00468BFF" "#ED0000FF" "#42B540FF" "#0099B4FF" "#925E9FFF" "#FDAF91FF"
#[7] "#AD002AFF" "#ADB6B6FF" "#1B1919FF"
legend.title = i,#基因名写在图例题目的位置
font.legend = 11,#图例的字体大小
#font.title = 12,font.x = 10,font.y = 10,#设置其他字体大小
legend.labs = c("high", "low"),
#在左下角标出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)
}
### 组图并保存
res <- arrange_ggsurvplots(p,
print = T,
ncol = 2,nrow = 2)#每页纸画几列几行
#保存到pdf文件
ggsave("SurvPlot.pdf",res,width = 12,height = 8)包装为批量绘制基因生存曲线的函数
可将数据清洗调整好,为OS_time, OS_event tibble格式第一列为TCGA id, 第三列开始为gene-到之后的gene名
head(mydata)## # A tibble: 6 x 22
## ID group CCL14 HBA1 CCL16 TUBB3 PAM50 OS_time OS_event RFS_time
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <fct> <int> <int> <int>
## 1 TCGA~ TNBC 3.37 0 0 5.94 Basal 967 1 NA
## 2 TCGA~ TNBC 4.97 3.39 0 6.00 Basal 584 0 NA
## 3 TCGA~ TNBC 4.77 0 0 5.04 Basal 2654 0 NA
## 4 TCGA~ TNBC 4.19 3.74 2.84 5.13 Basal 754 1 NA
## 5 TCGA~ TNBC 0 0 0 5.65 Basal 2048 0 2048
## 6 TCGA~ TNBC 5.34 0 4.07 6.04 Basal 1027 0 1027
## # ... with 12 more variables: RFS_event <int>, age <int>, ER <fct>,
## # PR <fct>, gender <fct>, HER2 <fct>, Margin_status <fct>, Node <int>,
## # M_stage <fct>, N_stage <fct>, T_stage <fct>, `Pathologic stage` <fct>genesurv<-function(mydata){
library(survival)
library(survminer)
mysurv<-Surv(mydata$OS_time,mydata$OS_event)
p<-list()
## 索引因为这里的基因是3:6，可调整格式固定下来
for(i in colnames(mydata)[3:6]){
## 分组
gene<-as.numeric(as.matrix(mydata[,i]))
grouplist<-factor(ifelse(gene>median(gene),"high","low"))
fit <- survfit(mysurv ~ grouplist, data = mydata)
## 分组比较
data.survdiff <- survdiff(mysurv ~ grouplist)
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]))
## 计算HR CI
HR <- paste("Hazard Ratio = ", round(HR,2), sep = "")
CI <- paste("95% CI: ", paste(round(low95,2), round(up95,2), sep = " - "), sep = "")
### 批量绘图
p[[i]]<-ggsurvplot(fit, data = mydata ,
#ggtheme = theme_bw(), #想要网格就运行这行
conf.int = F, #不画置信区间，想画置信区间就把F改成T
#conf.int.style = "step",#置信区间的类型，还可改为ribbon
censor = F, #不显示观察值所在的位置
palette = c("#00468BFF","#ED0000FF"), #线的颜色对应高、低
#pal_lancet("lanonc")(9)
#[1] "#00468BFF" "#ED0000FF" "#42B540FF" "#0099B4FF" "#925E9FFF" "#FDAF91FF"
#[7] "#AD002AFF" "#ADB6B6FF" "#1B1919FF"
legend.title = i,#基因名写在图例题目的位置
font.legend = 11,#图例的字体大小
#font.title = 12,font.x = 10,font.y = 10,#设置其他字体大小
legend.labs = c("high", "low"),
#在左下角标出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)
}
### 组图并保存
res <- arrange_ggsurvplots(p,
print = T,
ncol = 2,nrow = 2)#每页纸画几列几行
#保存到pdf文件
#ggsave("SurvPlot.pdf",res,width = 12,height = 8)
}
## 输入数据mydata-mydata需调整好格式
genesurv(mydata)
#ggsave("SurvPlot.pdf",res,width = 12,height = 8)包装可以绘制感兴趣的基因的函数
## 封装函数
## 参数1：genename
## 参数2：mydata，默认为mydata
genename_surv<-function(genename,data=mydata){
library(survival)
library(survminer)
mydata<-data
mysurv<-Surv(mydata$OS_time,mydata$OS_event)
gene<-as.numeric(as.matrix(mydata[,genename]))
grouplist<-factor(ifelse(gene>median(gene),"high","low"))
fit <- survfit(mysurv ~ grouplist, data = mydata)
## 分组比较
data.survdiff <- survdiff(mysurv ~ grouplist)
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]))
## 计算HR CI
HR <- paste("Hazard Ratio = ", round(HR,2), sep = "")
CI <- paste("95% CI: ", paste(round(low95,2), round(up95,2), sep = " - "), sep = "")
### 批量绘图
ggsurvplot(fit, data = mydata ,
#ggtheme = theme_bw(), #想要网格就运行这行
conf.int = F, #不画置信区间，想画置信区间就把F改成T
#conf.int.style = "step",#置信区间的类型，还可改为ribbon
censor = F, #不显示观察值所在的位置
palette = c("#00468BFF","#ED0000FF"), #线的颜色对应高、低
#pal_lancet("lanonc")(9)
#[1] "#00468BFF" "#ED0000FF" "#42B540FF" "#0099B4FF" "#925E9FFF" "#FDAF91FF"
#[7] "#AD002AFF" "#ADB6B6FF" "#1B1919FF"
legend.title = i,#基因名写在图例题目的位置
font.legend = 11,#图例的字体大小
#font.title = 12,font.x = 10,font.y = 10,#设置其他字体大小
legend.labs = c("high", "low"),
#在左下角标出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"))
}
### 输入获取感兴趣的基因
genename_surv("CCL14",mydata)
### 有了绘图函数可以批量绘制了
name<-colnames(mydata)[3:6]
p<-lapply(name,genename_surv)
## ggsurvplots组图
arrange_ggsurvplots(p,
print = T,
ncol = 2,nrow = 2)#每页纸画几列几行