Post date: Nov 02, 2014 10:43:42 PM
#Tried to use R script I used for other three taxa, but found 5 loci with 1 allele only, so the afreq files wouldn't read into R like I thought they would. So, to bypass this issue, we used cut to extract the first 2 columns only of each afreq mcmc file -- this is done in the terminal not in R:
cut -f 1,2 -d " " 25viii14_pop_snpfile_S-C-CS-PAc0afreq.txt > 2col_25viii14_pop_snpfile_S-C-CS-PAc0afreq.txt
cut -f 1,2 -d " " 25viii14_pop_snpfile_S-C-CS-PAc1afreq.txt > 2col_25viii14_pop_snpfile_S-C-CS-PAc1afreq.txt
cut -f 1,2 -d " " 25viii14_pop_snpfile_S-C-CS-R1c0afreq.txt > 2col_25viii14_pop_snpfile_S-C-CS-R1c0afreq.txt &
cut -f 1,2 -d " " 25viii14_pop_snpfile_S-C-CS-R1c1afreq.txt > 2col_25viii14_pop_snpfile_S-C-CS-R1c1afreq.txt &
cut -f 1,2 -d " " 25viii14_pop_snpfile_S-C-CS-R2c0afreq.txt > 2col_25viii14_pop_snpfile_S-C-CS-R2c0afreq.txt &
cut -f 1,2 -d " " 25viii14_pop_snpfile_S-C-CS-R2c1afreq.txt > 2col_25viii14_pop_snpfile_S-C-CS-R2c1afreq.txt &
cut -f 1,2 -d " " 3ix14_pop_snpfile_S-C-CS-SIc0afreq.txt > 2col_3ix14_pop_snpfile_S-C-CS-SIc0afreq.txt &
cut -f 1,2 -d " " 3ix14_pop_snpfile_S-C-CS-SIc1afreq.txt > 2col_3ix14_pop_snpfile_S-C-CS-SIc1afreq.txt &
#Then, I sourced code through R and ran in it in the background in greenhouse overnight to try to get pairwise Gst. Check greenhouse tomorrow morning to see how it did!
R CMD BATCH SpcalcFst_source.R &
#Here's the source code:
## concatenate files, and store mcmc samples in matrixes
nloci<-191678 ##Apparently this is the number of loci, not 191679. Fixed this error and started running again 3Nov14.
nsteps<-2000 * 2
afPA.mat<-matrix(c(afPAc0[,2],afPAc1[,2]),nrow=nsteps,ncol=nloci,byrow=T)
afR1.mat<-matrix(c(afR1c0[,2],afR1c1[,2]),nrow=nsteps,ncol=nloci,byrow=T)
afR2.mat<-matrix(c(afR2c0[,2],afR2c1[,2]),nrow=nsteps,ncol=nloci,byrow=T)
afSI.mat<-matrix(c(afSIc0[,2],afSIc1[,2]),nrow=nsteps,ncol=nloci,byrow=T)
##4Nov14 -- R got stuck here. Need to take a different approach. When get a chance, copy to my desktop and run at USU.
## remove single chain data objects
rm(afPAc0,afPAc1)
rm(afR1c0,afR1c1)
rm(afR2c0,afR2c1)
rm(afSIc0,afSIc1)
save(file="Sp_pairwisefst.Rwsp",list=ls())
af.mat<-list(afPA.mat,afR1.mat,afR2.mat,afSI.mat)
# calculate Fst = (Ht - Hs) / Ht
# first cacluate Hs p^2 + (1-p)^2
npop<-4
hs.mat<-af.mat
for(i in 1:npop){
hs.mat[[i]]<-af.mat[[i]]^2 + (1-af.mat[[i]])^2
}
ht<-matrix(NA,nrow=nsteps,ncol=nloci)
#ncomp = (n * n-1)/2
ncomp<-(4 * 3)/2
## must type ht, ncomp times
ht.mat<-list(ht,ht,ht,ht,ht,ht)
x<-1
for(i in 1:(npop-1)){
for(j in (i+1):npop){
atotal<-(af.mat[[i]] + af.mat[[j]])/2
ht.mat[[x]]<-atotal^2 + (1-atotal)^2
x<-x+1
}
}
save(file="Sp_pairwisefst.Rwsp",list=ls())
fst.mat<-list(ht,ht,ht,ht,ht,ht)
x<-1
for(i in 1:(npop-1)){
for(j in (i+1):npop){
hsmean<-(hs.mat[[i]] + hs.mat[[j]])/2
fst.mat[[x]]<-((1-ht.mat[[x]]) - (1-hsmean))/(1-ht.mat[[x]])
x<-x+1
}
}
# point estimates and ci
myest<-apply(mnf.mat,2,quantile,probs=c(0.5,0.025,0.975))
write.csv(myest,file="pairwisefst_Sp-CS.csv")
save(file="Sp_pairwisefst.Rwsp",list=ls())