零基础入门转录组数据分析——基因Wilcoxon秩和检验

呆猪儿

已于 2024-09-02 09:36:47 修改

阅读量1.4k

点赞数 29

分类专栏：生信之转录组——下游分析文章标签：数据分析数据挖掘算法 r语言笔记数据库大数据

于 2024-08-30 14:20:20 首次发布

本文链接：https://blog.csdn.net/weixin_49878699/article/details/141639844

版权

生信之转录组——下游分析专栏收录该内容

22 篇文章 142 订阅

订阅专栏

零基础入门转录组数据分析——基因Wilcoxon秩和检验

1. 单基因Wilcoxon秩和检验的基础知识

1.1 Wilcoxon秩和检验是什么？
Wilcoxon秩和检验（也称为Mann-Whitney U检验）是一种非参数检验，用于比较两个独立样本的中位数是否存在显著差异，这种检验不假设数据来自正态分布，因此非常适合于小样本或非正态分布的数据。

1.2 Wilcoxon秩和检验的假设？

原假设 —— 两个样本的总体中位数相等
备择假设 —— 两个样本的总体中位数不相等

1.3 p值的意义？
P值表示在零假设为真时，观察到当前或更极端结果的概率。通常，如果P值小于选定的显著性水平（如0.05），则拒绝原假设，认为两个总体的中位数存在显著差异。

1.4 Wilcoxon秩和检验和limma，DESeq2的区别是什么？
最主要的区别就在于Wilcoxon秩和检验不依赖于数据的分布形状，适用于提取出来小部分基因单独进行差异分析，而limma和DESeq2方法则是从宏观整体的层次来分析基因间的差异，会考虑基因彼此间的相互干扰。

综上所述： Wilcoxon秩和检验就是评估两组间基因表达中位值是否存在显著差异，适用于基因数目比较少的情况，例如：通过基因筛选之后仅剩10个左右基因的时候就可以用Wilcoxon秩和检验来比较组间差异了。

注意：做Wilcoxon秩和检验的时候如果数据中存在明显的异常值，需要先进行数据清洗或转换。

2. 基因Wilcoxon秩和检验（Rstudio）——代码实操

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

废话不多说，代码如下：

2. 1 数据处理

设置工作空间：

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

加载包：

library(rstatix)
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 (10个基因，这里用10个基因作为展示)

HubGene <- data.frame(symbol = c('VPS13D', 'MFF', 'ACSL1', 'VDAC1', 'PRELID1', 'BAK1',
                                 'CYCS', 'BCL2L10', 'MPV17L', 'PHB'))

HubGene 如下图所示，只有一列：10个基因的基因名
在这里插入图片描述
从TrainRawData中取出10个基因对应的表达矩阵，并且与之前准备的分组信息表TrainGroup进行合并

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

TrainData 如下图所示，第一列为样本名，第二列为分组情况，后面的都是基因表达量。。
在这里插入图片描述

之后将TrainData 转成长格式数据，关于长宽数据

TrainData <- TrainData %>% 
  pivot_longer(
    cols = -c("sample", "group"),
    names_to = "symbol",
    values_to = "Expression"
  )

转换后的TrainData 如下图所示第一列为样本名，第二列是对应的样本分组，第三列为基因的名称，第四列为不同基因对应的表达量。
在这里插入图片描述

2. 2 基因Wilcoxon秩和检验

接下来用处理好的TrainData 进行Wilcoxon秩和检验

group_by(symbol) —— 是dplyr包中的一个函数，用于按symbol列对数据进行分组。这意味着接下来的操作（如Wilcoxon检验）将针对每个不同的基因独立进行。
wilcox_test(Expression ~ group) —— 是rstatix包中的一个函数，用于执行Wilcoxon秩和检验。这里，它比较了Expression（表达量）在不同group（组别）之间的差异。
adjust_pvalue(method = ‘fdr’) —— 是rstatix包中的一个函数，用于调整p值以控制假发现率（False Discovery Rate, FDR）。使用FDR方法（也称为Benjamini-Hochberg方法）来调整p值，以减少由于多重测试而产生的假阳性结果

WilcoxonResults <- TrainData%>%
  group_by(symbol)%>%
  wilcox_test(Expression ~ group)%>%
  adjust_pvalue(method = 'fdr')

WilcoxonResults 如下图所示，symbol为基因名称；n1和n2分别对应对照和疾病的样本数量；最关键的就是最后两列，p和p.adj对应的是wilcoxon的显著性结果。
（* 表明 p < 0.05，** 表明 p < 0.01， *** 表明 p < 0.001，ns表明没有显著差异）
注意：这里p和p.adj选一个即可，看个人需求

在这里插入图片描述

2. 3 Wilcoxon秩和检验简单可视化

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

ggplot(TrainData, aes(x = symbol, y = Expression, fill = group)) +
  stat_boxplot(geom = "errorbar",
               width = 0.1,
               position = position_dodge(0.9)) +
  geom_boxplot(aes(x = symbol, y = Expression, fill = group),
               width = 0.2,
               position = position_dodge(0.9), 
               outlier.shape = NA, 
               outlier.colour = NA)+ 
  scale_fill_manual(values = c('#355783', "gold"), name = "Group")+
  labs(title = "", x = "", y = "Expression", size = 20) +
  stat_compare_means(data = TrainData,
                     mapping = aes(group = group),
                     label = "p.signif",
                     method = 'wilcox.test',
                     paired = F) +
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5, colour = "black", face = "bold", size = 18),
        axis.text.x = element_text(angle = 45, hjust=1, colour = "black", face = "bold", size = 10), 
        axis.text.y = element_text(hjust = 0.5, colour ="black", face="bold", size=12), 
        axis.title.x = element_text(size = 16, face = "bold"),
        axis.title.y = element_text(size = 16, face = "bold"),
        legend.text = element_text(face = "bold", hjust = 0.5, colour = "black", size = 12),
        legend.title = element_text(face = "bold", size = 12),
        legend.position = "top",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())