-
Notifications
You must be signed in to change notification settings - Fork 1
/
PAAC_BAC.R
100 lines (72 loc) · 2.52 KB
/
PAAC_BAC.R
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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#' Pseudo Amino Acid Composition (PseAAC) Descriptor
#'
#' This function calculates the Pseudo Amino Acid Composition (PseAAC)
extractPAAC_BAC <- function(
x, y, props = c("Hydrophobicity", "Hydrophilicity", "SideChainMass"),
lambda = 10, w = 0.05, customprops = NULL) {
protcheck <- function(x) {
AADict <- c(
"A", "R", "N", "D", "C", "E", "Q", "G", "H", "I",
"L", "K", "M", "F", "P", "S", "T", "W", "Y", "V"
)
all(strsplit(x, split = "")[[1]] %in% AADict)
}
if (protcheck(x) == FALSE) {
stop("x has unrecognized amino acid type")
}
if (nchar(x) <= lambda) {
stop('Length of the protein sequence must be greater than "lambda"')
}
AAidx <- read.csv(y, header = TRUE)
tmp <- data.frame(
AccNo = c("Hydrophobicity", "Hydrophilicity", "SideChainMass"),
A = c(0.62, -0.5, 15), R = c(-2.53, 3, 101),
N = c(-0.78, 0.2, 58), D = c(-0.9, 3, 59),
C = c(0.29, -1, 47), E = c(-0.74, 3, 73),
Q = c(-0.85, 0.2, 72), G = c(0.48, 0, 1),
H = c(-0.4, -0.5, 82), I = c(1.38, -1.8, 57),
L = c(1.06, -1.8, 57), K = c(-1.5, 3, 73),
M = c(0.64, -1.3, 75), F = c(1.19, -2.5, 91),
P = c(0.12, 0, 42), S = c(-0.18, 0.3, 31),
T = c(-0.05, -0.4, 45), W = c(0.81, -3.4, 130),
Y = c(0.26, -2.3, 107), V = c(1.08, -1.5, 43)
)
AAidx <- rbind(AAidx, tmp)
if (!is.null(customprops)) AAidx <- rbind(AAidx, customprops)
aaidx <- AAidx[, -1]
row.names(aaidx) <- AAidx[, 1]
n <- length(props)
# Standardize H0 to H
H0 <- as.matrix(aaidx[props, ])
H <- matrix(ncol = 20, nrow = n)
for (i in 1:n) {
H[i, ] <-
(H0[i, ] - mean(H0[i, ])) / (sqrt(sum((H0[i, ] - mean(H0[i, ]))^2) / 20))
}
AADict <- c(
"A", "R", "N", "D", "C", "E", "Q", "G", "H", "I",
"L", "K", "M", "F", "P", "S", "T", "W", "Y", "V"
)
dimnames(H) <- list(props, AADict)
# Compute (big) Theta
Theta <- vector("list", lambda)
xSplitted <- strsplit(x, split = "")[[1]]
N <- length(xSplitted)
for (i in 1:lambda) {
for (j in 1:(N - i)) {
Theta[[i]][j] <- mean((H[, xSplitted[j]] - H[, xSplitted[j + i]])^2)
}
}
# Compute (small) theta
theta <- sapply(Theta, mean)
# Compute first 20 features
fc <- summary(factor(xSplitted, levels = AADict), maxsum = 21)
Xc1 <- fc / (1 + (w * sum(theta)))
names(Xc1) <- paste("Xc1.", names(Xc1), sep = "")
# Compute last lambda features
Xc2 <- (w * theta) / (1 + (w * sum(theta)))
names(Xc2) <- paste("Xc2.lambda.", 1:lambda, sep = "")
# Combine (20 + lambda) features
Xc <- c(Xc1, Xc2)
return (Xc)
}