零基础入门转录组数据分析——基因Wilcoxon秩和检验
目录
- 零基础入门转录组数据分析——基因Wilcoxon秩和检验
- 1. 单基因Wilcoxon秩和检验的基础知识
- 2. 基因Wilcoxon秩和检验(Rstudio)——代码实操
- 2. 1 数据处理
- 2. 2 基因Wilcoxon秩和检验
- 2. 3 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())
Wilcoxon秩和检验结果如下图所示,横坐标为不同的基因名称,纵坐标为不同基因的表达水平,图中黄色的箱子表示疾病组,蓝色的箱子表示对照组,最上方的*表示显著性结果。
结语:
以上就是基因Wilcoxon秩和检验的所有过程,如果有什么需要补充或不懂的地方,大家可以私聊我或者在下方评论。
如果觉得本教程对你有所帮助,希望广大学习者能够花点自己的小钱支持一下(点赞旁的打赏按钮)作者创作(可以的话一杯蜜雪奶茶即可),感谢大家的支持~~~~~~ ^_^ !!!
祝大家能够开心学习,轻松学习,在学习的路上少一些坎坷~~~
- 目录部分跳转链接:零基础入门生信数据分析——导读