step09-gene_ROC

载入数据

library(plotROC)

## Warning: package 'plotROC' was built under R version 3.6.1

## Loading required package: ggplot2

library(ggplot2)
library(ggsci)
require(survivalROC)

## Loading required package: survivalROC

require(tidyverse)

## Loading required package: tidyverse

## -- Attaching packages ------------------------------------------- tidyverse 1.2.1 --

## √ tibble  2.1.3     √ purrr   0.3.2
## √ tidyr   0.8.3     √ dplyr   0.8.3
## √ readr   1.3.1     √ stringr 1.4.0
## √ tibble  2.1.3     √ 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>

mydata<-data %>% 
  filter(group=="TNBC") %>% 
  mutate(Os_time=Os_time/365)

ROC分析

Gene_roc<-function(genes=c("CCL16","CCL14","HBA1"),data=mydata,t=5,combine=T){
  ## 增加combine函数参数，将combine当成一个基因 
  if (combine==T){
    FML <- as.formula(paste0('OS_event~',paste(genes,collapse = '+')))
    multifit<-glm(FML,family=binomial(), data = data)
    data$combine<-predict.glm(multifit,type = "link")
    genes<-c(genes,"combine")
  }
  ## 将基因放到前面
  if(!require("tidyverse")) install.packages("tidyverse",update = F,ask = F)
  data<-data %>% select(genes,everything())
  ##
  sroc <- lapply(genes, function(i){
  stroc <- survivalROC(Stime =data$Os_time, status =data$OS_event, 
                       marker =data[[i]],
                       predict.time = t, method = "NNE", ## KM法或NNE法
                       span = .25 * 350^(-.2))
  data.frame(TPF = stroc[["TP"]], FPF = stroc[["FP"]], 
             c = stroc[["cut.values"]], 
             gene = rep(i, length(stroc[["FP"]])))
})
sroclong <- do.call(rbind, sroc)
class(sroclong)
head(sroclong)

## 整合到数据框中
sroclong <- do.call(rbind, sroc)
class(sroclong)
head(sroclong)

## 索引找出各个AUC值命名为lab
lab<-c()
stroc<-c()
for (i in 1:length(genes)){
 stroc[[i]] <- survivalROC(Stime =data$Os_time, status =data$OS_event, 
                       marker = data[[i]],##marker
                       predict.time = t, ##时间
                       method = "NNE", ## KM法或NNE法
                       span = .25 * 350^(-.2))
 ## 索引构建lab
 lab[i]<-paste(genes[i],":",round(stroc[[i]][[6]],3),sep = " ")
 
 }

### 绘制多基因图
pROC<-ggplot(sroclong, aes(x = FPF, y = TPF,label=c,color = gene)) + 
  geom_roc(n.cuts=0,labels = FALSE, stat = "identity")  ## n.cuts设置cut点数
pROC+theme_bw()+
    xlab("1-specificity")+ylab("sensitivity")+ #坐标轴标签
    ggtitle(paste(t,"years ROC curve",sep = " "))+
    geom_abline(intercept = 0, slope = 1, linetype = "dashed")+
    theme(legend.position = c(0.8,0.2),legend.background = element_blank())+
    scale_colour_lancet(name="AUC",labels=lab,limits=genes)

}
## 绘制感兴趣的多个基因的五年ROC,并combine
## 感兴趣的基因
genes<-colnames(mydata)[3:6]
Gene_roc(genes = genes,data=mydata,t=8,combine = T)

#ggsave(filename = "F:/Bioinfor_project/Breast/AS_research/AS/result/hubgene_ROC.pdf",width = 5,height = 5)