配列モチーフ探索

文献

Tompa et al. 2005 Nature Biotechnology: Phylogenetic footprinting、発現パターン、ChIP結果などの補助的情報を使わず、配列のみを使って解析する13のツールを fly, human, mouse, yeast, 人工データでテストした。ただし、エキスパートがパラメータをファインチューニングし、ベストヒットした1つのモチーフを取る、という条件での比較。
どのツールもyeastで特にハイスコア。flyはやや低め。人工データよりもrealのスコアが低いのは、今回正解としている以外のモチーフも配列に含まれててそいつをトップヒットとしてしまったことにペナルティがかかってしまったせい。
Hu et al. 2005 BMC Bioinformatics: あまりチューニングせずほぼデフォルト設定でベンチマーク。そうすると、どのツールも単体では意外とショボい。複数のツールを組み合わせて使うと良い。
GuhaThakurta 2006 Nucleic Acids Res: モチーフ探索の基本的な流れがわかりやすい
Das and Dai 2007 BMC Bioinformatics: 実験はしていないが網羅的なレビュー
Zambelli et al. 2013 Briefings in Bioinformatics: ChIP時代のツールも含めてざっくり

基本

配列モチーフとは

一般的に: だいたい 5–20bp くらいの生物学的に意味のある塩基配列パターン。 e.g. 転写因子結合サイト、タンパク質の機能ドメイン
palindoromic motifs: 相補鎖同士を5’から3’に読んで同じになってるような配列。 e.g. CACGTG, GCATGC
spaced dyad (gapped) motifs: 3–5bpくらいの保存的な結合サイトの間に非保存的なspacer/gapがある。結合タンパクがdimerな場合とか。 e.g. Yeast Gal4: 5’-CGGnnnnnnnnnnnCCG-3'
Shine-Dalgarno sequence: 原核生物mRNAの開始コドン上流にあるリボソーム呼び込みモチーフ。 16SリボソームにアンチSD配列がある。 e.g. E. coli AGGAGGA
Kozak sequence: 真核生物版SD配列。

見つけ方

Signature

共発現 coexpression/coregulation
高頻度 overrepresented
保存性 conserved -> phylogenetic footprinting

したがって、戦略としては

1ゲノムの中で共発現してる遺伝子のプロモーターを集めて高頻度な部分配列を探す
複数種のorthologousな遺伝子のプロモーターを集めて特に保存的な部分配列を探す
上記2つの合わせ技

表し方

Consensus

1つの文字列でバシッと表現。正規表現の文字集合 [] を使ってポジション毎の複数候補を表現したり、 IUPAC命名法に従って degenerate symbol で表すこともある。 https://www.bioinformatics.org/sms/iupac.html

e.g TATAAT, TATA[AT]A[AT], TATAWAW

PSPM: Position-Specific Probability Matrix

ポジションごとの塩基・アミノ酸の相対的な出現頻度をそのままの整数、あるいは合計1になるような実数[0, 1]の行列で表示。 Profile とも呼ばれる。 (Bucher 1990)

A     C     G     T
0.625 0.125 0.250 0.000
0.375 0.375 0.000 0.250
0.125 0.375 0.500 0.000
0.250 0.375 0.250 0.125
0.000 0.250 0.250 0.500
0.375 0.375 0.125 0.125
0.250 0.125 0.375 0.250
0.250 0.125 0.250 0.375

PSSM: Position-Specific Scoring Matrix

頻度からいろんなスコアに変換して行列にする。ツールによって異なるのでよく分からん。

PWM: Position Weight Matrix

一番良く見る呼び方だが、PSPMと同義だったり、PSSMと同義だったり。。。

Sequence Logo

ポジションごとの保存性・確実さと各塩基・アミノ酸の寄与をアルファベットの大きさで視覚的に表現。 (Schneider and Stephens 1990)

ポジション i における塩基・アミノ酸 a の高さは、相対頻度 f と情報量 R の積:
$\text{Height}(a, i) = f(a, i) \times R(i)$

ポジション i における情報量 (information content, IC) は、定数ひく不確実性ひく補正項:
$R(i) = \log_2(s) - H(i) - e(n)$
ただし s = 4 [DNA] or 20 [Protein]

ポジション i における不確実性 (Shannon entropy):
$H_i = - \sum _{k}{f(k, i)} \times \log_2 f(k, i)$

配列の数 n が少ない時のための補正項:
$e(n) = \frac{1}{\ln 2} \times \frac{s - 1}{2n}$
ただし s = 4 [DNA] or 20 [Protein]

Consensus Logo

コンセンサス配列の重み付け表示版、あるいはSequence Logoのトップヒット限定版。 1行の文字列なので図を使わず書式付きリッチテキストとして扱える。

アルゴリズムの分類

一覧

Das and Dai 2007 表を改変

MEME Suite: MEME, DREME, MEME-ChIP, MAST
MoD Tools: Weeder, WeederH, Pscan, PscanChIP
RSAT (Regulatory Sequence Analysis Tools): Oligo-Analysis, Dyad-Analysis, Consensus, peak-motifs; ローカル使用には作者へのメールが必要。
Brutlag Bioinformatics Group: BioProspector, MDscan
SCOPE: BEAM, PRISM, SPACER を組み合わせる。ローカル使用には作者へのメールが必要。

word-based (string-based)

オリゴヌクレオチドの数え上げと頻度比較 = exhaustive enumeration

全通りやるので global optimality にたどり着く

短いモチーフ配列(=真核生物)向き

Oligo-Analysis, Dyad-Analysis

Weeder

MITRA,

YMF

PSM: probabilistic sequence model

長いモチーフ配列(=原核生物)向き

探索の仕方によっては local optimum に陥りがち

EM (expectation maximization) algorithm

w 行のPWM $\theta$ を適当に作る (位置 j が塩基 a である確率は $\theta_{ja}$)
配列が n 本あり、それらの部分配列 S が $\theta$ から生成される確率(すなわち尤度)は

\[ L(\theta) = P(S_1, ..., S_n|\theta) = \prod_i^n \prod_j^w \theta_{i, S_{i j}} \]
部分配列がモチーフ $\theta$ から生成された場合と、そうじゃないただのバックグラウンド $\theta_0$ から生成された場合を

\[\begin{aligned} z_i &= 1 [\text{if} S_i \text{is generated by} \theta]\\ z_i &= 0 [\text{if} S_i \text{is generated by} \theta_0] \end{aligned}\]

のように missing parameter として表すと、尤度は

\[ L(z, \theta, \theta^0) = \prod_i^n [z_i P(S_i|\theta) + (1 - z_i)P(S_i|\theta^0)] \]
EMアルゴリズムでは以下の平均対数尤度を使う。

\[ \log \tilde{L}(z, \theta, \theta^0) = \sum_i^n [q(z_i = 1) \log P(S_i|\theta) + q(z_i = 0) \log P(S_i|\theta^0)] \]

ただし $q(z_i)$ はラベル変数の事後分布で、

\[\begin{aligned} q(z_i = 1) &\sim P(z_i = 1) P(S_i|\theta)\\ q(z_i = 0) &\sim P(z_i = 0) P(S_i|\theta_0) \end{aligned}\]
E (Expectation) ステップ: $\theta$ を使って $q(z_i)$ を更新
M (Maximization) ステップ: $q(z_i)$ を使って、尤度最大となるよう $\theta$ を更新
EMステップを繰り返す

MEME では各配列が持つモチーフの数が1以外の場合も扱えるのと、最初のPWMの作り方を工夫して大域最適解に行きやすくなっているという点でただのEMより改善されているらしい。

Gibbs sampling

AlignACE (Linux実行形式のみ配布), MotifSampler

N 本の各配列のモチーフ位置の初期値をランダムに与える
ランダムに配列をひとつ選び、それ以外の配列のモチーフを並べてPWMを計算
先に選んだ配列の上でPWMをスライドさせつつスコア (l-merがPWMから生成される確率など) を計算
スコアに比例した確率で新しいモチーフ位置を決める
各配列のモチーフ位置が安定するまで2–4を繰り返す

スコアには生の生成確率ではなく、バックグラウンドスコア(塩基出現頻度から予測される出現頻度)との比を考慮した相対エントロピーを用いることがある。

\[ \sum_i^l \sum_a^{ACGT} p_{ai} \log_2 \frac{p_{ai}}{b_a} \]

Ensemble

基本的にはどのアルゴリズムもsensitivity, accuracy共に低い。既存のアルゴリズムを複数ensemble的に組み合わせるべし。

EMD

(Hu et al. 2006)

配列 N 本のデータセットに対し M 個の異なるアルゴリズムをそれぞれ R 回ずつ走らせて K 個の予測サイトを得る
配列毎・アルゴリズム毎に予測サイトを集めてグループ化する
配列毎に全アルゴリズムの予測スコアを足しあわせて voting, smoothing, extracting

PF: Phylogenetic Footprinting

古くは Tagle et al. 1988 がこれを用いてグロブリンの調節領域を予測した。

ClustalW などによるアラインメントでは、近すぎると情報量が無く、遠すぎるとメチャクチャ。

複数種の配列で保存的な部分配列を探すだけなら上記のような1 genome用のツールを使ってもある程度はいけるが、そのまま使うのではなくて系統関係による重み付けが必要。

Backgroundの取り方

統計的仮説検証のためにランダム生成DNA配列を使ってしまうと、繰り返しが多く含まれていたりする現実のDNA配列と比べて “too null” になり、検出力を overestimate することにつながる。かといって、そういう配列を予めマスクしてしまっていいかというとそれも微妙 (Simcha et al. 2012)

ランダムではなくバックグラウンドと比べて多いかどうかを明示的に言いたい場合は discriminatory motif と呼んだりする。

Markov process: 各塩基の頻度を考慮して配列をランダム生成するのが0階(単純)マルコフ過程。 1つ前の塩基を考慮して次の塩基の出方が影響を受ける (つまり塩基ペアの頻度情報をつかう)のが1階マルコフ過程。て感じで配列を生成するマルコフ過程のオーダーを高くしていく。 Thijs et al. 2001 BioProspector, MotifSampler, YMF, MEME
バックグラウンド配列セットにおけるモチーフ出現頻度を考慮: Co-Bind 2001, WordSpy 2005, ANN-Spec 2000

見つけた後の評価

STAMP, TOMTOM, MotIV

ツールの性能評価

ツールを作るのも難しいが、性能を正しく評価するのも難しい。

ツール開発論文の多くは「MEMEよりxx%も精度が向上した」などと報告しているが、それはデータ・評価項目・コマンドオプションなどに依存している。

テストデータ

現実の配列を使う: 真の正解が分からない。つまり、未知のモチーフを発見すると不当なペナルティを食らう。
人工的に生成した配列を使う: 自然の配列ができあがる真の確率過程が分からないので、特定のアルゴリズムに有利になったりするような偏りが生じるかも。

評価項目

Raw Score (塩基単位): 予測されたモチーフを正解モチーフと塩基単位で比較し、正解 (True Positive, True Negative) と間違い (False Positive, False Negative) をカウントする。
nTP, nTN, nFP, nFN:
Raw Score (サイト単位): 予測されたモチーフが正解モチーフの1/4以上オーバーラップしてたら当たりとし、モチーフ(サイト)単位で TP, FP, FN をカウントする。 (TN は数えられない)
sTP, sFP, sFN
Sensitivity (塩基 or サイト単位): 正解モチーフ配列のうち、どれくらい拾えたか
$Sn = \frac{TP}{TP + FN}$
Specificity (塩基 or サイト単位): 本当はモチーフじゃない配列に、どれくらいノーと言えたか
$Sp = \frac{TN}{TN + FP}$
Positive Predictive Value (塩基 or サイト単位): モチーフだと予測されたサイトに、どれくらい正解が含まれていたか (これをSpecificityと呼ぶことも多い)
$PPV = \frac{TP}{TP + FP}$
Performance Coefficient (塩基単位): 正解率 (Pevzner and Sze 2000)
$nPC = \frac{TP}{TP + FN + FP}$
Correlation Coefficient (塩基単位): 予測と正解のピアソン相関係数 (Burset, M. and Guigó 1996)
$nCC = \frac{TP \times TN - FN \times FP}{\sqrt{(TP + FN)(TN + FP)(TP + FP)(TN + FN)}}$
Average Site Performance (塩基単位): 感度と精度の平均 (Burset, M. and Guigó 1996)
$nASP = \frac{Sn + PPV}{2}$

モチーフを含まないネガコンをデータに含めた場合は常に $TP = FN = 0$、アルゴリズムがモチーフ無しと予測した場合は常に $TP = FP = 0$ となり、いろんな尺度が undefined/noninformative になってしまうのが悩みどころ

評価の要約方法

ツールによってデータや尺度の得意・不得意があるので、複数のデータセットでいろんな尺度を計算することになる。それをどのように要約するか。

Averaged: 普通に算術平均
Normalized: 全アルゴリズムの平均からの残差を標準偏差で割りデータセットに対する平均を取る
Combined: それぞれのデータで尺度を計算してしまうのではなく、 7つの Raw Score をそれぞれ全データセットで足しあわせて、それから各尺度の計算をする。ゼロ除算を避けられる場合が多いので有用。

Database

DNA

TRANSFAC (Matys et al. 2003): 真核生物。情報量は多いが1つのTFが複数のプロファイルを持ってる。 Professional版と機能限定Public版がある。
JASPAR (Sandelin et al. 2004): 真核生物。情報量は少ないが1つのTFは1つのプロファイルだけ持ってる。無制限。

SCPD: The Promoter Database of Saccharomyces cerevisiae

DBTBS: Database of Transcriptional Regulation in Bacillus subtilis

RegulonDB:: Escherichia coli K-12 Transcriptional Regulatory Network

Protein

Heavy Watal