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

相关文章

来自专栏专知

Seq2seq强化学习实战 (Pytorch, Tensorflow, Theano)

【导读】本文是Kirti Bakshi在1月14日写的关于其强化学习课程的一个介绍,作者首先简单介绍了机器学习的缺点,以及为什么使用深度学习。然后讲述了其开设的...

3625
来自专栏数据派THU

这5个机器学习项目你不可错过!(附代码)

本文共2299字,建议阅读6分钟。 本文将给大家介绍五个十分可怕但还鲜为人知的机器学习项目,囊括了一些潜在的机器学习的新想法。

543
来自专栏月色的自留地

从锅炉工到AI专家(11)(END)

1317
来自专栏CVer

【CVPR 2018】979篇录用论文合集下载

CVPR 2018 共计录用979篇论文,现已将所有 PDF文件打包成一个文件夹,并提供检索表格。前几天,CVPR 2018大会上刚刚发布了最佳论文奖、学生最佳...

882
来自专栏美团技术团队

深度学习及AR在移动端打车场景下的应用

作为美团点评技术团队的传统节目,每年两次的Hackathon已经举办多年,产出很多富于创意的产品和专利,成为工程师文化的重要组成部分。本文就是2017年冬季Ha...

3229
来自专栏养码场

最适合练手30个的机器学习开源项目,赶紧收藏!

“ 场主,这篇文章炒鸡棒!内涵许多实战项目,很适合机器学习刚入门的小伙伴磨练来提升自己的技术水平。这些优质的开源项目都来自于GitHub上,排名十分靠前,反正很...

621
来自专栏AI研习社

终于!Supervise.ly 发布人像分割数据集啦(免费开源)

本文为雷锋字幕组编译的技术博客,原标题Releasing “Supervisely Person” dataset for teaching machines ...

822
来自专栏AI科技大本营的专栏

Facebook开源多款AI工具,支持游戏、翻译等

近日,Facebook 在年度开发者大会 F8 上宣布开源多款 AI 工具,除了 PyTorch、Caffe 等深度学习框架之外,此次开源的还包括 DenseP...

511
来自专栏大数据文摘

边玩边入门深度学习,我们帮你找了10个简易应用demo

1793
来自专栏新智元

AI解决密码学家终极挑战,600年未解伏尼契手稿有望破译

来源:gizmodo.com 编译:马文 【新智元导读】伏尼契手稿是一本内容不明的神秘书籍,里面充满着神秘的文字和插图。自从100多年前被发现以来,无数语言学家...

33610

扫描关注云+社区