Post date: May 20, 2015 12:0:59 AM
Did this on home laptop:
#17may15
sosopopinfo<-read.csv("sosorum_pops.csv",header=F)
sosogendata<-read.table("filtered_varEsosorumNoHdr.genest",header=F)
g<-(as.matrix(sosogendata))
gmn<-apply(g,1,mean)
nloci<-13999
nind<-258
gmnmat<-matrix(gmn,nrow=nloci,ncol=nind)
gprime<-g-gmnmat
gcovarmat<-matrix(NA,nrow=nind,ncol=nind)
for(i in 1:nind){
for(j in i:nind){
if (i==j){
gcovarmat[i,j]<-cov(gprime[,i],gprime[,j])
}
else{
gcovarmat[i,j]<-cov(gprime[,i],gprime[,j])
gcovarmat[j,i]<-gcovarmat[i,j]
}
}
}
pcgcov<-prcomp(x=gcovarmat,center=TRUE,scale=FALSE)
Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6
Standard deviation 0.02537 0.01593 0.01129 0.01079 0.008915 0.008829
Proportion of Variance 0.06415 0.02530 0.01269 0.01160 0.007920 0.007770
Cumulative Proportion 0.06415 0.08944 0.10213 0.11374 0.121660 0.129420
pops<-as.character(unique(sosopopinfo[,1]))
library(RColorBrewer)
pdf("sosorum_gcov_PCA_17may15.pdf",width=6,height=6)
par(mar=c(4.5,4.5,1.5,1.5))
plot(pcgcov$x[,1],pcgcov$x[,2],type="n",ylab="PC2 (2.5%)",xlab="PC1 (6.4%)")
title(main="A)",adj=0)
myc<-brewer.pal(5,"Accent")
for(i in 5:9){
a<-which(as.character(sosopopinfo[,1])==pops[i])
mn1<-mean(pcgcov$x[a,1])
mn2<-mean(pcgcov$x[a,2])
for(j in 1:length(a)){
lines(c(mn1,pcgcov$x[a[j],1]),c(mn2,pcgcov$x[a[j],2]),col=myc[i-4])
}
}
for(i in 5:9){
a<-which(as.character(sosopopinfo[,1])==pops[i])
mn1<-mean(pcgcov$x[a,1])
mn2<-mean(pcgcov$x[a,2])
text(mn1,mn2,pops[i])
}
legend(0.03,0.035,pops[5:9],col=myc,lty=1,lwd=1.5)
dev.off()
pdf("sosorum_gcov_PCA_17may15All.pdf",width=6,height=6)
par(mar=c(4.5,4.5,1.5,1.5))
plot(pcgcov$x[,1],pcgcov$x[,2],type="n",ylab="PC2 (2.5%)",xlab="PC1 (6.4%)")
title(main="A)",adj=0)
myc<-c(brewer.pal(4,"Accent"),rep("gray",5))
for(i in 1:9){
a<-which(as.character(sosopopinfo[,1])==pops[i])
mn1<-mean(pcgcov$x[a,1])
mn2<-mean(pcgcov$x[a,2])
for(j in 1:length(a)){
lines(c(mn1,pcgcov$x[a[j],1]),c(mn2,pcgcov$x[a[j],2]),col=myc[i])
}
}
for(i in 1:4){## change 9 to 4 to drop wild labels in second plot
a<-which(as.character(sosopopinfo[,1])==pops[i])
mn1<-mean(pcgcov$x[a,1])
mn2<-mean(pcgcov$x[a,2])
text(mn1,mn2,pops[i])
}
legend(-0.044,-0.01,c(pops[1:4],"wild"),col=myc[1:5],lty=1,lwd=1.5,cex=0.8)
dev.off()
## for colors see
help(brewer.pal) ## look at qualitative
colors()
#18may15
pdf("sosorum_gcov_PCA_18may15.pdf",width=6,height=7)
par(mfrow=c(2,1))
par(mar=c(2, 4, 1, 3.5))
plot(pcgcov$x[,1],pcgcov$x[,2],type="n",ylab="PC2 (2.5%)",xlab=NA,cex.lab=0.9,cex.axis=0.8)
myc<-brewer.pal(5,"Set2")
for(i in 5:9){
a<-which(as.character(sosopopinfo[,1])==pops[i])
mn1<-mean(pcgcov$x[a,1])
mn2<-mean(pcgcov$x[a,2])
for(j in 1:length(a)){
lines(c(mn1,pcgcov$x[a[j],1]),c(mn2,pcgcov$x[a[j],2]),col=myc[i-4])
}
}
for(i in 5:9){
a<-which(as.character(sosopopinfo[,1])==pops[i])
mn1<-mean(pcgcov$x[a,1])
mn2<-mean(pcgcov$x[a,2])
text(mn1,mn2,pops[i],cex=0.7)
}
legend(0.04,-0.014,pops[5:9],col=myc,lty=1,lwd=1.5,cex=0.8,bty="n")
par(mar=c(4, 4, 1, 3.5))
plot(pcgcov$x[,1],pcgcov$x[,2],type="n",ylab="PC2 (2.5%)",xlab="PC1 (6.4%)",cex.lab=0.9, cex.axis=0.8)
myc<-c(brewer.pal(4,"Set2"),rep("gray",5))
for(i in 1:9){
a<-which(as.character(sosopopinfo[,1])==pops[i])
mn1<-mean(pcgcov$x[a,1])
mn2<-mean(pcgcov$x[a,2])
for(j in 1:length(a)){
lines(c(mn1,pcgcov$x[a[j],1]),c(mn2,pcgcov$x[a[j],2]),col=myc[i])
}
}
for(i in 1:4){## change 9 to 4 to drop wild labels in second plot
a<-which(as.character(sosopopinfo[,1])==pops[i])
mn1<-mean(pcgcov$x[a,1])
mn2<-mean(pcgcov$x[a,2])
text(mn1,mn2,pops[i],cex=0.7)
}
legend(0.041,-0.01,c(pops[1:4],"wild"),col=myc[1:5],lty=1,lwd=1.5,cex=0.8,bty="n")
dev.off()
#gst
P<-matrix(NA,nrow=9,ncol=13999)
for(i in 1:9){
a<-which(as.character(sosopopinfo[,1])==pops[i])
P[i,]<-apply(g[,a],1,mean)/2
}
Hs<-2 * P * (1-P)
fst<-matrix(0,nrow=9,ncol=9)
for(i in 1:8){
for(j in i:9){
Pbar<-(P[i,] + P[j,])/2
Ht<-2 * Pbar * (1 - Pbar)
mHs<-(Hs[i,] + Hs[j,])/2
fst[i,j]<-mean(Ht-mHs)/mean(Ht)
fst[j,i]<-fst[i,j]
}
}
colnames(fst)<-pops
rownames(fst)<-pops
C-EL C-F1 C-F2 C-F3 W-CS W-EL
C-EL 0.000000000 0.002882524 0.002405027 0.008426096 0.008594560 0.005325068
C-F1 0.002882524 0.000000000 0.002067114 0.007811131 0.007139952 0.003238945
C-F2 0.002405027 0.002067114 0.000000000 0.007475998 0.008089740 0.004206818
C-F3 0.008426096 0.007811131 0.007475998 0.000000000 0.013978102 0.010039431
W-CS 0.008594560 0.007139952 0.008089740 0.013978102 0.000000000 0.007407459
W-EL 0.005325068 0.003238945 0.004206818 0.010039431 0.007407459 0.000000000
W-PO 0.007116060 0.004745318 0.005935812 0.011851950 0.008684956 0.003785814
W-SG 0.013513854 0.011243964 0.012472360 0.019100526 0.015521626 0.010559405
W-UP 0.018546821 0.015931961 0.017363999 0.023294731 0.019398915 0.014881842
W-PO W-SG W-UP
C-EL 0.007116060 0.01351385 0.01854682
C-F1 0.004745318 0.01124396 0.01593196
C-F2 0.005935812 0.01247236 0.01736400
C-F3 0.011851950 0.01910053 0.02329473
W-CS 0.008684956 0.01552163 0.01939892
W-EL 0.003785814 0.01055940 0.01488184
W-PO 0.000000000 0.01167396 0.01572172
W-SG 0.011673960 0.00000000 0.02290725
W-UP 0.015721718 0.02290725 0.00000000
## 95% CIs, lower-triangle = lower bound, upper triangle = upper bound
## bootstrapping to get CIs
fstCI<-matrix(0,nrow=9,ncol=9)
for(i in 1:8){
for(j in i:9){
Pbar<-(P[i,] + P[j,])/2
Ht<-2 * Pbar * (1 - Pbar)
mHs<-(Hs[i,] + Hs[j,])/2
bsGst<-rep(NA,1000)
for(k in 1:1000){
a<-sample(1:13999,13999,replace=TRUE)
bsGst[k]<-mean(Ht[a]-mHs[a])/mean(Ht[a])
}
fstCI[j,i]<-quantile(bsGst,probs=0.025)
fstCI[i,j]<-quantile(bsGst,probs=0.975)
}
}
colnames(fstCI)<-pops
rownames(fstCI)<-pops
C-EL C-F1 C-F2 C-F3 W-CS W-EL
C-EL 0.000000000 0.002956499 0.002472242 0.008693861 0.008810619 0.005462856
C-F1 0.002808502 0.000000000 0.002123471 0.008020406 0.007331873 0.003328373
C-F2 0.002339101 0.002012070 0.000000000 0.007699653 0.008305288 0.004313158
C-F3 0.008185666 0.007587838 0.007277009 0.000000000 0.014365632 0.010308079
W-CS 0.008361492 0.006941805 0.007876274 0.013576393 0.000000000 0.007625392
W-EL 0.005186631 0.003150116 0.004093806 0.009762587 0.007202771 0.000000000
W-PO 0.006907090 0.004612458 0.005783220 0.011528809 0.008424840 0.003685410
W-SG 0.013159005 0.010969716 0.012157092 0.018604101 0.015123570 0.010301390
W-UP 0.018061518 0.015485834 0.016880054 0.022645876 0.018838983 0.014457363
W-PO W-SG W-UP
C-EL 0.007300932 0.01389504 0.01901187
C-F1 0.004877172 0.01153466 0.01640919
C-F2 0.006097424 0.01281119 0.01785195
C-F3 0.012168332 0.01961758 0.02390328
W-CS 0.008911266 0.01595018 0.01993171
W-EL 0.003891808 0.01085709 0.01528939
W-PO 0.000000000 0.01200949 0.01619164
W-SG 0.011359641 0.00000000 0.02358603
W-UP 0.015245098 0.02229309 0.00000000