Code zu lang - Wie man beschleunigt

Question

Ich habe dieses kleine Stück Code in meinem Skript, das mehrere Annotation für Peak-Regionen berechnet. Der unten stehende Code ist ein Engpass für die Geschwindigkeit und dauert einige Stunden, da ich ungefähr 100.000 Regionen habe, für die ich dies ausführen muss, um die Anzahl der CpGs zu berechnen. Gibt es eine Möglichkeit, es zu beschleunigen?

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) 
} 

Answer 1

So etwas könnte funktionieren. Nicht sicher, wie viel schneller es sein wird.

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")) 

Answer 2

hier ist, was ich oben wie perc code erwähnt. Falls jemand es braucht. Der Code zählt nicht nur die Anzahl der CpGs, sondern auch deren Positionen und fügt GC Prozent hinzu.

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";