pyBigWig
0.3.23
Eine in C geschriebene Python -Erweiterung, um einen schnellen Zugriff auf Bigbed -Dateien und den Zugriff auf und auf die Erstellung von Bigwig -Dateien zu erstellen. Diese Erweiterung verwendet Libbigwig für den Zugriff auf lokale und Remotedatei.
Sie können diese Erweiterung direkt von GitHub mit:
pip install pyBigWig
oder mit Conda
conda install pybigwig -c conda-forge -c bioconda
Die folgenden Nicht-Python-Anforderungen müssen installiert werden:
curl-config -Konfiguration)Die Header und Bibliotheken für diese sind erforderlich.
Die Grundnutzung ist wie folgt:
>>> import pyBigWig
Dies funktioniert, wenn Ihr Arbeitsverzeichnis das Pybigwig -Quellcode -Verzeichnis ist.
>>> bw = pyBigWig.open("test/test.bw")
Beachten Sie, dass, wenn die Datei nicht vorhanden ist, eine Fehlermeldung angezeigt wird und None zurückgegeben wird. Seien Sie standardmäßig, alle Dateien werden zum Lesen und nicht zum Schreiben geöffnet. Sie können dies ändern, indem Sie einen Modus mit w : W:
>>> bw = pyBigWig.open("test/output.bw", "w")
Beachten Sie, dass eine zum Schreiben geöffnete Datei nicht für ihre Intervalle oder Statistiken abgefragt werden kann, sie kann nur geschrieben werden. Wenn Sie eine Datei zum Schreiben öffnen, müssen Sie als nächstes einen Header hinzufügen (siehe Abschnitt unten).
Lokaler und abgelegener Bigbed Read Access wird ebenfalls unterstützt:
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
Während Sie einen Modus für Bigbed -Dateien angeben können, wird er ignoriert. Das von pyBigWig.open() zurückgegebene Objekt ist das gleiche, unabhängig davon, ob Sie eine große oder große Datei eröffnen.
Da beide Bigwig- und Bigbed -Dateien geöffnet werden können, kann es erforderlich sein, festzustellen, ob ein bestimmtes bigWigFile -Objekt auf eine BigWig- oder Bigbed -Datei verweist. Zu diesem Zweck kann man die Funktionen isBigWig() und isBigBed() verwenden:
>>> bw = pyBigWig.open("test/test.bw")
>>> bw.isBigWig()
True
>>> bw.isBigBed()
False
bigWigFile -Objekte enthalten ein Wörterbuch, das die Chromosomenlängen hält und mit dem chroms() -Accessor zugegriffen werden kann.
>>> bw.chroms()
dict_proxy({'1': 195471971L, '10': 130694993L})
Sie können auch ein bestimmtes Chromosom direkt abfragen.
>>> bw.chroms("1")
195471971L
Die Längen werden den "langen" Ganzzahltyp gespeichert, weshalb es ein L -Suffix gibt. Wenn Sie ein nicht existierendes Chromosom angeben, wird nichts ausgegeben.
>>> bw.chroms("c")
>>>
Es ist manchmal nützlich, einen Bigwig -Header zu drucken. Dies wird hier als Python -Wörterbuch dargestellt, das: die Version (typischerweise 4 ), die Anzahl der Zoomniveaus ( nLevels ), die Anzahl der beschriebenen Grundlagen ( nBasesCovered ), der Mindestwert ( minVal ), der Maximalwert ( maxVal ), die Summe aller Werte ( sumData ) und die Summe aller quadratischen Werte ( sumSquared ). Die letzten beiden davon werden benötigt, um den Mittelwert und die Standardabweichung zu bestimmen.
>>> bw.header()
{'maxVal': 2L, 'sumData': 272L, 'minVal': 0L, 'version': 4L, 'sumSquared': 500L, 'nLevels': 1L, 'nBasesCovered': 154L}
Beachten Sie, dass dies auch für Bigbed -Dateien möglich ist und dieselben Wörterbuchschlüssel vorhanden sind. Einträge wie maxVal , sumData , minVal und sumSquared sind dann weitgehend nicht sinnvoll.
Bigwig -Dateien werden verwendet, um Werte zu speichern, die mit Positionen und Bereichen von ihnen verbunden sind. Normalerweise möchten wir schnell auf den Durchschnittswert über einen Bereich zugreifen, was sehr einfach ist:
>>> bw.stats("1", 0, 3)
[0.2000000054637591]
Nehmen wir anstelle des Mittelwerts an, wir wollten stattdessen den Maximalwert:
>>> bw.stats("1", 0, 3, type="max")
[0.30000001192092896]
Andere Optionen sind "min" (der Mindestwert), "Abdeckung" (der Bruch der abgedeckten Basen) und "std" (die Standardabweichung der Werte).
Es ist oft der Fall, dass wir stattdessen die Werte einer Reihe von gleichmäßig verteilten Behältern in einem bestimmten Intervall berechnen möchten, was ebenfalls einfach ist:
>>> bw.stats("1",99, 200, type="max", nBins=2)
[1.399999976158142, 1.5]
nBins standardmäßig 1, so wie es type standardmäßig mean .
Wenn die Start- und Endpositionen weggelassen werden, wird das gesamte Chromosom verwendet:
>>> bw.stats("1")
[1.3351851569281683]
Ein Hinweis für den Laienleser: Dieser Abschnitt ist eher technisch und nur aus Vollständigkeit enthalten. Die Zusammenfassung ist, dass, wenn Ihre Bedürfnisse einen genauen Mittelwert/max/etc erfordern. Zusammenfassende Werte für ein Intervall oder ein Intervalle und ein kleiner Kompromiss in der Geschwindigkeit ist akzeptabel, dass Sie die
exact=Truein der Funktion derstats()verwenden sollten.
Standardmäßig gibt es einige unintuitive Aspekte für die Berechnung von Statistiken zu Bereichen in einer Bigwig -Datei. Das Bigwig -Format wurde ursprünglich im Kontext von Genombrowsern erstellt. Dort ist die Berechnung der exakten Zusammenfassungsstatistiken für ein bestimmtes Intervall weniger wichtig als schnell in der Lage zu sein, eine ungefähre Statistik zu berechnen (schließlich müssen die Browser in der Lage sein, schnell eine Reihe von zusammenhängenden Intervallen anzuzeigen und Scrolling/Zooming zu unterstützen). Aus diesem Grund enthalten Bigwig-Dateien nicht nur Intervallwert-Assoziationen, sondern auch sum of values / sum of squared values / minimum value / maximum value / number of bases covered . Diese verschiedenen Größen werden als "Zoomebenen" bezeichnet. Die kleinste Zoomebene hat Mülleimer, die 16 -mal so groß sind wie die mittlere Intervallgröße in der Datei, und jede nachfolgende Zoomebene hat die viermal größeren Behälter als die vorherige. Diese Methodik wird in Kents Tools verwendet und daher wahrscheinlich in fast jeder derzeit vorhandenen Bigwig -Datei verwendet.
Wenn eine BigWig -Datei für eine zusammenfassende Statistik abgefragt wird, wird die Größe des Intervalls verwendet, um festzustellen, ob eine Zoomebene verwendet werden soll, und wenn ja, welches. Der optimale Zoom -Niveau ist das, was die größten Mülleimer nicht mehr als die Hälfte der Breite des gewünschten Intervalls aufweist. Wenn keine solche Zoomebene vorhanden ist, werden die ursprünglichen Intervalle stattdessen für die Berechnung verwendet.
Aus Gründen der Konsistenz mit anderen Tools nimmt Pybigwig dieselbe Methodik an. Da dies jedoch (a) in einigen Anwendungen unintuitiv und (b) unerwünscht ist, ermöglicht PyBigwig die Berechnung der exakten Zusammenfassungsstatistiken unabhängig von der Intervallgröße (dh ermöglicht das Ignorieren der Zoomniveaus). Dies wurde ursprünglich hier vorgeschlagen, und ein Beispiel ist unten:
>>> import pyBigWig
>>> from numpy import mean
>>> bw = pyBigWig.open("http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeCrgMapabilityAlign75mer.bigWig")
>>> bw.stats('chr1', 89294, 91629)
[0.20120902053804418]
>>> mean(bw.values('chr1', 89294, 91629))
0.22213841940688142
>>> bw.stats('chr1', 89294, 91629, exact=True)
[0.22213841940688142]
Während die stats() -Methode verwendet werden kann , um die ursprünglichen Werte für jede Basis (z. B. durch Einstellen nBins auf die Anzahl der Basen) abzurufen, ist es vorzuziehen, stattdessen den values() -Accessor zu verwenden.
>>> bw.values("1", 0, 3)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896]
Die erstellte Liste enthält immer einen Wert für jede Basis im angegebenen Bereich. Wenn eine bestimmte Basis keinen Wert in der BigWig -Datei hat, ist der zurückgegebene Wert nan .
>>> bw.values("1", 0, 4)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896, nan]
Manchmal ist es bequem, alle Einträge abzurufen, die einen Bereich überlappen. Dies kann mit der Funktion intervals() erfolgen:
>>> bw.intervals("1", 0, 3)
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896))
Was zurückgegeben wird, ist eine Liste von Tupeln, die enthalten sind: die Startposition, die End -End -Position und der Wert. Somit hat das obige Beispiel Werte von 0.1 , 0.2 und 0.3 an den Positionen 0 , 1 bzw. 2 .
Wenn die Start- und Endposition weggelassen wird, werden alle Intervalle auf dem angegebenen Chromosom zurückgegeben:
>>> bw.intervals("1")
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896), (100, 150, 1.399999976158142), (150, 151, 1.5))
Im Gegensatz zu Bigwig -Dateien enthalten Bigbed -Dateien Einträge, bei denen es sich um Intervalle mit einer zugehörigen Zeichenfolge handelt. Sie können diese Einträge mit der Funktion entries() zugreifen:
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
>>> bb.entries('chr1', 10000000, 10020000)
[(10009333, 10009640, '61035t130t-t0.026t0.42t404'), (10014007, 10014289, '61047t136t-t0.029t0.42t404'), (10014373, 10024307, '61048t630t-t5.420t0.00t2672399')]
Die Ausgabe ist eine Liste von Eintragstupeln. Die Tupelelemente sind die start und end jedes Eintrags, gefolgt von der zugehörigen string . Die Zeichenfolge wird genau so zurückgegeben, wie sie in der Bigbed -Datei aufbewahrt wird, sodass es Ihnen überlassen bleibt. Um festzustellen, was die verschiedenen Felder in dieser Zeichenfolge sind, wenden Sie sich an die SQL -Zeichenfolge:
>>> bb.SQL()
table RnaElements
"BED6 + 3 scores for RNA Elements data"
(
string chrom; "Reference sequence chromosome or scaffold"
uint chromStart; "Start position in chromosome"
uint chromEnd; "End position in chromosome"
string name; "Name of item"
uint score; "Normalized score from 0-1000"
char[1] strand; "+ or - or . for unknown"
float level; "Expression level such as RPKM or FPKM. Set to -1 for no data."
float signif; "Statistical significance such as IDR. Set to -1 for no data."
uint score2; "Additional measurement/count e.g. number of reads. Set to 0 for no data."
)
Beachten Sie, dass die ersten drei Einträge in der SQL -Zeichenfolge nicht Teil der Zeichenfolge sind.
Wenn Sie nur wissen müssen, wo Einträge sind und nicht ihre zugehörigen Werte, können Sie Speicher speichern, indem Sie zusätzlich withString=False in entries() angeben:
>>> bb.entries('chr1', 10000000, 10020000, withString=False)
[(10009333, 10009640), (10014007, 10014289), (10014373, 10024307)]
Wenn Sie eine Datei zum Schreiben geöffnet haben, müssen Sie ihm einen Kopfball geben, bevor Sie Einträge hinzufügen können. Der Header enthält alle Chromosomen in Ordnung und ihre Größen. Wenn Ihr Genom zwei Chromosomen, CHR1 und CHR2, mit den Längen 1 und 1,5 Millionen Basen hat, würde Folgendes einen geeigneten Kopfball hinzufügen:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)])
Bigwig-Header sind fälschlichem sensitiv, also sind chr1 und Chr1 unterschiedlich. Ebenso sind 1 und chr1 nicht gleich, sodass Sie Ensembl- und UCSC -Chromosomennamen nicht mischen können. Nach dem Hinzufügen eines Headers können Sie Einträge hinzufügen.
Standardmäßig werden bis zu 10 "Zoomebenen" für Bigwig -Dateien konstruiert. Sie können diese Standardnummer mit dem optionalen Argument maxZooms ändern. Eine häufige Verwendung davon besteht darin, eine BigWig -Datei zu erstellen, die einfach Intervalle und keine Zoomebenen enthält:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)], maxZooms=0)
Wenn Sie maxTooms=0 festlegen, beachten Sie bitte, dass IGV und viele andere Tools nicht funktionieren, da sie annehmen, dass mindestens eine Zoomebene vorhanden ist. Sie werden empfohlen, den Standard zu verwenden, es sei denn, Sie erwarten nicht, dass die BigWig -Dateien von anderen Paketen verwendet werden.
Angenommen, Sie haben eine Datei zum Schreiben geöffnet und einen Header hinzugefügt, können Sie dann Einträge hinzufügen. Beachten Sie, dass die Einträge in der Reihenfolge hinzugefügt werden müssen , da Bigwig -Dateien immer geordnete Intervalle enthalten. Es gibt drei Formate, die BigWig -Dateien intern verwenden können, um Einträge zu speichern. Das am häufigsten beobachtete Format ist identisch mit einer Bedgraph -Datei:
chr1 0 100 0.0
chr1 100 120 1.0
chr1 125 126 200.0
Diese Einträge würden wie folgt hinzugefügt:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [0, 100, 125], ends=[5, 120, 126], values=[0.0, 1.0, 200.0])
Jeder Eintrag nimmt 12 Bytes vor der Komprimierung ein.
Das zweite Format verwendet eine feste Span, aber eine variable Schrittgröße zwischen Einträgen. Diese können in einer Wackeldatei dargestellt werden wie:
variableStep chrom=chr1 span=20
500 -2.0
600 150.0
635 25.0
Die obigen Einträge beschreiben (1 basierte) Positionen 501-520, 601-620 und 636-655. Diese würden wie folgt hinzugefügt:
>>> bw.addEntries("chr1", [500, 600, 635], values=[-2.0, 150.0, 25.0], span=20)
Jeder Eintrag dieses Typs befindet sich vor der Komprimierung 8 Bytes.
Das endgültige Format verwendet für jeden Eintrag einen festen Schritt und eine Spanne, die dem FixedStep -Wiggle -Format entspricht:
fixedStep chrom=chr1 step=30 span=20
-5.0
-20.0
25.0
Die obigen Einträge beschreiben (1 basierte) Basen 901-920, 931-950 und 961-980 und würden wie folgt hinzugefügt:
>>> bw.addEntries("chr1", 900, values=[-5.0, -20.0, 25.0], span=20, step=30)
Jeder Eintrag dieses Typs befindet sich 4 Bytes.
Beachten Sie, dass Pybigwig versucht, Sie daran zu hindern, Einträge in einer falschen Reihenfolge hinzuzufügen. Dies erfordert jedoch zusätzliche Überköpfe. Sollte dies nicht akzeptabel sein, können Sie einfach validate=False angeben, wenn Sie Einträge hinzufügen:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [100, 0, 125], ends=[120, 5, 126], values=[0.0, 1.0, 200.0], validate=False)
Sie sind dann offensichtlich dafür verantwortlich, dass Sie keine Einträge außerhalb der Reihenfolge hinzufügen. Die resultierenden Dateien wären ansonsten nicht verwendbar.
Eine Datei kann mit einem einfachen bw.close() geschlossen werden, wie dies häufig mit anderen Dateitypen erfolgt. Für zum Schreiben geöffnete Dateien schreibt das Schließen einer Datei alle gepufferten Einträge in die Festplatte, konstruiert und schreibt den Dateiindex und konstruiert Zoomebenen. Folglich kann dies etwas Zeit dauern.
Ab Version 0.3.0 unterstützt Pybigwig die Eingabe von Koordinaten mit Numpy Inters und Vektoren in einigen Funktionen , wenn Numpy vor der Installation von Pybigwig installiert wurde . Um festzustellen, ob Pybigwig mit Numpy Support installiert wurde, indem Sie den numpy Accessor überprüfen:
>>> import pyBigWig
>>> pyBigWig.numpy
1
Wenn pyBigWig.numpy 1 ist, wurde Pybigwig mit Numpy -Unterstützung zusammengestellt. Dies bedeutet, dass addEntries() Numpy -Koordinaten akzeptieren kann:
>>> import pyBigWig
>>> import numpy
>>> bw = pyBigWig.open("/tmp/delete.bw", "w")
>>> bw.addHeader([("1", 1000)], maxZooms=0)
>>> chroms = np.array(["1"] * 10)
>>> starts = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90], dtype=np.int64)
>>> ends = np.array([5, 15, 25, 35, 45, 55, 65, 75, 85, 95], dtype=np.int64)
>>> values0 = np.array(np.random.random_sample(10), dtype=np.float64)
>>> bw.addEntries(chroms, starts, ends=ends, values=values0)
>>> bw.close()
Zusätzlich kann values() einen numpigen Vektor direkt ausgeben:
>>> bw = bw.open("/tmp/delete.bw")
>>> bw.values('1', 0, 10, numpy=True)
[ 0.74336642 0.74336642 0.74336642 0.74336642 0.74336642 nan
nan nan nan nan]
>>> type(bw.values('1', 0, 10, numpy=True))
<type 'numpy.ndarray'>
Wenn Sie keine Curl installiert haben, wird Pybigwig installiert, ohne auf Remotedateien zugreifen zu können. Sie können feststellen, ob Sie mit pyBigWig.remote auf Remotedateien zugreifen können. Wenn dies 1 zurückgibt, können Sie auf Remote -Dateien zugreifen. Wenn es 0 zurückgibt, können Sie nicht.
Ab Version 0.3.5 kann Pybigwig Bigwig -Dateien ohne Einträge lesen und schreiben. Bitte beachten Sie, dass solche Dateien im Allgemeinen nicht mit anderen Programmen kompatibel sind, da keine Definition darüber gibt, wie eine BigWig -Datei ohne Einträge aussehen sollte. Für eine solche Datei gibt der intervals() -Accessor None zurück, die Funktion stats() gibt eine Liste von None der gewünschten Länge zurück, und values() gibt [] (eine leere Liste) zurück. Dies sollte im Allgemeinen ermöglichen, dass Programme, die Pybigwig verwenden, ohne Probleme fortzusetzen.
Für diejenigen, die die Funktionalität von Pybigwig/libbigwig diesbezüglich nachahmen möchten, beachten Sie bitte, dass die Anzahl der abgedeckten Basen (wie im Dateikopf gemeldet) nach "leeren" Dateien untersucht wird.
Wackel-, Bigwig- und Bigbed-Dateien verwenden 0 basierte halbe Open-Koordinaten, die auch von dieser Erweiterung verwendet werden. Um auf den Wert für die erste Basis auf chr1 zuzugreifen, würde man die Startposition als 0 und die Endposition als 1 angeben. In ähnlicher Weise hätten die Basen 100 bis 115 einen Beginn von 99 und ein Ende von 115 . Dies ist lediglich um die Konsistenz mit der zugrunde liegenden Bigwig -Datei und kann sich in Zukunft ändern.
Pybigwig ist auch als Paket in Galaxy erhältlich. Sie finden es im Toolshed und die IUC hostet derzeit die XML -Definition davon auf GitHub.
