pyBigWig
0.3.23
Une extension Python, écrite en C, pour un accès rapide aux fichiers Bigbed et l'accès et la création de fichiers BigWig. Cette extension utilise LibbigWig pour un accès aux fichiers local et distant.
Vous pouvez installer cette extension directement à partir de GitHub avec:
pip install pyBigWig
ou avec conda
conda install pybigwig -c conda-forge -c bioconda
Les exigences non-python suivantes doivent être installées:
curl-config )Les en-têtes et les bibliothèques pour ceux-ci sont nécessaires.
L'utilisation de base est la suivante:
>>> import pyBigWig
Cela fonctionnera si votre répertoire de travail est le répertoire de code source PybigWig.
>>> bw = pyBigWig.open("test/test.bw")
Notez que si le fichier n'existe pas, vous verrez un message d'erreur et None ne sera retourné. Soyez par défaut, tous les fichiers sont ouverts pour la lecture et non l'écriture. Vous pouvez modifier cela en passant un mode contenant w :
>>> bw = pyBigWig.open("test/output.bw", "w")
Notez qu'un fichier ouvert pour l'écriture ne peut pas être interrogé pour ses intervalles ou statistiques, il ne peut être écrit que. Si vous ouvrez un fichier pour l'écriture, vous devrez ensuite ajouter un en-tête (voir la section à ce sujet ci-dessous).
L'accès à lecture locale et à distance en lecture est également pris en charge:
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
Bien que vous puissiez spécifier un mode pour les fichiers Bigbed, il est ignoré. L'objet renvoyé par pyBigWig.open() est le même que vous ouvriez un fichier bigwig ou bigbed.
Étant donné que les fichiers BigWig et Bigbed peuvent tous deux être ouverts, il peut être nécessaire de déterminer si un objet bigWigFile donné pointe vers un fichier bigwig ou bigbed. À cette fin, on peut utiliser les fonctions isBigWig() et isBigBed() :
>>> bw = pyBigWig.open("test/test.bw")
>>> bw.isBigWig()
True
>>> bw.isBigBed()
False
Les objets bigWigFile contiennent un dictionnaire contenant les longueurs chromosomiques, qui peuvent être accessibles avec l'accessoire chroms() .
>>> bw.chroms()
dict_proxy({'1': 195471971L, '10': 130694993L})
Vous pouvez également interroger directement un chromosome particulier.
>>> bw.chroms("1")
195471971L
Les longueurs sont stockées un type entier "long", c'est pourquoi il y a un suffixe L Si vous spécifiez un chromosome inexistant, rien n'est sorti.
>>> bw.chroms("c")
>>>
Il est parfois utile d'imprimer une en-tête BigWig. Ceci est présenté ici comme un dictionnaire Python contenant: la version (généralement 4 ), le nombre de niveaux de zoom ( nLevels ), le nombre de bases décrites ( nBasesCovered ), la valeur minimale ( minVal ), la valeur maximale ( maxVal ), la somme de toutes les valeurs ( sumData ), et la somme de toutes les valeurs carrés ( sumSquared ). Les deux derniers d'entre eux sont nécessaires pour déterminer la moyenne et l'écart type.
>>> bw.header()
{'maxVal': 2L, 'sumData': 272L, 'minVal': 0L, 'version': 4L, 'sumSquared': 500L, 'nLevels': 1L, 'nBasesCovered': 154L}
Notez que cela est également possible pour les fichiers Bigbed et que les mêmes clés de dictionnaire seront présentes. Les entrées telles que maxVal , sumData , minVal et sumSquared ne sont alors pas largement significatives.
Les fichiers BigWig sont utilisés pour stocker les valeurs associées aux positions et à la gamme d'entre elles. En règle générale, nous voulons accéder rapidement à la valeur moyenne sur une plage, ce qui est très simple:
>>> bw.stats("1", 0, 3)
[0.2000000054637591]
Supposons au lieu de la valeur moyenne, nous voulions plutôt la valeur maximale:
>>> bw.stats("1", 0, 3, type="max")
[0.30000001192092896]
Les autres options sont "min" (la valeur minimale), la "couverture" (la fraction des bases couvertes) et "STD" (l'écart type des valeurs).
C'est souvent le cas que nous aimerions plutôt calculer les valeurs d'un certain nombre de bacs uniformément espacés dans un intervalle donné, ce qui est également simple:
>>> bw.stats("1",99, 200, type="max", nBins=2)
[1.399999976158142, 1.5]
nBins est par défaut à 1, tout comme type par défaut mean .
Si les positions de début et de fin sont omises, tout le chromosome est utilisé:
>>> bw.stats("1")
[1.3351851569281683]
Une note au lecteur laïc: cette section est plutôt technique et incluse uniquement pour l'exhaustivité. Le résumé est que si vos besoins nécessitent une moyenne / max / etc. exacte. Valeurs de résumé pour un intervalle ou des intervalles et qu'un petit compromis de vitesse est acceptable, que vous devez utiliser l'option
exact=Truedans la fonctionstats().
Par défaut, il existe des aspects non intuitifs pour calculer les statistiques sur les plages dans un fichier bigwig. Le format BigWig a été créé à l'origine dans le contexte des navigateurs du génome. Là, le calcul des statistiques de résumé exact pour un intervalle donné est moins importante que de pouvoir calculer rapidement une statistique approximative (après tout, les navigateurs doivent être en mesure d'afficher rapidement un certain nombre d'intervalles contigus et de prendre en charge le défilement / zoom). Pour cette raison, les fichiers BigWig contiennent non seulement des associations d'intervalle, mais aussi sum of values / sum of squared values / minimum value / maximum value / number of bases covered pour les bacs de taille égale de différentes tailles. Ces différentes tailles sont appelées «niveaux de zoom». Le plus petit niveau de zoom a des bacs qui sont 16 fois la taille moyenne de l'intervalle dans le fichier et chaque niveau de zoom ultérieur a des bacs 4 fois plus grands que les précédents. Cette méthodologie est utilisée dans les outils de Kent et, par conséquent, probablement utilisée dans presque tous les fichiers BigWig existants.
Lorsqu'un fichier bigwig est interrogé pour une statistique de résumé, la taille de l'intervalle est utilisée pour déterminer s'il faut utiliser un niveau de zoom et, dans l'affirmative, lequel. Le niveau de zoom optimal est celui qui a les plus gros bacs, pas plus de la moitié de la largeur de l'intervalle souhaité. Si aucun niveau de zoom n'existe, les intervalles d'origine sont plutôt utilisés pour le calcul.
Par souci de cohérence avec d'autres outils, Pybigwig adopte cette même méthodologie. Cependant, comme il s'agit (a) non intuitif et (b) indésirable dans certaines applications, Pybigwig permet de calculer les statistiques de résumé exactes quelle que soit la taille de l'intervalle (c'est-à-dire qu'elle permet d'ignorer les niveaux de zoom). Cela a été initialement proposé ici et un exemple est ci-dessous:
>>> 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]
Bien que la méthode stats() puisse être utilisée pour récupérer les valeurs d'origine pour chaque base (par exemple, en définissant nBins sur le nombre de bases), il est préférable d'utiliser l'accessoire values() .
>>> bw.values("1", 0, 3)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896]
La liste produite contiendra toujours une valeur pour chaque base de la plage spécifiée. Si une base particulière n'a pas de valeur associée dans le fichier bigwig, la valeur renvoyée sera nan .
>>> bw.values("1", 0, 4)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896, nan]
Parfois, il est pratique de récupérer toutes les entrées qui chevauchent une plage. Cela peut être fait avec la fonction intervals() :
>>> bw.intervals("1", 0, 3)
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896))
Ce qui est renvoyé est une liste de tuples contenant: la position de début, la position de fin et la valeur. Ainsi, l'exemple ci-dessus a des valeurs de 0.1 , 0.2 et 0.3 aux positions 0 , 1 et 2 , respectivement.
Si la position de début et de fin est omise, tous les intervalles sur le chromosome spécifié sont retournés:
>>> bw.intervals("1")
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896), (100, 150, 1.399999976158142), (150, 151, 1.5))
Contrairement aux fichiers BigWig, les fichiers Bigbed contiennent des entrées, qui sont des intervalles avec une chaîne associée. Vous pouvez accéder à ces entrées à l'aide de la fonction 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 sortie est une liste de tuples d'entrée. Les éléments de tuple sont la position start et end de chaque entrée, suivie de sa string associée. La chaîne est retournée exactement car elle est maintenue dans le fichier Bigbed, donc l'analyse vous est laissée. Pour déterminer quels sont les différents champs dans ces chaînes, consultez la chaîne 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."
)
Notez que les trois premières entrées de la chaîne SQL ne font pas partie de la chaîne.
Si vous avez seulement besoin de savoir où se trouvent les entrées et non leurs valeurs associées, vous pouvez enregistrer la mémoire en spécifiant en outre withString=False dans entries() :
>>> bb.entries('chr1', 10000000, 10020000, withString=False)
[(10009333, 10009640), (10014007, 10014289), (10014373, 10024307)]
Si vous avez ouvert un fichier pour l'écriture, vous devrez lui donner un en-tête avant de pouvoir ajouter des entrées. L'en-tête contient tous les chromosomes, dans l'ordre , et leurs tailles. Si votre génome a deux chromosomes, CHR1 et CHR2, de longueurs de 1 et 1,5 million de bases, ce qui suit ajouterait un en-tête approprié:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)])
Les en-têtes BigWig sont sensibles à la casse, donc chr1 et Chr1 sont différents. De même, 1 et chr1 ne sont pas les mêmes, vous ne pouvez donc pas mélanger les noms de chromosomes Ensembl et UCSC. Après avoir ajouté un en-tête, vous pouvez ensuite ajouter des entrées.
Par défaut, jusqu'à 10 "niveaux de zoom" sont construits pour les fichiers BigWig. Vous pouvez modifier ce numéro par défaut avec l'argument facultatif maxZooms . Une utilisation courante de ceci est de créer un fichier bigwig qui contient simplement des intervalles et pas de niveaux de zoom:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)], maxZooms=0)
Si vous définissez maxTooms=0 , veuillez noter que IGV et de nombreux autres outils ne fonctionneront pas car ils supposent qu'au moins un niveau de zoom sera présent. Il est conseillé d'utiliser la valeur par défaut, sauf si vous ne vous attendez pas à ce que les fichiers BigWig soient utilisés par d'autres packages.
En supposant que vous ayez ouvert un fichier pour l'écriture et ajouté un en-tête, vous pouvez ensuite ajouter des entrées. Notez que les entrées doivent être ajoutées dans l'ordre, car les fichiers BigWig contiennent toujours des intervalles commandés. Il existe trois formats que les fichiers BigWig peuvent utiliser en interne pour stocker des entrées. Le format le plus couramment observé est identique à un fichier Bedgraph:
chr1 0 100 0.0
chr1 100 120 1.0
chr1 125 126 200.0
Ces entrées seraient ajoutées comme suit:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [0, 100, 125], ends=[5, 120, 126], values=[0.0, 1.0, 200.0])
Chaque entrée occupe 12 octets avant la compression.
Le deuxième format utilise une portée fixe, mais une taille de pas variable entre les entrées. Ceux-ci peuvent être représentés dans un fichier Wiggle comme:
variableStep chrom=chr1 span=20
500 -2.0
600 150.0
635 25.0
Les entrées ci-dessus décrivent les positions (basées sur 1) 501-520, 601-620 et 636-655. Ceux-ci seraient ajoutés comme suit:
>>> bw.addEntries("chr1", [500, 600, 635], values=[-2.0, 150.0, 25.0], span=20)
Chaque entrée de ce type occupe 8 octets avant la compression.
Le format final utilise une étape et une portée fixe pour chaque entrée, correspondant au format Wiggle fixe-step:
fixedStep chrom=chr1 step=30 span=20
-5.0
-20.0
25.0
Les entrées ci-dessus décrivent les bases (à 1 base) 901-920, 931-950 et 961-980 et seraient ajoutées comme suit:
>>> bw.addEntries("chr1", 900, values=[-5.0, -20.0, 25.0], span=20, step=30)
Chaque entrée de ce type occupe 4 octets.
Notez que PybigWig essaiera de vous empêcher d'ajouter des entrées dans un ordre incorrect. Cela nécessite cependant une sur-tête supplémentaire. Si cela ne soit pas acceptable, vous pouvez simplement spécifier validate=False lors de l'ajout d'entrées:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [100, 0, 125], ends=[120, 5, 126], values=[0.0, 1.0, 200.0], validate=False)
Vous êtes évidemment responsable de vous assurer que vous n'ajoutez pas les entrées en panne. Les fichiers résultants ne seraient autrement pas que Largyy ne serait pas utilisable.
Un fichier peut être fermé avec un simple bw.close() , comme cela est généralement fait avec d'autres types de fichiers. Pour les fichiers ouverts pour l'écriture, la fermeture d'un fichier écrit toutes les entrées tamponnées sur le disque, construit et écrit l'index de fichier et construit des niveaux de zoom. Par conséquent, cela peut prendre un peu de temps.
À partir de la version 0.3.0, PybigWig prend en charge l'entrée des coordonnées à l'aide de numpy entiers et vecteurs dans certaines fonctions si Numpy a été installé avant d'installer PybigWig . Pour déterminer si Pybigwig a été installé avec un support Numpy en vérifiant l'accessoire numpy :
>>> import pyBigWig
>>> pyBigWig.numpy
1
Si pyBigWig.numpy est 1 , alors Pybigwig a été compilé avec un support Numpy. Cela signifie que addEntries() peuvent accepter les coordonnées 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()
De plus, values() peuvent produire directement un vecteur 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 vous n'avez pas installé Curl, PybigWig sera installé sans la possibilité d'accéder à des fichiers distants. Vous pouvez déterminer si vous pourrez accéder à des fichiers distants avec pyBigWig.remote . Si cela renvoie 1, vous pouvez accéder aux fichiers distants. S'il revient 0, vous ne pouvez pas.
À partir de la version 0.3.5, Pybigwig est capable de lire et d'écrire des fichiers BigWig manquant d'entrées. Veuillez noter que ces fichiers ne sont généralement pas compatibles avec d'autres programmes, car il n'y a aucune définition de la façon dont un fichier BigWig sans entrées devrait être consulté. Pour un tel fichier, l'accessor intervals() ne renverra None , la fonction stats() renverra une liste de None longueur souhaitée, et values() renvoie [] (une liste vide). Cela devrait généralement permettre aux programmes utilisant Pybigwig de se poursuivre sans problème.
Pour ceux qui souhaitent imiter la fonctionnalité de Pybigwig / Libbigwig à cet égard, veuillez noter qu'il examine le nombre de bases couvertes (comme indiqué dans l'en-tête de fichier) pour vérifier les fichiers "vides".
Les fichiers Wiggle, Bigwig et Bigbed utilisent des coordonnées à demi-ouvertures basées sur 0, qui sont également utilisées par cette extension. Ainsi, pour accéder à la valeur de la première base sur chr1 , on spécifierait la position de départ comme 0 et la position finale comme 1 . De même, les bases de 100 à 115 auraient un début de 99 et une fin de 115 . Ceci est simplement pour la cohérence avec le fichier bigwig sous-jacent et peut changer à l'avenir.
Pybigwig est également disponible en tant que package en galaxie. Vous pouvez le trouver dans l'outil et l'IUC héberge actuellement la définition XML de ceci sur GitHub.
