2014-12-04 21 views
7

Jestem całkiem nowy w Pythonie i napisałem (prawdopodobnie bardzo brzydki) skrypt, który ma losowo wybierać podzbiór sekwencji z pliku fastq. Plik fastq przechowuje informacje w blokach po cztery wiersze. Pierwszy wiersz w każdym bloku zaczyna się od znaku "@". Plik fastq, którego używam jako pliku wejściowego, ma rozmiar 36 GB i zawiera około 14 000 000 linii.Jak mogę przyspieszyć mój skrypt Pythona?

Próbowałem przerobić istniejący skrypt, który wykorzystał zbyt dużo pamięci i udało mi się bardzo zmniejszyć zużycie pamięci. Ale skrypt trwa wiecznie i nie rozumiem dlaczego.

parser = argparse.ArgumentParser() 
parser.add_argument("infile", type = str, help = "The name of the fastq input file.", default = sys.stdin) 
parser.add_argument("outputfile", type = str, help = "Name of the output file.") 
parser.add_argument("-n", help="Number of sequences to sample", default=1) 
args = parser.parse_args() 


def sample(): 
    linesamples = [] 
    infile = open(args.infile, 'r') 
    outputfile = open(args.outputfile, 'w') 
    # count the number of fastq "chunks" in the input file: 
    seqs = subprocess.check_output(["grep", "-c", "@", str(args.infile)]) 
    # randomly select n fastq "chunks": 
    seqsamples = random.sample(xrange(0,int(seqs)), int(args.n)) 
    # make a list of the lines that are to be fetched from the fastq file: 
    for i in seqsamples: 
     linesamples.append(int(4*i+0)) 
     linesamples.append(int(4*i+1)) 
     linesamples.append(int(4*i+2)) 
     linesamples.append(int(4*i+3)) 
    # fetch lines from input file and write them to output file. 
    for i, line in enumerate(infile): 
     if i in linesamples: 
      outputfile.write(line) 

grep-krok trwa praktycznie w ogóle czasu, ale po ponad 500 minut, skrypt wciąż nie zaczął pisać do pliku wyjściowego. Więc przypuszczam, że jest to jeden z kroków pomiędzy grep a ostatnią pętlą for, która trwa tak długo. Ale nie rozumiem, który krok dokładnie i co mogę zrobić, aby przyspieszyć.

+7

Powinieneś rozważyć [profilowanie] (https://docs.python.org/2/library/profile.html) swojego programu, aby zobaczyć, które kroki się zawiesią.Czy próbowałeś też uruchomić swój kod na mniejszym pliku i sprawdzić, czy działa on do końca? Kolejnym krokiem, który rozważam później w optymalizacji, jest wykorzystanie wątków i procesów wieloprocesorowych do podziału pracy. –

+0

Nie wywoływać 'int' przez cały czas w pętli. Użyj również instrukcji 'with'. – Veedrac

Odpowiedz

2

W zależności od rozmiaru linesamples, od czasu, gdy przeszukujesz listę dla każdej iteracji, , przez infile. Możesz przekonwertować to na set, aby poprawić czas wyszukiwania. Ponadto, enumerate nie jest bardzo wydajny - zastąpiłem go konstrukcją line_num, którą zwiększamy w każdej iteracji.

def sample(): 
    linesamples = set() 
    infile = open(args.infile, 'r') 
    outputfile = open(args.outputfile, 'w') 
    # count the number of fastq "chunks" in the input file: 
    seqs = subprocess.check_output(["grep", "-c", "@", str(args.infile)]) 
    # randomly select n fastq "chunks": 
    seqsamples = random.sample(xrange(0,int(seqs)), int(args.n)) 
    for i in seqsamples: 
     linesamples.add(int(4*i+0)) 
     linesamples.add(int(4*i+1)) 
     linesamples.add(int(4*i+2)) 
     linesamples.add(int(4*i+3)) 
    # make a list of the lines that are to be fetched from the fastq file: 
    # fetch lines from input file and write them to output file. 
    line_num = 0 
    for line in infile: 
     if line_num in linesamples: 
      outputfile.write(line) 
     line_num += 1 
    outputfile.close() 
+0

Myślę, że ta odpowiedź poprawnie identyfikuje wąskie gardło "jeśli w liniach". Jawne zamknięcie uchwytów plików byłoby prawdopodobnie dobrym pomysłem, ale :) – cel

+0

@cel, uzgodniono - zaktualizowałem kod. – vikramls

+3

'wyliczanie' nie jest wydajne? Czy masz jakieś testy, aby to udowodnić? – Matthias

0

Spróbuj zrównoleglić swój kod. Mam na myśli to. Masz 14 000 000 linii wejściowych.

  1. Praca swoją grep i filtrować wiersze pierwszy i napisz go filteredInput.txt
  2. Splitu swój filteredInput do 10.000-100.000 liniach plików, takich jak filteredInput001.txt, filteredInput002.txt
  3. działa nasz kod na te podzielone pliki. Napisz swoje dane wyjściowe do różnych plików, takich jak output001.txt, output002.txt
  4. Połącz swoje wyniki jako ostatni krok.

Ponieważ Twój kod w ogóle nie działa. Możesz także uruchomić swój kod na tych przefiltrowanych danych wejściowych. Twój kod sprawdzi istnienie plików FilterInput i będzie wiedział, w którym kroku się znajduje, i wznowi od tego kroku.

Można również użyć wielu procesów Pythona w ten sposób (po kroku 1) za pomocą wątków powłoki lub python.

+1

Prawdopodobnie nie jest dobrym pomysłem zasugerowanie porównania przed optymalizacją algorytmu. Przy odpowiednim algorytmie IO będzie szyjką butelki, a nie procesorem. – cel

+0

@cel Jego kod nie działa w tej chwili, ale podzielenie problemu i parowanie nie jest dobrym pomysłem. –

1

Mówiłeś, że grep kończy działa dość szybko, więc w tym przypadku, a nie tylko przy użyciu grep do zliczania wystąpień @ mieć wyjście grep Przesunięcia bajcie każdego znaku @ to widzi (za pomocą opcji -b dla grep). Następnie użyj random.sample, aby wybrać blok, który chcesz zablokować. Po wybraniu przesunięć bajtów, które chcesz, użyj infile.seek, aby przejść do każdego przesunięcia bajtu i wydrukować 4 linie z tego miejsca.

0

Można użyć algorytmu Reservoir Sampling. Za pomocą tego algorytmu można odczytać dane tylko raz (nie trzeba z góry liczyć wierszy pliku), aby można było przesyłać dane przez skrypt. Istnieje przykładowy kod Pythona na stronie Wikipedii.

Istnieje również implementacja C do szybkiego próbkowania w Heng Li's seqtk.

Powiązane problemy