可以說我有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)
}
添加語言標籤可能? – sabisabi 2013-03-04 09:15:33
我得到:'.C(..)不可用於包tRap' – Arun 2013-03-04 09:17:29
爲什麼你每次都會覆蓋score?爲什麼不呢:'for(i in 1:length(matrix)){score [i] < - affinity(matrix [[i]],sequence)}' – Arun 2013-03-04 09:18:40