使用 d_ply 进行 For/while 循环优化答案

【问题标题】：For/while loop optimization with d_ply使用 d_ply 进行 For/while 循环优化
【发布时间】：2013-08-08 11:46:09
【问题描述】：

我正在运行一个自定义函数，它为 data.frame 的每一行生成并保存一些图：

zz <- "chr  start end   name  max
chr11 66184332 66197785 NPAS4_CBP20 90
chr11  62666002 62683613 BC047540_CBP20 100   
chr1 9824542  9828548  MIR3687_CBP20  500
chr6  33239767 33259282 B3GALT4_Pol2   1000
chr20  244996112   245029580   HNRNPU-AS1_Pol2   450
chr20 62487823 62525914 ABHD16B_Pol2   370
chr12 121146198   121179996   ACADS_Pol2  90"

my.genes <- read.table(text=zz, header = TRUE)

该函数有两个步骤，其中一个很慢（约 40 秒），但只需要为每个 chr 运行一个，然后为该 chr 的每一行运行一个。在伪代码中它会是这样的：

for each chr createData
   for each nameInChr do somethingWithData.

我的问题是，如何优化这个？嵌套 d_ply？按 chr 对 data.frame 进行排序，然后对 unique(mygenes$chr) 进行某种形式的应用？

如果这不包含实际功能，我很抱歉，但它很长，我认为这更像是一个很好的实践/理论问题。但是，如果需要，我可以添加适当的代码。

编辑

这是简化的函数 (plotGviz)。事实上，我认为 d_ply 不会起作用，因为它会为每一行运行 UcscTrack（我正确吗？）。最终目标是为每个独特的染色体运行 UcscTrack，然后使用该信息为每个基因（名称）构建图。这两个步骤（UcscTrack 和绘图）都很耗时，但基本原理是，首先运行代码的非唯一部分（适用于 chr）将节省整个时间 - 该表只是我所拥有的一小部分。

自定义 plotGviz 中的功能不能仅按照我设置的方式进行优化。

################
### libraries ##
library(Gviz)
library(GenomicRanges)
library(GenomicFeatures)
library(data.table)
library("RColorBrewer")
#######################


d_ply(my.genes, .(chr, name), plotGviz)

plotGviz <- function(gene) {
   chr <- gene$chr
   mygene.start <- gene$start
   mygene.end <- gene$end
   mygene <- gene$name
   max.cov <- gene$max

#################################################
## this part runs only once for each unique chr
## UcscTrack is what takes most time

   ## Gene annotations ##
   knownGenes <- UcscTrack(genome=gen, chromosome=chr,
   track="knownGene",  trackType="GeneRegionTrack",
   rstarts="exonStarts", rends="exonEnds", gene="name", symbol="name",
   transcript="name", strand="strand", fill="#FF7F00", name="UCSC Genes", geneSymbols = TRUE, showId = TRUE)

   refGenes <- UcscTrack(genome=gen, chromosome=chr,
   track="refGene", trackType="GeneRegionTrack", 
   rstarts="exonStarts", rends="exonEnds", gene="name", symbol="name",
   transcript="name", strand="strand", fill="#FF7F00", name="RefSeq Genes", geneSymbols = TRUE, showId = TRUE)

   ## axis scale ##
   ideoTrack <- IdeogramTrack(genome = gen, chromosome = chr)
   # plotTracks(ideoTrack, from = mygene.start, to = mygene.end)

   axisTrack <- GenomeAxisTrack(scale=0.25)


####################################
## now for each name in unique chr
######################

   ## construct tracks ##


   ## ChIP-seq coverage from BAM ##
   bamPol2 <- "~/bioinfo/srp_chip_seq/data/reads/merged_reads_CBP20line/Pol2.bam"
   bamInput <- "~/bioinfo/srp_chip_seq/data/reads/merged_reads_CBP20line/Input.bam"

  ## histogram
   bTrackPol2 <- DataTrack(range = bamPol2, genome = gen, chromosome = chr,
      name = "Pol2", type = "histogram", col.histogram="#984EA3",
      ylim=c(0,max.cov))

   bTrackInput <- DataTrack(range = bamInput, genome = gen, chromosome = chr,
      name = "Input", type = "histogram", col.histogram="#999999",
      ylim=c(0,max.cov))


   #################
   ## plot tracks ##

   pdf(paste("./plots/", mygene,"_cov.pdf", sep="")) 

   plotTracks(list(ideoTrack, axisTrack, bTrackCBP20, pTrackCBP20, bTrackPol2, pTrackPol2, bTrackInput, refGenes, bTrackRNAcov), from = mygene.start, to = mygene.end, fontfamily="Helvetica", background.title="white", col.title="black", col.axis="black")

   plotTracks(list(ideoTrack, axisTrack, bTrackCBP20, pTrackCBP20, bTrackPol2, pTrackPol2, bTrackInput, knownGenes, bTrackRNAcov), from = mygene.start, to = mygene.end, fontfamily="Helvetica", background.title="white", col.title="black", col.axis="black")
   dev.off()

   # plotTracks(list(axisTrack, refGenes, bTrackCBP20, pTrackCBP20, bTrackPol2, pTrackPol2, bTrackInput, bTrackRNA), from = mygene.start, to = mygene.end, fontfamily="Helvetica", background.title="white", col.title="black", col.axis="black")

   # head(displayPars(bTrackPol2))
}

edit2

我的临时解决方案为唯一的 chr 使用 for 循环，在为每个 chr 子设置所有名称之前创建 UcscTrack，并为绘图创建 d_pply。基本上在 2 个部分中吐出这么大的功能。

for (chrom in unique(my.genes$chr)) {
   print(paste("creating annotation for ", chrom, sep=""))

   ## Gene annotations ##
   knownGenes <- UcscTrack(genome=gen, chromosome=chrom,
   track="knownGene",  trackType="GeneRegionTrack",
   rstarts="exonStarts", rends="exonEnds", gene="name", symbol="name",
   transcript="name", strand="strand", fill="#FF7F00", name="UCSC Genes", geneSymbols = TRUE, showId = TRUE)

   refGenes <- UcscTrack(genome=gen, chromosome=chrom,
   track="refGene", trackType="GeneRegionTrack", 
   rstarts="exonStarts", rends="exonEnds", gene="name", symbol="name",
   transcript="name", strand="strand", fill="#FF7F00", name="RefSeq Genes", geneSymbols = TRUE, showId = TRUE)

   ## axis scale ##
   ideoTrack <- IdeogramTrack(genome = gen, chromosome = chrom)
   # plotTracks(ideoTrack, from = mygene.start, to = mygene.end)

   axisTrack <- GenomeAxisTrack(scale=0.25)

   ## select genes in chromosome
   df.g <- subset(my.genes, chr == chrom)
   d_ply(df.g, .(name), plotGviz)
}

【问题讨论】：

知道createData 和somethingWithData 是什么，这很难回答。我认为当问题是优化时，一般的答案是困难的，因为这将取决于正在发生的事情的具体情况
谢谢。我将编辑我的答案。
那么慢到无法接受吗？您可以尝试并行化。
第 1 步，UcscTrack 的 2 次调用，总共需要大约 40-50 秒，而绘图又增加了 20 秒（大约）。只要我不制作多于一把地块，这是可以接受的，但我正计划扩大它。关于并行化有什么建议吗？
plyr 让它变得非常简单。使用foreach 框架的后端之一，然后将d_ply 替换为d_ply(..., .parallel=TRUE, .paropts=list(.packages=c("Gviz", etc.)))。此 PDF 解释了第一步（以及与使用 foreach 和 plyr 不直接相关的更多内容）：cran.r-project.org/web/packages/doParallel/vignettes/…

标签： r loops split plyr

【解决方案1】：

试试data.table，它会比plyr 快。这是解决问题的通用方法：

library(data.table)
dt = data.table(my.genes)

dt[, .SD[, do_smth(), by = name], by = chr]

这是在另一个 SO 问题中使用上述方法的示例 - how to integrate properties defined on multiple rows using a data.frame or data.table long format approach

【讨论】：