step10-gene best separation

载入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':
## 
##     extract

##对数据集的基因进行bestSeparation统计
res.cut <- surv_cutpoint(mydata, time = "OS_time", 
                         event = "OS_event", 
                         variables = colnames(mydata)[3:6], 
                         minprop = 0.3) #默认组内sample不能低于30%

##按照bestSeparation分高低表达
res.cat <- surv_categorize(res.cut)
svdata<-mydata
##统计作图
my.surv <- Surv(res.cat$OS_time, res.cat$OS_event)
pl<-list()
## 按位置修改基因的数量
for (i in colnames(res.cat)[3:6]) {
  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("#00468BFF","#ED0000FF" ), #线的颜色对应高、低
           
           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],")")),
           
           #在左下角标出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)

## [1] 4

只保留p<=0.05的10个基因，筛掉了p>0.05的8个基因

批量出图

用survminer包自带的函数组图

res <- arrange_ggsurvplots(pl, 
                           print = T,
                           ncol = 2, nrow = 2)#每页纸画几列几行

#保存到pdf文件
ggsave("F:/Bioinfor_project/Breast/AS_research/AS/result/gene_bestSurvPlot.pdf",res,width = 8,height = 8)

包装为函数gene_bestsurv

gene_bestsurv<-function(gene=genes,data=mydata,col=3,row=3){
  ##加载必需的包
library(survival)
library(survminer)

##对数据集的基因进行bestSeparation统计
res.cut <- surv_cutpoint(data, time = "OS_time", 
                         event = "OS_event", 
                         variables = gene, 
                         minprop = 0.3) #默认组内sample不能低于30%

##按照bestSeparation分高低表达
res.cat <- surv_categorize(res.cut)

##统计作图
my.surv <- Surv(res.cat$OS_time, res.cat$OS_event)
pl<-list()
## 按位置修改基因的数量
for (i in gene) {
  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("#00468BFF","#ED0000FF" ), #线的颜色对应高、低
           
           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],")")),
           
           #在左下角标出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)
}
## 返回有几个图
print(length(pl))
res <- arrange_ggsurvplots(pl, 
                           print = T,
                           ncol = col, nrow = row)#每页纸画几列几行

}

出图

## 绘制多个基因
genes=colnames(mydata)[3:6]
gene_bestsurv(gene = genes,col = 2,row=2)

## [1] 4

#保存到pdf文件
#ggsave("F:/Bioinfor_project/Breast/AS_research/AS/result/gene_bestSurvPlot.pdf",res,width = 8,height = 8)

绘制感兴趣的基因

gene_bestsurv(gene="CCL16",col = 1,row =1)

## [1] 1

#ggsave("F:/Bioinfor_project/Breast/AS_research/AS/result/CCL16_bestSurvPlot.pdf",width = 4,height = 4)

step10-gene best separation

Chard Liu

2019年7月31日

载入hubgene数据

批量出图

包装为函数gene_bestsurv

出图

绘制感兴趣的基因