【直播】我的基因組58：用R包SNPRelate來對我的基因型跟hapmap計劃資料比較

阿新 • • 發佈：2022-05-03

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)