library(bio3d)
library(Hmisc)


prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/fluctuations/"
prefix129 <- "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/fluctuations/"
prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/fluctuations/"


f1 <- read.table( paste(prefix116, "10-50ns_fluct_fitEQ_1.out", sep="") )
f2 <- read.table( paste(prefix129, "10-50ns_fluct_fitEQ_1.out", sep="") )
f3 <- read.table( paste(prefix129, "10-50ns_fluct_fitEQ_2.out", sep="") )
f4 <- read.table( paste(prefix171, "10-50ns_fluct_fitEQ_1.out", sep="") )
f5 <- read.table( paste(prefix171, "10-50ns_fluct_fitEQ_2.out", sep="") )


#f3 <- read.table( paste(dir, "20-60ns_fluct_fitInt_1.out", sep="") )
#f4 <- read.table( paste(dir, "20-60ns_fluct_fitInt_2.out", sep="") )
#f5 <- read.table( paste(dir, "20-60ns_fluct_fitApi_1.out", sep="") )
#f6 <- read.table( paste(dir, "20-60ns_fluct_fitApi_2.out", sep="") )




pdb <- read.pdb("/net/lutefisk/slars/groel_md/bio3d_groel/raw_pdbs_new/split_chain/1XCK_A.pdb")
seq <- seq.pdb(pdb)
s <- array(seq)
sse <- dssp(pdb)


bs <- c(1:524)*0

bs[30] <- 2 ## Leu31
bs[31] <- 2
bs[32] <- 2
bs[52] <- 2

##bs[86] <- 2 ## Asp87
bs[86:90] <- 2 ## 

#bs[149] <- 2
#bs[150] <- 2

#bs[397] <- 2
#bs[414] <- 2
bs[453] <- 2
bs[478:480] <- 2
bs[492] <- 2
bs[494] <- 2

bs <- bs+1


cex=.7
cex.mtext=1.2
sse.cex=.5


pdf("fluct_all_10-50ns.pdf", w=9, h=5)
par(mfrow=c(2,1), mar=c(2.5,2.5,1.5,1), mgp=c(1.4,0.5,0))


## EQ domain
d<-abs(f2[,2]-f1[,2])
d2<-abs(f2[,2]-f3[,2])

d.inds <- which(d>.4)
d.inds = c(d.inds, which(d2>.4))


d.inds = d.inds[-which(d.inds<10 | (d.inds>36&d.inds<50) | (d.inds>136&d.inds<150))]
d.inds = d.inds[-which(d.inds==15 | d.inds==31 | d.inds==32 | d.inds==34 | d.inds==28 | d.inds==36 | d.inds==51 | d.inds==53 | d.inds==85 | d.inds==135)]
d.inds = d.inds[-which(d.inds==475 | d.inds==479 | d.inds==473)]
t=c(1:length(d.inds))*0
m <- cbind( f1[d.inds,2], f2[d.inds,2], f3[d.inds,2] )
maxcol <- max.col(m)
t[which(maxcol==1)]=m[which(maxcol==1),1]
t[which(maxcol==2)]=m[which(maxcol==2),2]
t[which(maxcol==3)]=m[which(maxcol==3),3]


xlim = c(1,135)
plot.bio3d(f2[,2], typ = "h", sse = sse, xlab = "", ylab = "RMSF (Å)", ylim = c(0,5), xlim=xlim, lwd=3.5, col=bs, axes=F) ## ATP
lines(f3[,2], typ="o", col="blue", lwd=1, cex=0.5) ## APO

lines(f1[,2], type='o', col="red", lwd=1, cex=0.5) ## APO
 
lines(f4[,2], type='o', col="grey70", lwd=1, cex=0.5) ## ATP & subs
lines(f5[,2], typ="o", col="green", lwd=1, cex=0.5) ## APO

abline(v=array(sse$helix$start), col="grey60", lty=2)
abline(v=array(sse$sheet$start), col="grey60", lty=2)
axis(1, at=sse$helix$start, labels=sse$helix$start)
axis(2)
box()
helices.inds=which(sse$helix$start>xlim[1] & sse$helix$start<xlim[2])
mtext(capitalize(letters[helices.inds]), side=3, line=0.3, at=(sse$helix$start[helices.inds]+sse$helix$end[helices.inds])/2, cex=sse.cex)


offset=1
mtext("A", side=3, at=-2, cex=cex.mtext, line=0.25)
text(d.inds, t+offset, paste( aa123(s[d.inds]), d.inds+1), srt=90, cex=cex)

xlim=c(410,524)
plot.bio3d(f2[,2], typ = "h", sse = sse, xlab = "", ylab = "RMSF (Å)", ylim = c(0,3.5), xlim=xlim, lwd=4.1, col=bs, axes=F) ## ATP
lines(f3[,2], typ="o", col="blue", lwd=1, cex=0.5) ## APO

lines(f1[,2], typ="o", col="red", lwd=1, cex=0.5)   ## APO


lines(f4[,2], type='o', col="grey70", lwd=1, cex=0.5)  ## ATP & subs
lines(f5[,2], typ="o", col="green", lwd=1, cex=0.5) ## APO

abline(v=array(sse$helix$start), col="grey60", lty=2)
abline(v=array(sse$sheet$start), col="grey60", lty=2)
axis(1, at=sse$helix$start, labels=sse$helix$start)
axis(2)
box()
helices.inds=which(sse$helix$start>xlim[1] & sse$helix$start<xlim[2])
mtext(capitalize(letters[helices.inds]), side=3, line=0.3, at=(sse$helix$start[helices.inds]+sse$helix$end[helices.inds])/2, cex=sse.cex)

offset=.7
mtext("B", side=3, at=xlim[1]-3, cex=cex.mtext, line=0.25)
text(d.inds, t+offset, paste( aa123(s[d.inds]), d.inds+1), srt=90, cex=cex)




dev.off()








pdf("fluct_ap_int_10-50ns.pdf", w=9, h=5)
par(mfrow=c(2,1), mar=c(2.5,2.5,1.5,1), mgp=c(1.4,0.5,0))


## Intermediate
d<-abs(f2[,2]-f1[,2])
d2<-abs(f2[,2]-f3[,2])

d.inds <- which(d>.4)
d.inds = c(d.inds, which(d2>.4))

d.inds = d.inds[-which((d.inds>195&d.inds<350) | d.inds==169 | d.inds==374)]
t=c(1:length(d.inds))*0
m <- cbind( f3[d.inds,2], f4[d.inds,2] )
maxcol <- max.col(m)
t[which(maxcol==1)]=m[which(maxcol==1),1]
t[which(maxcol==2)]=m[which(maxcol==2),2]

xlim=c(135,409)
plot.bio3d(f2[,2], typ = "h", sse = sse, xlab = "", ylab = "RMSF (Å)",
           ylim = c(0,5), xlim=xlim, lwd=1.5, col=bs, axes=F)

#axis(1, at=array(sse$helix$start), labels=array(sse$helix$start))
lines(f3[,2], typ="o", col="blue", lwd=1, cex=0.5) ## APO

lines(f1[,2], typ="o", col="red", lwd=1, cex=0.5)   ## APO


lines(f4[,2], type='o', col="grey70", lwd=1, cex=0.5)  ## ATP & subs
lines(f5[,2], typ="o", col="green", lwd=1, cex=0.5) ## APO


abline(v=array(sse$helix$start), col="grey60", lty=2)
abline(v=array(sse$sheet$start), col="grey60", lty=2)

axis(1, at=sse$helix$start, labels=sse$helix$start)
axis(2)
box()

helices.inds=which(sse$helix$start>xlim[1] & sse$helix$start<xlim[2])
mtext(capitalize(letters[helices.inds]), side=3, line=0.3, at=(sse$helix$start[helices.inds]+sse$helix$end[helices.inds])/2, cex=sse.cex)


offset=1
mtext("C", side=3, at=xlim[1]-6, cex=cex.mtext, line=0.25)
text(d.inds, t+offset, paste( aa123(s[d.inds]), d.inds+1), srt=90, cex=cex)



## Apical
d<-abs(f2[,2]-f1[,2])
d2<-abs(f2[,2]-f3[,2])

d.inds <- which(d>.4)
d.inds = c(d.inds, which(d2>.4))

t=c(1:length(d.inds))*0
m <- cbind( f1[d.inds,2], f2[d.inds,2] )
maxcol <- max.col(m)
t[which(maxcol==1)]=m[which(maxcol==1),1]
t[which(maxcol==2)]=m[which(maxcol==2),2]

 xlim=c(190,375)
plot.bio3d(f2[,2], typ = "h", sse = sse, xlab = "Residue No.", ylab = "RMSF (Å)",
           ylim = c(0,5), xlim=xlim, lwd=1.5, col=bs, axes=F)
lines(f3[,2], typ="o", col="blue", lwd=1, cex=0.5) ## APO

lines(f1[,2], typ="o", col="red", lwd=1, cex=0.5)   ## APO


lines(f4[,2], type='o', col="grey70", lwd=1, cex=0.5)  ## ATP & subs
lines(f5[,2], typ="o", col="green", lwd=1, cex=0.5) ## APO

abline(v=array(sse$helix$start), col="grey60", lty=2)
abline(v=array(sse$sheet$start), col="grey60", lty=2)

axis(1, at=sse$helix$start, labels=sse$helix$start)
axis(2)
box()

helices.inds=which(sse$helix$start>xlim[1] & sse$helix$start<xlim[2])
mtext(capitalize(letters[helices.inds]), side=3, line=0.3, at=(sse$helix$start[helices.inds]+sse$helix$end[helices.inds])/2, cex=sse.cex)


offset=1
mtext("D", side=3, at=xlim[1]-5, cex=cex.mtext, line=0.25)
text(d.inds, t+offset, paste( aa123(s[d.inds]), d.inds+1), srt=90, cex=cex)





dev.off()




## Whole
#plot.bio3d(f1[,2], typ = "h", sse = sse, xlab = "Residue No.", ylab = "RMSF (Å)", ylim = c(0,15), lwd=3.5, col=bs, axes=T)
#lines(f2[,2], type='l', col="red", lwd=1, cex=0.5)
#abline(v=array(sse$helix$start), col="grey60", lty=2)
#abline(v=array(sse$sheet$start), col="grey60", lty=2)