【直播】我的基因组55:简单的PCA分析千人基因组的人群分布

好久不见,我们的直播又开始啦!今天,我们主要讲的是人群分布,先用简单的PCA来分析一下千人基因组的人群分布吧!

PCA分析,就是主成分分析,我博客有讲过(点击最底部的阅读原文或复制链接http://www.bio-info-trainee.com/1232.html进行查看)。 PCA的原本目的是因为变量太多,想把它们合并成两三个变量,从而简化分析步骤。变量的多少代表维度的多少,一千维的数据已经无法想象了,但是二维和三维还是比较符合认知的。假设用PCA给千人基因组所有个体一个二维坐标,画在图上,就可以清清楚楚看到他们的分布了,是不是有规律的聚集成一个个种群。

主成分分析可以得到p个主成分,但是,由于各个主成分的方差是递减的,包含的信息量也是递减的,所以实际分析时,一般不是选取p个主成分,而是根据各个主成分累计贡献率的大小选取前k个主成分。这里贡献率就是指某个主成分的方差占全部方差的比重,实际也就是某个特征值占全部特征值总和的比重。贡献率越大,说明该主成分所包含的原始变量的信息越强。主成分个数k的选取,主要根据主成分的累积贡献率来决定,即一般要求累计贡献率达到85%以上,这样才能保证综合变量能包括原始变量的绝大多数信息。

先下载千人基因组计划的数据吧(http://www.1000genomes.org/category/population/)涉及到2504个人,5大人种,26个亚人种!信息都存储在integrated_call_samples_v3.20130502.ALL.panel文件里面。

下载千人基因组计划所有人的基因型的方法我以前讲过

  1. mkdir -p ~/annotation/variation/human/1000genomes
  2. cd ~/annotation/variation/human/1000genomes
  3. ##ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/
  4. nohup wget -c -r -nd -np -k -L -p ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502 &

我们随便打开一个染色体看看里面的数据是什么~

  1. zcat ALL.chr1.phase3_shapeit2_mvncall_integrated_v3plus_nounphased.rsID.genotypes.GRCh38_dbSNP_no_SVs.vcf.gz |cut -f 1-3

一般全基因组数据都是成千上万的位点,我没有看到教程告诉我如何挑选位点,比如http://online.cambridgecoding.com/notebooks/cca_admin/genetic-ancestry-analysis-python 我只能凭感觉挑选 allele frequency 在0.5附近的~

  1. zcat ALL.chr1.phase3_shapeit2_mvncall_integrated_v3plus_nounphased.rsID.genotypes.GRCh38_dbSNP_no_SVs.vcf.gz |perl -alne '{/AF=(.*?);/;print join("\t",$F[2],@F[9..$#F]) if $1<0.51 && $1>0.49}' >chr1.alf.test

即使这样,还是有9520个位点,这样太多了,画图会比较卡,我就挑选了前面的一千个位点来做的。

  1. dat[1:4,1:4]## 做好这个数据即可
  2. # apply PCA - scale. = TRUE is highly
  3. # advisable, but default is FALSE.
  4. dat.pca <- prcomp(dat,
  5. center = TRUE,
  6. scale. = TRUE)
  7. # print method
  8. print(dat.pca)
  9. # plot method
  10. plot(dat.pca, type = "l")
  11. summary(dat.pca)
  12. biplot (dat.pca , scale =0,var.axes =F)
  13. group_info <-read.table('integrated_call_samples_v3.20130502.ALL.panel',header =T,stringsAsFactors = F)
  14. head(group_info)
  15. pop = group_info[match(colname,group_info$sample ),'pop']
  16. super_pop = group_info[match(colname,group_info$sample),'super_pop']
  17. #library(devtools)
  18. #install_github("ggbiplot", "vqv")
  19. dat_tmp= as.data.frame(dat)
  20. dat_tmp$sample = rownames(dat_tmp)
  21. dat_df <- merge(dat_tmp,group_info,by='sample')
  22. library(ggbiplot)
  23. g <- ggbiplot(dat.pca, obs.scale = 1 ,var.scale = 1, var.axes=F,
  24. groups = super_pop, ellipse = TRUE,
  25. circle = TRUE)
  26. g <- g + scale_color_discrete(name = '')
  27. g <- g + theme(legend.direction = 'horizontal',
  28. legend.position = 'top')
  29. print(g)
  30. library(ggfortify)
  31. autoplot(prcomp(dat), data = dat_df, colour = 'super_pop')

还是老规矩,通过Google我找到了可视化方法!

用谷歌搜索来使用ggplot2做可视化(下)

就是上面代码中的ggbiplot和ggfortify包,很容易就把千人基因组按照5个种群给分开了,当然,如果按照26个亚种会很难看,我就不秀图片了!而且其实前两个主成分的贡献度都很低,说明它们两个的把人群分开的作用力不够。

首先是ggbiplot的图片!

然后是ggfortify 图片:

我只是随机挑选的千人基因组计划的1号染色体的1000个位点!我们看到,我们的数据区分的不是很明显,我挑选的1000个位点没办法把人群清晰分开(前两个主成分作用力太小了),刚开始我选择的是26个人种,更加麻烦,现在就标记5个超级人种,勉强还能看到规律。非洲人跟其他人区分挺好的。

  • AFR, African
  • AMR, Ad Mixed American
  • EAS, East Asian
  • EUR, European
  • SAS, South Asian

对于这个特定的分析,PCA一个最强大的部分是,一开始,我们无需指定寻找的集群数。

你们觉得哪一个好看呢?(投票ing)

参考文献:

https://www.r-bloggers.com/computing-and-visualizing-pca-in-r/

https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_pca.html

https://www.analyticsvidhya.com/blog/2016/03/practical-guide-principal-component-analysis-python/

http://stats.stackexchange.com/questions/72839/how-to-use-r-prcomp-results-for-prediction

文:Jimmy

图文编辑:吃瓜群众

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-02-13

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏木宛城主

曾今的代码系列——自己的分页控件+存储过程实现分页

项目里面的测试代码,仅供参考 LoginByAjax <title>Ajax登陆</title> <script src="Scripts/c...

1865
来自专栏菩提树下的杨过

基于sliverlight + wcf的web 文字版IM 示例

演示地址: http://task.24city.com/default.html 预览界面: ? 一、布局 采用Grid布局,5行2列 第一行:为登录/注册信...

3226
来自专栏跟着阿笨一起玩NET

从sql server 中读取二进制图片

401
来自专栏xingoo, 一个梦想做发明家的程序员

【插件开发】—— 6 SWT 复杂控件使用以及布局

前文回顾: 1 插件学习篇 2 简单的建立插件工程以及模型文件分析 3 利用扩展点,开发透视图 4 SWT编程须知 5 SWT简单控件的使用与布局搭...

2349
来自专栏张善友的专栏

弹出式模态窗体选择文本控件

2006年就要到来了,最近比较忙,很少更新blog,今天发一个模态窗体选择文本控件辞旧迎新.新年在发几个asp.net2.0 webPart控件同各位分享: ...

1907
来自专栏我和未来有约会

Silverlight制作逐帧动画 v2 - part2

Silverlight制作逐帧动画 v2 - part2 在这里完善了一下算法,加入了fps的机制进去。 private string[] ...

1896
来自专栏james大数据架构

CSS好看的按钮

好看的按钮 <style> .btn { BORDER-RIGHT: #7b9ebd 1px solid; PADDING-RIGHT: 2px; BORDE...

2007
来自专栏跟着阿笨一起玩NET

treeview 绑定文件夹和文件

461
来自专栏互联网开发者交流社区

STC-单片机控制系统

1113
来自专栏Golang语言社区

GO语言 TCP传输实例

package main import ( "net" "fmt" ) var ( maxRead = 1100 msgStop = []byt...

3406

扫码关注云+社区