BedTools被广泛应用于基因组数据的区间操作。而pybedtools则扩增bedtools的功能，并将其包装成python可用的特征水平操作。

全部的在线文档见：http://daler.github.io/pybedtools/

为什么选择pybedtools？

以下是一个如何获取基因间SNP在5kb以内的基因名称的例子。

from pybedtools import BedTool

snps = BedTool('snps.bed.gz') # [1]

genes = BedTool('hg19.gff') # [1]

intergenic_snps = snps.subtract(genes) # [2]

nearby = genes.closest(intergenic_snps, d=True, stream=True) # [2, 3]

for gene in nearby: # [4]

    if int(gene[-1]) < 5000: # [4]

    print gene.name # [4]

使用到的特征概括如下：

[1] 支持所有BEDTools支持的格式 (例子中的是gzipped BED和GFF)

[2] 对所有BEDTools的程序和参数(例子中的是subtract、closest 以及将e -d flag参数传递到closest);

[3] 传递结果(例如 Unix pipes, 在本例子中是stream=True)

[4] 通过参数例获取数据的特性声明结果中各种特征 (在此时 [-1] and .name)

以下，是使用shell脚本进行相同的分析。值得注意的是这个是需要perl，bash，awk的知识的，而运行时间与pybedtools相同。

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

.