2013-03-04 38 views
0

可以说我有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) 
} 
+0

添加语言标签可能? – sabisabi 2013-03-04 09:15:33

+0

我得到:'.C(..)不可用于包tRap' – Arun 2013-03-04 09:17:29

+1

为什么你每次都会覆盖score?为什么不呢:'for(i in 1:length(matrix)){score [i] < - affinity(matrix [[i]],sequence)}' – Arun 2013-03-04 09:18:40

回答

2

尝试lapply

lapply(matrix, affinity, sequence) 

PS:这是非常糟糕的主意,打电话给你的数据matrix

+0

我得到了你的意思,只是为了清楚我有一个更大的矩阵不被称为矩阵:p – 2013-03-04 09:19:34

+0

我认为'sapply'应该更容易在这里作为'score',我想是一个数值。否则,你将不得不用一个'do.call(...)'来包装结果作为一个向量。 – Arun 2013-03-04 09:26:36

+1

@SanderVanderZeeuw :)。如果此答案解决了您的问题,请考虑upvoting并接受answeR :) – 2013-03-04 09:26:58