搜索
查看: 2551|回复: 0

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

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-25 16:12:10 | 显示全部楼层 |阅读模式
【直播】我的基因组55:简单的PCA分析千人基因组的人群分布

好久不见,我们的直播又开始啦!今天,我们主要讲的是人群分布,先用简单的PCA来分析一下千人基因组的人群分布吧!


PCA分析,就是主成分分析,我博客有讲过(点击最底部的阅读原文或复制链接[url=]http://www.bio-info-trainee.com/1232.html[/url]进行查看)。 PCA的原本目的是因为变量太多,想把它们合并成两三个变量,从而简化分析步骤。变量的多少代表维度的多少,一千维的数据已经无法想象了,但是二维和三维还是比较符合认知的。假设用PCA给千人基因组所有个体一个二维坐标,画在图上,就可以清清楚楚看到他们的分布了,是不是有规律的聚集成一个个种群。




主成分分析可以得到p个主成分,但是,由于各个主成分的方差是递减的,包含的信息量也是递减的,所以实际分析时,一般不是选取p个主成分,而是根据各个主成分累计贡献率的大小选取前k个主成分。这里贡献率就是指某个主成分的方差占全部方差的比重,实际也就是某个特征值占全部特征值总和的比重。贡献率越大,说明该主成分所包含的原始变量的信息越强。主成分个数k的选取,主要根据主成分的累积贡献率来决定,即一般要求累计贡献率达到85%以上,这样才能保证综合变量能包括原始变量的绝大多数信息。




    先下载千人基因组计划的数据吧([url=]http://www.1000genomes.org/category/population/)[/url], 涉及到2504个人,5大人种,26个亚人种!信息都存储在integrated_call_samples_v3.20130502.ALL.panel文件里面。



下载千人基因组计划所有人的基因型的方法我以前讲过

  • [AppleScript] 纯文本查看 复制代码
    mkdir -p ~/annotation/variation/human/1000genomes
    cd ~/annotation/variation/human/1000genomes
    ##[url]ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/[/url]
    nohup wget  -c -r -nd -np -k -L -p  [url]ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502[/url] &
    





我们随便打开一个染色体看看里面的数据是什么~
  • [AppleScript] 纯文本查看 复制代码
    zcat ALL.chr1.phase3_shapeit2_mvncall_integrated_v3plus_nounphased.rsID.genotypes.GRCh38_dbSNP_no_SVs.vcf.gz |cut -f 1-3



一般全基因组数据都是成千上万的位点,我没有看到教程告诉我如何挑选位点,比如[url=]http://online.cambridgecoding.com/notebooks/cca_admin/genetic-ancestry-analysis-python[/url] 我只能凭感觉挑选 allele frequency 在0.5附近的~

  • [AppleScript] 纯文本查看 复制代码
    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个位点,这样太多了,画图会比较卡,我就挑选了前面的一千个位点来做的。


  • [AppleScript] 纯文本查看 复制代码
    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我找到了可视化方法!

用谷歌搜索来使用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)


参考文献:
[url=]https://www.r-bloggers.com/computing-and-visualizing-pca-in-r/[/url]
[url=]https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_pca.html[/url]
[url=]https://www.analyticsvidhya.com/blog/2016/03/practical-guide-principal-component-analysis-python/[/url]
[url=]http://stats.stackexchange.com/questions/72839/how-to-use-r-prcomp-results-for-prediction[/url]

文:Jimmy

图文编辑:吃瓜群众



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:【直播】我的基因组54:把我的variation跟dbSNP数据库相比较
下一篇:【直播】我的基因组56:探索遗传起源
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-11-12 10:15 , Processed in 0.041877 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.