CoolBoxを使ったHi-CやChIP-seqなどゲノムトラックのプロット

Hi-C解析やChIP-seq解析の結果を描画する可視化ライブラリ CoolBox の紹介と簡単な使い方。

Hi-Cのコンタクトマップを描画するツールはあんまり多くなくて、とくに他の実験で得られたゲノムトラックと同時に描画しようとするとなかなか苦労する。

探索的解析を兼ねてコンタクトマップを可視化するツールとしては JuiceBox や HiGlass がある。とくにHiGlassはおすすめで、いろいろな入力ファイルに対応してさまざまなプロットをインタラクティブに生成できる。だけど、自由度が大きいぶん使い方も複雑で細かなカスタマイズが難しい。

もっと手軽にプロットして、さくっと論文に出てくるような図を作りたい。という欲求にかなりいい感じで答えてくれる CoolBox というツールが最近出たので、ここでは簡単な使い方と実際のデータの描画例を紹介する。

CoolBoxはコンタクトマップとChIPやRNA-seqなどのゲノムカバレッジのトラックを簡単に描画するPythonライブラリ。

論文は以下。

Xu, W., Zhong, Q., Lin, D. et al. CoolBox: a flexible toolkit for visual analysis of genomics data. BMC Bioinformatics 22, 489 (2021). https://doi.org/10.1186/s12859-021-04408-w

Githubのレポジトリとドキュメントは以下。

インストール

Pythonライブラリなのでcondaやpipで簡単に導入できる。が、あまり環境は汚したくないし、dockerから使ってしまうのがおすすめ。

ready-to-useなdockerイメージを配布してくれているので、docker pullすればすぐに使える。

Linuxコマンドとして使うこともできるが、描画はいろいろと試行錯誤があるのでJupyter上で使ったほうが便利。

docker pull nanguage/coolbox
# 自分の描画したいデータが置いてあるディレクトリを/dataにマウントして、
# jupyterにアクセスするために適当なポートを8888にマッピングする
docker run -v $(pwd):/data -p 8888:8888 nanguage/coolbox:latest jupyter notebook --ip=0.0.0.0 --allow-root

使い方

まずはこのライブラリでどんなプロットができるのか、全部盛りの結果を見てみる。

import coolbox
from coolbox.api import *
coolbox.__version__

'0.3.7'

# ここで使うデータはCoolBoxのレポジトリにあるテストデータ
data_dir = './CoolBox/tests/test_data'
# 描画したいそれぞれのデータのパス
cool_file = f"{data_dir}/cool_chr9_4000000_6000000.mcool"
bed_file = f"{data_dir}/tad_chr9_4000000_6000000.bed"
bam_file = f"{data_dir}/bam_chr9_4000000_6000000.bam"
bedpe_file = f"{data_dir}/bedpe_chr9_4000000_6000000.bedpe"
gtf_file = f"{data_dir}/gtf_chr9_4000000_6000000.gtf"
bigwig_file = f"{data_dir}/bigwig_chr9_4000000_6000000.bw"
bedgraph_file = f"{data_dir}/bedgraph_chr9_4000000_6000000.bg"
bedannotation_file = f"{data_dir}/bed_chr9_4000000_6000000.bed"

frame = XAxis() + \
        Cool(cool_file, cmap="JuiceBoxLike", style='triangular', color_bar='vertical') + \
            TADCoverage(bed_file, border_only=True, alpha=1) + \
            Title("Hi-C with TADs") + \
        Spacer(0.1) + \
        BED(bed_file, border_only=True, alpha=1) + \
            Color("#ce00ce") +\
            Title("TADs (BED)") + \
        BAMCov(bam_file) + \
            MinValue(0.0) + MaxValue(70.0) + \
            Title("BAM Coverage") + \
        Spacer(0.1) + \
        Arcs(bedpe_file, line_width=1.5) + \
            Title("Arcs (BEDPE)") + \
        GTF(gtf_file, length_ratio_thresh=0.005) + \
            TrackHeight(6) + Title("Genes (GTF)") + \
        Spacer(0.1) + \
        BigWig(bigwig_file) + \
            Title("BigWig") + \
        BedGraph(bedgraph_file) + \
            Title("BedGraph") + \
        Spacer(0.1) + \
        BED(bedannotation_file) + \
            Feature(height=10, title="Genes (BED)")

frame.plot("chr9:4000000-6000000")

こんな感じの図が得られる。

使い方はかなり単純で、

1. 描画したいトラックのファイルパスをそれぞれの特別なオブジェクト（BED形式ファイル用、BAM用など）に入れて読み込む

2. 描画したい図の「上から下へ順番に」それらのオブジェクトを足し合わせる

3. 最後に、描画したい染色体領域を指定してプロットする

トラックごとにいちいち領域を指定しなくても最後に指定して一気に全部のトラックの位置を揃えて描画してくれる。

CoolBoxでは "+" 演算子が特別な役割を果たしていて、基本的には描画したいトラックを上から下に積み重ねていく。

だが中には、「直前で指定したトラックに修正を加える」ことを目的とした足し算もある。上の例でいうと、直前のトラックの右側にラベルを表示する Title とか、直前のトラックの描画の高さを変える TrackHeight など。Rのggplot2みたいな感じ。

以下、いくつかのトラックオブジェクトについて見てみる。パラメータの詳細や他のトラックについてはドキュメントを参照。

遺伝子アノテーショントラック

GTF形式ファイルが入力できる。

以下の例では、GTFオブジェクトにファイルのパスを入れて、さらに「座標の数値トラック」（ XAxis ）を積み重ねてから表示している。

gtf_file = f"{data_dir}/gtf_chr9_4000000_6000000.gtf"

frame = GTF(gtf_file) + \
        XAxis()
frame.plot("chr9:4000000-6000000")

また、RefSeqのBED形式で記述された遺伝子アノテーションを入力とすることもできる。

この場合は、一般的なBEDファイルを描画するBEDトラックに入れて描画する。

デフォルトの縦の高さでは描画がちょっとつぶれちゃったので、ここでは TrackHeight を足し合わせることで高さを広げている。

bed_file = f"{data_dir}/bed_chr9_4000000_6000000.bed"

frame = BED(bed_file) + \
            TrackHeight(12) + \
        XAxis()
frame.plot("chr9:4000000-6000000")

Hi-Cトラック

Hi-C解析の結果は、.cool 形式や .hic 形式が入力できる。

HiCExplorer で解析した場合は結果がHDF5形式で出てくるが、付属の hicConvertFormat コマンドで簡単にフォーマットを変換できる。

ここでは複数のコンタクトマップ解像度の情報を含んだ multi-cool (.mcool) 形式のファイルを描画する。

まずはよく見る三角形スタイルの描画。

cool形式ファイルは Cool で読み込んで style で描画スタイルを指定する。

cool_file = f"{data_dir}/cool_chr9_4000000_6000000.mcool"

frame = XAxis() + \
        Cool(cool_file, style='triangular')
frame.plot("chr9:4000000-6000000")

マトリックスのスタイルで描画する場合は以下のようにする。

なお、カラーマップはデフォルトではJuiceBoxのような白→赤の変化だが、matplotlibのカラーマップで好きなものを選択できる。

frame = XAxis() + \
        Cool(cool_file, style='matrix', cmap="YlOrRd")
frame.plot("chr9:4000000-6000000")

また、比較的近い領域のコンタクトの描画としてコンタクトマップの対角線周辺を四角く切り抜いた形式の描画もできる。

frame = XAxis() + \
        Cool(cool_file, style='window', depth_ratio=0.3, cmap="RdYlBu_r")
frame.plot("chr9:4500000-5750000")

論文でコンタクトマップの図を描くときは、検出されたTAD (topologically associating domain) を同時に描画することがある。

BED形式で記述されたTADは、通常のBEDオブジェクトでゲノム領域として描画することもできるが、TADとして描画するための特別なオブジェクトも用意されている。

TADCoverage オブジェクトをHi-Cトラックの直後に足し合わせると、コンタクトマップ上に直接、TADのエリアを表示してくれる。

以下では、コンタクトマップ上にTADを点線で描画すると同時に、同じコンタクトマップについて Inverted を作用させて上下反転したコンタクトマップを描画することで、TAD領域におけるコンタクトの大小を比較しやすくしている。

tad_file = f"{data_dir}/tad_chr9_4000000_6000000.bed"

frame = XAxis() + \
        Cool(cool_file, style='window', depth_ratio=0.3, cmap="Blues") + \
            TADCoverage(tad_file, border_only=True, border_color='red', alpha=0.8) + \
        Spacer(0.1) + \
        Cool(cool_file, style='window', depth_ratio=0.3, cmap="Blues") + \
            Inverted()
frame.plot("chr9:4500000-5750000")

さらに、BEDPE形式で記述されたHi-Cの「ループ」など、特定の2領域の相互作用の描画もコンタクトマップに重ねることができる。

peak_file = f"{data_dir}/peak_chr9_4000000_6000000.bedpe"

frame = XAxis() + \
        Cool(cool_file, style='window', depth_ratio=0.3, cmap="Blues") + \
            HiCPeaksCoverage(peak_file, color="red", line_width=5)
frame.plot("chr9:4500000-5750000")

以上はコンタクトマップの対角線周辺における比較的近接した位置のコンタクトの描画だったが、染色体上で遠く離れた位置の相互作用などが見たい場合もある。

そういった離れた位置にあるそれぞれのlocusを指定して描画する便利なオブジェクトとして、 JointView がある。

以下ではコンタクトマップをロードしたオブジェクトと、上側と左側に遺伝子アノテーションのトラックを指定したサブフレームを指定して JointView に食わせて、ふたつの領域を指定して描画している。

RANGE_1 = 'chr9:4750000-5250000'
RANGE_2 = 'chr9:5250000-5750000'

frame = Cool(cool_file)

sub_frames = {
    "left": GTF(gtf_file) + XAxis(),
    "top": XAxis() + GTF(gtf_file)
}

jv = JointView(frame, **sub_frames, space=0, padding_left=0)
jv.plot(RANGE_1, RANGE_2)

ゲノムカバレッジトラック

RNA-seq, ChIP-seqなど、ゲノムカバレッジ描画のためのオブジェクトとして BigWig がある。名前のとおり、基本的には事前にbigWig形式に変換したファイルを入力とする。

bigwig_file = f"{data_dir}/bigwig_chr9_4000000_6000000.bw"

frame = BigWig(bigwig_file) +\
        XAxis()
frame.plot("chr9:4000000-6000000")

描画の色や、最大値と最小値、トラックの高さなど細かく調整できる。

bigwig_file = f"{data_dir}/bigwig_chr9_4000000_6000000.bw"

frame = BigWig(bigwig_file) +\
            Color('blue') + MinValue(0.0) + MaxValue(1000.0) + \
        BigWig(bigwig_file) +\
            Color('red') + MinValue(0.0) + MaxValue(500.0) + \
        BigWig(bigwig_file) +\
            Color('green') + TrackHeight(6) + \
        XAxis()
frame.plot("chr9:4000000-6000000")

いちいち色を指定しなくても、Pythonのwith文を使うことで、複数のトラックの色を一気に指定することも可能。

with Color("#fd9c6b"):
    frame1 = BigWig(bigwig_file) + \
             BigWig(bigwig_file) + \
             BigWig(bigwig_file)

with Color("#66ccff"):
    frame2 = BigWig(bigwig_file) + \
             BigWig(bigwig_file) + \
             BigWig(bigwig_file)

frame = frame1 + frame2 + XAxis()
frame.plot("chr9:4000000-6000000")

また、これはゲノムカバレッジトラックに限った話ではないが、CoolBoxのどんなトラックも、locusを指定した縦棒（ Vlines ）や、特定のエリアの強調表示（ HighLights ）ができる。

ひとつひとつのトラックに Vlines や HighLights を足し合わせて描画することもできるけど、掛け算演算子を作用させることで全体に一気にオーバーレイすることもできる。

locus = ("chr9", 5280000)

frame *= Vlines(locus, line_width=2)

regions= ["chr9:4600000-5000000", "chr9:5750000-5950000"]
highlights = HighLights(regions, color="green", alpha=0.05)

frame *= highlights
frame.plot("chr9:4000000-6000000")

その他のトラック

その他にもいろいろなトラックに対応したオブジェクトが用意されている。Hi-C関連でよく使うのは、相互作用を記述したBEDPE形式やpairs形式を入力して、弧を描いてくれる Arcs オブジェクト。

bedpe_file = f"{data_dir}/bedpe_chr9_4000000_6000000.bedpe"
pairs_file = f"{data_dir}/pairs_chr9_4000000_6000000.pairs"

frame = Arcs(bedpe_file, linewidth=2, line_style="dotted") + \
            TrackHeight(3) + \
        XAxis() + \
        Arcs(pairs_file, linewidth=3) + \
            Inverted() + TrackHeight(6)
frame.plot("chr9:4000000-6000000")

実際のデータを使った論文の図の再現

ここまでを踏まえて、実際のHi-C解析の論文で掲載されていた図の再現をやってみる。

対象としたのは以下の論文。

Bonev, Boyan, et al. "Multiscale 3D genome rewiring during mouse neural development." Cell 171.3 (2017): 557-572.

解析済みのコンタクトマップは 4D Nucleome Data Portalで配布されているものを、またChIP-seq解析の結果に関してはGEO (GSE96107)からそれぞれダウンロードした。

この研究では、マウスの神経分化過程における3次元染色体構造を比較している。ここではそのうち、ES細胞の結果（ES）と神経前駆細胞の結果（NPC）を描画してみる。

まずはデータ読み込み。4DNはmulti-cool形式でデータを用意してくれているので、そのまま Cool オブジェクトに読み込む。十数GBの結構デカいデータだがかなり高速に読み込んで描画してくれるので快適。

mESC_cool_file = './Bonev_B_et_al/4DNESDXUWBD9/4DNFIC21MG3U.mcool'
NPC_cool_file = './Bonev_B_et_al/4DNESJ9SIAV5/4DNFIKK3QG34.mcool'

mESC_mat = Cool(mESC_cool_file, style='matrix', cmap="YlOrRd")
NPC_mat = Cool(NPC_cool_file, style='matrix', cmap="YlOrRd")

A/Bコンパートメントについては、普通のbigWig形式のゲノムトラックなので BigWig に食わせればいいのだが、Hi-Cで頻繁に使うトラックだけにもっと簡単に描画する特別なオブジェクトが用意されている。

ABCompartment は、bigWigファイルを入力するとプラス側（Aコンパートメント）とマイナス側（Bコンパートメント）でそれぞれ別の色で描画してくれる。

mESC_compartment_file = './Bonev_B_et_al/4DNESDXUWBD9/4DNFI92FVZRB.bw'
NPC_compartment_file = './Bonev_B_et_al/4DNESJ9SIAV5/4DNFIRPD9QXK.bw'

mESC_compartment = ABCompartment(mESC_compartment_file)
NPC_compartment = ABCompartment(NPC_compartment_file)

この二つを重ねて描画してみる。結果は論文におけるFig.1Bに対応している。まずはES細胞の結果。

TEST_RANGE = "chr3:10000-150000000"

frame = XAxis() + \
        mESC_mat + \
        mESC_compartment
frame.plot(TEST_RANGE)

つぎに神経前駆細胞の結果。

TEST_RANGE = "chr3:10000-150000000"

frame = XAxis() + \
        NPC_mat + \
        NPC_compartment
frame.plot(TEST_RANGE)

論文では、ひとつの結果として、細胞分化の過程で「Bコンパートメントどうしの相互作用が次第に強くなっていく」現象が報告されている。

ここでは JointView を使って、遠く離れたBコンパートメントどうしのコンタクトの値がESよりもNPCで強くなっていることを確認してみる。

また、同時にH3K9me3（ヘテロクロマチンに関連した転写抑制のヒストン修飾）のChIP-seq結果も描画してみる。

# マウスESのH3K9me3のChIP-seqデータはENCODEからとってきた
mESC_H3K9me3_file = './Bonev_B_et_al/chip_data/ENCODE_mESC/H3K9me3/ENCFF014CID.bigWig'
NPC_H3K9me3_file = './Bonev_B_et_al/chip_data/GSE96107_NPC_H3K9me3.bw'

mESC_H3K9me3 = BigWig(mESC_H3K9me3_file)
NPC_H3K9me3 = BigWig(NPC_H3K9me3_file)

JointView は以下のように、サブフレームに描画するトラックは好きにスタックできる。ここではA/BコンパートメントとChIP-seqのトラックを積み重ねて、下側と左側に描画する。

この図は論文のFig.1Hに対応する。まずはES細胞。

TEST_RANGE_1 = 'chr6:38000000-48000000'
TEST_RANGE_2 = 'chr6:126000000-136000000'

mESC_mat = Cool(mESC_cool_file, style='matrix', cmap="YlOrRd", color_bar=False)

frame = XAxis()+ \
        mESC_compartment + Title('A/B compartment') + \
        mESC_H3K9me3 + Title('H3K9me3')

sub_frames = {
    "left": frame,
    "bottom": frame
}

jv = JointView(mESC_mat, **sub_frames, space=0, padding_left=0)
jv.plot(TEST_RANGE_1, TEST_RANGE_2)

つぎに同じ領域の神経前駆細胞。コンパートメントのパターンは変化していないが、Bコンパートメントどうしのコンタクトがより強く出ているのがわかる。

TEST_RANGE_1 = 'chr6:38000000-48000000'
TEST_RANGE_2 = 'chr6:126000000-136000000'

NPC_mat = Cool(NPC_cool_file, style='matrix', cmap="YlOrRd", color_bar=False)

frame = XAxis()+ \
        NPC_compartment + Title('A/B compartment') + \
        NPC_H3K9me3 + Title('H3K9me3')

sub_frames = {
    "left": frame,
    "bottom": frame
}

jv = JointView(NPC_mat, **sub_frames, space=0, padding_left=0)
jv.plot(TEST_RANGE_1, TEST_RANGE_2)

最後に、ESからNPCへ分化する過程でZfp608遺伝子周辺で生じる、新たなTAD境界の出現（大きなひとつのTADだった領域がふたつのTADに分離）を見てみる。（論文のFig.2H）

ここではコンタクトマップに加えて、H3K27ac, CTCFのChIP-seq結果、さらに、論文には出ていない要素として Insulation score（コンタクトマップから計算できる量で、TAD境界の推定に利用できる）を描画する。

以下のようにトラックを積み重ねて、さらに特定位置のハイライトを全体に作用させて描画する。たしかにNPCでは、Zfp608のあたりでInsulation scoreが落ち込んでTADのsplitが生じているように見える。

mESC_H3K27ac_file = './Bonev_B_et_al/chip_data/ENCODE_mESC/H3K27ac/ENCFF676VTC.bigWig'
NPC_H3K27ac_file = './Bonev_B_et_al/chip_data/GSE96107_NPC_H3K27ac.bw'

mESC_H3K27ac = BigWig(mESC_H3K27ac_file, color="#0000ff")
NPC_H3K27ac = BigWig(NPC_H3K27ac_file, color="#00ff00")

mESC_CTCF_file = './Bonev_B_et_al/chip_data/GSE96107_ES_CTCF.bw'
NPC_CTCF_file = './Bonev_B_et_al/chip_data/GSE96107_NPC_CTCF.bw'

mESC_CTCF = BigWig(mESC_CTCF_file, color="#0000ff")
NPC_CTCF = BigWig(NPC_CTCF_file, color="#00ff00")

genes = GTF('./Bonev_B_et_al/gencode.vM25.chr_patch_hapl_scaff.annotation.gtf')

TEST_RANGE = 'chr18:53800000-56600000'
Zfp608 = 'chr18:54888048-54992555'

mESC_mat = Cool(mESC_cool_file, style='triangular', cmap='RdYlBu_r')
NPC_mat = Cool(NPC_cool_file, style='triangular', cmap="RdYlBu_r", orientation='inverted')

frame = mESC_mat + \
        genes + \
        XAxis() + \
        InsuScore(mESC_mat, color='#0000ff') + Title('Insulation scores') + \
        mESC_H3K27ac + Title('mESC:H3K27ac') + \
        mESC_CTCF + Title('mESC:CTCF') + \
        NPC_H3K27ac + Title('NPC:H3K27ac') + \
        NPC_CTCF + Title('NPC:CTCF') + \
        InsuScore(NPC_mat, color='#00ff00') + Title('Insulation scores') + \
        XAxis() + \
        NPC_mat

frame = frame * HighLights(Zfp608, color="#ff0000", alpha=0.05)

frame.plot(TEST_RANGE)

こんな感じで、他のTAD検出ツールやDifferential analysisのツールなどを使ってすでに描画したい領域が決まっている場合は、必要な構成要素をスタックしてさくさくとゲノムトラックの描画ができる。

なお、CoolBoxは特定領域を指定した静的な描画だけでなく、Jupyter上のipywidgetsを利用したインタラクティブな可視化もやってくれる。

だけど、正直あんまりレスポンスが速くなくて探索的な目的には使いにくいので、論文の図の生成よりも全体をぐりぐり動かして眺めたいときは、HiGlassなどを使った方がいいと思う。