【直播】我的基因组58：用R包SNPRelate来对我的基因型跟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 &

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)

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

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

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)