搜索
查看: 3233|回复: 0

【直播】我的基因组58:用R包SNPRelate来对我的基因型跟hapma...

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-25 16:54:52 | 显示全部楼层 |阅读模式
【直播】我的基因组58:用R包SNPRelate来对我的基因型跟hapmap计划数据比较

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

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

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



数据下载见前面的系列贴

  • [AppleScript] 纯文本查看 复制代码
    mkdir -p ~/annotation/variation/human/hapmap
    cd ~/annotation/variation/human/hapmap
    # [url]ftp://ftp.ncbi.nlm.nih.gov/hapmap/[/url]
    wget [url]ftp://ftp.ncbi.nlm.nih.gov/hapmap/phase_3/relationships_w_pops_051208.txt[/url]
    nohup wget -c -r -np -k -L -p  -nd -A.gz [url]ftp://ftp.ncbi.nlm.nih.gov/hapmap/phase_3/hapmap3_reformatted[/url] &


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


然后再学一下SNPRelate的用法:

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

只有一个核心函数,就是用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分析即可。



包说明书介绍的代码如下,我添加了注释,很简单就可以看懂!
  • [AppleScript] 纯文本查看 复制代码
    data(hapmap_geno)
    ## you need to create this data by yourself.
    # Create a gds file
    snpgdsCreateGeno("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数据库里面搜索即可。

  • [AppleScript] 纯文本查看 复制代码
    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软件的那两千多个位点吧!就对上面下载的数据进行批量处理:

  • [AppleScript] 纯文本查看 复制代码
    ls ~/annotation/variation/human/hapmap/*gz |while read id
    do
    echo $id
    file=$(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.genotype
    done


生成了11个种群的genotype文件,然后用下面的R代码处理。
  • [AppleScript] 纯文本查看 复制代码
    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 <- ASW
    colnames(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)-1
    dim(genotype);genotype[1:4,1:4]
    library(gdsfmt)
    library(SNPRelate)
    data(hapmap_geno)
    # Create a gds file
    snpgdsCreateGeno("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=pop
    table(super_pop)
    RV <- snpgdsPCA(genofile, num.thread=2)
    pc.percent <- RV$varprop*100
    head(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 file
    closefn.gds(genofile)


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

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


参考文献:

[url=]http://www.stat-gen.org/tut/tut_preproc.html[/url]

[url=]https://wurmlab.github.io/genomicscourse/2016-SIB/practicals/population_genetics/popgen[/url]




本帖子中包含更多资源

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

x



上一篇:【直播】我的基因组57:最简陋的祖源分析
下一篇:【直播】我的基因组59:把我的数据伪装成23andme或wegene的...
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-16 09:16 , Processed in 0.041539 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.