平成29年度
NGSハンズオン講習会
Hi-C解析
2017年9月1日
国立遺伝学研究所 東光一

1

本講義の内容n Hi-C解析とは
l Chromosome Conformation Capture の原理
l Hi-Cで何がわかるか？コンタクトマップの見方
n Hi-C解析の流れ（実習と並行）
l Hi-C解析のツール
l マッピング
l フィルタリング
l 正規化
l ピーク検出、TAD検出など
l 3D構造モデリング2

実習の内容n 目的：
Fastqファイル（公共データ）からスタートして、Hi-C解析論文でよく見るコンタクトマップや、
染色体3次元構造の構築までやってみる
n 複雑なコマンドは打ちません。
必要な操作はすべてpythonスクリプトにまとめてあるので、python 〇〇.py と打つだけ。
「ツールの使い方」よりも、Hi-Cデータの特徴や、データ解析で気をつけなければいけない点を理解する。

3

1_mapping_read_to_genome
2_filtering_reads
Bio-Linux-8.0.7_hm_kh.ova を起動。 3_normalization
4_convert_Juice
すべて、~/HiC の中で実行する。
5_detect_TADs
6_modeling_3D
解析のステップごとに、実行する
pythonスクリプトが入ったディレク
トリ（実行結果はそれぞれのResultsの中）data解析に使うfastqファイル
Index
ヒトゲノムのBowtie2インデックス
refヒトゲノム配列（fasta）
src
使用したライブラリのソースコード

実習の内容

4

Hi-C解析とはシーケンシングによって、ゲノム（染色体）の３次元空間内の立体構造を明らかにする手法。
構造は機能と密接に関わっている。
（ex. エンハンサー・プロモーターループ）
ゲノム配列そのものからはうかがい知れない染色体の立体構造を、ゲノム配列のシーケンシングによって明らかにするのがHi-C解析。

5

NGSの応用
n配列そのものを知りたい
l ゲノム
l メタゲノムl Reseqn

読み取られた配列データのパターンから、何か別のことを知りたい
l RNA-seql ChIP-seq
l ATAC-seql Hi-Cl iRep
精度の高いリファレンスゲノムがあることが前提。

6

Chromosome Conformation Capture (3C)

Dekker, Job, et al. "Capturing chromosome conformation." Science 295.5558 (2002): 1306-1311.

その後、次々に登場した3C-based methodすべての基礎となる手法。
鍵となるのは、DNAの空間的な近接性は、制限消化された断片間でのライゲーションの生じやすさで測れる、という考え方。

7

3C-based technologies

de Wit, Elzo, and Wouter de Laat. "A decade of 3C technologies: insights
into nuclear organization." Genes & development 26.1 (2012): 11-24.

8

4C: Chromosome conformation capture-on-chip

ゲノム上の特定の領域（viewpointと呼ばれる）と相互作用しているのはどの領域か？インバースPCRを利用してアレイで検出。
NGSと組み合わせた 4C-seq という拡張もある。

de Wit, Elzo, and Wouter de Laat. "A decade of 3C technologies: insights
into nuclear organization." Genes & development 26.1 (2012): 11-24.
9

4C: Chromosome conformation capture-on-chip
以下は4C-seqデータを使ってHi-Cデータのvalidationをした例

Ke, Yuwen, et al.
"3D chromatin structures of mature
gametes and structural reprogramming
during mammalian embryogenesis."
Cell 170.2 (2017): 367-381.
10

Hi-C

Lieberman-Aiden, Erez, et al.
"Comprehensive mapping of long-range interactions reveals
folding principles of the human genome." Science 326.5950 (2009): 289-293.

全領域 vs. 全領域の近接性をすべてシーケンシングで決定してしまう。
ビオチンプルダウンで、ライゲーションジャンクションが形成されている断片のみを濃縮。
ペアエンドでシーケンスしなければ意味がない。Forward, Reverse のリード
がリファレンスゲノムへマッピングされた位置を調べ、それらのゲノム上の領域がもともと空間的に近接していた、と解釈。

11

Hi-Cデータ解析のゴールのひとつ＝コンタクトマップ（接触確率行列）の生成
左図はコンタクトマップを、ヒートマップとして可視化した図。
コンタクトマップは対称行列。
タテとヨコに同じゲノム配列を並べて、
位置 i と位置 j にマップされるペアエ
ンドリードがあったら、行列の(i, j) のカウントをひとつ増やす。
したがって、行列で値が大きい要素は、
その領域間でマップされるペアがたくさん見つかる、すなわち接触確率が高い領域ペアであることを意味する。
Rao, Suhas SP, et al.
"A 3D map of the human genome at kilobase
resolution reveals principles of chromatin looping”
Cell 159.7 (2014): 1665-1680.
12コンタクトマップの見方
どのような３次元構造であったら、下図のようなコンタクトマップが得られるだろうか？

Rao, Suhas SP, et al.
"A 3D map of the human genome at kilobase
resolution reveals principles of chromatin looping”
Cell 159.7 (2014): 1665-1680.

13コンタクトマップの見方
どのような３次元構造であったら、下図のようなコンタクトマップが得られるだろうか？

L2
L1
L3
Rao, Suhas SP, et al.
"A 3D map of the human genome at kilobase
resolution reveals principles of chromatin looping”
Cell 159.7 (2014): 1665-1680.

14コンタクトマップの見方
Topologically Associated Domains (TADs)
コンタクトマップの対角線上に沿ってしばしば見られる正方形のかたち。
ヒストン修飾のパターンや複製タイミングなどと強く相関している構造。
エンハンサーの影響力をひとつのTADの内部に隔離するインシュレータとしての機能を持つ？
Akdemir, Kadir Caner, and Lynda Chin.
"HiCPlotter integrates genomic data withinteraction matrices." Genome biology 16.1
(2015): 198.
15

Topologically Associated Domains (TADs)

Ke, Yuwen, et al.
"3D chromatin structures of mature gametes and structural reprogramming
during mammalian embryogenesis."
Cell 170.2 (2017): 367-381.
16

Topologically Associated Domains (TADs)

Rao, Suhas SP, et al.
"A 3D map of the human genome at kilobase
resolution reveals principles of chromatin looping”
Cell 159.7 (2014): 1665-1680.

Fudenberg, Geoffrey, et al. "Formation ofchromosomal domains by loop extrusion.”
Cell reports 15.9 (2016): 2038-2049.

17

コンタクトマップの見方
どのような３次元構造であったら、下図のようなコンタクトマップが得られるだろうか？

Nagano, Takashi, et al.
“Cell-cycle dynamics ofchromosomal organization
at single-cell resolution.”
Nature 547 (2017): 61–67
18

Hi-Cデータの応用
1. ハプロタイプフェイジング
2. ゲノムアセンブリ
3. メタゲノム（meta3C）

Marbouty, Martial, et al.
"Metagenomic chromosome conformation
capture (meta3C) unveils the diversity ofchromosome organization in microorganisms.”
Elife 3 (2014): e03318.

19

Hi-Cのデメリット①
n 高い解像度のコンタクトマップを得るためには莫大な量のシーケンスが必要
l コンタクトマップを構成する要素の数は、ゲノムを分割するBinのサイズ（解像度）に対して二乗で大きくなる。そのため10倍の解像度のコンタクトマップを得る
ためには100倍のシーケンスが必要になる。

20

Hi-Cのデメリット②
n 集団の平均的構造であること（excl. single cell Hi-C）l 集団内の構造の多様性は3次元モデリングの際に「推論」するしかない

O'sullivan, Justin M., et al.
"The statistical-mechanics ofchromosome conformation capture.”
Nucleus 4.5 (2013): 390-398
21

Hi-Cのデメリット②n

集団の平均的構造であること（excl. single cell Hi-C）l
構造の多様性は相同染色体間でさえ観察される（フェイジングデータが利用できればいいが…）

Rao, Suhas SP, et al.
"A 3D map of the human genome at kilobase
resolution reveals principles of chromatin looping”
Cell 159.7 (2014): 1665-1680.

22

Hi-Cのデメリット③n3Cの原理的に、1対1の領域ペアの接触しか観測できない
l 複数の領域の同時接触はデータから間接的にわかるだけ。
l 新たな手法 Genome Architecture Mapping ならば克服できるかも。

Beagrie, Robert A., et al.
"Complex multi-enhancer contacts captured
by genome architecture mapping.”
Nature 543.7646 (2017): 519-524.
23

本講義の内容n Hi-C解析とは
l Chromosome Conformation Capture の原理
l Hi-Cで何がわかるか？コンタクトマップの見方
n Hi-C解析の流れ（実習と並行）
l Hi-C解析のツール
l マッピング
l フィルタリング
l 正規化
l ピーク検出、TAD検出など
l 3D構造モデリング24

Hi-C解析の流れ、利用可能なツール（一部）
クオリティコントロール
マッピング
フィルタリング
正規化

Trimmomatic, cutadapt, fastqc …

Bowtie2, BWA, Juicer, hiclib, HiCUP, HIPPIE…
Juicer, hiclib, HiCUP, HIPPIE, HOMER
Juicer, hiclib, HIPPIE, HOMER

ピーク（ループ）検出

Fit-Hi-C, GOTHiC, HOMER, HIPPIE, HiCCUPS

TAD 検出HiCseg, TADbit, Arrowhead, TADtree, Armatus

3次元モデリング

ChromSDE, ShRec3D, PASTIS …

25

1_mapping_read_to_genome
2_filtering_reads
Bio-Linux-8.0.7_hm_kh.ova を起動。 3_normalization
4_convert_Juice
すべて、~/HiC の中で実行する。
5_detect_TADs
6_modeling_3D
解析のステップごとに、実行する
pythonスクリプトが入ったディレク
トリ（実行結果はそれぞれのResultsの中）data解析に使うfastqファイル
Index
ヒトゲノムのBowtie2インデックス
refヒトゲノム配列（fasta）
src
使用したライブラリのソースコード

実習の内容

26

うまく動かない、どうしても失敗する場合
各ステップのディレクトリの中に、それぞれ “Results” というディレクト
リがあります。
その中に、そのステップで生成されるはずの正解データが入っているので、どうしても失敗する場合はそれをひとつ上の階層に mv していただ
ければ、次のステップ以降の解析も続けられます。

27

実習で使用するデータ
In situ Hi-CでKilobase解像度に達したヒトHi-Cのランドマーク的な論文。
100以上のサンプル、各サンプルあたり数億ペアエンド（１サンプルのfastqでも100GB近いファイルサイズ）
Rao, Suhas SP, et al.
"A 3D map of the human genome at kilobase resolution reveals principles
of chromatin looping” Cell 159.7 (2014): 1665-1680.

28

~/HiC/data
今回は Rao, et al. 2014 の公開データの中の１サンプル（Human B-Lymphocyte: GM12878）のみ、
さらに1000万リードにダウンサンプリングしたデータを扱う。
$cd ~/HiC/data$ls –l
すると、R1, R2 それぞれの fastq ファイルがある。

29

~/HiC/Ref
ヒトリファレンスゲノム（hg19）
http://hgdownload.cse.ucsc.edu/downloads.html からダウンロード
染色体ごとのFASTAファイル

30

~/HiC/Index~/HiC/Ref のヒトリファレンスゲノム（hg19）を
bowtie2-build したもの（Bowtie2のインデックスファイル）

31

Hi-C解析の流れ
クオリティコントロール

Trimmomatic などを使用してアダプター除去、低クオリティ除去など
今回は省略

マッピング
フィルタリング
正規化
ピーク（ループ）検出
TAD 検出
3次元モデリング32

Hi-C解析の流れ
クオリティコントロール
マッピング$cd ~/HiC/1_mapping_read_to_genome

フィルタリング
正規化
ピーク（ループ）検出
TAD 検出
3次元モデリング33

Illuminaのペアエンドシーケンス

http://assets.illumina.com/content/dam/illuminamarketing/images/technology/paired-end-sequencing-figure.gif

34

Hi-Cリードの特徴

Lieberman-Aiden, Erez, et al.
"Comprehensive mapping of long-range interactions reveals
folding principles of the human genome." Science 326.5950 (2009): 289-293.

35

Hi-Cライブラリの特徴
ライゲーションジャンクションは、インサートのどこにでも生じうる…R1, R2 それぞれ、リード全体がマッピング可能…R1がキメラリードとなっている…R2がキメラリードとなっている

36

Hi-Cリードマッピングの際の注意点
1. キメラリードを考慮する
2. ペアのマッピング方向や、インサートサイズを仮定するようなマッピ
ングはしない=> R1, R2 それぞれ個別に、キメラを考慮しつつマッピングする

Imakaev, Maxim, et al.
"Iterative correction of Hi-C data
reveals hallmarks of chromosome
organization.”
Nature methods 9.10 (2012): 999-1003.

37

マッピング戦略
1. R1, R2 個別にマッピングし、マッピング結果をパースして一対一の座
標ペア情報にまとめる。
R1, R2のマッピングのパターンは以下の3通り
I.
R1, R2それぞれリード全体がマップされる
それぞれのマッピング位置間でコンタクトがあった、とみなして座標ペア情報を記録する。
II. どちらかがキメラ
a. 一方がLocusA, LocusBのキメラ、もう一方がLocusB周辺の場合は、Locus A – B 間でコンタクトがあった、とみなして座標ペア情報を記録
b. 上記以外。破棄する。
III. 両方キメラ破棄する。
2. Iterative alignment method （今回はこっちの手法）

38

Iterative alignment method

Imakaev, Maxim, et al.
"Iterative correction of Hi-C data
reveals hallmarks of chromosome
organization.”
Nature methods 9.10 (2012): 999-1003.

39

$less mapping.py

R1をマッピング

R2をマッピング
40

初期ステップのトリミング長
次ステップで何bp延長するか

今回は、時間の関係上35bpまでで打ち切り。本当はリード全長に達するまでやる。
41

$python mapping.py

第一ラウンドBowtie2結果

第二ラウンドBowtie2結果

42

結果$ls –l ../data43

個別にマッピングした結果を統合する$less parse_results.py
統合結果を出力するファイル（HDF5形式）
ゲノムオブジェクトのロード。どの染色体を使うかを指定。

ゲノムの制限酵素消化断片にマッピングされたリードをアサインするため、実験で用いた制限酵素を指定する。
指定できる制限酵素は、BiopythonのRestrictionクラス
44

個別にマッピングした結果を統合する$python parse_results.py$ls -l

HDF5はバイナリファイルなので、中身を見たい場合はHDFViewなどのツールを使うか、pythonのHDF5モジュールなどで開く。
少なくともどちらかのリードがマッピングされたペアについて、座標情報などが格納されている。

45

Hi-C解析の流れ
クオリティコントロール
マッピング
フィルタリング$cd ~/HiC/2_filtering_reads

正規化
ピーク（ループ）検出
TAD 検出
3次元モデリング46

マッピングされたペア情報のフィルタリング
R1, R2 がマッピングされたすべてのペアが、Hi-Cとしての妥当な情報を持つわけではない。
以下のようなパターンでマッピングされたペアは「コンタクト」の情報を持たないため除去する。

Imakaev, Maxim, et al.
"Iterative correction of Hi-C data
reveals hallmarks of chromosome
organization.”
Nature methods 9.10 (2012): 999-1003.
47

$less filtering.py

ゲノムデータのロード

出力ファイルの指定
maximumMoleculeLengthは、イル
ミナライブラリの断片長の情報から設定（今回の実験では400bp）
隣接した制限断片が再びライゲー
ションした場合のペアを除去するために使う

さきほど作ったHDF5ファイ
ルのロード

48

$less filtering.py

追加で行うフィルタリング
filterRsiteStart():
おそらくライゲーションに失敗した
DNA断片
filterDuplicates():
ペアのどちらも同一の座標にマッピングされる２つのペアはPCR
duplicateの可能性が非常に高いため、除去する
filterLarge():10^5bp以上の制限断片にマップされるペアを除去。（リピート領域など、
アセンブル精度が低い領域）
filterExtreme():
マップされるリード数がトップ
0.5%の制限断片を除去。アーティ
ファクトの可能性が高い（元論文参照）

49

$less filtering.py

ゲノムを1MbpごとのBinに分割。
各制限断片上のマッピングされたリード数の情報から、raw read countのコンタ
クトマップ（ゲノム対称行列）としてデータをまとめる。
1Mbpで分割する場合、ヒトゲノムなら約3,000 約3,000 のサイズのマップとなる。
Binサイズの決定に定量的な基準はない。
引き出したい生物学的解釈によって適切に決めるしかない。
（「高品質」なコンタクトマップの基準としてたとえば：マトリックス内の90%以上の値が非ゼロであること、かつ、80%以上で1000以上のコンタクトがあること）

50

フィルタリングを実行する$python filtering.py$ls -l51

フィルタリングを結果を見る$less ./statistics.txt

52

Hi-C解析の流れ
クオリティコントロール
マッピング
フィルタリング
正規化$cd ~/HiC/3_normalization

ピーク（ループ）検出
TAD 検出
3次元モデリング53

データ正規化の必要性
1. サンプル間で比較する場合、サンプルごとにライブラリサイズが異なる。マッピングされたリードの数（サンプルのクオリティ）も異なる。
2. ゲノムの領域ごと、さらには領域間によっても、コンタクトが観測される確率が異なる
Hi-C実験は様々なバイアスの影響で、ある領域間のペアが観測されやすかったりされにくかったりする。
I. 制限酵素断片の長さ。両方とも長い断片の場合、両方とも短い場
合、あるいは長い断片と短い断片のペアはライゲーションが起きにくい。共に中間的な長さの場合にLigationされやすい。
II. 制限酵素断片のGC含量。シーケンシングのバイアス（読み取られやすさ）にばらつきがある。
III. Mappability. マッピングされるリードのゲノム中での「ユニークさ」。その領域がゲノム上でユニークな塩基配列であるかに依存する。
他の実験ではどうやって正規化しているか？
ChIP-seq: INPUTのデータで割り算
RNA-seq: そもそも1サンプル単独で評価しない。サンプル間の比較。
Hi-C実験にはコントロールがないことが問題。
54

Hi-C正規化の方法
1. Explicitにバイアスを仮定する手法
制限酵素断片長、GC含量、マッパビリティなど、バイアスを列挙（それぞれゲノム配列のみから計算可能）、領域ペアの観測確率をそれらのバイアスすべてをパラメータとした確率モデル（ポアソン、負の二項分布など）で表現し、観測値からバイア
スパラメータを学習する。
Yaffe and Tanay 2011 、HiCNormなど

2. Implicitにバイアスを仮定する手法
こちらの方が広く使われている。
Vanilla coverage, ICE, Knight and Ruiz 2012など

55理想的な（正規化された）コンタクトマップでは、ゲノム上のカバレッジが一定
Raw heatmap

Raw coverage

Normalized heatmap

Corrected coverage56各列の和を計算して割り算すると…
k

l

! 𝐴#$

! 𝐴#&

#

#

57

行列の対称性が崩れてしまう
k

l

×

1
∑# 𝐴#&kl
1
×∑# 𝐴#$58

そこで、行の和と列の和の積で割り算する
1
×∑# 𝐴#$ ∑# 𝐴#*k= Vanilla coverage normalization

1×∑# 𝐴#$ ∑# 𝐴#+

1
×∑# 𝐴#$ ∑# 𝐴#&

l

! 𝐴#$
#

59

Vanilla coverage normalizationの仮定
領域 i と領域 j のペアを観測する際のバイアスは、
領域 i のバイアスと、領域 j のバイアスの積に比例する。
つまり、それぞれのバイアスが観測に独立に影響する、と仮定している。
各領域のバイアスは、GC含量やマッパビリティなど、様々な要因が重ね合わさった結果として生じる複合的なバイアス（implicit bias ）強い仮定であるが、Explicit バイアスを仮定して推定した結果ときわめてよく一致する。

Imakaev, Maxim, et al.
"Iterative correction of Hi-C data
reveals hallmarks of chromosome
organization.”
Nature methods 9.10 (2012): 999-1003.

Explicit

Implicit60

Iterative correction (ICE method)
単独の Vanilla coverage normalizationは補正が強すぎる。
（和が非常に小さい列では、割り算結果が爆発する）
Þ Vanilla coverage normalization を何回も適用し、行列全体が収束するまで計算する
このような行列の補正手法は、 “matrix balancing” と呼ばれ、歴史的に何度も再発明されてきた。
ICEと同様の matrix balancing 手法だが、より収束の早いKnight and Ruiz 2012もよく使われる。

61

$less normalize.py
ゲノムデータのロード

Raw read count コンタクトマップのロード
正規化に悪影響を与える可能性のあるBinを除去する。

ICE正規化

正規化後のコンタクト
マップを出力62

正規化を実行$python normalize.py
以下のような結果ファイル（heatmap.pdf）ができているはず。

63

あとでTAD検出、3Dモデリングに使用するため、19番染色体の領域だけ切り出しておく。
$python submatrix.py$less norm_mat.txt

64JuiceBox で、コンタクトマップをインタラクティブに可視化して見る。
JuiceBoxは、独自形式で保存されたコンタクトマップデータを可視化するので、ここまでの実習で生成したデータをJuiceBoxの形式に変換する必要がある。
全ゲノムを見るのはメモリ的にきついので、ここでは一番染色体だけ見る$cd ~/4_convert_Juice$less convert_to_JuiceText.py$python convert_to_JuiceText.py$less ./forJuice.txt$./convert_to_JuiceHiC.sh
以上で、test.hic というバイナリファイル（Juice
形式のコンタクトマップを格納したファイル）ができているはず65

JuiceBoxの実行$./execute_Juicebox.shFile => Open => Local
から、いま作成した test.hic を開く。
Chromosomesで拡大。
Annotationsから
ENCODEデータとの比較なども可能。

66

Hi-C解析の流れ
クオリティコントロール
マッピング
フィルタリング
正規化
ピーク（ループ）検出

実習は省略。手法の紹介のみ。

TAD 検出
3次元モデリング67

コンタクトマップ上のピーク検出＝特に相互作用の強い領域ペアを特定する

L2
L1
L3

Rao, Suhas SP, et al.
"A 3D map of the human genome at kilobase
resolution reveals principles of chromatin looping”
Cell 159.7 (2014): 1665-1680.

68ピーク検出手法によって、得られるピークの数や位置は大きく異なる

Forcato, Mattia, et al. "Comparison of computational methods for Hi-C data
analysis." Nature methods 14.7 (2017): 679.
コンタクトマップの解像度も大きく影響する。
それぞれのピーク検出ツールが、どのように「バックグラウンド」を仮定しているかをちゃんと理解することが重要。

69

Fit-Hi-C (Global background)

Ay, Ferhat, Timothy L. Bailey, and William Stafford
Noble. "Statistical confidence estimation for Hi-C data
reveals regulatory chromatin contacts.”
Genome research 24.6 (2014): 999-1011.

ゲノム上の距離の関数として観測されたリードカウントをスプライン関数でモデリングする。
最初のスプラインは、外れ値を除去するために使われる。
その後、外れ値以外を使って、より洗練されたスプラインをモデリングする。これがヌルモデルとなる。
ヌルモデルの値（ある距離で期待されるリードカウント）を、ICE正規化手法で算出されたバイアスの値も加味して、ある距離のリードカウント観測期待値を計算する。
最後に、期待値と実際の観測値について、二項分布でp-valueを計算（多重検定補
正）する。
70

HiCCUPS (Local background)

周辺の相互作用強度と、K&R
正規化で算出されたバイアス
値からポアソン分布のパラメータを計算し、ピーク位置のp-valueを求める。
Rao, Suhas SP, et al.
"A 3D map of the human genome at kilobase
resolution reveals principles of chromatin looping”
Cell 159.7 (2014): 1665-1680.

71

Hi-C解析の流れ
クオリティコントロール
マッピング
フィルタリング
正規化
ピーク（ループ）検出
TAD 検出$cd ~/HiC/5_detect_TADs

3次元モデリング72

Topologically Associated Domains (TADs)の検出

Forcato, Mattia, et al. "Comparison of computational methods for Hi-C data
analysis." Nature methods 14.7 (2017): 679.
得られるTADのサイズや数はツールによってさまざま。
異なる解像度のコンタクトマップでも安定して同様のTADが得られるか、などを検討することが大事。

73

TADtreeCaleb Weinreb, Benjamin J. Raphael; “Identification ofhierarchical chromatin domains”, Bioinformatics,
Volume 32, Issue 11, 1 June 2016, Pages 1601–1609

ネストしたTAD, sub-TADを検出するPythonスクリプト。
インプットは正規化したコンタクトマップ。

74

TADtreeCaleb Weinreb, Benjamin J. Raphael; “Identification ofhierarchical chromatin domains”, Bioinformatics,
Volume 32, Issue 11, 1 June 2016, Pages 1601–1609

TAD内のコンタクトは距離に線形に増加していくが、他のTAD内部に含まれている場合、その増加率が大きくなる、という観察結果に基づいた手法。動的計画法でベストなTAD階層構造を特定する。
75

$less ./control_file.txt
TADの最大許容サイズ（Binいくつまでか）
高解像度データでは当然大きくすべき。
いじらないことが奨励されてるパラメータ

TAD境界に関わるチューニングパラメータ。
あまりに多くのTADが同一のポイントから生じている場合、おそらく値が大きすぎる。

検出するTADの最大数。

76

TADtreeを実行する。
$python TADtree.py ./control_file.txt
結果が、./output/chr19 以下に出力される。
設定したN以下の検出結果がすべて出力されるが、十分な数のTADを検出している結果が適切なため、proportion_duplicates.txt ファイルを参照
して、重複したTADが検出され始めた時点の結果をチェックする。
結果はBED形式で出力される（ただしポジションの数字の単位は塩基ではなくBinの数であることに注意）。

77

TAD検出後の解析検出されたTAD境界の位置で頻繁に見つかるモチーフや、特定のDNA結合タンパク質の結合パターンを見つける、など。
サンプル間でそれぞれTADを検出し、TADの変化と周辺の遺伝子発現の変化を比較する、など。

Forcato, Mattia, et al.
"Comparison of computational
methods for Hi-C data analysis.”Nature methods 14.7 (2017): 679.

Ke, Yuwen, et al.
"3D chromatin structures of mature
gametes and structural reprogramming
during mammalian embryogenesis."
Cell 170.2 (2017): 367-381.

78

Hi-C解析の流れ
クオリティコントロール
マッピング
フィルタリング
正規化
ピーク（ループ）検出
TAD 検出
3次元モデリング$cd ~/HiC/6_modeling_3D79

染色体3Dモデルの構築
Hi-C実験で得られたコンタクトマップのデータを、「空間制約」として３次元構造のモデリングに組み込む。
大きく分類して２つのアプローチがある。（Serra, et al. 2015）
1. コンタクトマップを直接３次元構造に変換する解析的アプローチ。
コンタクトマップが単一のコンセンサス的な構造を表現していると仮定する。したがって、実際はシングルセル実験のデータによりふさわしい。
2. コンタクトマップを制約として、モンテカルロサンプリングやベイズ法でそれを満たす構造を探す最適化アプローチ。これはさらに２つのカテゴリに分類でき、
a. シミュレーションがそれぞれ独立に、構造の「集合」を生成する方法。得られた構造の集合が集団内多様性を表現する。
b. アンサンブルベースの方法。多くの構造を同時にシミュレート
して、コンタクトマップ制約を満たす構造集合を探索する。
80

染色体3Dモデルの構築
解析的アプローチは次の２つのステップからなる。
1. コンタクトマップをユークリッド距離行列など、「距離の性質」を満たした行列データに変換する。
コンタクトマップで表現されているデータは、多くの異なる構造の「平均」であるため、必ずしも距離の性質のひとつ「三角不等式」
が満たされない。
https://ja.wikipedia.org/wiki/三角不等式

2. 距離行列を満たす直交座標系上の３次元構造をなんらかの最適化手法などで見つける。

81

コンタクトマップから距離行列への変換
コンタクトマップ（値が大きいほど接触頻度が高い）

距離行列（値が大きいほど距離が離れている）

82

コンタクトマップから距離行列への変換
コンタクトから距離へいかに変換するか？
単純なやり方は、単に逆数をとる。
𝐷# ,. =

1
𝐴# ,.

1

α=1とする場合が多い（が、それは自明ではない）
しかし、この変換だと Ai,j = 0 （コンタクトが観測されなかった）場合、
距離 Di,j が無限大になってしまう。
スパースなコンタクトマップでは致命的。
コンタクトマップのすべての要素に適当に小さな値を足して嵩上げするか？ => ほとんどの領域間で適当に足した値の逆数が支配的になってしまうÞ Shortest-path reconstruction
ShRec3Dの手法（ Lesne, et al. 2014）。
本実習で再現（公開ツールはMATLAB）。
83

Shortest-path reconstruction
グラフ理論を応用した距離行列構成手法。
まず、コンタクトマップを構成するBinそれぞれをノードとし、コンタ
クトマップの値の逆数を重みとしたエッジを持つグラフを構成する。
コンタクトマップの値がゼロの場合は、それらのノード間にエッジはひかれない。

84

Shortest-path reconstruction
エッジが直接ひかれているノード間の距
離は、単純にそのエッジの重みとする。

エッジが直接ひかれていないノード間の距離は、そのノード間の「最短パス」を構成するエッジ群の重みの和として定義する。
ノード間の最短パスはFloyd-Warshallア
ルゴリズム（動的計画法）で高速に検索できる。
都合のいいことに、こうして計算した「距離」は三角不等式を満たす。

85

コンタクトマップから距離行列への変換$less ./convert_contact_to_distance.py
Python の NetworkXを使用してグラフ構築、最短パス検索をしている。
今回は、正規化のセクションで生成した19番染色体のコンタクトマップを使って、19番染色体の距離行列を構築してみる。
$python ./convert_contact_to_distance.py ../3_normalization/norm_mat.txt
ディレクトリ内に、dist.npy という距離行列を格納したバイナリファイルができるはず。

86

距離行列から3次元構造への変換
距離行列から空間構造への変換は、もっとも簡単なやり方としては、多次元尺度構成法（Multi-dimensional scaling; MDS）を使う。

メタ16Sの論文で頻繁に出てくるPCoAもMDSの一種。
上図やメタ16S論文では距離行列から二次元平面上の配置への変換に
MDSを用いているが、本来は変換先の空間は何次元でもOK。
今回の実習では3次元。
87距離行列から3次元構造への変換$less ./modeling_3d.py
さっき作ったdist.npy（距離行列データ）を内部でロードして、MDS計
算、得られた３次元座標に基づく染色体の可視化までを行なっている。
実行方法は以下。
$python modeling_3d.py

マウスドラッグでぐりぐり動かすことができる。

88

参考文献
•

Dekker, Job, et al. "Capturing chromosome conformation." science 295.5558 (2002): 1306-1311.

•

de Wit, Elzo, and Wouter de Laat. "A decade of 3C technologies: insights into nuclear organization." Genes &development 26.1 (2012): 11-24.

•

Ke, Yuwen, et al. "3D chromatin structures of mature gametes and structural reprogramming during
mammalian embryogenesis." Cell 170.2 (2017): 367-381.

•

Lieberman-Aiden, Erez, et al. "Comprehensive mapping of long-range interactions reveals folding principles of
the human genome." science 326.5950 (2009): 289-293.

•

Rao, Suhas SP, et al. "A 3D map of the human genome at kilobase resolution reveals principles of chromatin
looping." Cell 159.7 (2014): 1665-1680.

•

Akdemir, Kadir Caner, and Lynda Chin. "HiCPlotter integrates genomic data with interaction matrices." Genome
biology 16.1 (2015): 198.

•

Fudenberg, Geoffrey, et al. "Formation of chromosomal domains by loop extrusion." Cell reports 15.9 (2016):
2038-2049.

•

Nagano, Takashi, et al. “Cell-cycle dynamics of chromosomal organization at single-cell resolution.” Nature 547(2017): 61–67•

Marbouty, Martial, et al. "Metagenomic chromosome conformation capture (meta3C) unveils the diversity ofchromosome organization in microorganisms." Elife 3 (2014): e03318.

•

O'sullivan, Justin M., et al. "The statistical-mechanics of chromosome conformation capture." Nucleus 4.5(2013): 390-398.

•

Beagrie, Robert A., et al. "Complex multi-enhancer contacts captured by genome architecture
mapping." Nature 543.7646 (2017): 519-524.
89•Imakaev, Maxim, et al. "Iterative correction of Hi-C data reveals hallmarks of chromosome
organization." Nature methods 9.10 (2012): 999-1003.

•

Yaffe, Eitan, and Amos Tanay. "Probabilistic modeling of Hi-C contact maps eliminates systematic biases tocharacterize global chromosomal architecture." Nature genetics 43.11 (2011): 1059-1065.

•

Hu, Ming, et al. "HiCNorm: removing biases in Hi-C data via Poisson regression." Bioinformatics 28.23 (2012):
3131-3133.

•

Knight, Philip A., and Daniel Ruiz. "A fast algorithm for matrix balancing." IMA Journal of Numerical
Analysis 33.3 (2013): 1029-1047.

•

Forcato, Mattia, et al. "Comparison of computational methods for Hi-C data analysis." Nature methods 14.7
(2017): 679.

•

Ay, Ferhat, Timothy L. Bailey, and William Stafford Noble. "Statistical confidence estimation for Hi-C data
reveals regulatory chromatin contacts." Genome research 24.6 (2014): 999-1011.

•

Caleb Weinreb, Benjamin J. Raphael; Identification of hierarchical chromatin domains, Bioinformatics, Volume32, Issue 11, 1 June 2016, Pages 1601–1609•

Serra, François, et al. "Restraint-based three-dimensional modeling of genomes and genomic domains." FEBS
letters 589.20PartA (2015): 2987-2995.

•

Lesne, Annick, et al. "3D genome reconstruction from chromosomal contacts." Nature methods 11.11 (2014):
1141-1143.

90