论文
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包
代码语言：javascript复制devtools::install_github("hhwagner1/LandGenCourse")
devtools::install_github("hhwagner1/LandGenCourseData")
读取数据集
代码语言：javascript复制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个位点
生成每个样本属于哪个群体
代码语言：javascript复制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
主成分分析
代码语言：javascript复制pcaS<-prcomp(data,center = T)
获取每个主成分所解释的变异占比
代码语言：javascript复制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个颜色
代码语言：javascript复制colors <- RColorBrewer::brewer.pal(9, "Paired")
作图
代码语言：javascript复制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再做一个图
代码语言：javascript复制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
最后是拼图
代码语言：javascript复制library(patchwork)
p1 p2 plot_layout(guides = "collect")&theme(legend.position = "top")
image.png
