a=c(0,1,2,0,0,1,2,1,0,1,2,2) b=matrix(a,3,4,byrow=T) b snps=b snpMean= apply(snps,2,mean) #get mean for each snp snpMean snps=t(snps)-snpMean #operation on matrix and vector goes in direction of column snps K=crossprod(snps, snps) K library(compiler) #required for cmpfun source("http://www.zzlab.net/GAPIT/gapit_functions.txt") library(gplots) myGD=read.table(file="http://zzlab.net/GAPIT/data/mdp_numeric.txt",head=T) taxa=myGD[,1] favorite=c("33-16", "38-11", "B73", "B73HTRHM", "CM37", "CML333", "MO17", "YU796NS") B73=taxa%in%c("B73") #taxa=taxa[!B73] index=taxa%in%favorite snps=myGD[,-1] #snps=snps[!B73,] #K=GAPIT.kinship.loiselle(t(myGD[,-1]), method="additive", use="all") #K[index,index] K1=GAPIT.kinship.VanRaden(snps) K1[index,index] K2=GAPIT.kinship.ZHANG(snps) K2[index,index] #heatmap.2(K1, cexRow =.2, cexCol = 0.2, col=rev(heat.colors(256)), scale="none", symkey=FALSE, trace="none") #quartz() #heatmap.2(K2, cexRow =.2, cexCol = 0.2, col=rev(heat.colors(256)), scale="none", symkey=FALSE, trace="none") n=nrow(snps) ind.a=seq(1:(n*n)) i =1:n j=(i-1)*n ind.d=i+j par(mfrow=c(1,3)) plot(K2[ind.a],K1[ind.a],main="All elements",xlab="Zhang",ylab="VanRaden") lines(K2[ind.d],K1[ind.d],main="All elements",xlab="Zhang",ylab="VanRaden",col="red",type="p") plot(K2[ind.d],K1[ind.d],main="Diagonals",xlab="Zhang",ylab="VanRaden") plot(K2[-ind.d],K1[-ind.d],main="Off diag",xlab="Zhang",ylab="VanRaden")