ホーム > 特許ランキング > 国立大学法人 東京大学
知財求人 - 知財ポータルサイト「IP Force」
▶ 国立大学法人 東京大学の特許一覧 ▶ トヨタ自動車株式会社の特許一覧 ▶ 独立行政法人農業・食品産業技術総合研究機構の特許一覧 ▶ 学校法人 龍谷大学の特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2022-02-15
(45)【発行日】2022-02-24
(54)【発明の名称】トランスクリプトーム解析装置及び解析方法
(51)【国際特許分類】
G16B 40/00 20190101AFI20220216BHJP
【ＦＩ】
G16B40/00
【請求項の数】
12
(21)【出願番号】P 2018003697
(22)【出願日】2018-01-12
(65)【公開番号】P2019125045
(43)【公開日】2019-07-25
【審査請求日】2020-08-07
(73)【特許権者】
【識別番号】504137912
【氏名又は名称】国立大学法人 東京大学
(73)【特許権者】
【識別番号】000003207
【氏名又は名称】トヨタ自動車株式会社
(73)【特許権者】
【識別番号】501203344
【氏名又は名称】国立研究開発法人農業・食品産業技術総合研究機構
(73)【特許権者】
【識別番号】597065329
【氏名又は名称】学校法人 龍谷大学
(74)【代理人】
【識別番号】110002572
【氏名又は名称】特許業務法人平木国際特許事務所
(72)【発明者】
【氏名】近藤 聡
(72)【発明者】
【氏名】大音 徳
(72)【発明者】
【氏名】阿部 円佳
(72)【発明者】
【氏名】青木 直大
(72)【発明者】
【氏名】福田 あかり
(72)【発明者】
【氏名】廣瀬 竜郎
(72)【発明者】
【氏名】永野 惇
【審査官】山内 裕史
(56)【参考文献】
【文献】特開２００４－３５５１７４（ＪＰ，Ａ）
【文献】再公表特許第２００３／０８５５４８（ＪＰ，Ａ１）
【文献】特開２０１６－０４５８０１（ＪＰ，Ａ）
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ１６Ｂ ５／００ － ９９／００
(57)【特許請求の範囲】
【請求項1】
目的変数データと遺伝子発現量データとを含む複数のデータセットに対して、遺伝子発現量データをランダムに削減した第1～第mのサブデータセット（m≧2）を生成するデータセット生成手段と、第1～第mのサブデータセットのそれぞれに対して正則化項を有する回帰分析法を適用して、目的変数データを目的変数とし遺伝子発現量データを説明変数とする第1～第mの予測式を算出する予測式算出手段と、
第1～第mの予測式に含まれる遺伝子発現量データに対応する遺伝子のリストを生成する遺伝子リスト生成手段と
を備えるトランスクリプトーム解析装置。
【請求項2】
上記予測式算出手段は、上記回帰分析法としてLASSO（least absolute shrinkage and selection operator）を適用することを特徴とする請求項１記載のトランスクリプトーム解析装置。
【請求項3】
上記データセット生成手段は、1000～20000通りのサブデータセット（m=1000～20000）を生成することを特徴とする請求項１記載のトランスクリプトーム解析装置。
【請求項4】
上記遺伝子リスト生成手段は、第1～第mの予測式に基づいて遺伝子の出現確率を算出し、算出した出現確率と関連付けて遺伝子のリストを生成することを特徴とする請求項１記載のトランスクリプトーム解析装置。
【請求項5】
上記遺伝子リスト生成手段は、遺伝子のアノテーション情報が格納されたデータベースから、リストに含まれる遺伝子のアノテーション情報を読み出し、読み出したアノテーション情報と関連づけて遺伝子のリストを生成することを特徴とする請求項１記載のトランスクリプトーム解析装置。
【請求項6】
上記遺伝子リスト生成手段により生成したリストに含まれる複数の遺伝子について、上記複数のデータセットに含まれる目的変数データと遺伝子発現量データとを用いた重回帰分析により、所定の目的変数に関する予測モデル式を生成する予測モデル式生成手段を更に有することを特徴とする請求項１記載のトランスクリプトーム解析装置。
【請求項7】
中央処理装置が、目的変数データと遺伝子発現量データとを含む複数のデータセットに対して、遺伝子発現量データをランダムに削減したサブデータセットを生成するサブデータセット生成工程と、中央処理装置が、サブデータセットに対して正則化法を適用して、目的変数データを目的変数とし遺伝子発現量データを説明変数とする予測式を算出する予測式算出工程と、
記憶装置が、予測式に含まれる遺伝子発現量データに対応する遺伝子を記録する遺伝子記録工程と、
上記サブデータセット生成工程、上記予測式算出工程及び上記遺伝子記録工程をm回（m≧2）繰り返し、中央処理装置が、記録された遺伝子のリストを生成する遺伝子リスト生成工程と
を含むトランスクリプトーム解析方法。
【請求項8】
上記予測式算出工程では、中央処理装置が、上記正則化法としてLASSO（least absolute shrinkage and selection operator）を適用することを特徴とする請求項７記載のトランスクリプトーム解析方法。
【請求項9】
上記サブデータセット生成工程では、中央処理装置が、1000～20000通りのサブデータセット（n=1000～20000）を生成することを特徴とする請求項７記載のトランスクリプトーム解析方法。
【請求項10】
上記遺伝子リスト生成工程では、中央処理装置が、第1～第m回の繰り返しで生成した第1～第mの予測式に基づいて遺伝子の出現確率を算出し、算出した出現確率と関連付けて遺伝子のリストを生成することを特徴とする請求項７記載のトランスクリプトーム解析方法。
【請求項11】
上記遺伝子リスト生成工程では、遺伝子のアノテーション情報が格納されたデータベースから、中央処理装置が、リストに含まれる遺伝子のアノテーション情報を読み出し、読み出したアノテーション情報と関連づけて遺伝子のリストを生成することを特徴とする請求項７記載のトランスクリプトーム解析方法。
【請求項12】
上記遺伝子リスト生成工程の後、生成したリストに含まれる複数の遺伝子について、上記複数のデータセットに含まれる目的変数データと遺伝子発現量データとを用いた重回帰分析により、中央処理装置が、所定の目的変数に関する予測モデル式を生成する予測モデル式生成工程を更に有することを特徴とする請求項７記載のトランスクリプトーム解析方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、トランスクリプトームに関する情報を解析するトランスクリプトーム解析装置及び解析方法に関する。
【背景技術】
【0002】
遺伝子発現に基づいて生物の表現型を予測する試みとして、遺伝子発現データと表現型データとから重回帰分析する方法が知られている（非特許文献１及び特許文献１）。非特許文献１に開示された方法では、遺伝子発現データの重複をなくすため、同じオペロンについては最も発現レベルが高いデータのみを適用するなどして、遺伝子発現データを限定していた。
【0003】
ところで、トランスクリプトームは、一般的に、所定の状態や条件下における組織或いは細胞内に存在する全ての転写産物を意味する。トランスクリプトームは、ゲノム上のコーディング領域からの転写産物（すなわち、mRNA）と、非コーディング領域からの転写産物（いわゆるncRNA）とを含む。トランスクリプトームを解析することによって、環境要因による遺伝子発現の変動、表現型に関連して発現する遺伝子の同定など、遺伝子の発現状態に基づいた新たな知見を得ることができる。
【0004】
トランスクリプトームを解析する際には、例えば、組織や細胞内に存在する転写産物をマイクロアレイ技術や次世代シーケンス技術を適用して網羅的に計測する。計測されたデータは、大量の塩基配列データであり典型的なビックデータである。
【0005】
得られたデータを統計的に解析する方法としては、特許文献２に開示されるように、多変量解析の一手法である主成分分析を適用する方法が知られている。当該方法では、分析により得られた塩基配列データではなくトレーニングデータについて主成分分析することで、条件の異なるサンプル間で比較可能な結果を導くことができる。
【0006】
また、トランスクリプトーム解析法としては、特許文献３に開示されるように、遺伝子発現情報（状態変数）と形質情報（特性変数）から、解析対象の特性変数推定モデルを生成する方法が知られている。特許文献３に開示された方法では、特性変数を目的変数（従属変数）、状態変数のそれぞれを説明変数として、正則化項を有する回帰モデルを生成している。回帰モデルの算出式として、LASSO回帰（Least Absolute Shrinkage and Selection Operator）が例示されている。
【0007】
ところでLASSO回帰とは、統計学や機械学習の分野における過剰適合を防ぐために用いられる正則化の一手法（L1型正則化法）であり、大量のデータのうち重要でないデータのパラメータを０としてデータから削除する、スパース正則化法に基づく回帰モデリングである（非特許文献２）。
【先行技術文献】
【特許文献】
【0008】
【文献】ＷＯ２０１６／１４８１０７
【文献】特許第５８５４３４６号
【文献】特開２０１７－５１１１８号公報
【非特許文献】
【0009】
【文献】Nature Communications 5, Article number: 5792 (2014)
【文献】Robert Tibshirani, Journal of the Royal Statistical Society. Series B (Methodological) Vol. 58, No. 1 (1996), pp. 267-288
【発明の概要】
【発明が解決しようとする課題】
【0010】
ところで、上述したトランスクリプトーム解析においては、解析対象のサンプル数と比較して、塩基配列データが得られる転写産物の数が極めて大きいため、非特許文献１に開示された方法では十分に意味のある解析結果を得ることが困難であった。また、特許文献３に開示されたLASSO回帰分析を適用した解析方法については、解析対象のサンプル数と比較して、塩基配列データが得られる転写産物の数が極めて大きい場合であっても良好な解析結果が期待される。しかしながら、トランスクリプトーム解析においては、解析結果の更なる精度向上が求められていた。
【0011】
そこで、本発明は、上述した実情に鑑み、転写産物に関する塩基配列データを用いて、より高精度なトランスクリプトーム解析を行うことができるトランスクリプトーム解析装置及び解析方法を提供することを目的とする。
【課題を解決するための手段】
【0012】
上述した目的を達成した本発明は以下を包含する。
（１）目的変数データと遺伝子発現量データとを含む複数のデータセットに対して、遺伝子発現量データをランダムに削減した第1～第mのサブデータセット（m≧2）を生成するデータセット生成手段と、第1～第mのサブデータセットのそれぞれに対して正則化項を有する回帰分析法を適用して、目的変数データを目的変数とし遺伝子発現量データを説明変数とする第1～第mの予測式を算出する予測式算出手段と、第1～第mの予測式に含まれる遺伝子発現量データに対応する遺伝子のリストを生成する遺伝子リスト生成手段とを備えるトランスクリプトーム解析装置。
（２）上記予測式算出手段は、上記回帰分析法としてLASSO（least absolute shrinkage and selection operator）を適用することを特徴とする（１）記載のトランスクリプトーム解析装置。
（３）上記データセット生成手段は、1000～20000通りのサブデータセット（m=1000～20000）を生成することを特徴とする（１）記載のトランスクリプトーム解析装置。
（４）上記遺伝子リスト生成手段は、第1～第mの予測式に基づいて遺伝子の出現確率を算出し、算出した出現確率と関連付けて遺伝子のリストを生成することを特徴とする（１）記載のトランスクリプトーム解析装置。
（５）上記遺伝子リスト生成手段は、遺伝子のアノテーション情報が格納されたデータベースから、リストに含まれる遺伝子のアノテーション情報を読み出し、読み出したアノテーション情報と関連づけて遺伝子のリストを生成することを特徴とする（１）記載のトランスクリプトーム解析装置。
（６）上記遺伝子リスト生成手段により生成したリストに含まれる複数の遺伝子について、上記複数のデータセットに含まれる目的変数データと遺伝子発現量データとを用いた重回帰分析により、所定の目的変数に関する予測モデル式を生成する予測モデル式生成手段を更に有することを特徴とする（１）記載のトランスクリプトーム解析装置。
（１）目的変数データと遺伝子発現量データとを含む複数のデータセットに対して、遺伝子発現量データをランダムに削減した第1～第mのサブデータセット（m≧2）を生成するデータセット生成手段と、第1～第mのサブデータセットのそれぞれに対して正則化項を有する回帰分析法を適用して、目的変数データを目的変数とし遺伝子発現量データを説明変数とする第1～第mの予測式を算出する予測式算出手段と、第1～第mの予測式に含まれる遺伝子発現量データに対応する遺伝子のリストを生成する遺伝子リスト生成手段とを備えるトランスクリプトーム解析装置。
（２）上記予測式算出手段は、上記回帰分析法としてLASSO（least absolute shrinkage and selection operator）を適用することを特徴とする（１）記載のトランスクリプトーム解析装置。
（３）上記データセット生成手段は、1000～20000通りのサブデータセット（m=1000～20000）を生成することを特徴とする（１）記載のトランスクリプトーム解析装置。
（４）上記遺伝子リスト生成手段は、第1～第mの予測式に基づいて遺伝子の出現確率を算出し、算出した出現確率と関連付けて遺伝子のリストを生成することを特徴とする（１）記載のトランスクリプトーム解析装置。
（５）上記遺伝子リスト生成手段は、遺伝子のアノテーション情報が格納されたデータベースから、リストに含まれる遺伝子のアノテーション情報を読み出し、読み出したアノテーション情報と関連づけて遺伝子のリストを生成することを特徴とする（１）記載のトランスクリプトーム解析装置。
（６）上記遺伝子リスト生成手段により生成したリストに含まれる複数の遺伝子について、上記複数のデータセットに含まれる目的変数データと遺伝子発現量データとを用いた重回帰分析により、所定の目的変数に関する予測モデル式を生成する予測モデル式生成手段を更に有することを特徴とする（１）記載のトランスクリプトーム解析装置。
【0013】
（７）目的変数データと遺伝子発現量データとを含む複数のデータセットに対して、遺伝子発現量データをランダムに削減したサブデータセットを生成するサブデータセット生成工程と、サブデータセットに対して正則化法を適用して、目的変数データを目的変数とし遺伝子発現量データを説明変数とする予測式を算出する予測式算出工程と、予測式に含まれる遺伝子発現量データに対応する遺伝子を記録する遺伝子記録工程と、上記サブデータセット生成工程、上記予測式算出工程及び上記遺伝子記録工程をm回（m≧2）繰り返し、記録された遺伝子のリストを生成する遺伝子リスト生成工程とを含むトランスクリプトーム解析方法。
（８）上記予測式算出工程では、上記正則化法としてLASSO（least absolute shrinkage and selection operator）を適用することを特徴とする（７）記載のトランスクリプトーム解析方法。
（９）上記サブデータセット生成工程では、1000～20000通りのサブデータセット（n=1000～20000）を生成することを特徴とする（７）記載のトランスクリプトーム解析方法。
（１０）上記遺伝子リスト生成工程では、第1～第m回の繰り返しで生成した第1～第mの予測式に基づいて遺伝子の出現確率を算出し、算出した出現確率と関連付けて遺伝子のリストを生成することを特徴とする（７）記載のトランスクリプトーム解析方法。
（１１）上記遺伝子リスト生成工程では、遺伝子のアノテーション情報が格納されたデータベースから、リストに含まれる遺伝子のアノテーション情報を読み出し、読み出したアノテーション情報と関連づけて遺伝子のリストを生成することを特徴とする（７）記載のトランスクリプトーム解析方法。
（１２）上記遺伝子リスト生成工程の後、生成したリストに含まれる複数の遺伝子について、上記複数のデータセットに含まれる目的変数データと遺伝子発現量データとを用いた重回帰分析により、所定の目的変数に関する予測モデル式を生成する予測モデル式生成工程を更に有することを特徴とする（７）記載のトランスクリプトーム解析方法。
（８）上記予測式算出工程では、上記正則化法としてLASSO（least absolute shrinkage and selection operator）を適用することを特徴とする（７）記載のトランスクリプトーム解析方法。
（９）上記サブデータセット生成工程では、1000～20000通りのサブデータセット（n=1000～20000）を生成することを特徴とする（７）記載のトランスクリプトーム解析方法。
（１０）上記遺伝子リスト生成工程では、第1～第m回の繰り返しで生成した第1～第mの予測式に基づいて遺伝子の出現確率を算出し、算出した出現確率と関連付けて遺伝子のリストを生成することを特徴とする（７）記載のトランスクリプトーム解析方法。
（１１）上記遺伝子リスト生成工程では、遺伝子のアノテーション情報が格納されたデータベースから、リストに含まれる遺伝子のアノテーション情報を読み出し、読み出したアノテーション情報と関連づけて遺伝子のリストを生成することを特徴とする（７）記載のトランスクリプトーム解析方法。
（１２）上記遺伝子リスト生成工程の後、生成したリストに含まれる複数の遺伝子について、上記複数のデータセットに含まれる目的変数データと遺伝子発現量データとを用いた重回帰分析により、所定の目的変数に関する予測モデル式を生成する予測モデル式生成工程を更に有することを特徴とする（７）記載のトランスクリプトーム解析方法。
【発明の効果】
【0014】
本発明に係るトランスクリプトーム解析装置及び解析方法によれば、トランスクリプトームに関する高精度な解析が可能となる。したがって、本発明に係るトランスクリプトーム解析装置及び解析方法を適用することによって、例えば、所定の状態や条件といった要因による遺伝子発現の変動解析、表現型に関連する遺伝子の発現解析、或いは、遺伝子発現に基づく形質の予測解析等を高精度に行うことができる。
【図面の簡単な説明】
【0015】
【図1】本発明に係るトランスクリプトーム解析装置の一実施形態を示す機能ブロック図である。
【図2】本発明に係るトランスクリプトーム解析方法の一実施形態を示すフローチャートである。
【図3】トランスクリプトーム解析装置及び解析方法で出力される遺伝子のリストの一例を示す特性図である。
【図4】トランスクリプトーム解析装置及び解析方法で出力される遺伝子のリストの他の例を示す特性図である。
【図5】本発明に係るトランスクリプトーム解析装置の他の実施形態を示す機能ブロック図である。
【図6】本発明に係るトランスクリプトーム解析方法の他の実施形態を示すフローチャートである。
【図7】トランスクリプトーム解析装置及び解析方法で出力される予測値と実測値との関係を示す特性図である
【図8】本発明を適用した予測評価システムの構成を示すブロック図である。
【図9】Arroz da TerraとOuu365の発芽14日後を撮像した写真である。
【図10】BIL104系統を用いて地上部乾物重のQTL解析を行った結果を示す特性図である。
【図11】実施例で作出した系統の数と地上部乾燥重量との関係を示す特性図である。
【図12】実施例で作出した系統について、地上部生鮮重を示す特性図である。
【図13-1】実施例において、発現量バイオマーカーとして高頻度で選出された158遺伝子のリストを示す特性図である。
【図13-2】実施例において、発現量バイオマーカーとして高頻度で選出された158遺伝子のリストを示す特性図である。
【図13-3】実施例において、発現量バイオマーカーとして高頻度で選出された158遺伝子のリストを示す特性図である。
【図13-4】実施例において、発現量バイオマーカーとして高頻度で選出された158遺伝子のリストを示す特性図である。
【図13-5】実施例において、発現量バイオマーカーとして高頻度で選出された158遺伝子のリストを示す特性図である。
【図13-6】実施例において、発現量バイオマーカーとして高頻度で選出された158遺伝子のリストを示す特性図である。
【図13-7】実施例において、発現量バイオマーカーとして高頻度で選出された158遺伝子のリストを示す特性図である。
【図13-8】実施例において、発現量バイオマーカーとして高頻度で選出された158遺伝子のリストを示す特性図である。
【図13-9】実施例において、発現量バイオマーカーとして高頻度で選出された158遺伝子のリストを示す特性図である。
【図14】実施例で作出したRNA-seq解析に用いたBIL系統と親品種について、qLTG3-1発現量と地上部生鮮重との関係を示す特性図である。
【図15】実施例で作出した、RNA-seq解析に用いたBIL系統と親品種について、SG-1発現量と地上部生鮮重との関係を示す特性図である。
【図16】実施例で作出した、BIL104系統すべてと親品種について、SG-1発現量と地上部生鮮重との関係を示す特性図である。
【図17】実施例１で作成したリストに含まれる遺伝子に関する発現量データ及び地上部生鮮重データとから算出した地上部生鮮重の予測値と、地上部生鮮重の実測値との関係を示す特性図である。
【発明を実施するための形態】
【0016】
以下、本発明に係るトランスクリプトーム解析装置及び/又は解析方法を図面を参照して詳細に説明する。
【0017】
〔第１の実施形態〕
本発明に係るトランスクリプトーム解析装置１は、図１に示すように、所定の目的変数データについて多数の遺伝子発現量データ（p次元、但しpは転写産物の数に相当する）を含むデータセットから第1～第mのデータセット（2≦m≦p-1）を生成するデータセット生成部２と、データセット生成部２で生成した第1～第mのデータセットのそれぞれに対して正則化法を適用して、目的変数データを目的変数とし遺伝子発現量データを説明変数とする第1～第mの予測式を算出する予測式算出部３と、予測式算出部３で算出した第1～第mの予測式に含まれる遺伝子発現量データに対応する遺伝子のリストを生成する遺伝子リスト生成部４とを備えている。また、トランスクリプトーム解析装置１は、遺伝子のアノテーション情報が格納された外部のデータベース５にアクセスできるものであっても良い。
本発明に係るトランスクリプトーム解析装置１は、図１に示すように、所定の目的変数データについて多数の遺伝子発現量データ（p次元、但しpは転写産物の数に相当する）を含むデータセットから第1～第mのデータセット（2≦m≦p-1）を生成するデータセット生成部２と、データセット生成部２で生成した第1～第mのデータセットのそれぞれに対して正則化法を適用して、目的変数データを目的変数とし遺伝子発現量データを説明変数とする第1～第mの予測式を算出する予測式算出部３と、予測式算出部３で算出した第1～第mの予測式に含まれる遺伝子発現量データに対応する遺伝子のリストを生成する遺伝子リスト生成部４とを備えている。また、トランスクリプトーム解析装置１は、遺伝子のアノテーション情報が格納された外部のデータベース５にアクセスできるものであっても良い。
【0018】
トランスクリプトーム解析装置１に入力するデータセットは、所定の目的変数データと、遺伝子発現量データ（p次元）とを含んでいる。ここで、目的変数データとは、量的形質或いは質的形質を含む表現型に関する数値データ、周辺環境条件の有無並びに程度に関する数値データ、解析対象生物に対する処理の有無並びに程度に関する数値データを含む意味である。より具体的に、目的変数データは、植物体といった解析対象生物の生育量に関するデータ（例えば、地上部重量、根部重量、葉面積、種子収量等）、解析対象生物に負荷するストレスに関するデータ（例えば、高温度処理時間、低温度処理時間、薬剤処理濃度、病害虫ストレス時間等）を挙げることができる。
【0019】
また、遺伝子発現量データとは、観察される転写産物について、発現量の相対量を示す数値データである。より具体的に、遺伝子発現量データとしては、市販されている遺伝子発現解析用（トランスクリプトーム解析用）マイクロアレイや、市場で提供されている遺伝子発現解析受託サービス等を利用して得られるマイクロアレイデータや、次世代シーケンス装置を用いた発現解析（RNA-Seq）を利用して得られたシーケンスデータ等を挙げることができる。特に、遺伝子発現量データとしては、次世代シーケンス装置を用いた発現解析（RNA-Seq）を利用して得られたシーケンスデータとすることが好ましい。次世代シーケンス装置を用いた発現解析（RNA-Seq）を利用して得られたシーケンスデータには、解析対象生物における転写産物が網羅されているからである。
【0020】
トランスクリプトーム解析装置１によれば、上記データセットを解析することで、目的変数データを説明できる遺伝子のリストを生成することができる。なお、遺伝子のリストとは、タンパク質をコードする狭義の遺伝子のリストに限定されず、非コーディング領域からの転写産物のリストも含む意味である。
【0021】
例えば、植物体の初期生育量（所定の期間の植物重量）を目的変数データとした場合、トランスクリプトーム解析装置１によれば、初期生育量を説明できる遺伝子リストを生成することができる。また、植物体に処理する薬剤濃度を目的変数データとした場合、トランスクリプトーム解析装置１によれば、処理する薬剤濃度に関連して発現する遺伝子のリストを生成することができる。さらに、サンプリング時の気温を目的変数データとした場合、トランスクリプトーム解析装置１によれば、生育温度に関連して発現する遺伝子のリストを生成することができる。
【0022】
【0023】
先ず、マイクロアレイ装置や次世代シーケンス装置等から出力された遺伝子発現量データと目的変数データとを入力する（ステップＳ１）。ここで、入力された遺伝子発現量データをp次元とし、p次元説明変数ベクトルx = {x1,……,xp}とする。また入力された目的変数をyとする。なお、本例では、p次元説明変数ベクトルxと目的変数yとからなるn組のデータセット（{(yi,xi)| i=1, ……,n}）が入力されるものとする。
【0024】
次に、データセット生成部２において、入力されたn組のデータセットに含まれる遺伝子発現量データをランダムにサンプリングすることで、p-1次元以下のサブデータセットを生成する（ステップＳ２）。なお、本ステップでは、初期値をm=1とするm番目のサブデータセットを生成する。言い換えると、生成する第mのサブデータセットは、入力されたデータセットに含まれる遺伝子発現量データをランダムに削減して、入力されたデータセットより少ない数の遺伝子発現量データを含むデータセットとして定義される。
【0025】
ここで、生成される第mのサブデータセットは、入力されたデータセットに含まれる遺伝子発現量データ（p次元）の一部であれば良く、例えば、p次元の遺伝子発現量データのうち５～９０％のデータとすることができ、５～７０％のデータとすることができ、５～５０％のデータとすることができ、１０～５０％のデータとすることができ、１０～２５％のデータとすることができ、１０～１５％のデータとすることができる。
【0026】
例えば、遺伝子発現量データ数が30000である場合（すなわちP=30000、30000個の転写産物）、データセット生成部２で生成する第mのサブデータセットは、ランダムに選択された1000～20000個の遺伝子発現量データを含むことができ、好ましくは1500～15000個の遺伝子発現量データを含むことができ、より好ましくは1500～7500個の遺伝子発現量データを含むことができ、更に好ましくは1500～4500個の遺伝子発現量データを含むことができる。
【0027】
次に、予測式算出部３において、データセット生成部２で生成した第mのサブデータセットに対して正則化項を有する回帰分析法を適用して、目的変数データを目的変数とし遺伝子発現量データを説明変数とする第mの予測式を算出する（ステップＳ３）。ここで、正則化項を有する回帰分析法とは、正則化回帰モデルとも呼ばれ、最小二乗法に制約（罰則）を付け加えて推定量を縮小させる解析法である。具体的に、正則化項を有する回帰分析法としては、LASSO回帰分析法、Ridge回帰分析法及びelastic net回帰分析法を挙げることができる。特に、LASSO回帰分析法を適用して予測式を算出することが好ましい。本ステップにおいて算出される予測式は、特にLASSO回帰分析法を適用した場合、目的変数を説明するときに重要でない遺伝子発現量データのパラメータを０とした予測式となる。
【0028】
なお、LASSO回帰分析法を適用して予測式を算出する際には、Friedman et al., Regularization Paths for Generalized Linear Models via Coordinate Descent, Journal of Statistical Software, January 2010, Volume 33, Issue 1を参照することができる。
【0029】
次に、遺伝子リスト生成部４において、予測式算出部３で算出した予測式に含まれる遺伝子発現量データを抽出し、抽出した遺伝子発現量データに対応する遺伝子を記録する（ステップＳ４）。すなわち、予測式は正則化項を有する回帰分析法により算出されているため、目的変数を説明するときに重要な遺伝子発現量データのみを抽出することができる。例えば、正則化項を有する回帰分析法としてLASSO回帰分析法を適用した場合、パラメータを０とした遺伝子発現量データ以外の遺伝子発現量データを抽出する。
【0030】
次に、ステップＳ５において、上記ステップＳ２～Ｓ４の処理を予め規定していた回数（m回）を繰り返したか判断する。例えば、繰り返し回数として10000回（m=10000）を予め規定していた場合、初期値を１とした第１のサブデータセットに対して上記ステップＳ２～Ｓ４の処理を実行した後、ステップＳ５においてm=1の値を10000とを比較し、ステップＳ６に進む。ステップＳ６においてm値を一つ増やし、ステップＳ２～Ｓ５をm値が10000になるまで繰り返す。
【0031】
以上のステップＳ２～Ｓ６をm回繰り返すことによって、ステップＳ１で入力されたn組のデータセットについて、第1～mのサブデータセットのそれぞれに対し第1～第mの予測式を算出することができる。
【0032】
次に、ステップＳ７では、遺伝子リスト生成部４において、第1～mの予測式についてステップＳ４で記録した遺伝子をリストとして出力する。遺伝子リスト生成部４で生成する遺伝子のリストは、特に限定されないが、遺伝子リスト生成部４で抽出した遺伝子発現量データに対応する遺伝子を列挙した形式でも良いし、抽出した遺伝子発現量データに対応する遺伝子と当該遺伝子の出現確率とを関連づけた形式でも良い。ここで出現確率とは、第1～第mの予測式に含まれる全ての遺伝子数に対する特定の遺伝子が含まれる回数として算出することができる。
【0033】
また、遺伝子リスト生成部４で生成する遺伝子のリストは、上述のように算出した出現確率が所定の値を超えた遺伝子のみを含む形式であっても良いし、上述のように算出した出現確率が高いものから順に列挙する形式であっても良い。
【0034】
一例として、遺伝子リスト生成部４で生成する遺伝子のリストの出力例を図３に示す。遺伝子リスト生成部４で生成する遺伝子のリストは、図３に示すように、転写産物毎に割り振られたIDと、転写産物が由来する遺伝子に関するシンボルと、転写産物毎に算出された出現確率と、目的変数データとの相関係数とを含んでいる。
【0035】
さらに、図１に示すように、トランスクリプトーム解析装置１が外部のデータベース５にアクセスし、リストに含まれる遺伝子のアノテーション情報等を検索し、得られたアノテーション情報等を遺伝子に関連づけた形式としても良い。また、トランスクリプトーム解析装置１は、検索したアノテーション情報に基づいて遺伝子を群分けし、群毎に遺伝子をリストとする形式としても良い。
【0036】
一例として、遺伝子リスト生成部４で生成する遺伝子のリストの出力例を図４に示す。遺伝子リスト生成部４で生成する遺伝子のリストは、図４に示すように、転写産物毎に割り振られたIDと、転写産物が由来する遺伝子に関するシンボルと、当該遺伝子シンボルで特定される遺伝子の機能に関する情報と、転写産物毎に算出された出現確率と、目的変数データとの相関係数とを含んでいる。
【0037】
図３及び/又は４に示した遺伝子リストによれば、所定の目的変数データに関して解析した結果として、当該目的変数データを説明できる遺伝子群を理解することができる。特にこれら遺伝子リストに上記出現確率が関連づけられている場合、リストに挙げられた各遺伝子について、当該出現確率に基づいて目的変数データとの関連性の強さを判断することができる。さらに、これら遺伝子リストにアノテーション情報が関連づけられている場合、リストに挙げられた各遺伝子について、当該アノテーション情報に基づいて目的変数データとの関連性について生物学的な意味を理解することができる。
【0038】
〔第２の実施形態〕
ところで、本発明に係るトランスクリプトーム解析装置及び解析方法は、上述した第１の実施形態に限定されず、図５及び６に示すように、所定の目的変数データに関して作成した遺伝子のリストを利用して、当該目的変数データに関する予測モデル式を作成するものであっても良い。なお、図５及び６に示すトランスクリプトーム解析装置１０及び解析方法において、図１及び２に示したトランスクリプトーム解析装置及び解析方法と同一の構成及び工程に対しては、図１及び２と同一の符号を付すことにより、その詳細な説明は省略する。
ところで、本発明に係るトランスクリプトーム解析装置及び解析方法は、上述した第１の実施形態に限定されず、図５及び６に示すように、所定の目的変数データに関して作成した遺伝子のリストを利用して、当該目的変数データに関する予測モデル式を作成するものであっても良い。なお、図５及び６に示すトランスクリプトーム解析装置１０及び解析方法において、図１及び２に示したトランスクリプトーム解析装置及び解析方法と同一の構成及び工程に対しては、図１及び２と同一の符号を付すことにより、その詳細な説明は省略する。
【0039】
図５に示したトランスクリプトーム解析装置１０は、遺伝子リスト生成部４で生成したリストに含まれる遺伝子に基づいて予測モデル式を生成する予測モデル式生成部１１を備えている。予測モデル式生成部１１を備えるトランスクリプトーム解析装置は、遺伝子リスト生成部４で生成する遺伝子のリスト（例えば図３及び４）に加えて、所定の目的変数データを説明する説明変数を含む予測モデル式を生成できる。
【0040】
トランスクリプトーム解析装置１０では、図６に示すように、上述した第１の実施形態と同様に、ステップＳ１～Ｓ６にて遺伝子リストを生成する。その後、トランスクリプトーム解析装置１０では、予測モデル式生成部１１において、リストに含まれる遺伝子について、ステップＳ１で入力したn個のデータセットから当該遺伝子に関する説明変数データ及び目的変数データをそれぞれ読み出す。そして、各遺伝子に関する目的変数y及び説明変数xの値を用いた重回帰分析や機械学習により所定の目的変数データを説明する予測モデル式を構築することができる。
【0041】
また、予測モデル式生成部１１では、遺伝子リスト生成部４で生成したリストに含まれる全ての遺伝子に関する目的変数y及び説明変数xの値を用いた重回帰分析や機械学習により所定の目的変数データを説明する予測モデル式を生成しても良いし、当該リストに含まれる一部の遺伝子に関する目的変数y及び説明変数xの値を用いた重回帰分析や機械学習により所定の目的変数データを説明する予測モデル式を生成しても良い。リストに含まれる一部の遺伝子としては、例えば、出現頻度の値が閾値を超える範囲の遺伝子することができ、所定のアノテーション情報が関連づけられている遺伝子としても良い。
【0042】
予測式を構築するための方法としては、特に限定されないが、例えば、LASSO回帰解析法、Ridge回帰解析法及びelastic net解析法などから選ばれる重回帰法や、Random forest法及びDeep learningなどの機械学習法を挙げることができる。
【0043】
一例として、予測モデル式生成部１１においてRandom forest法を適用して予測モデル式を作成することができる。このRandom forest法にて作成した予測モデル式は、所定の目的変数yを算出する決定木の形式のモデル式であり、遺伝子リスト生成部４で生成したリストに含まれる遺伝子の遺伝子発現量データxの関数として生成される。予測モデル式生成部１１で生成した予測モデル式によれば、所定の生物について取得した遺伝子発現量データに基づいて、当該生物に関して目的変数の予測値を算出することができる。
【0044】
ここで、遺伝子リスト生成部４で生成したリストに含まれる目的変数y（実測値）と、Random forest法を適用して作成した予測モデル式に基づいて算出した予測値との関係を図７に示す。図７に示すように、Random forest法を適用して作成した予測モデル式によれば、算出した予測値が実測値と非常に高い適合度を示すことがわかる。なお、図７に示したグラフは、後述の実施例に記載したデータを用いたものである。
【0045】
例えば植物の種子収量を目的変数として上記予測モデル式を得た場合、検査対象の植物から所得した遺伝子発現量データを用いることで、当該植物の種子収量を予測することができる。すなわち、検査対象の植物に関して栽培試験を経ずとも、次世代シーケンサーにより簡易に取得できる遺伝子発現量データから、上述した種子収量等の目的変数を推測することができる。
【0046】
以上のように、本実施の形態に示すトランスクリプトーム解析装置によれば、所定の目的変数について予測モデル式を作成することができる。作成した予測モデル式を用いることで、例えば図８に示すように、検査対象生物の特性評価システム２０を構築することができる。
【0047】
図８に示した特性評価システムは、本実施の形態に示すトランスクリプトーム解析装置で作成した所定の目的変数に関する予測モデル式を格納した記憶部２１と、検査対象生物に関する遺伝子発現データに基づいて目的変数を予測する予測部２２とを備えている。記憶部２１は、検査対象の生物毎に様々な目的変数について予測モデル式を格納している。例えば記憶部２１は、検査対象である植物について、地上部重量、根部重量、葉面積、種子収量、高温度処理時間、低温度処理時間、薬剤処理濃度、病害虫ストレス時間等の目的変数についてそれぞれ予測モデル式を格納することができる。
【0048】
予測部２２は、検査対象の生物に関する遺伝子発現量データが入力されると、記憶部２１に格納された予測モデル式のそれぞれに遺伝子発現量データを代入し、種々の目的変数について予測値を算出することができる。
【0049】
このように、特性評価システム２０は、様々な目的変数についてそれぞれ予測モデル式を記憶部２１に格納しておくことによって、検査対象生物の遺伝子発現量データに基づいて、これら目的変数に関する予測値を出力することができる。例えば、所定の植物について遺伝子発現量データが入力されると、地上部重量、根部重量、葉面積、種子収量、高温度処理時間、低温度処理時間、薬剤処理濃度、病害虫ストレス時間等の目的変数について一括して或いは選択した範囲で予測値を得ることができる。
【0050】
以上で説明した第１の実施形態及び第２の実施形態に係るトランスクリプトーム解析装置及び解析方法は、例えば、中央処理装置（CPU）、主記憶装置、補助記憶装置、出力装置及び入力装置を備えるコンピュータによって実現することができる。すなわち、例えば、目的変数データ及び遺伝子発現量データは、中央処理装置の制御のもと、入力装置を介して入力することができ、主記憶装置或いは補助記憶装置に記憶することができる。また、例えば、第mのサブデータセットは、中央処理装置の制御のもと、所定のアルゴリズムに従って生成することができる。さらに、第mのサブデータセットに基づく第mの予測式は、中央処理装置の制御のもと、所定のアルゴリズムに従って生成することができる。このように、以上で説明した第１の実施形態及び第２の実施形態に係るトランスクリプトーム解析装置及び解析方法は、中央処理装置の制御のもとで実現することができる。
【0051】
ただし、以上で説明した第１の実施形態及び第２の実施形態に係るトランスクリプトーム解析装置及び解析方法は、いわゆるクラウドコンピューティングにより実現することもできる。クラウドコンピューティングでは、例えば、クラウドサーバーに格納した目的変数データ及び遺伝子発現量データを利用することができ、また、生成した予測式や遺伝子リストをクラウドサーバーに格納することもできる。
【実施例】
【0052】
以下、実施例により本発明をより詳細に説明するが、本発明の技術的範囲はこれら実施例に限定されるものではない。
【0053】
〔実施例１〕
1．材料および方法
1－1．実験材料イネ系統と栽培条件
本実施例において、Ouu 365/Arroz da Terra//Ouu 365戻し交配自殖系統(BILs)は、Fukuda et al., 2014, Plant Production Science 17:41-46.に記述した系統を使用した。系統種子を50倍希釈次亜塩素酸で消毒し、水道水で3回洗浄したのち、30℃水中で2日間浸漬して発芽させた。1系統あたり24粒の発芽種子を水耕栽培用フローターに播種し(Fukuda et al., 2012, Plant Production Science 15:183-191.)、水耕栽培用溶液上で生育させた(Hayashi and Chino, 1986, Plant and Cell Physiology 27:1387-1393.)。水耕液は2日おきに作り替え、20℃、12時間明暗周期のグロースチャンバー内で14日間生育させた。
1．材料および方法
1－1．実験材料イネ系統と栽培条件
本実施例において、Ouu 365/Arroz da Terra//Ouu 365戻し交配自殖系統(BILs)は、Fukuda et al., 2014, Plant Production Science 17:41-46.に記述した系統を使用した。系統種子を50倍希釈次亜塩素酸で消毒し、水道水で3回洗浄したのち、30℃水中で2日間浸漬して発芽させた。1系統あたり24粒の発芽種子を水耕栽培用フローターに播種し(Fukuda et al., 2012, Plant Production Science 15:183-191.)、水耕栽培用溶液上で生育させた(Hayashi and Chino, 1986, Plant and Cell Physiology 27:1387-1393.)。水耕液は2日おきに作り替え、20℃、12時間明暗周期のグロースチャンバー内で14日間生育させた。
【0054】
1－2．方法
1－2－1．QTL解析
発芽14日後のBIL104系統と親2系統の苗を採取し、乾燥機で80℃2日間乾燥させた後、種子と根部分を取り除き、秤量した。実験は3反復のBiological replicatesにて行い、苗地上部乾物重量の平均値をQTL解析に用いた。BILの遺伝子型は124種のSSRマーカーを用いて解析し(Fukuda et al., 2014, Plant Production Science 17:41-46.)、MAPMAKER/EXP 3.0 (Lander et al., 1987, Genomics 1:174-181. doi:10.1016/0888-7543(87)90010-3)とQTL Cartographer 2.5 (Wang et al., 2010, Statistical Genetics & Bioinformatics, North Carolina State Universityにて提供)を用いてQTL解析を行った。
1－2－1．QTL解析
発芽14日後のBIL104系統と親2系統の苗を採取し、乾燥機で80℃2日間乾燥させた後、種子と根部分を取り除き、秤量した。実験は3反復のBiological replicatesにて行い、苗地上部乾物重量の平均値をQTL解析に用いた。BILの遺伝子型は124種のSSRマーカーを用いて解析し(Fukuda et al., 2014, Plant Production Science 17:41-46.)、MAPMAKER/EXP 3.0 (Lander et al., 1987, Genomics 1:174-181. doi:10.1016/0888-7543(87)90010-3)とQTL Cartographer 2.5 (Wang et al., 2010, Statistical Genetics & Bioinformatics, North Carolina State Universityにて提供)を用いてQTL解析を行った。
【0055】
1－2－2．RNAの単離とRNA-seq
親品種のOuu365とArroz da Terra、ならびにそれぞれ初期生育量の異なるBIL20系統を選出し、RNA-seq解析に用いた。発芽14日後の苗について、種子と根部分を取り除き、苗地上部の生鮮重量を測定した後、液体窒素中で凍結し、解析に用いるまで-80℃で保存した。RNeasy mini Kit (Qiagen社製)を用いてRNAを抽出した後、RNA-seq解析を行った。RNAの定量・定性を2100-Bioanalyzer (Agilent Technologies社製)を用いて行った後、TruSeq RNA LT Sample Preparation Kit v2 (Illumina Inc社製)を用いてシークエンス用ライブラリーを作成した。Illumina Hiseq 2000により、100bp, single-end readにて、ライブラリーのシークエンスを行った。シークエンス結果のFastqファイルをDDBJ Sequence Read Archive (DRA) 、 accession no. DRA006312に示した。
シークエンスデータはOryza sativa-Nipponbare-Reference-IRGSP-1.0 genome (Oryza sativa.IRGSP-1.0.24.dna.toplevel.fa.gz, ftp://ftp.ensemblgenomes.org/pub/release-24/plants/fasta/oryza_sativa/dna/) およびgene set (Oryza sativa.IRGSP-1.0.24.gtf.gz, ftp://ftp.ensemblgenomes.org/pub/release-24/plants/gtf/oryza_sativa/) を参照配列として、 TopHat2 (Kim et al., 2013, Genome Biology 14:13. doi:10.1186/gb-2013-14-4-r36; Trapnell et al., 2009, Bioinformatics 25:1105-1111. doi:10.1093/bioinformatics/btp120)を用いてマッピングを行った。各遺伝子の発現量について、 FPKM (Fragments Per Kilobase Million)値として算出した。
親品種のOuu365とArroz da Terra、ならびにそれぞれ初期生育量の異なるBIL20系統を選出し、RNA-seq解析に用いた。発芽14日後の苗について、種子と根部分を取り除き、苗地上部の生鮮重量を測定した後、液体窒素中で凍結し、解析に用いるまで-80℃で保存した。RNeasy mini Kit (Qiagen社製)を用いてRNAを抽出した後、RNA-seq解析を行った。RNAの定量・定性を2100-Bioanalyzer (Agilent Technologies社製)を用いて行った後、TruSeq RNA LT Sample Preparation Kit v2 (Illumina Inc社製)を用いてシークエンス用ライブラリーを作成した。Illumina Hiseq 2000により、100bp, single-end readにて、ライブラリーのシークエンスを行った。シークエンス結果のFastqファイルをDDBJ Sequence Read Archive (DRA) 、 accession no. DRA006312に示した。
シークエンスデータはOryza sativa-Nipponbare-Reference-IRGSP-1.0 genome (Oryza sativa.IRGSP-1.0.24.dna.toplevel.fa.gz, ftp://ftp.ensemblgenomes.org/pub/release-24/plants/fasta/oryza_sativa/dna/) およびgene set (Oryza sativa.IRGSP-1.0.24.gtf.gz, ftp://ftp.ensemblgenomes.org/pub/release-24/plants/gtf/oryza_sativa/) を参照配列として、 TopHat2 (Kim et al., 2013, Genome Biology 14:13. doi:10.1186/gb-2013-14-4-r36; Trapnell et al., 2009, Bioinformatics 25:1105-1111. doi:10.1093/bioinformatics/btp120)を用いてマッピングを行った。各遺伝子の発現量について、 FPKM (Fragments Per Kilobase Million)値として算出した。
【0056】
1－2－3．発現量バイオマーカーと遺伝子選出頻度の算出
苗地上部生鮮重を表す発現量バイオマーカーと、遺伝子の選出頻度について、以下の方法で算出した。発現量の平均値が0.01以上の遺伝子37043種について、以下のように解析に用いた。各遺伝子の発現量について、FPKM値に0.01を加えた後にLog2値に変換した。発現量バイオマーカーについて、LASSO法を用い、L1線形回帰モデルにより選出を行った(Tibshirani, 1996, Journal of the Royal Statistical Society Series B-Methodological 58:267-288. )。バイオマーカー遺伝子の選出頻度を計算するため、トランスクリプトームの部分集団(subset)を用いてのバイオマーカーの選出を繰り返し行った。37043遺伝子の中から10％の遺伝子をランダムに選択し、変数としてLASSO解析に用いた。インプットした変数の中から8遺伝子を、適切な、係数がゼロでない説明変数として選出し、発現量バイオマーカーとした。部分集団(subset)の選出と、発現量バイオマーカーの算出は10000回繰り返し行った。各遺伝子の選出頻度を、10000回のトライアルでバイオマーカーに使用された割合として決定した。解析はRのglmnet package (R Core Team, 2015, R: A language and environment for statistical computing. https://www.R-project.org/)を用いて行った。
苗地上部生鮮重を表す発現量バイオマーカーと、遺伝子の選出頻度について、以下の方法で算出した。発現量の平均値が0.01以上の遺伝子37043種について、以下のように解析に用いた。各遺伝子の発現量について、FPKM値に0.01を加えた後にLog2値に変換した。発現量バイオマーカーについて、LASSO法を用い、L1線形回帰モデルにより選出を行った(Tibshirani, 1996, Journal of the Royal Statistical Society Series B-Methodological 58:267-288. )。バイオマーカー遺伝子の選出頻度を計算するため、トランスクリプトームの部分集団(subset)を用いてのバイオマーカーの選出を繰り返し行った。37043遺伝子の中から10％の遺伝子をランダムに選択し、変数としてLASSO解析に用いた。インプットした変数の中から8遺伝子を、適切な、係数がゼロでない説明変数として選出し、発現量バイオマーカーとした。部分集団(subset)の選出と、発現量バイオマーカーの算出は10000回繰り返し行った。各遺伝子の選出頻度を、10000回のトライアルでバイオマーカーに使用された割合として決定した。解析はRのglmnet package (R Core Team, 2015, R: A language and environment for statistical computing. https://www.R-project.org/)を用いて行った。
【0057】
1－2－4．SG1遺伝子のシークエンシング
Ouu365とArroz da TerraのSG1遺伝子のコード領域、および上流-2108bpの領域について、PCRにより、以下のプライマーを用いて増幅した（5’-GGGACGTGATAACCGACTCA-3’（配列番号１）および5’-CCCCACTGTACGTTCTCTCC-3’（配列番号２））。PCR産物をillustra ExoProStar kitを用いて精製し、Fasmac社に送付してシークエンシングを行った。
翻訳開始点より-1948bp上流の1塩基置換について検出するため、以下のプライマーを用いてPCR差物を増幅し（5’-GGGACGTGATAACCGACTCA-3’（配列番号３）及び5’-TTCAGGTCACCTAGCCCATC-3’（配列番号４））、制限酵素HaeIIIにより切断を行った。Arroz da Terra型の配列GGCCはHaeIIIにより切断されるが、Ouu365型の配列AGCCは切断されなかった。
Ouu365とArroz da TerraのSG1遺伝子のコード領域、および上流-2108bpの領域について、PCRにより、以下のプライマーを用いて増幅した（5’-GGGACGTGATAACCGACTCA-3’（配列番号１）および5’-CCCCACTGTACGTTCTCTCC-3’（配列番号２））。PCR産物をillustra ExoProStar kitを用いて精製し、Fasmac社に送付してシークエンシングを行った。
翻訳開始点より-1948bp上流の1塩基置換について検出するため、以下のプライマーを用いてPCR差物を増幅し（5’-GGGACGTGATAACCGACTCA-3’（配列番号３）及び5’-TTCAGGTCACCTAGCCCATC-3’（配列番号４））、制限酵素HaeIIIにより切断を行った。Arroz da Terra型の配列GGCCはHaeIIIにより切断されるが、Ouu365型の配列AGCCは切断されなかった。
【0058】
1－2－5．定量リアルタイムPCR
苗地上部から、トータルRNAを、上記の通り抽出した。1μgのトータルRNAを用いてPrimeScript RT reagent Kit with gDNA Eraser (Takara Bio社)によって、cDNAの合成を行った。Thermal Cycler Dice Real Time System IIIを用い、SYBR Premix Ex TaqとプライマーセットOA045647 (Takara Bio社)によって、SG1 cDNA量の定量をリアルタイムPCRにより行った。リアルタイムPCRの測定は3反復のtechnical replicatesにより行った。SG1mRNAのコピー数の算出のため、Ouu365のｃDNAを鋳型として、SG1のPCR産物を以下のプライマーを用いて増幅し（5’-CGACCAGCTGATCTCCAA-G3’（配列番号５）及び5’-CATTTTTACTGGCCCTTCCA-3’（配列番号６））、リアルタイム定量PCRのスタンダードとして用いた。スタンダード用PCR産物は、Qubit fluorometer (Thermo Fisher Scientific社)を用いて定量を行い、その分子量からコピー数を算出した。SG1発現量(copies per ng RNA)はLog2値に変換した後、QTL解析に用いた。
苗地上部から、トータルRNAを、上記の通り抽出した。1μgのトータルRNAを用いてPrimeScript RT reagent Kit with gDNA Eraser (Takara Bio社)によって、cDNAの合成を行った。Thermal Cycler Dice Real Time System IIIを用い、SYBR Premix Ex TaqとプライマーセットOA045647 (Takara Bio社)によって、SG1 cDNA量の定量をリアルタイムPCRにより行った。リアルタイムPCRの測定は3反復のtechnical replicatesにより行った。SG1mRNAのコピー数の算出のため、Ouu365のｃDNAを鋳型として、SG1のPCR産物を以下のプライマーを用いて増幅し（5’-CGACCAGCTGATCTCCAA-G3’（配列番号５）及び5’-CATTTTTACTGGCCCTTCCA-3’（配列番号６））、リアルタイム定量PCRのスタンダードとして用いた。スタンダード用PCR産物は、Qubit fluorometer (Thermo Fisher Scientific社)を用いて定量を行い、その分子量からコピー数を算出した。SG1発現量(copies per ng RNA)はLog2値に変換した後、QTL解析に用いた。
【0059】
2．結果
2－1．戻し交雑自殖系統（BIL）のQTL解析
Arroz da TerraとOuu365の発芽14日後の地上部乾物重の平均はそれぞれ5.11mg、2.91mgであり、Arroz da Terraが有意に重かった（t-test, 5%水準）。BIL104系統の地上部乾物重量は、2.52から5.47mgの間に分布した（図９）。BIL104系統を用いて地上部乾物重のQTL解析を行った結果、第3，7及び10染色体上にArroz da Terra型で地上部乾物重を増加させるQTLが検出された(表１、図１０)。なお、図１０中、黒四角は、地上部乾物重を増加させるQTLの位置を示している。白抜き楕円は、SG1発現量を低下させるeQTLの位置を示している。
2－1．戻し交雑自殖系統（BIL）のQTL解析
Arroz da TerraとOuu365の発芽14日後の地上部乾物重の平均はそれぞれ5.11mg、2.91mgであり、Arroz da Terraが有意に重かった（t-test, 5%水準）。BIL104系統の地上部乾物重量は、2.52から5.47mgの間に分布した（図９）。BIL104系統を用いて地上部乾物重のQTL解析を行った結果、第3，7及び10染色体上にArroz da Terra型で地上部乾物重を増加させるQTLが検出された(表１、図１０)。なお、図１０中、黒四角は、地上部乾物重を増加させるQTLの位置を示している。白抜き楕円は、SG1発現量を低下させるeQTLの位置を示している。
【0060】
【表1】
【0061】
2－2．RNA-seq 解析とバイオマーカー遺伝子の選出
初期生育量と関連を持つ転写産物の探索のため、親品種2品種と、BILの中から異なる初期生育量を持つ20系統を用いて(図１１)、発芽14日後の苗地上部からRNAを抽出し、RNA-seq解析に用いた。なお、図１１中の白抜き三角形は、RNA-seq分析に使用したBIL系統それぞれについて地上部乾物重の平均値を示している。苗地上部生鮮重は、図１２に示した。
サンプル当たり平均41.6Mのリード数が得られ、96.1％にあたる40.0Mリード/サンプルがOs-Nipponbare-Reference-IRGSP-1.0 genome上にマッピングされた。遺伝子発現量はFPKM値（fragments per kilobase of coding sequence per million reads）として算出した。苗地上部生鮮重を表すバイオマーカーとなる遺伝子の選出頻度について、上記「1－2－3．発現量バイオマーカーと遺伝子選出頻度の算出」に示した方法を用いて以下のように決定した。全発現遺伝子の中から10％の遺伝子をランダムに選択して部分集団（subset）とし、LASSO解析を用いて部分集団内から8遺伝子を、苗地上部生鮮重を表す説明変数（発現量バイオマーカー）として算出した。部分集団(subset)の選出と、発現量バイオマーカーの算出を10000回繰り返し、各遺伝子が発現量バイオマーカーとして選出される頻度を決定した。高頻度で選出された遺伝子は、その発現量が苗地上部生鮮重と連動していることを示す。高頻度（１％以上の確立）で158遺伝子が選出された。これら選出された158遺伝子のリストを図１３に示した。選出された158遺伝子の発現量は、すべて地上部生鮮重と有意な相関を持っていた(5%水準)。
初期生育量と関連を持つ転写産物の探索のため、親品種2品種と、BILの中から異なる初期生育量を持つ20系統を用いて(図１１)、発芽14日後の苗地上部からRNAを抽出し、RNA-seq解析に用いた。なお、図１１中の白抜き三角形は、RNA-seq分析に使用したBIL系統それぞれについて地上部乾物重の平均値を示している。苗地上部生鮮重は、図１２に示した。
サンプル当たり平均41.6Mのリード数が得られ、96.1％にあたる40.0Mリード/サンプルがOs-Nipponbare-Reference-IRGSP-1.0 genome上にマッピングされた。遺伝子発現量はFPKM値（fragments per kilobase of coding sequence per million reads）として算出した。苗地上部生鮮重を表すバイオマーカーとなる遺伝子の選出頻度について、上記「1－2－3．発現量バイオマーカーと遺伝子選出頻度の算出」に示した方法を用いて以下のように決定した。全発現遺伝子の中から10％の遺伝子をランダムに選択して部分集団（subset）とし、LASSO解析を用いて部分集団内から8遺伝子を、苗地上部生鮮重を表す説明変数（発現量バイオマーカー）として算出した。部分集団(subset)の選出と、発現量バイオマーカーの算出を10000回繰り返し、各遺伝子が発現量バイオマーカーとして選出される頻度を決定した。高頻度で選出された遺伝子は、その発現量が苗地上部生鮮重と連動していることを示す。高頻度（１％以上の確立）で158遺伝子が選出された。これら選出された158遺伝子のリストを図１３に示した。選出された158遺伝子の発現量は、すべて地上部生鮮重と有意な相関を持っていた(5%水準)。
【0062】
2－3．苗地上部重QTL内に含まれる高頻度選出バイオマーカー遺伝子
選出された高頻度遺伝子と苗地上部重量QTLとを比較すると、第3、7及び10染色体上QTL内に含まれる遺伝子が、それぞれ5個、6個及び4個あった。そのうち第3染色体上遺伝子の中に、既存の低温発芽遺伝子、qLTG3-1が含まれていた(RAP ID: Os03g0103300, Fujino et al., 2008, Theoretical and Applied Genetics 108:794-799. doi:10.1007/s00122-003-1509-4)。RNA-seqに用いた系統のqLTG3-1発現量と地上部生鮮重との間には、有意な正の相関が見られた（図１４）。親品種のひとつArroz da Terraは、機能型のqLTG3-1遺伝子を有していることが報告されているが(Fujino and Iwata, 2011, Theoretical and Applied Genetics 123:1089-1097. doi:10.1007/s00122-011-1650-4)、もう一方の親品種Ouu365は、qLTG3-1遺伝子コード領域内に71bpの欠損を持ち、その機能を失っていることが確認されている(Fukuda et al., 2014, Plant Production Science 17:41-46)。RNA-seqに用いたBIL系統のqLTG3-1遺伝子型を調査した結果、Arroz da Terra型のqLTG3-1を持つ系統の地上部生鮮重とqLTG3-1発現量は、Ouu365型のqLTG3-1を持つ系統に比べ、有意に高かった(t-test, 1% level)。
選出された高頻度遺伝子と苗地上部重量QTLとを比較すると、第3、7及び10染色体上QTL内に含まれる遺伝子が、それぞれ5個、6個及び4個あった。そのうち第3染色体上遺伝子の中に、既存の低温発芽遺伝子、qLTG3-1が含まれていた(RAP ID: Os03g0103300, Fujino et al., 2008, Theoretical and Applied Genetics 108:794-799. doi:10.1007/s00122-003-1509-4)。RNA-seqに用いた系統のqLTG3-1発現量と地上部生鮮重との間には、有意な正の相関が見られた（図１４）。親品種のひとつArroz da Terraは、機能型のqLTG3-1遺伝子を有していることが報告されているが(Fujino and Iwata, 2011, Theoretical and Applied Genetics 123:1089-1097. doi:10.1007/s00122-011-1650-4)、もう一方の親品種Ouu365は、qLTG3-1遺伝子コード領域内に71bpの欠損を持ち、その機能を失っていることが確認されている(Fukuda et al., 2014, Plant Production Science 17:41-46)。RNA-seqに用いたBIL系統のqLTG3-1遺伝子型を調査した結果、Arroz da Terra型のqLTG3-1を持つ系統の地上部生鮮重とqLTG3-1発現量は、Ouu365型のqLTG3-1を持つ系統に比べ、有意に高かった(t-test, 1% level)。
【0063】
2－4．苗地上部重QTL外にある高頻度選出バイオマーカー遺伝子
苗地上部重量QTL内に含まれなかった高頻度選出遺伝子の中に、既存の組織伸長抑制遺伝子SG1(Short Grain 1, RAP ID: Os09g0459200, Nakagawa et al., 2012, Plant Physiology 158:1208-1219. doi:10.1104/pp.111.187567)が含まれていた。RNA-seqに用いた系統のSG1遺伝子発現量と、地上部生鮮重は、有意な負の相関を持っていた（図１５）。SG1は、過剰発現形質転換体において、植物ホルモンのブラシノステロイドへの応答性を低下させ、植物体を矮化させることが知られている(Nakagawa et al., 2012, Plant Physiology 158:1208-1219. doi:10.1104/pp.111.187567)。しかし、SG1が自然変異を持つかどうかは、今まで報告されていない。親品種のArroz da TerraとOuu365のSG1遺伝子の塩基配列を比較した結果、コード領域内に塩基置換や欠失・挿入変異は無かった。翻訳開始点の上流-1948bと-2038bの位置に単塩基置換があったが、RNA-seq解析に用いた系統のSG1遺伝子発現量は、この位置の遺伝子型によって差は見られなかった。
苗地上部重量QTL内に含まれなかった高頻度選出遺伝子の中に、既存の組織伸長抑制遺伝子SG1(Short Grain 1, RAP ID: Os09g0459200, Nakagawa et al., 2012, Plant Physiology 158:1208-1219. doi:10.1104/pp.111.187567)が含まれていた。RNA-seqに用いた系統のSG1遺伝子発現量と、地上部生鮮重は、有意な負の相関を持っていた（図１５）。SG1は、過剰発現形質転換体において、植物ホルモンのブラシノステロイドへの応答性を低下させ、植物体を矮化させることが知られている(Nakagawa et al., 2012, Plant Physiology 158:1208-1219. doi:10.1104/pp.111.187567)。しかし、SG1が自然変異を持つかどうかは、今まで報告されていない。親品種のArroz da TerraとOuu365のSG1遺伝子の塩基配列を比較した結果、コード領域内に塩基置換や欠失・挿入変異は無かった。翻訳開始点の上流-1948bと-2038bの位置に単塩基置換があったが、RNA-seq解析に用いた系統のSG1遺伝子発現量は、この位置の遺伝子型によって差は見られなかった。
【0064】
2－5．BIL104系統のSG1発現量の定量リアルタイムPCR解析
RNA-seq解析に用いた以外のBIL系統においても、SG1発現量と苗地上部重量とが相関を持つか確認するため、BIL104系統すべてと親品種について、定量リアルタイムPCRによるSG1発現量の測定を行った。その結果、BIL104系統と親品種のSG1発現量と地上部生鮮重は、有意な負の相関を示した（図１６）。翻訳開始点の上流-1948bの遺伝子型によるSG1の発現量の違いは見られなかった。SG1発現量に影響する染色体領域を調査するため、BIL104系統のSG1発現量を用いて発現量QTL解析(eQTL解析)を行った結果、第3染色体上と第7染色体上の2か所に、Arroz da Terra型でSG1発現量を低下させるeQTLが検出された(表２及び図１０)。このうち、作用力の強い第7染色体上のeQTLは、苗重量QTLと同位置にあった(図１０)。一方で、SG1遺伝子が存在する第9染色体上には、eQTLは検出されなかった。
RNA-seq解析に用いた以外のBIL系統においても、SG1発現量と苗地上部重量とが相関を持つか確認するため、BIL104系統すべてと親品種について、定量リアルタイムPCRによるSG1発現量の測定を行った。その結果、BIL104系統と親品種のSG1発現量と地上部生鮮重は、有意な負の相関を示した（図１６）。翻訳開始点の上流-1948bの遺伝子型によるSG1の発現量の違いは見られなかった。SG1発現量に影響する染色体領域を調査するため、BIL104系統のSG1発現量を用いて発現量QTL解析(eQTL解析)を行った結果、第3染色体上と第7染色体上の2か所に、Arroz da Terra型でSG1発現量を低下させるeQTLが検出された(表２及び図１０)。このうち、作用力の強い第7染色体上のeQTLは、苗重量QTLと同位置にあった(図１０)。一方で、SG1遺伝子が存在する第9染色体上には、eQTLは検出されなかった。
【0065】
【表2】
【0066】
本実施例によって、BIL20系統および親系統を用いたRNA-seq解析により、初期生育を表す指標となるバイオマーカー遺伝子候補が選出できることが明らかとなった。また、その中に、既存の組織伸長抑制遺伝子SG1が含まれていた。SG1が組織伸長抑制の作用を持つことは、activation-tagによる過剰発現形質転換体により確認されているが(Nakagawa et al., 2012, Plant Physiology 158:1208-1219. doi:10.1104/pp.111.187567)、自然状態において、SG1発現量に系統間に違いがあるかは不明であった。本実施例のトランスクリプトーム解析により、BIL系統のSG1発現量と苗地上部生鮮重が負の相関を持つことが明らかとなった(図１５)。さらに、定量リアルタイムPCR解析によって、RNA-seqに使用された系統のみでなく、104のBIL系統すべてにおいて、SG1発現量と苗地上部生鮮重とが負の相関を持つことが明らかとなり（図１６）、SG1が苗の初期生育量に影響していることが示唆された。これらの結果から、22系統のRNA-seqデータを用いた、本実施例のトランスクリプトーム解析は、初期生育に関わる転写産物を検出するために有効な手段であると考えられた。
【0067】
2－6．本実施例で示したトランスクリプトーム解析の有用性
転写産物の網羅解析（トランスクリプトーム解析）は、様々な形態的・生理的性質に影響する転写産物を検出できる強力な手段であるが、一方で、転写産物は多くの環境要因・遺伝要因の影響を複雑に受ける。そのため、特定の性質を表す発現量バイオマーカーを統計的に選出するには、ノイズを取り除くため、数百以上の多数のサンプル数を用いることが望ましいと考えられている。しかしながら、数百以上の多数のサンプル数を準備し、RNA-Seq等の遺伝子発現解析を行うことは困難な場合が多い。
本実施例に示したトランスクリプトーム解析においては、BIL20系統および親品種2系統の22系統という、比較的小さいサンプルサイズを用いて苗重量をあらわす発現量バイオマーカーの検出を試みた。その結果、本実施例に示したトランスクリプトーム解析によれば、候補バイオマーカーとして、qLTG3-1とSG1という、ゲノム変異を持つものと持たないもの、2種の既存の遺伝子を検出することができた。この結果より、本実施例に示したトランスクリプトーム解析は、比較的小さいサンプルサイズの解析でも、効果的に発現量バイオマーカーの選出を行える可能性を示している。
転写産物の網羅解析（トランスクリプトーム解析）は、様々な形態的・生理的性質に影響する転写産物を検出できる強力な手段であるが、一方で、転写産物は多くの環境要因・遺伝要因の影響を複雑に受ける。そのため、特定の性質を表す発現量バイオマーカーを統計的に選出するには、ノイズを取り除くため、数百以上の多数のサンプル数を用いることが望ましいと考えられている。しかしながら、数百以上の多数のサンプル数を準備し、RNA-Seq等の遺伝子発現解析を行うことは困難な場合が多い。
本実施例に示したトランスクリプトーム解析においては、BIL20系統および親品種2系統の22系統という、比較的小さいサンプルサイズを用いて苗重量をあらわす発現量バイオマーカーの検出を試みた。その結果、本実施例に示したトランスクリプトーム解析によれば、候補バイオマーカーとして、qLTG3-1とSG1という、ゲノム変異を持つものと持たないもの、2種の既存の遺伝子を検出することができた。この結果より、本実施例に示したトランスクリプトーム解析は、比較的小さいサンプルサイズの解析でも、効果的に発現量バイオマーカーの選出を行える可能性を示している。
【0068】
〔実施例２〕
本実施例では、実施例１で作成した高頻度遺伝子リスト（図１３）を使用して、苗地上部生鮮重の予測値を算出した。
本実施例では、実施例１で作成した高頻度遺伝子リスト（図１３）を使用して、苗地上部生鮮重の予測値を算出した。
【0069】
1．方法
実施例1で作成した高頻度遺伝子リスト（図１３）158遺伝子のうち上位100遺伝子の遺伝子発現量及び苗地上部生鮮重（図１２）を用いてrandom forest法（Breiman, L., 2001, Machine Learning 45: 5-32）により、遺伝子発現量から苗地上部生鮮重を予測した。random forestでは、これら100遺伝子に関する、実施例１で測定した発現量データと苗地上部生鮮重量を入力値として決定木の形式で予測モデル式を作成し、当該予測モデル式に基づいて上記100遺伝子に関する発現量データから予測値を算出するものである。
実施例1で作成した高頻度遺伝子リスト（図１３）158遺伝子のうち上位100遺伝子の遺伝子発現量及び苗地上部生鮮重（図１２）を用いてrandom forest法（Breiman, L., 2001, Machine Learning 45: 5-32）により、遺伝子発現量から苗地上部生鮮重を予測した。random forestでは、これら100遺伝子に関する、実施例１で測定した発現量データと苗地上部生鮮重量を入力値として決定木の形式で予測モデル式を作成し、当該予測モデル式に基づいて上記100遺伝子に関する発現量データから予測値を算出するものである。
【0070】
2．結果
５分割交差検証（cross validation）を20回繰り返し、苗地上部生鮮重の予測値を求めた。横軸を苗地上部生鮮重の実測値とし、縦軸を上記予測モデル式により算出された予測値（平均値）としてデータをプロットしたグラフを図１７に示した。図１７に示したデータについてR2（自由度調整済決定係数）を算出したところ0.8554となり、非常に高い適合度を示した。すなわち、実施例１で作成したリストに含まれる遺伝子に関する遺伝子発現量データ及び苗地上部生鮮重を用いて策した予測モデル式は、実際のデータに当てはまっていることを表しており、説明変数（遺伝子発現量データ）が目的変数（苗地上部生鮮重）をよく説明していると言える。
５分割交差検証（cross validation）を20回繰り返し、苗地上部生鮮重の予測値を求めた。横軸を苗地上部生鮮重の実測値とし、縦軸を上記予測モデル式により算出された予測値（平均値）としてデータをプロットしたグラフを図１７に示した。図１７に示したデータについてR2（自由度調整済決定係数）を算出したところ0.8554となり、非常に高い適合度を示した。すなわち、実施例１で作成したリストに含まれる遺伝子に関する遺伝子発現量データ及び苗地上部生鮮重を用いて策した予測モデル式は、実際のデータに当てはまっていることを表しており、説明変数（遺伝子発現量データ）が目的変数（苗地上部生鮮重）をよく説明していると言える。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13-1】
【図13-2】
【図13-3】
【図13-4】
【図13-5】
【図13-6】
【図13-7】
【図13-8】
【図13-9】
【図14】
【図15】
【図16】
【図17】
【配列表】