【直播】我的基因组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 条评论
登录 后参与评论

相关文章

来自专栏DHUtoBUAA

粗略的物体碰撞预测及检测

  该博客实时更新于我的Github。

3386
来自专栏机器学习、深度学习

人群运动--Scene-Independent Group Profiling in Crowd

Scene-Independent Group Profiling in Crowd CVPR2014 http://www.ee.cuhk.edu.hk/...

1719
来自专栏大数据挖掘DT机器学习

玩玩文本挖掘-wordcloud、主题模型与文本分类

本文主要介绍文本挖掘的常见方法,主要包括词频分析及wordcloud展现、主题模型、文本分类、分类评价等。分类主要包括无监督分类(系统聚类、KMeans...

3116
来自专栏生信技能树

lncRNA实战项目-第六步-WGCNA相关性分析

WGCNA将lncRNA分成18个模块(3635个lncRNA),空间模块中lncRNA表达呈现明显的组织区域特异性,如:CB (M1, 794个lncRNAs...

92810
来自专栏机器之心

教程 | 使用深度学习进行医疗影像分析:文件格式篇

选自Medium 作者:Taposh Dutta-Roy 机器之心编译 参与:Nurhachu Null、李泽南 今年 3 月,英伟达的 GTC 2017 大会...

2846
来自专栏机器人网

PID控制原理:看完这个故事你就明白了

小明接到这样一个任务:有一个水缸漏水,且漏水的速度是不定的,但要求水面高度维持在某个位置,一旦发现水面高度低于要求位置,就要往水缸里加水。 ? 开始小明用瓢加水...

3205
来自专栏大数据挖掘DT机器学习

用R语言对城管事件数据分析

作者:夏尔康 https://ask.hellobi.com/blog/xiaerkang/3975 这次使用主成分分析主要目的并不是降维,而是分析城管数据中的...

3169
来自专栏PPV课数据科学社区

如何用R语言对城管事件数据分析?

这次使用主成分分析主要目的并不是降维,而是分析城管数据中的事件类别之间是否存在关系,当然,城管事件类型有好几百,这里就只选取从去年九月到目前发生量前十的事件类别...

3328
来自专栏AI研习社

用金庸、古龙群侠名称训练 LSTM,会生成多么奇葩的名字?

AI 研习社按:本文转载自 Magicly 博客,获作者授权。阅读原文请见:http://magicly.me/2017/04/07/rnn-lstm-gene...

35411
来自专栏iOSDevLog

Scikit-Learn教程:棒球分析 (一)

一个scikit-learn教程,通过将数据建模到KMeans聚类模型和线性回归模型来预测MLB每赛季的胜利。

1232

扫码关注云+社区