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 <- "~/lutefisk/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="") )


f6 <- read.table( paste(prefix116, "10-50ns_fluct_fitINT_1.out", sep="") )
f7 <- read.table( paste(prefix129, "10-50ns_fluct_fitINT_1.out", sep="") )
#f6 <- read.table( paste(prefix116, "10-50ns_fluct_fitINT_1.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_EQ_10-50ns_1.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(f4[,2]-f1[,2])

dcut=0.6
d.inds <- which(d>dcut)
#d.inds <- c(d.inds, which(d2>dcut))


#dup=duplicated(d.inds)
#d.inds = d.inds[dup]


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==31 | d.inds==32 | d.inds==34 | d.inds==36 | d.inds==51 | d.inds==53)]
#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]) ##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]
#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(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


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_EQ_10-50ns_2.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(f3[,2]-f1[,2])
d2<-abs(f5[,2]-f1[,2])

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


dup=duplicated(d.inds)
d.inds = d.inds[dup]


t=c(1:length(d.inds))*0
m <- cbind( f1[d.inds,2], f3[d.inds,2], f5[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(f1[,2], typ = "h", sse = sse, xlab = "", ylab = "RMSF (Å)", ylim = c(0,5), xlim=xlim, lwd=3.5, col=bs, axes=F) ## APO

lines(f3[,2], typ="o", col="blue", lwd=1, cex=0.5) ## APO
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(f1[,2], typ = "h", sse = sse, xlab = "", ylab = "RMSF (Å)", ylim = c(0,3.5), xlim=xlim, lwd=4.1, col=bs, axes=F) ## APO

lines(f3[,2], typ="o", col="blue", lwd=1, cex=0.5) ## APO
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()












###############3












pdf("fluct_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(f6[,2]-f7[,2])

d.inds <- which(d>.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( f6[d.inds,2], f7[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(f7[,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(f6[,2], typ="o", col="red", lwd=1, cex=0.5)   ## APO


#slines(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)


dev.off()



## 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)