pyBigWig
0.3.23
Una extensión de Python, escrita en C, para acceso rápido a archivos Bigbed y acceso y creación de archivos BigWig. Esta extensión utiliza libBigWig para el acceso local y remoto de archivos.
Puede instalar esta extensión directamente desde GitHub con:
pip install pyBigWig
o con conda
conda install pybigwig -c conda-forge -c bioconda
Se deben instalar los siguientes requisitos que no son de Python:
curl-config )Se requieren los encabezados y las bibliotecas para estos.
El uso básico es el siguiente:
>>> import pyBigWig
Esto funcionará si su directorio de trabajo es el directorio del código fuente de PybigWig.
>>> bw = pyBigWig.open("test/test.bw")
Tenga en cuenta que si el archivo no existe, verá un mensaje de error y no se devolverá None . Sea predeterminado, todos los archivos se abren para leer y no escribir. Puede alterar esto pasando un modo que contiene w :
>>> bw = pyBigWig.open("test/output.bw", "w")
Tenga en cuenta que un archivo abierto para la escritura no se puede consultar para sus intervalos o estadísticas, solo se puede escribir. Si abre un archivo para escribir, entonces deberá agregar un encabezado (consulte la sección sobre esto a continuación).
El acceso de lectura Bigbed local y remoto también es compatible:
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
Si bien puede especificar un modo para archivos Bigbed, se ignora. El objeto devuelto por pyBigWig.open() es el mismo independientemente de si está abriendo un archivo BigWig o Bigbed.
Dado que los archivos BigWig y Bigbed se pueden abrir, puede ser necesario determinar si un objeto bigWigFile dado apunta a un archivo BigWig o Bigbed. Para ese fin, se puede usar las funciones isBigWig() e isBigBed() :
>>> bw = pyBigWig.open("test/test.bw")
>>> bw.isBigWig()
True
>>> bw.isBigBed()
False
Los objetos bigWigFile contienen un diccionario que contiene las longitudes del cromosoma, a las que se puede acceder con el accesor chroms() .
>>> bw.chroms()
dict_proxy({'1': 195471971L, '10': 130694993L})
También puede consultar directamente un cromosoma particular.
>>> bw.chroms("1")
195471971L
Las longitudes se almacenan como el tipo entero "largo", por lo que hay un sufijo L Si especifica un cromosoma inexistente, entonces no se emite nada.
>>> bw.chroms("c")
>>>
A veces es útil imprimir un encabezado de Bigwig. Esto se presenta aquí como un diccionario de Python que contiene: la versión (típicamente 4 ), el número de niveles de zoom ( nLevels ), el número de bases descritas ( nBasesCovered ), el valor mínimo ( minVal ), el valor máximo ( maxVal ), el suma de todos los valores ( sumData ), y la suma de todos los valores cuadrados ( sumSquared ). Los dos últimos son necesarios para determinar la media y la desviación estándar.
>>> bw.header()
{'maxVal': 2L, 'sumData': 272L, 'minVal': 0L, 'version': 4L, 'sumSquared': 500L, 'nLevels': 1L, 'nBasesCovered': 154L}
Tenga en cuenta que esto también es posible para archivos grandes y las mismas claves de diccionario estarán presentes. Las entradas como maxVal , sumData , minVal y sumSquared no son en gran medida significativas.
Los archivos BigWig se utilizan para almacenar valores asociados con posiciones y rangos de ellos. Por lo general, queremos acceder rápidamente al valor promedio en un rango, lo cual es muy simple:
>>> bw.stats("1", 0, 3)
[0.2000000054637591]
Supongamos que en lugar del valor medio, en su lugar queríamos el valor máximo:
>>> bw.stats("1", 0, 3, type="max")
[0.30000001192092896]
Otras opciones son "min" (el valor mínimo), "cobertura" (la fracción de bases cubiertas) y "STD" (la desviación estándar de los valores).
A menudo, es el caso del que nos gustaría calcular los valores de un número de contenedores espaciados uniformemente en un intervalo dado, que también es simple:
>>> bw.stats("1",99, 200, type="max", nBins=2)
[1.399999976158142, 1.5]
nBins predeterminado es 1, al igual que type se mean .
Si se omiten las posiciones de inicio y final, se usa todo el cromosoma:
>>> bw.stats("1")
[1.3351851569281683]
Una nota para el lector laico: esta sección es bastante técnica e incluida solo en aras de la integridad. El resumen es que si sus necesidades requieren una media exacta/max/etc. Valores resumidos para un intervalo o intervalos y que una pequeña compensación en la velocidad es aceptable, que debe usar la opción
exact=Trueen la funciónstats().
Por defecto, hay algunos aspectos no intuitivos para calcular estadísticas sobre rangos en un archivo BigWig. El formato BigWig se creó originalmente en el contexto de los navegadores del genoma. Allí, calcular estadísticas de resumen exacta para un intervalo dado es menos importante que poder calcular rápidamente una estadística aproximada (después de todo, los navegadores deben poder mostrar rápidamente una serie de intervalos contiguos y soportar desplazamiento/zoom). Debido a esto, los archivos BigWig contienen no solo asociaciones de valor de intervalo, sino también sum of values / sum of squared values / minimum value / maximum value / number of bases covered para contenedores de tamaño igualmente de varios tamaños. Estos diferentes tamaños se denominan "niveles de zoom". El nivel de zoom más pequeño tiene contenedores que son 16 veces el tamaño medio del intervalo en el archivo y cada nivel de zoom posterior tiene contenedores 4 veces más grandes que el anterior. Esta metodología se utiliza en las herramientas de Kent y, por lo tanto, probablemente se usa en casi todos los archivos bigwig actualmente existentes.
Cuando se consulta un archivo BigWig para una estadística de resumen, el tamaño del intervalo se usa para determinar si usar un nivel de zoom y, de ser así, cuál. El nivel de zoom óptimo es el que tiene los contenedores más grandes no más de la mitad del ancho del intervalo deseado. Si no existe dicho nivel de zoom, los intervalos originales se utilizan para el cálculo.
En aras de la consistencia con otras herramientas, Pybigwig adopta esta misma metodología. Sin embargo, dado que esto es (a) no intuitivo y (b) indeseable en algunas aplicaciones, Pybigwig permite el cálculo de las estadísticas de resumen exactas independientemente del tamaño del intervalo (es decir, permite ignorar los niveles de zoom). Esto se propuso originalmente aquí y un ejemplo está a continuación:
>>> 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]
Mientras que el método stats() se puede usar para recuperar los valores originales para cada base (por ejemplo, configurando nBins en el número de bases), es preferible usar el accesor values() .
>>> bw.values("1", 0, 3)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896]
La lista producida siempre contendrá un valor para cada base en el rango especificado. Si una base en particular no tiene un valor asociado en el archivo BigWig, el valor devuelto será nan .
>>> bw.values("1", 0, 4)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896, nan]
A veces es conveniente recuperar todas las entradas que superponen algún rango. Esto se puede hacer con la función intervals() :
>>> bw.intervals("1", 0, 3)
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896))
Lo que se devuelve es una lista de tuplas que contienen: la posición de inicio, la posición final final y el valor. Por lo tanto, el ejemplo anterior tiene valores de 0.1 , 0.2 y 0.3 en las posiciones 0 , 1 y 2 , respectivamente.
Si se omiten la posición de inicio y finalización, se devuelven todos los intervalos en el cromosoma especificados:
>>> bw.intervals("1")
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896), (100, 150, 1.399999976158142), (150, 151, 1.5))
A diferencia de los archivos BigWig, los archivos Bigbed contienen entradas, que son intervalos con una cadena asociada. Puede acceder a estas entradas utilizando la función entries() :
>>> 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')]
La salida es una lista de tuplas de entrada. Los elementos de la tupla son la posición start y end de cada entrada, seguido de su string asociada. La cadena se devuelve exactamente como se mantiene en el archivo Bigbed, por lo que lo que la deja se deja. Para determinar cuáles son los diversos campos en estas cadenas, consulte la cadena SQL:
>>> 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."
)
Tenga en cuenta que las primeras tres entradas en la cadena SQL no son parte de la cadena.
Si solo necesita saber dónde están las entradas y no sus valores asociados, puede guardar la memoria especificando adicionalmente withString=False en entries() ::
>>> bb.entries('chr1', 10000000, 10020000, withString=False)
[(10009333, 10009640), (10014007, 10014289), (10014373, 10024307)]
Si ha abierto un archivo para escribir, entonces deberá darle un encabezado antes de poder agregar cualquier entrada. El encabezado contiene todos los cromosomas, en orden y sus tamaños. Si su genoma tiene dos cromosomas, CHR1 y CHR2, de longitudes 1 y 1.5 millones de bases, entonces lo siguiente agregaría un encabezado apropiado:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)])
Los encabezados de BigWig son sensibles a los casos, por lo que chr1 y Chr1 son diferentes. Del mismo modo, 1 y chr1 no son los mismos, por lo que no puede mezclar los nombres de cromosomas EnsemBL y UCSC. Después de agregar un encabezado, puede agregar entradas.
Por defecto, se construyen hasta 10 "niveles de zoom" para archivos BigWig. Puede cambiar este número predeterminado con el argumento opcional maxZooms . Un uso común de esto es crear un archivo BigWig que simplemente contenga intervalos y no hay niveles de zoom:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)], maxZooms=0)
Si establece maxTooms=0 , tenga en cuenta que IGV y muchas otras herramientas no funcionarán, ya que suponen que al menos un nivel de zoom estará presente. Se recomienda utilizar el valor predeterminado a menos que no espere que los archivos BigWig sean utilizados por otros paquetes.
Suponiendo que haya abierto un archivo para escribir y agregó un encabezado, puede agregar entradas. Tenga en cuenta que las entradas deben agregarse en orden, ya que los archivos BigWig siempre contienen intervalos ordenados. Hay tres formatos que los archivos BigWig pueden usar internamente para almacenar entradas. El formato más comúnmente observado es idéntico a un archivo Bedgraph:
chr1 0 100 0.0
chr1 100 120 1.0
chr1 125 126 200.0
Estas entradas se agregarían de la siguiente manera:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [0, 100, 125], ends=[5, 120, 126], values=[0.0, 1.0, 200.0])
Cada entrada ocupa 12 bytes antes de la compresión.
El segundo formato utiliza un tramo fijo, pero un tamaño de paso variable entre las entradas. Estos pueden representarse en un archivo de maniobra como:
variableStep chrom=chr1 span=20
500 -2.0
600 150.0
635 25.0
Las entradas anteriores describen posiciones (basadas en 1) 501-520, 601-620 y 636-655. Estos se agregarían de la siguiente manera:
>>> bw.addEntries("chr1", [500, 600, 635], values=[-2.0, 150.0, 25.0], span=20)
Cada entrada de este tipo ocupa 8 bytes antes de la compresión.
El formato final utiliza un paso y un tramo fijo para cada entrada, correspondiente al formato de manejo fijo:
fixedStep chrom=chr1 step=30 span=20
-5.0
-20.0
25.0
Las entradas anteriores describen bases (basadas en 1) 901-920, 931-950 y 961-980 y se agregarían de la siguiente manera:
>>> bw.addEntries("chr1", 900, values=[-5.0, -20.0, 25.0], span=20, step=30)
Cada entrada de este tipo ocupa 4 bytes.
Tenga en cuenta que Pybigwig intentará evitar que agregue entradas en un orden incorrecto. Esto, sin embargo, requiere una sobrecarga adicional. Si eso no es aceptable, simplemente puede especificar validate=False al agregar entradas:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [100, 0, 125], ends=[120, 5, 126], values=[0.0, 1.0, 200.0], validate=False)
Obviamente, usted es responsable de asegurarse de no agregar entradas fuera de servicio. Los archivos resultantes de otra manera no serían utilizables.
Un archivo se puede cerrar con un simple bw.close() , como se hace comúnmente con otros tipos de archivos. Para los archivos abiertos para escribir, cerrar un archivo escribe cualquier entrada buffada en el disco, construye y escribe el índice de archivos, y construye los niveles de zoom. En consecuencia, esto puede llevar un poco de tiempo.
A partir de la versión 0.3.0, PybigWig admite la entrada de coordenadas que usan enteros y vectores numpy en algunas funciones si Numpy se instaló antes de instalar PybigWig . Para determinar si PybigWig se instaló con soporte Numpy revisando el accesor numpy :
>>> import pyBigWig
>>> pyBigWig.numpy
1
Si pyBigWig.numpy es 1 , entonces Pybigwig se compiló con soporte numpy. Esto significa que addEntries() puede aceptar coordenadas numpy:
>>> 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()
Además, values() pueden emitir directamente un vector numpy:
>>> 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'>
Si no tiene Curl instalado, PyBigWig se instalará sin la capacidad de acceder a archivos remotos. Puede determinar si podrá acceder a archivos remotos con pyBigWig.remote . Si eso devuelve 1, puede acceder a archivos remotos. Si devuelve 0, entonces no puedes.
A partir de la versión 0.3.5, Pybigwig puede leer y escribir archivos BigWig que carecen de entradas. Tenga en cuenta que dichos archivos generalmente no son compatibles con otros programas, ya que no hay una definición de cómo debe verse un archivo BigWig sin entradas. Para dicho archivo, el accesor intervals() no devolverá None , la función stats() no devolverá una lista de None de las longitud deseadas, y values() devolverá [] (una lista vacía). Esto generalmente debería permitir que los programas que utilicen PybigWig continúen sin problemas.
Para aquellos que deseen imitar la funcionalidad de Pybigwig/Libbigwig a este respecto, tenga en cuenta que analiza la cantidad de bases cubiertas (como se informa en el encabezado del archivo) para verificar los archivos "vacíos".
Los archivos de Wiggle, BigWig y Bigbed usan coordenadas medias basadas en 0, que también son utilizadas por esta extensión. Entonces, para acceder al valor para la primera base en chr1 , se especificaría la posición inicial como 0 y la posición final como 1 . Del mismo modo, las bases 100 a 115 tendrían un comienzo de 99 y un final de 115 . Esto es simplemente por consistencia con el archivo BigWig subyacente y puede cambiar en el futuro.
Pybigwig también está disponible como un paquete en Galaxy. Puede encontrarlo en la plataforma de herramientas y el IUC actualmente está alojando la definición XML de esto en GitHub.
