# Analysis of GiardinB coiled-coil prediction
mat=read.delim("PMatrix.txt")
dim(mat)
per<- seq(0,56, by=7)
bgr<- mat[,5]
a<-apply ( mat[,6+per],1,sum)
b<-apply ( mat[,7+per],1,sum)
c<-apply ( mat[,8+per],1,sum)
d<-apply ( mat[,9+per],1,sum)
e<-apply ( mat[,10+per],1,sum)
f<-apply ( mat[,11+per],1,sum)
g<-apply ( mat[,12+per],1,sum)
poststate <- data.frame(a,b,c,d,e,f,g,bgr)
dim(poststate)

aa=145; bb=155
plot ( aa:bb, (((aa:bb)-aa)/(bb-aa)), type="n" )
for (i in 1:8)
{ lines (aa:bb, poststate[aa:bb,i], col=i)}

maPalette <-
function (low = "white", high = c("green", "red"), mid = NULL,
    k = 50)
{
    low <- col2rgb(low)/255
    high <- col2rgb(high)/255
    if (is.null(mid)) {
        r <- seq(low[1], high[1], len = k)
        g <- seq(low[2], high[2], len = k)
        b <- seq(low[3], high[3], len = k)
    }
    if (!is.null(mid)) {
        k2 <- round(k/2)
        mid <- col2rgb(mid)/255
        r <- c(seq(low[1], mid[1], len = k2), seq(mid[1], high[1],
            len = k2))
        g <- c(seq(low[2], mid[2], len = k2), seq(mid[2], high[2],
            len = k2))
        b <- c(seq(low[3], mid[3], len = k2), seq(mid[3], high[3],
            len = k2))
    }
    rgb(r, g, b)
}

aa<-140;bb<-162
values <- (as.matrix(poststate[aa:bb,1:cc] ))^0.6 
round(values,2);max<-0.6;values[values>max] <- max

colnb = 20
cols = maPalette(low = "white", high = "blue",  k=colnb) 
br<-unique(order(as.numeric(values)))
par(lab=c(11,1,3))
par(las=0)
image(x=aa:bb, y=1:cc, z=values , zlim=c(0,max(values)), 
	col = cols,  xaxs = "i", yaxs = "i", yaxt="n", xaxt="n",
	xlab="amino acid sequence", ylab="heptad position")
axis(2, at = seq(1, 7, by = 1), labels=c("a","b","c","d","e","f","g"))
axis(1, at = seq(aa, bb, by = 1), labels=as.character(mat[aa:bb,3],cex=0.7))
box(lty=7)

pdf(file ="FigGiardin140-162blue.pdf",
         width = 8.5, height = 4, onefile = TRUE, family = "Helvetica",
         title = "", fonts = "mono", paper = "a4")
par(lab=c(11,1,3))
par(las=0)
image(x=aa:bb, y=1:cc, z=values , zlim=c(0,max(values)), 
	col = cols,  xaxs = "i", yaxs = "i", yaxt="n",xaxt="n",
	xlab="amino acid sequence", ylab="heptad position")
axis(2, at = seq(1, 7, by = 1), labels=c("a","b","c","d","e","f","g"))
axis(1, at = seq(aa, bb, by = 1), labels=as.character(mat[aa:bb,3], cex.axis=0.7)  )
box(lty=7)
dev.off()   


colnb = 30;cols = maPalette(low = "white", high = "black",  k=colnb) 
par(lab=c(11,1,3))
par(las=0)
image(x=aa:bb, y=1:cc, z=values , zlim=c(0,max(values)), 
	col = cols,  xaxs = "i", yaxs = "i", yaxt="n", xaxt="n",
	xlab="amino acid sequence", ylab="heptad position")
axis(2, at = seq(1, 7, by = 1), labels=c("a","b","c","d","e","f","g"))
axis(1, at = seq(aa, bb, by = 1), labels=as.character(mat[aa:bb,3],cex=0.7))
box(lty=7)

pdf(file ="FigGiardin140-162bw.pdf",
         width = 8.5, height = 4, onefile = TRUE, family = "Helvetica",
         title = "", fonts = "mono", paper = "a4")
par(lab=c(11,1,3))
par(las=0)
image(x=aa:bb, y=1:cc, z=values , zlim=c(0,max(values)), 
	col = cols,  xaxs = "i", yaxs = "i", yaxt="n",xaxt="n",
	xlab="amino acid sequence", ylab="heptad position")
axis(2, at = seq(1, 7, by = 1), labels=c("a","b","c","d","e","f","g"))
axis(1, at = seq(aa, bb, by = 1), labels=as.character(mat[aa:bb,3], cex.axis=0.7)  )
box(lty=7)
dev.off()