前言I have a Python, I have a BEDTools.
Ah,pybedtools!
pybedtools 是采用 Python 对BEDTools 软件的封装和扩展。为啥要用pybedtools ，而不是直接使用BEDTools 呢？ 当然是我们想在Python 使用 BEDTools 啦（）
Why pybedtools? 作者也是给了一个例子
找出在基因间隔区的SNPs, 然后获取离SNPs距离<5kb以内的基因名字。
若用pybedtools的话：
import pybedtools
snps = pybedtools.example_bedtool('snps.bed.gz')
genes = pybedtools.example_bedtool('hg19.gff')
intergenic_snps = snps.subtract(genes)
nearby = genes.closest(intergenic_snps, d=True, stream=True)
for gene in nearby:
if int(gene[-1]) < 5000:
print(gene.name)
而若用shell 脚本的话,你可能需要这些代码：
snps=snps.bed.gz
genes=hg19.gff
intergenic_snps=/tmp/intergenic_snps
snp_fields=`zcat $snps | awk '(NR == 2){print NF; exit;}'`
gene_fields=9
distance_field=$(($gene_fields + $snp_fields + 1))
intersectBed -a $snps -b $genes -v > $intergenic_snps
closestBed -a $genes -b $intergenic_snps -d \
awk '($'$distance_field' < 5000){print $9;}' \
perl -ne 'm/[ID|Name|gene_id]=(.*?);/; print "$1\n"'
rm $intergenic_snps
通过比较，若用pybedtools，你只需要熟悉Python 和pybedtools。 而用shell实现相同功能的话，你可能需要Perl, bash, awk 和bedtools。
对于我而言，主要还是因为使用pybedtools，可以让我全程使用Python代码得到和bedtools同样的结果。
使用安装
通过conda
$ conda install --channel conda-forge --channel bioconda pybedtools
通过pip 安装
$ pip install pybedtools -i https://pypi.tuna.tsinghua.edu.cn/simple
注：不能在ｗｉｎｄｏｗ系统下安装
Create a BedTool
使用pybedtools ， 第一步就是导入pybedtools 包和创建BedTool对象。BedTool对象封装了BedTools程序所有可用的程序功能，使它们可以更好的在Python中使用。所以，大多情况，我们用pybedtools ，即在操作BedTool对象。BedTool对象的创建可以使用interval file (BED, GFF, GTF, VCF, SAM, or BAM format)。
由文件创建BedTool对象：
# 其中'test.bed' 是文件路径。
test = pybedtools.BedTool('test.bed')
此外，也可以从字符串里创建：
s = '''
chrX 1 100
chrX 25 800
'''
## 由参数from_string控制
a = pybedtools.BedTool(s, from_string=True)
pybedtools 也提供了一些测试文件，下文将使用这些文件作为演示。
# list_example_files会列出示例文件
>>> pybedtools.list_example_files()
['164.gtf',
...
'a.bed',
'a.bed.gz',
'b.bed',
'bedpe.bed',
'bedpe2.bed',
'c.gff',
'd.gff',
...
'venn.b.bed',
'venn.c.bed',
'x.bam',
'x.bed',
'y.bam']
使用示例文件来创建BedTool对象
## 传入文件名'a.bed'即可
a = pybedtools.example_bedtool('a.bed')
查看前几行数据
>>> a.head()
chr1 1 100 feature1 0 +
chr1 100 200 feature2 0 +
chr1 150 500 feature3 0 -
chr1 900 950 feature4 0 +
Intersections
BEDTools 常用来取交集，来看看在pybedtools怎样操作
a = pybedtools.example_bedtool('a.bed')
b = pybedtools.example_bedtool('b.bed')
a_and_b = a.intersect(b)
>>> a.head()
chr1 1 100 feature1 0 +
chr1 100 200 feature2 0 +
chr1 150 500 feature3 0 -
chr1 900 950 feature4 0 +
>>> b.head()
chr1 155 200 feature5 0 -
chr1 800 901 feature6 0 +
>>> a_and_b.head()
chr1 155 200 feature2 0 +
chr1 155 200 feature3 0 -
chr1 900 901 feature4 0 +
若是查看，是否重叠，和intersectBed 用法一样，添加参数u：
a_with_b = a.intersect(b, u=True)
>>> a_with_b.head()
chr1 100 200 feature2 0 +
chr1 150 500 feature3 0 -
chr1 900 950 feature4 0 +
运行intersect方法后，返回的是一个新的 BedTool 示例，其会将内容暂时保存在一个硬盘临时地址上。临时地址获取：
>>> a_and_b.fn
'/tmp/pybedtools.kfnxvuuz.tmp'
>>> a_with_b.fn
'/tmp/pybedtools.qre024y9.tmp'
Saving the results
BedTool.saveas() 方法可以复制BedTool实例指向的临时文件到一个新的地址。同时可以添加trackline，用于直接上传UCSC Genome Browser,，而不需要再次打开文件添加 trackline。
c = a_with_b.saveas('intersection-of-a-and-b.bed', trackline='track name="a and b"')
print(c.fn)
# intersection-of-a-and-b.bed
BedTool.saveas() 方法返回一个新的BedTool实例指向新的地址。
而BedTool.moveto()用于直接移动文件到一个新的地址，若你不需要添加trackline，文件又很大时，这个方法会很快。
d = a_with_b.moveto('another_location.bed')
print(d.fn)
# 'another_location.bed'
既然已经移动的了，也即老的文件，不存在了,若再次查看其内容会报错.
>>> a_with_b.head()
FileNotFoundError Traceback (most recent call last)
/tmp/ipykernel_2075/544676037.py in <module>
.........
FileNotFoundError: [Errno 2] No such file or directory: '/tmp/pybedtools.4uqthcml.tmp'
Default arguments
BEDTools 使用中，我们会指定-a 和-b参数，而在之前的操作中，我们并没有指定。
当这是因为pybedtools给出了默认设置，进行方法调用的实例作为-a, 传入的另一个实例作为-b，即exons对应-a ("a.bed"), snps对应-b("b.bed")
import pybedtools
exons = pybedtools.example_bedtool('a.bed')
snps = pybedtools.example_bedtool('b.bed')
exons_with_snps = exons.intersect(snps, u=True)
$ intersectBed -a a.bed -b b.bed -u > tmpfile
上面两者是一样的。
不过也可以显示的指定a,b
# these all have identical results
x1 = exons.intersect(snps)
x2 = exons.intersect(a=exons.fn, b=snps.fn)
x3 = exons.intersect(b=snps.fn)
x4 = exons.intersect(snps, a=exons.fn)
x1 == x2 == x3 == x4
# True
Chaining methods together (pipe)
方法的链式调用，类似linux shell下的管道操作
例如：
a 和b查看交集与否，然后合并坐标区间，分开运行是这样
x1 = a.intersect(b, u=True)
x2 = x1.merge()
链式调用可以这样：
x3 = a.intersect(b, u=True).merge()
再来个例子
x4 = a.intersect(b, u=True).saveas('a-with-b.bed').merge().saveas('a-with-b-merged.bed')
Operator overloading
操作符重载， 一个amazing 的方法！
x5 = a.intersect(b, u=True)
x6 = a + b
x5 == x6
# True
x7 = a.intersect(b, v=True)
x8 = a - b
x7 == x8
# True
即 + 操作符表示了 intersectBed 使用 -u 蚕食, - 操作符表示了 intersectBed 使用 -v 参数。
使用了操作符重载还是可以继续链式调用的
x9 = (a + b).merge()
Intervals
在pybedtools中， 以Interval对象来表示BED，GFF，GTF或VCF文件中的一行数据。注上文及下文都是是在Python3版本进行演示（不会还有人用Python2吧）
还是继续创建一个BedTool对象作为例子，
a = pybedtools.example_bedtool('a.bed')
feature = a[0]
features = a[1:3]
BedTool 支持切片操作, 获取单行元素是一个Interval对象
>>> type(feature)
pybedtools.cbedtools.Interval
Common Interval attributes
打印Interval对象，会返回其所代表行数据。
>>> print(feature)
chr1 1 100 feature1 0 +
所有feature（指Interval），不管由什么文件类型创建而来BedTool，都具有chrom, start, stop, name, score, and strand 属性 。其中，start, stop是整数，而其他所有其他（包括score）都是字符串。
不过，我们应该经常会有，只有chrom, start, stop 的数据，我们看看pybedtools如何处理其余缺失属性。
>>> feature2 = pybedtools.BedTool('chrX 500 1000', from_string=True)[0]
>>> print(feature2)
chrX 500 1000
貌似，print(feature2)打印的原始行数据。
>>> feature2.chrom
'chrX'
>>> feature2.start
500
>>> feature2.stop
1000
>>> feature2.name
'.'
>>>feature2.score
'.'
>>>feature2.strand
'.'
缺失默认值是“.”.
Indexing into Interval objects
Interval 可以通过list 和 dict 的方式进行索引。
>>> feature2[:]
['chrX', '500', '1000']
>>> feature2[0]
'chrX'
>>> feature2["chrom"]
'chrX'
不过以字典方式进行索引有个好处，对缺失值可以自动返回默认值，而通过整数索引会报错
>>> feature2["name"]
'.'
>>> feature2[3]
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
/tmp/ipykernel_2075/566460392.py in <module>
----> 1 feature2[3]
/opt/conda/lib/python3.9/site-packages/pybedtools/cbedtools.pyx in pybedtools.cbedtools.Interval.__getitem__()
IndexError: field index out of range
Fields
Interval 有一个 Interval.fields 属性， 将原始行分割成一个list. 使用整数索引时，对这个fields list 进行索引。
>>> f = pybedtools.BedTool('chr1 1 100 asdf 0 + a b c d', from_string=True)[0]
>>> f.fields
['chr1', '1', '100', 'asdf', '0', '+', 'a', 'b', 'c', 'd']
>>> len(f.fields)
10
BED is 0-based, others are 1-based
多格式使用时，一个麻烦的事，BED文件的坐标系统不同于GFF/GTF/VCF文件.
- BED 文件是0-based，（染色体上第一个位置为0），特征不包含stop 位置
- GFF, GTF, and VCF 文件是 1-based （即染色体上第一个位置为1）特征包含stop位置
为了方便，pybedtools 做了这些约定：
Interval.start是0-based start position，不管什么文件。即使是GFF,或者其他1-based feature。- 使用
len()获取Interval的长度，总是返回Interval.stop - Interval.start.这样不管什么文件格式，都能保证长度正确。也简化了潜在代码。我们可以人为所有Intervals是一样的 Interval.fields内容全是字符串，是对原始行的分割。- 所以对于GFF feature， Interval.fields[3] 和 Interval[3] 与 Interval.start 是不一样的。即Interval.fields[3] = Interval.start +1.因为Interval.fields[3]是原始的1-based 数据，而Interval.start 采用0-based
Worked example
下面举两个例子
创建一个GFF Interval
gff = ["chr1",
"fake",
"mRNA",
"51", # <- start is 1 greater than start for the BED feature below
"300",
".",
"+",
".",
"ID=mRNA1;Parent=gene1;"]
gff = pybedtools.create_interval_from_list(gff)
创建一个和之前 gff 坐标一样的BED Interval。但它们的坐标的系统不一样。
bed = ["chr1",
"50",
"300",
"mRNA1",
".",
"+"]
bed = pybedtools.create_interval_from_list(bed)
确认下各自的文件格式，
>>> gff.file_type
'gff'
>>> bed.file_type
'bed'
各自的原始表示
>>> print(gff)
chr1 fake mRNA 51 300 . + . ID=mRNA1;Parent=gene1;
>>> print(bed)
chr1 50 300 mRNA1 . +
>>> bed.start == gff.start == 50
True
>>> bed.end == gff.end == 300
300
>>> len(bed) == len(gff) == 250
True
GFF features have access to attributes
GFF and GTF files have lots of useful information in their attributes field (the last field in each line). 这些attributes 可以通过Interval.attrs访问。
>>> print(gff)
chr1 fake mRNA 51 300 . + . ID=mRNA1;Parent=gene1;
>>> gff.attrs
{'ID': 'mRNA1', 'Parent': 'gene1'}
>>> gff.attrs['Awesomeness'] = "99"
>>> gff.attrs['ID'] = 'transcript1'
可以直接修改attrs
gff.attrs['Awesomeness'] = "99"
gff.attrs['ID'] = 'transcript1'
print(gff.attrs)
# {'ID': 'transcript1', 'Parent': 'gene1', 'Awesomeness': '99'}
print(gff)
# chr1 fake mRNA 51 300 . + . ID=transcript1;Parent=gene1;Awesomeness=99;
Filtering
BedTool.filter() 可以对BedTool 对象进行过滤。传递一个函数，其接收的第一个参数是一个Interval. 返回True/False来进行过滤。
挑选>100bp的特征
a = pybedtools.example_bedtool('a.bed')
b = a.filter(lambda x: len(x) > 100)
print(b)
# chr1 150 500 feature3 0 -
也可以定义的更加通用。挑选长度由传入参数决定
def len_filter(feature, L):
"Returns True if feature is longer than L"
return len(feature) > L
a = pybedtools.example_bedtool('a.bed')
print(a.filter(len_filter, L=200))
# chr1 150 500 feature3 0 -
也是支持链式调用
a.intersect(b).filter(lambda x: len(x) < 100).merge()
Fast filtering functions in Cython
featurefuncs 模块里包含了一些用Cython实现的功能，相比用纯Python,速度大约提升70%左右。
from pybedtools.featurefuncs import greater_than
def L(x,width=100):
return len(x) > 100
### Cython 的长度比较大于实现
test.filter(greater_than, 100)
### 纯Python的大于实现
test.filter(L, 100)
Each
相似BedTool.filter(), 作用函数应用于每个Interval,根据返回布尔值来判断保留与否。BedTool.each()也是将函数应用于每个Interval， 但主要是对Interval进行修改。
a = pybedtools.example_bedtool('a.bed')
b = pybedtools.example_bedtool('b.bed')
## intersect" with c=True, 重复特征的计数
with_counts = a.intersect(b, c=True)
然后定义一个函数，进行计数的归一化，
def normalize_count(feature, scalar=0.001):
"""
assume feature's last field is the count
"""
counts = float(feature[-1])
normalized = round(counts / (len(feature) * scalar), 2)
# need to convert back to string to insert into feature
feature.score = str(normalized)
return feature
>>> with_counts.head()
chr1 1 100 feature1 0 + 0
chr1 100 200 feature2 0 + 1
chr1 150 500 feature3 0 - 1
chr1 900 950 feature4 0 + 1
>>> normalized = with_counts.each(normalize_count)
>>> print(normalized)
chr1 1 100 feature1 0.0 + 0
chr1 100 200 feature2 10.0 + 1
chr1 150 500 feature3 2.86 - 1
chr1 900 950 feature4 20.0 + 1
链式调用：
a.intersect(b, c=True).each(normalize_count).filter(lambda x: float(x[4]) < 1e-5)
