2013-03-04 47 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