前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >跟着Molecular Ecology学数据分析:使用R语言对群体SNP数据做主成分分析

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

作者头像
用户7010445
发布2021-12-01 18:15:26
8930
发布2021-12-01 18:15:26
举报

论文

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

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2021-11-26,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 小明的数据分析笔记本 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 论文
  • 读取数据集
  • 生成每个样本属于哪个群体
  • 主成分分析
  • 获取每个主成分所解释的变异占比
  • 用ggplot2来做个图
  • 用主成分3,4再做一个图
  • 最后是拼图
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档