ゲ ノ ム 情 報 学 基 礎 論 (2)連鎖不平衡と多因子性

ゲノム情報学基礎論
(2)連鎖不平衡と多因子性
大阪大学 大学院医学系研究科 ゲノム情報学共同研究講座
中谷明弘
Akihiro NAKAYA, Ph.D.
Department of Genome Informatics
Graduate School of Medicine, Osaka University
全体の概要
n 担当 中谷明弘（医学系研究科ゲノム情報学・特任教授）
n 会場 医学系研究科講義棟
講義室2
n 日程 通年・偶数月全5回 18:00∼20:00
1. 
6月16日（火）計算ゲノミクスのための情報処理
2. 
8月25日（火）連鎖不平衡と多因子性（変異データの特性と特徴）
3.  10月27日（火）分類と層別化（検体や変異の階層的クラスタリング）
4.  12月15日（火）特徴量の構築（属性値のモデル化と予測のための学習）
5. 
2月16日（火）数え上げ（形質に関連する因子の組み合わせ最適化）
n 情報ページ
http://www.gi.med.osaka-u.ac.jp/lectures/GIF2015/
本日の内容
n エクソーム解析用VCFファイル内の変異情報
n ローカス間の連鎖不平衡
n 形質関連ハプロタイプと常染色体劣性遺伝
n HDR法による特異的な遺伝子型パタンの探索
n 注釈付き変異情報ファイル
n 計算機実験
1. 
2. 
候補変異の位置情報を含んだファイルの生成
n 注釈付き変異情報ファイルのフィルタリングする
n Perlで作成する
HDRで特異的な領域の探索
n 上記で生成した候補変異のデータを使用する
n 既存のJavaアプリで処理する
エクソーム解析のデータ
n FASTQ（リード情報）→
VCF: Variant Call Format
FASTQ: acronym?, one of FAST-Xs, Q=quality, cf, FASTA
VCF（変異情報）→ 下流の解析
エクソーム解析用VCFファイル内の変異情報
n 変異データ（SNVとInDel）
WT(=hg19)
n 変異はhg19
ヘテロ
ホモ
(WT) を基準
n 日本人の参照配列ではない
n 検体当たり10万弱の変異
n アノテーション（変異の特性情報付加）
n 絞り込み e.g., 低頻度(< 0.1%) && 非同義置換 && ･･･ → 数十∼数百変異
（コントロール）検体群
注目している検体
CTRL群内で低頻度
CTRL群内で高頻度
変異無し＝hg19
CTRL群内で高頻度
1000人ゲノムの日本人検体群
n 1000G
NA18940∼NA18949
赤：ヘテロ
青：ホモ （hg19が基準）
http://www.1000genomes.org からデータ取得
1000人ゲノムの日本人検体群
青・緑：遺伝子領域
赤：エクソン領域
赤：ヘテロ
青：ホモ
1番染色体短腕のテロメア付近
（約50kbpの領域）
NA18949のみ「変異無し」
→ 実はこちらが「変異有り」?
VCFファイルには含まれないので見落としてしまう?
全検体が「変異有り」
→ 実はこれが「参照配列」?
単一検体の変異情報のみでは
恐らく十分ではない
ここまでのまとめ
n NGSによるエクソーム解析のデータ（FASTQ:リード
→ VCF:変異）
n 特定の検体での変異の有無だけ見ても十分ではない
『注目している検体に変異があった！』
1. 
日本人に共通した「変異」かも知れない（hg19基準）
→ 実はそれが日本人の「参照配列」かも知れない
2. 
変異が検出されていない位置が変異かも知れない
検体群内あるいはデータベースでの変異の頻度情報との突き合わせが必要
ローカス間の連鎖不平衡
n 連鎖不平衡（LD:
linkage disequilibrium）
同一染色体上の近接した領域が一緒に世代間で伝達
n 形質に関連した（原因の）染色体領域
n 領域内の変異自体が形質の原因
n 領域内の変異が「原因」と連鎖
形
質
の
原
因
ロ
ー
カ
ス
近隣のローカスと連鎖している？
エクソン領域の変異が
タンパクに影響
遺伝的要因がエクソン領域にあるとは限らない
形
質
の
原
因
ロ
ー
カ
ス
近接の遺伝子型パタンと連鎖
直接は影響しない可能性 → マーカー
形質関連ハプロタイプ
n 未知の遺伝性疾患が無限に現れる訳ではない
n 疾患関連ハプロタイプは出尽くしている？
n LDブロックの崩壊と融合が釣り合った定常状態？
n 致死ハプロタイプ
n 有史以来の全ての
日本人 （2000年分くらい？）
n 少なくとも120世代くらい れば･･･
n 特定の遺伝性疾患の
原因 はどこかにあるはず
n 疾患家系のデータ → 遺伝様式
n 常染色体劣性
n 常染色体優性
n ･･･
1000人ゲノムの日本人検体群
青・緑：遺伝子領域
赤：エクソン領域
赤：ヘテロ
青：ホモ
1番染色体短腕のテロメア付近
（約50kbpの領域）
例えばこの領域の変異パタン
バラバラではなく恐らくある程度の
塊 で祖先から引き継がれている
常染色体劣性遺伝の場合
n 遺伝的背景を想定
変異
n IBDとIBS
祖先の誰かに入ってしまった疾患関連変異
→ 近隣の変異パタンも引き連れて受け継ぐ
候補ローカス
※IBD: identical by descent
IBS: identical by state
家系内の疾患検体に特異的かつ共通な変異位置
→ 候補ローカス（頻度等で絞り込み）
→ 周辺領域に特異的な遺伝子型パタン？
周辺領域
n 幾つかの仮定/前提
n エクソン領域の稀少変異
n 疾患検体に特異的
n ホモ接合
非血族婚
血族婚
常染色体劣性ならば発症者でホモ接合しており
尚且つ、近隣もホモ化しているという仮定
n 近接の遺伝子型パタンと連鎖
絞り込まない
変異も使用
必要十分条件？
ホモ
ホモ
ホモ
IBD
ホモ
ホモ位置の相互一致率
ホモ
ホモ
ホモ
ホモ
IBS
形質に対する遺伝的要因の多因子性
n 形質に関連する変異が染色体上で連続しているとは限らない
n 染色体を跨がった
連鎖
n シスとトランス
n 相互作用・交互作用・エピスタシス
→ 組み合わせ論的探索
ここまでとこれからのまとめ
1. 
複数の変異の間の関係
n 染色体上で近隣･･･連鎖不平衡（LDブロック）
n 染色体間で関連･･･交互作用（エピスタシス）
2. 
常染色体劣性の場合 n ホモ接合
3. 
稀少なホモ変異を候補とした形質関連変異の評価
n 候補の絞り込み
n 候補の特異性評価
特異的な遺伝子型パタンの探索
n HDR
method (Hamming distance ratio method)
n Imai A., Nakaya A., and Ott J., et al.,
n Scientific Reports, 2015 Jul 6;5:12028. doi:10.1038/srep12028.
www.nature.com/scientificreports
OPEN
Received: 11 November 2014
Accepted: 11 June 2015
Published: 06 July 2015
Beyond Homozygosity Mapping:
Family-Control analysis based on
Hamming distance for prioritizing
variants in exome sequencing
Atsuko Imai1,2, Akihiro Nakaya2, Somayyeh Fahiminiya3, Martine Tétreault3,
Jacek Majewski3, Yasushi Sakata1, Seiji Takashima1,4, Mark Lathrop3 & Jurg Ott5
HDR 法
n ２つの最適化パラメータ
検体1
n 領域の場所（中心の位置）p
n 領域のサイズ
検体2
v/v
d
d/2
!v/v
n1 ← (v/v,v/v)
n 最適化関数：HDRスコア
n ハミング距離(Hamming
distance)
n 異なった変異ローカスの数
n 全数との比(Ratio)を取って正規化
n2 ← (v/v,!v/v)
p
d/2
HDR = (n2+n3)/(n1+n2+n3)
n3 ← (!v/v,v/v)
(v/v,v/v)が無いとき (n1=0) 最大値 1.0
n = n1+n2+n3
アノテーションパイプライン
VCFに対して解析に必要な情報を蓄積/追加
PLにより変異検出（VCF）
PLにより検体当たり8∼9万変異検出
##fileform at=VC Fv4.1
#C H RO M PO S
chr1
888639
chr1
888659
chr1
899928
chr1
899938
ID
.
.
.
.
REF
T
T
G
G
ALT
C
C
C
C
Q U AL
1105.77
932.77
93.03
11.83
FILTER
IN FO
FO RM AT N A19091_ERR250253_b
PASS
AC =2;AF=1.G00;
T:AN
AD =2;
:D P:
DGP=30;
1|
Q 1:
:O0,
GD30:
:PL
el30:
s=0.
60:
00;
.:1134,
FS=0.90,
000;
0 H aplotypeScore=0.0000;M LEAC =2;M LEAF=1.00;M Q =56.00
PASS
AC =2;AF=1.G00;
T:AN
AD =2;
:D P:
DGP=26;
1|
Q 1:
:O0,
GD26:
:PL
el26:
s=0.
60:
00;
.:961,
FS=0.
78,
000;
0 H aplotypeScore=0.0000;M LEAC =2;M LEAF=1.00;M Q =56.00
Low C overage
AC =2;AF=1.G00;
T:AN
AD =2;
:D P:
DGP=4;
1/1:
Q :PL
D0,
el
4:
s=0.
4:12:
00;
121,
FS=0.
12,0000;H aplotypeScore=0.0000;M LEAC =2;M LEAF=1.00;M Q =39.12;M
Low C overage;
ACVeryLow
=2;AF=1.
QG00;
ual
T:AN
AD =2;
:D P:
DGP=1;
1/1:
Q :PL
D0,
el
1:
s=0.
1:3:38,
00;3,
FS=0.
0 000;H aplotypeScore=0.0000;M LEAC =2;M LEAF=1.00;M Q =60.00;M
Annovarによりアノテーション, 遺伝子名, Exon/Intron情報, 1000G/ESP/ExAC頻度情報を付加
Annovarによる情報付加
フィルタリング
公共DBからの情報抽出
情報追加
情報蓄積
フィルタリング
候補遺伝子抽出
稀少変異抽出
#C hr
chr1
chr1
chr1
chr1
chr1
chr1
Start
End
978604
978605
980460
980460
980537
980537
981087
981087
981931
981931
982444
982444
Ref
CT
G
G
A
A
A
Alt
A
A
G
G
G
Func.refG ene
intronic
intronic
intronic
intronic
exonic
intronic
G ene.refG ene
ExonicFunc.
esp6500si
refG ene _al
1000g2014oct_al
l
ExAC
l _ALL ExAC _EAS Filter
N A19091_ERR250253_b_G
N A19091_ERR250253_b_RefD
N A19091_ERR250
T
AG RN
0.399263
0.479233
0.5609
0.7856 PASS
hom
0
21
AG RN
0.911941
Low C overage
hom
0
4
AG RN
Low Q D ;Lowhet
Q ual;VeryLow Q ual10
2
AG RN
0.832078
0.869409
0.9128
0.9999 PASS
hom
0
16
AG RN
synonym ous 0.
S780216
NV
0.797724
0.8767
0.9649 PASS
hom
0
5
AG RN
0.821885
PASS
hom
0
5
Exonic/Splicing, アレル頻度(e.g., 0.1％未満), Hetero/Homo, ･･･等の条件
→ 数百程度に絞り込み
#C hr
Start
Ref
Alt
Func.refG ene
G ene.refG ene
ExonicFunc.refG ene
esp6500si_all
1000g2014oct_all
ExAC _EAS
snp135
snp137
snp138
chr14
24728975
G
A
exonic
TG M 1
nonsynonym ous S N V
chr14
24879204
G
A
exonic
N YN RIN
nonsynonym ous S N V
0.000399361
0.000399361
.
.
rs121918731
rs182330062
rs121918731
rs182330062
rs121918731
rs182330062
C LIN SIG =pathogenic;C LN D B N =Autos
clinvar_20150330
om al_recessive_congenital_ichthyosis_
Q U AL
222.77
91.77
Filter
PASS
PASS
N A19091_ERR250253_b_G
hetT
het
rs121918731:[W ith P athogenic allele], rs182330062:N A,
snp-info
rs121918731:[W ith P athogenic allele], rs182330062:N A,
rs121918731:[W ith P athogenic allele] rs182330062:N A
pub
H G VD
0.6026
アレル頻度情報：ESP6500, 1000G, ExAC
SNP情報：dbSNP（SNP135/137/138)
疾患関係情報：ClinVar
non-­‐pathogenic probable-­‐non-­‐pathogenic probable-­‐pathogenic pathogenic histocompa2bility 臨床関連性情報：dbSNP
文献情報：PubMed
日本人固有変異情報：HGVD
注釈付き変異情報ファイル
n タブ区切りファイル
0:#Chr, 1:Start, 2:End, 3:Ref, 4:Alt, 5:Func.refGene, 6:Gene.refGene, …
27:esp6500si_all, 28:1000g2012apr_all, 29:snp135, 30:snp137, …
48:Filter, 49:NA18940_ERR034596_b_GT, …
計算機実験
1. 
候補変異の位置情報を含んだファイルの生成
n 注釈付き変異情報ファイルのフィルタリングする
n Perlで作成する
2. 
HDRで特異的な領域の探索
n 上記で生成した候補変異のデータを使用する
n 既存のJavaアプリで処理する 正確にはScalaアプリ
変異情報ファイルのフィルタリング
#!/usr/bin/perl
$filename = $ARGV[0];
open($FILE, $filename);
<$FILE>; # skip the first line
こ
の
ブ
ロ
ッ
ク
を
１
行
読
み
出
す
毎
に
実
行
$count = 0; 条件を満たす行の数を保持する変数。
while(<$FILE>) {
my $line = $_;
chomp $line;
my @vals = split(/\t/, $line);
my
my
my
my
my
my
my
$chr = $vals[0];
$start = $vals[1];
$end = $vals[2];
$func = $vals[5];
$g1000 = $vals[28];
$filter = $vals[48];
$gt = $vals[49];
$ARGV[0]に1つめのコマンドライン引数が入っている。
それを$filenameという変数に保存。
ファイルを開いて$FILEに関連づけ（オマジナイ）。
<$FILE>でファイルから1行読み出せる。
最初の1行（列の名前）は要らないので捨てている。
while(<$FILE>)で1行の読み出しを繰り返すことができる。
読み出した1行は$_に保持される･･･$line変数に保存。
chompで行末の改行文字を削除。
タブ文字（\t）で1行を分解して配列@valsに保存。
0∼n-1番目の要素は$vals[0]∼$vals[n-1]で参照可能。
各列の値を変数に保持（見易さのため）。
各列の値でフィルタリング。
eqは文字列の一致の判定。&&は複数条件の論理積（ かつ ）。
if ($func eq "exonic" && $g1000 < 0.001 && $filter eq "PASS" && $gt eq "hom") {
printf "$chr\t$start\n";
$count++;
条件を満たせば$chrと$startの値をタブ区切りで出力して改行。
} 条件を満たす行の数を保持する変数を1増やす。
}
printf STDERR "$count\n";
変異情報ファイルのフィルタリング
n テキストエディタで入力して
.pl で終わるファイルに保存
n コマンドラインから以下の通りに打ち込んでリターンで実行
% perl FilterAnno.pl NA18940_ERR034596_hg19_multianno.txt > pos.txt
n 末尾で
> pos.txt とすることで結果をファイルに保存可能（リダイレクト）
※ Windowsでは C:¥Perl64¥bin¥perl.exe
HDRによる特異的領域の探索
n 必要な入力ファイル
n ケース検体のVCFファイル（1個）
n コントロール検体のVCFファイル（1個∼複数個）
n 注目するローカスの位置情報ファイル
n 必要な設定パラメータ
n 領域のサイズの範囲（と刻み幅）
http://www.gi.med.osaka-u.ac.jp/software/hdr/
HDRの起動
n 以下のJARファイルをダブルクリック
http://www.gi.med.osaka-u.ac.jp/lectures/GIF2015/dat/JARfiles/hdrfx.jar
n 或いは以下をコマンドラインから実行
Windows
C:¥Program Files¥Java¥jre1.8.0_45¥bin¥java.exe -Xmx8G ‒jar hdrfx.jar
% java -Xmx8G ‒jar hdrfx.jar
Unix系
ファイルを読み込む
Data メニューから
ケースVCF → 赤色タブ
コントロールVCF → 青色タブ
位置情報 → 緑色タブ
1000G NA18940
NA18941 ∼ NA18949
※ 1000G NA18940 ∼ NA18949のVCFファイルは未フィルタリング
→ VCFファイルに当該行があるか否かで遺伝子型のパタンを評価する。
なので、クオリティーの低い物、ホモ、ヘテロ、高頻度、低頻度、
エクソン領域の内外、･･･ の全てを含んだままになっている。
（必要に応じてVCFファイルのフィルタリングを行う）。
いずれのファイルも複数を読み込み可能
ケースVCFと位置情報ファイルを選択
ケースVCF: Case プルダウンメニューからタブ名で選択
この例では1つずつしか読み込んでいないので一択
位置情報: Positions プルダウンメニューからタブ名で選択
※コントロールVCFは読み込んだ全てを使用
探索範囲のパラメータを設定する
探索する最小の領域サイズ（e.g., 10kbp）
探索する最大の領域サイズ（e.g., 30kbp）
刻み幅（e.g., 10kbp）
この設定の場合：10kbp, 20kbp, 30kbp と 10kbpから10kbp 刻みで探索する
探索処理の実行
Exec & Save ボタンで実行 → 結果を保存するファイル名を指定
黄色タブに実行結果のサマリ
パーメータ等を変更して複数回実行可能
探索処理の結果
n 未加工のVCFファイル･･･エクソン領域
( 100bp) の変異
n 結果ファイル
n 前半：各ポジション・領域サイズでの全検体間のスコア
n 後半： 〃 での全検体間の
n1+n2+n3
chr6 32489856
黄色タブの実行結果のサマリ
検体1 検体2
d
chr6
chr6
chr6
chr6
chr6
chr6
32489856
32489856
32489856
32489856
32489856
32489856
n1+n2
スコア
10
10
10
10
10
10
1
1
1
1
1
1
n1+n3
n2+n3
2 0.919 16 64 68 74
3 0.950 16 47 57 60
6 0.901 16 73 73 81
7 0.923 16 54 60 65
9 1.000 16 23 39 39
10 0.891 16 55 57 64
n1+n2+n3
変異の順位と有意さ
n maxstatRSプログラム
http://www.jurgott.org/linkage/maxstatRS.html
まとめ
n 連鎖不平衡や交互作用のような複数の要因の扱い
n 手作業や目視の限界
n 定型化できない処理
n お薦めのパタン
n Perlでのプログラミング 入力/出力ファイルの整形やジョブ管理
n 既存のソフトウェア活用 面倒な部分はまずは他力本願 n HDR
n NGSデータの解析処理のテンプレート
n エクソーム解析PLと連携（コマンドライン版）
n 下流の統計解析との統合
とはいうものの複雑なことも十分可能
# 実はこれもPerlプログラムの出力