0

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)

その関数には 2 つのステップがあり、そのうちの 1 つは低速 (~40 秒) ですが、chr ごとに 1 つ、次にその chr の行ごとに実行するだけで済みます。擬似コードでは、次のようになります。

for each chr createData
   for each nameInChr do somethingWithData.

私の質問は、これをどのように最適化するのですか? ネストされた d_ply? data.frame を chr で並べ替えてから、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))
}

編集2

一意の 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)
}
4

1 に答える 1

0

を試してみてください。data.tableより速くなりplyrます。問題に対する一般的なアプローチは次のとおりです。

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

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

そして、別のSOの質問での上記のアプローチの例を次に示します-data.frameまたはdata.tableロングフォーマットアプローチを使用して、複数の行で定義されたプロパティを統合する方法

于 2013-08-08T16:40:21.593 に答える