TOP > クイック検索 > 国内特許検索 > トピックモデルを用いた遺伝子情報推定装置 > 明細書

明細書 :トピックモデルを用いた遺伝子情報推定装置

発行国 日本国特許庁(JP)
公報種別 公開特許公報(A)
公開番号 特開2018-139043 (P2018-139043A)
公開日 平成30年9月6日(2018.9.6)
発明の名称または考案の名称 トピックモデルを用いた遺伝子情報推定装置
国際特許分類 G06F  19/24        (2011.01)
C12N  15/09        (2006.01)
C12Q   1/68        (2018.01)
FI G06F 19/24
C12N 15/00 A
C12Q 1/68 Z
請求項の数または発明の数 8
出願形態 OL
全頁数 17
出願番号 特願2017-033381 (P2017-033381)
出願日 平成29年2月24日(2017.2.24)
発明者または考案者 【氏名】岩山 幸治
【氏名】永野 惇
出願人 【識別番号】597065329
【氏名又は名称】学校法人 龍谷大学
個別代理人の代理人 【識別番号】110000914、【氏名又は名称】特許業務法人 安富国際特許事務所
審査請求 未請求
テーマコード 4B063
Fターム 4B063QA01
4B063QA13
4B063QQ42
4B063QS39
要約 【課題】低コストかつ高精度で遺伝子情報を推定できる装置を提供する。
【解決手段】既存のトピックモデルに基づき、遺伝子情報に適した新しいトピックモデルを用いた遺伝子情報を推定する。このモデルでは単語の生成分布を従来の多項分布から負の二項分布に置き換える。
【選択図】図2
特許請求の範囲 【請求項1】
(i)複数のサンプルに対応する遺伝子情報の測定データを読み込む手順、
(ii)遺伝子情報の測定データに基づき遺伝子の平均出現頻度を算出する手順、
(iii)各サンプルの各遺伝子にトピックが割り当てられる確率を初期化する手順、
(iv)遺伝子へのトピック割り当て確率を更新する手順、
(v)トピック情報を更新する手順、
(vi)トピック情報の分散を更新する手順、
(vii)負の二項分布のdispersionパラメータを更新する手順、
および、
(viii)トピック情報およびトピック出現頻度から遺伝子情報を推定する手順、
を有する、
トピックモデルを用いた遺伝子情報推定装置。
【請求項2】
手順(iv)、(v)、および(vii)において、過分散の程度を表すパラメータをφとし、サンプル(s)においてトピック(k)を割り当てられた遺伝子(i)の期待出現数をμskiとしたとき、サンプル(s)における遺伝子(i)の出現数(rsi)が、
【数1】
JP2018139043A_000040t.gif
という確率に従う、請求項1に記載の装置。
【請求項3】
サンプル(s)における全遺伝子のカウントの総和(ν)、全サンプルでの全遺伝子の出現数の総和に対する遺伝子(i)の出現数の割合の対数(m)、指数分布あるいはジェフリーズ事前分布に従う分散(τki)、および平均が0で分散がτkiの正規分布に従うトピック(k)における遺伝子(i)の期待カウントの平均からのずれを表す量が与えられたときに、サンプル(s)においてトピック(k)を割り当てられた遺伝子(i)の期待出現数(μski)が
μski=νexp(ηki+m
と表される、請求項2に記載の装置。
【請求項4】
請求項1~3のいずれか一項に記載の遺伝子情報推定装置を構成する各手段としてコンピュータを機能させるための、遺伝子情報推定プログラム。
【請求項5】
請求項4に記載の遺伝子情報推定プログラムが記録されたことを特徴とするコンピュータ読み取り可能な記録媒体。
【請求項6】
(i)複数のサンプルに対応する遺伝子情報の測定データを読み込むステップ、
(ii)遺伝子情報の測定データに基づき遺伝子の平均出現頻度を算出するステップ、
(iii)各サンプルの各遺伝子にトピックが割り当てられる確率を初期化するステップ、
(iv)遺伝子へのトピック割り当て確率を更新するステップ、
(v)トピック情報を更新するステップ、
(vi)トピック情報の分散を更新するステップ、
(vii)負の二項分布のdispersionパラメータを更新するステップ、
および、
(viii)トピック情報およびトピック出現頻度から遺伝子情報を推定するステップ、
を有する、トピックモデルを用いた遺伝子情報推定方法。
【請求項7】
ステップ(iv)、(v)、および(vii)において、過分散の程度を表すパラメータをφとし、サンプル(s)においてトピック(k)を割り当てられた遺伝子(i)の期待出現数をμskiとしたとき、サンプル(s)における遺伝子(i)の出現数(rsi)が、
【数2】
JP2018139043A_000041t.gif
という確率に従う、請求項6に記載の方法。
【請求項8】
サンプル(s)における全遺伝子のカウントの総和(ν)、全サンプルでの全遺伝子の出現数の総和に対する遺伝子(i)の出現数の割合の対数(m)、指数分布あるいはジェフリーズ事前分布に従う分散(τki)、および平均が0で分散がτkiの正規分布に従うトピック(k)における遺伝子(i)の期待カウントの平均からのずれを表す量が与えられたときに、サンプル(s)においてトピック(k)を割り当てられた遺伝子(i)の期待出現数(μski)が
μski=νexp(ηki+m
と表される、請求項7に記載の方法。
発明の詳細な説明 【技術分野】
【0001】
本発明は、トピックモデルを用いた遺伝子情報推定装置に関する。
【背景技術】
【0002】
近年、核酸塩基配列を次世代シーケンサで決定することで、遺伝子情報を網羅的に解析する技術が開発されている。例えば、DNA上の遺伝子は、転写、翻訳などを経てその機能を発現するが、この遺伝子発現の状態を解析する目的で、転写産物であるRNAの配列を次世代シーケンサで決定して遺伝子の発現を網羅的に定量するRNA-Seq法が用いられている(非特許文献1)。
【0003】
RNA-Seq等の遺伝子情報には、数千から数万の遺伝子の情報が含まれるため、出力されるデータは高次元のカウントデータとなる。同様の高次元カウントデータを扱うことの多い自然言語処理の分野では、文書の潜在的な意味を扱う手法としてトピックモデルが提案されている。トピックモデルは、文書内の単語の共起に基づき、文書を構成するトピックと各トピックにおける単語の出現頻度を同時に推定する。トピックモデルをRNA-Seqで得られた遺伝子発現データに適用することで、類似した発現パターンを持つ遺伝子群をトピックとして抽出できるのではないかと考えられる。実際、トピックモデルの一つであるLatent Dirichle Allocation(LDA;非特許文献2)を、単一細胞のRNA-Seqデータに適用することで、推定されたトピックから細胞間の階層的な構造を推定できることが示されている(非特許文献3)。
【先行技術文献】
【0004】

【非特許文献1】Wang et al.Nature Reviews Genetics 2009 Jan;10(1):57
【非特許文献2】Blei et al.Journal of Machine Learning Research 2003 3:993
【非特許文献3】DuVerle et al.BMC Bioinformatics 2016 Sep 13;17(1):363
【発明の概要】
【発明が解決しようとする課題】
【0005】
RNA-Seq等の遺伝子情報においては、読んだ配列の数(総リード数;depth)が多いほどノイズの少ない質の高いデータが得られる反面、コストが上がるために扱えるサンプルの数が減るというトレードオフが存在する。総リード数の少ないデータから精度よく発現量を推定できれば、より多くのサンプルの遺伝子情報発現量を定量することができる。RNA-Seqの場合、実際の発現データを見ると、発現パターンの相関した遺伝子が多数存在する。こうした遺伝子間の発現パターンの関係性をトピックモデルにより抽出することで、ノイズの大きな定量データからでも他の遺伝子の情報を利用して高精度な発現量の推定が可能になると期待される。しかし、RNA-Seqデータには、過分散、つまりポアソン分布や二項分布で予想される分散よりも大きな分散を持つという性質がある。従来のトピックモデルでは、多項分布を用いているため、このような大きな分散を説明することができない。
【課題を解決するための手段】
【0006】
本発明では、既存のトピックモデルに基づき、過分散を有する遺伝子情報の測定データにも適用可能な新しいトピックモデル及びその推定方法を提案する。このモデルでは、従来、単語の生成分布に用いられていた多項分布を、負の二項分布に置き換える。
【0007】
すなわち、本発明は、
(i)複数のサンプルに対応する遺伝子情報の測定データを読み込む手順、
(ii)遺伝子情報の測定データに基づき遺伝子の平均出現頻度を算出する手順、
(iii)各サンプルの各遺伝子にトピックが割り当てられる確率を初期化する手順、
(iv)遺伝子へのトピック割り当て確率を更新する手順、
(v)トピック情報を更新する手順、
(vi)トピック情報の分散を更新する手順、
(vii)負の二項分布のdispersionパラメータを更新する手順、
および、
(viii)トピック情報およびトピック出現頻度から遺伝子情報を推定する手順、
を有する、
トピックモデルを用いた遺伝子情報推定装置に関する。
【0008】
手順(iv)、(v)、および(vii)において、過分散の程度を表すパラメータをφとし、サンプル(s)においてトピック(k)を割り当てられた遺伝子(i)の期待出現数をμskiとしたとき、サンプル(s)における遺伝子(i)の出現数(rsi)が、
【数1】
JP2018139043A_000003t.gif
という確率に従うことが好ましい。
【0009】
サンプル(s)における全遺伝子のカウントの総和(ν)、全サンプルでの全遺伝子の出現数の総和に対する遺伝子(i)の出現数の割合の対数(m)、指数分布あるいはジェフリーズ事前分布に従う分散(τki)、および平均が0で分散がτkiの正規分布に従うトピック(k)における遺伝子(i)の期待カウントの平均からのずれを表す量が与えられたときに、サンプル(s)においてトピック(k)を割り当てられた遺伝子(i)の期待出現数(μski)が
μski=νexp(ηki+m
と表される。
【0010】
また、本発明は、前記遺伝子情報推定装置を構成する各手段としてコンピュータを機能させるための、遺伝子情報推定プログラムに関する。
【0011】
また、本発明は、前記遺伝子情報推定プログラムが記録されたことを特徴とするコンピュータ読み取り可能な記録媒体に関する。
【0012】
また、本発明は、
(i)複数のサンプルに対応する遺伝子情報の測定データを読み込むステップ、
(ii)遺伝子情報の測定データに基づき遺伝子の平均出現頻度を算出するステップ、
(iii)各サンプルの各遺伝子にトピックが割り当てられる確率を初期化するステップ、
(iv)遺伝子へのトピック割り当て確率を更新するステップ、
(v)トピック情報を更新するステップ、
(vi)トピック情報の分散を更新するステップ、
(vii)負の二項分布のdispersionパラメータを更新するステップ、
および、
(viii)トピック情報およびトピック出現頻度から遺伝子情報を推定するステップ、
を有する、トピックモデルを用いた遺伝子情報推定方法に関する。
【0013】
ステップ(iv)、(v)、および(vii)において、過分散の程度を表すパラメータをφとし、サンプル(s)においてトピック(k)を割り当てられた遺伝子(i)の期待出現数をμskiとしたとき、サンプル(s)における遺伝子(i)の出現数(rsi)が、
【数2】
JP2018139043A_000004t.gif
という確率に従うことが好ましい。
【0014】
サンプル(s)における全遺伝子のカウントの総和(ν)、全サンプルでの全遺伝子の出現数の総和に対する遺伝子(i)の出現数の割合の対数(m)、指数分布あるいはジェフリーズ事前分布に従う分散(τki)、および平均が0で分散がτkiの正規分布に従うトピック(k)における遺伝子(i)の期待カウントの平均からのずれを表す量が与えられたときに、サンプル(s)においてトピック(k)を割り当てられた遺伝子(i)の期待出現数(μski)が
μski=νexp(ηki+m
と表される。
【発明の効果】
【0015】
過分散のある遺伝子情報を扱うための新しいトピックモデルを提案し、シミュレーションによって推定精度の検証を行った。低コストを想定した総リード数が10程度のデータから本発明のモデルによって推定した遺伝子発現量は、よりコストのかかる総リード数10程度で定量した発現量よりも真の値に近かった(図8)。このことは、1サンプル当たりの総リード数を増やすことで質の良いデータを得るよりも、総リード数を減らして多くのサンプルの定量を行ったうえで本発明のモデルによる発現量の推定を行うことが有効であることを示している。
【図面の簡単な説明】
【0016】
【図1】本発明の実施形態に係る遺伝子情報推定装置1の構成を示すブロック図を示す。
【図2】遺伝子発現推定装置1による処理の流れを示すフローチャートを示す。
【図3A】トピックモデル推定プログラム41の一実施形態による処理の流れを示すフローチャートを示す。
【図3B】トピックモデル推定プログラム41の他の実施形態による処理の流れを示すフローチャートを示す。
【図4】遺伝子情報推定プログラム43による処理の流れを示すフローチャートを示す。
【図5】実施例における対数周辺尤度の下限の近似値を示す。
【図6】実施例において推定された各サンプルにおけるトピックの出現確率を示す。
【図7】実施例における各遺伝子の出現確率のトピック間の差を示す。
【図8】実施例におけるリード数の期待値と推定値との誤差の分布を示す。
【発明を実施するための形態】
【0017】
[1]遺伝子情報推定装置
本発明は、
(i)複数のサンプルに対応する遺伝子情報の測定データを読み込む手順、
(ii)遺伝子情報の測定データに基づき遺伝子の平均出現頻度を算出する手順、
(iii)各サンプルの各遺伝子にトピックが割り当てられる確率を初期化する手順、
(iv)遺伝子へのトピック割り当て確率を更新する手順、
(v)トピック情報を更新する手順、
(vi)トピック情報の分散を更新する手順、
(vii)負の二項分布のdispersionパラメータを更新する手順、
および、
(viii)トピック情報およびトピック出現頻度から遺伝子情報を推定する手順、
を有する、
トピックモデルを用いた遺伝子情報推定装置に関する。

【0018】
手順(i)では複数のサンプルに対応する遺伝子情報の測定データを読み込む。手順(ii)では遺伝子情報の測定データに基づき遺伝子の平均出現頻度を算出する。遺伝子の平均出現頻度の算出法については後述する。手順(iii)では各サンプルの各遺伝子にトピックが割り当てられる確率を初期化する。初期化方法については後述する。

【0019】
遺伝子へのトピック割り当て確率を更新する手順(iv)では、各サンプルにおいて各遺伝子のカウントがどのトピックから生成されたかについて事後分布の推定を行う。

【0020】
トピック情報を更新する手順(v)では、割り当てられたトピックに基づき、トピック情報、すなわち各トピックにおける遺伝子の出現頻度の平均からのずれを推定する。

【0021】
トピック情報の分散を更新する手順(vi)では、推定されたトピック情報に基づき、トピック情報の分散を更新する。

【0022】
負の二項分布のdispersionパラメータを更新する手順(vii)では、トピックごとに各遺伝子のカウントデータを生成する負の二項分布のパラメータを推定する。すなわち、手順(iv)~(v)で割り当てられたトピックと推定されたトピック情報から算出される期待値と測定データの差分に基づき、負の二項分布のdispersionパラメータを推定する。

【0023】
トピック情報およびトピック出現頻度から遺伝子情報を推定する手順(viii)では、トピックモデル推定プログラムにより推定されたトピック情報と平均出現頻度から各トピックにおける遺伝子の出現頻度を算出し、同じくトピックモデル推定プログラムで推定された各サンプルにおけるトピックの出現頻度と掛け合わせて得られる各サンプルにおける遺伝子の出現頻度を遺伝子情報の推定として出力する。

【0024】
[2]遺伝子情報の測定データ
遺伝子情報の測定データは、試料から核酸を調製する工程、核酸からライブラリーを調製し増幅する工程、および次世代シークエンサーによりライブラリーの配列決定を行う工程を経ることにより得られる。各工程における作業手順および条件は特に限定されず、通常の方法を用いることが可能である。試料は限定されず、植物細胞、微生物細胞、動物細胞等を使用可能である。核酸はDNAであってもよく、RNAであってもよい。遺伝子情報は、複数の遺伝子の発現や、複数の遺伝子の存在に関連する情報であれば特に限定されないが、本発明の遺伝子情報推定装置は、過分散を持つ遺伝子情報に好適に適用できる。

【0025】
遺伝子情報の測定データとして、RNA-Seqデータ、ChIP-Seqデータが挙げられる。RNA-Seqデータは、RNAを次世代シークエンサーで定性・定量して得られ、トランスクリプトーム解析に用いられる。ChIP-Seqデータは、クロマチン免疫沈降法と次世代シークエンサーによる配列決定を組み合わせて得られ、DNA結合タンパク質が結合するDNA配列を示す。ChIP-Seqデータの取得時に対象となるDNA結合タンパク質としては、ヒストン、転写因子、DNA修飾酵素が挙げられる。

【0026】
遺伝子情報の測定データとして、DNAのメチル化を示すデータも挙げられる。DNAのメチル化を示すデータは、バイサルファイト変換を行ったゲノムDNAを次世代シークエンサーで解析してメチル化シトシンを検出する方法や、抗メチルシトシン抗体やメチルCpGタンパク質によるDNAの沈降と次世代シークエンサーによる配列決定を組み合わせる方法により得ることができる。

【0027】
本発明の装置は、メタゲノム解析において遺伝子情報を推定するためにも使用可能である。メタゲノム解析は、複数の微生物を含む試料から核酸を抽出し、その核酸塩基配列をまとめて解読し、試料中に存在した微生物の構成を解析する手法である。

【0028】
遺伝子情報の測定データとして、RNA-Seqデータを用いた場合、出力される遺伝子情報は、リード数、すなわち遺伝子発現量の期待値である。これにより、試料におけるトランスクリプトームを高い精度で推定できる。

【0029】
通常、遺伝子情報の推定のためにはサンプルあたりのデータ量が多い必要があるが、本発明の装置によれば、RNA-Seqの場合にはサンプルあたり100,000から1000,000リード程度という少ないデータ量でも高い精度で遺伝子情報を推定できる。

【0030】
[3]遺伝子情報推定装置の全体構成
図1は、本発明の実施形態に係る遺伝子情報推定装置1の構成を示すブロック図である。遺伝子情報推定装置1は、演算手段2と、入力手段3と、記憶手段4と、出力手段5とを備えている。各手段2~5はバスライン6に接続されている。

【0031】
演算手段2は、例えば、CPU(Central Processing Unit)およびRAM(Random Access Memory)から構成される主制御装置である。演算手段2は、トピックモデル推定部21と、遺伝子情報推定部22と、メモリ23を含む。

【0032】
入力手段3は、例えば、キーボード、マウス、ディスクドライブ装置などから構成される。記憶手段4は、例えば、一般的なハードディスク装置などから構成され、プログラム格納部40と、データ格納部44とを含む。

【0033】
プログラム格納部40には、演算手段2で用いられるプログラムとして、トピックモデル推定プログラム41、および遺伝子情報推定プログラム43とを記憶させておくことが可能である。演算手段2は、記憶手段4のプログラム格納部40から、トピックモデル推定プログラム41、および遺伝子情報推定プログラム43をそれぞれ読み込み、メモリ23に格納し実行することで、トピックモデル推定部21と、遺伝子情報推定部22とをそれぞれ実現する。

【0034】
データ格納部44には、演算手段2で用いられる各種データとして、遺伝子情報の測定データ45、遺伝子の平均出現頻度46、トピック情報47、トピック出現頻度48とを記憶する。遺伝子情報の測定データ45は入力手段3を介して入力され、記憶手段4のデータ格納部44に記憶される構成とすることが可能である。遺伝子の平均出現頻度46、トピック情報47、トピック出現頻度48は、演算手段2の演算処理結果を示すデータである。

【0035】
出力手段5は、例えば、グラフィックボード(出力インタフェース)およびそれに接続されたモニタである。モニタは、例えば、液晶ディスプレイ等から構成され、可視化を行った結果等を表示する。

【0036】
本発明の課題は、前述の遺伝子情報推定装置を構成する各手段としてコンピュータを機能させるための遺伝子情報推定プログラム、前記遺伝子情報推定プログラムが記録されたコンピュータ読み取り可能な記録媒体によっても解決できる。また、本発明の課題は、前述の遺伝子情報推定装置の手順と同様のステップ含む、トピックモデルを用いた遺伝子情報推定方法によっても解決できる。

【0037】
遺伝子情報推定装置は、一般的なコンピュータを、遺伝子情報推定プログラムにより動作させることで実現できる。このプログラムは、通信回線を介して提供することも可能であるし、CD-ROM等のコンピュータ読み取り可能な記録媒体に書き込んで配布することも可能である。コンピュータは、コンピュータ読み取り可能な記録媒体から遺伝子情報推定プログラムを読み出し、このプログラムに基づいた各機能を実現することができる。また、遺伝子情報推定プログラムをインストールされたコンピュータは、CPUが、ROM等に格納されたこのプログラムをRAMに展開することにより、遺伝子情報推定装置としての効果を奏することができる。

【0038】
[4]遺伝子情報推定装置による処理の流れ
図2は、図1に示した遺伝子情報推定装置1による処理の流れを示すフローチャートである。S11において、遺伝子情報推定装置1は遺伝子情報の測定データ45を読み込む。S12において、遺伝子情報推定装置1は遺伝子情報の測定データ45から遺伝子の平均出現頻度46を算出する。

【0039】
遺伝子の平均出現頻度46は、例えば遺伝子の全サンプルを通しての出現回数を、その全ての遺伝子についての総和で割ることで算出できる。

【0040】
S13において、遺伝子情報推定装置1はトピックモデル推定プログラム41によりトピック情報47およびトピック出現頻度48を推定する。S14において、遺伝子情報推定装置1はトピック情報47およびトピック出現頻度48から遺伝子情報を推定する。S15において、遺伝子情報推定装置1は推定された遺伝子情報を出力手段5から出力する。

【0041】
[5]トピックモデル推定プログラムによる処理の流れ
図3Aは、図1に示したトピックモデル推定プログラム41による処理の流れを示すフローチャートである。

【0042】
S201において、遺伝子情報推定装置1は、トピックモデル推定プログラム41の初期化を行う。S202において、遺伝子情報推定装置1は遺伝子へのトピック割り当て確率の更新を行う。S203において、遺伝子情報推定装置1はトピック情報(η)の更新1を行う。S204において、遺伝子情報推定装置1はトピック情報の分散(τ)の更新を行う。S205において、遺伝子情報推定装置1は負の二項分布のdispersionパラメータ(φ)の更新を行う。S206において、遺伝子情報推定装置1は対数周辺尤度の下限が収束しているか否かを評価し、収束していない場合にはS202に戻る。対数周辺尤度の下限が収束している場合にはS211に進む。S211において、遺伝子情報推定装置1はトピック情報47およびトピック出現頻度48を出力手段5から出力する。

【0043】
以下、トピックモデル推定プログラムにおける処理の流れの詳細を、必要に応じて図3Aを参照しながら説明する。過分散なデータを単一のトピックから発生させるために、先行技術であるSparse Additive Generative model(SAGE;Eisenstein et al.Proc.Int.Conf.Mach.Learn.2011)における単語あるいは遺伝子の出現確率を負の二項分布に置き換える。サンプルで各遺伝子に割り当てられるトピックの確率は、SAGEと同様にDirichlet分布から生成される。

【0044】
【数3】
JP2018139043A_000005t.gif
【数4】
JP2018139043A_000006t.gif

【0045】
ここで、θはサンプルsにおいて各トピックを割り当てる確率、zsiはサンプルsのi番目の遺伝子に割り当てられたトピック、sとGはそれぞれサンプルと遺伝子の数を表す。先行技術と同様に、すべてのトピックにおける平均的な遺伝子の出現頻度の対数をm、各トピックにおける出現頻度がこのmからどの程度ずれているかを表す量をηとし、その事前分布は以下の正規分布とする。

【0046】
【数5】
JP2018139043A_000007t.gif

【0047】
分散τkiは指数分布
【数6】
JP2018139043A_000008t.gif
あるいはJeffreys事前分布
【数7】
JP2018139043A_000009t.gif
を用いる。こうすることで、ηは多くの要素が0となることが知られている。
サンプルsにおける遺伝子iにトピックkが割り当てられたとき、すなわちzsi=kのとき、そのリード数は、
【数8】
JP2018139043A_000010t.gif
という負の二項分布に従うものとする。ここで、φは負の二項分布において過分散の度合いを決めるdispersionというパラメータである。この分布の期待出現数μskiは、サンプルsにおける全遺伝子のカウントの総和をνと表記すると、
【数9】
JP2018139043A_000011t.gif
となる。

【0048】
トピックモデル推定プログラムの初期化処理(S201)では、各サンプルの各遺伝子にトピックの割り当てられる確率をランダムに初期化する。また、入力されたカウントデータから平均的な出現頻度の対数mを算出し、上記のトピックの割り当てられる確率と合わせて平均出現頻度からのずれηkiの初期値を決める。

【0049】
遺伝子へのトピック割り当て確率の更新(S202)は下記の手順により行う。すなわち、RNA-Seqデータに基づき、z、η、τの事後分布の推定を行う。これらの事後分布を解析的に求めることはできないため、事後分布の近似q(z)、q(η)、q(τ)を考え、対数周辺尤度の下限
【数10】
JP2018139043A_000012t.gif
を最大化する。ここで
【数11】
JP2018139043A_000013t.gif
は確率分布q()のもとでの期待値を意味する。同様に確率分布q()のもとでの分散を
【数12】
JP2018139043A_000014t.gif
、サンプルsの遺伝子iを除いたzをz\s,i、サンプルsにおいて遺伝子i以外でトピックkを割り当てられた遺伝子の数をns,k\s,iと表記し、
【数13】
JP2018139043A_000015t.gif
とおく。

【0050】
トピックモデルの推定は、τkiの事前分布を指数分布として推定を行う。最初に遺伝子へのトピックの割り当て確率の更新(S202)を行う。トピックの割り当て確率q(z)の更新式は、
【数14】
JP2018139043A_000016t.gif
となる。

【0051】
次に、トピック情報(η)の更新1(S203)を行う。q(η)はラプラス近似を行う。関数
【数15】
JP2018139043A_000017t.gif
を定義し、
【数16】
JP2018139043A_000018t.gif
と表記すると
【数17】
JP2018139043A_000019t.gif
となる。ここで、
【数18】
JP2018139043A_000020t.gif
の推定値はニュートン法による更新を1ステップにつき1回行う、
【数19】
JP2018139043A_000021t.gif

【0052】
次に、トピック情報の分散(τ)の更新(S204)を行う。τの事後分布q(τ)はガンマ分布
【数20】
JP2018139043A_000022t.gif
とし、そのパラメータak,iとbk,i
【数21】
JP2018139043A_000023t.gif
【数22】
JP2018139043A_000024t.gif
で更新する。

【0053】
以上の推定値から算出される期待出現頻度と測定されたカウントデータrsiとの差分に基づき、負の二項分布のdispersionパラメータφの更新を行う(S205)。以上の更新(S202~S205)を対数周辺尤度の下限が収束するまで繰り返す(S206)。対数周辺尤度の下限は各更新ステップで導出された近似に基づき算出される。

【0054】
より高精度の推定を行うために、トピックモデルの推定は複数段階行ってもよい。トピックモデルの推定を二段階行う場合のトピック推定プログラム41による処理の流れを図3Bに示す。

【0055】
図3Bにおいて、S201~S206およびS211は図3Aと同じである。S206において、遺伝子情報推定装置1は対数周辺尤度の下限が収束しているか否かを評価し、収束していない場合にはS202に戻る。対数周辺尤度の下限が収束している場合にはS207に進む。S207において、遺伝子情報推定装置1は遺伝子へのトピック割り当て確率の更新を行う。S208において、遺伝子情報推定装置1はトピック情報(η)の更新2を行う。S209において、遺伝子情報推定装置1は、負の二項分布のdispersionパラメータ(φ)の更新を行う。S210において、遺伝子情報推定装置1は対数周辺尤度の下限が収束しているか否かを評価し、収束していない場合にはS207に戻る。対数周辺尤度の下限が収束している場合にはS211に進む。S211において、遺伝子情報推定装置1はトピック情報47およびトピック出現頻度48を出力手段5から出力する。

【0056】
トピックモデルの推定を二段階行う場合、まず、τkiの事前分布を指数分布として前述の方法で推定を行う。対数周辺尤度の下限が収束した後、よりスパースなηを得るために、τkiの事前分布をJeffreys事前分布として遺伝子へのトピック割り当て確率(S207)、トピック情報(S208)、負の二項分布のdispersionパラメータ(S209)の更新を、対数周辺尤度の下限がもう一度収束するまで繰り返す。q(z)は一段階目と同じ式で更新できる。また、q(τ)はキャンセルされて消えるため、二段階目では更新を行わない。q(η)の更新のみ、
【数23】
JP2018139043A_000025t.gif
中の
【数24】
JP2018139043A_000026t.gif
が、
【数25】
JP2018139043A_000027t.gif
となるために更新式が異なる(S208)。

【0057】
遺伝子へのトピック割り当て確率の更新、トピック情報(η)の更新2、および負の二項分布のdispersionパラメータの更新(S207~S209)を、対数周辺尤度の下限が収束するまで繰り返す(S210)。対数周辺尤度の下限の収束は、S206と同じ方法により評価できる。対数周辺尤度の下限が収束した後、トピック情報47およびトピック出現頻度48が出力される(S211)。

【0058】
S202~S205、およびS207~S209において、dispersionφは負の二項分布の分散と期待値とdispersionの関係、
【数26】
JP2018139043A_000028t.gif
から、リード数の期待値と実測の二乗誤差に基づいて推定する。

【0059】
[5]遺伝子情報推定プログラムによる処理の流れ
図4は、図1に示した遺伝子情報推定プログラム43による処理の流れを示すフローチャートである。S31において、遺伝子情報推定装置1は、遺伝子の平均出現頻度46、トピック情報47、およびトピック出現頻度48を読み込む。遺伝子の平均出現頻度46は図2のS12で算出されたものであり、トピック情報47、およびトピック出現頻度48は図3Aおよび図3BのS211で出力されたものである。S32において、遺伝子情報推定装置1は遺伝子情報を推定する。S33において、遺伝子情報推定装置1は、推定された遺伝子情報を出力手段5から出力する。

【0060】
S32における遺伝子情報の推定において、サンプル(s)における遺伝子(i)の期待出現数
【数27】
JP2018139043A_000029t.gif
は、
サンプル(s)における全遺伝子のカウントの総和(ν)、トピックモデル推定プログラム(S13)により推定された全サンプルでの全遺伝子の出現数の総和に対する遺伝子(i)の出現数の割合の対数(m)、および、以上で推定された各サンプルにおける各遺伝子へのトピックの割り当て確率、トピック(k)における遺伝子(i)の期待カウントの平均からのずれを表す量(η)から、
【数28】
JP2018139043A_000030t.gif
と推定される。
【実施例】
【0061】
(実施例1)トピックの出現確率の推定
先行技術との比較を行うために、疑似的なRNA-Seqデータを作成して解析を行った。統制群と実験群がそれぞれ20サンプルある状況を想定した。遺伝子の総数を10,000とし、そのうち500遺伝子は実験群で発現量が増加し、別の500遺伝子は発現量が減少する。残りの9,000遺伝子は二群間で発現量が変わらないとする。統制群における各遺伝子のリード数は以下の負の二項分布に従う。
【実施例】
【0062】
【数29】
JP2018139043A_000031t.gif
【数30】
JP2018139043A_000032t.gif
【数31】
JP2018139043A_000033t.gif
【実施例】
【0063】
【数32】
JP2018139043A_000034t.gif
は総リード数の期待値で、低コストで総リード数の少ない条件
【数33】
JP2018139043A_000035t.gif
と高コストで総リード数の多い条件
【数34】
JP2018139043A_000036t.gif
の二条件を想定した。dispersionφは、公開されている実測のデータ(Pickerell et al.Nature 2010)に基づきedgeR(Zhou et al.Nucleic Acids Res.2014)で推定を行った。実験群の期待リード数μは以下のように決定した。
【実施例】
【0064】
【数35】
JP2018139043A_000037t.gif
【数36】
JP2018139043A_000038t.gif
【実施例】
【0065】
ここで、ρは実験群で発現量が増加する遺伝子群については、ρ=+1、減少する遺伝子群についてはρ=-1、残りの変化しない遺伝子群については、ρ=0とした。生成した疑似データについて、LDA、SAGE、本発明の提案モデルの3つのモデルを適用した。
【実施例】
【0066】
異なるトピック数で提案モデルの推定を行った時の、対数周辺尤度の下限の近似値を図5に示す。図5の横軸はトピック数を表し、縦軸は変分下限の近似値を表す。トピック数が2の時に最大となることから、以降の解析ではトピック数を2とした。疑似データは統制群と実験群の二群からなるため、この結果はトピック数を適切に選択できていることを示している。
【実施例】
【0067】
3つのモデルで推定された各サンプルにおけるトピックの出現確率を図6に示す。図6において、横軸はサンプルを表し、縦軸はトピックの出現確率を表す。統制群と実験群を左右に分けて表示している。すべてのモデルで、二群間でトピックの出現確率が変化しているが、LDAとSAGEに対して、本発明のモデル(proposed)では二群間の違いが顕著になっている。
【実施例】
【0068】
(実施例2)トピックにおける遺伝子の出現確率の比較
各遺伝子について推定された二つのトピックにおける出現確率の差を比較した。LDAについては、各トピックにおける出現確率の差を、SAGEと本発明のモデルについては、
【数37】
JP2018139043A_000039t.gif
を計算し、ヒストグラムを作成した(図7)。図7において、実験群で発現量が増加する500遺伝子(up)を白色、減少する500遺伝子(down)を黒色、それ以外の9,000遺伝子(const)を灰色で表示している。LDAでは、ほとんどの遺伝子が0付近に集中しており、発現量の増加する遺伝子と減少する遺伝子とそれ以外の遺伝子を区別することが困難である。3種類の遺伝子は、SAGEよりも本発明のモデル(proposed)の方がより明瞭に分離している。
【実施例】
【0069】
(実施例3) 発現量の推定精度の比較
発現量の推定精度を比較するために、リード数の期待値μsiと推定値の誤差を算出した(図8)。図8ではリード数の期待値の大きさに応じて分割してプロットした。期待総リード数10で生成した疑似データ(raw)、コストをかけた状態を想定して期待総リード数を10とした疑似データ(deep)、SAGEの推定(sage)、本発明のモデルの推定(proposed)を、それぞれ黒色の実線、灰色の実線、黒色の破線、灰色の破線でプロットした。
【実施例】
【0070】
総リード数の多い状況を想定したデータ(deep)は総リード数をそろえるために、生成したデータを100で割った値と期待値を比較した。LDAは、元のデータ(raw)よりも誤差が大きくなっており、発現量の推定には使用できない。sageと本発明のモデルを比較すると、約60%の観測で本発明のモデルの方が期待値に近い推定を与えている。
【実施例】
【0071】
本発明のモデルは、今回のシミュレーションの二群比較のように統制のとれた実験のサンプルだけではなく、edgeRやDESeqといった負の二項分布を用いた既存の解析法が利用できない未知の構成のサンプルから定量したRNA-Seqデータに対しても、有効と考えられる。
【符号の説明】
【0072】
1 遺伝子情報推定装置
2 演算手段
21 トピックモデル推定部
22 遺伝子情報推定部
23 メモリ
3 入力手段
4 記憶手段
40 プログラム格納部
41 トピックモデル推定プログラム
43 遺伝子情報推定プログラム
44 データ格納部
45 遺伝子情報の測定データ
46 遺伝子の平均出現頻度
47 トピック情報
48 トピック出現頻度
5 出力手段
6 バスライン

図面
【図1】
0
【図2】
1
【図3A】
2
【図3B】
3
【図4】
4
【図5】
5
【図6】
6
【図7】
7
【図8】
8