1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
|
#
# computePI: To Compute the Theoretical Isoelectric Point
#
# seq is a protein sequence as a vector of single characters in upper case,
# note that there is no argument checking here.
#
computePI <- function(seq){
#
# Remove stop codons translated as the character '*' :
#
if(length(which(seq == "*")) != 0) seq <- seq[- which(seq == "*")]
compoAA <- table(factor(seq, levels = LETTERS))
nTermR <- which(LETTERS == seq[1])
cTermR <- which(LETTERS == seq[length(seq)])
computeCharge <- function(pH, compoAA, pK, nTermResidue, cTermResidue){
cter <- 10^(-pK[cTermResidue,1]) / (10^(-pK[cTermResidue,1]) + 10^(-pH))
nter <- 10^(-pH) / (10^(-pK[nTermResidue,2]) + 10^(-pH))
carg <- as.vector(compoAA['R'] * 10^(-pH) / (10^(-pK['R',3]) + 10^(-pH)))
chis <- as.vector(compoAA['H'] * 10^(-pH) / (10^(-pK['H',3]) + 10^(-pH)))
clys <- as.vector(compoAA['K'] * 10^(-pH) / (10^(-pK['K',3]) + 10^(-pH)))
casp <- as.vector(compoAA['D'] * 10^(-pK['D',3]) /(10^(-pK['D',3]) + 10^(-pH)))
cglu <- as.vector(compoAA['E'] * 10^(-pK['E',3]) / (10^(-pK['E',3]) + 10^(-pH)))
ccys <- as.vector(compoAA['C'] * 10^(-pK['C',3]) / (10^(-pK['C',3]) + 10^(-pH)))
ctyr <- as.vector(compoAA['Y'] * 10^(-pK['Y',3]) / (10^(-pK['Y',3]) + 10^(-pH)))
charge <- carg + clys + chis + nter - (casp + cglu + ctyr + ccys + cter)
return(charge)
}
critere <- function( p1, p2, p3, p4, p5){
computeCharge(pH = p1, compoAA = p2, pK = p3, nTermResidue = p4, cTermResidue = p5)^2
}
nlmres <- suppressWarnings(nlm(critere, 7, p2 = compoAA, p3 = SEQINR.UTIL$pk, p4 = nTermR, p5 = cTermR))
#
# If minimum is not zero, try whith a different guess:
#
while( ! identical(all.equal( nlmres$minimum, 0 ), TRUE))
{
nlmres <- suppressWarnings(nlm(critere, runif(1, 0, 14), p2 = compoAA, p3 = SEQINR.UTIL$pk, p4 = nTermR, p5 = cTermR))
}
return(nlmres$estimate)
}
|
