示例数据来源于论文
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
image.png
Chip-seq ATAC-seq CUT&RUN 单细胞转录组 空间转录组
有时间可以仔细看看
今天的推文介绍一下论文中fig1k右侧的小图的R代码
示例数据截图
image.png
读取数据
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)
image.png
这里有两个警告信息
第二个 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)
image.png
发现一个小问题，论文中的图展示的是G的基因型的值是大于A的基因型的，但是提供的原始数据里G是小于A的。不知道是我理解的有问题还是这个论文里提供的原始作图数据的一个小错误，大家找到论文和原始数据来对照着看一下。
image.png
欢迎大家关注我的公众号
小明的数据分析笔记本
声明：文中观点不代表本站立场。本文传送门：https://eyangzhen.com/417094.html
