Heavy Watal

配列モチーフ探索

文献

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

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

表し方

Consensus
1つの文字列でバシッと表現。 正規表現の文字集合 [] を使ってポジション毎の複数候補を表現したり、 IUPAC命名法に従って degenerate symbol で表すこともある。 http://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]

http://schneider.ncifcrf.gov/paper/logopaper/

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
  1. w 行のPWM $\theta$ を適当に作る (位置 j が塩基 a である確率は $\theta_{ja}$)

  2. 配列が n 本あり、それらの部分配列 S が $\theta$ から生成される確率(すなわち尤度)は

    \[ L(\theta) = P(S_1, ..., S_n|\theta) = \prod_i^n \prod_j^w \theta_{i, S_{i j}} \]
  3. 部分配列がモチーフ $\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)] \]
  4. 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}\]
  5. E (Expectation) ステップ: $\theta$ を使って $q(z_i)$ を更新

  6. M (Maximization) ステップ: $q(z_i)$ を使って、尤度最大となるよう $\theta$ を更新

  7. EMステップを繰り返す

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

Gibbs sampling
AlignACE (Linux実行形式のみ配布), MotifSampler
  1. N 本の各配列のモチーフ位置の初期値をランダムに与える
  2. ランダムに配列をひとつ選び、それ以外の配列のモチーフを並べてPWMを計算
  3. 先に選んだ配列の上でPWMをスライドさせつつスコア (l-merがPWMから生成される確率など) を計算
  4. スコアに比例した確率で新しいモチーフ位置を決める
  5. 各配列のモチーフ位置が安定するまで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)
  1. 配列 N 本のデータセットに対し M 個の異なるアルゴリズムをそれぞれ R 回ずつ走らせて K 個の予測サイトを得る
  2. 配列毎・アルゴリズム毎に予測サイトを集めてグループ化する
  3. 配列毎に全アルゴリズムの予測スコアを足しあわせて 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

Pfam

PROSITE

ProDom

PRINTS

References