library(bio3d)
load("swissprot_blast.RData")

pdb <- read.pdb("1XCK_A.pdb")
sse <- dssp(pdb)
blast <- blast.pdb( seq.pdb(pdb), database = "swissprot" )
hits <- plot.blast(blast, cutoff=305)
seqs <- get.seq(hits$gi.id, outfile = "seqs_swissprot_cutoff305.fasta" )

#save(blast, seqs, file="swissprot_blast.RData")

ide <- identity(seqs$ali)
excl.rows <- which(ide[1,]<0.30)

seqs2 <- NULL
seqs2$ali <- seqs$ali[-excl.rows,]
seqs2$id <- seqs$id[-excl.rows]

ide <- identity(seqs2$ali)
k <- ide.filter(seqs2$ali, ide=ide, cutoff=0.9)
ide <- identity(seqs2$ali[k$ind,])

aln <- seqaln(seqs2$ali[k$ind,], id = seqs2$id[k$ind])

gaps <- gap.inspect(aln$ali)

ent10 <- conserv(aln, method="entropy10")
ent22 <- conserv(aln, method="entropy22")
ident <- conserv(aln, method="identity")
simil <- conserv(aln, method="similarity")

groel.inds <- which(aln$ali[1,]!="-")
col=(seq(1:length(groel.inds))*0)+1
#tmp <- c(134,475,34,480,58,209,92,78,284,367)
tmp <- c(83,327,177,322,412,475,179,390,176,376,7,11,23,75,177,400,17,101,479,488,85,499,412,489,231,310,65,524,31,91,84,506,77,507,92,54,78)
col[tmp]=3

pdf("consv.pdf", w=15, h=10)
par(mfcol=c(3,1))
plot.bio3d(ent10[groel.inds], sse=sse, ylab="Entropy", xlab="Residue Number", col=col)
#lines(ent22[groel.inds], lty=3, col="red")
mtext(3, at=-10, line=-2, text="A", cex=1.4)
mtext(3, at=sse$helix$start, text=LETTERS[1:18], adj=0)

plot.bio3d(c(freq/2, seq(1:24)*0), type="h", xlab="Residue Number",
     sse=sse, xaxt="n",ylab="Number of Contacts" )
mtext(3, at=-10, line=-2, text="B", cex=1.4)
mtext(3, at=sse$helix$start, text=LETTERS[1:18], adj=0)

plot.bio3d(c(rf, seq(1:22)*0), type="h", xlab="Residue Number",
     sse=sse, xaxt="n",ylab="RMSF (Å)" )
lines(as.numeric(pdb$atom[pdb$calpha,"b"])/10, lty=2, col="red")

mtext(3, at=-10, line=-2, text="C", cex=1.4)
mtext(3, at=sse$helix$start, text=LETTERS[1:18], adj=0)
dev.off()




#plot(c(freq, seq(1:24)*0), ident[groel.inds])

cont.avg <- NULL
mypdbs <- c("1XCK_A.pdb", "1SVT_A.pdb", "1SVT_N.pdb", "1KP8_A.pdb")
for ( i in 1:2 ) {

  pdb <- read.pdb(mypdbs[i])
  
  l <- dm(pdb,selection="//A/1:525////")
  #l <- dm.xyz(pdbs$xyz[i,1:1572], mask.lower=FALSE)
  l[upper.tri(l)]=NA  # make top diagonal NA

  ## Find residues with contacting atoms (<=5 Angstrom)
  inds.stru <- which(l<=5, arr.ind=TRUE)   

  ## Find non-consecutive residues (>5 residues sequence separation)
  seq.sep  <- abs(as.numeric(pdb$atom[inds.stru[,1],"resno"]) -
                  as.numeric(pdb$atom[inds.stru[,2],"resno"]))
  inds.seq <- which(seq.sep>5) # seperated by > 5 residues

  ## All-atom contacts (indices)
  inds <- inds.stru[inds.seq,]

  ## All-atom contacts (now in terms of residue numbers)
  tmp <- unique( paste(pdb$atom[inds.stru[,1],"resno"],
                       pdb$atom[inds.stru[,2],"resno"], sep="#") )
  
  contacts <- matrix(as.numeric(unlist(strsplit(tmp,split="#"))),
                     ncol=2, byrow=TRUE )

  cont.avg=rbind(cont.avg, contacts)
  
}

#freq  <- table(contacts)
freq  <- table(cont.avg)
xaxis <-as.vector(bounds(as.numeric(names(freq)))[,c(1,2)])
  
  
# Plot residue contacts

#x11()
pdf("contacts.pdf")
plot(freq, typ="h", xlab="Residue Number",
     xaxt="n",ylab="Number of Contacts" )
axis(1,at=xaxis,labels=xaxis)
dev.off()



load("../../../bio3d_groel/xrayPCA.RData")
rf <- rmsf(xyz)