跟着Science正刊论文学数据分析:R语言做Kruskal-Wallis秩和检验

推荐的理由：没有太复杂的数学公式，都是一些简单的小例子，配有一些有趣的插图。涉及到的算法有排序、递归、哈希表、动态规划等等，内容还挺多的。动态规划在生物信息学的序列比对中是非常常用的算法。闲来无事时可以拿出来翻翻。

论文

Reconstruction of the human amylase locus reveals ancient duplications seeding modern-day variation
https://www.science.org/doi/10.1126/science.adn0609
论文作者的推特 https://x.com/evobioclio
之前做过估计基因拷贝数，但是审稿人对估计的准确性一直提出问题，最后也没有给出太好的回复。想看一下这篇论文里具体是怎么研究基因拷贝数的。论文的内容很长，估计需要看好长时间。看的过程中争取把能复现的部分都记录下来。
还有一篇今年发的Nature论文也是研究这个基因位点
https://www.nature.com/articles/s41586-024-07911-1
Recurrent evolution and selection shape structural diversity at the amylase locus
这个论文也好好看看
今天的推文是记录 Kruskal-Wallis秩和检验
论文里写
To evaluate whether the AMY1 gene copy numbers exhibit population-specific patterns, we performed Kruskal-Wallis rank sum test using the copy number information
Kruskal-Wallis秩和检验 主要作用

This non-parametric test was used to determine if there are statistically significant differences between the medians of three or more independent groups
非参数检验，检验3组或者更多组的中位数是否有显著的差异
优势

We chose this test due to its robustness in handling non-normally distributed data
在处理非正态分布数据时具有稳健性。
之前有一篇推文介绍过这个统计检验方法
跟着Nature Genetics学作图：R语言ggplot2箱线图展示不同类别基因家族核苷酸多样性
代码

library(tidyverse)
library(readxl)

dat<-read_excel(“D:/R_4_1_0_working_directory/env001/2024.data/20241018/science.adn0609_tables_s1_to_s15.xlsx”,
sheet = “S3”,
skip = 3)

dat %>% pull(SuperPopulation) %>%
table()

kruskal.test(AMY1 Copy Number~SuperPopulation,
data = dat)

这里总共5组数据，Kruskal-Wallis秩和检验 给出一个P值，想看下两两之间的差异用 Nemenyi test
代码
dat %>%
mutate(SuperPopulation=factor(SuperPopulation)) -> dat

PMCMRplus::kwAllPairsNemenyiTest(AMY1 Copy Number~SuperPopulation,
data = dat)

下三角矩阵对应的应该是P值
画个云雨图

library(smplot2)
ggplot(dat, aes(x = SuperPopulation, y = AMY1 Copy Number, fill = SuperPopulation)) +
sm_raincloud() +
theme(text = element_text(size = 13),
axis.text.x = element_text(angle=0,hjust=0.5,vjust=1))+
labs(x=NULL,y=NULL)

欢迎大家关注我的公众号