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)
參考
