跳转到主要内容

使用libBigWig访问bigWig文件的包

项目描述

PyPI version Travis-CI status bioconda-badge DOI

pyBigWig

一个用C语言编写的Python扩展,用于快速访问bigBed文件以及访问和创建bigWig文件。此扩展使用libBigWig进行本地和远程文件访问。

目录

安装

您可以直接使用以下命令从github安装此扩展:

pip install pyBigWig

或者使用conda

conda install pybigwig -c conda-forge -c bioconda

需求

以下非Python依赖项必须安装:

  • libcurl(以及curl-config配置)
  • zlib

需要这些的头文件和库。

用法

基本用法如下:

加载扩展

>>> import pyBigWig

打开bigWig或bigBed文件

如果您的当前工作目录是pyBigWig源代码目录,这将工作。

>>> bw = pyBigWig.open("test/test.bw")

请注意,如果文件不存在,您将看到错误消息,并且将返回None。默认情况下,所有文件都是为读取而打开的,而不是写入。您可以通过传递包含w的模式来更改此设置。

>>> bw = pyBigWig.open("test/output.bw", "w")

请注意,打开用于写入的文件不能查询其区间或统计信息,它只能写入。如果您打开一个文件用于写入,那么您接下来需要添加一个头信息(请参阅下面的章节)。

支持本地和远程bigBed读取访问

>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")

虽然您可以为大Bed文件指定模式,但它将被忽略。pyBigWig.open()返回的对象在打开bigWig或bigBed文件时是相同的。

确定文件类型

由于bigWig和bigBed文件都可以打开,可能需要确定给定的bigWigFile对象指向的是bigWig还是bigBed文件。为此,可以使用isBigWig()isBigBed()函数。

>>> bw = pyBigWig.open("test/test.bw")
>>> bw.isBigWig()
True
>>> bw.isBigBed()
False

访问染色体列表及其长度

bigWigFile对象包含一个字典,其中包含染色体长度,可以通过chroms()访问器访问。

>>> bw.chroms()
dict_proxy({'1': 195471971L, '10': 130694993L})

您还可以直接查询特定的染色体。

>>> bw.chroms("1")
195471971L

长度存储为“长”整数类型,这就是为什么有一个L后缀。如果您指定一个不存在的染色体,则不会输出任何内容。

>>> bw.chroms("c")
>>> 

打印头信息

有时打印bigWig的头信息很有用。这里以Python字典的形式呈现,包含:版本(通常是4)、缩放级别数(nLevels)、描述的基数字符数(nBasesCovered)、最小值(minVal)、最大值(maxVal)、所有值的总和(sumData)以及所有平方值的总和(sumSquared)。后两个值用于确定平均值和标准差。

>>> bw.header()
{'maxVal': 2L, 'sumData': 272L, 'minVal': 0L, 'version': 4L, 'sumSquared': 500L, 'nLevels': 1L, 'nBasesCovered': 154L}

请注意,这也适用于bigBed文件,并将出现相同的字典键。例如,maxValsumDataminValsumSquared等条目基本上没有意义。

计算范围上的摘要信息

bigWig文件用于存储与位置和范围相关的值。通常我们想要快速访问某个范围的平均值,这非常简单。

>>> bw.stats("1", 0, 3)
[0.2000000054637591]

假设我们想要的是最大值,而不是平均值

>>> bw.stats("1", 0, 3, type="max")
[0.30000001192092896]

其他选项是“min”(最小值)、“coverage”(覆盖的基分数)和“std”(值的平均值)。

通常我们想要计算给定区间内均匀分布的若干个bins的值,这也很简单。

>>> bw.stats("1",99, 200, type="max", nBins=2)
[1.399999976158142, 1.5]

nBins默认为1,就像type默认为mean一样。

如果省略起始和结束位置,则使用整个染色体。

>>> bw.stats("1")
[1.3351851569281683]

关于统计和缩放级别的说明

对普通读者的说明:本节相当技术性,仅为完整起见而包含。总结来说,如果您的需求需要区间或区间的精确平均/最大等汇总值,并且可以接受在速度上的微小妥协,那么您应在stats()函数中使用exact=True选项。

默认情况下,在bigWig文件中对范围进行统计计算有一些不直观的方面。bigWig格式最初是在基因组浏览器的背景下创建的。在那里,对给定区间的精确汇总统计不如快速计算近似统计重要(毕竟,浏览器需要能够快速显示多个连续区间并支持滚动/缩放)。因此,bigWig文件不仅包含区间-值关联,还包含各种大小等大小箱子的值的总和/平方值的总和/最小值/最大值/覆盖的碱基数。这些不同的大小被称为“缩放级别”。最小的缩放级别具有16倍文件中平均区间大小的箱子,每个后续缩放级别的箱子比前一个箱子大4倍。这种方法在Kent的工具中使用,因此,很可能在几乎每个现有的bigWig文件中都使用。

当查询bigWig文件以获取汇总统计时,使用区间的大小来确定是否使用缩放级别,如果是,则确定使用哪个级别。最佳缩放级别是具有最大箱子的宽度不超过所需区间宽度一半的那个级别。如果不存在这样的缩放级别,则改用原始区间进行计算。

为了与其他工具保持一致,pyBigWig采用同样的方法。然而,由于这种方法(A)不直观,并且(B)在某些应用中不受欢迎,pyBigWig允许无论区间大小如何都进行精确的汇总统计计算(即,它允许忽略缩放级别)。这最初在这里提出,下面是一个示例

>>> 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]

检索范围中单个碱基的值

虽然可以使用stats()方法检索每个基对的原始值(例如,通过将nBins设置为碱基数),但最好使用values()访问器。

>>> bw.values("1", 0, 3)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896]

产生的列表将始终包含指定范围内每个基对的一个值。如果某个基对在大Wig文件中没有关联的值,则返回的值将是nan

>>> bw.values("1", 0, 4)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896, nan]

检索范围中的所有区间

有时检索与某些范围重叠的所有条目很方便。这可以使用intervals()函数完成

>>> bw.intervals("1", 0, 3)
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896))

返回的是包含:起始位置、结束位置和值的元组的列表。因此,上面的例子在位置012的值分别为0.10.20.3

如果省略起始和结束位置,则返回指定染色体上的所有区间

>>> bw.intervals("1")
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896), (100, 150, 1.399999976158142), (150, 151, 1.5))

检索bigBed条目

与bigWig文件不同,bigBed文件包含条目,它们是与字符串关联的区间。您可以使用entries()函数访问这些条目

>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
>>> bb.entries('chr1', 10000000, 10020000)
[(10009333, 10009640, '61035\t130\t-\t0.026\t0.42\t404'), (10014007, 10014289, '61047\t136\t-\t0.029\t0.42\t404'), (10014373, 10024307, '61048\t630\t-\t5.420\t0.00\t2672399')]

输出是一个条目元组的列表。元组元素是每个条目的startend位置,后面跟着其关联的string。字符串以它在bigBed文件中保持的方式返回,因此解析它留给您。要确定这些字符串中的各种字段,请查阅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."
    )

请注意,SQL字符串中的前三个条目不属于字符串。

如果您只需要知道条目的位置而不是它们的关联值,则可以在entries()中指定withString=False以节省内存

>>> bb.entries('chr1', 10000000, 10020000, withString=False)
[(10009333, 10009640), (10014007, 10014289), (10014373, 10024307)]

向bigWig文件添加头信息

如果您已打开文件进行写入,则您在添加任何条目之前需要提供标题。标题包含所有染色体,按顺序以及它们的尺寸。如果您有长度为1和1.5百万碱基对的两个染色体chr1和chr2,则以下将添加适当的标题

>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)])

bigWig标题区分大小写,因此chr1Chr1不同。同样,1chr1也不同,因此您不能混合Ensembl和UCSC染色体名称。添加标题后,您就可以添加条目了。

默认情况下,为bigWig文件构建了多达10个“缩放级别”。您可以使用maxZooms可选参数更改此默认数量。这种用途的常见方法是创建一个只包含区间而没有缩放级别的bigWig文件。

>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)], maxZooms=0)

如果您设置maxTooms=0,请注意,IGV和许多其他工具将无法正常工作,因为它们假设至少存在一个缩放级别。除非您不期望bigWig文件被其他包使用,否则建议使用默认设置。

向bigWig文件添加条目

假设您已经打开了用于写入的文件并添加了标题,然后您可以添加条目。请注意,条目必须按顺序添加,因为bigWig文件始终包含有序区间。bigWig文件可以内部使用三种格式来存储条目。最常见观察到的格式与bedGraph文件相同。

chr1	0	100	0.0
chr1	100	120	1.0
chr1	125	126	200.0

这些条目将按照以下方式添加

>>> bw.addEntries(["chr1", "chr1", "chr1"], [0, 100, 125], ends=[5, 120, 126], values=[0.0, 1.0, 200.0])

每个条目在压缩前占用12字节。

第二种格式使用固定的跨度,但条目之间的步长是可变的。这些可以用wiggle文件表示为

variableStep chrom=chr1 span=20
500	-2.0
600	150.0
635	25.0

上述条目描述的是(基于1)位置501-520、601-620和636-655。它们将按照以下方式添加

>>> bw.addEntries("chr1", [500, 600, 635], values=[-2.0, 150.0, 25.0], span=20)

这种类型的每个条目在压缩前占用8字节。

最后一种格式为每个条目使用固定的步长和跨度,对应于固定步长的wiggle格式

fixedStep chrom=chr1 step=30 span=20
-5.0
-20.0
25.0

上述条目描述的是(基于1)碱基901-920、931-950和961-980,将按照以下方式添加

>>> bw.addEntries("chr1", 900, values=[-5.0, -20.0, 25.0], span=20, step=30)

这种类型的每个条目在压缩前占用4字节。

请注意,pyBigWig会尝试防止您添加顺序错误的条目。但这需要额外的开销。如果不能接受,可以在添加条目时指定validate=False

>>> bw.addEntries(["chr1", "chr1", "chr1"], [100, 0, 125], ends=[120, 5, 126], values=[0.0, 1.0, 200.0], validate=False)

因此,您显然需要负责确保不添加顺序错误的条目。否则,生成的文件将很大程度上不可用。

关闭bigWig或bigBed文件

文件可以通过简单的bw.close()关闭,就像通常处理其他文件类型一样。对于打开用于写入的文件,关闭文件会将任何缓冲条目写入磁盘,构建并写入文件索引,并构建缩放级别。因此,这可能需要一些时间。

Numpy

截至版本0.3.0,pyBigWig支持在安装pyBigWig之前安装numpy的情况下使用numpy整数和向量的某些函数输入坐标。要确定pyBigWig是否通过检查numpy访问器安装了numpy支持

>>> import pyBigWig
>>> pyBigWig.numpy
1

如果pyBigWig.numpy1,则pyBigWig是用numpy支持编译的。这意味着addEntries()可以接受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()

此外,values()可以直接输出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'>

远程文件访问

如果您没有安装curl,pyBigWig将无法访问远程文件。您可以使用pyBigWig.remote来确定您是否能够访问远程文件。如果它返回1,则可以访问远程文件。如果它返回0,则不能。

空文件

截至版本0.3.5,pyBigWig能够读取和写入缺少条目的bigWig文件。请注意,这类文件通常与其他程序不兼容,因为没有定义没有条目的bigWig文件应该是什么样子。对于这类文件,intervals()访问器将返回Nonestats()函数将返回一个长度为所需长度的None列表,而values()将返回一个空列表(空列表)。这通常允许使用pyBigWig的程序继续运行而不会出现任何问题。

对于那些希望在这方面模仿pyBigWig/libBigWig的功能的用户,请注意,它通过查看文件标题中报告的覆盖碱基数来检查“空”文件。

关于坐标的说明

Wiggle、bigWig 和 bigBed 文件使用基于0的半开坐标,这个扩展也使用这些坐标。因此,要访问 chr1 上的第一个碱基的值,需要指定起始位置为 0,结束位置为 1。同样,碱基100到115的起始位置为 99,结束位置为 115。这样做是为了与底层 bigWig 文件保持一致性,将来可能会发生变化。

星系

pyBigWig 也可在 Galaxy 中作为软件包使用。您可以在工具库中找到它,并且 IUC 目前在该 GitHub 上托管了此软件包的 XML 定义。

项目详情


下载文件

下载适合您平台的文件。如果您不确定选择哪个,请了解有关 安装软件包 的更多信息。

源分布

pybigwig-0.3.23.tar.gz (75.0 kB 查看哈希值)

上传时间

构建分布

pyBigWig-0.3.23-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (188.3 kB 查看哈希值)

上传时间 CPython 3.12 manylinux: glibc 2.17+ x86-64

pyBigWig-0.3.23-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (187.0 kB 查看哈希值)

上传时间 CPython 3.11 manylinux: glibc 2.17+ x86-64

pyBigWig-0.3.23-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (184.6 kB 查看哈希值)

上传时间 CPython 3.10 manylinux: glibc 2.17+ x86-64

pyBigWig-0.3.23-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (182.3 kB 查看哈希值)

上传时间 CPython 3.9 manylinux: glibc 2.17+ x86-64

由以下提供支持