Post date: Jun 17, 2015 4:36:57 PM
Here is a PCA plot for the seed beetles based on modal genotypes and a genetic covariance matrix. This is part of the evidence that L2 and L3 samples come from two different sub-lines. Here is the code:
estpost_popmod -p g -s 0 -w 0 -o g_L1-F91.txt popaf_L1-F91_chain0.hdf5 popaf_L1-F91_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L1-F100.txt popaf_L1-F100_chain0.hdf5 popaf_L1-F100_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L1R-F35.txt popaf_L1R-F35_chain0.hdf5 popaf_L1R-F35_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L1R-F46.txt popaf_L1R-F46_chain0.hdf5 popaf_L1R-F46_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L2-F78.txt popaf_L2-F78_chain0.hdf5 popaf_L2-F78_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L2-F87.txt popaf_L2-F87_chain0.hdf5 popaf_L2-F87_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L2R-F35.txt popaf_L2R-F35_chain0.hdf5 popaf_L2R-F35_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L2R-F45.txt popaf_L2R-F45_chain0.hdf5 popaf_L2R-F45_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L3-F76.txt popaf_L3-F76_chain0.hdf5 popaf_L3-F76_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L3-F85.txt popaf_L3-F85_chain0.hdf5 popaf_L3-F85_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L3R-F35.txt popaf_L3R-F35_chain0.hdf5 popaf_L3R-F35_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_L3R-F45.txt popaf_L3R-F45_chain0.hdf5 popaf_L3R-F45_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_M-13.txt popaf_M-13_chain0.hdf5 popaf_M-13_chain1.hdf5
estpost_popmod -p g -s 0 -w 0 -o g_M-14.txt popaf_M-14_chain0.hdf5 popaf_M-14_chain1.hdf5
cat g_* > genotypes.csv
### R code ####
N<-576
L<-51371
G<-matrix(scan("genotypes.csv",n=N*L,sep=","),nrow=N,ncol=L,byrow=TRUE)
Ghat<-round(G,0)
GhatCent<-Ghat
for(i in 1:L){
GhatCent[,i]<-GhatCent[,i]-mean(GhatCent[,i])
}
GhatCov<-matrix(NA,nrow=N,ncol=N)
for(i in 1:N){
for(j in i:N){
GhatCov[i,j]<-cov(GhatCent[i,],GhatCent[j,])
GhatCov[j,i]<-GhatCov[i,j]
}}
ids<-c(rep(c("L1-F100","L1-F91","L1R-F35","L1R-F46","L2-F78","L2-F87","L2R-F35","L2R-F45","L3-F76","L3-F85","L3R-F35","L3R-F45"),each=40),rep(c("M-13","M-14"),each=48))
## pca
pcout<-prcomp(GhatCov,center=FALSE,stand=FALSE)
mn1<-tapply(X=pcout$x[,1],INDEX=ids,mean)
mn2<-tapply(X=pcout$x[,2],INDEX=ids,mean)
mn3<-tapply(X=pcout$x[,3],INDEX=ids,mean)
uids<-unique(ids)
P<-length(uids)
myc<-rep(c("darkblue","lightblue","darkgreen","limegreen","red","orange","black"),each=2)
myl<-rep(c(1,2),7)
myl[1]<-2
myl[2]<-1
pdf("cmacIndPca.pdf",width=6,height=12)
par(mfrow=c(2,1))
par(mar=c(4,5,1.5,0.5))
plot(pcout$x[,1],pcout$x[,2],type='n',xlab="PC1 (42%)",ylab="PC2 (20%)",cex.lab=1.4)
for(j in 1:P){
x<-which(ids==uids[j])
n<-length(x)
for(a in 1:n){
lines(c(mn1[j],pcout$x[x[a],1]),c(mn2[j],pcout$x[x[a],2]),lwd=0.7,lty=myl[j],col=myc[j])
}
}
legend(-0.9,0.48,uids,lty=myl,col=myc,cex=0.5)
plot(pcout$x[,1],pcout$x[,3],type='n',xlab="PC1 (42%)",ylab="PC3 (13%)",cex.lab=1.4)
for(j in 1:P){
x<-which(ids==uids[j])
n<-length(x)
for(a in 1:n){
lines(c(mn1[j],pcout$x[x[a],1]),c(mn3[j],pcout$x[x[a],3]),lwd=0.7,lty=myl[j],col=myc[j])
}
}
dev.off()