Post date: Nov 08, 2015 5:38:24 PM
I fit the HMM for Doro's Fst estimates for N1 R x GIS, N1 A x C, and FHA R x GIS. I can do others if we want, but these were the ones in the most recent plot she sent, and are probably most relevant. The plots are largely as expected, high states on LG8 (only) form the R x GIS comparisons, and nothing for A x C (the bump on LG8 wasn't picked up by the HMM). All of the files are in /labs/evolution/projects/timema_speciation/tcristinae/fstdoro/ (the perl script to make the R scripts is in the scripts file). Here is the R code for one of the comparisons:
library(HiddenMarkov)
nl<-42687
nei<-read.table("HudsonFst_s090p001_binnedSNPs_FHA_N1_morphshosts.txt",header=TRUE)
Mstep.normc <- function(x, cond, pm, pn){
nms <- sort(names(pm))
n <- length(x)
m <- ncol(cond$u)
if (all(nms==c("mean", "sd"))){
mean <- pm$mean
#mean <- as.numeric(matrix(x, nrow = 1) %*% cond$u)/apply(cond$u,
#MARGIN = 2, FUN = sum)
sd <- pm$sd
return(list(mean=mean, sd=sd))
}
if (all(nms == "mean")) {
mean <- pm$mean
#mean <- as.numeric(matrix(x, nrow = 1) %*% cond)/apply(cond$u,
#MARGIN = 2, FUN = sum)
return(list(mean = mean))
}
if (all(nms == "sd")) {
sd <- pm$sd
return(list(sd = sd))
}
stop("Invalid specification of parameters")
}
rnormc <- rnorm
dnormc <- dnorm
pnormc <- pnorm
qnormc <- qnorm
lgwin<-which(is.na(nei[,1])==FALSE)
for(parm in 12){
lparm<-log(nei[lgwin,parm]/(1-nei[lgwin,parm]))
lparm[is.finite(lparm)==FALSE]<-log(0.000001/0.999999)
nl<-length(lgwin)
delta<-c(0.05,0.95)
PiA<-matrix(c(0.9,0.1,0.1,0.9),nrow=2,byrow=T)
PiB<-matrix(c(0.8,0.2,0.01,0.99),nrow=2,byrow=T)
Pi<-list(PiA,PiB)
mn<-mean(lparm,na.rm=TRUE)
s<-sd(lparm,na.rm=TRUE)
q90<-quantile(lparm,prob=0.9)
params<-vector("list",2)
fit<-vector("list",2)
for(i in 1:2){
init<-dthmm(lparm,Pi[[i]],delta,"normc",list(mean=c(q90,mn),sd=rep(s,2)),discrete=FALSE)
params[[i]]<-BaumWelch(init,bwcontrol(maxiter = 500, tol = 1e-04, prt = TRUE, posdiff = FALSE,
converge = expression(diff < tol)))
fit[[i]]<-Viterbi(params[[i]])
}
save(list=ls(),file=paste("hmm2s_N1_RxGISpar",parm,".rwsp",sep=""))
mycolors<-c("orange","darkgray")
png(file=paste("hmm2stPlotN1_RxGISpar",parm,".png",sep=""),width=2400,height=800)
eparm<-1/(1+exp(-1 * lparm))
plot(eparm,col=mycolors[fit[[1]]],pch=19,ylim=c(0,1),xlab="snp number",ylab="Fst",cex.lab=1.5,cex.axis=1.2,main="")
lgs<-nei[lgwin,1]
bnds<-which(lgs[1:(nl-1)] != lgs[2:nl]) + 0.5
abline(v=bnds,lty=2)
dev.off()
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=paste("hmm2stSummaryN1_RxGISpar",parm,".txt",sep=""),row.names=F,col.names=F)
cor.test(fit[[1]],fit[[2]])
}
rm(list=ls())