【直播】我的基因组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

图文编辑：吃瓜群众