简介

ROC（Receiver Operating Characteristic）曲线分析是一种广泛用于分类模型性能评估的工具，尤其适用于处理二分类问题。在基因研究中，ROC曲线通常用于评估不同基因在预测疾病状态、临床分型或其他生物学现象中的敏感性和特异性。它通过比较基因的表达水平与实际的类别标签（如健康与疾病）之间的关系，评估该基因在区分不同类别的能力。

AUC（Area Under the Curve）是ROC曲线下方的面积值，它可以作为模型性能的一个量化指标。AUC值越接近1，表示模型的预测能力越强。AUC = 1：表示完美分类器。AUC = 0.5：表示分类器没有任何区分能力，相当于随机预测。AUC < 0.5：表示分类器性能较差，预测结果不可靠。

今天生信分析实习生来分享如何使用R语言进行多基因的ROC 分析，并计算AUC值绘制ROC曲线，后台私信发送“转录组分析五”可领取示例数据与源代码。

step1 加载必要的包和读取数据

数据主要有三个文件：

“deseq2_result_diff_count.csv”文件为转录组分析一里面获取的差异基本丰度数据，该节内容也包含了“group.csv”文件。

“string_node_degrees.tsv”文件为转录组分析四获得的基因互作强度文件

## 加载必要的包
library(ggplot2)
library(tidyverse)
library(pROC)

# 读取degree数据和所有差异基因的表达数据
degree_data <- read.csv('string_node_degrees.tsv',comment.char = '',header = T,sep = '\t',check.names = F)
exp_data <- read.csv('deseq2_result_diff_count.csv',header = T,row.names = 1)
group <- read.csv('group.csv',header = T)

step2 提取基因丰度合并数据

根据基因的degree强度提取5个Hub基因名，然后提取这些基因的丰度数据绘制ROC曲线。

# 按照node_degree降序排序，并选取前5个节点,然后提取对应的基因表达数据
top5_nodes <- degree_data %>%
  arrange(desc(node_degree)) %>%
  head(5) %>%
  pull(`#node`)

top5_nodes <- as.data.frame(t(exp_data[top5_nodes,]))

## 将基因表达数据的样本顺序排序设置为group里面的样本顺序，然后合并分组信息
top5_nodes <- top5_nodes[group$sample,]
rownames(top5_nodes)==group$sample
top5_nodes$Group <- group$group

step3 将分组信息转为二分类变量

# 将Group列转换为二元变量，Control = 0, 其他组 = 1
top5_nodes$Group_binary <- ifelse(top5_nodes$Group == "Control", 0, 1)

step4 进行多物种ROC分析，并提取数据

# 初始化一个空数据框，用于存储所有基因的ROC结果
roc_results <- data.frame(Gene = character(),
                          Sensitivity = numeric(),
                          Specificity = numeric(),
                          AUC = numeric(),
                          stringsAsFactors = FALSE)

# 遍历每个基因进行ROC分析
for (gene in names(top5_nodes)[1:(ncol(top5_nodes)-2)]) { # 遍历基因列
  # 计算ROC
  roc_obj <- roc(top5_nodes$Group_binary, top5_nodes[[gene]], smooth=T)

  # 提取每个基因的多行数据
  temp_data <- data.frame(
    Gene = rep(gene, length(roc_obj$sensitivities)), # 每行对应同一个基因
    Sensitivity = roc_obj$sensitivities,            # 敏感性
    Specificity = roc_obj$specificities,            # 特异性
    AUC = rep(auc(roc_obj), length(roc_obj$sensitivities)) # 重复AUC值
  )

  # 合并到总数据框
  roc_results <- rbind(roc_results, temp_data)
}

step5 绘图

# 修改Gene列以包含AUC值
roc_results$Gene <- paste0(roc_results$Gene, " (AUC=", round(roc_results$AUC, 2), ")")

## 绘制多基因ROC曲线

ROC_plot <- ggplot(roc_results, aes(x = 1 - Specificity, y = Sensitivity, color = Gene)) +
  geom_line(size=1.5) + # 绘制基因ROC曲线
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray",size=1.5) + # 添加对角线
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))+
  theme_bw() +
  labs(
    title = "ROC Curves with AUC for Hub Genes",
    x = "1 - Specificity",
    y = "Sensitivity",
    color = "Gene (AUC)"
  ) +
  theme(legend.position = "right",
        axis.title = element_text(size = 25,face = 'bold'),
        axis.text = element_text(size = 20,color = 'black'),
        legend.title = element_text(size = 25,face = 'bold'),
        axis.ticks = element_blank(),
        legend.text = element_text(size = 20),
        plot.margin = unit(c(4,4,4,4),'cm'),
        plot.title = element_text(hjust = 0.5,size = 30,face = 'bold'))

ggsave(plot = ROC_plot,'ROC_plot.png',width = 12,height = 8,dpi = 300)
ggsave(plot = ROC_plot,'ROC_plot.pdf',width = 12,height = 8)