论文
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、生物信息学入门学习资料及自己的学习笔记！
后续还要仔细看看这篇论文，看看其中对结果的解释
