使用pyBigWig模块查看bigwig文件中的内容

bam, bedgraph, bigwig是3种常见的存储测序深度信息的文件,都可以方便的导入IGV浏览器进行查看,其中bigwig最为常用。在chip_seq, atac_seq中,通常都会提供该种格式的文件,来来可视乎测序深度的分布。

bigwig是一种二进制格式的文件,常规情况下,无法直接浏览其内容。在python中,通过pyBigWig模块,可以方便的查看其文本内容,该模块的基本用法如下

1. 打开文件

该模块支持bigbed和bigwig两种文件格式,打开文件的代码如下

>>> bw = pyBigWig.open('ZM24TRK4.bigwig')
>>> bw.isBigBed()
False
>>> bw.isBigWig()
True
另外,还提供了相应的函数,来帮助我们判断文件格式。

2. 读取内容

测序深度的统计,有固定窗口和变长窗口两种方式,这两种都是针对染色体进行统计,通过如下方式可以查看文件中包含的染色体以及长度

>>> bw.chroms()
{'D10': 64331360L, 'D11': 69233432L, 'D12': 60589447L, 'D13': 60180343L, 'A11': 115914562L, 'A10': 109403145L, 'A13': 107345356L, 'A12': 102340594L, 'D07': 55165807L, 'D06': 63144374L, 'D05': 62317630L, 'D04': 54327443L, 'D03': 53049694L, 'D02': 68262679L, 'D01': 61090159L, 'D09': 50656052L, 'D08': 66058605L, 'A08': 121158745L, 'A09': 79884170L, 'A02': 104266397L, 'A03': 105393268L, 'A01': 112577161L, 'A06': 121900087L, 'A07': 93248268L, 'A04': 80868428L, 'A05': 107448258L}
>>> bw.chroms('D10')
64331360L

再结合intervals方法,就可以遍历所有染色体的测序深度信息了,代码如下

>>> for chrom in bw.chroms():
...     for entry in bw.intervals(chrom):
...             start, end, pos = entry
...             print(chrom, start, end, pos)
...

3. 提取子集

intervals方法有多个参数,必选参数为染色体名称,除此之外,还可以添加起始和终止位置,返回特定区间内的测序深度,用法如下

>>> bw.intervals('D10', 0, 200)
((0, 50, 1.0), (50, 100, 2.0), (100, 150, 3.0), (150, 200, 4.0))

4. 关闭文件

文件读取完后,要记得关闭文件,代码如下

>>> bw.close()

通过该模块,可以将bigwig的内容转换为纯文本,帮助我们更加直观的了解bigwig中存储的信息。

·end·
(0)

相关推荐