【直播】我的基因组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文件里面。
下载千人基因组计划所有人的基因型的方法我以前讲过
mkdir -p ~/annotation/variation/human/1000genomes
cd ~/annotation/variation/human/1000genomes
##ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/
nohup wget -c -r -nd -np -k -L -p ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502 &
我们随便打开一个染色体看看里面的数据是什么~
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附近的~
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个位点，这样太多了，画图会比较卡，我就挑选了前面的一千个位点来做的。
dat[1:4,1:4]## 做好这个数据即可
# apply PCA - scale. = TRUE is highly
# advisable, but default is FALSE.
dat.pca <- prcomp(dat,
center = TRUE,
scale. = TRUE)
# print method
print(dat.pca)
# plot method
plot(dat.pca, type = "l")
summary(dat.pca)
biplot (dat.pca , scale =0,var.axes =F)
group_info <-read.table('integrated_call_samples_v3.20130502.ALL.panel',header =T,stringsAsFactors = F)
head(group_info)
pop = group_info[match(colname,group_info$sample ),'pop']
super_pop = group_info[match(colname,group_info$sample),'super_pop']
#library(devtools)
#install_github("ggbiplot", "vqv")
dat_tmp= as.data.frame(dat)
dat_tmp$sample = rownames(dat_tmp)
dat_df <- merge(dat_tmp,group_info,by='sample')
library(ggbiplot)
g <- ggbiplot(dat.pca, obs.scale = 1 ,var.scale = 1, var.axes=F,
groups = super_pop, ellipse = TRUE,
circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal',
legend.position = 'top')
print(g)
library(ggfortify)
autoplot(prcomp(dat), data = dat_df, colour = 'super_pop')
还是老规矩，通过Google我找到了可视化方法！
就是上面代码中的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
图文编辑：吃瓜群众