pyBigWig
0.3.23
Uma extensão Python, escrita em C, para acesso rápido a arquivos Bigbed e acesso e criação de arquivos bigwig. Esta extensão usa Libbigwig para acesso a arquivos locais e remotos.
Você pode instalar esta extensão diretamente do GitHub com:
pip install pyBigWig
ou com conda
conda install pybigwig -c conda-forge -c bioconda
Os requisitos de seguidores não-python devem ser instalados:
curl-config )Os cabeçalhos e bibliotecas são necessários.
O uso básico é o seguinte:
>>> import pyBigWig
Isso funcionará se o seu diretório de trabalho for o diretório de código -fonte do Pybigwig.
>>> bw = pyBigWig.open("test/test.bw")
Observe que, se o arquivo não existir, você verá uma mensagem de erro e None será retornado. Seja padrão, todos os arquivos são abertos para leitura e não escrita. Você pode alterar isso passando um modo contendo w :
>>> bw = pyBigWig.open("test/output.bw", "w")
Observe que um arquivo aberto para gravação não pode ser consultado para seus intervalos ou estatísticas, ele só pode ser gravado. Se você abrir um arquivo para gravar, será necessário adicionar um cabeçalho (consulte a seção abaixo).
O acesso de leitura local e remoto também é suportado:
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
Embora você possa especificar um modo para arquivos bigbed, ele é ignorado. O objeto retornado por pyBigWig.open() é o mesmo, independentemente de você estar abrindo um arquivo bigwig ou bigbed.
Como os arquivos bigwig e bigbed podem ser abertos, pode ser necessário determinar se um determinado objeto bigWigFile aponta para um arquivo bigwig ou bigbed. Para esse fim, pode -se usar as funções isBigWig() e isBigBed() :
>>> bw = pyBigWig.open("test/test.bw")
>>> bw.isBigWig()
True
>>> bw.isBigBed()
False
Os objetos bigWigFile contêm um dicionário que mantém os comprimentos do cromossomo, que podem ser acessados com o acessador chroms() .
>>> bw.chroms()
dict_proxy({'1': 195471971L, '10': 130694993L})
Você também pode consultar diretamente um cromossomo específico.
>>> bw.chroms("1")
195471971L
Os comprimentos são armazenados do tipo inteiro "longo", e é por isso que há um sufixo L Se você especificar um cromossomo inexistente, nada será produzido.
>>> bw.chroms("c")
>>>
Às vezes é útil imprimir o cabeçalho de um bigwig. Isso é apresentado aqui como um dicionário python contendo: a versão (normalmente 4 ), o número de níveis de zoom ( nLevels ), o número de bases descritas ( nBasesCovered ), o valor mínimo ( minVal ), o valor máximo ( maxVal ), o soma de todos os valores ( sumData ) e a soma de todos os valores quadrados ( sumSquared ). Os dois últimos são necessários para determinar a média e o desvio padrão.
>>> bw.header()
{'maxVal': 2L, 'sumData': 272L, 'minVal': 0L, 'version': 4L, 'sumSquared': 500L, 'nLevels': 1L, 'nBasesCovered': 154L}
Observe que isso também é possível para arquivos Bigbed e as mesmas teclas de dicionário estarão presentes. Entradas como maxVal , sumData , minVal e sumSquared não são em grande parte significativas.
Os arquivos bigwig são usados para armazenar valores associados a posições e faixas deles. Normalmente, queremos acessar rapidamente o valor médio em um intervalo, o que é muito simples:
>>> bw.stats("1", 0, 3)
[0.2000000054637591]
Suponha que, em vez do valor médio, queríamos o valor máximo:
>>> bw.stats("1", 0, 3, type="max")
[0.30000001192092896]
Outras opções são "min" (o valor mínimo), "cobertura" (a fração das bases cobertas) e "std" (o desvio padrão dos valores).
Muitas vezes, é o caso de que gostaríamos de calcular valores de algum número de caixas uniformemente espaçadas em um determinado intervalo, o que também é simples:
>>> bw.stats("1",99, 200, type="max", nBins=2)
[1.399999976158142, 1.5]
nBins padrão é 1, assim como os padrões type mean .
Se as posições de início e final forem omitidas, todo o cromossomo será usado:
>>> bw.stats("1")
[1.3351851569281683]
Uma nota para o leitor leigo: Esta seção é bastante técnica e incluída apenas por uma questão de completude. O resumo é que, se suas necessidades exigirem uma média exata/max/etc. Valores resumidos para um intervalo ou intervalos e que uma pequena troca de velocidade é aceitável, que você deve usar a opção
exact=Truena funçãostats().
Por padrão, existem alguns aspectos não intuitivos nas estatísticas de computação em intervalos em um arquivo bigwig. O formato Bigwig foi criado originalmente no contexto de navegadores do genoma. Lá, a computação de estatísticas de resumo exata para um determinado intervalo é menos importante do que poder rapidamente calcular uma estatística aproximada (afinal, os navegadores precisam ser capazes de exibir rapidamente vários intervalos contíguos e suportar rolagem/zoom). Por esse motivo, os arquivos bigwig contêm não apenas associações de valor de intervalo, mas também sum of values / sum of squared values / minimum value / maximum value number of bases covered para caixas de tamanho igual de vários tamanhos. Esses tamanhos diferentes são chamados de "níveis de zoom". O menor nível de zoom possui caixas que são 16 vezes o tamanho médio do intervalo no arquivo e cada nível de zoom subsequente possui caixotes 4 vezes maior que o anterior. Essa metodologia é usada nas ferramentas de Kent e, portanto, provavelmente usada em quase todos os arquivos bigwig atualmente existentes.
Quando um arquivo bigwig é consultado para uma estatística resumida, o tamanho do intervalo é usado para determinar se deve usar um nível de zoom e, em caso afirmativo, qual. O nível ideal de zoom é o que possui as maiores caixas não mais que metade da largura do intervalo desejado. Se não existir esse nível de zoom, os intervalos originais serão usados para o cálculo.
Por uma questão de consistência com outras ferramentas, o Pybigwig adota essa mesma metodologia. No entanto, como isso é (a) não intuitivo e (b) indesejável em algumas aplicações, o pybigwig permite a computação de estatísticas de resumo exatas, independentemente do tamanho do intervalo (ou seja, permite ignorar os níveis de zoom). Isso foi proposto originalmente aqui e um exemplo está abaixo:
>>> 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]
Embora o método stats() possa ser usado para recuperar os valores originais para cada base (por exemplo, definindo nBins para o número de bases), é preferível usar o acessador values() .
>>> bw.values("1", 0, 3)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896]
A lista produzida sempre conterá um valor para cada base do intervalo especificado. Se uma base específica não tiver valor associado no arquivo bigwig, o valor retornado será nan .
>>> bw.values("1", 0, 4)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896, nan]
Às vezes, é conveniente recuperar todas as entradas que se sobrepõem a algum alcance. Isso pode ser feito com a função intervals() :
>>> bw.intervals("1", 0, 3)
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896))
O que é retornado é uma lista de tuplas que contêm: a posição inicial, a posição final e o valor. Assim, o exemplo acima possui valores de 0.1 , 0.2 e 0.3 nas posições 0 , 1 e 2 , respectivamente.
Se a posição inicial e final for omitida, todos os intervalos do cromossomo especificados serão retornados:
>>> bw.intervals("1")
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896), (100, 150, 1.399999976158142), (150, 151, 1.5))
Ao contrário dos arquivos bigwig, os arquivos bigbed mantêm entradas, que são intervalos com uma string associada. Você pode acessar essas entradas usando a função 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')]
A saída é uma lista de tuplas de entrada. Os elementos da tupla são a posição start e end de cada entrada, seguida por sua string associada. A string é retornada exatamente como é mantida no arquivo Bigbed, então analisar a ela é deixado para você. Para determinar quais são os vários campos nessas string, consulte a string 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."
)
Observe que as três primeiras entradas na sequência SQL não fazem parte da string.
Se você só precisar saber onde estão as entradas e não seus valores associados, poderá salvar a memória especificando adicionalmente withString=False nas entries() :
>>> bb.entries('chr1', 10000000, 10020000, withString=False)
[(10009333, 10009640), (10014007, 10014289), (10014373, 10024307)]
Se você abriu um arquivo para gravar, precisará dar um cabeçalho antes de adicionar quaisquer entradas. O cabeçalho contém todos os cromossomos, em ordem , e seus tamanhos. Se o seu genoma tiver dois cromossomos, Chr1 e Chr2, de comprimentos 1 e 1,5 milhão de bases, o seguinte adicionaria um cabeçalho apropriado:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)])
Os cabeçalhos Bigwig são sensíveis ao minúsculas, portanto chr1 e Chr1 são diferentes. Da mesma forma, 1 e chr1 não são os mesmos, então você não pode misturar nomes de cromossomos Ensembl e UCSC. Depois de adicionar um cabeçalho, você pode adicionar entradas.
Por padrão, até 10 "níveis de zoom" são construídos para arquivos bigwig. Você pode alterar esse número padrão com o argumento opcional maxZooms . Um uso comum disso é criar um arquivo bigwig que simplesmente contém intervalos e não níveis de zoom:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)], maxZooms=0)
Se você definir maxTooms=0 , observe que o IGV e muitas outras ferramentas não funcionarão, pois assumem que pelo menos um nível de zoom estará presente. Você é aconselhado a usar o padrão, a menos que não espere que os arquivos bigwig sejam usados por outros pacotes.
Supondo que você tenha aberto um arquivo para gravar e adicionar um cabeçalho, você pode adicionar entradas. Observe que as entradas devem ser adicionadas em ordem, pois os arquivos bigwig sempre contêm intervalos ordenados. Existem três formatos que os arquivos bigwig podem usar internamente para armazenar entradas. O formato mais comumente observado é idêntico a um arquivo de grão de cama:
chr1 0 100 0.0
chr1 100 120 1.0
chr1 125 126 200.0
Essas entradas seriam adicionadas da seguinte forma:
>>> 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 da compressão.
O segundo formato usa uma extensão fixa, mas um tamanho de etapa variável entre as entradas. Estes podem ser representados em um arquivo Wiggle como:
variableStep chrom=chr1 span=20
500 -2.0
600 150.0
635 25.0
As entradas acima descrevem posições (1 baseado em 1) 501-520, 601-620 e 636-655. Estes seriam adicionados da seguinte maneira:
>>> bw.addEntries("chr1", [500, 600, 635], values=[-2.0, 150.0, 25.0], span=20)
Cada entrada desse tipo ocupa 8 bytes antes da compressão.
O formato final usa uma etapa e extensão fixa para cada entrada, correspondendo ao formato FILLSTEP Wiggle:
fixedStep chrom=chr1 step=30 span=20
-5.0
-20.0
25.0
As entradas acima descrevem as bases (1 baseado em 1) 901-920, 931-950 e 961-980 e seriam adicionadas da seguinte forma:
>>> bw.addEntries("chr1", 900, values=[-5.0, -20.0, 25.0], span=20, step=30)
Cada entrada desse tipo ocupa 4 bytes.
Observe que o Pybigwig tentará impedir que você adicione entradas em um pedido incorreto. Isso, no entanto, requer excesso de cabeça adicional. Se isso não for aceitável, você pode simplesmente especificar validate=False ao adicionar entradas:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [100, 0, 125], ends=[120, 5, 126], values=[0.0, 1.0, 200.0], validate=False)
Obviamente, você é responsável por garantir que não adicione entradas fora de ordem. Caso contrário, os arquivos resultantes não seriam utilizáveis.
Um arquivo pode ser fechado com um simples bw.close() , como geralmente é feito com outros tipos de arquivo. Para arquivos abertos para gravar, o fechamento de um arquivo grava qualquer entradas em buffer em disco, constrói e grava o índice de arquivo e constrói os níveis de zoom. Consequentemente, isso pode levar um pouco de tempo.
A partir da versão 0.3.0, o Pybigwig suporta a entrada de coordenadas usando números inteiros e vetores Numpy em algumas funções se Numpy foi instalado antes da instalação do Pybigwig . Para determinar se o Pybigwig foi instalado com suporte Numpy, verificando o acessador numpy :
>>> import pyBigWig
>>> pyBigWig.numpy
1
Se pyBigWig.numpy for 1 , o pybigwig foi compilado com suporte numpy. Isso significa que addEntries() pode aceitar 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()
Além disso, values() podem gerar diretamente um vetor 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'>
Se você não possui o CURL instalado, o Pybigwig será instalado sem a capacidade de acessar arquivos remotos. Você pode determinar se poderá acessar arquivos remotos com pyBigWig.remote . Se isso retornar 1, você poderá acessar arquivos remotos. Se ele retornar 0, então você não pode.
A partir da versão 0.3.5, o Pybigwig é capaz de ler e escrever arquivos bigwig sem entradas. Observe que esses arquivos geralmente não são compatíveis com outros programas, pois não há definição de como um arquivo bigwig sem entradas deve parecer. Para esse arquivo, o acessador intervals() retornará None , a função stats() retornará uma lista de None duração desejada e values() retornarão [] (uma lista vazia). Isso geralmente deve permitir que programas utilizando o pybigwig continuem sem problemas.
Para aqueles que desejam imitar a funcionalidade do pybigwig/libbigwig a esse respeito, observe que ele analisa o número de bases cobertas (conforme relatado no cabeçalho do arquivo) para verificar arquivos "vazios".
Os arquivos Wiggle, Bigwig e Bigbed usam coordenadas meio abertas baseadas em 0, que também são usadas por esta extensão. Portanto, para acessar o valor da primeira base no chr1 , especificaria -se a posição inicial como 0 e a posição final como 1 . Da mesma forma, as bases 100 a 115 teriam um início de 99 e um final de 115 . Isso é simplesmente por uma questão de consistência com o arquivo bigwig subjacente e pode mudar no futuro.
O Pybigwig também está disponível como um pacote no Galaxy. Você pode encontrá -lo no galpão de ferramentas e a IUC está hospedando a definição XML disso no Github.
