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

生信技能树

发布于 2018-03-08 11:00:34

1.9K0

发布于 2018-03-08 11:00:34

好久不见，我们的直播又开始啦！今天，我们主要讲的是人群分布，先用简单的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')