String vergleiche und Zählen von sehr großer Datenmenge

[nork]

Moin @ll,

ich habe gerade mit der Python Programmierung angefangen und habe gleich ein Performance Problem.

Ich habe eine riesige Datei mit mehr als 18.526.903 Zeilen. Nun möchte ich diese Datei aufschlüsseln und auszählen.
In dieser Datei sind Basen, Chromosomnamen und die Coverage der Base festgehalten. Ich möchte nun ein Histogramm für die Coverage erstellen, die Verteilung der Basen und die Verteilung der Basen pro Chromosom.

Die gute Nachricht ist das es soweit funktioniert nur leider ewig dauert. Das Programm läuft später auf einem Clusterserver und es könnten daher auch mehrere Prozessoren benutzt werden.

Ich hoffe das jemand von euch vielleicht eine Idee hat wie es besser/schneller geht.

Vielen Dank für eure Hilfe!

Nork

Code: Alles auswählen

class Count(object):
    def __init__(self,vcfFile):
        #threading.Thread.__init__(self)
        self.vcfFile = vcfFile
    
    def count(self):
        print(self.vcfFile)
        fileOne = open(self.vcfFile,'rb')
        fileOneOne = open(self.vcfFile,'rb')

        #Den Header in der VCF Datei auslassen
        iSetRowOne = int(SetFirstRow(fileOneOne))

        #Die Spaltennamen aus der VCF Datei holen
        readerVCF = csv.reader(fileOne, delimiter="\t")
        for i in range(iSetRowOne):
            strTitlesFirst = readerVCF.next()

        readerVCFtoDict = csv.DictReader(fileOne, strTitlesFirst, delimiter="\t")
        
        #VCF Datei Zeile für Zeile durchgehen.
        try:
            for rowBase in readerVCFtoDict:
                #Die Zeile Info enthält mehrere Informationen und diese Spalte ich auf
                dictData = rowBase['INFO'][:].strip()
                if rowBase['ALT'] == 'X':
                    if rowBase['REF'] != 'N':
                        #Statistiken in die jeweiligen Dict schreiben
                        ddRowBaseInfo = constructStrToDict(dictData,';','=')                 
                        dictHistBase[rowBase['REF']] += 1
                        dictHistCoverage[int(ddRowBaseInfo['DP'])] += 1
                        dictHistCoveragePerChr[rowBase['#CHROM']][int(ddRowBaseInfo['DP'])] += 1
                        iMergeTrue += 1
        finally:
            print("Statistik generiert!")

Code: Alles auswählen

def constructStrToDict(strInput,seperatorOne,seperatorTwo):
    ddInfoRow = {}
    strItems = strInput.split(seperatorOne)
    for item in strItems:
        #print(item)
        a = 0
        a = item.find(seperatorTwo)
        if a > 0 :
            key,value = item.split(seperatorTwo)
        else:
            key,value = item.split('I')
        ddInfoRow[key] = value
    return ddInfoRow

Bei Fragen stehe natürlich zur Verfügung.

[sma]

Ein paar allgemeine Anmerkungen, ohne das ich den Tipp geben kann, wie's sofort schneller wird.

Ein CPython-Programm kann immer nur einen Prozesserkern nutzen, mehr CPUs bringen dir leider nichts, wenn du dein Problem nicht in mehrere Programme parallelisieren kannst. Paralleles zählen und auswerten klingt nach map/reduce, doch ob's da was fertiges für Python gibt, weiß ich nicht. Die große Java-Keule wäre vielleicht Hadoop oder so.

Doch zum Python-Stil: Dein "vcfFile" ist eher ein "vcfFilename", kein File-Objekt. Warum öffnest du es zweimal? Und warum binär? Das ist doch eine Textdatei, oder? Allgemein finde ich die Idee nicht gut, große Dateien komplett in den Hauptspeicher zu saugen. Lade und verarbeite lieber Zeile für Zeile. Das mag aber kein Problem bei heutigen Rechnern sein, wo ich problemlos mehrere GB im RAM halten kann. Ob der eingebaute CSV-Reader auch streamen kann, müsste ich nachgucken, im zweifel sollte es doch aber auch so möglich sein, eine Datei Zeile für Zeile an allen TABs zu splitten. Benutze nicht range() für irgend etwas großes. Oder ist iSetRowOne nicht groß? Was macht denn SetFirstRow() mit der zweiten Datei? Sollte es viel laden, mache doch die Datei auch wieder zu, damit das Betriebssystem Speicher freigeben kann. Ressourcenverwaltung ist deine Aufgabe, nutze sie so sparsam wie möglich. Der "with"-Befehl bei Python hilft bei Dateien, das schließen nicht zu vergessen. Das try-finally sorgt nur dafür, dass auch bei Fehlern die Erfolgsmeldung ausgegeben wird. Willst du das? Will man diese nicht nur sehen, wenn's durchgelaufen ist? Ich vestehe so leider gar nicht, was der eigentliche Algorithmus macht. Namen wie "dictData" oder "ddRowBaseInfo" helfen da auch nicht wirklich. Der Datentyp (siehe auch "iSetRowOne") hat eigentlich nie etwas im Namen einer Variablen zu suchen. Benenne sie danach, was sie für eine Aufgabe hat.

Ich würde sagen, der Algorithmus ist O(N) und nicht O(N^2) oder so, d.h. da ist wenig Optimierungspotential. Du benutzt viele Dictionaries, wo Exemplare von Klassen vielleicht besser (und effizienter) wären, aber das hilft wohl nicht, um bei 18 Mio Zeilen das ganze wirklich zu beschleunigen. Ich würde meinen, dass [:] zum Kopieren vor dem strip() ist überflüssig, aber ansonsten überlasse ich die Analyse mal jemandem, der sich vielleicht damit auskennt, worum es bei Chromosomen und Basen geht.

Stefan

[pillmuncher]

@nork: Zeig doch mal einen (aussagekräftigen) Auszug aus deinem File.

[BlackJack]

@nork: Du bist mit Informationen ein wenig geizig. Die Gesamtgrösse der Datei und eine genauere Zeitangabe als „ewig” wären Beispielsweise nützlich das mal auf konkrete Zahlen von Zeilen und Bytes pro Zeit umzurechnen.

Ansonsten gilt bei Optimierungen eigentlich immer, dass man erst einmal messen sollte wo die meiste Zeit verbraucht wird, damit man nicht an Stellen schraubt, die kaum bis keinen Einfluss auf die Geschwindigkeit haben. Bei solchen Aktionen, die grosse Dateien verarbeiten, sollte man zusätzlich die Zeit bestimmen, die ein Einlesen ohne weitere Verarbeitung braucht um den Anteil zu erfahren, den man sowieso nicht weiter beeinflussen kann.

Warum ist `Count` eine Klasse und keine Funktion? Wo kommen `dictHistBase`, `dictHistCoverage`, `dictHistCoveragePerChr`, und `iMergeTrue` her? `iMergeTrue` wird bei dem Quelltext ausserdem zu einem `UnboundLocalError` führen, das kann also nicht das sein, was Du tatsächlich verwendest.

Ansonsten sind offensichtliche Sachen zusätzlich von den Anmerkungen von sma IMHO:

`dictData` wird für jeden Datensatz berechnet und gebunden, egal ob der Datensatz nun verwendet wird oder nicht.

Die beiden folgenden ``if``-Anweisungen kann man zusammen fassen. Den numerischen Wert vom 'DP'-Schlüssel kann man auch nur einmal berechnen. Beziehungsweise könnte man `dictHistCoverage` am Ende aus `dictHistCoveragePerChr` berechnen.

Die Zuweisung ``a = 0`` in `constructStrToDict()` ist überflüssig. Der ``in``-Operator wäre vielleicht auch schneller als die `find()`-Methode.

Ungetestet:

Code: Alles auswählen

import csv
from collections import defaultdict
from functools import partial
from itertools import islice

# ...

def count(vcf_filename):
    header_length = get_header_length(vcf_filename)
    vcf_file = open(vcf_filename, 'rb')
    for _ in islice(csv.reader(vcf_file, delimiter='\t'), header_length):
        pass

    rows = (
        row for row in csv.DictReader(vcf_file, delimiter='\t')
        if row['ALT'] == 'X' and row['REF'] != 'N'
    )
    base_histogram = defaultdict(int)
    coverage_per_chromosome_histogram = defaultdict(partial(defaultdict, int))
    for row in rows:
        dp_value = int(parse_info(row['INFO'])['DP'])
        base_histogram[row['REF']] += 1
        coverage_per_chromosome_histogram[row['#CHROM']][dp_value] += 1

    coverage_histogram = defaultdict(int)
    for histogram in coverage_per_chromosome_histogram.itervalues():
        for key, value in histogram.iteritems():
            coverage_histogram[key] += value
    
    return (
        sum(base_histogram.itervalues()),
        base_histogram,
        coverage_histogram,
        coverage_per_chromosome_histogram
    )


def parse_info(string, separator_a=';', separator_b='='):
    return dict(
        item.split((separator_b if separator_b in item else 'I'), 1)
        for item in string.split(separator_a)
    )

[nork]

Moin,

danke sma für eine erste Analyse. Bei der Variablennamen vergabe hab ich sicher noch nachhol Potential, aber sein altes Schema raus zu bekommen ist nicht so einfach. Wobei ich gerne den Datentyp vor eine Variable setzte um zu wissen was da eigentlich drin sein soll.

Ich hatte über Zeilenweises einlesen nachgedacht, aber mich dagegen entschieden weil ich dachte die Daten aus dem Ram zu ziehen schneller sei als jedes mal die Festplatte zu bemühen. Der Clusterserver hat um die 64gb und die Datei ist ~10gb groß sollte also kein Problem sein.

@pillmuncher

Code: Alles auswählen

##fileformat=VCFv4.1
##samtoolsVersion=0.1.18 (r982:295)
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	1_Testrun
chr1	114813	.	G	X	0	.	DP=1;I16=1,0,0,0,36,1296,0,0,0,0,0,0,13,169,0,0	PL	0,3,4
chr1	572208	.	T	X	0	.	DP=1;I16=1,0,0,0,41,1681,0,0,1,1,0,0,22,484,0,0	PL	0,3,4

so sehen die Dateien aus.

Mit freundlichen Grüßen
Norman

[nork]

@nork: Du bist mit Informationen ein wenig geizig. Die Gesamtgrösse der Datei und eine genauere

Also die Größe beträgt 9,2 GB (9889235283 Bytes). Die Größe mag aber potentiel eher größer sein.

Zeitangabe als „ewig” wären Beispielsweise nützlich das mal auf konkrete Zahlen von Zeilen und Bytes pro Zeit umzurechnen.

Ich habe das Skript auf meinem Rechner noch nicht durchlaufen lassen. Werde es heute über Nacht mal versuchen. Dauert aber mehr als 3Stunden zu diesem Zeitpunkt brach ich ab.

Ansonsten gilt bei Optimierungen eigentlich immer, dass man erst einmal messen sollte wo die meiste Zeit verbraucht wird, damit man nicht an Stellen schraubt, die kaum bis keinen Einfluss auf die Geschwindigkeit haben. Bei solchen Aktionen, die grosse Dateien verarbeiten, sollte man zusätzlich die Zeit bestimmen, die ein Einlesen ohne weitere Verarbeitung braucht um den Anteil zu erfahren, den man sowieso nicht weiter beeinflussen kann.

Danke werde dann gleich mal Testen.

Warum ist `Count` eine Klasse und keine Funktion? Wo kommen `dictHistBase`, `dictHistCoverage`, `dictHistCoveragePerChr`, und `iMergeTrue` her? `iMergeTrue` wird bei dem Quelltext ausserdem zu einem `UnboundLocalError` führen, das kann also nicht das sein, was Du tatsächlich verwendest.

Entschuldige die habe ich beim einfügen wohl oder bearbeiten im Forum wohl verschwinden lassen.

`dictData` wird für jeden Datensatz berechnet und gebunden, egal ob der Datensatz nun verwendet wird oder nicht.

Jeder Datensatz hat auch eine Info Zeile die ich zerlegen möchte.

Die beiden folgenden ``if``-Anweisungen kann man zusammen fassen.

Stimmt.

Den numerischen Wert vom 'DP'-Schlüssel kann man auch nur einmal berechnen. Beziehungsweise könnte man `dictHistCoverage` am Ende aus `dictHistCoveragePerChr` berechnen.

Der DP-Wert wechselt von Zeile zu Zeile und schwankt zwischen 0-10000. Wüsste jetzt nicht direkt wie ich den Wert vorher berechnen lassen sollte.

Die Zuweisung ``a = 0`` in `constructStrToDict()` ist überflüssig. Der ``in``-Operator wäre vielleicht auch schneller als die `find()`-Methode.

Ungetestet:

Code: Alles auswählen

import csv
from collections import defaultdict
from functools import partial
from itertools import islice

# ...

def count(vcf_filename):
    header_length = get_header_length(vcf_filename)
    vcf_file = open(vcf_filename, 'rb')
    for _ in islice(csv.reader(vcf_file, delimiter='\t'), header_length):
        pass

    rows = (
        row for row in csv.DictReader(vcf_file, delimiter='\t')
        if row['ALT'] == 'X' and row['REF'] != 'N'
    )
    base_histogram = defaultdict(int)
    coverage_per_chromosome_histogram = defaultdict(partial(defaultdict, int))
    for row in rows:
        dp_value = int(parse_info(row['INFO'])['DP'])
        base_histogram[row['REF']] += 1
        coverage_per_chromosome_histogram[row['#CHROM']][dp_value] += 1

    coverage_histogram = defaultdict(int)
    for histogram in coverage_per_chromosome_histogram.itervalues():
        for key, value in histogram.iteritems():
            coverage_histogram[key] += value
    
    return (
        sum(base_histogram.itervalues()),
        base_histogram,
        coverage_histogram,
        coverage_per_chromosome_histogram
    )


def parse_info(string, separator_a=';', separator_b='='):
    return dict(
        item.split((separator_b if separator_b in item else 'I'), 1)
        for item in string.split(separator_a)
    )

Danke für die vielen Hinweise werde es gleich mal testen.

Gruß
Norman

[BlackJack]

@nork: Du verwechselst die Begriffe Zeile und Spalte anscheinend teilweise. Du hast keine Info-*Zeile*, sondern eine Info-*Spalte*. Und den Wert in dieser Spalte zerlegt Dein Code für *jeden* Datensatz; auch für solche, die Du dann gar nicht weiter verarbeitest. Das ist unnötige Arbeit.

Beim 'DP'-Schlüssel meinte ich einmal pro Datensatz ermitteln und nicht zweimal. Denn innerhalb eines Datensatzes bleibt der Wert ja gleich.

[nork]

Code: Alles auswählen

Anweisung: ('countBlackJack', {'vcfFile': '2011_08_19_Testrun_01-1-Idx_11-7.bam/2011_08_19_Testrun_01-1-Idx_11-7.bam.small.vcf'})
Programm beendet in 0.0160757025083 min.
Anweisung: ('countNork', {'vcfFile': '2011_08_19_Testrun_01-1-Idx_11-7.bam/2011_08_19_Testrun_01-1-Idx_11-7.bam.small.vcf'})
Programm beendet in 0.0105002363523 min.

Ich habe die Time Methode verwendet hoffe die ist genau genug. Beide Codes sind unverändert werde jetzt mal versuchen das beste aus beiden zu vereinen. @BlackJack der Code ist auf jeden Fall übersichtlicher. Vielen Dank für den Hinweis!

Die Test für die unveränderlichen Faktoren steht noch aus aber werde das noch machen.

http://bioinet.org/2011_08_19_Testrun_0 ... .small.vcf
Ich habe eine verkleinerte Datei (nur jede 10000 Zeile) erstellt zum testen.

Vielen Dank schon mal soweit.

Gruß
Norman