Post date: Oct 13, 2014 5:26:45 PM
I summarized the HMM results for each population pair using the plotHmm*R scripts in /labs/evolution/projects/timema_radiation/popgen/fsts/. I ran all of these with the wrap_qsub_rc_hmm.pl script. Here is the script:
mycolors<-c("orange","darkgray","lightblue")
###########################################
load("hmmVP_CxVP_Q.rwsp")
png(file="hmmPlotVP_CxVP_Q.png",width=1500,height=500)
plot(kept,(1/(1+exp(-1 * lfst))),col=mycolors[fit[[1]]],pch=19,xlab="snp number",ylab=expression(paste(F[ST])),cex.lab=1.5,cex.axis=1.2,main="")
dev.off()
## first row gives mean Fst, next is prop in each, next 3 are TPM
out<-round(rbind(1/(1+ exp(-1*params[[1]]$pm$mean)),table(fit[[1]])/length(fit[[1]]),params[[1]]$Pi),4)
write.table(out,file="hmmSummaryVP_CxVP_Q.txt",row.names=F,col.names=F)
cor.test(fit[[1]],fit[[2]])
rm(list=ls())
The results are in /labs/evolution/projects/timema_radiation/popgen/fsts/ (at least for now). There is one *png and one *txt file per population pair. The *png files show Fst along the genome (with low-coverage SNPs dropped) with different colors for the three HMM states (1=orange, 2=gray, 3=lightblue). The *txt files have five lines. The first line gives mean Fst for each state (in order, 1, 2, 3). The second line shows the proportion of loci in each state. Lines 3-5 give the transition probability matrix. That is the probability of transition from one state to another state between SNPs (in order from 1-3); note that this can be useful for thinking about the amount of structure in Fst along chromosomes (lots of structure would mean high transition probabilities along the diagonal such that you tend to stay in the same state).
In every case state 1 is the high Fst state (with Fst quite a bit higher than the others). States 2 and 3 are more or less similar depending on the comparison and while state 2 is usually the intermediate Fst state, this is sometimes state 3 instead. Thus, for the moment, we might want to focus on state 1 (high Fst=orange) vs. the other two states together (I can re-color them like this if we want). Just like before we see relatively rapid transitions among states, with more or less structure depending on the comparison (and somewhat it seems on the mean Fst).
HFTP_CxHFRS_Q doesn't seem to have converged and it only used one of the three states. It has a much lower correlation between independent runs (~0.3 vs ~1.0) than other population pairs. I am re-making a plot based on the second run to see what I get.