R语言做双侧威尔克森（Wilcoxon）配对符号秩检验简单小例子

小明的数据分析笔记本

示例数据来源于论文

A disease-associated gene desert directs macrophage inflammation through ETS2

https://www.nature.com/articles/s41586-024-07501-1

这个论文里提供了很多组学数据处理流程的代码

代码链接

https://github.com/JamesLeeLab/ETS2manuscript_Stankey_CT_et_al_2024/tree/main

Chip-seq ATAC-seq CUT&RUN 单细胞转录组 空间转录组

有时间可以仔细看看

今天的推文介绍一下论文中fig1k右侧的小图的R代码

示例数据截图

读取数据

library(readxl)
dat<-read_excel("2024.data/20240607/fig1k.xlsx")
dat

双侧威尔克森（Wilcoxon）配对符号秩检验，论文中图注部分的写法是

Statistical analysis was performed using two-sided Wilcoxon matched-pair tests

代码

wilcox.test(dat$`rs2836882: G`,dat$`rs2836882: A`,
            alternative = "two.side",
            paired = TRUE)

这里有两个警告信息

第二个 cannot compute exact p-value with zeroes，我在通义千问上问了一下这个原因，给出的回答是

在使用R语言的wilcox.test函数进行配对检验时，如果遇到提示“cannot compute exact p-value with zeroes”，这意味着在您的数据中存在零值（0），这导致无法计算精确的p值。 Wilcoxon符号秩检验（也称为Wilcoxon配对检验）是一种非参数统计方法，用于比较配对数据集的中位数是否有显著差异。该检验在计算过程中依赖于数据的排序和符号，而零值在这一上下文中可能引起问题，尤其是在尝试进行精确概率计算时，因为它影响了秩的分配。

这个数据集里确实存在两组数据一样的情况，如果把数据一样的去掉就不会有这个警告信息

new.dat<-dat %>% filter(Donor!="Heterozygote7")
wilcox.test(new.dat$`rs2836882: G`,new.dat$`rs2836882: A`,
            alternative = "two.side",
            paired = TRUE)

第一个警告信息

cannot compute exact p-value with ties

通义千问给出的回答是

在使用R语言的wilcox.test函数进行配对检验时，如果遇到提示“cannot compute exact p-value with ties”，这是因为Wilcoxon符号秩检验在计算过程中遇到了“相同值”（ties）的问题。 Wilcoxon检验是基于数据的排序来工作的，它将每个配对的差值按照绝对值大小排序，然后赋予每个差值一个秩（rank）。当存在两个或多个相同的差值（即“ties”）时，这些差值会获得相同的秩。这种情况下，传统的精确概率计算方法可能无法直接应用，因为相同的秩会导致某些排列的计算变得复杂或不可行，特别是在样本量较小或ties较多的情况下。

就是如果有两组数据的差值是一样的就会有这个警告信息

new.dat01<-dat %>% filter(Donor!="Heterozygote7"&Donor!="Heterozygote11"&Donor!="Heterozygote14")
wilcox.test(new.dat01$`rs2836882: G`,new.dat01$`rs2836882: A`,
            alternative = "two.side",
            paired = TRUE)

或者在函数里加上exact=FALSE参数就没有警告信息了，如果把差值是0或者差值相等的数据去掉一组，p 值会变小

作图代码

dat %>% 
  pivot_longer(!Donor) %>% 
  mutate(x=case_when(
    name == "rs2836882: G" ~ 1,
    name == "rs2836882: A" ~ 2
  )) -> dat.plot01


dat %>% 
  mutate(xA=2,xG=1) -> dat.plot02
ggplot()+
  geom_segment(data=dat.plot02,
               aes(x=xG,xend=xA,y=`rs2836882: G`,yend=`rs2836882: A`))+
  geom_point(dat=dat.plot01,
             aes(x=x,y=value,color=name),size=5,
             show.legend = FALSE)+
  theme_bw(base_size = 15)+
  theme(panel.grid = element_blank())+
  scale_color_manual(values = c("red","blue"))+
  scale_x_continuous(breaks = c(1,2),
                     limits = c(0.5,2.5),
                     labels = c("G","A"))+
  labs(x="rs2836882",y="ATAC-seq reads\n(individual)")+
  theme(panel.border = element_blank(),
        axis.line = element_line())+
  annotate(geom = "segment",x=1,xend=2,y=48,yend=48)+
  annotate(geom = "text",x=1.5,y=50,label="P=0.0001",size=8)

发现一个小问题，论文中的图展示的是G的基因型的值是大于A的基因型的，但是提供的原始数据里G是小于A的。不知道是我理解的有问题还是这个论文里提供的原始作图数据的一个小错误，大家找到论文和原始数据来对照着看一下。

欢迎大家关注我的公众号 小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！

微信扫一扫分享文章