置換モデル

生物学では、 置換モデルは一連のシンボルが別の特性セットに変化するプロセスを表します。たとえば、分類学では、シーケンス内の各位置は、存在または非存在の種の特性に対応する場合があります。アルファベットは、不在の場合は「0」、存在の場合は「1」で構成されます。その場合、シーケンス00110は、たとえば、種に羽がないか卵を産み、毛皮があり、温血で、水中で呼吸できないことを意味します。別のシーケンス11010は、種が羽を持ち、卵を産み、毛皮がなく、温血で、水中で呼吸できないことを意味します。系統発生学では、最初にヌクレオチドまたはタンパク質配列のアライメントを取得し、次にアライメントの対応する位置の塩基またはアミノ酸を文字として取得することにより、シーケンスが取得されることがよくあります。これによって達成されるシーケンスは、AGGCGAGCTTAおよびGCCGTAGACGCのように見える場合があります。

置換モデルは、多くのことに使用されます。

進化系統樹を系統発生学または分類学で構築する。
シーケンスをシミュレートして、他の方法とアルゴリズムをテストします。

中立で独立した有限サイトモデル

これまで使用されてきたほとんどの代替モデルは、中立で独立した有限サイトモデルです。

ニュートラル選択は置換に対して機能しないため、制約はありません。あるサイトでの独立した変更は、別のサイトでの変更の確率に影響しません。有限サイト有限数のサイトがあるため、進化の過程で、1つのサイトを複数回変更できます。これは、たとえば、時間0および時間tで文字の値が0である場合、変更が発生していないか、1に戻って0に戻った、または1に変わって、 0に戻ってから1に戻ってから0に戻って、というように続きます。

分子時計と時間の単位

通常、系統樹の枝の長さは、サイトごとの予想される置換数として表されます。進化モデルが、先祖シーケンス内の各サイトが特定の子孫シーケンスに進化するまでに通常xの置換を受けることを示す場合、先祖と子孫は分岐長xで分離されていると見なされます。

枝の長さは地質年で測定される場合があります。たとえば、化石記録により、祖先の種と子孫の種の間の年数を決定できる場合があります。一部の種は他の種よりも速い速度で進化するため、枝の長さのこれらの2つの測定値は常に正比例するわけではありません。 1サイトあたりの年間置換の予想数は、ギリシャ文字のmu（μ）で示されることがよくあります。

1年あたりの予想される置換数μが、どの種の進化が調査されているかに関係なく一定である場合、モデルは厳密な分子時計を持つと言われます。厳密な分子時計の重要な意味は、祖先の種とその現代の子孫の間で予想される置換の数は、どの子孫種を調べるかとは無関係でなければならないということです。

厳密な分子時計の仮定は、特に進化の長期にわたって特に非現実的であることが多いことに注意してください。たとえば、げっ歯類は霊長類と遺伝的に非常に似ていますが、ゲノムの一部の領域での発散以来、推定時間内にはるかに多くの置換を受けています。これは、生成時間が短いこと、代謝率が高いこと、集団構造が増加していること、種分化率が高いこと、または体のサイズが小さいことによる可能性があります。分子時計の仮定の下でカンブリア紀の爆発のような古代の出来事を研究するとき、クラッドと系統発生データの間の不十分な同時発生がしばしば観察されます。進化の速度を可変にできるモデルに関するいくつかの研究がありました（例とを参照）。

系統発生における異なる進化系統間の分子時計の速度の変動を考慮に入れることができるモデルは、「厳格」とは反対に「緩和」と呼ばれます。このようなモデルでは、レートは先祖と子孫の間で相関するかしないかを想定でき、系統間のレート変動は多くの分布から引き出すことができますが、通常は指数分布と対数正規分布が適用されます。系統発生が少なくとも2つのパーティション（系統のセット）に分割され、それぞれに厳密な分子時計が異なるレートで適用される場合、「ローカル分子時計」と呼ばれる特別なケースがあります。

時間可逆モデルと定常モデル

多くの有用な代替モデルは時間可逆です。数学に関しては、他のすべてのパラメーター（2つのシーケンス間で予想されるサイトごとの置換数など）が一定である限り、モデルはどのシーケンスが祖先で、どのシーケンスが子孫であるかを気にしません。

実際の生物学的データの分析が実行されると、一般に祖先の種のシーケンスにはアクセスできず、現在の種にのみアクセスできます。ただし、モデルが時間可逆である場合、どの種が祖先種であったかは関係ありません。代わりに、系統樹は、いずれかの種を使用して根付かせたり、新しい知識に基づいて後で根付かしたり、根付かなかったりすることができます。これは、「特別な」種が存在せず、すべての種が最終的に同じ確率で互いに派生するためです。

モデルが時間可逆であるのは、それがプロパティを満たす場合のみです（表記については以下で説明します）

πiQij=πjQji{\ displaystyle \ pi _ {i} Q_ {ij} = \ pi _ {j} Q_ {ji}}

または、同等に、詳細なバランスプロパティ、

πiP（t）ij =πjP（t）ji {\ displaystyle \ pi _ {i} P（t）_ {ij} = \ pi _ {j} P（t）_ {ji}}

すべてのi 、 j 、およびtに対して 。

時間可逆性と定常性を混同しないでください。 Qが時間とともに変化しない場合、モデルは静止しています。以下の分析では、定常モデルを想定しています。

代用モデルの数学

定常、中立、独立、有限サイトモデル（一定の進化速度を仮定）には、2つのパラメーターπ 、ベース（またはキャラクター）周波数の平衡ベクトル 、およびレートマトリックス Qがあります。これは、1つのタイプのベース別のタイプのベースに変更します。 i ≠ jの要素Qij {\ displaystyle Q_ {ij}}は、ベースiがベースjに移動するレートです。 Q行列の対角線は、行の合計がゼロになるように選択されます。

Qii = −∑ {j∣j≠i} Qij、{\ displaystyle Q_ {ii} =-{\ sum _ {\ lbrace j \ mid j \ neq i \ rbrace} Q_ {ij}} \ ,,}

平衡行ベクトルπは、レート行列Qによって消滅する必要があります。

πQ= 0。{\ displaystyle \ pi \、Q = 0 \ ,.}

遷移行列関数は、分岐の長さ（時間単位、場合によっては置換）から条件付き確率の行列への関数です。 P（t）{\ displaystyle P（t）}と表記されます。 i番目の列とj番目の行のエントリPij（t）{\ displaystyle P_ {ij}（t）}は、時間tの後、条件jが与えられた位置にベースjがある確率です。時間0でその位置にベースiがあります。モデルが時間可逆である場合、2つのシーケンス間の分岐長の合計がわかっていれば、一方が他方の祖先でない場合でも、これを2つのシーケンス間で実行できます。

P ij （t）の漸近特性は、 P ij （0）= δijのようになります。ここで、 δijはクロネッカーデルタ関数です。つまり、シーケンスとそれ自体の間で基本構成に変更はありません。もう一方の極端な場合、limt→∞Pij（t）=πj、{\ displaystyle \ lim _ {t \ rightarrow \ infty} P_ {ij}（t）= \ pi _ {j} \ ,,}またはその他つまり、時間が無限になると、元のベースに関係なく、その位置にベースiが存在する場合にその位置にベースjが見つかる確率は、その位置にベースjがあるという平衡確率になります。また、πP（T）=π{\ displaystyleの\のPI P（T）= \ PI}すべてのtについてことになります。

遷移行列は、行列のべき乗を介してレート行列から計算できます。

P（t）= eQt = ∑n =0∞Qntnn！、{\ displaystyle P（t）= e ^ {Qt} = \ sum _ {n = 0} ^ {\ infty} Q ^ {n} {\ frac {t ^ {n}} {n！}} \ ,,}

ここで、 Q nは、 n番目のべき乗を与えるのに十分な回数だけ乗算された行列Qです。

Qは対角化可能である場合、行列指数を直接計算することができると共に、Qは = U -1ΛU Qの対角であるとします

Λ=（λ1…0⋮⋱⋮0…λ4）、{\ displaystyle \ Lambda = {\ begin {pmatrix} \ lambda _ {1}＆\ ldots＆0 \\\ vdots＆\ ddots＆\ vdots \\ 0＆\ ldots＆\ lambda _ {4} \ end {pmatrix}} \ ,,}

ここで、Λは対角行列で、{λi} {\ displaystyle \ lbrace \ lambda _ {i} \ rbrace}はQの固有値で、それぞれが多重度に従って繰り返されます。それから

P（t）= eQt = eU−1（Λt）U = U−1eΛtU、{\ displaystyle P（t）= e ^ {Qt} = e ^ {U ^ {-1}（\ Lambda t）U} = U ^ {-1} e ^ {\ Lambda t} \、U \ ,,}

対角行列Eの λTがによって与えられます

eΛt=（eλ1t…0⋮⋱⋮0…eλ4t）。{\ displaystyle e ^ {\ Lambda t} = {\ begin {pmatrix} e ^ {\ lambda _ {1} t}＆\ ldots＆0 \\\ vdots ＆\ ddots＆\ vdots \\ 0＆\ ldots＆e ^ {\ lambda _ {4} t} \ end {pmatrix}} \ ,.}

GTR：一般化された時間可逆

GTRは、最も一般的な中立、独立、有限サイト、時間可逆モデルです。 1986年にサイモン・タヴァレによって一般的な形で最初に記述されました。

ヌクレオチドのGTRパラメーターは、平衡基本周波数ベクトルπ→=（π1、π2、π3、π4）{\ displaystyle {\ vec {\ pi}} =（\ pi _ {1}、\ pi _ {2 }、\ pi _ {3}、\ pi _ {4}）}、各サイトで各ベースが発生する頻度、およびレートマトリックスを提供します

Q = 1 {-（x_ {1} + x_ {2} + x_ {3}）}＆x_ {1}＆x_ {2}＆x_ {3} \\ {\ pi _ {1} x_ {1} \ pi _ { 2}}＆{-（{\ pi _ {1} x_ {1} \ pi _ {2}} + x_ {4} + x_ {5}）}＆x_ {4}＆x_ {5} \\ { \ pi _ {1} x_ {2} \ over \ pi _ {3}}＆{\ pi _ {2} x_ {4} \ over \ pi _ {3}}＆{-（{\ pi _ {1 } x_ {2} \ over \ pi _ {3}} + {\ pi _ {2} x_ {4} \ over \ pi _ {3}} + x_ {6}）}＆x_ {6} \\ {\ pi _ {1} x_ {3} \ over \ pi _ {4}}＆{\ pi _ {2} x_ {5} \ over \ pi _ {4}}＆{\ pi _ {3} x_ {6 } \ over \ pi _ {4}}＆{-（{\ pi _ {1} x_ {3} \ over \ pi _ {4}} + {\ pi _ {2} x_ {5} \ over \ pi _ {4}} + {\ pi _ {3} x_ {6} \ over \ pi _ {4}}）} \ end {pmatrix}}}

モデルは時間可逆性でなければならず、長時間平衡ヌクレオチド（塩基）の周波数に近づく必要があるため、対角線より下の各速度は、2塩基の平衡比を掛けた対角線より上の逆数に等しくなります。そのため、ヌクレオチドGTRには6つの置換率パラメーターと4つの平衡基本周波数パラメーターが必要です。 4つの周波数パラメーターの合計は1である必要があるため、空き周波数パラメーターは3つしかありません。多くの場合、合計9個の無料パラメーターは、8個のパラメーターと単位時間あたりの置換の総数であるμ{\ displaystyle \ mu}にさらに削減されます。置換（μ{\ displaystyle \ mu} = 1）で時間を測定する場合、8つの空きパラメーターのみが残ります。

一般に、パラメータの数を計算するには、マトリックスの対角線より上のエントリの数をカウントします。つまり、サイトごとのn個の特性値n2-n2 {\ displaystyle {{n ^ {2} -n} \ over 2} }、次に平衡周波数にn-1を加算し、μ{\ displaystyle \ mu}が固定されているため1を減算します。あなたが得る

n2−n2 +（n−1）−1 = 12n2 + 12n−2。{\ displaystyle {{n ^ {2} -n} \ over 2} +（n-1）-1 = {1 \ over 2} n ^ {2} + {1 \ over 2} n-2。}

たとえば、アミノ酸配列（タンパク質を構成する「標準」アミノ酸は20個あります）には、208個のパラメーターがあります。ただし、ゲノムのコーディング領域を研究する場合、コドン置換モデルを使用することが一般的です（コドンは3塩基で、タンパク質の1つのアミノ酸をコードします）。 43 = 64 {\ displaystyle 4 ^ {3} = 64}のコドンがあり、2078の自由なパラメーターが得られますが、複数の塩基が異なるコドン間の遷移率がゼロであると想定される場合、20 ×19×32 + 63−1 = 632 {\ displaystyle {{20 \ times 19 \ times 3} \ over 2} + 63-1 = 632}パラメーター。

機械的モデルと経験的モデル

進化モデルの主な違いは、検討中のデータセットに対して毎回推定されるパラメーターの数と、大規模なデータセットで一度推定されるパラメーターの数です。機構モデルは、すべての置換を、できれば最尤法を使用して、分析されるすべてのデータセットに対して推定される多数のパラメーターの関数として記述します。これには、特定のデータセットの特性に合わせてモデルを調整できるという利点があります（たとえば、DNAの異なる構成バイアス）。特に相互に補正できる場合は、使用するパラメーターが多すぎると問題が発生する可能性があります。次に、データセットが小さすぎて、すべてのパラメーターを正確に推定するのに十分な情報が得られない場合がよくあります。

経験的モデルは、大規模なデータセットから多くのパラメーター（通常、レートマトリックスと文字頻度のすべてのエントリ、上記のGTRモデルを参照）を推定することで作成されます。これらのパラメーターは固定され、すべてのデータセットで再利用されます。これには、これらのパラメーターをより正確に推定できるという利点があります。通常、現在のデータセットのみから置換行列のすべてのエントリを推定することはできません。マイナス面として、推定されたパラメータは一般的すぎる可能性があり、特定のデータセットに十分に適合しない可能性があります。

大規模なゲノム配列決定により、非常に大量のDNAおよびタンパク質配列が生成されるため、任意の数のパラメーターを使用して経験的モデルを作成するのに十分なデータが利用可能です。上記の問題のため、2つのアプローチは、大規模なデータでほとんどのパラメーターを一度推定することで組み合わされることが多く、残りのいくつかのパラメーターは検討中のデータセットに合わせて調整されます。次のセクションでは、DNA、タンパク質、またはコドンベースのモデルで採用されているさまざまなアプローチの概要を説明します。

DNA置換のモデル

DNAモデルのより正式な説明については、主要記事「 DNA進化のモデル 」を参照してください。

DNA進化のモデルは、JukesとCantorによって1969年に最初に提案されました。すべての塩基について等しい遷移率と等しい平衡周波数を仮定しています。 1980年、木村は2つのパラメーターを持つモデルを導入しました。1つは遷移用、もう1つは転換率用で、1981年にFelsensteinは、置換率が標的ヌクレオチドの平衡周波数に対応する4パラメーターモデルを提案しました。長谷川、岸野、矢野（HKY）は、最後の2つのモデルを5パラメーターモデルに統合しました。 1990年代に、HKYに似たモデルが複数の研究者によって開発および改良されました。

DNA置換モデルの場合、主に機構モデル（上記の説明）が使用されます。推定するパラメーターの数が少ないため、これは実現可能ですが、DNAは多くの場合、生物や遺伝子の種類に応じて特定の目的（高速発現や安定性など）に高度に最適化されるため、これらの状況に合わせてモデルを調整する必要があります。

アミノ酸置換のモデル

多くの分析、特により長い進化距離の場合、進化はアミノ酸レベルでモデル化されます。すべてのDNA置換がエンコードされたアミノ酸を変更するわけではないため、ヌクレオチド塩基ではなくアミノ酸を見ると情報が失われます。しかし、いくつかの利点はアミノ酸情報の使用に有利です：DNAはアミノ酸よりも組成バイアスを示す傾向がはるかに高く、DNA内のすべての位置が同じ速度で進化するわけではありません（非同義の変異は同族のものよりも人口）、しかしおそらく最も重要なのは、それらの急速に進化する位置と限られたアルファベットサイズ（4つの可能な状態のみ）のために、DNAは逆置換によってはるかに苦しみ、より長い距離を正確に推定することを難しくしています。

DNAモデルとは異なり、アミノ酸モデルは伝統的に経験的モデルです。彼らは1970年代にDayhoffと同僚によって先駆者となり、少なくとも85％の同一性を持つタンパク質配列からの置換率を推定しました。これにより、サイトで複数の置換が観察される可能性が最小限になりました。推定レートマトリックスから、PAM250などの名前で知られている一連の置換確率マトリックスが導出されました。 Dayhoffモデルは、相同性検索結果の重要性を評価するために使用されましたが、系統解析にも使用されました。 Dayhoff PAMマトリックスは比較的少数のアラインメントに基づいていました（その時点で利用できなかったため）が、1990年代には、ほぼ同じ方法論を使用して新しいマトリックスが推定されましたが、そのとき利用可能な大きなタンパク質データベースに基づいています（後者は「JTT」行列として知られています）。