【直播】我的基因组58:用R包SNPRelate来对我的基因型跟hapmap计划数据比较

hapmap计划的人群分布结果和千人基因组计划的分布结果来分析是一样的!【直播】我的基因组55:简单的PCA分析千人基因组的人群分布

这两个计划里面收集的样本的种群信息都比较完善,而且每个样本的基因型信息很容易就下载了。

但是hapmap收集的样本要比千人基因组计划少一些,如下:

数据下载见前面的系列贴

mkdir -p ~/annotation/variation/human/hapmapcd ~/annotation/variation/human/hapmap# ftp://ftp.ncbi.nlm.nih.gov/hapmap/wget ftp://ftp.ncbi.nlm.nih.gov/hapmap/phase_3/relationships_w_pops_051208.txtnohup wget -c -r -np -k -L -p  -nd -A.gz ftp://ftp.ncbi.nlm.nih.gov/hapmap/phase_3/hapmap3_reformatted &

这样就得到了hapmap计划涉及到的所有样本的基因型文件。

然后再学一下SNPRelate的用法:

说明书还比较好读:http://corearray.sourceforge.net/tutorials/SNPRelate/

只有一个核心函数,就是用snpgdsPCA来对包含了GDS格式的基因型信息的文件做分析!

所以重点就是创建GDS格式的基因型文件!

有两种方式来创建GDS文件,被R包作者包装成了两个函数:分别是snpgdsCreateGeno和snpgdsVCF2GDS

其中snpgdsCreateGeno需要自己导入6个数据,比较复杂,第一个是genmat,每个样本在每个位点的基因型(0,1,2)矩阵,然后是sample.id(共279个)和snp.id(共1000个)看名字就知道是样本编号和位点的编号,然后是snp.chromosome和snp.position记录着那1000个snp位点的染色体及坐标信息,最后是snp.allele说明该位点是由什么突变到什么的。

而snpgdsVCF2GDS可以直接读取多样本的VCF文件,一般来说需要自己把多个样本的vcf文件合并成一个,稍微简单一点!

创建好的GDS文件,可以用openfn.gds,index.gdsn,read.gdsn,closefn.gds函数来操作,但是意义不大,我们只需要做PCA分析即可。

包说明书介绍的代码如下,我添加了注释,很简单就可以看懂!

data(hapmap_geno)## you need to create this data by yourself.# Create a gds filesnpgdsCreateGeno("test.gds", genmat = hapmap_geno$genotype,sample.id = hapmap_geno$sample.id, snp.id = hapmap_geno$snp.id,snp.chromosome = hapmap_geno$snp.chromosome,snp.position = hapmap_geno$snp.position,snp.allele = hapmap_geno$snp.allele, snpfirstdim=TRUE)# Open the GDS file(genofile <- snpgdsOpen("test.gds"))## 需要详细理解 genofile 这个对象里面包含的数据内容RV <- snpgdsPCA(genofile, num.thread=2)## 做PCA分析的时候不需要样本的种群信息,但是画图的时候需要,可以看看聚类是否符合认知。pop <- read.gdsn(index.gdsn(genofile,path="sample.annot/pop.group"))## 如果你没有在 snpgdsCreateGeno 里面添加 sample.anno详细,那么上面这个代码是无效的,不过你可以直接赋值pop,就是一个向量,指明你的sample.id(共279个)所属种群即可。plot(RV$eigenvect[,2], RV$eigenvect[,1],col=as.integer(factor(pop)),xlab="PC 2", ylab="PC 1")legend("topleft", legend=levels(factor(pop)), pch="o", col=1:4)

我就基于前面对千人基因组计划数据的探索来使用这个包:

根据我对这个包的学习,目前我只有我挑选的snp位点的dbSNP的ID,并没有保留它们的染色体坐标以及突变形式,我需要重新再写个程序,支持直接去dbSNP数据库里面搜索即可。

zcat ~/annotation/variation/human/dbSNP/All_20160601.vcf.gz |perl -alne 'BEGIN{open FH,"/home/jianmingzeng/biosoft/fastpop/FastPop/snp.txt";while(<FH>){chomp;$h{$_}=1};close FH}{print "$F[2]\t$F[0]\t$F[1]\t$F[3]/$F[4]" if exists $h{$F[2]}}' >fastpop.dbSNP

还是挑选前面的fastpop软件的那两千多个位点吧!就对上面下载的数据进行批量处理:

ls ~/annotation/variation/human/hapmap/*gz |while read iddoecho $idfile=$(basename $id )pop=${file%%.*}zcat $id |perl -alne 'BEGIN{open FH,"/home/jianmingzeng/biosoft/fastpop/FastPop/snp.txt";while(<FH>){chomp;$h{$_}=1};close FH}{print join("\t",$F[0],@F[11..$#F]) if exists $h{$F[0]}}' >$pop.choose.genotypedone

生成了11个种群的genotype文件,然后用下面的R代码处理。

listFiles=list.files("./","*genotype")ASW <- read.table(listFiles[1],stringsAsFactors = F);ASW[1:4,1:4]sample_list<-paste("ASW",1:(ncol(ASW)-1),sep = '_')pop <- rep("ASW",ncol(ASW)-1)for (f in listFiles[2:length(listFiles)] ){this_pop=strsplit(f,'\\.')[[1]][1];tmp <- read.table( f ,stringsAsFactors = F);tmp[1:4,1:4]pop <- c(pop,rep(this_pop,ncol(tmp)-1))sample_list<-c(sample_list,paste(this_pop,1:(ncol(tmp)-1),sep = '_'))ASW <- merge(ASW,tmp,by='V1')}exprSet <- ASWcolnames(exprSet)=c('rsID',sample_list)exprSet[1:4,1:4]snp_info=read.table('fastpop.dbSNP',stringsAsFactors = F)head(snp_info)snp_info <- snp_info[match(as.character(exprSet$rsID),snp_info[,1]),]genotype <- exprSet[,-1]genotype <- apply(genotype,1,function(x){as.numeric(as.factor(x))})genotype <- t(genotype)-1dim(genotype);genotype[1:4,1:4]library(gdsfmt)library(SNPRelate)data(hapmap_geno)# Create a gds filesnpgdsCreateGeno("hapmap.gds", genmat = genotype,sample.id = as.character(sample_list), snp.id = as.character(exprSet[,1]),snp.chromosome = snp_info[,2],snp.position = snp_info[,3],snp.allele = snp_info[,4], snpfirstdim=TRUE)# Open the GDS file(genofile <- snpgdsOpen("hapmap.gds"))table(pop);super_pop=poptable(super_pop)RV <- snpgdsPCA(genofile, num.thread=2)pc.percent <- RV$varprop*100head(round(pc.percent, 2))plot(RV$eigenvect[,2], RV$eigenvect[,1], col=as.integer(factor(super_pop)),xlab="PC 2", ylab="PC 1")legend("bottomleft", y.intersp=0.3,legend=levels(factor(super_pop)), pch="o", col=1:length(super_pop))# close the genotype fileclosefn.gds(genofile)

人种太多了,上色就很麻烦,我也懒得把我自己的基因型放进去了,比较千人基因组计划的分析结果挺好的。

这个hapmap首先基因型就是通过芯片得到的,准确性没有千人基因组计划的测序数据好。

参考文献:

http://www.stat-gen.org/tut/tut_preproc.html

https://wurmlab.github.io/genomicscourse/2016-SIB/practicals/population_genetics/popgen

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-02-18

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏点滴积累

geotrellis使用(三十一)使用geotrellis直接将GeoTiff发布为TMS服务

前言 传统上我们需要先将Tiff中存储的影像等数据先切割成瓦片,而后再对外提供服务。这样的好处是服务器响应快,典型的用空间来换时间的操作。然而这样造成的问题是空...

2689
来自专栏SeanCheney的专栏

Python解答LeetCode

古有科举八股,今有LeetCode。 八股定格式而取文采心意,LeetCode定题目且重答案背诵。 美其名曰:"practice makes perfec...

47816
来自专栏大数据挖掘DT机器学习

中文分词实践(基于R语言)

背景:分析用户在世界杯期间讨论最多的话题。 思路:把用户关于世界杯的帖子拉下来,然后做中文分词+词频统计,最后将统计结果简单做个标签云. 后续:中文分词是中文...

3386
来自专栏华章科技

大白话讲解遗传算法

种群(Population):生物的进化以群体的形式进行,这样的一个群体称为种群。

661
来自专栏从零开始学 Web 前端

从零开始学 Web 之 CSS3(八)CSS3三个案例

而分辨率则一般用像素来度量 px,表示屏幕水平和垂直方向的像素数,例如 1920*1080 指的是屏幕垂直方向和水平方向分别有1920和1080个像素点而构成。

911
来自专栏ATYUN订阅号

用香蕉也能玩电脑游戏—Tensorflow对象检测接口的简单应用

Tensorflow最近发布了用于对象检测的对象检测接口(Object Detection API),能够定位和识别图像中的对象。它能够快速检测图像允许从视频帧...

3284
来自专栏新智元

【重磅】谷歌TensorFlow 1.0发布,智能手机也能玩转深度学习

【新智元导读】 近日,谷歌开源深度学习框架 TensorFlow 发布了完整的1.0版本,不仅改进了库中的机器学习功能,而且对 Python 和 Java 用户...

2777
来自专栏james大数据架构

创建支持多种屏幕尺寸的Android应用

Android涉及各种各样的支持不同屏幕尺寸和密度的设备。对于应用程序,Android系统通过设备和句柄提供了统一的开发环境,大部分工作是校正每一个应用程序的用...

1826
来自专栏Golang语言社区

数据说话:Go语言的Switch和Map性能实测

在开发pgx(一个针对Go语言的PostgreSQL driver)的时候,有好几次我都需要在20多个代码分支间跳转。通常我会选用switch语句。还有个更加可...

3555
来自专栏恰同学骚年

借助 Lucene.Net 构建站内搜索引擎(上)

前言:最近翻开了之前老杨(杨中科)的Lucene.Net站内搜索项目的教学视频,于是作为老杨脑残粉的我又跟着复习了一遍,学习途中做了一些笔记也就成了接下来您看到...

1072

扫码关注云+社区