pyBigWig
0.3.23
cで記述されたPython拡張機能は、Bigbedファイルへの迅速なアクセスとBigwigファイルへのアクセスと作成。この拡張機能は、ローカルおよびリモートファイルアクセスにLibbigwigを使用します。
この拡張機能は、次のようにGitHubから直接インストールできます。
pip install pyBigWig
またはコンドラと
conda install pybigwig -c conda-forge -c bioconda
フォロー以外の非パイソン要件をインストールする必要があります。
curl-config config）これらのヘッダーとライブラリが必要です。
基本的な使用法は次のとおりです。
>>> import pyBigWig
これは、作業ディレクトリがPybigwigソースコードディレクトリである場合に機能します。
>>> bw = pyBigWig.open("test/test.bw")
ファイルが存在しない場合、エラーメッセージが表示され、返されNoneに注意してください。デフォルトでは、すべてのファイルが読み取りのために開かれ、書き込みではありません。 w ：を含むモードを渡すことでこれを変更できます。
>>> bw = pyBigWig.open("test/output.bw", "w")
執筆用に開かれたファイルは、その間隔や統計を照会することはできず、書き込みのみであることに注意してください。書き込み用のファイルを開くと、次にヘッダーを追加する必要があります（以下のセクションを参照）。
ローカルおよびリモートビッグベッドの読み取りアクセスもサポートされています。
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
Bigbedファイルのモードを指定できますが、無視されます。 pyBigWig.open()によって返されるオブジェクトは、Bigwigファイルを開いているかBigbedファイルを開いているかに関係なく同じです。
BigwigとBigbedファイルは両方を開くことができるため、特定のbigWigFileオブジェクトがBigwigまたはBigbedファイルを指しているかどうかを判断する必要がある場合があります。そのために、 isBigWig()およびisBigBed()関数を使用できます。
>>> bw = pyBigWig.open("test/test.bw")
>>> bw.isBigWig()
True
>>> bw.isBigBed()
False
bigWigFileオブジェクトにはchroms()アクセターでアクセスできる染色体の長さを保持する辞書が含まれています。
>>> bw.chroms()
dict_proxy({'1': 195471971L, '10': 130694993L})
特定の染色体を直接照会することもできます。
>>> bw.chroms("1")
195471971L
長さは「長い」整数タイプに保存されているため、 L接尾辞があります。存在しない染色体を指定すると、出力はありません。
>>> bw.chroms("c")
>>>
Bigwigのヘッダーを印刷するのが便利な場合があります。これは、バージョン（通常4 ）、ズームレベルの数（ nLevels ）、記述されているベースの数（ nBasesCovered ）、最小値（ minVal ）、最大値（ maxVal ）、すべての値（ sumData ）の合計、およびすべての2乗値（ sumSquared ）の合計。これらの最後の2つは、平均と標準偏差を決定するために必要です。
>>> bw.header()
{'maxVal': 2L, 'sumData': 272L, 'minVal': 0L, 'version': 4L, 'sumSquared': 500L, 'nLevels': 1L, 'nBasesCovered': 154L}
これはビッグベッドファイルでも可能であり、同じ辞書キーが存在する可能性があることに注意してください。 maxVal 、 sumData 、 minVal 、 sumSquaredなどのエントリは、ほとんど意味がありません。
Bigwigファイルは、ポジションとそれらの範囲に関連する値を保存するために使用されます。通常、範囲を超えて平均値にすばやくアクセスしたいと考えています。これは非常に簡単です。
>>> bw.stats("1", 0, 3)
[0.2000000054637591]
平均値の代わりに、代わりに最大値が必要だとします。
>>> bw.stats("1", 0, 3, type="max")
[0.30000001192092896]
その他のオプションは、「最小値」（最小値）、「カバレッジ」（カバーされた塩基の割合）、および「STD」（値の標準偏差）です。
代わりに、特定の間隔で均等に間隔を空けたビンの数の値を計算したい場合がよくありますが、これも簡単です。
>>> bw.stats("1",99, 200, type="max", nBins=2)
[1.399999976158142, 1.5]
nBinsデフォルトで1で、 type meanであると同じように。
開始位置と終了位置が省略されている場合、染色体全体が使用されます。
>>> bw.stats("1")
[1.3351851569281683]
レイリーダーへのメモ：このセクションはかなり技術的で、完全性のためにのみ含まれています。要約は、ニーズが必要な場合は正確な平均/max/etcが必要な場合です。間隔または間隔の要約値、および速度の小さなトレードオフが許容できることは、
stats()関数のexact=Trueオプションを使用する必要があることを許容します。
デフォルトでは、Bigwigファイルの範囲に関する統計を計算するには、直感的ではない側面がいくつかあります。 Bigwig形式は、もともとゲノムブラウザーのコンテキストで作成されました。そこでは、特定の間隔の正確な要約統計を計算することは、おおよその統計を迅速に計算できるよりも重要ではありません（結局のところ、ブラウザは多くの連続した間隔をすばやく表示し、スクロール/ズームをサポートする必要があります）。このため、Bigwigファイルには、間隔値の関連性だけでなく、2乗sum of values / sum of squared values / minimum value / maximum value /さまざまなサイズの等しくサイズのビンのnumber of bases coveredも含まれています。これらのさまざまなサイズは、「ズームレベル」と呼ばれます。最小のズームレベルには、ファイルの平均インターバルサイズの16倍のビンがあり、その後の各ズームレベルには、前のビンが4倍大きくなっています。この方法論は、Kentのツールで使用されているため、現在既存のほとんどすべてのBigwigファイルで使用される可能性があります。
概要統計のためにBigwigファイルが照会された場合、間隔のサイズはズームレベルを使用するかどうか、もしそうならどちらを使用するかを決定するために使用されます。最適なズームレベルは、目的の間隔の幅の半分以下の最大のビンを持つものです。そのようなズームレベルが存在しない場合、元の間隔は計算に使用されます。
他のツールとの一貫性のために、Pybigwigはこの同じ方法論を採用しています。ただし、これは（a）直感的ではなく、（b）一部のアプリケーションでは望ましくないため、Pybigwigはインターバルサイズに関係なく正確な要約統計の計算を可能にします（つまり、ズームレベルを無視できます）。これは最初にここで提案され、例は以下にあります。
>>> import pyBigWig
>>> from numpy import mean
>>> bw = pyBigWig.open("http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeCrgMapabilityAlign75mer.bigWig")
>>> bw.stats('chr1', 89294, 91629)
[0.20120902053804418]
>>> mean(bw.values('chr1', 89294, 91629))
0.22213841940688142
>>> bw.stats('chr1', 89294, 91629, exact=True)
[0.22213841940688142]
stats()メソッドは、各ベースの元の値を取得するために使用できますが（たとえば、 nBinsベースの数に設定することにより）、代わりにvalues()アクセサを使用することが望ましいです。
>>> bw.values("1", 0, 3)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896]
生成されたリストには、指定された範囲内のすべてのベースに対して常に1つの値が含まれます。特定のベースにBigwigファイルに関連する値がない場合、返された値はnanになります。
>>> bw.values("1", 0, 4)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896, nan]
ある範囲の重複するすべてのエントリを取得するのが便利な場合があります。これはintervals()関数で実行できます。
>>> bw.intervals("1", 0, 3)
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896))
返されたのは、開始位置、エンドエンド位置、値を含むタプルのリストです。したがって、上記の例は、それぞれ0 1および2で0.1 、および0.3の値を持ってい0.2 。
開始位置と終了位置が省略されている場合、指定された染色体上のすべての間隔が返されます。
>>> bw.intervals("1")
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896), (100, 150, 1.399999976158142), (150, 151, 1.5))
Bigwigファイルとは対照的に、Bigbedファイルはエントリを保持します。これは、関連する文字列の間隔です。 entries()関数を使用してこれらのエントリにアクセスできます。
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
>>> bb.entries('chr1', 10000000, 10020000)
[(10009333, 10009640, '61035t130t-t0.026t0.42t404'), (10014007, 10014289, '61047t136t-t0.029t0.42t404'), (10014373, 10024307, '61048t630t-t5.420t0.00t2672399')]
出力は、エントリタプルのリストです。タプル要素は、各エントリのstart位置とend位置であり、その後に関連するstringが続きます。文字列は、Bigbedファイルに保持されているとまったく同じように返されるため、解析することはあなたに残されます。これらの文字列にさまざまなフィールドが何であるかを判断するには、SQL文字列を参照してください。
>>> bb.SQL()
table RnaElements
"BED6 + 3 scores for RNA Elements data"
(
string chrom; "Reference sequence chromosome or scaffold"
uint chromStart; "Start position in chromosome"
uint chromEnd; "End position in chromosome"
string name; "Name of item"
uint score; "Normalized score from 0-1000"
char[1] strand; "+ or - or . for unknown"
float level; "Expression level such as RPKM or FPKM. Set to -1 for no data."
float signif; "Statistical significance such as IDR. Set to -1 for no data."
uint score2; "Additional measurement/count e.g. number of reads. Set to 0 for no data."
)
SQL文字列の最初の3つのエントリは、文字列の一部ではないことに注意してください。
関連する値ではなく、エントリがどこにあるかを知る必要がある場合、 entries()でwithString=Falseを指定することでメモリを保存できます。
>>> bb.entries('chr1', 10000000, 10020000, withString=False)
[(10009333, 10009640), (10014007, 10014289), (10014373, 10024307)]
執筆用のファイルを開いた場合は、エントリを追加する前にヘッダーを与える必要があります。ヘッダーには、すべての染色体とそのサイズが含まれています。ゲノムには、長さ1および150万の塩基の2つの染色体、Chr1とChr2がある場合、次の場合、適切なヘッダーが追加されます。
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)])
Bigwigヘッダーは大文字と小文字が区別されるため、 chr1とChr1は異なります。同様に、 1とchr1同じではないため、EnsemblとUCSC染色体の名前を混ぜることはできません。ヘッダーを追加した後、エントリを追加できます。
デフォルトでは、Bigwigファイル用に最大10 "ズームレベル"が構築されています。 maxZoomsオプションの引数を使用して、このデフォルト番号を変更できます。これの一般的な使用は、間隔を保持し、ズームレベルを保持しないBigwigファイルを作成することです。
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)], maxZooms=0)
maxTooms=0を設定した場合、IGVおよび他の多くのツールは、少なくとも1つのズームレベルが存在すると仮定するため、機能しないことに注意してください。 Bigwigファイルが他のパッケージで使用されるとは思わない限り、デフォルトを使用することをお勧めします。
書き込み用のファイルを開いてヘッダーを追加したと仮定すると、エントリを追加できます。 Bigwigファイルには常に順序付けられた間隔が含まれているため、エントリを順番に追加する必要があることに注意してください。 Bigwigファイルが内部で使用できる3つの形式があります。エントリを保存します。最も一般的に観察される形式は、ベッドグラフファイルと同一です。
chr1 0 100 0.0
chr1 100 120 1.0
chr1 125 126 200.0
これらのエントリは次のように追加されます。
>>> bw.addEntries(["chr1", "chr1", "chr1"], [0, 100, 125], ends=[5, 120, 126], values=[0.0, 1.0, 200.0])
各エントリは、圧縮の前に12バイトを占有します。
2番目の形式では、固定スパンを使用しますが、エントリ間の可変ステップサイズです。これらは、次のようにwiggleファイルで表現できます。
variableStep chrom=chr1 span=20
500 -2.0
600 150.0
635 25.0
上記のエントリは、（1ベースの）位置501-520、601-620、および636-655を説明しています。これらは次のように追加されます：
>>> bw.addEntries("chr1", [500, 600, 635], values=[-2.0, 150.0, 25.0], span=20)
このタイプの各エントリは、圧縮の前に8バイトを占有します。
最終的な形式では、固定ステップの小刻み形式に対応する各エントリの固定ステップとスパンを使用します。
fixedStep chrom=chr1 step=30 span=20
-5.0
-20.0
25.0
上記のエントリは（1ベース）ベース901-920、931-950、および961-980を説明し、次のように追加されます。
>>> bw.addEntries("chr1", 900, values=[-5.0, -20.0, 25.0], span=20, step=30)
このタイプの各エントリは4バイトを占めています。
Pybigwigは、誤った順序でエントリを追加できないようにしようとすることに注意してください。ただし、これには追加のオーバーヘッドが必要です。それが受け入れられない場合は、エントリを追加するときにvalidate=Falseを指定できます。
>>> bw.addEntries(["chr1", "chr1", "chr1"], [100, 0, 125], ends=[120, 5, 126], values=[0.0, 1.0, 200.0], validate=False)
その後、あなたは明らかに、エントリを順番に追加しないことを保証する責任があります。それ以外の場合、結果のファイルは使用できません。
他のファイルタイプで一般的に行われるように、ファイルは単純なbw.close()で閉じることができます。書き込みのために開かれたファイルの場合、ファイルを閉じると、バッファーされたエントリがディスクに書き込み、ファイルインデックスを構築して書き込み、ズームレベルを作成します。その結果、これには少し時間がかかる場合があります。
バージョン0.3.0の時点で、Pybigwigは、Pybigwigを設置する前にnumpyがインストールされている場合、いくつかの関数のnumpy整数とベクトルを使用して座標の入力をサポートしています。 PybigwigがnumpyアクセサをチェックしてNumpyサポートで設置されたかどうかを判断するには：
>>> import pyBigWig
>>> pyBigWig.numpy
1
pyBigWig.numpyが1場合、PybigwigはNumpyサポートで編集されました。これはaddEntries() numpy座標を受け入れることができることを意味します。
>>> import pyBigWig
>>> import numpy
>>> bw = pyBigWig.open("/tmp/delete.bw", "w")
>>> bw.addHeader([("1", 1000)], maxZooms=0)
>>> chroms = np.array(["1"] * 10)
>>> starts = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90], dtype=np.int64)
>>> ends = np.array([5, 15, 25, 35, 45, 55, 65, 75, 85, 95], dtype=np.int64)
>>> values0 = np.array(np.random.random_sample(10), dtype=np.float64)
>>> bw.addEntries(chroms, starts, ends=ends, values=values0)
>>> bw.close()
さらに、 values() numpyベクターを直接出力できます。
>>> bw = bw.open("/tmp/delete.bw")
>>> bw.values('1', 0, 10, numpy=True)
[ 0.74336642 0.74336642 0.74336642 0.74336642 0.74336642 nan
nan nan nan nan]
>>> type(bw.values('1', 0, 10, numpy=True))
<type 'numpy.ndarray'>
Curlがインストールされていない場合、リモートファイルにアクセスすることなくPybigwigがインストールされます。 pyBigWig.remoteを使用してリモートファイルにアクセスできるかどうかを判断できます。それが1を返す場合、リモートファイルにアクセスできます。 0を返す場合、できません。
バージョン0.3.5の時点で、PybigwigはエントリのないBigwigファイルを読み書きすることができます。このようなファイルは、一般に他のプログラムと互換性がないことに注意してください。これは、エントリのないBigwigファイルがどのように見えるかについての定義がないためです。このようなファイルの場合、 intervals()アクセサはNone返され、 stats()関数は目的の長さのNoneを返し、 values() [] （空のリスト）を返します。これにより、通常、Pybigwigを利用して問題なく継続するプログラムが可能になります。
この点に関して、Pybigwig/Libbigwigの機能を模倣したい人のために、「空の」ファイルをチェックするためにカバーされているベースの数（ファイルヘッダーで報告されているように）を調べることに注意してください。
Wiggle、Bigwig、およびBigbedファイルは、0ベースのハーフオープン座標を使用します。これは、この拡張機能でも使用されます。したがって、 chr1の1 baseの値にアクセスするには、開始位置を0として、最終位置は1として指定します。同様に、100〜115のベースには99の開始と115の終了があります。これは、単に基礎となるBigwigファイルとの一貫性のためであり、将来的に変化する可能性があります。
Pybigwigは、Galaxyのパッケージとしても入手できます。ツールシェッドで見つけることができ、IUCは現在、GitHubでこれのXML定義をホストしています。
