对数据框列表应用 sapply 或其他应用函数，而不是嵌套 for 循环

Question

我有两个数据框/数据列表“humanSplitandratSplit”，它们的格式为

> ratSplit$Kidney_F_GSM1328570
  ratGene        ratReplicate alignment RNAtype
1    Crot Kidney_F_GSM1328570         7     REV
2    Crot Kidney_F_GSM1328570        12     REV
3    Crot Kidney_F_GSM1328570         4     REV

和

> humanSplit$Fetal_Brain_408_AGTCAA_L009_R1_report.txt
   humanGene                            humanReplicate alignment RNAtype
53     ZFP28 Fetal_Brain_408_AGTCAA_L009_R1_report.txt         5     reg
55     RC3H1 Fetal_Brain_408_AGTCAA_L009_R1_report.txt         9     reg
56     IFI27 Fetal_Brain_408_AGTCAA_L009_R1_report.txt         4     reg

下面使用的另一个文件是geneList的形式：

ABAT,Abat
ABCA1,Abca1
ABCA12,Abca12
ABCA2,Abca2
ABCA3,Abca17
ABCA4,Abca4
ABCA5,Abca5

现在我想在经过一些数据处理后在ratSplit 和humanSplit 之间的所有元素对组合之间进行费希尔精确测试。并最终想将 Fisher 测试的结果写入csv 文件中。现在我正在做双重循环。但我想知道如何使用sapply 或其他相关的东西来提高效率。

目前我正在做以下事情：在这里我首先制作一个data.frameresult，在其中我存储/附加从费舍尔测试中获得的所有信息在每个步骤中成对。最后，当整个循环完成时，我将最终的result 写入csv 文件。我的理解是使用sapply我需要将循环内部转换成一个函数，然后调用sapply。但我不确定优化它的最佳方法是什么。任何帮助将不胜感激

result <- data.frame(humanReplicate = "human_replicate", ratReplicate = "rat_replicate", pvalue = "p-value", alternative = "alternative_hypothesis", 
                     Conf.int1 = "conf.int1", Conf.int2 ="conf.int2", oddratio = "Odd_Ratio")
for(i in 1:length(ratSplit)) {
  for(j in 1:length(humanSplit)) {
    ratReplicateName <- names(ratSplit[i])
    humanReplicateName <- names(humanSplit[j])

    #merging above two based on the one-to-one gene mapping as in geneList defined above.
    mergedHumanData <-merge(geneList,humanSplit[[j]], by.x = "human", by.y = "humanGene")
    mergedRatData <- merge(geneList, ratSplit[[i]], by.x = "rat", by.y = "ratGene")

    mergedHumanData <- mergedHumanData[,c(1,2,4,5)] #rearrange column
    mergedRatData <- mergedRatData[,c(2,1,4,5)]  #rearrange column
    mergedHumanRatData <- rbind(mergedHumanData,mergedRatData) #now the columns are "human", "rat", "alignment", "RNAtype"

    agg <- aggregate(RNAtype ~ human+rat, data= mergedHumanRatData, FUN=getGeneType) #agg to make HmYn form
    HmRnTable <- table(agg$RNAtype) #table of HmRn ie RNAtype in human and rat.

    #now assign these numbers to variables HmYn. Consider cases when some form of HmRy is not present in the table. That's why
    #is.integer0 function is used
    HyRy <- ifelse(is.integer0(HmRnTable[names(HmRnTable) == "HyRy"]), 0, HmRnTable[names(HmRnTable) == "HyRy"][[1]])
    HnRn <- ifelse(is.integer0(HmRnTable[names(HmRnTable) == "HnRn"]), 0, HmRnTable[names(HmRnTable) == "HnRn"][[1]])
    HyRn <- ifelse(is.integer0(HmRnTable[names(HmRnTable) == "HyRn"]), 0, HmRnTable[names(HmRnTable) == "HyRn"][[1]])
    HnRy <- ifelse(is.integer0(HmRnTable[names(HmRnTable) == "HnRy"]), 0, HmRnTable[names(HmRnTable) == "HnRy"][[1]])

    contingencyTable <- matrix(c(HnRn,HnRy,HyRn,HyRy), nrow = 2)

    fisherTest <- fisher.test(contingencyTable) 
    newLine <- data.frame(t(c(humanReplicate = humanReplicateName, ratReplicate = ratReplicateName, pvalue = fisherTest$p,
                              alternative = fisherTest$alternative, Conf.int1 = fisherTest$conf.int[1], Conf.int2 =fisherTest$conf.int[2], 
                              oddratio = fisherTest$estimate[[1]])))


    result <-rbind(result,newLine)
  }
}

write.table(result, file = "newData5.csv", row.names = FALSE, append = FALSE, col.names = TRUE, sep = ",")

Answer 1

由于我们缺少geneList，因此很难对此进行测试，但我推断代码可以正常工作，因此您只想加快速度。以下是一些帮助提示：

1. 不要像这样预定义result。通过将每列的第一个条目设置为列名的字符串，您可以确保所有后续条目都将被强制转换为字符串。虽然这并不一定会咬你，因为无论如何你最终都会写入 CSV，但这是一种糟糕的形式。如果您打算将其读回 R，它会咬您一口，因为第一行将用作列名，但第二行将全部是字符串，迫使后续数据也为字符串。 （然后你必须清理自己的数据，浪费。）

2. 在脚本结束时，您调用rbind。这可能会在一段时间内没问题，但重复调用rbind 将导致每次都复制整个data.frame。随着更多的行被追加，这将导致您的代码显着变慢。这可以通过下面列出的两种方法之一来解决。

3. 由于您使用每个names(HmRnTable) == "HyRy" 两次，我的技术是首先将其保存到一个向量（或标量，如果您使用which(...)），然后在HmRnTable 的子集中使用此变量。它可能会加快速度一点，但也可能使代码更易于阅读。事实上，您可以将这四个作业中的每一个都缩短为（未经测试）：

   idx <- is.integer0(HmRnTable[names(HmRnTable) == 'HyRy'])
   HyRy <- HmRnTable[idx][[1]]
   HyRy[idx] <- 0
   ## repeat for HyRn, HnRy, and HnRn
   

4. 我强烈建议您将大部分代码放入一个函数中，该函数接受两个参数（i 和 j）或四个参数（list1、index1、list2、index2）。 （你做什么取决于你有多少关于变量范围引用的强迫症。）我假设这个函数会从fisher.test返回结果，而不需要按摩。这将使创建函数时的测试更容易，并在稍后包含在此脚本中。我将在下面将函数引用为myfisher(i,j)。

5. 我推断您要运行大量比较（因为每次迭代确实不应该花费那么长时间）。 @Floo0 关于outer 的评论可以工作，expand.grid 也可以。无论哪种方式，您都在将 list1 的每个元素与 list2 的每个元素进行比较。如果你开始：

   (eg <- expand.grid(i = 1:length(ratSplit),
                      j = 1:length(humanSplit)))
   ##    i j
   ## 1  1 1
   ## 2  2 1
   ## 3  3 1
   ## 4  4 1
   ## ...
   

   这给了我们一个简单的data.frame，我们可以在上面使用apply。不过，老实说，我喜欢ddply 在这种情况下的优雅，因为它很容易从 data.frame 到 data.frame。

   library(plyr)
   ret <- ddply(eg, .(i, j), function(df) {
       with(myfisher(df$i, df$j),
            data.frame(pv = p.value, ci1 = conf.int[1], ci2 = conf.int[2],
                       alt = alternative))
   })
   

   请注意，我明确没有包括 humanReplicateName 或 ratReplicateName，因为它们可以添加到 ddply 之前（或之后，引用 ret$ 而不是 eg$）：

   eg$ratReplicateName <- names(ratSplit[ eg$i ])
   eg$humanReplicateName <- names(humanSplit[ eg$i ])
   

   名称也会神奇地出现在输出中。在循环内部处理的更少。

   到目前为止，这将产生一个 data.frame，然后您可以将其保存到 CSV。

   我将再提出一项建议，根据运行时间长短，这可能有点矫枉过正。我有时会发现我的长时间运行的脚本被打断了，也许是因为我不得不即时调整一些东西；我发现了一个错误；或者如果计算机必须重新启动。没有简单的方法可以让中游继续，但我已经采用了一些技术来缓解这种情况。

6. 使用不假定返回值类型的d_ply 而不是ddply。不要返回（单行）data.frame，而是立即将该行（有或没有标题，您的调用）保存到文件中。尽管以下代码并不是真正的“原子”代码，因此可能会出现竞争条件，但它足够强大，可以满足我们的大部分需求：

   library(plyr)
   savedir <- './output'
   ret <- ddply(eg, .(i, j), function(df) {
       fn <- file.path(savedir, sprintf('%s-%s.csv', df$i, df$j))
       if (! file.exists(fn)) {
           ret <- with(myfisher(df$i, df$j),
                       data.frame(pv = p.value, ci1 = conf.int[1], ci2 = conf.int[2],
                                  alt = alternative))
           write.table(ret, file = fn, sep = ",", append = FALSE, 
                       row.names = FALSE, col.names = TRUE)
       }
   })
   

   这样做的一个好处是，如果/当它被中断时，最您要做的就是查看“./output/”中的所有文件，删除 0 字节的文件，然后重新运行，它只会在丢失的文件上执行。哦，它变得更好了。

7. 并行化它。如果您需要做到这一点（并且某些功能的改进不如其他功能），您可以在系统上使用多个内核。你可以这样做：

   library(parallel)
   cl <- makeCluster(detectCores() - 1) # I like to keep one free
   clusterEvalQ(cl, {
       load('myTwoLists.rda') # with ratSplit and humanSplit variables
       source('mycode.R') # with myfisher(i,j)
       ## any libraries you may want/need to add, if you get more advanced
   })
   eg <- expand.grid(i = 1:length(ratSplit),
                     j = 1:length(humanSplit))
   eg$ratReplicateName <- names(ratSplit[ eg$i ])
   eg$humanReplicateName <- names(humanSplit[ eg$i ])
   ign <- parApply(cl, eg, 1, function(r) {
       i <- r[1] ; j <- r[2]
       fn <- file.path(savedir, sprintf('%s-%s.csv', i, j)
       if (! file.exists(fn)) {
           ret <- with(myfisher(i, j),
                       data.frame(ratName = r[3], humanName = r[4],
                                  pv = p.value, ci1 = conf.int[1], ci2 = conf.int[2],
                                  alt = alternative))
           write.table(ret, file = fn, sep = ",", append = FALSE, 
                       row.names = FALSE, col.names = TRUE)
       }
   })
   stopCluster(cl)
   

   请注意，我们不再将行引用为 data.frame。如果 Hadley 的代码可以在本地并行工作，我会很高兴的。 （据我所知，它确实存在，但我还没有找到它！） 编辑：我过去使用parallel 的项目都没有使用ddply，反之亦然，所以我从来没有玩过使用foreach的ddply(..., .parallel = TRUE)。

Caveat Emptor：此代码尚未在此上下文中进行测试，尽管它几乎是从工作代码中复制/粘贴的。由于编辑的原因，我可能会不合时宜或缺少括号。希望对您有所帮助！