2016-05-05 3 views
-3

ピーク領域のいくつかの注釈を計算しているこの小さなコードがスクリプトにあります。以下のコードはスピードのボトルネックであり、CpGカウントを計算するためにこれを実行する必要がある約10万のリージョンがあるので、数時間かかる。それをスピードアップする方法はありますか?コードが長すぎます - スピードアップする方法

for (i in 1:nrow(dataMtx)){ 
peakCord<-gsub("chr", "", peakCord) 
peakSeq<-system(sprintf("samtools faidx genome.fa %s", peakCord[i]), intern=T) 
peakSeq<-gsub(">.*$", "", peakSeq) 
peakSeq<-paste(peakSeq, collapse='') 
dataMtx$CpGCount[i] <- sum(str_count(peakSeq, "CG")) 
print(i) 
} 
+0

あなたは 'apply()'関数とそれに関連するファミリーを試してみることができます。おそらくdata.tableまたはdplyrパッケージを使用したいと思うかもしれません。または、関数をコンパイルしてコンパイルしてください - 私は個人的にそれを試しませんでしたが。 – zacdav

答えて

0

このようなものが動作する可能性があります。どのくらい速くなるかわからない。

library(dplyr) 
library(stringi) 

result = 
    data_frame(peakCord = peakCord) %>% 
    rowwise %>% 
    mutate(peakCord.replace = 
      peakCord %>% 
      stri_replace_all_fixed("chr", ""), 
     peakSeq = 
      peakCord.replace %>% 
      sprintf("samtools faidx genome.fa %s", .) %>% 
      system(intern = T) %>% 
      stri_replace_all(">.*$", "") %>% 
      paste(collapse=''), 
     CpGCount = peakSeq %>% stri_count_fixed("CG")) 
+0

おかげさまでbramtayl ..私はさらに多くの地域(約500,000)を持っていると思うし、forループでRでやっているのは本当に遅いです。だから、私は同じためにperlスクリプトを書いてしまった。 Rループで撮影した4日間とは異なり、約30万リージョンで約1時間で仕事をすることができました。とにかく素晴らしい提案のおかげでありがとう。 – Sudhir

0

ここで私は上記のperkコードで述べたとおりです。誰かがそれを必要とする場合。コードでは、CpGの数だけでなくその位置もカウントされ、GCの割合も加算されます。

use strict; 
use warnings; 
BEGIN { our $start_run = time(); } 

open(POSITIONS,"mergedPeaks.bed"); # "mergedPeaks.bed"); 
my $filename='outfile.txt'; 
open(my $fh, '>', $filename) or die "Could not open file '$filename' $!"; 
my $string1="CG"; 
my @positions=(); 
while(<POSITIONS>){ 
chomp; 
unless($_=~ m/^chrom/){ 
my ($seqName,$begin,$end) = split(/\t/); 
$seqName=~s/chr//; 
open(SAMTOOLS,"samtools faidx /n/meissnerfs2/Everyone/sthakurela/annotationFiles/genomeFASTA/hs/hg19/genome.fa $seqName:$begin-$end |"); 
my @data = <SAMTOOLS>; 
chop(@data); 
my $seq=join("", @data); 
$seq =~ s/\d+|\:|\-//g; 
while ($seq =~ /$string1/gi){ 
push(@positions, pos($seq)- length($string1)); 
} 
my $length=scalar @positions; 
my $seqLen=length($seq); 
my $GC_count=($seq=~tr/GC/GC/); 
my $GCper=sprintf("%.2f", ($GC_count/$seqLen)*100); 
print $fh $_, "\t", $length,"\t", (join(",",@positions)), "\t", $GCper, "\n"; 
@positions=(); 
@data=(); 
$GC_count=0; 
$GCper=0; 
$seqLen=0; 
close(SAMTOOLS); 
}} 
close(POSITIONS); 
close $fh; 

my $end_run = time(); 
my $run_time = $end_run - our $start_run; 
print "Job took $run_time seconds\n"; 
関連する問題