Rscript/get_value_beta.R

In Solipugids/microarray1433: Affinity-Based Prediction of Phosphopeptide-Binding Specificity

setwd("/Users/wangchao/Development/microarray1433/Rscript")
orig_data <- read.table("../data/array.dat",header=T,colClasses=c("integer","character","character","numeric","integer","factor"))
library(Biostrings)
data(BLOSUM62)
library(reshape2)

#normorlised BLOSUM62
matrix_transform <- function(x){
	output <- input <- x[-(nrow(x)-4):-nrow(x),-(ncol(x)-4):-ncol(x)]
	for (i in 1:nrow(input)) {
		for (j in 1:ncol(input)){
			output[i,j] = (input[i,j]/sqrt(input[i,i]*input[j,j])+1)/2
		}
	}
	return(output)
}

nBLOSUM <- matrix_transform(BLOSUM62)

aff.dist <- function(x){
	dist <- log2(max(x[1],x[2])/min(x[1],x[2]))
	return(dist)
}

seq.dist <- function(x) {
	seq <- gsub(pattern = "XXX|pS/T",replacement = "", x)
	seq.mat <- strsplit(seq, "")
	a <- nBLOSUM[seq.mat[[1]][1],seq.mat[[2]][1]]
	b <- nBLOSUM[seq.mat[[1]][2],seq.mat[[2]][2]]
	c <- nBLOSUM[seq.mat[[1]][3],seq.mat[[2]][3]]
	single <- matrix(c( 1-a, 1-b, 1-c, sqrt((1-a) * (1-b)), sqrt((1-a) * (1-c)), sqrt((1-b) * (1-c)), ((1-a) * (1-b) * (1-c))^(1/3)),byrow=T,nrow=1)
	return(single)
}

cal.lm <- function(row, m) {
	seq.list <- m[(row-5):(row+5), 3]
	affinity.list <- log2(m[(row-5) : (row+5), 4]+1)
	seq.mat <- matrix(c(seq.list,rep(seq.list[6], length(seq.list))), ncol = 2)
	seq.input <- t(apply(seq.mat, 1, seq.dist))
	affinity.mat <- matrix(c(affinity.list, rep(affinity.list[6], length(affinity.list))), ncol = 2)
	affinity.mat <- affinity.mat + 1 
	affinity.input <- t(apply(affinity.mat, 1, aff.dist))
	seq.input <- seq.input[-6,]
	affinity.input <- affinity.input[1,-6]
	output <- lm(affinity.input ~ seq.input, x = T, y = T)
	return(output$coefficient)
}

mainfun <- function(isoforms) {
	two.sides <- orig_data[orig_data$Isoforms == isoforms,]
	left.side <- two.sides[grep("pS/TXXX",two.sides$Sequence),]
	right.side <- two.sides[grep("XXXpS/T",two.sides$Sequence),]

	left.side <- left.side[order(left.side$RANK),]
	right.side <- right.side[order(right.side$RANK),]

	left.nrow <- 6 : (nrow(left.side) - 5)
	left.exp <- lapply(left.nrow, cal.lm, left.side)
	left.exp <-do.call(rbind, left.exp)

	right.nrow <- 6 : (nrow(right.side) -5)
	right.exp <- lapply(right.nrow, cal.lm, right.side)
	right.exp <- do.call(rbind, right.exp)
	biaoti <- paste("LEFT", isoforms,sep=" ")
	boxplot(left.exp, names = F, axes = F, main = biaoti, outline = F)
	axis(2)
	biaoti <- paste("RIGHT", isoforms,sep=" ")
	boxplot(right.exp, names = F, axes = F, main = biaoti, outline = F)
	axis(2)
	#image(1:ncol(right.exp), 1:nrow(right.exp), t(right.exp), axes = F, xlab = '', ylab = biaoti)
}

layout(matrix(1:14, byrow = F, nrow = 2))
lapply(levels(orig_data$Isoforms), mainfun)