BEDファイルは、ゲノム上の位置情報を表現するのによく使われるファイルフォーマットであり、ゲノムブラウザによる可視化などで使われます。
本記事では、BEDファイルフォーマットの構造とともに、pythonかBEDファイルを操作することができるpybedtoolsの使い方を解説します。
BEDファイルとは
BEDファイルはシンプルなタブ区切りのテキスト形式でありながら、柔軟にゲノム領域情報や注釈情報を扱えるため、多くのバイオインフォマティクスツール（例：BEDTools、pybedtools、UCSC Genome Browserなど）で広く利用されています。
以下は、簡単なBEDファイルの例になります。非常にシンプルですね。
chr1 100 200
chr2 300 400
chr3 500 600BEDファイルは、3列（BED3）から最大12列（BED12）まで定義できますが、BED12形式の各フィールドについて下表にまとめました。
| カラム番号 | カラム名 | 説明 | 例 | 備考 |
|---|---|---|---|---|
| 1 | chrom | 対象となる染色体やスキャフォールドの名前。 | chr1 | 必須フィールド |
| 2 | chromStart | 区間の開始位置（0-based）。開始位置は含まれる。 | 11873 | 必須フィールド |
| 3 | chromEnd | 区間の終了位置（0-based）。終了位置は含まれない。 | 14409 | 必須フィールド |
| 4 | name | 区間の名称（遺伝子名やエクソン番号など）。 | uc001aaa.3 | オプションフィールド。 |
| 5 | score | 数値によるスコア。一般には0～1000が推奨されるが任意の値が使用可能。 | 0 | 可視化時などに濃淡表現で利用。 |
| 6 | strand | ストランド情報。プラス（+）またはマイナス（-）を指定。 | + | オプション。 |
| 7 | thickStart | 太く表示される（例：CDS部分）区間の開始位置。 | 12000 | 可視化用情報。 |
| 8 | thickEnd | 太く表示される区間の終了位置。 | 12500 | 可視化用情報。 |
| 9 | itemRgb | 表示時の色をRGB形式（カンマ区切り）で指定。 | 255,0,0 | UCSC Genome Browserなどで利用。 |
| 10 | blockCount | 区間内のブロック（例：エクソン）の数。 | 3 | 複数ブロックがある場合に記載。 |
| 11 | blockSizes | 各ブロックのサイズをカンマ区切りのリストで指定。 | 100,200,150 | ブロックごとの長さ。 |
| 12 | blockStarts | 各ブロックの開始位置を、区間全体の開始位置からの相対位置で指定。 | 0,500,900 | ブロックのオフセット値。 |
pybedtoolsとは
pybedtoolsは、BEDToolsの機能をPython上で簡単に利用できるライブラリです。BEDTools自体は、ゲノム上の区間（interval）同士の交差、統合、差集合などを効率的に計算するための定番ツールとして広く使われています。
BEDToolsのインストール
pybedtoolsを使用するために、本家のBEDToolsをインストールする必要があります。
Macの場合は、Homebrewを使って以下のコマンドでインストールできます。
brew install bedtoolspybedtoolsのインストール
pybedtoolsはpipコマンドを使用して簡単にインストールすることができます。
pip install pybedtoolsファイルの読み込み
サンプルとなるBEDファイルは、pybedtoolsの中にあるものを、一旦読み込んでファイル出力します。
import pybedtools
# サンプルデータの読み込み
a = pybedtools.example_bedtool('a.bed')
b = pybedtools.example_bedtool('b.bed')
# sample1, sample2として保存
a.saveas("sample1.bed")
b.saveas("sample2.bed")これにより、2つのBEDファイルが保存されます。
sample1
chr1 1 100 feature1 0 +
chr1 100 200 feature2 0 +
chr1 150 500 feature3 0 -
chr1 900 950 feature4 0 +sample2
chr1 155 200 feature5 0 -
chr1 800 901 feature6 0 +用意した、2つのBEDファイルをBedToolオブジェクトとして読み込んでみます。
sample1 = pybedtools.BedTool('sample1.bed')
sample2 = pybedtools.BedTool('sample2.bed')BEDファイルの操作
ここでは、pybedtoolsを利用してBEDファイルの基本的な操作を行う例を紹介します。
交差部分の抽出
異なる実験で得られたピーク情報（BED形式）同士の共通領域を抽出するには、intersect()メソッドを使用します。
s1_and_s2 = sample1.intersect(sample2)
print("---sample1----")
sample1.head()
print("---sample2----")
sample2.head()
print("---sample1 and sample2 intersected----")
s1_and_s2.head()実行すると、下記のように共通領域が出力されます。
---sample1----
chr1 1 100 feature1 0 +
chr1 100 200 feature2 0 +
chr1 150 500 feature3 0 -
chr1 900 950 feature4 0 +
---sample2----
chr1 155 200 feature5 0 -
chr1 800 901 feature6 0 +
---sample1 and sample2 intersected----
chr1 155 200 feature2 0 +
chr1 155 200 feature3 0 -
chr1 900 901 feature4 0 +区間の統合
区間を統合するには、merge()メソッドを使用します。
# 100bp以内に位置する区間を統合（ストランド情報も考慮）
sample1.merge(d=100, s=True).head()このように、区間が統合されます。
chr1 1 200
chr1 150 500
chr1 900 950ソートする
順序を並び替えるには、sort()メソッドを使用します。
デフォルトでは、sortメソッドは染色体名（chrom）と開始位置（start）に基づいて昇順でソートを行います。ただし、BEDToolsのsortコマンドが提供するオプションを引数として渡すことで、ソートの挙動をカスタマイズできます。主なオプションは以下のとおりです：
-
-sizeA: 特徴量のサイズ（end - start）に基づいて昇順でソートします。 -
-sizeD: 特徴量のサイズに基づいて降順でソートします。 -
-chrThenSizeA: 染色体名でソートした後、特徴量のサイズに基づいて昇順でソートします。 -
-chrThenSizeD: 染色体名でソートした後、特徴量のサイズに基づいて降順でソートします。 -
-chrThenScoreA: 染色体名でソートした後、スコア（5列目）に基づいて昇順でソートします。 -
-chrThenScoreD: 染色体名でソートした後、スコアに基づいて降順でソートします。
以下のコードは-sizeDオプションを使用して、特徴量のサイズが大きい順にソートします。
sample1.sort(sizeD=True).head()以下のように、特徴量のサイズによって並び替えられます。
chr1 150 500 feature3 0 -
chr1 100 200 feature2 0 +
chr1 1 100 feature1 0 +
chr1 900 950 feature4 0 +
5. まとめ
本記事では、pybedtoolsを利用したBEDファイル操作の基本と、BEDファイルフォーマットの概要を解説しました。今回は、非常に小さなBEDファイルを操作しましたが、ぜひ実際の分析でBEDファイルを操作するときに使ってみてください。
