library(bio3d)

pdb1 <- read.pdb("1XCK_reference.pdb")
ca.inds <- atom.select(pdb1, "///9:515///CA/")

prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/"
prefix129 <- "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/"
prefix117 <- "/net/gulrotkake/slars/groel_md/1XCK/117_1XCK_hepta/results/"
prefix128 <- "/net/gulrotkake/slars/groel_md/1XCK/128_1XCK_hepta_MGATP/results/"
prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/"

sim116 <- NULL
sim129 <- NULL
sim117 <- NULL
sim128 <- NULL
sim171 <- NULL

## monomers
for ( i in 1:7 ) {
  tmprmsd <- read.table(paste(prefix116, "rmsd_monomer/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim116$monomer$rmsd=cbind(sim116$monomer$rmsd, tmprmsd[,"V2"])
}

for ( i in 1:14 ) {
  tmprmsd <- read.table(paste(prefix129, "rmsd_monomer/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim129$monomer$rmsd=cbind(sim129$monomer$rmsd, tmprmsd[,"V2"])
}

for ( i in 1:7 ) {
  tmprmsd <- read.table(paste(prefix171, "rmsd_monomer/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim171$monomer$rmsd=cbind(sim171$monomer$rmsd, tmprmsd[,"V2"])
}

## dimers
for ( i in 1:7 ) {
  tmprmsd <- read.table(paste(prefix116, "rmsd_dimers/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim116$dimer$rmsd=cbind(sim116$dimer$rmsd, tmprmsd[,"V2"])
}

for ( i in 1:7 ) {
  tmprmsd <- read.table(paste(prefix129, "rmsd_dimers/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim129$dimer$rmsd=cbind(sim129$dimer$rmsd, tmprmsd[,"V2"])
}

for ( i in 1:7 ) {
  tmprmsd <- read.table(paste(prefix171, "rmsd_dimers/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim171$dimer$rmsd=cbind(sim171$dimer$rmsd, tmprmsd[,"V2"])
}

## trimers
for ( i in c(1,2,3,4,5,7) ) {
  tmprmsd <- read.table(paste(prefix116, "rmsd_trimer/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim116$trimer$rmsd=cbind(sim116$trimer$rmsd, tmprmsd[,"V2"])
}

for ( i in c(1,2,3,4,5,7) ) {
  tmprmsd <- read.table(paste(prefix129, "rmsd_trimer/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim129$trimer$rmsd=cbind(sim129$trimer$rmsd, tmprmsd[,"V2"])
}


## heptamers
for ( i in 1:2 ) {
  tmprmsd <- read.table(paste(prefix116, "rmsd_heptamer/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim116$heptamer$rmsd=cbind(sim116$heptamer$rmsd, tmprmsd[,"V2"])
}

for ( i in 1:2 ) {
  tmprmsd <- read.table(paste(prefix129, "rmsd_heptamer/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim129$heptamer$rmsd=cbind(sim129$heptamer$rmsd, tmprmsd[,"V2"])
}

for ( i in 1:2 ) {
  tmprmsd <- read.table(paste(prefix171, "rmsd_heptamer/50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim171$heptamer$rmsd=cbind(sim171$heptamer$rmsd, tmprmsd[,"V2"])
}

## complex
sim116$complex$rmsd=read.table(paste(prefix116, "rmsd_complex/rmsd.dat", sep=""))
sim129$complex$rmsd=read.table(paste(prefix129, "rmsd_complex/rmsd.dat", sep=""))
#sim171$complex$rmsd=read.table(paste(prefix171, "rmsd_complex/rmsd.dat", sep=""))


## calculating means
sim116$monomer$rmsd.mean <- rowMeans(sim116$monomer$rmsd)
sim129$monomer$rmsd.mean.atp <- rowMeans(sim129$monomer$rmsd[,1:7])
sim129$monomer$rmsd.mean.free <- rowMeans(sim129$monomer$rmsd[,8:14])
#sim171$monomer$rmsd.mean.atp <- rowMeans(sim171$monomer$rmsd[,1:7])
#sim171$monomer$rmsd.mean.free <- rowMeans(sim171$monomer$rmsd[,8:14])

sim116$dimer$rmsd.mean <- rowMeans(sim116$dimer$rmsd)
sim129$dimer$rmsd.mean.atp <- rowMeans(sim129$dimer$rmsd[,1:7])
#sim129$dimer$rmsd.free.mean <- rowMeans(sim129$dimer$rmsd[,8:14])
#sim171$dimer$rmsd.mean.atp <- rowMeans(sim171$dimer$rmsd[,1:7])
#sim171$dimer$rmsd.free.mean <- rowMeans(sim171$dimer$rmsd[,8:14])

sim116$trimer$rmsd.mean <- rowMeans(sim116$trimer$rmsd)
sim129$trimer$rmsd.mean.atp <- rowMeans(sim129$trimer$rmsd[,1:6])
#sim171$trimer$rmsd.mean.atp <- rowMeans(sim171$trimer$rmsd[,1:6])


## calculating means (with inds)
#inds <- seq(1,5000,20)
#sim116$monomer$rmsd.mean <- rowMeans(sim116$monomer$rmsd[inds,])
#sim129$monomer$rmsd.mean.atp <- rowMeans(sim129$monomer$rmsd[inds,1:7])
#sim129$monomer$rmsd.mean.free <- rowMeans(sim129$monomer$rmsd[inds,8:14])

pdf("rmsd_all.pdf", w=12, h=16)
par(mfcol=c(5,1), mar=c(5,4,4,2), mgp=c(3,1,0))


## plotting monomers
plot(sim116$monomer$rmsd[,1], type='l', col="grey80", ylim=c(1, 4.5),
     ylab="RMSD (Å)", xlab="Time (ns)", axes=FALSE )
  

for ( i in 2:7 ) {
  lines(sim116$monomer$rmsd[,i], lty=2, col="grey80")
  lines(sim129$monomer$rmsd[,i], lty=2, col="grey80")
}

lines(sim116$monomer$rmsd.mean, lty=1, col="blue")
lines(sim129$monomer$rmsd.mean.atp, lty=1, col="red")
lines(sim129$monomer$rmsd.mean.free, lty=1, col="orange")
#lines(sim171$monomer$rmsd.mean.atp, lty=2, col="black")
#lines(sim171$monomer$rmsd.mean.free, lty=2, col="green")

axis(2)
at=c(0,200,400,600,800,1000)
axis(1, at=at, labels=c(0,10,20,30,40,50))
box()
legend("topright", legend=c("GroEL (ATP-free)", "GroEL-ATP (cis)", "GroEL-ATP (trans)"), col=c("blue", "red", "orange"), lty=1, ncol=3)
mtext("(A) Monomers", side=3, line=0.7, at=0, cex=1.1, adj=0)



## plotting dimers
plot(sim116$dimer$rmsd[,1], type='l', col="grey80", ylim=c(1, 4.5),
     ylab="RMSD (Å)", xlab="Time (ns)", axes=FALSE )
  

for ( i in 2:7 ) {
  lines(sim116$dimer$rmsd[,i], lty=2, col="grey80")
  lines(sim129$dimer$rmsd[,i], lty=2, col="grey80")
}

#for ( i in 1:7 ) {
#  lines(sim171$dimer$rmsd[,i], lty=2, col="grey80")
#}


lines(sim116$dimer$rmsd.mean, lty=1, col="blue")
lines(sim129$dimer$rmsd.mean.atp, lty=1, col="red")
#lines(sim129$dimer$rmsd.mean.free, lty=2, col="orange")
#lines(sim171$dimer$rmsd.mean.atp, lty=2, col="black")
#lines(sim171$dimer$rmsd.mean.free, lty=2, col="green")

axis(2)
axis(1, at=at, labels=c(0,10,20,30,40,50))
box()
legend("bottomright", legend=c("GroEL (ATP-free)","GroEL-ATP"), col=c("blue", "red"), lty=1, ncol=2)
mtext("(B) Dimers", side=3, line=0.7, at=0, cex=1.1, adj=0)



## plotting trimers
plot(sim116$trimer$rmsd[,1], type='l', col="grey80", ylim=c(1, 4.5),
     ylab="RMSD (Å)", xlab="Time (ns)", axes=FALSE )
  

for ( i in 2:6 ) {
  lines(sim116$trimer$rmsd[,i], lty=2, col="grey80")
  lines(sim129$trimer$rmsd[,i], lty=2, col="grey80")
}

lines(sim116$trimer$rmsd.mean, lty=1, col="blue")
lines(sim129$trimer$rmsd.mean.atp, lty=1, col="red")
#lines(sim129$trimer$rmsd.mean.free, lty=2, col="orange")
#lines(sim171$trimer$rmsd.mean.atp, lty=2, col="black")
#lines(sim171$trimer$rmsd.mean.free, lty=2, col="green")

axis(2)
axis(1, at=at, labels=c(0,10,20,30,40,50))
box()
legend("bottomright", legend=c("GroEL (ATP-free)","GroEL-ATP"), col=c("blue", "red"), lty=1, ncol=2)
mtext("(C) Trimers (1x cis + 2x trans)", side=3, line=0.7, at=0, cex=1.1, adj=0)



## plotting heptamers
plot(sim116$heptamer$rmsd[,1], type='l', col="blue", ylim=c(1, 4.5),
     ylab="RMSD (Å)", xlab="Time (ns)", axes=FALSE )
  
#lines(sim116$heptamer$rmsd[,2], lty=2, col="blue")
lines(sim129$heptamer$rmsd[,1], lty=1, col="red")
#lines(sim129$heptamer$rmsd[,2], lty=2, col="orange")
#lines(sim171$heptamer$rmsd[,1], lty=2, col="black")
#lines(sim171$heptamer$rmsd[,2], lty=2, col="green")

axis(2)
axis(1, at=at, labels=c(0,10,20,30,40,50))
box()

legend("bottomright", legend=c("GroEL (ATP-free)", "GroEL-ATP"), col=c("blue", "red"), lty=1, ncol=2)
mtext("(D) Heptameric complex", side=3, line=0.7, at=0, cex=1.1, adj=0)

## plotting complex
plot(sim116$complex$rmsd, type='l', col="blue", ylim=c(1, 4.5),
     ylab="RMSD (Å)", xlab="Time (ns)", axes=FALSE )
  
lines(sim129$complex$rmsd, lty=1, col="red")
#lines(sim171$complex$rmsd, lty=2, col="black")

axis(2)
axis(1, at=at, labels=c(0,10,20,30,40,50))
box()
legend("bottomright", legend=c("GroEL (ATP-free)","GroEL-ATP"), col=c("blue", "red"), lty=1, ncol=2)
mtext("(E) Double ring complex", side=3, line=0.7, at=0, cex=1.1, adj=0)
dev.off()