假设我有 4 个矩阵：
matrix<-list(
MA0275.1 = structure(c(0, 76, 0, 24, 0, 100, 0, 0, 0, 0,
100, 0, 0, 0, 100, 0, 72, 11, 16, 0, 53, 0, 0, 47), .Dim = c(4L,
6L), .Dimnames = list(c("A", "C", "G", "T"), NULL), id = "MA0275.1", accession = "ASG1"),
MA0276.1 = structure(c(0, 220, 8, 35, 0, 291, 0, 3, 61, 21,
133, 10, 58, 54, 101, 12, 130, 0, 54, 0, 0, 11, 8, 147, 33,
150, 8, 35, 80, 0, 92, 26, 0, 8, 249, 19, 0, 0, 256, 18), .Dim = c(4L,
10L), .Dimnames = list(c("A", "C", "G", "T"), NULL), id = "MA0276.1", accession = "ASH1"),
MA0277.1 = structure(c(63, 13, 13, 13, 100, 0, 0, 0, 100,
0, 0, 0, 88, 13, 0, 0, 75, 0, 25, 0, 0, 0, 100, 0, 78, 16,
3, 3, 81, 6, 6, 6, 63, 13, 13, 13), .Dim = c(4L, 9L), .Dimnames = list(
c("A", "C", "G", "T"), NULL), id = "MA0277.1", accession = "AZF1"),
MA0278.1 = structure(c(64, 217, 425, 292, 104, 552, 150,
192, 484, 111, 114, 288, 78, 401, 186, 333, 455, 51, 370,
122, 248, 34, 670, 46, 98, 724, 143, 33, 52, 918, 7, 20,
348, 346, 280, 24, 12, 3, 977, 6, 966, 5, 23, 4, 26, 6, 962,
4, 9, 10, 4, 975, 47, 930, 7, 15, 892, 42, 16, 49, 487, 123,
320, 68, 288, 140, 317, 254, 373, 110, 81, 434, 178, 367,
184, 268, 402, 140, 341, 114, 435, 229, 241, 94), .Dim = c(4L,
21L), .Dimnames = list(c("A", "C", "G", "T"), NULL), id = "MA0278.1", accession = "BAS1"))
这些矩阵用于对序列进行评分（亲和力）。我拥有的功能使我能够给我 1 个亲和力分数，因此只使用第一个矩阵。
但如果分数是在 data.frame() 中为每个矩阵定义的，那就太好了。我试图用一个 for 循环来做到这一点。
sequence<-"GCCTTTCCTTCTCTTCTCCGCGTGTGGAGGGAGCCAGCGCTTAGGCCGGAGCGAGCCTGGGGGCCGCCCGCCGTGAAGACATCGCGGGGACCGATTCACC"
for (i in matrix) {
score<-affinity(i,sequence)
}
这给了我一个 1 矩阵的数值。所以 for 循环不能正常工作。我希望它给我每个矩阵的所有亲和力分数。
亲和力函数：
affinity<-function (pwm, seq, Rmax = NULL, lambda = 0.7, pseudo.count = 1,
gc.content = 0.5, slide = FALSE)
{
if (is.null(seq) || is.na(seq) || mode(seq) != "character") {
stop("sequence must be a character string of length >= ncol(pwm)")
}
gap.pos = sapply(1:nchar(seq), function(i) substr(seq, i,
i) == "-")
seq = gsub("-", "", seq)
if (nchar(seq) < ncol(pwm)) {
stop("sequence must be a character string of length >= ncol(pwm)")
}
Rmax = ifelse(is.null(Rmax), exp(0.584 * ncol(pwm) - 5.66),
Rmax)
pwm = pwm + pseudo.count
at.content = 1 - gc.content
pwm = apply(pwm, 2, function(p) {
maxAT = max(p[c(1, 4)])
maxCG = max(p[c(2, 3)])
if (maxAT > maxCG) {
transformed = c(log(maxAT/p[1])/lambda, log((maxAT/at.content) *
(gc.content/p[2]))/lambda, log((maxAT/at.content) *
(gc.content/p[3]))/lambda, log(maxAT/p[4])/lambda)
}
else {
transformed = c(log((maxCG/gc.content) * (at.content/p[1]))/lambda,
log(maxCG/p[2])/lambda, log(maxCG/p[3])/lambda,
log((maxCG/gc.content) * (at.content/p[4]))/lambda)
}
if (maxAT == maxCG) {
transformed = log(maxAT/p)/lambda
}
return(transformed)
})
if (slide) {
z = .C("R_affinity", as.double(pwm), ncol(pwm), as.character(seq),
as.integer(nchar(seq)), as.double(Rmax), as.double(lambda),
double(length = nchar(seq) - ncol(pwm) + 1), PACKAGE = "tRap")
res = z[[7]]
gapped = numeric(length = nchar(seq) - ncol(pwm) + 1)
gapped[!gap.pos] = res
res = gapped
}
else {
z = .C("R_affinity_sum", as.double(pwm), ncol(pwm), as.character(seq),
as.integer(nchar(seq)), as.double(Rmax), as.double(lambda),
double(length = 1), PACKAGE = "tRap")
res = z[[7]]
}
return(res)
}