零基础入门转录组数据分析—

零基础入门转录组数据分析——单基因ROC分析

1. ROC分析的基础知识

1.1 ROC分析是什么？
ROC（Receiver Operating Characteristic）分析是一种用于评估判断准确性的统计方法。它通过将灵敏度和特异度结合起来，以图示的方式展示在不同临界值下，是否能正确判断阳性和阴性样本的能力。

1.2 ROC分析的基本原理？
ROC曲线是通过一系列不同的二分类方式（临界值或决定阈），计算诊断结果的真阳性率和假阳性率，以真阳性率（灵敏度）为纵坐标，假阳性率（1-特异度）为横坐标绘制的曲线，用曲线下面积（AUC）作为衡量预测准确性的指标。

1.3 ROC分析的主要用途？

评估诊断性能： 在生信领域，ROC分析用于评估诊断的性能。通过ROC曲线和曲线下面积（AUC），可以量化区分不同分组的准确性，帮助选择出相对重要的基因。
模型评估： 在机器学习和统计建模中，ROC分析用于评估分类模型的性能。通过比较不同模型的ROC曲线和AUC值，可以选择出最优的模型。
阈值选择： ROC分析可以帮助确定分类或诊断的最佳阈值。通过观察ROC曲线上的不同点，并结合实际需求（如灵敏度和特异度的平衡），可以选择出最合适的阈值。

1.4 ROC分析的优点？

综合评估： ROC曲线将灵敏度和特异度结合在一起，以图示的方式展示诊断或分类方法的整体性能，便于直观理解和评估。
不固定阈值： ROC分析不依赖于特定的分类阈值，而是展示了在不同阈值下的性能变化，有利于使用者根据实际需求选择合适的阈值。
标准化评估： AUC值作为ROC曲线下的面积，提供了一个标准化的评估指标，便于不同测试或方法之间的比较。
适用于不平衡数据集： 在正负样本不平衡的数据集中，ROC分析仍然能够提供一个相对稳定的评估结果，因为AUC值不受样本比例的影响。

举个栗子： 现在有8个样本，1个基因，每个样本对应该基因都有表达量。其中4个样本为疾病组，另外4个样本为对照组。4个疾病样本对应该基因的表达量分别为1，1.1，1.2，1.2；4个对照样本对应该基因的表达量分别为0.5，0.6，0.5，0.7这样可以明显看出疾病样本中该基因为高表达，那么用该基因的表达水平进行ROC分析，它的AUC值就是1，可以完美区分疾病样本和对照样本。

综上所述： ROC分析就是评估区分准确性的方法，不仅能评估单个变量的区分准确性，还能评估模型的区分准确性。

在本章节中仅用单个基因为例来展示如何进行ROC分析，对于模型的ROC分析不做展示

注意：ROC曲线对极端值（如极少数的极端高或低预测值）较为敏感。这些极端值可能会扭曲ROC曲线的形状，从而影响判断的准确性。

2. 单基因ROC分析（Rstudio）——代码实操

本项目以TCGA——肺腺癌为例展开分析
物种：人类（Homo sapiens）
R版本：4.2.2
R包：tidyverse, pROC

废话不多说，代码如下：

2. 1 数据处理

设置工作空间：

rm(list = ls()) # 删除工作空间中所有的对象
setwd('/XX/XX/XX') # 设置工作路径
if(!dir.exists('./19_ROC')){dir.create('./19_ROC')
} 
setwd('./19_ROC/')

加载包：

library(pROC)
library(tidyverse)

导入要分析的表达矩阵TrainRawData，并对TrainRawData的列名进行处理（这是因为在读入的时候系统会默认把样本id中的“-”替换成“.”，所以要给替换回去）

TrainRawData <- read.csv("./data_fpkm.csv", row.names = 1, check.names = F)  # 行名为全部基因名，每列为样本名
colnames(TrainRawData) <- gsub('.', '-', colnames(TrainRawData), fixed = T)

TrainRawData如下图所示，行为基因名（symbol），列为样本名
在这里插入图片描述
导入分组信息表TrainGroup

TrainGroup <- read.csv("./data_group.csv", row.names = 1) # 为每个样本的分组信息（tumor和control）
colnames(TrainGroup) <- c('sample', 'group')

TrainGroup 如下图所示，第一列sample为样本名，第二列为样本对应的分组 （分组为二分类变量：disease和control）
在这里插入图片描述
导入要用于分析的基因HubGene (1个基因，这里只用一个基因作为展示)

HubGene <- data.frame(symbol = 'VDAC1')

HubGene 如下图所示，只有一列：1个基因的基因名

在这里插入图片描述

从TrainRawData中取出这1个基因对应的表达矩阵，并且与之前准备的分组信息表TrainGroup进行合并

TrainData <- TrainRawData[HubGene$symbol, ] %>% t() %>% as.data.frame()
TrainData <- merge(TrainData, TrainGroup, by.x = "row.names", by.y = 'sample')
TrainData <- column_to_rownames(TrainData, var = 'Row.names')

TrainData 如下图所示，行为样本名，第一列为基因的表达量，第二列为分组信息（这里是disease和control）。
在这里插入图片描述

2. 2 单基因ROC分析

直接通过roc函数进行单基因ROC分析

response = TrainData$group —— 这个参数指定了真实的目标变量或类别标签，简单说就是想要预测的分组
predictor = TrainData[, 1] —— 这个参数指定了要参与预测的变量，这里是TrainData的第一列，也就是前面提到的VDAC1表达量
levels = c(‘control’, ‘disease’) —— 这个参数指定了response变量中类别的顺序，用于ROC曲线的计算。在这里，它告诉roc函数’control’是负类（或参考类），而’disease’是正类。（注意：这个levels非常重要，设置参考类如果反了，那么相应的结果也会是反的）

roc <- roc(response = TrainData$group, predictor = TrainData[, 1],levels = c('control', 'disease'))
roc$auc

roc$auc如下图所示，表明曲线下面积（Area under the curve, AUC）为0.8867（AUC越接近于1表明预测准确性越高）。

在这里插入图片描述

2. 3 ROC曲线简单可视化

接下来一步就是要对ROC分析结果进行简单可视化，毕竟文字的展示效果不如图片更加直观。

plot(roc, main = paste0("ROC Curve for Disease Prediction of ", colnames(TrainData)[1]), # 设置主标题col = "#DD7123", # 设置曲线颜色lwd = 2,  # 线条宽度print.auc = T, # 打印AUC的值print.auc.x = 0.4, # AUC值得位置（x轴）print.auc.y = 0.5, # AUC值得位置（y轴）print.auc.pattern = 'AUC=%.3f', # AUC值得格式，表明保留三位小数print.auc.cex = 1.2,  # AUC值字体大小grid = c(0.5,0.2), # 网格线设置grid.col = "gray", # 网格线颜色xlab = "False Positive Rate (1 - Specificity)",  # x轴标签  ylab = "True Positive Rate (Sensitivity)",  # y轴标签 font.lab = 2,  # 轴标签字体样式  cex.axis = 1.2,  # 轴标签字体大小  cex.main = 1.5  # 主标题字体大小  )

ROC曲线如下图所示

横坐标 —— 假阳性率（False Positive Rate, FPR），也称为1-特异性或误报率。它表示在所有实际为阴性的样本中，被错误地判断为阳性的比例。FPR的值越接近0，说明模型的误报率越低，即模型在判断为阳性的样本中，真正为阴性的样本占比越少。
纵坐标 —— 真阳性率（True Positive Rate, TPR）：也称为敏感度或真正率。它表示在所有实际为阳性的样本中，被正确地判断为阳性的比例。TPR的值越接近1，说明模型的敏感度越高，即模型能够准确地识别出更多的真正阳性样本。
曲线位置 —— ROC曲线越靠近左上角（FPR越小，TPR越大），说明模型的预测准确率越高，这是因为左上角的点代表了最低的误报率和最高的真正率。
曲线下面积（AUC） —— AUC值用于量化地表示模型的预测准确性。AUC值越高，说明模型的预测性能越好。AUC值在0.5到1之间，当AUC值等于0.5时，表示模型没有预测价值（即随机猜测）；当AUC值大于0.7时，通常认为具有较高的诊断价值。

在这里插入图片描述