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

跟着Science正刊论文学数据分析:R语言做Kruskal-Wallis秩和检验 by 小明的数据分析笔记本

首先推荐一本算法相关的书

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

论文

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)

欢迎大家关注我的公众号

小明的数据分析笔记本

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

ixxmu / mp_duty