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

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

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

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

数据下载见前面的系列贴

1. mkdir -p ~/annotation/variation/human/hapmap

2. cd ~/annotation/variation/human/hapmap

3. # ftp://ftp.ncbi.nlm.nih.gov/hapmap/

4. wget ftp://ftp.ncbi.nlm.nih.gov/hapmap/phase_3/relationships_w_pops_051208.txt

5. nohup wget -c -r -np -k -L -p  -nd -A.gz ftp://ftp.ncbi.nlm.nih.gov/hapmap/phase_3/hapmap3_reformatted &

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

然后再学一下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分析即可。

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

1. data(hapmap_geno)

2. ## you need to create this data by yourself.

3. # Create a gds file

4. snpgdsCreateGeno("test.gds", genmat = hapmap_geno$genotype,

5. sample.id = hapmap_geno$sample.id, snp.id = hapmap_geno$snp.id,

6. snp.chromosome = hapmap_geno$snp.chromosome,

7. snp.position = hapmap_geno$snp.position,

8. snp.allele = hapmap_geno$snp.allele, snpfirstdim=TRUE)

9. # Open the GDS file

10. (genofile <- snpgdsOpen("test.gds"))

11. ## 需要详细理解 genofile 这个对象里面包含的数据内容

12. RV <- snpgdsPCA(genofile, num.thread=2)

13. ## 做PCA分析的时候不需要样本的种群信息，但是画图的时候需要，可以看看聚类是否符合认知。

14. pop <- read.gdsn(index.gdsn(genofile,path="sample.annot/pop.group"))

15. ## 如果你没有在 snpgdsCreateGeno 里面添加 sample.anno详细，那么上面这个代码是无效的，不过你可以直接赋值pop，就是一个向量，指明你的sample.id(共279个)所属种群即可。

16. plot(RV$eigenvect[,2], RV$eigenvect[,1],col=as.integer(factor(pop)),

17. xlab="PC 2", ylab="PC 1")

18. legend("topleft", legend=levels(factor(pop)), pch="o", col=1:4)

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

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

1. 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软件的那两千多个位点吧！就对上面下载的数据进行批量处理：

1. ls ~/annotation/variation/human/hapmap/*gz |while read id

2. do

3. echo $id

4. file=$(basename $id )

5. pop=${file%%.*}

6. 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

7. done

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

1. listFiles=list.files("./","*genotype")

2. ASW <- read.table(listFiles[1],stringsAsFactors = F);ASW[1:4,1:4]

3. sample_list<-paste("ASW",1:(ncol(ASW)-1),sep = '_')

4. pop <- rep("ASW",ncol(ASW)-1)

5. for (f in listFiles[2:length(listFiles)] ){

6. this_pop=strsplit(f,'\\.')[[1]][1];

7. tmp <- read.table( f ,stringsAsFactors = F);tmp[1:4,1:4]

8. pop <- c(pop,rep(this_pop,ncol(tmp)-1))

9. sample_list<-c(sample_list,paste(this_pop,1:(ncol(tmp)-1),sep = '_'))

10. ASW <- merge(ASW,tmp,by='V1')

11. }

12. exprSet <- ASW

13. colnames(exprSet)=c('rsID',sample_list)

14. exprSet[1:4,1:4]

15. snp_info=read.table('fastpop.dbSNP',stringsAsFactors = F)

16. head(snp_info)

17. snp_info <- snp_info[match(as.character(exprSet$rsID),snp_info[,1]),]

18. genotype <- exprSet[,-1]

19. genotype <- apply(genotype,1,function(x){

20. as.numeric(as.factor(x))

21. })

22. genotype <- t(genotype)-1

23. dim(genotype);genotype[1:4,1:4]

24. library(gdsfmt)

25. library(SNPRelate)

26. data(hapmap_geno)

27. # Create a gds file

28. snpgdsCreateGeno("hapmap.gds", genmat = genotype,

29. sample.id = as.character(sample_list), snp.id = as.character(exprSet[,1]),

30. snp.chromosome = snp_info[,2],

31. snp.position = snp_info[,3],

32. snp.allele = snp_info[,4], snpfirstdim=TRUE)

33. # Open the GDS file

34. (genofile <- snpgdsOpen("hapmap.gds"))

35. table(pop);

36. super_pop=pop

37. table(super_pop)

38. RV <- snpgdsPCA(genofile, num.thread=2)

39. pc.percent <- RV$varprop*100

40. head(round(pc.percent, 2))

41. plot(RV$eigenvect[,2], RV$eigenvect[,1], col=as.integer(factor(super_pop)),

42. xlab="PC 2", ylab="PC 1")

43. legend("bottomleft", y.intersp=0.3,legend=levels(factor(super_pop)), pch="o", col=1:length(super_pop))

44. # close the genotype file

45. closefn.gds(genofile)

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

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

参考文献：