library(haplo.stats) hgdpURL <- "http://www.biostat.umn.edu/~cavanr/HGDP_AKT1.txt" hgdp <- read.delim(file=hgdpURL) attach(hgdp) geno=cbind(substr(AKT1.C6024T,1,1), substr(AKT1.C6024T,2,2), substr(AKT1.C0756A,1,1), substr(AKT1.C0756A,2,2), substr(AKT1.G2375A,1,1), substr(AKT1.G2375A,2,2), substr(AKT1.G2347T,1,1), substr(AKT1.G2347T,2,2)) SNPnames <- c("C6024T","C0756A","G2375A", "G2347T") i=1 id=levels(Population)[i] dat=geno[Population==id,] h1=haplo.em(dat, locus.label=SNPnames, control=haplo.em.control(min.posterior=1.0e-4)) h1 names(h1) h1$hap.prob h1$haplotype maxNoHap=40 hapProb=matrix(NA, 52, maxNoHap) hapMat=vector("list",52) for(i in 1:52){ id=levels(Population)[i] dat=geno[Population==id,] h1=haplo.em(dat, locus.label=SNPnames, control=haplo.em.control(min.posterior=1.0e-4)) nh=dim(h1$haplotype)[1] hapProb[i,1:nh]=h1$hap.prob hapMat[[i]]=h1$haplotype } dat=numeric(0) for(i in 1:52){ nh=dim(hapMat[[i]])[1] for(j in 1:nh) dat=c(dat,hapMat[[i]][j,]) } length(dat)/4 ddat1=matrix(dat, byrow=T, ncol=4) dat2=rep(NA, 155) for(i in 1:155) dat2[i]=paste(ddat1[i,1],ddat1[i,2],ddat1[i,3],ddat1[i,4],sep="") table(dat2) hapProb[1:5,1:10] t(hapProb[1:5,1:10]) c(t(hapProb[1:5,1:10])) hapProbPop=c(t(hapProb))[!is.na(c(t(hapProb)))] length(hapProbPop) i=1 names(table(dat2))[i] hid=names(table(dat2))[i] hapProbPop[dat2==hid] par(mfrow=c(3,3)) for(i in 1:9){ hid=names(table(dat2))[i] hist(hapProbPop[dat2==hid]) } par(mar=c(2,2,2,1)) for(i in 1:9){ hid=names(table(dat2))[i] hist(hapProbPop[dat2==hid], xlab="", main=paste("hap:", hid)) }