JP6423426B2 - 転写産物判定方法 - Google Patents

Description

本発明は、次世代シークエンシング（ＮＧＳ：ｎｅｘｔ ｇｅｎｅｒａｔｉｏｎ ｓｅｑｕｅｎｃｉｎｇ）リードに基づいて転写産物（例えば、ｍＲＮＡ）量の情報を提供することに関する。

次世代シークエンシング技術は、核酸サンプルを配列決定するときに大量のショートリードを作り出す。次世代シークエンシングに不可欠なステップは、ライブラリ調製または略してライブラリプレップ（ｌｉｂｒａｒｙ ｐｒｅｐ）である。このプロセスは、入力としてｍＲＮＡまたはｃＤＮＡを取り、それぞれがｍＲＮＡ分子の区分に対応する短いｃＤＮＡフラグメントのライブラリを作り出す。これらのフラグメントは、次にＮＧＳシーケンサーによって、通常はそれらの全体ではなくそれらの開始および／またはそれらの終結部において部分的に配列決定される。これは、ヌクレオチドの短い配列を生じ、この短い配列は、リードと称され、遺伝コードの核酸塩基を表すＡ、Ｃ、Ｇ、Ｔまたは０、１、２、３のような４つのＡＳＣＩＩ文字の一群の配列として、最も一般にはＮＧＳシーケンサーによって記憶される。元のサンプル中にどのｍＲＮＡ分子が存在したかを推測するために、リードがリファレンスゲノム上へマッピングされる。

次世代シークエンシングは、様々なゲノム・マッピング手順（特許文献１）または、例えば、配列リードをある生物バリアントへ関連付けるためにマッピングされたゲノムを用いることによるＤＮＡ同定方法（特許文献２）において利用されてきた。

特許文献３は、生物のトランスクリプトームのプロファイルを得るための方法を記載し、この方法は、シークエンシング・リードを得るために１つ以上のｃＤＮＡ分子を配列決定するステップと、それぞれのシークエンシング・リードをリファレンス配列へアライメントするステップとを備える。しかしながら、従来の方法には知られていない、短い配列リードを用いたトランスクリプトーム解析の根底にある主要な問題は、配列のずれが異なるアイソフォームのような、例えば、遺伝子の系統差、点変異、または１つのタンパク質のスプライスバリアントなどの転写バリアントの場合におけるアライメント・ステップである。通常、短い配列リードを１つの転写バリアントへ正しくアライメントすることは困難である。

配列リードに基づいて転写産物シークエンシング・データをアセンブルする最も一般的な方法は、「Ｃｕｆｆｌｉｎｋｓ」法である（非特許文献１）。Ｃｕｆｆｌｉｎｋｓは、ＲＮＡ−Ｓｅｑ実験において観測されたリードを「説明する（“ｅｘｐｌａｉｎ”）」転写産物の簡潔なセットを構築する。Ｃｕｆｆｌｉｎｋｓは、比較アセンブリ問題を二部グラフにおける最大マッチングの問題へ帰着させることによってこれを行う。基本に、Ｃｕｆｆｌｉｎｋｓは、リード・アライメントについて被覆関係を構築し、この関係に係る有向非巡回グラフ上において最小パス被覆を見出すことによって、ディルワースの定理の構成的な証明を実施する。この統計的手法を用いて、Ｃｕｆｆｌｉｎｋｓは、既知のリファレンス・アノテーションを用いるかあるいはリファレンスゲノムのみを用いた転写産物の非経験的なアセンブリ後に、サンプルに存在する転写産物アイソフォームの存在比を推定することができる。Ｃｕｆｆｌｉｎｋｓは、フラグメントのセットを所与として、ペアエンド・シークエンシング実験の統計モデルを用い、転写産物のセットにおける存在比について尤度を導出する。この尤度関数が一意的な最大値を有することを示すことができて、Ｃｕｆｆｌｉｎｋｓは、数値的最適化アルゴリズムを用いてこれを見出す。プログラムは、次に、転写産物に関して提案された存在比を所与として、これらの確率を乗算し、実験においてフラグメントを観測するであろう総合的な尤度を計算する。Ｃｕｆｆｌｉｎｋｓの統計モデルが線形であるために、尤度関数は、一意的な最大値を有し、数値的最適化アルゴリズムをもつＣｕｆｆｌｉｎｋｓによってこれが見出される。

非特許文献２は、フラグメント・バイアスを補正することによりＲＮＡ−Ｓｅｑ発現推定値を改善する方法に関する。非特許文献３は、自然個体群における転写産物量のクラスの混合モデリングを記載する。

従来の方法は、転写バリアントを正しく識別して他の転写産物との比較において正しい転写産物量または存在比を得ることができなかった。本明細書において比較により示されるように、Ｃｕｆｆｌｉｎｋｓ法でさえ、いくつかの実験では正しい転写産物量情報に到達することができなかった。

本発明の目標は、転写産物量のより正確な評価を可能にする改善された方法を提供することである。

米国特許出願公開第２０１３／１１０４１０（Ａ１）号 国際公開第２００９／０８５４１２（Ａ１）号 国際公開第２００９／０９１７９８（Ａ１）号

コール・トラップネル（Ｃｏｌｅ Ｔｒａｐｎｅｌｌ）、ブライアン・Ａ・ウィリアムズ（Ｂｒｉａｎ Ａ Ｗｉｌｌｉａｍｓ）、ジオ・ペルテア（Ｇｅｏ Ｐｅｒｔｅａ）、アリ・モルタザヴィ（Ａｌｉ Ｍｏｒｔａｚａｖｉ）、ゴードン・クワン（Ｇｏｒｄｏｎ Ｋｗａｎ）、マレイケ・Ｊ・ヴァン・バレン（Ｍａｒｉｊｋｅ Ｊ ｖａｎ Ｂａｒｅｎ）、スチーブン・Ｌ・ザルツバーグ（Ｓｔｅｖｅｎ Ｌ Ｓａｌｚｂｅｒｇ）、バーバラ・Ｊ・ウォルド（Ｂａｒｂａｒａ Ｊ Ｗｏｌｄ）およびリオル・パクター（Ｌｉｏｒ Ｐａｃｈｔｅｒ）転写産物アセンブリおよびＲＮＡ−Ｓｅｑによる定量化が細胞分化の間のアノテーションされていない転写産物およびアイソフォームのスイッチングを明らかにする（Ｔｒａｎｓｃｒｉｐｔ ａｓｓｅｍｂｌｙ ａｎｄ ｑｕａｎｔｉｆｉｃａｔｉｏｎ ｂｙ ＲＮＡ−Ｓｅｑ ｒｅｖｅａｌｓ ｕｎａｎｎｏｔａｔｅｄ ｔｒａｎｓｃｒｉｐｔｓ ａｎｄ ｉｓｏｆｏｒｍ ｓｗｉｔｃｈｉｎｇ ｄｕｒｉｎｇ ｃｅｌｌ ｄｉｆｆｅｒｅｎｔｉａｔｉｏｎ）Ｎａｔ Ｂｉｏｔｅｃｈｎｏｌ ２８（５）：５１１−５１５、２０１０年５月 アダム・ロバーツ（Ａｄａｍ Ｒｏｂｅｒｔｓ）、コール・トラップネル（Ｃｏｌｅ Ｔｒａｐｎｅｌｌ）、ジュリー・ ドナヒー（Ｊｕｌｉｅ Ｄｏｎａｇｈｅｙ）、ジョン・Ｌ・リン（Ｊｏｈｎ Ｌ Ｒｉｎｎ）およびリオル・パクター（Ｌｉｏｒ Ｐａｃｈｔｅｒ）フラグメント・バイアス補正によるｒｎａ−ｓｅｑ発現推定値の改善（Ｉｍｐｒｏｖｉｎｇ ｒｎａ−ｓｅｑ ｅｘｐｒｅｓｓｉｏｎ ｅｓｔｉｍａｔｅｓ ｂｙ ｃｏｒｒｅｃｔｉｎｇ ｆｏｒ ｆｒａｇｍｅｎｔ ｂｉａｓ）Ｇｅｎｏｍｅ Ｂｉｏｌ １２（３）：Ｒ２２ ２０１１年３月 ウェン−ピン（Ｗｅｎ−Ｐｉｎｇ）ら：ゲノム生物学（Ｇｅｎｏｍｅ Ｂｉｏｌｏｇｙ）８（６）（２００７）：Ｒ９８

本発明は、転写産物量を推定する方法を提供し、方法は、ａ）対象となる遺伝子座内の転写産物の潜在的な混合物から転写産物フラグメントシークエンシングデータを得るステップと、ｂ）前記フラグメントシークエンシングデータを対象となる前記遺伝子座の遺伝子座標へ割り当て、それによってフラグメント遺伝子座標カバレッジのデータセットを得るステップであって、遺伝子座標ごとの前記カバレッジは、結合されて（全カバレッジ・ヒストグラムまたは全体のヒストグラムとも呼ばれる）カバレッジ包絡曲線を形成する、ステップと、ｃ）前記混合物の転写産物の数をセットするステップと、ｄ）転写産物ｉごとにモデリングされた遺伝子カバレッジの確率分布関数を予めセットするステップであって、ｉは、転写産物のための数値識別子を示し、前記確率分布関数は、前記転写産物ｉの重み係数α_ｉと、少なくとも２つの確率サブ関数ｊの和との数学的積からなり、ｊは、確率サブ関数のための数値識別子を示し、各確率サブ関数ｊは、重み係数β_ｉ，ｊによって独立に重み付けされる、ステップと、ｅ）サム関数を得るために転写産物ごとの確率分布関数を加算するステップと、ｆ）サム関数を前記カバレッジ包絡曲線へフィッティングし、それによってフィットを向上させるために、α_ｉおよびβ_ｉ，ｊに関する値を最適化するステップと、ｇ）予めセットされた収束判定基準が満たされるまでステップｅ）およびｆ）を繰り返し、それによって収束判定基準が満たされた後に最適化されるような重み係数α_ｉによって与えられる、混合物の転写産物ごとの推定転写産物量を得るステップとを備える。

本発明は、前記方法を利用した、例えば、前記方法およびステップをコンピュータ上で行うかまたは補助するための機械コードを含む、コンピュータプログラム製品をさらに提供する。コンピュータプログラム製品は、任意の種類のメモリ装置上に提供できる。さらに提供されるのは、システム、例えば、本発明の方法のステップを行うことを補助するようにプログラムされたコンピュータ装置である。算出ステップは、通常、オペレータの補助なしに行われる。入力および設定ステップは、プログラムまたはシステムによって、例えば、ステップｄ）において確率サブ関数の数およびタイプに関する選択肢の提案を示唆することによって補助されてもよい。もちろん、プログラムまたはシステムは、オペレータからのさらなる入力なしにデフォルト・パラメータを用いて行われてもよい。

以下の詳細な記載および好ましい実施形態は、本発明のすべての態様に適用され、明示的に示された場合を除いて、制限なしに互いに組み合わせることができる。好ましい実施形態および態様は、特許請求の範囲において定義される。

発明の詳細な開示
本発明は、転写産物フラグメント配列のサンプルから転写産物量情報を得るために数値的方法を利用する。

この方法は、遺伝子カバレッジ情報を得るために、一般に転写産物フラグメント配列と呼ばれる（ＮＧＳ）リードをリファレンスゲノムのようなリファレンス配列へアライメントする（ステップｂ）。このために用いられる従来の統計学的ツールは、観測データの特質についてしばしば非現実的な仮定を行い、それゆえに転写産物濃度の推定値が不正確である。最も広く用いられるツールのいくつか、例えば、Ｃｕｆｆｌｉｎｋｓは、リードの分布が転写産物に沿って均一であると仮定し、これは現在のｍＲＮＡ−Ｓｅｑプロトコルと矛盾する。本発明は、転写産物に沿ったリード分布のバイアスと、転写産物量とを同時に知ることが可能な統計学的モデルを提供する。このために、転写産物のリードまたはフラグメント分布を混合関数によってモデリングして、この混合関数を転写産物量とともにフィッティング・ステップにおいてトレーニングする。フィッティング・ステップに用いる方法は、期待値最大化アルゴリズムを用いた最大尤度の枠組みなど、従来の最大化または最小化手順から推測できる。本発明のモデルにおけるリードの全確率分布は、混合関数を混合した関数なので、このモデルをＭｉｘ^２（Ｍｉｘｑｕａｒｅと読む）モデルと称する。以下に示すのは、Ｍｉｘ^２モデルが非常に多用途であり、随意的なパラメータの連結によって、データに本来備わった様々な構造へ調整できることである。特に、転写産物量を得るために用いる方法は、転写産物に係る確率分布に適しうる。実験において示すのは、Ｍｉｘ^２モデルが転写産物量に関してＣｕｆｆｌｉｎｋｓプログラムに用いる統計学的モデルよりかなり良好な推定値を達成することである。不正確な転写産物アノテーションから開始したときでも、Ｍｉｘ^２モデルは、正しいアノテーションをデータから知ることができ、先行技術よりはるかに優れた存在比の推定値を生み出す。フィッテング・ステップの間の優れた学習能力ゆえに、例えば、割り当てステップａ）の間またはステップｄ）における（例えば、ランダム）確率分布関数の選択の間に選択される開始パラメータが決定的に重要ではない。仮定する転写産物の数も異なってもよい。誤った転写産物アノテーションまたは転写産物数の仮定は、例えば、１つ以上の転写産物に関する統計分布関数を存在比ゼロへフィッティングすることによって、おそらく補正されるであろう。重み係数アルファは、確率としてモデリングでき、収束後に転写産物の存在比を表すことになろう。

本明細書において、「対象となる遺伝子座（“ｇｅｎｅｔｉｃ ｌｏｃｕｓ ｏｆ ｉｎｔｅｒｅｓｔ”）」は、本方法を染色体上の遺伝子の１つの連続的な配列伸長鎖には限定しない。「対象となる遺伝子座」は、一般に、遺伝子配列の１つ以上の区分を指す。遺伝子座は、ゲノム座標と関連付けることができて、配列リード（転写産物フラグメント配列）に関する位置情報を提供する。ゲノム「位置（“ｐｏｓｉｔｉｏｎ”）」または「座標（“ｃｏｏｒｄｉｎａｔｅ”）」は、本明細書では、リファレンス配列の開始からある距離を隔てたリファレンス配列上で数値により識別されるヌクレオチドを指すために用いる。遺伝子座標は、遺伝子座標が転写産物と適合すれば、ゲノム座標または転写産物座標で表すことができ、遺伝子座標がいずれかの転写産物エクソン内に位置すれば、これが当てはまる。適合するゲノム座標は、転写産物開始からの相対距離を算出して、ゲノム座標に先行するイントロンの長さを減算することによって転写産物座標へ変換される。転写産物座標は、このプロセスを逆戻りすることによってゲノム座標へ変換できる。１つの転写産物の座標上で定義された確率分布が確率サブ関数の縮小、伸張およびシフトによって別の転写産物の転写産物座標上の確率分布へ変換されてもよい。

核酸塩基タイプ（例えば、Ａ、Ｔ／Ｕ、Ｇ、Ｃ）がゲノム座標または位置と関連付けられても、関連付けられなくてもよい。通常、各転写産物に対して１つのタイプの核酸塩基が各遺伝子座標と関連付けられる。しかしながら、本発明は、点変異を識別するために用いることができるので、遺伝子座標ごとに、異なる重複転写産物の核酸塩基構成が異なってもよく、すなわち、特にサンプルが異なる生物または異なる対立遺伝子の核酸分子を含んだ場合には、異なる転写産物間で１つ以上の核酸塩基が異なってもよい。対象となる遺伝子座内の（ａｎｄ ａ ｌｏｃｕｓ ｏｆ ｉｎｔｅｒｅｓｔ）転写産物、それらのフラグメント間で生じうるミスマッチの場合においても共通の遺伝子座標を配分することができる。フラグメント配列を遺伝子座標へ割り当てるステップｂ）が遺伝子座標を提供するこれらの配列間の配列比較またはアライメントを備えてもよい。ヌクレオチドの比較を備える（例えば、前述のオープンソースＣｕｆｆｌｉｎｋｓ法による）配列比較は、当技術分野においてよく知られている。本発明の方法は、対象となる所与の遺伝子座上にずれた配列をもつ異なる系統または生物のデータに対して用いることができ、或いは、異なるスプライスバリアント、すなわち、エクソン配列の異なる組み合わせによって区別できる異なる転写産物を識別するために本発明の方法を用いてもよい。このように、割り当てステップａ）の間のミスマッチは、許容することも却下することもできる。かかるミスマッチは、好ましくは、１００塩基当たり多くとも１、２、３、４、５、６、７、８、９、１０または１１個のミスマッチである。

本発明は、転写産物の混合物における個々の転写産物の（例えば、確率に相当する相対量のような）存在比をモデリングできる。転写産物フラグメントシークエンシングデータは、少なくとも２、好ましくは３、４、５、６、７、８、９、１０、あるいは少なくともまたは多くとも１５、少なくともまたは多くとも２０、少なくとも２５、少なくともまたは多くとも３０、少なくともまたは多くとも４０個の転写産物配列を備えてもよい。転写産物の数は、例えば、核酸を例えば増幅または除去することにより核酸サンプルに対して行う選択ステップにおいて選択できる。このように除去または増幅される核酸が共通の配列伸長鎖、例えば、選択したオリゴヌクレオチド・プローブに対応する配列、アンカー配列またはプライマー配列を含んでもよい。１つ以上のゲノム遺伝子座がこのような方法で選択されてもよい。本発明の特に好ましい実施形態では、１つの遺伝子の複数の転写バリアントが検討される。とは言え一般に、「転写産物」は、本明細書において任意の遺伝子または遺伝子の組み合わせおよびその任意のバリアントの任意の核酸またはその配列、特にそのｍＲＮＡまたはｃＤＮＡ配列バリアントを指してもよい。

好ましくは、前記転写産物フラグメントシークエンシングデータの転写産物フラグメント配列は、５から８００ヌクレオチド、好ましくは６から６００ヌクレオチド、７から４００ヌクレオチド、８から２００ヌクレオチドまたは９から１５０ヌクレオチド、なおさらに好ましくは１０から１００ヌクレオチド、特に好ましくは１２から７０ヌクレオチドの長さを有する。

転写産物フラグメント配列の数は、少なくともまたは多くとも１００、少なくともまたは多くとも５００、少なくともまたは多くとも１０００、少なくともまたは多くとも５０００、少なくともまたは多くとも１００００であってもよい。好ましくは、転写産物ごとに、少なくともまたは多くとも１０、或いは少なくともまたは多くとも２０、少なくともまたは多くとも５０であってもよい。組み合わせてまたは代わりとして、転写産物フラグメント配列の数は、多くとも４０００００、多くとも３０００００、多くとも２０００００、多くとも１０００００または多くとも５００００であってもよい。

少なくとも１つ以上、例えば、すべての転写産物の転写産物長は、例えば、１００から１００００００、好ましくは１０００から１０００００ヌクレオチドまたは２０００から１００００ヌクレオチドであってもよい。

好ましい実施形態において、対象となる遺伝子座は、１つ以上の遺伝子または遺伝因子の、例えば、転写産物配列をコードする、１つ以上のアイソフォームを備え、好ましくは１つの遺伝子または遺伝因子の少なくとも２、３、４個またはそれ以上のプライスバリアントを備える。対象となる遺伝子座は、別の遺伝子または遺伝因子の１つ以上の他のプライスバリアントを備えてもよい。スプライスバリアントに加えてまたは代わりに、対象となる遺伝子座は、異なる対立遺伝子を備えてもよい。好ましい実施形態において、遺伝子または遺伝因子は、タンパク質（例えば、ｍＲＮＡ）をコードするが、タンパク質をコードしない転写産物、例えば、ｍｉｃｒｏＲＮＡ、ｓｎｏＲＮＡもしくはｒＲＮＡ、ならびにそれらの前駆体、特にｐｒｅ−ｍｉｃｒｏＲＮＡまたはｐｒｅ−ｒＲＮＡを含む、調節または触媒ＲＮＡもコードする。

本明細書では、「遺伝子」および「遺伝因子」は、１つ以上の転写産物を形成するために転写される配列をもつ遺伝子ヌクレオチドに関する。

本明細書では、「アイソフォーム」は、転写産物の特定のバリアントに関係して用いる。１つの「遺伝子」または「遺伝因子」の転写産物は、例えば、スプライスバリアントの場合、異なってもよく、従って異なるアイソフォームを生じる。他のアイソフォーム変化は、異なる対立遺伝子、または、例えば異なる微生物もしくは系統の混合における転写産物物質の異なるソースによって生じうる。

転写産物の数をセットするステップは、予めアノテーションされている配列データを対象となる遺伝子座から得ることと、転写産物の数を、対象となる前記遺伝子座から予想される、異なる転写産物配列としてカウントするスプライスバリアントを含む、異なる転写産物配列の少なくとも数にセットすることとを備えてもよい。前述のように、誤った転写産物は、サム関数のフィッティングの間に除去されてもよく、ゼロへ収束する１つ以上の重み係数アルファを結果として生じるため、転写産物のセット数は、（厳密に分かっても分からなくてもよい）転写産物の実際の数を超えてもよい。通常、対象となる遺伝子座は、例えば、核酸選択ステップのうちのいずれか１つのｃＤＮＡ生成ステップから知られる。アノテーションされている遺伝子データを用いて、転写産物の量に係る開始数ｉ_ｍａｘに到達することが可能である。転写産物ごとに、ステップｄ）において確率分布関数をセットすることによって存在比をモデリングすることになる。この関数は、重み係数アルファを含み、重み係数アルファは、和および確率分布フィッティングの収束、反復プロセス後に混合物における転写産物の存在比に対応する。転写産物の混合物は、本発明のモデルの第１の混合物である。各アルファは、フィッティング・プロセスの間に別々に修正できる。従って、転写産物の数のセッティングは、粗い推定値であってもよい。例えば、この数は、２、３、４、５、６、７、８、９、１０、１１、１２、１３、１４、１５であってもよく、１６、２０、３０、４０またはそれ以上のうちの少なくともいずれか１つであってもよい。

第２の混合物は、各確率分布関数に数学的に含まれる。確率分布関数は、２つ以上（例えば、２、３、４、５、６、７、８、９、１０またはそれ以上）の確率サブ関数または「ブロック」からなる。転写産物の長さ（転写産物座標範囲）にわたって非対称的なリード分布のモデリングを可能にするためには、少なくとも２つの確率サブ関数を有することが不可欠である。転写産物当たりの確率サブ関数の総量は、ｊ_ｍａｘである。ｊ_ｍａｘは、通常、およそ４〜８である。あまりに多い量の確率サブ関数は、過剰なトレーニングおよび計算効率の低下に繋がりかねない。

確率分布関数は、確率サブ関数の和である。注目すべきは、係数アルファが確率分布関数全体を重み付けし、従って、各基礎をなす確率サブ関数を均等に重み付けすることである。確率分布関数内では、確率サブ関数が重み係数β_ｉ，ｊ（ベータ_ｉ，ｊ）によって別々に重み付けされる。注目すべきは、転写産物ごとに、すなわち確率分布関数ごとに、ベータが独立係数または従属係数（「連結される」）であってもよいことである。独立係数の場合、各ベータは、フィッティング・ステップの間に個別に修正できる。さらに可能なのは、例えば、第１のベータ、または一般に任意のベータｊが各確率分布関数に対して同じになるように、ベータを異なる確率分布関数間で連結することである（例えば、ベータ_１，ｊ＝ベータ_２，ｊ）。これによって、モデリングされる重み係数ベータの数がｉ_ｍａｘ×ｊ_ｍａｘではなく、単にｊ_ｍａｘへ減少して、フィッティング・プロセスが簡略になり、計算リソースの使用量が減少する。この簡略化は、転写産物ごとに同様のフラグメント・カバレッジ分布が予想できる場合、例えば−２つ以上の、例えば、すべての転写産物（「連結された」転写産物）に対して−所定の配列部分では他の部分より常に大きい存在比が予想される場合に最もよく機能する。この仮定は、通常、およそ同じかまたは同様の長さの転写産物に関して、例えば、連結された別の転写産物の０．３から３．２、好ましくは０．５から２．１の（核酸塩基単位の）長さを有する転写産物に関して真である。もちろん、以下でより明確になるように、確率サブ関数の他のパラメータも、異なる転写産物間で同様に、または（例えば、ある１つの転写産物に関する確率分布関数全体を遺伝子座標方向に横へシフトさせるときに）所与の転写産物の確率サブ関数間でやはり同様に連結できる。また、ステップｄ）において、転写産物の確率サブ関数を、互いに（通常、確率分布関数の全体に対する部分の確率または確率寄与としてモデリングされる）およそ同じ最大高さまたは最大値に正規化することが可能である。

パラメータを連結するためのさらなる選択肢によれば、すべての転写産物、すなわち、各確率分布関数にわたってではなく、転写産物、すなわち、各確率分布関数の一群において、１つ以上のパラメータ、例えば、１つ以上のベータのみを連結することも可能である。かかる群は、確率分布関数のサム関数の同様の予想形状によって定義されてもよい。例えば、同様の長さおよび／または同様のＧＣ含量によって群が定義されてもよい。同様の長さとは、例えば、群のすべてのメンバーがすべてのメンバーの平均サイズの＋／−５０％以内、好ましくは＋／−３５％以内、より好ましくは＋／−３０％以内の長さを有するときである。同様のＧＣ含量とは、例えば、すべてのメンバーの平均ＧＣ含量の＋／−１０％以内、好ましくは＋／−５％以内であってもよい。

「およそ（“ａｂｏｕｔ”）」は、本明細書では、所定の値と同じ値もしくは所定の値と±１０％異なる値を指す。

所定の重み係数を用いた転写産物ｉに関する確率サブ関数の和は、転写産物（ｔｒａｎｓｃｒｉｐ）Ｉに関する確率分布関数を構成する。確率サブ関数としては、これらの確率サブ関数の和が確率関数を形成できることを必要条件として、任意の数学関数を選択できる。確率サブ関数は、フィッティングおよび最適化が可能な計算モデルのための基礎をなすため、少ないパラメータを用いた簡略な関数が好ましい。確率サブ関数は、遺伝子座標に応じて関数の形状を確定する、例えば、１、２、３、４、５、６、７、８、９または１０個の関数パラメータを備えるか、またはそれらから成ってもよい。確率サブ関数ｊは、好ましくは遺伝子座標ごとに正値で構成される。好ましくは、確率サブ関数は、非周期関数であるか、および／または特に好ましくは、確率サブ関数は、密度関数または確率関数であるが、（所定の転写産物に関する）確率分布関数と呼ばれる、確率サブ関数の和と区別するために本明細書では確率サブ関数または「ブロック」と呼ばれる。確率サブ関数は、通常、ある遺伝子座標において最大値を含み、前記最大値から離れた正または負の遺伝子座標では着実に減少する。好ましくは、確率サブ関数は、唯一の最大値を有する。確率サブ関数ｊの例は、ガウス関数、正方形関数、三角形関数、好ましくはガウス関数から選択される。

遺伝子座標またはその範囲は、好ましくは各転写産物に対して同じである。この目的のために、転写産物に関する確率分布関数を遺伝子座標方向に縮小もしくは拡張、またはシフトさせることができる。もちろん、リファレンス配列、例えば、ゲノムにおける正しいソース座標に再び到達するために、かかる変換を反転させることができる。しかしながら、かかる変換は、存在比モデリング・アルゴリズムにとって有益であろう。

この縮小、拡張およびシフトは、正しい（または尤度がより高い）転写産物長の知識を得るため、または高めるために用いることができる。プリファーメントにおいて、転写産物長の推定値は、このステップによって得られ、フラグメントの数（例えば、量または濃度）を確定するためにこれを用いる。フラグメント数は、関係：ＦＰＫＭ＝フラグメントの数／転写産物長（×重み係数）に基づいて、ＦＰＫＭ（１００万リード当たり、１０００塩基対当たりのフラグメント：ｆｒａｇｍｅｎｔｓ ｐｅｒ ｋｉｌｏ ｂａｓｅ ｐａｉｒ ｐｅｒ ｍｉｌｌｉｏｎ ｒｅａｄｓ）から推定できる。従って、フラグメントの数は、ＦＰＫＭ×転写産物長と相関する。

転写産物長の推定値は、ステップｅ）の収束後のサム関数を随意的に縮小もしくは拡張および／またはシフトさせた確率分布関数に従って、転写産物開始および終結位置をセットすることによって得ることができる。定義された開始および終結位置は有さないが、それに関する確率値のみを有しうる、複数の確率分布関数からなるサム関数の特質に起因して、開始および終結は、例えば、（開始に関しては）遺伝子座標までの部分面積、または（終結に関しては）遺伝子座標から始まる部分面積が、（サム関数または最初もしくは最後の確率分布関数の）曲線下の全面積のある割合である前記遺伝子座標として推定できる。この割合は、１％と１０％との間とすることができる。当業者は、用いる確率（ｐｒｏｂａｂｌｙ）分布関数の形状に応じて、かかる面積カットオフを確定するのに適した値を容易に試験できる。これは、既知の開始および終結位置をもつモデル核酸を用い、他の開始および終結位置をもつ本発明の関数を設定し、記載するように、本発明のアルゴリズムに改善された開始および終結位置をトレーニングさせて容易に試験できる。カットオフ値に用いる値の例は、０．５％、１％、２％、３％、４％、５％、６％、７％、８％、９％または１０％である。

さらなる選択肢によれば、フラグメント長分布または平均フラグメント長は、ステップｅ）の収束後に本発明のモデル関数を用いて算出される。

特に、１つ以上の転写産物、好ましくは各転写産物のための遺伝子座標（およびその確率分布関数）は、対象とならない遺伝子領域を削除するように随意的に変換した、ゲノムにおけるヌクレオチド位置に対応してもよく、好ましくは−リードまたは配列フラグメントが何もアライメントしないイントロンのような−前記対象とならない遺伝子領域は、前記転写産物フラグメントシークエンシングデータによるカバレッジを含まない。このように、連続的な確率分布関数を間歇的なゼロ分布を伴わずにモデリングできる。かかる対象とならない領域は、個々の確率分布関数から除去でき、加えてまたは代わりに（各反復ステップｇ）の間に様々な確率分布関数のサム関数をフィッティングする対象となる）カバレッジ包絡曲線からも除去できる。包絡曲線は、結合された配列データ・カバレッジのヒストグラムを表す。

このように、好ましくは、本発明の方法は、イントロンをもつ遺伝子座標位置を前記カバレッジ包絡曲線から除去することを備えたステップｂ２）をさらに備える。これは、もちろん、元のフラグメントシークエンシングデータまたは対象となる遺伝子座の遺伝子もしくはゲノム座標に戻って確率分布関数を再び参照する場合を除いて、ステップｆ）および好ましくはすべての他のステップにおけるカバレッジ包絡曲線の使用を置き換える、修正されたカバレッジ包絡曲線を得ることにつながる。

存在比をモデリングするために、遺伝子座標カバレッジがリードの配列情報全体を用いる必要はない。さらに、リード配列の限られた数のヌクレオチドのみ、例えば、単に最初の１つまたは開始点の情報を用いれば十分である。フラグメント遺伝子座標カバレッジは、遺伝子座標に割り当てられたフラグメント配列ごとに少なくとも１つのヌクレオチドのカウントを含むことができ、好ましくは、少なくとも１つのヌクレオチドは、フラグメント開始点またはフラグメント配列全体を備える。「備える（“ｃｏｍｐｒｉｓｅｓ”）」は、本明細書では、含む（“ｃｏｎｔａｉｎｉｎｇ”）におけるようにさらなるメンバーを許容する、開いた定義として理解するものとする。他方、「成る（”ｃｏｎｓｉｓｔｉｎｇ”）」は、成るの定義の特徴のさらなる要素を伴わない閉じた定義と見做される。このように、「備える」は、より広い定義であり、「成る」の定義を包含する。「備える」の語を用いた本明細書における任意の定義は、本発明の特別の実施形態では成るの制限とともに読まれてもよい。

確率分布関数の成分である、２つ以上の確率サブ関数は、（転写産物の）前記確率分布関数のゲノム座標範囲のうちに好ましくは均等に分布する。一般に、好ましくは、転写産物に関する確率サブ関数は、それぞれ異なる遺伝子座標に極大を備える。本発明の特に好ましい実施形態において、転写産物ごとの確率サブ関数の極大は、随意的に（イントロンの除去後のように、修正された遺伝子カバレッジを提供しうる）上述のような変換後に、少なくともおよそ、ステップａ）において転写産物に割り当てられた遺伝子座標の最初と最後のヌクレオチドとの間の差に１／ｊ_ｍａｘを乗じたものだけ分離される。このように、１つの転写産物（それゆえにカバレッジをモデリングする確率密度関数）によってカバーされる遺伝子座標に沿って複数の確率サブ関数が均等に分布することができる。もちろん、最適な均等分布からずれることもありうる。従って、確率サブ関数ごとに、遺伝子座標上の最大位置は、均等分布の最大値から、転写産物長またはｊ_ｍａｘで除した転写産物長の、例えば、０％〜５０％、好ましくは１０％〜４０％ずれてもよい。

このように、ステップｄ）では正値を用いて転写産物の全長をカバーするために、遺伝子座標において転写産物に関する確率サブ関数を配置するか、またはシフトさせることができる。

好ましい実施形態において、本方法は、少なくとも１つの転写産物、好ましくはｍＲＮＡまたはｃＤＮＡの配列リードを確定するステップを備え、前記リードは、前記転写産物フラグメントシークエンシングデータを提供するために前記転写産物のフラグメントの配列を備える。確定ステップは、当該技術分野で知られた任意の方法、例えば、ＰＣＲシークエンシングによって行うことができる。かかる方法は、マクサム−ギルバート・シークエンシング、チェーンターミネーション法、ショットガン・シークエンシング、ブリッジＰＣＲ、大規模並列処理特徴シークエンシング（ＭＰＳＳ：ｍａｓｓｉｖｅｌｙ ｐａｒａｌｌｅｌ ｓｉｇｎａｔｕｒｅ ｓｅｑｕｅｎｃｉｎｇ）、ポロニー・シークエンシング、ピロシークエンシング、イルミナ（ソロクサ）シークエンシング、ＳＯＬｉＤシークエンシング、イオン半導体シークエンシング、ＤＮＡナノボール・シークエンシング、ヘリスコープ一分子シークエンシング、一分子リアルタイム（ＳＭＲＴ：ｓｉｎｇｌｅ ｍｏｌｅｃｕｌｅ ｒｅａｌ ｔｉｍｅ）シークエンシング、ナノポアＤＮＡシークエンシング、ハイブリダイゼーションによるシークエンシング、質量分光法を用いたシークエンシング、マイクロ流体サンガー・シークエンシング、顕微鏡法ベース技術、ＲＮＡＰシークエンシング、インビトロウイルス（ｉｎ ｖｉｔｒｏ ｖｉｒｕｓ）高スループット・シークエンシングを含む。

本発明の方法は、代わりにまたは加えて、例えば、予め用意された配列データを用いる場合に、本発明の方法を行うことが可能な演算装置へ転写産物フラグメントシークエンシングデータを入力するステップを備えてもよい。前記演算装置は、次に、入力したシークエンシング・データを用いて存在比情報に到達するために本発明の方法を行うことになる。

存在比情報は、例えば、以下にさらに記載されるように、ビデオスクリーン、プリンタまたはコンピュータ可読媒体などの出力装置上に表示できる。

ステップｆ）におけるフィッティングは、期待値最大化アルゴリズムによって行うことができる。期待値最大化アルゴリズムは、デンプスターら、１９７７、に記載されている。もちろん、任意の他の最大化または最小化アルゴリズムを適用して本発明のサム関数が包絡曲線にフィッティングするようにできる。基本的に、重みアルファｉおよびベータｉ，ｊ−ならびに随意的に、確率サブ関数を定義する任意のさらなるパラメータは、サム関数（モデル）とカバレッジ包絡曲線（フラグメント配列に基づく実際のデータ）との間の差を、例えば、確率サブ関数のパラメータを変化させることによって最小化するように修正される。特に、算出ステップｅ）およびｆ）の繰り返しによる数回の反復後に、最小値を見出すべくかかる最小化または変更を行うために多くのアルゴリズムが存在する。

確率サブ関数のかかる追加のパラメータは、例えば、半値全幅値、または、特にガウス関数の場合には、標準偏差値（シグマｉ，ｊ）である。これらは、関数幅のパラメータである。好ましくは、転写産物ｉに関する確率サブ関数ごとに、確率サブ関数の最大値とは別のいずれかの遺伝子座標における関数幅（例えば、好ましくは半値全幅値または標準偏差）は、およそ同一であり、好ましくは同一である。代わりに、幅対最大値の比が転写産物ｉに関する確率サブ関数ごとにおよそ同一であり、好ましくは同一である。これは、転写産物の異なる確率分布関数間のベータについて上述したのと同様に関数幅を「連結する（“ｔｉｅｓ”）」。

一般に、好ましい実施形態において、フィッティング・ステップｆ）の間にひとつの転写産物の確率サブ関数の群内および／または複数の転写産物に対して同じ識別子ｊをもつ確率サブ関数の群内で確率サブ関数の少なくとも１、２、３、４、５、または６個のパラメータを連結し、例えば、互いに組み合わせて修正する。

最後のステップは、フィッティングの収束である。収束判定基準は、オペレータによってセットできる。これは、パラメータα_ｉもしくはβ_ｉ，ｊのうちのいずれか１つまたはそれらの組み合わせ、例えば、すべてのαｓもしくはすべてのβｓ、またはすべてのαｓおよびβｓの最大値の調整、あるいは、例えば、実施例の式（１５）のような、確率分布関数における最大値の調整とすることができる。例えば、収束判定基準は、式（２１）による対数尤度の増加が０．５未満である場合とすることができる。所望の品質を達成するために、オペレータによって任意の収束判定基準を選択できる。

プリファーメントにおいて、本発明の方法を他のバイアス低減方法とさらに組み合わせる。例えば、ステップｅ）で得た推定転写産物量をバイアス係数によってさらに重み付けし、それによって重み付けした推定転写産物量を得る。このバイアス係数は、ある転写産物が湿式化学的な転写産物またはフラグメント生成の間に増加または減少することを考慮に入れることができる。いくつかのフラグメント、例えば、Ｃで開始するフラグメントは、過剰発現することがあり、好ましくは１未満のバイアス係数によって重み付けし、それによって重み付けした推定転写産物量をステップｅ）の推定転写産物量に比べて減少させる。Ｃで開始するフラグメントが過剰発現する化学的な理由は、プライマーの動態にある、すなわち、ＧまたはＴで開始するテンプレートと比較して、Ｃで開始するテンプレートを用いるとプライマーがより良好にアニールする。それゆえに、Ｇおよび／またはＴで開始するフラグメントは、Ｃで開始するフラグメントのバイアス係数より大きいバイアス係数によって重み付けするとよい。一般に、当業者は、任意の多核酸特性に依存しうる系統的なバイアスが存在する任意のパラメータに関してかかるバイアス係数を取得し、適合したバイアスを含めてステップ３）の結果を相応に重み付けする（ａｎ ｗｅｉｇｈ ｔｈｅ ｒｅｓｕｌｔ ｏｆ）ことができる。

本発明は、さらに、本発明による方法をコンピュータ上で行うため、または本発明による方法をコンピュータによってサポートするためのコンピュータプログラム製品を備えたコンピュータ可読メモリ装置に関し、特に、ステップａ）、ｂ）、ｂ２）、ｃ）、ｄ）、ｅ）、ｆ）およびｇ）のいずれか１つをコンピュータ上で行うことができる。本発明に関する任意の方法またはステップをコンピュータに実装した方法として行うことができる。配列リードを確定する通常は湿式化学的なステップも、例えば、自動化または半自動化配列リーダを制御してそこからデータを得るためにコンピュータによって補助されてもよい。コンピュータプログラム製品またはメモリ装置にはサンプルからショートリードを得るリード生成コンポーネント、例えば、シーケンサー、好ましくは、コンピュータ・コンポーネントを備えるシーケンサーがさらに設けられてもよい。例えば、コンピュータ可読媒体は、磁気記憶装置（例えば、ハードディスク、フロッピディスク、磁気ストリップ、．．．）、光ディスク（例えば、コンパクトディスク（ＣＤ：ｃｏｍｐａｃｔ ｄｉｓｋ）、デジタル多用途ディスク（ＤＶＤ：ｄｉｇｉｔａｌ ｖｅｒｓａｔｉｌｅ ｄｉｓｋ）、スマートカードならびにフラッシュメモリ装置（例えば、カード、スティック、キーデバイス、．．．）を含みうるが、これらには限定されない。

本発明を以下の図面および例によってさらに説明するが、本発明のこれらの実施形態に限定されることはなく、各要素を本発明の任意の他の実施形態と組み合わせることができる。特に、本発明の方法におけるステップを記載または定義するために、以下に示す式のいずれか１つを別々に、個別に、または組み合わせて用いることができる。

ｃＤＮＡ生成、フラグメンテーション、シークエンシング、配列リードのリファレンス配列へのマッピング、およびフラグメント配列のマッピング・データセットの解析のステップに従うＮＧＳワークフローを示す。 ２つの転写産物の混合カバレッジ・ヒストグラムを示す。 参考図：Ｃｕｆｆｌｉｎｋｓモデルによって仮定されたフラグメント開始点の分布を示す。 参考図：Ｃｕｆｆｌｉｎｋｓモデルによって仮定されたカバレッジを示す。 Ｕｑｃｒｑ遺伝子上の転写産物リードのカバレッジを示す。 遺伝子座標変換を示す。平均確率ｐ_{（ａｖｇ）}（ｒ）のスケーリングおよびシフト・バージョンによる全確率ｐ_{（ｔｏｔａｌ）}（ｒ）の分解；全存在量包絡線（上部）、転写産物１に関する存在量（中央）、転写産物１に関して修正した遺伝子カバレッジを示すためのジャンクションの除去（シフト）およびスケーリングによる転写産物１に関する初期遺伝子カバレッジの変換（下部）を示す。 遺伝子座標上の存在比が等しい（それぞれが全体の１／３）３つのモデル転写産物のエクソン構造およびジャンクションを示す。 図７の転写産物に関するジャンクションの除去（シフト）後の転写産物座標における開始点確率分布ｐ（ｓ（ｒ）｜ｔ＝ｉ）を示す。 図７の転写産物に関するゲノム座標における開始点確率分布ｐ（ｓ（ｒ）｜ｔ＝ｉ）を示す。 ｉ＝１、２、３に対して存在比α_ｉ＝１／３が等しい図７および８の転写産物に関する開始点の全確率分布を示す。 図７の転写産物の拡張した開始および終結点を示す。モデルは、スケールおよびシフト・パラメータのラムダｉおよびニューｉを用いて正しい開始および終結点分布を推定する。 シグマの連結を示す：シグマ、確率サブ関数における随意的なさらなる標準偏差（幅）パラメータを転写産物ｉに関する各確率サブ関数内で一定に保ち、示されるのは、サム関数（下部）をもつ２つの転写産物１（中央）および２（上部）である。破線は、欠落した確率サブ関数であり、欠落していなければ１つの転写産物の遺伝子カバレッジ範囲にわたって均等に分布した確率サブ関数（ここではガウス関数）からなるサム関数をモデリングするための確率サブ関数である。 図１０における全フラグメント開始点分布からサンプリングした開始点のヒストグラムを示す。 サンプリングしたフラグメント長のヒストグラムを示す。 混合ガウシアンの成分（確率サブ関数）（点線）、実線：転写産物に関する確率分布関数を示す。 １ｔｉｅ−Ｍｉｘ^２モデルに対する重み係数アルファの収束を示す。 １ｔｉｅ−Ｍｉｘ^２モデルに対する重み係数ベータの収束を示す。 １ｔｉｅ−Ｍｉｘ^２モデルに対する最終的な重み係数ベータを示す。 １ｔｉｅ−Ｍｉｘ^２モデルに対する転写産物２の最終的な確率分布関数（ｐｄｆ：ｐｒｏｂａｂｉｌｉｔｙ ｄｉｓｔｒｉｂｕｔｉｏｎ ｆｕｎｃｔｉｏｎ）［ｐ（ｓ（ｒ）｜ｔ＝２）］を示す。 １ｔｉｅ−Ｍｉｘ^２モデルに対するサム関数［ｐ_{（ｔｏｔａｌ）}（ｓ（ｒ））］の収束を示す。 アルファ_（３）＝０．２を用いた、収束後の１ｔｉｅ−Ｍｉｘ^２モデル、Ｃｕｆｆｌｉｎｋｓモデルに対する推定および真のアルファ_ｉ（存在比、全体のうちの比率）を示す。 アルファ_（３）＝０．４を用いた、収束後の１ｔｉｅ−Ｍｉｘ^２モデル、Ｃｕｆｆｌｉｎｋｓモデルに対する推定および真のアルファ_ｉ（存在比、全体のうちの比率）を示す。 アルファ_（３）＝０．６を用いた、収束後の１ｔｉｅ−Ｍｉｘ^２モデル、Ｃｕｆｆｌｉｎｋｓモデルに対する推定および真のアルファ_ｉ（存在比、全体のうちの比率）を示す。 アルファ_（３）＝０．８を用いた、収束後の１ｔｉｅ−Ｍｉｘ^２モデル、Ｃｕｆｆｌｉｎｋｓモデルに対する推定および真のアルファ_ｉ（存在比、全体のうちの比率）を示す。 正しいおよび誤った転写産物アノテーションのエクソン構造を示す。 ５ｔｉｅ−Ｍｉｘ^２モデルに対するアルファ_ｉの収束を示す。 ５ｔｉｅ−Ｍｉｘ^２モデルに対するベータ_ｊの収束を示す。 ５ｔｉｅ−Ｍｉｘ^２モデルに対するミュー_ｊ（確率サブ関数ｊをシフトさせるためのパラメータ）の収束を示す。 ５ｔｉｅ−Ｍｉｘ^２モデルに対するシグマ_ｊ（確率サブ関数ｊのための幅パラメータ）の収束を示す。 ５ｔｉｅ−Ｍｉｘ^２モデルに対するニュー_ｉ（転写産物ｉに関する確率分布関数のためのシフトまたは並進係数）の収束を示す。 ５ｔｉｅ−Ｍｉｘ^２モデルに対するラムダ_ｉ（転写産物ｉに関する確率分布関数のためのスケーリング係数）の収束を示す。 転写産物１および５ｔｉｅ−Ｍｉｘ^２モデルに対する確率分布関数ｐ（ｒ｜ｔ＝ｉ）の収束を示す。垂直実線は、転写産物の正しい開始および終結点を示す。 転写産物２および５ｔｉｅ−Ｍｉｘ^２モデルに対する確率分布関数ｐ（ｒ｜ｔ＝ｉ）の収束を示す。垂直実線は、転写産物の正しい開始および終結点を示す。 転写産物３および５ｔｉｅ−Ｍｉｘ^２モデルに対する確率分布関数ｐ（ｒ｜ｔ＝ｉ）の収束を示す。垂直実線は、転写産物の正しい開始および終結点を示す。 ５ｔｉｅ−Ｍｉｘ^２モデルに対するサム関数ｐ_{（ｔｏｔａｌ）}（ｒ）の収束を示す。 アルファ_（３）＝０．２を用いた、収束後の５ｔｉｅ−Ｍｉｘ^２モデル、Ｃｕｆｆｌｉｎｋｓモデルに対する推定および真のアルファ_ｉ（存在比、全体のうちの比率）を示す。 アルファ_（３）＝０．４を用いた、収束後の５ｔｉｅ−Ｍｉｘ^２モデル、Ｃｕｆｆｌｉｎｋｓモデルに対する推定および真のアルファ_ｉ（存在比、全体のうちの比率）を示す。 アルファ_（３）＝０．６を用いた、収束後の５ｔｉｅ−Ｍｉｘ^２モデル、Ｃｕｆｆｌｉｎｋｓモデルに対する推定および真のアルファｉ（存在量、全体のうちの比率）を示す。 アルファ_（３）＝０．８を用いた、収束後の５ｔｉｅ−Ｍｉｘ^２モデル、Ｃｕｆｆｌｉｎｋｓモデルに対する推定および真のアルファ_ｉ（存在比、全体のうちの比率）を示す。 アルファ_（３）＝０．２を用いた、収束後の５ｔｉｅ−Ｍｉｘ^２モデル、１ｔｉｅ−Ｍｉｘ^２モデルに対する推定および真のアルファ_ｉを示す。 アルファ_（３）＝０．４を用いた、収束後の５ｔｉｅ−Ｍｉｘ^２モデル、１ｔｉｅ−Ｍｉｘ^２モデルに対する推定および真のアルファ_ｉを示す。 アルファ_（３）＝０．６を用いた、収束後の５ｔｉｅ−Ｍｉｘ^２モデル、１ｔｉｅ−Ｍｉｘ^２モデルに対する推定および真のアルファ_ｉを示す。 アルファ_（３）＝０．８を用いた、収束後の５ｔｉｅ−Ｍｉｘ^２モデル、１ｔｉｅ−Ｍｉｘ^２モデルに対する推定および真のアルファ_ｉを示す。

実施例１：ＮＧＳ方法への導入
元のサンプルにどのｍＲＮＡ分子が存在したかを推測するために、バローズ・ホィーラー変換のような既知の方法を用いてＮＧＳリードをリファレンスゲノム上へマッピングする。これは、リードごとにスプライス部位についての情報を潜在的に含んだ遺伝子座標のセットを与える。マッピング・プロセスを図１に視覚化する。ここではシーケンサーによって作り出されたショートリードの位置がリファレンスゲノム内で同定される。このプロセスをシーケンサーによって生成されたすべてのリードに対して繰り返し、図１において黒く塗りつぶした曲線下の短い直線によって示されるような、多数の短い配列を遺伝子軸上に生じさせる。マッピングしたリードを組み合わせた統計データは、遺伝子軸上に異なるタイプのヒストグラム（すなわち、異なるタイプのカバレッジ包絡曲線）をもたらす。図１における黒く塗りつぶした曲線は、例として、カバレッジ（包絡線）を示す。遺伝子軸上の所定の位置における曲線の値は、その位置をカバーするリードの数である。フラグメント開始点ヒストグラム、別のタイプのカバレッジ包絡曲線のような他のヒストグラムも同様に調べる。各遺伝子位置におけるこのヒストグラムの値は、この位置で開始するリードの数である。本発明の方法は、ヒストグラムの特定のタイプには依存せず、遺伝子軸上および遺伝子座におけるフラグメントのセット上の両方で任意のタイプのヒストグラムに適用可能である。

遺伝子のような遺伝子座は、それぞれがそれら自体のｍＲＮＡを生成する多数の異なる、おそらく重なった領域を含むことができる。かかる領域は、転写産物と称される。遺伝子座におけるヒストグラムｈ_{（ｌｏｃｕｓ）}（ｒ）は、それゆえに、遺伝子座における各個別の転写産物と関連付けられた混合ヒストグラム、すなわち

であり、ここでｈ（ｒ｜ｔ＝ｉ）は、転写産物ｔ＝ｉと関連付けられたヒストグラムであり、ｉは転写産物の名称、例として、Ａｐｏｅである。図２は、図１におけるヒストグラムの２つの転写産物のヒストグラムへの可能な分解を示す。転写産物１は、遺伝子座の全長にわたって広がり、転写産物２の左側で開始する。加えて、転写産物１は、２つのエクソンおよび１つのジャンクションを有する。他方、転写産物２は、転写産物１の第１のエクソンの末端で終結する単一のエクソンを有する。図２における黒く塗りつぶした曲線は、転写産物１のヒストグラムであり、その上の曲線から転写産物１のヒストグラムを差し引いたものが、（生のフラグメントシークエンシングデータからは知られない）転写産物２のヒストグラムである。転写産物ｔ＝ｉの重みα_ｉは、

によって与えられる。これは、全ヒストグラムｈ_{ｌｏｃｕｓ}（ｒ）中で転写産物ｔ＝ｉに属するカウントの割合である。ｈ_{ｌｏｃｕｓ}（ｒ）が、例として、フラグメント開始点のヒストグラムであれば、α_ｉは、転写産物ｔ＝ｉによって生成されたフラグメントのパーセンテージである。この場合、α_ｉは、転写産物ｔ＝ｉの存在比と称され、ｍＲＮＡサンプル内の転写産物ｔ＝ｉの濃度と直接に相関する。図２において、転写産物の重みは、ｈ_{ｌｏｃｕｓ}（ｒ）下の全面積に対するヒストグラムｈ（ｒ｜ｔ＝ｉ）下の面積の割合である。両方の転写産物に対してヒストグラムの面積がおよそ同じなので、転写産物ごとの重みは、約０．５である。従って、各転写産物は、遺伝子座における位置をカバーする約５０％のチャンスを有する。

式（１）におけるヒストグラムｈ_{ｌｏｃｕｓ}（ｒ）の分解は、通常、知られておらず、α_ｉの導出は、精緻な数学的機械を必要とする。分解（１）を見出す問題を数学的に扱いやすくするために（１）を確率論的枠組み、すなわち、

に再定式化し、ここでｐ_{ｌｏｃｕｓ}（ｒ）は、遺伝子座においてフラグメントｒを観測する全確率であり、ｐ（ｒ｜ｔ＝ｉ）は、転写産物ｔ＝ｉを所与として、遺伝子座においてフラグメントｒを観測する確率であり、α_ｉは、遺伝子座において転写産物ｔ＝ｉを観測する確率である。合計して１になるように確率分布ｐ_{ｌｏｃｕｓ}（ｒ）およびｐ（ｒ｜ｔ＝ｉ）を正規化することによってそれらをヒストグラムから導出する場合、（３）が（１）および（２）の直接の結果である。以下に展開する方法が次世代シークエンシングの分野だけでなく、一般的な確率論的設定において適用可能であるという事実を強調するために、下付き文字「ｌｏｃｕｓ」を「ｔｏｔａｌ」によって置き換えることにする。従って
ｐ_{（ｌｏｃｕｓ）}（ｒ）＝ｐ_{（ｔｏｔａｌ）}（ｒ） （４）
である。

転写産物確率α_ｉの推定においては、転写産物分布の形状が基本的に重要である。しかしながら、Ｃｕｆｆｌｉｎｋｓ（非特許文献１）のような現在の方法は、これらの分布を正確にモデリングしない。その代わりに、フラグメント開始点分布に関して、Ｃｕｆｆｌｉｎｋｓは、フラグメント長の分布のみに依存するモデルを用いる。Ｃｕｆｆｌｉｎｋｓにおけるフラグメント長のデフォルト・モデルは、平均値が２００および標準偏差が８０のガウシアンである。Ｃｕｆｆｌｉｎｋｓにおける他の仮定と併せて、これが示唆するのは、２００ｂｐのリード長および３０００ｂｐの転写産物に対してフラグメント開始点分布および転写産物のカバレッジが図３および４に示す関数によってモデリングされることである。

比較として、図５は、図４においてＣｕｆｆｌｉｎｋｓによって仮定された分布から著しくずれたＳｑｕａｒｅライブラリプレップのカバレッジを示す。カバレッジにおける転写産物の５’および／または３’末端へのバイアスは、現在のＮＧＳライブラリの間で共通の特徴であり、Ｃｕｆｆｌｉｎｋｓの仮定はそれゆえに妥当でない。結果として、Ｃｕｆｆｌｉｎｋｓによって推定される転写産物量および濃度は、それらの真の値とは実質的に異なる可能性がある。本明細書に記載する方法は、転写産物分布およびそれらの確率を同時に推定し、それゆえに転写産物量の推定においてかなりより正確である。

遺伝子位置において観測される分布は、どの転写産物が遺伝子位置に含まれるか、およびそれらの確率についてしばしば視覚的な手がかりを呈示する。再び、図２における例を考える。転写産物の平均分布が図６における分布と同様の形状を有し、転写産物確率ｐ（ｒ｜ｔ＝ｉ）がスケーリングおよびシフトによってｐ_ａｖｇ（ｒ）から近似的に導出されることが分かっていれば、図２における重ね合わせの成分を見出すタスクは、全確率分布ｐ_{ｔｏｔａｌ}（ｒ）において平均分布ｐ_ａｖｇ（ｒ）の形状と同様の形状を探すことによって解決できる。図６ではこのプロセスを視覚化し、ｐ_ａｖｇ（ｒ）のスケーリングおよびシフト・バージョンの重ね合わせによるｐ_{ｔｏｔａｌ}（ｒ）へのベスト・フィットが図２におけるのと同様の重ね合わせに繋がり、従って転写産物確率α_ｉの正しい推定をもたらすことを示す。

ｐ（ｒ｜ｔ＝ｉ）は、ｐ_ａｖｇ（ｒ）のスケーリングおよびシフト・バージョンなので、それらは、
ｐ（ｒ｜ｔ＝ｉ）＝ｐ_ａｖｇ（λ_ｉｒ−ν_ｉ） （５）
によって与えられ、ここでλ_ｉおよびν_ｉは、転写産物ｔ＝ｉのスケールおよびシフト・パラメータである。これが意味するのは、転写産物ｔ＝ｉがこの転写産物に固有の２つのパラメータλ_ｉおよびν_ｉのみを有することである。すべての他のパラメータ、すなわちｐ_ａｖｇ（ｒ）のパラメータは、異なるｐ（ｒ｜ｔ＝ｉ）間で連結される。これは、ｐ（ｒ｜ｔ＝ｉ）のパラメータを然るべく連結することによってｐ_ａｖｇ（ｒ）のような異なる転写産物に共通の構造を推定できるという、以下に記載する方法の中心となる考えを際立たせる。以下では、ｐ（ｒ｜ｔ＝ｉ）を混合関数によってモデリングすることになり、それゆえに（３）においてｐ_{ｔｏｔａｌ}（ｒ）から分解されるのは、混合関数のうちの混合関数である。読み易くするために、このモデルを本発明のモデルと呼ぶことにする。以下のセクションは、Ｍｉｘ^２モデルとも称する本発明のモデルを一般的に紹介し、そのいくつかの変形を考察する。加えて、本発明のモデル、Ｍｉｘ^２モデルが確率分布Ｐ（ｒ｜ｔ＝ｉ）を確実に推定できて、Ｃｕｆｆｌｉｎｋｓモデルよりかなり正確なα_ｉに関する推定値を生み出すことを実験が示す。

実施例２：座標変換
２．１．ゲノムにおける位置および転写産物座標における位置
遺伝子軸は、生物に関して配列決定された塩基対の配列であり、通常、０または１で開始し、生物の複雑さによっては数百万塩基対の長さにまで達しうる。加えて、遺伝子軸は、通常、染色体またはコンティグへさらに分割される。遺伝子軸を図５の上部に視覚化し、このグラフィックは、およそ塩基対５３，２４２，５００と５３，２４４，２００との間の染色体１１上のゲノムの選択を表す。転写産物は、通常、遺伝子軸上のエクソンの配列（エクソン_１，．．．，エクソン_Ｎ）として定義し、ｉ番目のエクソンは、ｓ（エクソン_ｉ）で開始し、ｅ（エクソン_ｉ）で終結する遺伝子軸上の区間［ｓ（エクソン_ｉ），．．．，ｅ（エクソン_ｉ）］である。２つの連続するエクソン間のギャップ［ｅ（エクソン_ｉ）＋１，ｓ（エクソン_ｉ＋１）−１］をイントロンと称し、イントロンに先行する最後のヌクレオチドからイントロンの後に続く最初のヌクレオチドへの接続をジャンクションと称する。転写産物の３つの例を図７に示す。この図上のｘ軸は、１０００から５５００に及ぶ遺伝子軸を示し、一方でｙ軸は転写産物ｉｄを示す。このように、転写産物１は、位置１０００で開始し、位置２５００で終結する単一のエクソンからなる。それに対して、転写産物２は、エクソン配列（［１５００，３２００］，［４０００，５０００］）によって定義され、一方で転写産物３は、エクソン配列（［２７００，３２００］，［４０００，５５００］）よって定義される。従って、転写産物２および３は、図７では破線矢印で示す同じジャンクションを有する。以下では、転写産物の長さをｌ（ｔ）によって示すことにする。従って、ｌ（転写産物１）＝１５０１、ｌ（転写産物２）＝２７０２およびｌ（転写産物３）＝２００２である。図７におけるエクソン配列および転写産物の長さをさらに表１にまとめる。

転写産物のエクソン内にある遺伝子軸上の点は、この転写産物の座標系内でも参照できる。転写産物座標は、転写産物の開始からの距離から先行するイントロンの長さを差し引いたものであり、それゆえに１とｌ（ｔ）との間の数である。このように、転写産物のエクソンｅ_ｉにおける位置Ｐに対して、転写産物およびゲノム座標Ｐ_{ｔｒａｎｓ}は、次のように変換できる。

例として、図７の遺伝子座標における位置４５００を考える。この位置は、第２および第３の両方の転写産物のエクソン_２内にある。転写産物２では、この位置の転写産物座標は、４５００−１５００＋１−８０１＝２２００であり、一方で転写産物３では、転写産物座標は、４５００−２７００＋１−８０１＝１０００である。位置４５００は、転写産物１のエクソン内にはないので、この位置を転写産物１のための転写産物座標へ変換することはできない。転写産物のエクソン内にある遺伝子軸上の位置は、転写産物と適合すると称することになる。転写産物と適合する遺伝子（例えば、ゲノム）座標の転写産物の座標への変換は、Ｔによって示すことになる。Ｔが関連付けられる転写産物は、文脈から明らかになるであろう。転写産物座標における確率分布ｐ_{ｔｒａｎｓ}（Ｐ_{ｔｒａｎｓ}｜ｔ＝ｉ）は、遺伝子またはゲノム座標における確率分布へ次のように変換できる。

表記法を簡略化するために、以下では都合がよいときはいつでも、「ｇｅｎｏｍｅ」および「ｔｒａｎｓ」を落とすことにする。図８および９は、転写産物座標およびゲノム座標における確率分布の間の関係を示す。図８は、図７における３つの転写産物のそれぞれに対してスケーリングした確率分布を示す。混合関数における成分としてのそれらの使用を予想して、各確率分布に係数１／３を乗じた。実線は、３つの転写産物のうちで最も短い転写産物１に属し、１５０１塩基対の長さを有する。それに対して、転写産物２および３は、２００２および２７０２塩基対長であり、転写産物座標におけるそれらの確率分布は、それゆえに、１からそれぞれ２００２および２７０２へ及ぶ。

図９は、ゲノム座標における転写産物確率分布を示す。ここでは転写産物の異なる開始点およびジャンクションを考慮に入れた。

ゲノム座標における転写産物確率の和をとると、図１０に示す対象となる遺伝子座全体における位置に関する全確率分布ｐ_{ｔｏｔａｌ}（ｒ）を得る。図１０における曲線は、ＮＧＳデータから導出した遺伝子座における分布の平滑化バージョンと同様である。後のセクションでは、遺伝子座におけるフラグメントの開始点と解釈されるサンプルを生成するために、図１０におけるのと同様の分布を用いることになる。次に、本発明のモデルの確率密度関数のパラメータをトレーニングするために、これらのサンプルを用いることになり、結果として生じたｐ（ｒ｜ｔ＝ｉ）およびｐ_{ｔｏｔａｌ}（ｒ）が図９および１０における曲線の正確な推定値であり、続いて重み係数α_ｉとしての存在比をもたらすことを示すことになる。

２．２ ゲノムおよび転写産物座標におけるフラグメント
フラグメントは、転写産物内の連続的な配列である。転写産物と同様に、フラグメントｒは、それゆえに遺伝子軸上の区間の配列、ｒ＝（ｒｉｎｔ_１，．．．，ｒｉｎｔ_ｋ）からなり、ここでｒｉｎｔ_ｉ＝［ｓ（ｒｉｎｔ_ｉ），ｅ（ｒｉｎｔ_ｉ）］は、開始ｓ（ｒｉｎｔ_ｉ）および終結ｅ（ｒｉｎｔ_ｉ）をもつｉ番目の区間である。フラグメントｒは、その開始および終結が転写産物のエクソン内に位置し、かつ隣接する区間の間のギャップ［ｅ（ｒｉｎｔ_ｉ）＋１，ｓ（ｒｉｎｔ_ｉ＋１）−１］が転写産物のイントロンであれば、すなわち、いくつかのｉ≦ｋに対して、
ｓ（ｉｎｔ_１）∈ｅｘｏｎｉ （９）
ｅ（ｉｎｔ_Ｋ）∈ｅｘｏｎｋ （１０）
であり、
［ｅ（ｉｎｔ_ｋ）＋１，ｓ（ｉｎｔ_ｋ＋１）−１］∈｛［ｅ（ｅｘｏｎ_ｉ）＋１，ｓ（ｅｘｏｎ_ｉ＋１）−１］：ｉ＝１，．．．，Ｎ−１｝∀ｋ＝１，．．．，Ｋ−１ （１１）
であれば、その転写産物と適合する。

フラグメントｒが転写産物と適合すれば、それを転写産物座標における区間［Ｔ（ｓ（ｒｉｎｔ_１）），Ｔ（ｅ（ｒｉｎｔ_Ｋ））］へ変換できる。フラグメントの区間の間のギャップは、それゆえに転写産物座標への変換において除去される。前セクションにおけるように、フラグメントｒの転写産物座標への変換は、Ｔによって示すことになる。例として、フラグメント（［２０００，３０００］，［４０００，４５００］）および図７における３つの転写産物を考える。このフラグメントの開始および終結は、それぞれ転写産物２の第１および第２のエクソン中にある。同時に、フラグメントのギャップ［３００１，３９９９］は、転写産物２のイントロンと合致する。このように、フラグメントは、転写産物２と適合し、転写産物座標における区間［５０１，２２００］へ変換される。それに対して、転写産物１および転写産物３に対しては、フラグメントの開始または終結のいずれかがそのエクソンの１つの内にはなく、フラグメントは、それゆえにこれらの転写産物のいずれとも適合しない。

転写産物座標におけるフラグメント上の確率分布ｐ_{ｔｒａｎｓ}（ｒ_ｔａｎｓ｜ｔ＝ｉ）は、それゆえにゲノム座標におけるフラグメント上の確率分布へ次のように変換できる。

前の通り、「ｇｅｎｏｍｅ」および「ｔｒａｎｓ」の下付き文字は、都合が良いときはいつでも落とすことにする。その開始および終結を用いる代わりに、区間をその開始ｓ（ｒ）および長さｌ（ｒ）によって表すこともできる。これは、転写産物座標における確率分布の次の因数分解を可能にする。

フラグメント長の包括的な分布は、通常、ライブラリのバイオアナライザ・トレースから推測できるので、これらの因数分解は、ＮＧＳデータにとって便利である。（１３）における因数分解は、Ｃｕｆｆｌｉｎｋｓによって用いられ、Ｃｕｆｆｌｉｎｋｓは、所定の長さのフラグメントを転写産物上のどこにでも等しい確率で配置できることを追加的に仮定する。対照的に、セクション４の実験において、Ｍｉｘ^２モデルは、（１４）における因数分解を用い、それによってフラグメント開始点の転写産物特有のばらつきのより効率的な推定が可能になる。

実施例３：転写産物確率の推定および転写産物特有の確率分布
以下に記載するモデルは、混合関数を混合した関数を用い、それゆえにＭｉｘ^２モデルと称することにする。

３．１ Ｍｉｘ^２モデルの数学的基礎
以下では、ｒは、フラグメントおよび位置の両方を表すことができる。しかしながら、便宜上、ｒを常にフラグメントと呼ぶことにする。遺伝子座において特定のフラグメントｒを観測する確率ｐ_{ｔｏｔａｌ}（ｒ）は、転写産物に対してフラグメントを観測する確率に、その転写産物がフラグメントを生成する確率を重み付けしたものの和である。従って、ｐ_{ｔｏｔａｌ}（ｒ）は、次の混合確率分布によって与えられる。

セクション２において記載したように、ｒがｔ＝ｉと適合すればｐ（ｒ｜ｔ＝ｉ）＝ｐ_{ｔｒａｎｓ}（Ｔ（ｒ）｜ｔ＝ｉ）であり、そうでない場合にはｐ（ｒ｜ｔ＝ｉ）＝０である。本方法は、確率分布ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）が混合確率分布、すなわち、

であることを仮定する。

ここでＭ_ｉは混合成分の数であり、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）は確率分布である。β_ｉ，ｊ≧０は、和が１になる重みであり、すなわち、

である。

前の通り、変数ｔが転写産物を表すものとし、一方で新たに導入した変数ｂは、一般に確率サブ関数とも呼ばれる「ビルディングブロック」を表すものとする。このモデルでは、フラグメントｒを所与として転写産物ｔ＝ｉおよびブロックｂ＝ｊが観測された事後確率を算出することが可能である。この事後確率は、

によって与えられる。

フラグメントｒが転写産物ｔ＝ｉと適合しなければ、ｐ（ｔ＝ｉ，ｂ＝ｊ）＝０、従ってｐ（ｒ｜ｔ＝ｉ，ｂ＝ｊ｜ｒ）＝０であることに留意すべきである。ｔおよびｂに関する事後確率は、

によって与えられる。

そのうえ、モデルは、ｐ_{ｔｏｔａｌ}（ｒ）のいくつかのパラメータが異なる転写産物（確率分布関数）およびブロック（確率サブ関数）の間で連結されることを仮定する。これは、混合重みβ_ｉ，ｊと、確率分布ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ），ｉ＝１，．．．，Ｎ，ｊ＝１，．．．，Ｍの確率サブ関数のセットを確定する任意のパラメータΘとの両方を含むことができる。連結は、組み合わせたすべての転写産物およびブロックに対して適用されてもよく、または転写産物およびブロックの群内で適用されてもよい。確率分布ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）のパラメータを連結することは、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）間のある類似性を示唆するため、パラメータ連結は、パラメータ連結のタイプによって示唆されるような類似性を呈示する転写産物に対してのみ適用されるべきである。例として、導入において提案したように、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）が単一のｐ_ａｖｇ（ｒ）からスケーリングおよびシフトによって導出されたことをパラメータ連結が示唆するならば、パラメータ連結は、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）が相互のスケーリングおよびシフト・バージョンである転写産物に対してのみ適用されるべきである。いくつかのライブラリプレップでは、後者は、ある範囲の長さをもつ転写産物に対してのみ有効である。それゆえに、その長さが異なる範囲内にある転写産物には、ｐ（ｒ｜ｔ＝ｉ）が相互のスケーリングおよびシフト・バージョンであることを必要とするパラメータ連結を適用すべきではない。

遺伝子領域におけるフラグメントｒのデータセットＲに関して、データセットＲの尤度をモデルｐ_{ｔｏｔａｌ}（ｒ）の下で最大化することによって、Ｍｉｘ^２モデルの転写産物特有の分布ｐ（ｒ｜ｔ＝ｉ）および転写産物確率α_ｉを推定する。パラメータα_ｉ、β_ｉ，ｊ、Θの尤度関数は、それゆえに

によって与えられる。

（２１）の最大化は、多数の異なる最適化方法を用いて近似解を見出すことができる制約付き非線形最適化問題である。以下で用いる最適化方法は、期待値最大化（ＥＭ：ｅｘｐｅｃｔａｔｉｏｎ ｍａｘｉｍｉｚａｔｉｏｎ）アルゴリズムである。ＥＭアルゴリズムは、尤度関数の極大を見出す反復法である。大域的最適に近い極大を得るために、多分必要なのは、モデルパラメータに対して異なる初期化を用いて、ＥＭアルゴリズムを数回試みることである。次に、最大尤度をもたらす結果が最適化問題の解として選ばれる。ＥＭアルゴリズムは、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）の形状またはそれらのパラメータの連結についての仮定なしに、α_ｉを推定するために用いることができる。Ｃｕｆｆｌｉｎｋｓでも用いられる、α_ｉに関するＥＭ更新式は、

によって与えられる。

ここで｜Ｒ｜は、データセットＲにおけるフラグメントの数であり、α_ｉ ^{（ｎ＋１）}は、ｎ＋１回目の反復後のα_ｉに関する推定値であり、ｐ^（ｎ）（ｔ＝ｉ｜ｒ）は、フラグメントｒを所与として、ｎ回目の反復からのパラメータを用いて推定した転写産物ｔ＝ｉを観測する事後確率である。この場合、最適化問題は凹なので、異なる初期化を用いてＥＭアルゴリズムを繰り返す必要はない。凹問題では、ＥＭアルゴリズムが常に収束する単一の大域的最適が存在する。モデルパラメータは、モデルｐ_{ｔｏｔａｌ}（ｒ）の観測データＲへのフィットを最適化することによって推定する。しかしながら、最終目標は、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）への良好なフィット、延いては転写産物確率α_ｉの良好な推定値を見出すことである。ｐ_{ｔｏｔａｌ}（ｒ）の最適化それ自体は、α_ｉの良好な推定値を示唆しない。ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）のパラメータが適切に連結されている場合に限り、ｐ_{ｔｏｔａｌ}（ｒ）の最適化は、転写産物確率α_ｉの良好な推定値をもたらすであろう。次のセクションは、異なる連結方策を用いたＭｉｘ^２モデルのいくつかの変形を紹介する。セクション４ではそれらの変形のうちの２つを調べて、スケーリングおよびシフトした転写産物確率分布ｐ（ｒ｜ｔ＝ｉ）の場合、２つの変形がＣｕｆｆｌｉｎｋｓモデルよりかなり正確なα_ｉに関する推定値をもたらすことを示すことになろう。

３．２ Ｍｉｘ^２モデルの変形
先のセクションにおけるＭｉｘ^２モデルの数学的基礎は、かなり一般的である。このセクションは、多くの具体的に実現したＭｉｘ^２モデルを考察して、その様々な利点および限界を際立たせる。

３．２．１ 単一のパラメータ群内の連結
このセクションにおいて考察する最も簡単なＭｉｘ^２モデルは、異なる転写産物間のｉ＝１，．．．，Ｎに対して重みβ_ｉ，ｊのみを連結する。従って、このモデルのパラメータのセットは、｛α_ｉ，β_ｊ：ｉ＝１，．．．，Ｎ，ｊ＝１，．．．，Ｍ｝であり、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）は、

によって与えられる。

このＭｉｘ^２モデルは、単一のパラメータ群内のパラメータ、すなわちβ_ｊのみを連結するので、以下ではこのモデルを１ｔｉｅ−Ｍｉｘ^２と呼ぶことにする。１ｔｉｅ−Ｍｉｘ^２モデルの下でフラグメントｒを観測する確率は、

によって与えられる。

式（２４）では和の順序を交換することができて、

である。前述のように、（２６）における被加数は、フラグメントｒが転写産物ｔ＝ｉと適合しなければゼロである。ｐ（ｒ｜ｂ＝ｊ）が確率分布なので、式（２５）は、α_ｉを所与として、β_ｊを推定する問題がα_ｉの推定と概念的に同じであることを示す。このように、ＥＭアルゴリズムを適用することができて、β_ｊに関する次の更新式を与える。

これは、式（２２）および（２７）の反復的な適用によって、１ｔｉｅ−Ｍｉｘ^２モデルのα_ｉおよびβ_ｊの両方を同時にトレーニングできることを示唆する。Ｃｕｆｆｌｉｎｋｓモデルと比較して、１ｔｉｅ−Ｍｉｘ^２モデルの尤度関数が凹ではないため、原理的に、１ｔｉｅ−Ｍｉｘ^２モデルに対するＥＭアルゴリズムは、α_ｉおよびβ_ｊの初期化に依存して異なる解へ収束する可能性がある。しかしながら、セクション４における実験では満足すべき結果へ収束するのにα_ｉおよびβ_ｊの単一の初期化で十分であった。データをフィッティングするためには、１ｔｉｅ−Ｍｉｘ^２モデルに係る確率分布ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）を選ばなければならない。従って、データの構造について予備知識が必要である。観測したｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）が異なるｔ＝ｉに対する相互のスケーリング・バージョンであれば、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）も同様に異なるｔ＝ｉに対するスケーリング・バージョンでなければならない。加えて、同じｂ＝ｊ、しかし異なるｔ＝ｉに対して、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）は、同じスケーリング係数β_ｊによって調節される転写産物座標系の領域内になければならない。後者は、転写産物の正しい開始および長さが分かっていれば達成できる。この場合、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）をｔ＝ｉの転写産物座標系に沿って等距離の位置に配置できる。正しい開始および長さが分からなければ、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）をどのように置くかが直ちに明確ではない。次のセクションにおけるモデルは、転写産物の正しい開始および長さをデータから知ることによってこの制限を回避する。

３．２．２ ５つのパラメータ群内の連結
先のセクションにおけるモデルは、転写産物の正しい開始点および長さの知識に依存する。このセクションは、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）を自動的に配置およびスケーリングするモデルを考察する。このモデルは、５つのパラメータ群内のパラメータを連結するので、以下では５ｔｉｅ−Ｍｉｘ^２モデルと称することにする。５ｔｉｅ−Ｍｉｘ^２モデルは、分布ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）に対してガウシアンを用い、その内部パラメータ、すなわち、それらの平均値μ_ｉ，ｊおよび標準偏差σ_ｉ，ｊは、パラメータのセット

から導出する。特に、μ_ｉ，ｊおよびσ_ｉ，ｊは、

によって与えられる。

転写産物ｔ＝ｉおよびブロックｂ＝ｊに関するｐｄｆは、それゆえに

によって与えられる。前の通り、ブロックｊに対してβ_ｉ，ｊを転写産物ｉ＝１，．．．，Ｎ間で連結し、従ってβ_ｉ，ｊ＝β_ｊである。確率分布

は、転写産物ｔ＝ｉに関するｐｄｆ ｐ（ｒ｜ｔ＝ｉ）を得るためにν_ｉによってシフトさせてλ_ｉによってスケーリングした、遺伝子座における転写産物の平均ｐｄｆ ｐ_ａｖｇ（ｒ）として解釈でき、ここでｇ_{μｊ，σｊ}（ｒ）は、平均値μｊおよび標準偏差σｊをもつガウシアンである。モデルの能力を制限することなく、選択した転写産物に関してλ_ｉ＝１およびν_ｉ＝０とセットすることが可能である。式（３２）におけるｐｄｆは、転写産物ｔ＝ｉの座標系ではなく実際の連続空間における確率分布である。従って、結果として生じるｐｄｆ ｐ（ｒ｜ｔ＝ｉ）およびｐ_{ｔｏｔａｌ}（ｒ）も連続的なｐｄｆである。この点を強調するために、以下ではこれらの確率分布を

によって示し、

の重みを

によって示すことになる。これらの連続的な確率分布から、転写産物ｔ＝ｉに関して

を、可能なフラグメントＦ（ｔ）のセットに制限して、それらの値を和が１になるように、すなわち、

に正規化することによって、ゲノムおよび転写産物座標におけるｐ_{ｔｏｔａｌ}（ｒ）およびｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）を算出することができ、ここで（３３）の分母の和における和集合は、遺伝子座におけるすべての転写産物に及ぶ。次に、α_ｉを

から次のように導出できる。

式（３５）は、連続的なモデル重みアルファ_{（ｉ，Ｒ）}のリスケーリングになる。リスケーリングの効果は、ゆっくりと変化する確率分布ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）、およびすべての転写産物がかなりの長さを有する遺伝子座に対してはあまり顕著ではない。従って、後で記載する実験ではリスケーリング式の効果を無視できる。遺伝子座およびｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）がこれらの必要条件に適合する状況では、それゆえに転写産物確率α_ｉとしてアルファ_{（ｉ，Ｒ）}を直接に用いることができ、セクション４においてこれを行うことになる。λ_ｉは、転写産物の異なる長さを考慮するスケーリング・パラメータであり、一方でν_ｉは、転写産物開始点の転写産物アノテーションとは異なる相対的なオフセットを考慮する。これら２つのパラメータは、転写産物アノテーションにおける誤った開始点および長さを補正できるため、転写産物アノテーションから開始および終結の定義を一斉に除去することが可能である。例として、すべての転写産物の開始点を遺伝子座における転写産物のすべての開始点より小さい同じ値にセットできる。同様に、すべての転写産物の終結点を遺伝子座における転写産物のすべての終結点より大きい同じ値にセットできる。転写産物アノテーションのこの拡張を図１１に示す。転写産物アノテーションを拡張することは、転写産物開始または終結と適合しなかったためにこれまで転写産物に対して無効であったフラグメントがこの転写産物に対して今や有効になるという効果を有する。図１１では、例として、転写産物２に対して有効なすべてのフラグメントが転写産物３に対しても有効である。加えて、イントロンに架からないか、または入らないすべてのフラグメントは、すべての３つの転写産物に対して有効である。

このセクションにおけるモデルのパラメータは、ＥＭアルゴリズムを用いて効率的にトレーニングできる。平均値、分散およびオフセット・パラメータに関するＥＭ更新式を以下に示す。

スケール・パラメータλ_ｉに対してＥＭアルゴリズムの補助関数の微分を算出すると、以下の２次式が得られる。

ほとんどの実際的な環境では、式（３９）は、１つの正および１つの負の解を有する。λ_ｉは、スケール・パラメータなので正でなければならず、それゆえに（３９）の正の解のみが対象となる。セクション４における実験は、５ｔｉｅ−Ｍｉｘ^２モデルが、実際、不正確な転写産物アノテーションを補正でき、この場合、正しい転写産物アノテーションに依存するモデルよりはるかに正確なα_ｉに関する推定値をもたらすことを示す。

３．２．３ ２つのパラメータ群内の連結
前モデルのように、このセクションにおけるモデルは、ブロックｊに対してβ_ｉ，ｊを転写産物にわたって連結し、それゆえにβ_ｉ，ｊ＝β_ｊである。加えて、このセクションにおけるモデルは、５ｔｉｅ−Ｍｉｘ^２モデルからの

を連結するが、所定のブロックに対して転写産物にわたって連結するのではなく、より一般的に（ｉ，ｊ），ｉ＝１，．．．，Ｎ，ｊ＝１，．．．，Ｍの対にわたって連結する。例として、Ｎ≧２，Ｍ≧３であれば、可能なセットは、Ｌ＝｛（１，１），（１，２），（２，３）｝であり、従って、

である。このセクションにおけるモデルは、２つの群内のパラメータのみを連結するので、２ｔｉｅ−Ｍｉｘ^２モデルと称することにする。このように、２ｔｉｅ−Ｍｉｘ^２モデルに関して、
σ_ｉ，ｊ＝λ_ｉσ_Ｌ（ｉ，ｊ）∈Ｌ （４０）
である。

２ｔｉｅ−Ｍｉｘ^２モデルに対する動機付けは、以下の通りである。図１２に示す状況を考える。この図では、１０００ｂｐで開始して４０００ｂｐで終結する領域内にジャンクションが何もない２つの転写産物がある。第２の転写産物は、全領域に及び、一方で第１の転写産物は、位置１５００で開始して位置３５００で終結する。各転写産物のｐｄｆ（ｐ（ｒ｜ｔ＝ｉ））が図１２におけるガウシアンより上の２つの曲線の形状を有すると想定する。しかしながら、図１２の下部に示す全ｐｄｆ ｐ_{ｔｏｔａｌ}（ｒ）は、厳密には２つの転写産物に関するｐ（ｒ｜ｔ＝ｉ）の和ではない。その代わりに、ｐ_{ｔｏｔａｌ}（ｒ）は、２２５０ｂｐと２７５０ｂｐの間の領域にくぼみを有する。これは、ＮＧＳシークエンシングの共通の特徴であり、ある領域内のヌクレオチド配列がこの領域と交差しているすべての転写産物に係るフラグメントの生成に影響を及ぼしうる。シークエンシングを開始するシークエンシング・プライマーにある配列に結合する能力が欠如していると、このような結果になりうる。図１２における例では、例として、くぼみの領域内の配列がシークエンシング・プライマーにはより不利であり、それゆえに、この領域における全ｐｄｆ ｐ_{ｔｏｔａｌ}（ｒ）が著しく減少する。この問題に対する解は、くぼみの領域と交差していない他の転写産物によって共有されるかもしれないβ_ｊの連結ではなく、その領域との交差を有するブロックｐ（ｒ｜ｔ＝ｉ，ｂ＝ｊ）に属する

の連結である。図１２では、これらのブロックを破線で示す。すべてのこれらのブロックの

を非常に大きい値にセットすれば、連続的なｘ軸全体にわたってガウシアンｐ（ｒ｜ｔ＝ｉ，ｂ＝ｊ）が広がり、ゲノムの遺伝子座における全ｐｄｆ ｐ_{ｔｏｔａｌ}（ｒ）にはほとんど何も寄付しない。ｐ_{ｔｏｔａｌ}（ｒ）においてくぼみが生じうる領域が予め知られていなければ、互いに近接したブロックに属するすべての

を連結できる。本例では、転写産物１のための２つの隣接ブロックの

を、それらの上にある転写産物２のブロックの

と連結してもよく、例えば、Ｌ＝｛（１，１），（１，２），（２，３）｝であり、従ってこのセクションの初めにおけるように、

である。σ_Ｌに関するＥＭ更新式は、式（３７）から導出でき、

によって与えられる。

モデルの他のパラメータに関するＥＭ更新式は、前セクションにおけるのと同じである。図１２の状況では、β_ｊを更新することが賢明である。他のパラメータの更新は、（ｉ，ｊ）∈Ｌに対するｐ（ｒ｜ｔ＝ｉ，ｂ＝ｊ）の空間的な関連性を破壊しかねず、それゆえに潜在的に有害である。

３．２．４ ６つのパラメータ群内の連結
５ｔｉｅ−Ｍｉｘ^２の直接的な拡張は、

をスケーリングする別のパラメータ

の導入によって与えられる。特に、

である。

このように、５ｔｉｅ−Ｍｉｘ^２モデルと比較して、このセクションにおけるモデルは、

を独立にスケーリングし、パラメータを６つのパラメータ群内で連結するので、６ｔｉｅ−Ｍｉｘ^２モデルと称する。６ｔｉｅ−Ｍｉｘ^２モデルでは、

を、同じ転写産物に対して異なるブロック間で連結し、従って

である。６ｔｉｅ−Ｍｉｘ^２モデルにおいて、λ_ｉおよびκ_ｉは、平均値μ_ｊおよび標準偏差σ_ｊと同様の役割を果たし、式（３６）、（３７）と同様の構造を有する以下の式を用いてトレーニングできる。

５ｔｉｅ−Ｍｉｘ^２モデルを凌ぐ６ｔｉｅ−Ｍｉｘ^２モデルの利点は、異なる転写産物間の平滑度の違いをモデリングできることである。しかしながら、パラメータ数の増加が不十分なモデル収束に繋がる可能性もあり、それゆえに６ｔｉｅ−Ｍｉｘ^２モデルによって示唆される平滑度における分散をＮＳＧデータが呈示する場合にのみ用いるべきである。

３．２．５ 高次Ｍｉｘ^２モデル
前セクションでは、式（３０）および（３１）によって与えられるλの転写産物特有のアフィン線形変換を通じてμ_ｊおよびσ_ｊから、ガウシアンｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）の平均値μ_ｉ，ｊおよび標準偏差σ_ｉ，ｊを導出した。このセクションは、この概念を一般化して、μ_ｉ，ｊおよびσ_ｉ，ｊをλにおける多項式、すなわち、

によって生成することを可能にする。

μ_ｊ（ｋ_０）に関するＥＭ更新式は、以下の表現によって与えられる。

残りのモデルパラメータに関するＥＭ更新式は、非線形方程式をもたらし、これらの式は、反復法によってのみそれら自体を解決できる。

実施例４：スケーリングおよびシフトした転写産物特有の確率分布に関する実験
このセクションにおける実験では、確率ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）を次のように因数分解する。
ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）＝ｐ_{ｔｒａｎｓ}（ｓ（ｒ）｜ｔ＝ｉ，ｂ＝ｊ）ｐ_{ｔｒａｎｓ}（ｌ（ｒ）｜ｔ＝ｉ，ｂ＝ｊ，ｓ（ｒ）） （４８）
ここでｓ（ｒ）およびｌ（ｒ）は、フラグメントｒの開始および長さである。加えて、確率ｐ_{ｔｒａｎｓ}（ｌ（ｒ）｜ｔ＝ｉ，ｂ＝ｊ，ｓ（ｒ））がｓ（ｒ）および転写産物の長さｌ（ｔ）にのみ依存すると仮定する。それゆえに（４８）は、
ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）＝ｐ_{ｔｒａｎｓ}（ｓ（ｒ）｜ｔ＝ｉ，ｂ＝ｊ）ｐ_{ｔｒａｎｓ}（ｌ（ｒ）｜ｌ（ｔ），ｓ（ｒ）） （４９）に帰着する。

結果として、混合モデル（１６）は、開始点分布ｐ_{ｔｒａｎｓ}（ｓ（ｒ）｜ｔ＝ｉ）にフラグメント長に関する確率分布を乗じた混合モデル、すなわち、

である。

このように、モデルの正しいｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）への収束は、分布

の正しい開始点分布への収束およびα_ｉの正しい存在比への収束をチェックすることによって評価できる。

このセクションでは２つのタイプの実験を考察する。第１のタイプは、１００００個のフラグメントのセットに対するＭｉｘ^２モデルの収束の詳細な解析であり、重みα_１＝０．２８、α_２＝０．３２およびα_３＝０．４を用いた図９におけるｐ（ｒ｜ｔ＝ｉ）の重ね合わせからフラグメント開始点をランダムに取り出した。加えて、フラグメント開始点ｓ（ｒ）ごとに、ｓ（ｒ）、すなわち、１，．．．，ｌ（ｔ）−ｓ（ｒ）＋１で開始するフラグメントの可能な長さへ正規化した平均値２００および標準偏差８０のガウシアンから、フラグメント長をランダムに取り出した。サンプリングした開始点およびフラグメント長のヒストグラムを図１３および１４に示す。図からわかるように、開始点ヒストグラムは、図１０における分布の形状におおよそ従い、一方でフラグメント長ヒストグラムは、ガウシアンの形状に従う。

このセクションにおける第２のタイプの実験は、重みα_１，α_２，α_３の６０個の異なるセットに対して第１のタイプを繰り返し、Ｍｉｘ^２モデルによって推定したα_１，α_２，α_３を真の重みと比較する。前の実験は、Ｍｉｘ^２モデルの収束およびそれらのパラメータについて詳細な知見を提供し、一方で後の実験は、Ｍｉｘ^２モデルによって推定する重みα_１，α_２，α_３の総合的な精度に注目する。α_３をそれぞれ０．２、０．４、０．６、０．８の値にセットすることによって６０個の重みのセットを選んだ（ｗｈｅｒｅ ｃｈｏｓｅｎ）。α_３の値ごとに、α_１の１５個の異なる値を０と１−α_３との間に等距離間隔で選び、α_２は、α_２＝１−α_１−α_３にセットした。重みのセットごとに、真の重みと推定した重みとの間のカルバック・ライブラー（ＫＬ：Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ）ダイバージェンスを算出し、重みのすべてのセットにわたるＫＬダイバージェンスの平均値を尺度として用いて、ＣｕｆｆｌｉｎｋｓおよびＭｉｘ^２モデルの精度を定量化した。

４．１ １ｔｉｅ−Ｍｉｘ^２モデル
４．１．１ α_１＝０．２８、α_２＝０．３２、α_３＝０．４に対する１ｔｉｅ−Ｍｉｘ^２モデルの収束
前述のように、このセクションにおいてモデルを推定するために用いたデータは、重みα_１＝０．２８、α_２＝０．３２、α_３＝０．４を用いた図９におけるｐ（ｒ｜ｔ＝ｉ）の重ね合わせからフラグメント開始点をサンプリングすることによって生成し、一方でフラグメント長は、再正規化した平均値２００および標準偏差８０のガウシアンからサンプリングした。これらのサンプリング・データセットのヒストグラムを図１３および図１４に示す。このセクションにおける１ｔｉｅ−Ｍｉｘ^２モデルは、ガウシアンであるように選んだ８つのビルディングブロックｂ_ｉ，ｊ（ｓ（ｒ））を用いる。ガウシアンの平均値は、転写産物の長さにわたって均等に分布し、

の位置で開始して、

のステップで進み、ここで両方の数を最も近い整数へ四捨五入した。このように、例として、長さが２７０２ｂｐの転写産物２に関して、このモデルは、１６９、５０７、８４５、１１８３、１５２１、１８５９、２１９７、２５３５における平均値を生じた。各ガウシアンの標準偏差は、第１の平均値に等しいようにセットした。従って、転写産物２に関して、各ガウシアンの標準偏差は、１６９であった。モデルにおける初期のβ_ｊは、１／８にセットした。転写産物２に関して結果として生じたｐ（ｒ｜ｔ＝ｉ）の初期化を、重み付けしたブロックβ_ｊｐ（ｓ（ｒ）｜ｔ＝ｉ，ｂ＝ｊ）とともに図１５に示す。最初および最後のブロックは、転写産物内のブロックより若干高い。これは、これらのガウシアンがそのすその部分を失っており、それらの正規化定数がそれゆえに他のガウシアンより高いという事実に起因する。

モデルにおける初期の存在比は、α_１＝α_２＝α_３＝１／３にセットした。後続の反復間のα_ｉおよびβ_ｊの間の差が０．００１未満か、または後続の反復間の対数尤度における増加が０．５未満となるまでＥＭアルゴルズムを適用した。これらの終了条件を用いると、ＥＭアルゴリズムが２０回の反復後に収束した。

図１６は、α_ｉの収束を示す。この図のｘ軸はＥＭアルゴリズムの反復を表し、一方でｙ軸はα_ｉの値を表す。図１６における３つの曲線は、すべて初期値１／３から始まり、図１６に３つの水平線で示す正解に近い値へ収束する。最終値は、α_１＝０．２９、α_２＝０．３、α_３＝０．４１であった。比較として、Ｃｕｆｆｌｉｎｋｓモデルを用いて導出した推定値は、α_１＝０．３７、α_２＝０．１２、α_３＝０．５１である。これらの数をα_ｉの真の値と推定値との間のカルバック・ライブラー（ＫＬ）ダイバージェンスとともに表２にまとめる。ＫＬダイバージェンスは、同一の分布に対して０であり、分布間の差の増加とともに増加する。表２は、Ｃｕｆｆｌｉｎｋｓモデルに対する真のα_ｉと推定されたα_ｉとの間のＫＬダイバージェンスが１ｔｉｅ−Ｍｉｘ^２モデルのＫＬダイバージェンスより２桁優ることを示す。

図１７は、８個のブロックに関する８つのベータの収束を示す。前の通り、ｘ軸はＥＭアルゴリズムの反復を表し、ｙ軸はβ_ｊの値を表す。β_ｊはすべて１／８に初期化したので、曲線は、同じ位置から始まる。α_ｉの場合のように、β_ｊに関するグランドトルースは何もないので、それらの品質は、結果として生じた分布ｐ（ｒ｜ｔ＝ｉ）をチェックすることによってのみ評価できる。最終的なβ_ｊをそれらのブロックｉｄに対してプロットした図１８に示す。これは、図１９から分かるように、β_ｊがｐ（ｒ｜ｔ＝ｉ）の形状におおよそ従い、次には正しいｐ（ｒ｜ｔ＝ｉ）の良好な近似に繋がることを示す。

最後に、図２０は、全確率分布ｐ_{ｔｏｔａｌ}（ｒ）の収束を示す。これは、総合的な尤度の観点から、ＥＭアルゴリズムが最初の反復において正解に向けた最も大きなステップを進むことを示す。その後の反復は、総合的な尤度の増加に対して比較的わずかな効果を有する。しかしながら、初期の反復では、α_ｉはそれらの正解からまだはるかに遠い。これは、総合的な尤度におけるわずかな変化もパラメータの有意な変化に繋がりうることを示す。結果として、ＥＭアルゴリズムの終了条件は、全尤度の小さい増加に対して寛容であるべきである。

４．１．２ 重みの６０個のセットに関する１ｔｉｅ−Ｍｉｘ^２およびＣｕｆｆｌｉｎｋｓモデルの間の比較
このセクションは、セクション４の初めに記載した手順に従って選んだ、重みα_１、α_２、α_３の６０個の異なるセットを用いた実験を考察する。従って、α_３は、０．２、０．４、０．６、０．８の値をとり、α_１およびα_２は、０と１−α_３との間に等距離間隔で分布する。重みのセットごとに、図１５におけるように転写産物２に対してｐ（ｒ｜ｔ＝ｉ）を初期化した。他の２つの転写産物に関する分布ｐ（ｒ｜ｔ＝ｉ）をそれに応じて初期化した。収束判定基準を満たすまでＥＭアルゴリズムを実行した。最終的な反復で得られたα_１を１ｔｉｅ−Ｍｉｘ^２モデルの結果として選んだ。同様に、Ｃｕｆｆｌｉｎｋｓモデルに対して収束判定基準のうちの１つを満たすまでＥＭアルゴリズムを行い、最終的な反復からのα_ｉをＣｕｆｆｌｉｎｋｓモデルの結果として選んだ。図２１は、α_３＝０．２に対するこれらの実験の結果を示す。このグラフのｘ軸はα_１の真の値を示し、ｙ軸は真のα_ｉならびにＣｕｆｆｌｉｎｋｓおよび１ｔｉｅ−Ｍｉｘ^２モデルによって推定したα_ｉを示す。図２１における一点鎖線は、真のα_ｉを指し、一方で点線および破線は、それぞれＣｕｆｆｌｉｎｋｓおよび１ｔｉｅ−Ｍｉｘ^２の推定値を指す。図２１は、α_３＝０．２に対して、１ｔｉｅ−Ｍｉｘ^２の推定値がα_ｉの真の値と非常によく一致し、一方でＣｕｆｆｌｉｎｋｓの推定値がむしろ不十分であることを示す。α_２＝０に対してのみ、Ｃｕｆｆｌｉｎｋｓおよび１ｔｉｅ−Ｍｉｘ^２モデルの推定値が符合する。この状況では、図７からわかるように、完全に分離した転写産物１および３のみが存在する。この場合には、それゆえにＥＭアルゴリズムは、転写産物１および３に割り当てられたフラグメントの数を単にカウントし、ｐ（ｒ｜ｔ＝ｉ）の形状には依存しない。図２２は、α_３＝０．４に対して、Ｃｕｆｆｌｉｎｋｓおよび１ｔｉｅ−Ｍｉｘ２モデルの推定値を示す。α_３＝０．２については、１ｔｉｅ−Ｍｉｘ^２モデルの推定値は、Ｃｕｆｆｌｉｎｋｓモデルと比較して非常に正確である。α_２＝０に対してのみ、両方の推定値が符合する。同様の描像がα_３＝０．６およびα_３＝０．８に関して浮かび上がる。

表３から分かるように、これは、総合的にＣｕｆｆｌｉｎｋｓモデルでは平均ＫＬダイバージェンス０．１２３６８および１ｔｉｅ−Ｍｉｘ^２モデルでは平均ＫＬダイバージェンス３．６３６９ｅ−０４に繋がる。このように、平均ＫＬダイバージェンスの観点から、１ｔｉｅ−Ｍｉｘ^２モデルの精度は、Ｃｕｆｆｌｉｎｋｓモデルより３桁上回る。

４．２ ５ｔｉｅ−Ｍｉｘ^２モデル
このセクションにおけるモデルは、転写産物特有のオフセットおよびスケーリング・パラメータν_ｉおよびλ_ｉを推定し、それゆえに、正しい転写産物アノテーションには依存しない。このように、このモデルの可能性を実証するために、３つの誤った転写産物アノテーションを用いてこのモデルをトレーニングする。これらのアノテーションを図２５に見ることができる。この図における実線は、図７および表１と同じ正しい転写産物アノテーションを示す。図２５における点線は、このセクションにおけるモデルをトレーニングするために用いた誤った転写産物アノテーションを示す。誤った転写産物アノテーションにおけるエクソンの厳密な開始および終了位置ならびに正しいアノテーションとのそれらの差を表４に示す。このセクションにおける５ｔｉｅ−Ｍｉｘ^２モデルは、位置１で開始し、位置１００００で終結する拡張した転写産物アノテーションを用いる。従って、転写産物２および３の拡張したアノテーションは、同一である。５ｔｉｅ−Ｍｉｘ^２モデルのν_ｉおよびλ_ｉを誤ったアノテーションと適合するように初期化して表５に示す。転写産物２のλ_ｉを１に選び、それゆえに転写産物１および３のλ_ｉは、１２７５／２９００＝０．４３９７および２３００／２９００＝０．７９３１によって与えられる。５ｔｉｅ−Ｍｉｘ^２モデルのパラメータの複雑な相互作用に起因して、このモデルの尤度面は、準最適な極大を有する。これらの極大のうちの１つに捕らわれることを回避し、それゆえにα_ｉに関する複数の準最適な推定値を得るために、モデルパラメータを適切に初期化する必要がある。このセクションにおける方策は、実施例４．１．１において１ｔｉｅ−Ｍｉｘ^２モデルに関して得た初期値を誤ったアノテーションに対して用い、それらをランダムにある量変化させることであった。このような方法で、２００個の異なる初期パラメータセットを生成し、それらに対して収束判定基準のうちの１つを満たすまで、すなわち、後続の反復間の対数尤度における差が０．５未満か、あるいはα_ｉおよびβ_ｊの間の差が０．００１未満となるまでＥＭアルゴリズムを行った。得られた２００個の結果から、最大尤度をもつ結果を推定値として選んだ。

このセクションにおいて調べた他のモデル、すなわち、Ｃｕｆｆｌｉｎｋｓおよび１ｔｉｅ−Ｍｉｘ^２モデルでは、誤ったアノテーションを補正できず、転写産物２は、そのシフトした開始に起因して失うよりも、その増加した長さによってより多くのフラグメントを獲得するため、推定したα_２が高過ぎることが予想できる。同様に、転写産物１は、転写産物３がその開始において失うよりも、その終結においてより多くのフラグメントを失うため、α_１がα_３よりさらに強く過小評価されることが予想できる。

４．２．１ α_１＝０．２８、α_２＝０．３２、α_３＝０．４に対する５ｔｉｅ−Ｍｉｘ^２モデルの収束
図２６は、重みα_１＝０．２８、α_２＝０．３２、α_３＝０．４に対して、ＥＭアルゴリズムの収束後に最大尤度をもつ初期パラメータセットに関するα_ｉの収束を示す。この初期パラメータセットではＥＭアルゴリズムが１４９回の反復後に収束した。図２６におけるｘ軸はＥＭアルゴリズムの反復を示し、ｙ軸は対応するαを示す。図２６における破線は、ＥＭアルゴリズムの経過中のα_ｉを示し、一方で水平一点鎖線は、真のα_ｉを示す。図２６は、α_ｉが真の値に非常に近い値へ収束することを示す。これは、表６にも反映され、同表は、推定したα_ｉとともにＥＭアルゴリズムの最終的な反復後の正しい重みからのそれらのＫＬダイバージェンスを示す。ＥＭアルゴリズムの間の他のモデルパラメータの値を図２７、２８、２９、３０および３１に示す。前述のように、β_ｊ、μ_ｊ、σ_ｊの取得値の品質は、結果として生じたｐ（ｒ｜ｔ＝ｉ）を検討することによってのみ評価できる。比較として、転写産物特有のシフトおよびスケール・パラメータν_ｉおよびλ_ｉは、最終的なモデルが誤った初期仮定を補正するかどうかについて指標を与える。ν_ｉの増加は、右へのシフトを示し、一方でλ_ｉの増加は、長さにおける増加を示す。図３０および３１は、それゆえに５ｔｉｅ−Ｍｉｘ^２モデルが転写産物１の開始を右へシフトさせて、その長さを増加させ、一方では転写産物２および３の開始を左へシフトさせて、その長さを減少させることを示す。これは、表５に示すような、誤ったアノテーションの正しいアノテーションからのずれと合致する。

ＥＭアルゴリズムの間の誤った転写産物アノテーションの補正は、ｐ^（ｎ）（ｒ｜ｔ＝ｉ）の収束を示す図３２、３３および３４においてもわかる。これらの図中の縦線は、問題の転写産物の転写産物座標における正しい開始および終結を示す。図３２は、転写産物１に関するｐ^（ｎ）（ｒ｜ｔ＝ｉ）の収束を示す。この図中の一点鎖線は、初期の転写産物特有のｐｄｆ ｐ^（０）（ｒ｜ｔ＝ｉ）を示し、一方で実線は、ｐ^{（１４９）}（ｒ｜ｔ＝ｉ）、１４９回の反復後のＥＭアルゴリズムの結果を示す。図３２は、図３０および３１と合致して、ｐ^（ｎ）（ｒ｜ｔ＝ｉ）がＥＭアルゴリズムの間にわずかに右へ動き、長さが増加することを示す。ｐ^{（１４９）}（ｒ｜ｔ＝ｉ）において生じた転写産物開始および転写産物長の値は、ほとんど完全である。図３２は、正しいアノテーションへの主要な一歩がＥＭアルゴリズムの最初の数回の反復において生じることも示す。これは、ν_ｉ ^（ｎ）およびλ_ｉ ^（ｎ）が２０回のＥＭ反復後に実質的に一定のままであることを示す図３０および３１と合致する。図３３および３４は、転写産物２および３に関するｐ^（ｎ）（ｒ｜ｔ＝ｉ）の収束を示す。これらの図が示すのは、ＥＭアルゴリズムが両方の転写産物を左へシフトさせ、それらの長さを減少させることによって、転写産物１に対するように、初期の誤った転写産物アノテーションを補正することを示す。最後に、図３５は、

の収束を示す。図３５における実線からわかるように、最終的な

は、図１３におけるヒストグラムの形状に対する良好な近似である。

４．２．２ 重みの６０個のセットに関する５ｔｉｅ−Ｍｉｘ^２、１ｔｉｅ−Ｍｉｘ^２およびＣｕｆｆｌｉｎｋｓモデル間の比較
このセクションは、重みα_１、α_２、α_３の６０個の異なるセット、ならびに５ｔｉｅ−Ｍｉｘ^２、１ｔｉｅ−Ｍｉｘ^２およびＣｕｆｆｌｉｎｋｓモデルを用いた実験を考察する。５ｔｉｅ−Ｍｉｘ^２モデルのパラメータを初期化するために、かつ１ｔｉｅ−Ｍｉｘ^２およびＣｕｆｆｌｉｎｋｓモデルのための固定参照枠として図２５における誤った転写産物アノテーションを用いた。図３６は、真の値α_３＝０．２に対して、５ｔｉｅ−Ｍｉｘ^２およびＣｕｆｆｌｉｎｋｓモデルに関して推定したα_ｉを示す。この図は、真のα_ｉと５ｔｉｅ−Ｍｉｘ^２モデルによって推定したα_ｉとの間の良好な一致を示し、一方でＣｕｆｆｌｉｎｋｓモデルによって推定したα_ｉは、正解から著しくずれる。真のα_１の小さい値に対してのみ、Ｃｕｆｆｌｉｎｋｓモデルは、５ｔｉｅ−Ｍｉｘ^２モデルに近づく推定値を作り出す。この場合には、転写産物１が非常に低い濃度を有し、それゆえに転写産物２および３からのフラグメントのみがＥＭアルゴリズムにおいて役割を果たす。転写産物１がその誤った末端に起因して転写産物２に対して失うフラグメントの数は無視できるため、これらのフラグメントは、ＥＭアルゴリズムによって転写産物２および３の間でほとんど正しく配分される。しかしながら、α_１の増加とともに、Ｃｕｆｆｌｉｎｋｓモデルは、α_２を著しく過大評価し、α_１およびα_３の過小評価がより顕著になる。この影響は、図３７ではさらに強く、短い初期の減少後に、α_２に関するＣｕｆｆｌｉｎｋｓの推定値は、０．５においてほとんど一定であるように見える。比較として、５ｔｉｅ−Ｍｉｘ^２モデルによって推定したα_ｉは、重みの全範囲にわたってやはり非常に正確である。

図３８および３９において、５ｔｉｅ−Ｍｉｘ^２モデルの推定値は、やはり非常に正確であり、一方でＣｕｆｆｌｉｎｋｓモデルは、α_１およびα_３を著しく過小評価する。加えて、図３８および３９ではＣｕｆｆｌｉｎｋｓモデルは、それぞれアルファ_２ ０．４および０．２８に関してほとんど一定の推定値をもたらす。表７からわかるように、完全な６０個のパラメータセットにわたってＣｕｆｆｌｉｎｋｓモデルは、平均ＫＬダイバージェンス０．２１９７７を得て、一方で５ｔｉｅ−Ｍｉｘ^２モデルは、平均ＫＬダイバージェンス０．０１４４８２を得る。このように、平均ＫＬダイバージェンスの観点から、５ｔｉｅ−Ｍｉｘ^２モデルの精度は、Ｃｕｆｆｌｉｎｋｓモデルの精度を１５倍上回る。事実、５ｔｉｅ−Ｍｉｘ^２モデルは、表３によれば、平均ＫＬダイバージェンス０．１２３６８をもたらす正しいアノテーションを用いたＣｕｆｆｌｉｎｋｓモデルよりさらにいっそう正確である。このように、平均ＫＬダイバージェンスの観点から、誤った転写産物アノテーションを用いた５ｔｉｅ−Ｍｉｘ^２モデルの精度は、正しいアノテーションを用いたＣｕｆｆｌｉｎｋｓモデルの精度を８倍上回る。

４．２．３ １ｐａｒ−Ｍｉｘ^２モデルとの比較
１ｔｉｅ−Ｍｉｘ^２モデルは、セクション４．１．２における正しい転写産物アノテーションを用いた実験ではＣｕｆｆｌｉｎｋｓモデルよりはるかに良好であったので、誤った転写産物アノテーションを両方が用いた５ｔｉｅ−Ｍｉｘ^２モデルと１ｔｉｅ−Ｍｉｘ^２モデルとを比較する価値がある。図４０は、真のα_３の値０．２に対して、１ｔｉｅ−Ｍｉｘ^２および５ｔｉｅ−Ｍｉｘ^２モデルに関して推定したα_ｉを示す。図３６におけるように、α_２は、固定した誤った転写産物アノテーションを用いる１ｔｉｅ−Ｍｉｘ^２モデルによって著しく過大評価される。図４０における１ｔｉｅ−Ｍｉｘ^２モデルのα_ｉの推定値は、図３６におけるＣｕｆｆｌｉｎｋｓモデルの推定値に比べてわずかにより正確であるが、５ｔｉｅ−Ｍｉｘ^２モデルより依然としてかなり不十分である。この傾向は、図４１、４２および４３にも呈示され、これらの図は、１ｔｉｅ−Ｍｉｘ^２モデルに関する推定値がＣｕｆｆｌｉｎｋｓモデルに関する推定値を少し改善するが、５ｔｉｅ−Ｍｉｘ^２モデルの推定値よりかなり不十分であることを示す。総じて、これは、表８における１ｔｉｅ−Ｍｉｘ^２モデルに関する平均ＫＬダイバージェンスをもたらし、同表は、ＫＬダイバージェンスの観点から、５ｔｉｅ−Ｍｉｘ^２モデルの精度が１ｔｉｅ−Ｍｉｘ^２モデルの精度より約１２倍高いことを示す。要約すれば、このセクションにおける実験は、正しい転写産物アノテーションがないときにα_ｉに関する信頼性の高い推定値を得るためには、正しい転写産物位置を知ることが可能なモデルの使用が重要であることを示す。

実施例５：結論
先のセクションにおける実験は、正しいおよび誤った転写産物アノテーションの両方を用いて１ｔｉｅ−Ｍｉｘ^２および５ｔｉｅ−Ｍｉｘ^２モデルをＣｕｆｆｌｉｎｋｓモデルと比較した。６０個のパラメータセットα_１，α_２，α_３のセットに関する実験の結果を表９にまとめる。これらの結果は、正しい転写産物アノテーションを用いると１ｔｉｅ−Ｍｉｘ^２モデルは、Ｃｕｆｆｌｉｎｋｓモデルよりはるかに優れており、一方で誤った転写産物アノテーションを用いると５ｔｉｅ−Ｍｉｘ^２モデルが１ｔｉｅ−Ｍｉｘ^２およびＣｕｆｆｌｉｎｋｓモデルの両方に勝ることを示す。これは、適切なＭｉｘ^２モデルの使用が存在比の推定値の精度をかなり改善することを示唆する。

実施例６：先行技術との比較
Ｃｕｆｆｌｉｎｋｓは、転写産物アセンブリおよび転写産物量の推定のための方法を実装したプログラムである。実装の詳細な記載は、非特許文献１に見ることができる。Ｃｕｆｆｌｉｎｋｓでは、本発明によるα_ｉに対応する転写産物量を式（１３）における因数分解を通じて推定し、式（１３）においてｐ_{ｔｒａｎｓ}（ｓ（ｒ）｜ｔ＝ｉ，ｌ（ｒ））は均一であり、ｐ_{ｔｒａｎｓ}（ｌ（ｒ）｜ｔ＝ｉ）は、転写産物ｔ＝ｉに依存しない。Ｍｉｘ^２モデルとは対照的に、Ｃｕｆｆｌｉｎｋｓは、フラグメント開始点ｓ（ｒ）の分布をデータから知ることがなく、そのうえ、正しい転写産物アノテーションの利用可能性に依存する。

Ｃｕｆｆｌｉｎｋｓの存在比推定への拡張を実装するロバーツら、２０１１、では、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）が式（１３）におけるように因数分解される。ｐ_{ｔｒａｎｓ}（ｓ（ｒ）｜ｔ＝ｉ，ｌ（ｒ））は、配列特有の位置的な重みの正規化された積であり、これらの重みが転写産物座標における位置ごとに定義される。この結果、単一の転写産物に関してさえトレーニングが必要な多数のパラメータが生じる。参考文献９の実験ではモデルを計算上扱いやすくするために位置的な重みをステップ関数に制限しなければならず、そのうちの５つのみが推定されて、すべての転写産物間で共有される。参考文献９とは対照的に、本発明のモデルは、すべての確率ｐ_{ｔｒａｎｓ}（ｓ（ｒ）｜ｔ＝ｉ）を個別には推定しないが、混合確率分布（１６）のパラメータを推定する。結果として、本発明のモデルを推定するために必要なパラメータは、はるかにより少なく、それゆえに、そのパラメータ推定が参考文献９におけるモデルよりロバストであり、計算上扱いやすい。加えて、参考文献９におけるモデルは、転写産物量に関する正確な開始推定値を必要とし、それゆえに、その重みが単一のアイソフォーム遺伝子上で推定される。従って、参考文献９におけるモデルは、本明細書に示される実験に用いるデータには適用できない。これに対して、本発明のモデルは、転写産物量に関する正確な開始推定値を必要とせず、それゆえに、複数のアイソフォーム遺伝子上でトレーニングすることができる。加えて、参考文献９におけるモデルは、本発明のモデルと異なり、正しい転写産物アノテーションを必要とする。

ウーら、２０１１、では、各エクソンおよび転写産物に単一の重みを割り当てることによって、フラグメントの分布におけるバイアスが遺伝子のエクソンに関してモデリングされる。本発明のモデルとは対照的に、これは、ｐ_{ｔｒａｎｓ}（ｓ（ｒ）｜ｔ＝ｉ）が各エクソン上で一定であるように制限する。加えて、ウーら、２０１１、における重みは、部分的に確率論的枠外の発見的方法によって推定される。ウーら、２０１１、における研究が本発明のモデルとさらに異なるのは、ウーら、２０１１、における研究がエクソン上のリード・カウントの確率にポアソン分布を用いるのに対して、本発明のモデルが、転写産物を所与として、フラグメントの確率に混合モデルを用いる点である。加えて、ウーら、２０１１、における研究は、正しい転写産物アノテーションに依存する。

リーら、２０１０、は、配列特有のバイアスを補正し、それゆえに、参考文献９における配列特有の重みと同様である。リーら、２０１０、は、カウント・データをモデリングするためにポアソン分布を用いる点ではウーら、２０１１、における研究にも類似する。ウーら、２０１１、とは対照的に、リーら、２０１０、におけるカウント・データは、エクソンには基づかず、カバレッジによって与えられる。リーら、２０１０、における研究は、転写産物量を推定しないという点で、本発明のモデルならびに参考文献９およびウーら、２０１１、におけるモデルとは異なる。

グラウスら、２０１２、では、転写産物によって生成されたリードを観測する確率に関してベイジアン・モデルが用いられる。このモデルは、転写産物量を推定するために参考文献９と同じ配列特有の位置的な重みおよびギブスサンプリング手順を用いる。このように、本発明のモデルと比較して、グラウスら、２０１２、におけるパラメータ推定は、最大尤度の枠組みではなくベイジアン内で行われる。グラウスら、２０１２、におけるモデルは、ｐ（ｒ｜ｔ＝ｉ）に混合確率分布を何も用いず、正しい転写産物アノテーションを必要とするという点においても本発明のモデルとは異なる。

リーら、２０１０および２０１１、におけるモデルは、グラウスら、２０１２、におけるモデルと同様の構造をもつベイジアン・モデルである。しかしながら、グラウスら、２０１２、とは対照的に、リーら、２０１０および２０１１、における転写産物量は、最大尤度の枠内で、特にＥＭアルゴリズムを用いてトレーニングされる。ｐ（ｓ（ｒ）｜ｔ＝ｉ）のような、残りのモデルパラメータは、事前に導出されるか、または発見的方法を用いて推定されるかいずれかである。これは、すべてのそのパラメータを最大尤度の枠内で推定する本発明のモデルとは対照的である。加えて、リーら、２０１０および２０１１、におけるモデルは、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）をモデリングするために混合確率分布を用いず、正しい転写産物アノテーションを必要とする。

本発明に従って、本発明のモデルを推定するために期待値最大化アルゴリズムを用いた。期待値最大化アルゴリズムの一般的な枠組みは、デンプスターら、１９７７、において開発され、一方で本発明のモデルのための具体的なＥＭ更新式は、本明細書に記載するように導出した。

実施例７：Ｍｉｘ ^２ モデルを用いた転写産物群の連結
式（２０）に続く段落において、Ｍｉｘ^２モデルのパラメータの連結は、異なる転写産物ｔ＝ｉに関するｐ（ｒ｜ｔ＝ｉ）間のある類似性を示唆することを述べた。それゆえに、Ｍｉｘ^２モデルのパラメータは、この類似性を呈示する転写産物ｔ＝ｉ間でのみ連結されるべきである。群内の転写産物のみがそれらのパラメータを共有するような異なる群に転写産物が分けられるならば、Ｍｉｘ^２モデルのＥＭ更新式（２７）、（３６）および（３７）は、修正される必要がある。以下では、各転写産物ｔ＝ｉは、関数Ｇ（ｉ）＝ｋを通じてリトリーブできる関連付けられた群ｇ＝ｋを有する。その場合、群ｇ＝ｋ内のパラメータβ_ｋ，ｊ，μ_ｋ，ｊおよびσ_ｋ，ｊのＥＭ更新式は、以下のように与えられる。

制約条件

を実行するためにラグランジェ法を用いて、β_ｋ，ｊに関する微分をとると、

を結果として生じ、ここで

である。転写産物間で連結される残りのパラメータ、すなわち、μ_ｋ，ｊおよびσ_ｋ，ｊに関しては、転写産物ｔ＝１，．．．，Ｎの完全なセットにわたる和を群ｇ＝ｋ内の転写産物にわたる和で置き換えることによって更新式（３６）および（３７）を修正する必要があり、すなわち、

である。遺伝子座における各転写産物をそれ自体の群へ配置すると、すべてのパラメータ連結が解除されて、前述のように、フラグメント分布および転写物量に関する不正確な推定値に繋がる。転写産物を異なる群へ配置するときには、それゆえに、同様でないフラグメント分布ｐ（ｒ｜ｔ＝ｉ）の分離と、Ｍｉｘ^２モデルの安定性を確保するために十分な数のパラメータの連結の維持との間で適正なバランスをとることが重要である。賢明な必要条件は、例として、各群が少なくとも２つの転写産物を含むことである。代わりに、３つ以上の群が存在するときには、多くとも１つの群が単一の転写産物を含むことを必要とするのが賢明である。

実施例８：Ｍｉｘ ^２ モデルを用いたフラグメント長分布のトレーニング
セクション４において考察した実験では、転写産物特有のフラグメント確率ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ，ｂ＝ｊ）を（４８）におけるように因数分解し、フラグメント長の確率分布ｐ_{ｔｒａｎｓ}（ｌ（ｒ）｜ｔ＝ｉ，ｂ＝ｊ，ｓ（ｒ））は、フラグメント開始ｓ（ｒ）および転写産物ｔ＝ｉの長さｌ（ｔ＝ｉ）のみに依存すると仮定した。加えて、ｐ_{ｔｒａｎｓ}（ｌ（ｒ）｜ｌ（ｔ＝ｉ），ｓ（ｒ））は、所与であると仮定した。これらの仮定は、ｐ_{ｔｒａｎｓ}（ｌ（ｒ）｜ｔ＝ｉ，ｂ＝ｊ，ｓ（ｒ））をデータセットＲから推定するのであれば必要なく、Ｍｉｘ^２モデルの枠内でこの推定を行うことができる。このために、フラグメント開始ｓ（ｒ）の分布については、ｐ_{ｔｒａｎｓ}（ｌ（ｒ）｜ｔ＝ｉ，ｂ＝ｊ，ｓ（ｒ））を混合確率分布、すなわち、

として書き、ｂｓ＝ｊは、先にｂ＝ｊによって示した隠れ変数である。ここではｂｓは、「フラグメント開始のビルディングブロック」に対する簡略記憶記号であり、一方で隠れ変数ｂｌは、「フラグメント長のビルディングブロック」に対する簡略記憶記号である。ｐ_{ｔｒａｎｓ}（ｌ（ｒ）｜ｔ＝ｉ，ｂｓ＝ｊ，ｓ（ｒ））は、ｂｓ＝ｊに依存しないと仮定するのが賢明であり、従って（５７）は、

に帰着する。（１５）、（４８）および（５８）を組み合わせると、フラグメントの確率に関する次の表現をもたらし、

ここで以下が成り立つ。

従って、（５９）は、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）の混合重みがβ_ｊとγ_ｋとの積であるＭｉｘ^２モデルであり、ｐ_{ｔｒａｎｓ}（ｒ｜ｔ＝ｉ）の混合成分は、（６０）における条件付き確率分布の積である。

（１５）および（２５）に加えて、（５９）から以下を導出できる。

結果として、（２２）および（２７）と同様に、ＥＭアルゴリズムを用いてγ_ｋを次のように推定でき、

である。フラグメント長のビルディングブロックｐ_{ｔｒａｎｓ}（ｌ（ｒ）｜ｔ＝ｉ，ｂｌ＝ｋ，ｓ（ｒ））に関しては、例として、その平均値がｓ（ｒ）とｌ（ｔ＝ｉ）との間、または１とｌ（ｔ＝ｉ）との間のいずれかに等距離に分布して、離散的または連続的な１次元確率空間内のいずれかにおいて正規化したガウシアンを用いることができる。連続的な１次元確率空間上の分布を選ぶならば、それらの内部パラメータ、例えば、平均値、標準偏差、シフトおよびスケール・パラメータは、フラグメント開始ｐ_{ｔｒａｎｓ}（ｓ（ｒ）｜ｔ＝ｉ）の確率分布に関するこれらのパラメータと同様に推定できる。

実施例９：固定した転写産物エンドポイントを用いたν _ｉ の評価
５ｔｉｅ−Ｍｉｘ^２モデルは、シフトおよびスケール・パラメータν_ｉおよびλ_ｉを推定する。ν_ｉを０にセットしてＥＭアルゴリズムの間に更新しなければ、転写産物ｔ＝ｉの開始は、変化しないままである。転写産物のエンドポイントを固定することになれば、５ｔｉｅ−Ｍｉｘ^２モデルを少し修正する必要がある。

Ｌ_ｉがアノテーションによる転写産物ｔ＝ｉの長さであれば、転写産物ｔ＝ｉの固定したエンドポイントに関して、以下が成り立つ。

これは、ν_ｉに関する次のＥＭ更新式をもたらす。

これから、λ_ｉ ^{（ｎ＋１）}を（６５）に従って導出でき、すなわち、

である。

実施例１０：Ｍｉｘ ^２ モデルの他のバイアス・モデルとの組み合わせ
第４のセクションの実験において考察したＭｉｘ^２パラメータの連結は、位置的なフラグメンテーション・バイアス、すなわち、転写産物内のフラグメント開始に関係するバイアスに関するモデルを実装する。配列特有のバイアスのような、他の種類のバイアスを他のモデル、例えば、非特許文献２における可変長隠れマルコフ・モデル（ＶＬＭＭ：ｖａｒｉａｂｌｅ ｌｅｎｇｔｈ ｈｉｄｄｅｎ Ｍａｒｋｏｖ ｍｏｄｅｌ）を用いて記述してもよい。典型的に、非位置的なバイアスに関するモデルは、ヌクレオチド配列の観測度数をバイアスされていないデータの帰無仮説下におけるそれらの度数と比較する。バイアスされたデータ中でｒの単一のコピーを観測することを考慮すれば、バイアスが何もないときのフラグメントｒの多重度ｍ＝ｃにわたる確率分布ｐ（ｍ＝ｃ｜ｒ）を導出するためにこの比較を用いることができる。次に、分布ｐ（ｍ＝ｃ｜ｒ）は、Ｍｉｘ^２モデルのＥＭ更新式における各フラグメントｒを予想されるｒの多重度により重み付けすることによって、非位置的なバイアスをデータから計算上除去するために用いることができる。存在比α_ｉのＥＭ更新式に関して、例として、これは、

をもたらす。留意すべきは、ｐ（ｍ＝ｃ｜ｒ）の期待値がそれにより無限大になりかねないため、Ｍｉｘ^２モデルの下でデータセットの尤度Ｒを最大化することによって分布ｐ（ｍ＝ｃ｜ｒ）を推定することができないことである。Ｍｉｘ^２モデルの他のバイアス・モデルとの組み合わせについては、このセクションにおいて考察したように、ｐ（ｍ＝ｃ｜ｒ）をＭｉｘ^２モデルの最大尤度の枠外で推定することがそれゆえに重要である。

実施例１１：ＦＰＫＭおよびＲＰＫＭ値を補正するためのＭｉｘ ^２ モデルの使用
ＲＮＡ−Ｓｅｑにおける転写産物の濃度は、通常、ＦＰＫＭ（１００万リード当たり、１０００塩基対当たりのフラグメント数）またはＲＰＫＭ（１００万リード当たり、１０００塩基対当たりのリード数）尺度を用いて測定され、転写産物ｔ＝ｉに関して、後者は、

によって与えられる。ＦＰＫＭ尺度では、（６９）における転写産物ｔ＝ｉの長さｌ（ｔ＝ｉ）が調整した転写産物長

によって置き換えられる（非特許文献１）。このように、正確なＦＰＫＭおよびＲＰＫＭ値を得るためには（６９）において正確な転写産物長ｌ（ｔ＝ｉ）を用いることが重要である。５ｔｉｅ−Ｍｉｘ^２モデルは、不正確な転写産物アノテーションを補正でき、それゆえにより高い精度の転写産物長をもたらし、改善された精度のＲＰＫＭおよびＦＰＫＭ値を結果として生じる。５ｔｉｅ−Ｍｉｘ^２モデルに基づいて転写産物長に関する推定値を得るために種々の方法を利用できる。転写産物開始は、例として、

ここで０＜ｃ_ｓ＜１である位置ｘとして推定できる。典型的に、ｃ_ｓは、ｃ_ｓ＝０．０５のような小さい正値となるであろう。同様に、転写産物終結は、ｃ_ｓを置き換えた１に近いであろう０＜ｃ_ｅ＜１のｃ_ｅ値、例えばｃ_ｅ＝０．９５によって（７０）が成り立つ位置として推定できる。同様に、（７０）は、ｃ_ｓおよびｃ_ｅとともに、任意の混合成分ｐ（ｒ｜ｔ＝ｉ，ｂ＝ｊ）に適用できる。この場合、転写産物開始を確定するために用いる混合成分は、典型的に、小さいｓ（ｒ）の領域に集中するのに対して、転写産物終結を確定するために用いる混合成分は、大きいｓ（ｒ）の領域に集中するであろう。５ｔｉｅ−Ｍｉｘ^２モデルから転写産物ｔ＝ｉの長さに関する推定値を得るための別の方法は、転写産物アノテーションの長さにＭｉｘ^２モデルの収束後のλ_ｉを乗じる。

参考文献
［１］Ａ．Ｐ．デンプスター（Ｄｅｍｐｓｔｅｒ）、Ｎ．Ｍ．レアード（Ｌａｉｒｄ）、およびＤ．Ｂ．ルービン（Ｒｕｂｉｎ）ｅｍアルゴリズムを用いた不完全データからの最大尤度（Ｍａｘｉｍｕｍ ｌｉｋｅｌｉｈｏｏｄ ｆｒｏｍ ｉｎｃｏｍｐｌｅｔｅ ｄａｔａ ｖｉａ ｔｈｅ ｅｍ ａｌｇｏｒｉｔｈｍ）Ｊｏｕｒｎａｌ ｏｆ ｔｈｅ Ｒｏｙａｌ Ｓｔａｔｉｓｔｉｃａｌ Ｓｏｃｉｅｔｙ、Ｓｅｒｉｅｓ Ｂ、３９（１）：１ ３８、１９７７
［２］ピーター・グラウス（Ｐｅｔｅｒ Ｇｌａｕｓ）、アンティ・ホンケラ（Ａｎｔｔｉ Ｈｏｎｋｅｌａ）およびマグナス・ラトリ（Ｍａｇｎｕｓ Ｒａｔｔｒａｙ）生物学的多様性をもつＲＮＡ−Ｓｅｑデータからの発現変動転写産物の同定（Ｉｄｅｎｔｉｆｙｉｎｇ ｄｉｆｆｅｒｅｎｔｉａｌｌｙ ｅｘｐｒｅｓｓｅｄ ｔｒａｎｓｃｒｉｐｔｓ ｆｒｏｍ ＲＮＡ−Ｓｅｑ ｄａｔａ ｗｉｔｈ ｂｉｏｌｏｇｉｃａｌ ｖａｒｉａｔｉｏｎ）Ｂｉｏｉｎｆｏｒｍａｔｉｃｓ ２８（１３）：１７２１−１７２８、２０１２
［３］ボー・リー（Ｂｏ Ｌｉ）およびコリン・デューイ（Ｃｏｌｉｎ Ｄｅｗｅｙ）ＲＳＥＭ：リファレンスゲノム有無のｒｎａ−ｓｅｑデータからの正確な転写産物定量化（ａｃｃｕｒａｔｅ ｔｒａｎｓｃｒｉｐｔ ｑｕａｎｔｉｆｉｃａｔｉｏｎ ｆｒｏｍ ｒｎａ−ｓｅｑ ｄａｔａ ｗｉｔｈ ｏｒ ｗｉｔｈｏｕｔ ａ ｒｅｆｅｒｅｎｃｅ ｇｅｎｏｍｅ）ＢＭＣ Ｂｉｏｉｎｆｏｒｍａｔｉｃｓ １２（１）：３２３、２０１１
［４］ボー・リー（Ｂｏ Ｌｉ），ビクター・ルオーティ（Ｖｉｃｔｏｒ Ｒｕｏｔｔｉ）、ロン・Ｍ・スチュワート（Ｒｏｎ Ｍ Ｓｔｅｗａｒｔ）、ジェームズ・Ａ・トムソン（Ｊａｍｅｓ Ａ Ｔｈｏｍｓｏｎ）、コリン・Ｎ・デューイ（Ｃｏｌｉｎ Ｎ Ｄｅｗｅｙ）リード・マッピングの不確実性を伴うＲｎａ−ｓｅｑ遺伝子発現推定（Ｒｎａ−ｓｅｑ ｇｅｎｅ ｅｘｐｒｅｓｓｉｏｎ ｅｓｔｉｍａｔｉｏｎ ｗｉｔｈ ｒｅａｄ ｍａｐｐｉｎｇ ｕｎｃｅｒｔａｉｎｔｙ）Ｂｉｏｉｎｆｏｒｍａｔｉｃｓ ２６（４）：４９３−５００、２０１０年２月
［５］ジュン・リー（Ｊｕｎ Ｌｉ）、フェイ・ジアン（Ｈｕｉ Ｊｉａｎｇ）およびウィン・ワン（Ｗｉｎｇ Ｗｏｎｇ）ＲＮＡ−Ｓｅｑデータにおけるショートリードの割合の不均一性モデリング（Ｍｏｄｅｌｉｎｇ ｎｏｎ−ｕｎｉｆｏｒｍｉｔｙ ｉｎ ｓｈｏｒｔ−ｒｅａｄ ｒａｔｅｓ ｉｎ ＲＮＡ ｓｅｑ ｄａｔａ）Ｇｅｎｏｍｅ Ｂｉｏｌｏｇｙ １１（５）：Ｒ５０＋、２０１０年
［６］アダム・ロバーツ（Ａｄａｍ Ｒｏｂｅｒｔｓ）、コール・トラップネル（Ｃｏｌｅ Ｔｒａｐｎｅｌｌ）、ジュリー・ ドナヒー（Ｊｕｌｉｅ Ｄｏｎａｇｈｅｙ）、ジョン・Ｌ・リン（Ｊｏｈｎ Ｌ Ｒｉｎｎ）およびリオル・パクター（Ｌｉｏｒ Ｐａｃｈｔｅｒ）フラグメント・バイアス補正によるｒｎａ−ｓｅｑ発現推定値の改善（Ｉｍｐｒｏｖｉｎｇ ｒｎａ−ｓｅｑ ｅｘｐｒｅｓｓｉｏｎ ｅｓｔｉｍａｔｅｓ ｂｙ ｃｏｒｒｅｃｔｉｎｇ ｆｏｒ ｆｒａｇｍｅｎｔ ｂｉａｓ）Ｇｅｎｏｍｅ Ｂｉｏｌ １２（３）：Ｒ２２ ２０１１年３月
［７］コール・トラップネル（Ｃｏｌｅ Ｔｒａｐｎｅｌｌ）、ブライアン・Ａ・ウィリアムズ（Ｂｒｉａｎ Ａ Ｗｉｌｌｉａｍｓ）、ジオ・ペルテア（Ｇｅｏ Ｐｅｒｔｅａ）、アリ・モルタザヴィ（Ａｌｉ Ｍｏｒｔａｚａｖｉ）、ゴードン・クワン（Ｇｏｒｄｏｎ Ｋｗａｎ）、マレイケ・Ｊ・ヴァン・バレン（Ｍａｒｉｊｋｅ Ｊ ｖａｎ Ｂａｒｅｎ）、スチーブン・Ｌ・ザルツバーグ（Ｓｔｅｖｅｎ Ｌ Ｓａｌｚｂｅｒｇ）、バーバラ・Ｊ・ウォルド（Ｂａｒｂａｒａ Ｊ Ｗｏｌｄ）およびリオル・パクター（Ｌｉｏｒ Ｐａｃｈｔｅｒ）転写産物アセンブリおよびＲＮＡ−Ｓｅｑによる定量化が細胞分化の間のアノテーションされていない転写産物およびアイソフォームのスイッチングを明らかにする（Ｔｒａｎｓｃｒｉｐｔ ａｓｓｅｍｂｌｙ ａｎｄ ｑｕａｎｔｉｆｉｃａｔｉｏｎ ｂｙ ＲＮＡ−Ｓｅｑ ｒｅｖｅａｌｓ ｕｎａｎｎｏｔａｔｅｄ ｔｒａｎｓｃｒｉｐｔｓ ａｎｄ ｉｓｏｆｏｒｍ ｓｗｉｔｃｈｉｎｇ ｄｕｒｉｎｇ ｃｅｌｌ ｄｉｆｆｅｒｅｎｔｉａｔｉｏｎ）Ｎａｔ Ｂｉｏｔｅｃｈｎｏｌ ２８（５）：５１１−５１５、２０１０年５月
［８］ジュヨンプオン・ウー（Ｚｈｅｎｇｐｅｎｇ Ｗｕ）、シー・ワーン（Ｘｉ Ｗａｎｇ）およびシュエゴン・ジャーン（Ｘｕｅｇｏｎｇ Ｚｈａｎｇ）、ＲＮＡ−Ｓｅｑにおけるアイソフォーム発現推測を改善するための不均一リード分布モデルの使用（Ｕｓｉｎｇ ｎｏｎ−ｕｎｉｆｏｒｍ ｒｅａｄ ｄｉｓｔｒｉｂｕｔｉｏｎ ｍｏｄｅｌｓ ｔｏ ｉｍｐｒｏｖｅ ｉｓｏｆｏｒｍ ｅｘｐｒｅｓｓｉｏｎ ｉｎｆｅｒｅｎｃｅ ｉｎ ＲＮＡ−Ｓｅｑ） Ｂｉｏｉｎｆｏｒｍａｔｉｃｓ ２７（４）：５０２−５０８、２０１１年２月
［９］ロバーツ（Ｒｏｂｅｒｔｓ）ら：ゲノム生物学（Ｇｅｎｏｍｅ Ｂｉｏｌｏｇｙ）１２（３）（２０１１）：Ｒ２２
［１０］ウェン−ピン（Ｗｅｎ−Ｐｉｎｇ）ら：ゲノム生物学（Ｇｅｎｏｍｅ Ｂｉｏｌｏｇｙ）８（６）（２００７）：Ｒ９８

Claims (15)

1. コンピュータを用いて転写産物量を推定する方法であって、
   ａ）対象となる遺伝子座の転写産物の潜在的な混合物から転写産物フラグメントシークエンシングデータを得るステップと、
   ｂ）前記フラグメントシークエンシングデータを対象となる前記遺伝子座の遺伝子座標に割り当て、それによってフラグメント遺伝子座標カバレッジのデータセットを得るステップであって、遺伝子座標ごとの前記カバレッジは、結合されてカバレッジ包絡曲線を形成する、ステップと、
   ｃ）前記混合物の転写産物の数をセットするステップと、
   ｄ）転写産物ｉごとにモデリングされた遺伝子カバレッジの確率分布関数を予めセットするステップであって、ｉは、転写産物のための数値識別子を示し、前記確率分布関数は、少なくとも２つの確率サブ関数ｊの和を乗じた前記転写産物ｉの重み係数α_ｉによって定義され、ｊは、確率サブ関数のための数値識別子を示し、各確率サブ関数ｊは、重み係数β_ｉ，ｊによって独立に重み付けされる、ステップと、
   ｅ）サム関数を得るために各転写産物の前記確率分布関数を加算するステップと、
   ｆ）前記サム関数を前記カバレッジ包絡曲線へフィッティングし、それによってフィットを向上させるために、α_ｉおよびβ_ｉ，ｊに関する値を最適化するステップと、
   ｇ）予めセットされた収束判定基準が満たされるまでステップｅ）およびｆ）を繰り返し、それによって前記収束判定基準が満たされた後に最適化されるような前記重み係数α_ｉによって与えられる、前記混合物の転写産物ごとの推定転写産物量を得るステップと
   を備える方法であって、
   前記コンピュータが、前記ステップａ）における転写産物フラグメントシークエンシングデータを入力する手段、前記ステップｂ）〜ｇ）を行う演算手段、及び前記ステップｇ）における前記混合物の転写産物ごとの推定転写産物量を出力する手段を備える、方法。
2. 転写産物フラグメントシークエンシングデータは、少なくとも５つの転写産物フラグメント配列を備える、請求項１に記載の方法。
3. 対象となる前記遺伝子座は、１つ以上の遺伝子または遺伝因子の１つ以上のアイソフォームを備え、好ましくは１つの遺伝子または遺伝因子の少なくとも２つのスプライスバリアントを備える、請求項１または２に記載の方法。
4. 転写産物の数をセットする前記ステップは、予めアノテーションされている配列データを対象となる前記遺伝子座から得ることと、転写産物の前記数を、対象となる前記遺伝子座から予想される、異なるアイソフォームとしてカウントするスプライスバリアントを含む、異なるアイソフォームの少なくとも前記数にセットすることとを備える、請求項１から３のいずれか一項に記載の方法。
5. 前記確率サブ関数ｊは、各遺伝子座標ごとに正値から構成され、好ましくは密度関数である、請求項１から４のいずれか一項に記載の方法。
6. 前記確率サブ関数ｊは、非周期関数であり、好ましくはガウス関数、正方形関数、三角形関数、特に好ましくはガウス関数である、請求項１から５のいずれか一項に記載の方法。
7. 前記遺伝子座標は、対象とならない遺伝子領域を削除するように随意的に変換された、ゲノムにおけるヌクレオチド位置に対応し、好ましくは前記対象とならない遺伝子領域は、前記転写産物フラグメントシークエンシングデータによるカバレッジを含まない、請求項１から６のいずれか一項に記載の方法。
8. スプライスジャンクションをもつ遺伝子座標位置を前記カバレッジ包絡曲線から除去することを備えたステップｂ２）をさらに備える、請求項１から７のいずれか一項に記載の方法。
9. 前記フラグメント遺伝子座標カバレッジは、遺伝子座標に割り当てられたフラグメント配列ごとに少なくとも１つのヌクレオチドのカウントを含み、好ましくは、前記少なくとも１つのヌクレオチドは、フラグメント開始点またはフラグメント配列全体を備える、請求項１から８のいずれか一項に記載の方法。
10. 転写産物に関する前記確率サブ関数は、それぞれ異なる遺伝子座標における極大を備える、請求項１から９のいずれか一項に記載の方法。
11. ステップｄ）において、転写産物に関する前記確率サブ関数は、正値を用いて転写産物の全長をカバーするために前記遺伝子座標において配置またはシフトされる、請求項１から１０のいずれか一項に記載の方法。
12. 少なくとも１つの転写産物、好ましくはｍＲＮＡの配列リードを確定するステップを備え、前記リードは、前記転写産物フラグメントシークエンシングデータを提供するために前記転写産物のフラグメントの配列を備える、請求項１から１１のいずれか一項に記載の方法。
13. 前記転写産物フラグメントシークエンシングデータの前記転写産物フラグメント配列は、５から８００ヌクレオチド、好ましくは６から６００ヌクレオチド、より好ましくは７から４００ヌクレオチド、なおさらに好ましくは８から２００ヌクレオチド、とりわけ好ましくは９から１５０ヌクレオチド、特に好ましくは１０から１００ヌクレオチド、最も好ましくは１２から７０ヌクレオチドの長さを有する、請求項１から１２のいずれか一項に記載の方法。
14. 転写産物ｉに関する確率サブ関数ごとの半値全幅値は、およそ同一である、請求項１から１３のいずれか一項に記載の方法。
15. コンピュータ上で請求項１から１４のいずれか一項の方法を行うためのコンピュータプログラム製品を備える、コンピュータ可読メモリ装置。
    

Images