跟着Molecular Ecology学数据分析：使用R语言对群体SNP数据做主成分分析

论文

The population structure and recent colonization history of Oregon threespine stickleback determined using restriction-site associated DNA-sequencing

本地pdf文件 nihms465650.pdf

image.png

这个论文对应的数据是可以公开下载的

image.png

找到了一本电子书 https://bookdown.org/hhwagner1/LandGenCourse_book/ 里面用到这篇文章的数据做了群体PCA，今天的推文我们试着重复一下这本电子书中的代码

如果要用这个数据的话首先得安装R包

devtools::install_github("hhwagner1/LandGenCourse")
devtools::install_github("hhwagner1/LandGenCourseData")

读取数据集

require("LandGenCourse")

example_df<-system.file("extdata", "stickleback_data.txt", 
            package = "LandGenCourseData")
data <- read.table(example_df, 
                   sep="\t", 
                   as.is=T,
                   check.names=F)

dim(data)
data[1:5,1:10]

看到了书里介绍library()和require()函数的区别：

library() 加载的包即使之前已经加载过了还是会加载一遍
require() 如果之前加载过就不会再加载了

数据集应该是行是样本，列是位点，总共571个样本，10000个位点

生成每个样本属于哪个群体

pops <- unique(unlist(lapply(rownames(data),
                             function(x){y<-c();y<-c(y,unlist(strsplit(x,"_")[[1]][1]))})))  

pops

sample_sites <- rep(NA,nrow(data))
for (i in 1:nrow(data)){
  sample_sites[i] <- strsplit(rownames(data),"_")[[i]][1]}
N <- unlist(lapply(pops,function(x){length(which(sample_sites==x))}))
names(N) <- pops
N

主成分分析

pcaS<-prcomp(data,center = T)

获取每个主成分所解释的变异占比

perc <- round(100*(pcaS$sdev^2 / sum(pcaS$sdev^2))[1:10],2)
names(perc) <- apply(array(seq(1,10,1)), 1, function(x){paste0("PC", x)})
perc 

用ggplot2来做个图

首先是生成9个颜色

colors <- RColorBrewer::brewer.pal(9, "Paired") 

作图

library(ggplot2)

ggplot(data=pca.df,aes(x=PC1,y=PC2))+
  geom_point(aes(color=pop))+
  scale_color_manual(values = colors)+
  theme_bw()+
  labs(x="PC1 (37.47%)",
       y="PC2 (3.78%)")

image.png

用主成分3,4再做一个图

pca.df34<-data.frame(pcaS$x[,3:4],
                   pop=sample_sites)
head(pca.df34)

ggplot(data=pca.df34,aes(x=PC3,y=PC4))+
  geom_point(aes(color=pop),size=5)+
  scale_color_manual(values = colors)+
  theme_bw()+
  labs(x="PC1 (2.55%)",
       y="PC2 (0.88%)")+
  theme(legend.position = "top",
        legend.title = element_text(hjust=0.5))+
  guides(color=guide_legend(nrow=1,
                            title.position = "top"))

image.png

最后是拼图

library(patchwork)

p1+p2+plot_layout(guides = "collect")&theme(legend.position = "top")

image.png

推文的示例数据和代码大家可以直接找到开头提到的在线电子书，或者直接给推文打赏1元，入股打赏了没有收到我的回复，可以加我的微信mingyan24催我

欢迎大家关注我的公众号

小明的数据分析笔记本

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

后续还要仔细看看这篇论文，看看其中对结果的解释