知財求人 - 知財ポータルサイト「IP Force」
▶ 中国水産科学研究院黄海水産研究所の特許一覧 ▶ 南方海洋科学与工程广東省実験室湛江の特許一覧
特表2024-521656東星斑の抗病性優良品種を育成するためのゲノム選択方法
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公表特許公報(A)
(11)【公表番号】
(43)【公表日】2024-06-04
(54)【発明の名称】東星斑の抗病性優良品種を育成するためのゲノム選択方法
(51)【国際特許分類】
A01K 67/02 20060101AFI20240528BHJP
C12N 15/09 20060101ALI20240528BHJP
C12Q 1/6888 20180101ALI20240528BHJP
G06Q 50/02 20240101ALI20240528BHJP
G16B 30/00 20190101ALI20240528BHJP
C12N 15/11 20060101ALN20240528BHJP
C12Q 1/6809 20180101ALN20240528BHJP
C12Q 1/6869 20180101ALN20240528BHJP
【ＦＩ】
A01K67/02
C12N15/09 Z
C12Q1/6888 Z
G06Q50/02
G16B30/00
C12N15/11 Z ZNA
C12Q1/6809 Z
C12Q1/6869 Z
【審査請求】有
【予備審査請求】未請求
(21)【出願番号】P 2023570117
(86)(22)【出願日】2022-06-02
(85)【翻訳文提出日】2023-11-10
(86)【国際出願番号】 CN2022096720
(87)【国際公開番号】W WO2023103303
(87)【国際公開日】2023-06-15
(31)【優先権主張番号】202111478938.0
(32)【優先日】2021-12-06
(33)【優先権主張国・地域又は機関】CN
(81)【指定国・地域】
【公序良俗違反の表示】
（特許庁注：以下のものは登録商標）
１．Ｌｉｎｕｘ
(71)【出願人】
【識別番号】520057601
【氏名又は名称】中国水産科学研究院黄海水産研究所
(71)【出願人】
【識別番号】523427593
【氏名又は名称】南方海洋科学与工程广東省実験室湛江
【氏名又は名称原語表記】ＧＵＡＮＧＤＯＮＧ ＰＲＯＶＩＮＣＩＡＬ ＬＡＢＯＲＡＴＯＲＹ ＯＦ ＳＯＵＴＨＥＲＮ ＭＡＲＩＮＥ ＳＣＩＥＮＣＥ ＡＮＤ ＥＮＧＩＮＥＥＲＩＮＧ （ＺＨＡＮＪＩＡＮＧ）
【住所又は居所原語表記】Ｎｏ．１ Ｗｅｎｔｉ Ｒｏａｄ， Ｘｉａｓｈａｎ Ｄｉｓｔｒｉｃｔ（Ｌｉｂｒａｒｙ） Ｚｈａｎｊｉａｎｇ， Ｇｕａｎｇｄｏｎｇ ５２４０００， Ｃｈｉｎａ
(74)【代理人】
【識別番号】100178434
【弁理士】
【氏名又は名称】李 じゅん
(72)【発明者】
【氏名】陳 松林
(72)【発明者】
【氏名】盧 昇
(72)【発明者】
【氏名】劉 洋
(72)【発明者】
【氏名】周 茜
(72)【発明者】
【氏名】王 磊
(72)【発明者】
【氏名】朱 春華
(72)【発明者】
【氏名】張 天時
(72)【発明者】
【氏名】陳 亜東
(72)【発明者】
【氏名】徐 文騰
【テーマコード（参考）】
4B063
5L050
【Ｆターム（参考）】
4B063QA07
4B063QA13
4B063QA20
4B063QQ02
4B063QQ42
4B063QR32
4B063QR72
4B063QX02
5L050CC01
(57)【要約】
【課題】本発明は、東星斑の抗病性優良品種を育成するためのゲノム選択方法に関するものである。
【解決手段】前記方法において、異なる地理由来の非家族系グループの東星斑の幼魚を使用して抗病性の参照グループを生成し、抗病性参照グループを生成するための東星斑の幼魚は、抗病性のスクリーニング後に取得されたものである。抗病性の表現型を有するが参照グループには含まれていない個体を検証グループとして生成して、本発明で提供されたＳＮＰ部位によりゲノム推定育種価を計算する効果を検証する。これらのＳＮＰ部位には東星斑の抗病性性状に関連するＳＮＰ部位が含まれる。本発明より提供されるＳＮＰ部位によって候補グループのゲノム推定育種価を計算し、ゲノム推定育種価が高い個体を抗病性が強いと認めて、東星斑の抗病性優良品種を育成するために利用する。本発明により提供されるゲノム選択技術に基づく東星斑の抗病性優良品種を育成するための方法は、抗病性の強い親をスクリーニングするために利用することができる。選択された親は、東星斑の選択育種に使用することができる。これにより、本発明は東星斑の繁殖グループの生存率を向上させるための効率的な技術的な手段を提供した。
【解決手段】前記方法において、異なる地理由来の非家族系グループの東星斑の幼魚を使用して抗病性の参照グループを生成し、抗病性参照グループを生成するための東星斑の幼魚は、抗病性のスクリーニング後に取得されたものである。抗病性の表現型を有するが参照グループには含まれていない個体を検証グループとして生成して、本発明で提供されたＳＮＰ部位によりゲノム推定育種価を計算する効果を検証する。これらのＳＮＰ部位には東星斑の抗病性性状に関連するＳＮＰ部位が含まれる。本発明より提供されるＳＮＰ部位によって候補グループのゲノム推定育種価を計算し、ゲノム推定育種価が高い個体を抗病性が強いと認めて、東星斑の抗病性優良品種を育成するために利用する。本発明により提供されるゲノム選択技術に基づく東星斑の抗病性優良品種を育成するための方法は、抗病性の強い親をスクリーニングするために利用することができる。選択された親は、東星斑の選択育種に使用することができる。これにより、本発明は東星斑の繁殖グループの生存率を向上させるための効率的な技術的な手段を提供した。
【特許請求の範囲】
【請求項1】
ゲノム選択に基づく東星斑の抗病性優良品種を育成するための方法であって、１）異なる地理由来の東星斑の幼魚により東星斑の抗病性の参照グループを生成し、且つ、
方法の選択効果を検証するため、抗病性の表現型を有するが参照グループには含まれていない個体により検証グループを生成するステップと、
２）参照グループと検証グループに対して全ゲノムリシーケンシングを行い、ゲノム範囲内の高品質の一塩基多型（ＳＮＰ）の部位を捕獲するステップと、
３）抗病性性状ゲノム推定育種価ＧＥＢＶを予測する精度がＳＮＰ数に応じて変化する傾向に関して、ゲノム最良線形不偏予測Ｇ－ＢＬＵＰとＢａｙｅｓＣπとの２つの方法を比較し、これにより、候補グループのＧＥＢＶを計算するためのＳＮＰ数と方法を選択するステップと、
４）ステップ３）で取得した結果に基づいて、東星斑の候補グループのＧＥＢＶを計算するための、東星斑の抗病性性状に関連するＳＮＰ部位が含まれているＳＮＰマーカーを選択するステップと、
５）ステップ３）で選択されたゲノム選択方法と、ステップ４）で選出したＳＮＰマーカーとを利用して、候補グループのＧＥＢＶを計算し、ＧＥＢＶが高い個体の抗病性が強いと認めて子孫を繁殖する親として使用するステップとを含む、
ことを特徴とするゲノム選択に基づく東星斑の抗病性優良品種を育成するための方法。
【請求項2】
前記ステップ１）における、抗病性の参照グループを生成するための東星斑の幼魚は、抗病性のスクリーニングを行って得られたものである、ことを特徴とする請求項１に記載の方法。
【請求項3】
前記方法では、検証グループを利用して、ステップ３）で選択されたゲノム選択方法とステップ４）で選出されたＳＮＰマーカーの予測効果を検証する、ことを特徴とする請求項１に記載の方法。
【請求項4】
前記ＳＮＰマーカーは、ＳＥＤ ＩＤ ＮＯ：１～ＳＥＤ ＩＤ ＮＯ：１００００のいずれかの配列の第３６位と第１０７位に位置する、ことを特徴とする請求項１に記載の方法。
【請求項5】
東星斑の抗病性性状に関連する前記ＳＮＰ部位は、ＳＥＤ ＩＤ ＮＯ：１～ＳＥＤ ＩＤ ＮＯ：１５のいずれかの配列の第３６位と第１０７位に位置する、ことを特徴とする請求項１に記載のゲノム選択方法。
【請求項6】
ＳＮＰマーカーセットであって、上記ＳＮＰマーカーセットにおけるＳＮＰ部位は、ＳＥＤ ＩＤ ＮＯ：１～ＳＥＤ ＩＤ ＮＯ：１００００のいずれかの配列の第３６位と第１０７位に位置する、ことを特徴とするＳＮＰマーカーセット。
【請求項7】
方法であって、東星斑の抗病性優良品種のスクリーニングにおいて、請求項６に記載のＳＮＰマーカーセットの利用する方法。【発明の詳細な説明】
【技術分野】
【0001】
本発明は、水産遺伝品種育成の技術分野に属し、特に、東星斑の抗病性優良品種を育成するためのゲノム選択方法に関する。
【背景技術】
【0002】
東星斑（トンシンバーン）の学名はニセスジハタ（Ｐｌｅｃｔｒｏｐｏｍｕｓ ｌｅｏｐａｒｄｕｓ）であり、スズキ目、ハタ科、ハタ亜科、スヂハタ属に属する。東星斑は肉質が緻密で美味しく、全身が鮮やかな赤色で、経済的価値や観賞価値が高く、貴重な海洋養殖魚類である。市場での需要が高いため、乱獲により野生個体数が大幅に減少され、東星斑は市場で希少な魚となっている。近年、東星斑の人工養殖技術の進歩と発展につれて、養殖範囲と規模は徐々に拡大されている。しかしながら、高密度な養殖モデルの普及に加え、東星斑の養殖環境の悪化や生殖質の劣化により、養殖過程で病害が頻繁に発生され、多大な経済的損失を引き起こした。その中でも、ビブリオ・ハーヴェイ（Ｖｉｂｒｉｏ ｈａｒｖｅｙｉ）による「腐体病」が養殖業者にとって大きな問題となっている。例えば、２０１１年に中国海南省でビブリオ・ハーヴェイを優勢菌とする細菌性疾患が発生し、発病した魚が大量に死亡し、尾部潰瘍や全身の白化などの症状が現れた。通常、養殖農家は病気に対処する際は予防に重点を置き、病気発生後はプール全体にオキシテトラサイクリン、エンロフロキサシン、ノルフロキサシンなどの抗生物質を撒いて治療する。しかし、抗生物質の大規模な使用は水産物の品質に影響を与えるだけでなく、食品の安全性の問題を引き起こし、市販魚の価値に影響を与え、消費者に潜在的な損害をもたらす。したがって、東星斑の抗病性優良品種を育成し、東星斑養殖グループ自体の抗病性を向上させることが非常に必要である。
【0003】
繁殖能力のある個体をいかに正確に選び出すかは、優良品種を育種するための核心的な問題である。魚類の抗病性育種において、通常、大規模に確立された家系と人工感染実験を利用して抗病性表現型を取得した後、最良線形不偏予測（ＰＢＬＵＰ）から推定された推定育種価（ＥＢＶ）によって選択育種する。ただし、病原体の垂直感染を避けるため、人工感染実験で生存された個体を繁殖親として使用するのではなく、ＥＢＶの高い家系から未感染の健康な個体を選択して繁殖させるのが一般的である。したがって、選択効率が高くなく、毎年得られる遺伝的進歩は限られている。選択育種された魚類が家系を確立できない場合、表現型や経験に基づいて選択することになり、選択育種の効率が大幅に低下してしまう。現在、東星班の人工繁殖は主に自然産卵や自然受精に依存しており、家族系を確立できない。そのため、東星斑のような家系確立が難しい魚類では有効な選択育種方法が不足している。本発明は、非家族系グループを利用して東星斑の抗病性優良品種を育成するための方法を詳細に記載し、東星斑の優良品種を育成するための技術方法を提供し、家族系確立が難しい魚類の選択育種のために新たなアイデアを提供することを目的としている。
【発明の概要】
【0004】
本発明の目的は、非家族系グループに基づく東星斑の抗病性優良品種を育成するためのゲノム選択方法を確立し、従来の選択育種における家系確立の必要性と選択精度の低さの問題を克服し、東星斑の抗病性優良品種を育成するための分子育種方法を提供し、東星斑養殖グループの抗病性を向上させ、東星斑の抗病性優良品種の育種プロセスを促進することにある。
【0005】
本発明により提供される東星斑の抗病性優良品種を育成するためのゲノム選択方法は、
１）異なる地理由来の東星斑の幼魚を使用して、東星斑の抗病性参照グループを生成し、
さらに、本発明の上記方法による選択効果の検証を行うために、抗病性表現型を有するが参照グループには含まれていない一部の個体を検証グループとして生成し、
ここで、抗病性参照グループを生成するための東星斑の幼魚は、抗病性のスクリーニングを行った後に取得されるものである、ステップと、
２）参照グループと検証グループに対して全ゲノムリシーケンシングを行い、ゲノム範囲内の高品質の一塩基多型（ＳＮＰ）部位を捕獲するステップと、
３）抗病性性状ゲノム推定育種価（ＧＥＢＶ）を予測する精度がＳＮＰ数に応じて変化する傾向に関して、ゲノム最良線形不偏予測（Ｇ－ＢＬＵＰ）とＢａｙｅｓＣπとの２つの方法を比較し、これにより、候補グループのＧＥＢＶを計算するためのＳＮＰ数と方法を選択するステップと、
４）ステップ３）で取得した結果に基づいて、東星斑の候補グループのＧＥＢＶを計算するための、東星斑の抗病性性状に関連するＳＮＰ部位が含まれているＳＮＰマーカーを選択するステップと、
５）検証グループを利用してステップ３）で選択されたゲノム選択方法とステップ４）で選び出されたＳＮＰマーカーとの予測効果を検証し、生存された検証個体の平均ＧＥＢＶが死亡個体の平均ＧＥＢＶよりも大きい場合、ＧＥＢＶと個体の抗病性との間に正の相関があることを示し、ＧＥＢＶが高い個体が抗病性も強いため、ＧＥＢＶを利用して候補の親を選択して、抗病性が強い子孫を育種することができるステップと、
６）ステップ３）で選択されたゲノム選択方法と、ステップ４）で選出したＳＮＰマーカーとを利用して、候補グループのＧＥＢＶを計算し、ＧＥＢＶが高い個体の抗病性が強いと認めて子孫を繁殖する親として使用するステップと、を含む。
１）異なる地理由来の東星斑の幼魚を使用して、東星斑の抗病性参照グループを生成し、
さらに、本発明の上記方法による選択効果の検証を行うために、抗病性表現型を有するが参照グループには含まれていない一部の個体を検証グループとして生成し、
ここで、抗病性参照グループを生成するための東星斑の幼魚は、抗病性のスクリーニングを行った後に取得されるものである、ステップと、
２）参照グループと検証グループに対して全ゲノムリシーケンシングを行い、ゲノム範囲内の高品質の一塩基多型（ＳＮＰ）部位を捕獲するステップと、
３）抗病性性状ゲノム推定育種価（ＧＥＢＶ）を予測する精度がＳＮＰ数に応じて変化する傾向に関して、ゲノム最良線形不偏予測（Ｇ－ＢＬＵＰ）とＢａｙｅｓＣπとの２つの方法を比較し、これにより、候補グループのＧＥＢＶを計算するためのＳＮＰ数と方法を選択するステップと、
４）ステップ３）で取得した結果に基づいて、東星斑の候補グループのＧＥＢＶを計算するための、東星斑の抗病性性状に関連するＳＮＰ部位が含まれているＳＮＰマーカーを選択するステップと、
５）検証グループを利用してステップ３）で選択されたゲノム選択方法とステップ４）で選び出されたＳＮＰマーカーとの予測効果を検証し、生存された検証個体の平均ＧＥＢＶが死亡個体の平均ＧＥＢＶよりも大きい場合、ＧＥＢＶと個体の抗病性との間に正の相関があることを示し、ＧＥＢＶが高い個体が抗病性も強いため、ＧＥＢＶを利用して候補の親を選択して、抗病性が強い子孫を育種することができるステップと、
６）ステップ３）で選択されたゲノム選択方法と、ステップ４）で選出したＳＮＰマーカーとを利用して、候補グループのＧＥＢＶを計算し、ＧＥＢＶが高い個体の抗病性が強いと認めて子孫を繁殖する親として使用するステップと、を含む。
【0006】
前記ＳＮＰ部位は、合計２０，０００（２０ｋ）個があり、ＳＥＤ ＩＤ ＮＯ：１～ＳＥＤ ＩＤ ＮＯ：１００００のいずれかの配列の第３６位と第１０７位に位置し、部位前後３５ｂｐ以内の配列はフランキング配列（flanking sequence）である。
【0007】
前記ＳＮＰ部位には、東星斑の抗病性に関連する３０個のＳＮＰ部位が含まれ、ＳＥＤ ＩＤ ＮＯ：１～ＳＥＤ ＩＤ ＮＯ：１５のいずれかの配列の第３６位と第１０７位に位置し、部位前後３５ｂｐ以内の配列はフランキング配列である。
【0008】
本発明により提供されるゲノム選択技術に基づく東星斑の抗病性優良品種を育成するための方法と２０，０００個のＳＮＰ部位は、抗病性の強い親をスクリーニングするために利用することができる。選択された親は、東星斑の選択育種に直接使用することができ、これにより、本発明は、東星斑繁殖グループの生存率を向上させるために効率的な技術的な手段を提供した。
【図面の簡単な説明】
【0009】
【図1】東星斑参照グループ及び検証グループの共通部位欠失率の箱ひげ図である。
【図2】Ｇ－ＢＬＵＰとＢａｙｅｓＣπとにより東星斑の抗病性性状ＧＥＢＶを予測する精度の変化を示すグラフである。
【図3】検証グループＧＥＢＶを推定するための２０ｋのＳＮＰの密度分布図である。
【図4】検証グループにおける生存個体と死亡個体の抗病性ＧＥＢＶ平均値の比較図である。
【図5】候補グループにおける抗病性が高い個体と感染され易い個体の抗病性ＧＥＢＶ平均値の比較図である。
【発明を実施するための形態】
【0010】
以下、本発明の抗ビブリオ・ハーヴェイの東星斑養殖グループを育成するための方法について、詳細に説明する。
【実施例1】
【0011】
実施例１：東星斑抗病性参照グループを生成する
中国海南省と山東省の多数の大規模な養殖場から東星斑の幼魚を収集し、収集された幼魚を人工感染実験が可能な場所へ輸送し、ビブリオ・ハーヴェイ菌液を腹腔内注射することによりこれらの幼魚を人為的に感染させる。菌液を注入した時刻を開始点（０時）とし、一度注入後の幼魚の生存状況を８時間ごとに観察し、適時に死亡個体を除去し、個体情報（全長、体重、死亡時刻など）を記録し、尾鰭の鰭条を切り取り、無水エタノール中に保存する。実験は合計５日間行われ、終了後、生存個体の尾鰭の鰭条を採集し、無水エタノール中に保存し、個体出所、全長、体重、死亡時刻などのデータを記録する。その中から一部の個体を選択して全ゲノムリシーケンシングを行い、東星斑の抗ビブリオ・ハーヴェイ参照グループを生成する（表１）。人工感染実験における被験された東星斑の幼魚の生存状態を抗病性表現型とし、その後の分析に使用する。抗病性表現型を二値性状として定義する。即ち、０は感染実験で死亡した幼魚を表し、１は実験期間で生存した幼魚を表す。
中国海南省と山東省の多数の大規模な養殖場から東星斑の幼魚を収集し、収集された幼魚を人工感染実験が可能な場所へ輸送し、ビブリオ・ハーヴェイ菌液を腹腔内注射することによりこれらの幼魚を人為的に感染させる。菌液を注入した時刻を開始点（０時）とし、一度注入後の幼魚の生存状況を８時間ごとに観察し、適時に死亡個体を除去し、個体情報（全長、体重、死亡時刻など）を記録し、尾鰭の鰭条を切り取り、無水エタノール中に保存する。実験は合計５日間行われ、終了後、生存個体の尾鰭の鰭条を採集し、無水エタノール中に保存し、個体出所、全長、体重、死亡時刻などのデータを記録する。その中から一部の個体を選択して全ゲノムリシーケンシングを行い、東星斑の抗ビブリオ・ハーヴェイ参照グループを生成する（表１）。人工感染実験における被験された東星斑の幼魚の生存状態を抗病性表現型とし、その後の分析に使用する。抗病性表現型を二値性状として定義する。即ち、０は感染実験で死亡した幼魚を表し、１は実験期間で生存した幼魚を表す。
【0012】
実験では、採集された幼魚を６つのグループに分け、注射後、それぞれ同じ仕様の繊維強化プラスチック（ＦＲＰ）水槽に投入する。中国山東省産からの東星斑の幼魚の平均体重と平均全長は、いずれも中国海南省産からの幼魚より高かった。一元配置分散検定では、６つの実験グループ間で感染生存率に顕著な差がないことを示した（ｐ＝０．０５）。感染後の生存率範囲は３１．８５％～４３．４１％であり、合計５３４匹の東星斑の幼魚を選択して参照グループを生成した（表１）。
【0013】
【表1】
【実施例2】
【0014】
実施例２：東星斑抗病性参照グループのゲノムリシーケンシングと高品質ＳＮＰのスクリーニング
１、参照グループの全ゲノムリシーケンシング及び変異の検出
ＤＮＡ抽出キット（ＴＩＡＮＧＥＮ、北京）により提供される標準手順を使用し、東星斑抗病性参照グループのゲノムＤＮＡを抽出及び精製した。参照グループのゲノムＤＮＡに加えて、本発明の上記方法の効果を検証するために、抗病性表現型を有する２９８匹の東星斑のゲノムＤＮＡも抽出した。精製されたゲノムＤＮＡを使用してイルミナ（ｌｌｕｍｉｎａ）ペアエンドライブラリを生成した後、ＩｌｌｕｍｉｎａＨｉＳｅｑ２０００シーケンシングプラットフォームを使用してシーケンシングし、低品質のシーケンシングリード（ｒｅａｄｓ）を除外した後、ソフトウェアＢＷＡを使用してこれらのリードを東星斑の参照ゲノムと比較した。ソフトウェアＧＡＴＫを使用し、変異情報（ＳＮＰ及びＩＮＤＥＬ）を検出した。最終的に、参照グループにおいて、４９８個の個体のシーケンシングが成功し、個体のシーケンシングデータ量は１０Ｇで、合計３４，０２６，７１９個の変異が取得され、染色体上の変異は２８，７８９，３８７個で、非染色体上の変異は５，２３７，３３２個で、個体平均シーケンシング深度は７．６６Ｘで、変異結果はファイル「ｖａｒｉａｎｃｅｓ．Ｒｅｆ．ｖｃｆ．ｇｚ」に保存される。検証グループにおいて、２９８個の個体のシーケンシングが成功し、個体のシーケンシングデータ量は５Ｇで、合計１２，４５２，４８３個の変異が取得され、染色体上の変異は１１，３５５，３７３個で、非染色体上の変異は１，０９７，１１０個で、個体平均シーケンシング深度は３．６６Ｘで、変異結果はファイル「ｖａｒｉａｎｃｅｓ．Ｖａｌ．ｖｃｆ．ｇｚ」に保存される。
１、参照グループの全ゲノムリシーケンシング及び変異の検出
ＤＮＡ抽出キット（ＴＩＡＮＧＥＮ、北京）により提供される標準手順を使用し、東星斑抗病性参照グループのゲノムＤＮＡを抽出及び精製した。参照グループのゲノムＤＮＡに加えて、本発明の上記方法の効果を検証するために、抗病性表現型を有する２９８匹の東星斑のゲノムＤＮＡも抽出した。精製されたゲノムＤＮＡを使用してイルミナ（ｌｌｕｍｉｎａ）ペアエンドライブラリを生成した後、ＩｌｌｕｍｉｎａＨｉＳｅｑ２０００シーケンシングプラットフォームを使用してシーケンシングし、低品質のシーケンシングリード（ｒｅａｄｓ）を除外した後、ソフトウェアＢＷＡを使用してこれらのリードを東星斑の参照ゲノムと比較した。ソフトウェアＧＡＴＫを使用し、変異情報（ＳＮＰ及びＩＮＤＥＬ）を検出した。最終的に、参照グループにおいて、４９８個の個体のシーケンシングが成功し、個体のシーケンシングデータ量は１０Ｇで、合計３４，０２６，７１９個の変異が取得され、染色体上の変異は２８，７８９，３８７個で、非染色体上の変異は５，２３７，３３２個で、個体平均シーケンシング深度は７．６６Ｘで、変異結果はファイル「ｖａｒｉａｎｃｅｓ．Ｒｅｆ．ｖｃｆ．ｇｚ」に保存される。検証グループにおいて、２９８個の個体のシーケンシングが成功し、個体のシーケンシングデータ量は５Ｇで、合計１２，４５２，４８３個の変異が取得され、染色体上の変異は１１，３５５，３７３個で、非染色体上の変異は１，０９７，１１０個で、個体平均シーケンシング深度は３．６６Ｘで、変異結果はファイル「ｖａｒｉａｎｃｅｓ．Ｖａｌ．ｖｃｆ．ｇｚ」に保存される。
【0015】
２、高品質ＳＮＰ変異部位をスクリーニングする
以下の手順に従って、上述の変異をフィルタリングし、高品質ＳＮＰ部位を保留した後、欠失している部位情報を充填する。Ｌｉｎｕｘ環境で次のコマンドを実行する：
以下の手順に従って、上述の変異をフィルタリングし、高品質ＳＮＰ部位を保留した後、欠失している部位情報を充填する。Ｌｉｎｕｘ環境で次のコマンドを実行する：
【0016】
#ＶＣＦｔｏｏｌｓを使用して染色体上のＳＮＰ部位を出力する
vcftools --gzvcf variances.Ref.vcf.gz
--remove-indels
--chr chr1 --chr chr2 --chr chr3 --chr chr4 --chr chr5 --chr chr6 --chr chr7
--chr chr8 --chr chr9 --chr chr10 --chr chr11 --chr chr12 --chr chr13
--chr chr14 --chr chr15 --chr chr16 --chr chr17--chr chr18 --chr chr19
--chr chr20 --chr chr21 --chr chr22 --chr chr23 --chr chr24
--recode-INFO-all --recode --stdout | bgzip -c -f > raw.snps.Ref.vcf.gz
vcftools --gzvcf variances.Val.vcf.gz
--remove-indels
--chr chr1 --chr chr2 --chr chr3 --chr chr4 --chr chr5 --chr chr6 --chr chr7
--chr chr8 --chr chr9 --chr chr10 --chr chr11 --chr chr12 --chr chr13
--chr chr14 --chr chr15 --chr chr16 --chr chr17--chr chr18 --chr chr19
--chr chr20 --chr chr21 --chr chr22 --chr chr23 --chr chr24
--recode-INFO-all --recode --stdout | bgzip -c -f > raw.snps.Val.vcf.gz
vcftools --gzvcf variances.Ref.vcf.gz
--remove-indels
--chr chr1 --chr chr2 --chr chr3 --chr chr4 --chr chr5 --chr chr6 --chr chr7
--chr chr8 --chr chr9 --chr chr10 --chr chr11 --chr chr12 --chr chr13
--chr chr14 --chr chr15 --chr chr16 --chr chr17--chr chr18 --chr chr19
--chr chr20 --chr chr21 --chr chr22 --chr chr23 --chr chr24
--recode-INFO-all --recode --stdout | bgzip -c -f > raw.snps.Ref.vcf.gz
vcftools --gzvcf variances.Val.vcf.gz
--remove-indels
--chr chr1 --chr chr2 --chr chr3 --chr chr4 --chr chr5 --chr chr6 --chr chr7
--chr chr8 --chr chr9 --chr chr10 --chr chr11 --chr chr12 --chr chr13
--chr chr14 --chr chr15 --chr chr16 --chr chr17--chr chr18 --chr chr19
--chr chr20 --chr chr21 --chr chr22 --chr chr23 --chr chr24
--recode-INFO-all --recode --stdout | bgzip -c -f > raw.snps.Val.vcf.gz
【0017】
実行後、参照グループと検証グループの染色体上のＳＮＰ部位は、それぞれファイル「ｒａｗ．ｓｎｐｓ．Ｒｅｆ．ｖｃｆ．ｇｚ」と「ｒａｗ．ｓｎｐｓ．Ｖａｌ．ｖｃｆ．ｇｚ」に保存される。
【0018】
#ＰＬＩＮＫ２を使用して上述のＶＣＦファイルを読み取り、二対立遺伝子部位を保留する
plink2 --vcf raw.snps.Ref.vcf.gz --autosome-num 24
--max-alleles 2 --min-alleles 2
--make-bpgen --out raw.Ref
plink2 --vcf raw.snps.Val.vcf.gz --autosome-num 24
--max-alleles 2 --min-alleles 2
--make-bpgen --out raw.Val
plink2 --vcf raw.snps.Ref.vcf.gz --autosome-num 24
--max-alleles 2 --min-alleles 2
--make-bpgen --out raw.Ref
plink2 --vcf raw.snps.Val.vcf.gz --autosome-num 24
--max-alleles 2 --min-alleles 2
--make-bpgen --out raw.Val
【0019】
#Ｒを使用して「ｒａｗ．Ｒｅｆ．ｂｉｍ」ファイルと「ｒａｗ．Ｖａｌ．ｂｉｍ」ファイルを修正する
library(data.table)
bim <- fread(“raw.Ref.bim”)
bim$V2 <- paste(“rs”, bim$V1, ‘:’, bim$V4, sep = “”)
fwrite(bim, “raw.Ref.bim”, sep = “＼t”, col.names = F, quote = F)
bim <- fread(“raw.Val.bim”)
bim$V2 <- paste(“rs”, bim$V1, ‘:’, bim$V4, sep = “”)
fwrite(bim, “raw.Can.bim”, sep = “＼t”, col.names = F, quote = F)
library(data.table)
bim <- fread(“raw.Ref.bim”)
bim$V2 <- paste(“rs”, bim$V1, ‘:’, bim$V4, sep = “”)
fwrite(bim, “raw.Ref.bim”, sep = “＼t”, col.names = F, quote = F)
bim <- fread(“raw.Val.bim”)
bim$V2 <- paste(“rs”, bim$V1, ‘:’, bim$V4, sep = “”)
fwrite(bim, “raw.Can.bim”, sep = “＼t”, col.names = F, quote = F)
【0020】
#ＰＬＩＮＫ２を使用して各部位の欠失率をカウントする
plink2 --bpfile raw.Ref --autosome-num 24 --missing --out Ref
plink2 --bpfile raw.Val --autosome-num 24 --missing --out Val
plink2 --bpfile raw.Ref --autosome-num 24 --missing --out Ref
plink2 --bpfile raw.Val --autosome-num 24 --missing --out Val
【0021】
#Ｒを使用して２つのファイル内の共通部位の欠失率をカウントし、欠失率が０．０５未満の部位を出力する
library(data.table)
vmiss.ref <- fread(“Ref.vmiss”)
vmiss.val <- fread(“Val.vmiss”)
comm.loci <- intersect(vmiss.ref$ID, vmiss.val$ID)
ref.comm <- vmiss.ref[ID %in% comm.loci]
val.comm <- vmiss.val[ID %in% comm.loci]
new.vmiss <- data.table(SNP = comm.loci, MISSING = (ref.comm$MISSING_CT + val.comm$MISSING_CT), OBS = (ref.comm$OBS_CT + val.comm$OBS_CT))
new.vmiss$F_MISS <- round(new.vmiss$MISSING / new.vmiss$OBS, 3)
kept.loci <- new.vmiss[F_MISS < 0.05, SNP]
write.table(kept.loci, “kept.loci.txt”, sep = “＼t”, col.names = F, row.names = F, quote = F)
実行後、ファイル「ｋｅｐｔ．ｌｏｃｉ．ｔｘｔ」には、２つのファイル中で欠失率が０．０５未満の部位が保存される（図１）。
library(data.table)
vmiss.ref <- fread(“Ref.vmiss”)
vmiss.val <- fread(“Val.vmiss”)
comm.loci <- intersect(vmiss.ref$ID, vmiss.val$ID)
ref.comm <- vmiss.ref[ID %in% comm.loci]
val.comm <- vmiss.val[ID %in% comm.loci]
new.vmiss <- data.table(SNP = comm.loci, MISSING = (ref.comm$MISSING_CT + val.comm$MISSING_CT), OBS = (ref.comm$OBS_CT + val.comm$OBS_CT))
new.vmiss$F_MISS <- round(new.vmiss$MISSING / new.vmiss$OBS, 3)
kept.loci <- new.vmiss[F_MISS < 0.05, SNP]
write.table(kept.loci, “kept.loci.txt”, sep = “＼t”, col.names = F, row.names = F, quote = F)
実行後、ファイル「ｋｅｐｔ．ｌｏｃｉ．ｔｘｔ」には、２つのファイル中で欠失率が０．０５未満の部位が保存される（図１）。
【0022】
#２つのファイルを合併する
##ＰＬＩＮＫ２を使用してＶＣＦファイルを生成する
plink2 --bpfile raw.Ref --autosome-num 24 --extract kept.loci.txt
--export vcf-4.2 --out snps.ref
plink2 --bpfile raw.Val --autosome-num 24 --extract kept.loci.txt
--export vcf-4.2 --out snps.val
##ＶＣＦファイルを圧縮する
bgzip -c -f snps.ref.vcf > snps.ref.vcf.gz
bgzip -c -f snps.val.vcf > snps.val.vcf.gz
##ｂｃｆｔｏｏｌｓを使用してインデックスを作成する
bcftools index -t snps.ref.vcf.gz
bcftools index -t snps.val.vcf.gz
##ｂｃｆｔｏｏｌｓを使用してＶＣＦファイルを合併する
bcftools merge -m snps -O snps.ref.vcf.gz snps.val.vcf.gz -o snps.merged.vcf.gz
##ＰＬＩＮＫ２を使用してＶＣＦファイルを生成する
plink2 --bpfile raw.Ref --autosome-num 24 --extract kept.loci.txt
--export vcf-4.2 --out snps.ref
plink2 --bpfile raw.Val --autosome-num 24 --extract kept.loci.txt
--export vcf-4.2 --out snps.val
##ＶＣＦファイルを圧縮する
bgzip -c -f snps.ref.vcf > snps.ref.vcf.gz
bgzip -c -f snps.val.vcf > snps.val.vcf.gz
##ｂｃｆｔｏｏｌｓを使用してインデックスを作成する
bcftools index -t snps.ref.vcf.gz
bcftools index -t snps.val.vcf.gz
##ｂｃｆｔｏｏｌｓを使用してＶＣＦファイルを合併する
bcftools merge -m snps -O snps.ref.vcf.gz snps.val.vcf.gz -o snps.merged.vcf.gz
【0023】
#ＰＬＩＮＫ２を使用し、合併後のＳＮＰ部位について品質管理を行う
##合併後のファイルを読み取る
plink2 --vcf snps.merged.vcf.gz --autosome-num 24
--max-alleles 2 --min-alleles 2
--make-bpgen --out snps.merged
##品質管理
plink2 --bpfile snps.merged --autosome-num 24
--maf 0.05 --hwe 0.001
--make-bpgen --out snps.qc
##品質管理後のデータをＶＣＦに変換する
plink2 --bpfile snps.qc --autosome-num 24 --export vcf --out snps4impute
##ＶＣＦファイルを圧縮する
bgzip -c -f snps4impute.vcf > snps4impute.vcf.gz
##欠失しているＳＮＰ部位を充填する
java -Xmx20g -jar beagle.05Apr21.9b7.jar gt=snps4impute.vcf.gz ne=500 iterations=30 burnin=10 out=snps.imputed
実行後、合計１，２１３，４９６個の高品質ＳＮＰがその後の分析のために残され、充填された遺伝子型はファイル「ｓｎｐｓ．ｉｍｐｕｔｅｄ．ｖｃｆ．ｇｚ」に保存される。
##合併後のファイルを読み取る
plink2 --vcf snps.merged.vcf.gz --autosome-num 24
--max-alleles 2 --min-alleles 2
--make-bpgen --out snps.merged
##品質管理
plink2 --bpfile snps.merged --autosome-num 24
--maf 0.05 --hwe 0.001
--make-bpgen --out snps.qc
##品質管理後のデータをＶＣＦに変換する
plink2 --bpfile snps.qc --autosome-num 24 --export vcf --out snps4impute
##ＶＣＦファイルを圧縮する
bgzip -c -f snps4impute.vcf > snps4impute.vcf.gz
##欠失しているＳＮＰ部位を充填する
java -Xmx20g -jar beagle.05Apr21.9b7.jar gt=snps4impute.vcf.gz ne=500 iterations=30 burnin=10 out=snps.imputed
実行後、合計１，２１３，４９６個の高品質ＳＮＰがその後の分析のために残され、充填された遺伝子型はファイル「ｓｎｐｓ．ｉｍｐｕｔｅｄ．ｖｃｆ．ｇｚ」に保存される。
【実施例3】
【0024】
実施例３：２つのゲノム選択方法で東星斑の抗病性ＧＥＢＶの精度を予測する
東星斑の抗ビブリオ・ハーヴェイ病ＧＥＢＶを推定するための方法及び適切なＳＮＰマーカー数を検討するために、本発明は、生成された参照グループの表現型及び遺伝子型データを利用し、５倍交差検証法によって、SNP数の減少に応じてＧ－ＢＬＵＰとＢａｙｅｓＣπが東星斑の抗ビブリオ・ハーヴェイＧＥＢＶを予測する精度の変化傾向を比較した。予測精度は被検者動作特性曲線下面積（ＡＵＣ）を評価指標とし、合計２５回計算し、２５回の平均値で最終精度を示す。計算前にＲパッケージ：ｄａｔａ．ｔａｂｌｅ、ＡＳＲｅｍｌ、ＢＧＬＲ、Ｒｃｐｐ、ｐＲＯＣ及びＰａｒａｌｌｅｌをインストールする必要がある。
東星斑の抗ビブリオ・ハーヴェイ病ＧＥＢＶを推定するための方法及び適切なＳＮＰマーカー数を検討するために、本発明は、生成された参照グループの表現型及び遺伝子型データを利用し、５倍交差検証法によって、SNP数の減少に応じてＧ－ＢＬＵＰとＢａｙｅｓＣπが東星斑の抗ビブリオ・ハーヴェイＧＥＢＶを予測する精度の変化傾向を比較した。予測精度は被検者動作特性曲線下面積（ＡＵＣ）を評価指標とし、合計２５回計算し、２５回の平均値で最終精度を示す。計算前にＲパッケージ：ｄａｔａ．ｔａｂｌｅ、ＡＳＲｅｍｌ、ＢＧＬＲ、Ｒｃｐｐ、ｐＲＯＣ及びＰａｒａｌｌｅｌをインストールする必要がある。
【0025】
１、東星斑抗病性参照グループから異なる数のＳＮＰセット（ＳＮＰサブセット）を抽出する
実施例２で取得された高品質ＳＮＰから異なる数のＳＮＰサブセットを抽出する。抗病性参照グループの個体番号は、ファイル「ｉｎｄｉｖ．ｒｅｆ．ｔｘｔ」に保存される。Ｌｉｎｕｘ環境で次のコマンドを実行し、ＳＮＰサブセットを抽出する：
#ＶＣＦファイルを読み込み、関連データを準備する
##ＰＬＩＮＫ２を使用してＶＣＦファイルを読み取り、品質管理を再度実行する
plink2 --vcf snps.imputed.vcf.gz --autosome-num 24 --maf 0.05 --make-bpgen --out snps.imputed
##ＰＬＩＮＫ２を使用して対立遺伝子頻度をカウントする
plink2 --bpfile snps.imputed --autosome-num 24 --freq --out snps.imputed
##Ｒで次のコードを実行し、遺伝子型の出力に必要なファイルを準備する
library(data.table)
freq <- fread(“snps.imputed.afreq”)
alleles <- rbind(freq[ALT_FREQS <= 0.5, c(“ID”, “ALT”)], freq[ALT_FREQS > 0.5, c(“ID”, “REF”)], use.names = F)
write.table(alleles, “export.A.txt”, sep = “＼t”, col.names = F, quote = F)
実施例２で取得された高品質ＳＮＰから異なる数のＳＮＰサブセットを抽出する。抗病性参照グループの個体番号は、ファイル「ｉｎｄｉｖ．ｒｅｆ．ｔｘｔ」に保存される。Ｌｉｎｕｘ環境で次のコマンドを実行し、ＳＮＰサブセットを抽出する：
#ＶＣＦファイルを読み込み、関連データを準備する
##ＰＬＩＮＫ２を使用してＶＣＦファイルを読み取り、品質管理を再度実行する
plink2 --vcf snps.imputed.vcf.gz --autosome-num 24 --maf 0.05 --make-bpgen --out snps.imputed
##ＰＬＩＮＫ２を使用して対立遺伝子頻度をカウントする
plink2 --bpfile snps.imputed --autosome-num 24 --freq --out snps.imputed
##Ｒで次のコードを実行し、遺伝子型の出力に必要なファイルを準備する
library(data.table)
freq <- fread(“snps.imputed.afreq”)
alleles <- rbind(freq[ALT_FREQS <= 0.5, c(“ID”, “ALT”)], freq[ALT_FREQS > 0.5, c(“ID”, “REF”)], use.names = F)
write.table(alleles, “export.A.txt”, sep = “＼t”, col.names = F, quote = F)
【0026】
実行後、１，２１１，２５９個の高品質ＳＮＰと７９６個の個体が残され、その後の分析に使用される。その中で、４９８個の個体が参照グループで、２９８個の個体が検証グループである。
【0027】
#Ｒを使用し、参照グループから異なるＳＮＰ数を含むサブセットを抽出し、バイナリ形式で保存する。ＳＮＰ数はそれぞれ５００、７００、１ｋ、２ｋ、８ｋ、１０ｋ、２０ｋ、５０ｋ、８０ｋ、１００ｋである。
##ＰＬＩＮＫ２を使用し、０／１／２モードでエンコードされた遺伝子型ＡＡ／Ａａ／ａａファイルを生成する
plink2 --bpfile snps.imputed --autosome-num 24 --export-allele export.A.txt --keep indiv.ref.txt --export A --out geno.ref
##個体番号を出力する
awk ‘{print $2}’ geno.ref.raw | sed ‘1d’ > iid.ref.txt
##Ｒで次のスクリプトを実行してＳＮＰサブセットを抽出し、バイナリ形式で保存する
library(data.table)
Den.snp <- c(5e2, 7e2, 1e3, 2e3, 8e3, 1e4, 2e4, 5e4, 8e4, 1e5)
geno <- fread(paste(“geno.ref.raw”))
geno[, c(1: 6):= NULL]
geno <- as.matrix(geno)
for (i in 1: length(Den.snp)) {
coor <- floor(seq(from = 1, to = ncol(geno), length.out = Den.snp[i]))
M <- geno[, coor]
save(M, file = paste(“ref.snp.s”, i, “.bin”, sep = “”))
}
実行後、すべてのＳＮＰサブセットは、後で簡単に読み取ることができるようにバイナリ形式で保存される。
##ＰＬＩＮＫ２を使用し、０／１／２モードでエンコードされた遺伝子型ＡＡ／Ａａ／ａａファイルを生成する
plink2 --bpfile snps.imputed --autosome-num 24 --export-allele export.A.txt --keep indiv.ref.txt --export A --out geno.ref
##個体番号を出力する
awk ‘{print $2}’ geno.ref.raw | sed ‘1d’ > iid.ref.txt
##Ｒで次のスクリプトを実行してＳＮＰサブセットを抽出し、バイナリ形式で保存する
library(data.table)
Den.snp <- c(5e2, 7e2, 1e3, 2e3, 8e3, 1e4, 2e4, 5e4, 8e4, 1e5)
geno <- fread(paste(“geno.ref.raw”))
geno[, c(1: 6):= NULL]
geno <- as.matrix(geno)
for (i in 1: length(Den.snp)) {
coor <- floor(seq(from = 1, to = ncol(geno), length.out = Den.snp[i]))
M <- geno[, coor]
save(M, file = paste(“ref.snp.s”, i, “.bin”, sep = “”))
}
実行後、すべてのＳＮＰサブセットは、後で簡単に読み取ることができるようにバイナリ形式で保存される。
【0028】
２、Ｇ－ＢＬＵＰの予測精度を評価する
ＡＳＲｒｅｍｌ－Ｒで、次のモデルを使用して表現型データをフィッティングさせ、Ｇ－ＢＬＵＰの予測精度を評価する。
yijk=Φ(Ori.Locationi+b.weightj+ak)
上述のモデルで、yijkは表現型値であり、場所ｉから採集された個体ｋを体重ｊのときに感染させた場合の表現型を示す。Ori.Locationiは固定効果であり、被検個体の採集場所を示す。b.weightjは共変量であり、感染実験時の被検個体の体重を示す。akは個体のランダムな相加効果である。Φは標準正規分布の累積関数を示す。Ｌｉｎｕｘ環境で次のコマンドを実行する。
##Ｒで次のスクリプトを実行し、３列Ｇ行列の逆行列を生成する
library(data.table)
library(Rcpp)
sourceCpp(“matrix_multiplication.cpp”)
ginv <- function(invG) {
Ginv <- data.frame(row = rep(1 : nrow(invG), nrow(invG)),
column = rep(1 : nrow(invG), each = nrow(invG)),
value = as.numeric(invG),
lower.mat = as.logical(lower.tri(invG, diag = T)))
Ginv <- Ginv[Ginv$lower.mat == T, c("row", "column", "value")]
Ginv <- Ginv[order(Ginv$row, Ginv$column), ]
return(Ginv)
}
No.subsets <- 10
IID <- fread(“iid.ref.txt”, header = F)$V1
A22 <- diag(1, length(IID))
w <- 0.95
for (i in 1: No.subsets) {
load(paste(“ref.snp.s”, i, “.bin”, sep = “”))
pi <- colSums(M) / (2 * nrow(M))
P <- matrix((2 * pi), nrow = nrow(M), ncol = ncol(M), byrow = T)
Z <- as.matrix(M - P)
sum2pq <- 2 * (sum(pi * (1- pi)))
Zt <- t(Z)
ZZt <- eigenMatMult(Z, Zt)
G <- (ZZt / sum2pq)
Ga <- (w * G) + ((1 - w) * A22)
invG <- solve(Ga)
Ginv <- ginv(invG)
attr(Ginv, “rowNames”) <- IID
attr(Ginv, “INVERSE”) <- TRUE
save(Ginv, file = paste(“Ginv.ref.s”, i, “.bin”, sep = “”))
}
##実行前に、次のスクリプトをファイル「ｍａｔｒｉｘ＿ｍｕｌｔｉｐｌｉｃａｔｉｏｎ．ｃｐｐ」に保存し、作業パスに配置する。
// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]
#include <RcppArmadillo.h>
#include <RcppEigen.h>
// [[Rcpp::export]]
SEXP armaMatMult(arma::mat A, arma::mat B){
arma::mat C = A * B;
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMatMult(Eigen::MatrixXd A, Eigen::MatrixXd B){
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMapMatMult(const Eigen::Map<Eigen::MatrixXd> A, Eigen::Map<Eigen::MatrixXd> B){
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
##Ｒで次のスクリプトを実行し、５倍交差検証に必要なファイル「ｃｏｏｒ４ｃｖ．ｒｅｆ．ｃｓｖ」を生成する。
library(data.table)
fold <- 5
No.cv <- 5
No.indiv <- 500
res <- matrix(nrow = (No.indiv / fold), ncol = (fold * No.cv))
for (i in 1: fold) {
coor <- data.table(seq = c(1: No.indiv), or = rnorm(No.indiv))
setorder(coor, -or)
res[, (i * 5 - 4): (i * 5)] <- coor$seq
}
for (i in 1: ncol(res))
res[, i] <- sort(res[ ,i])
res[which(res == 499)] <- NA
res[which(res == 500)] <- NA
write.table("coor4cv.ref.csv", sep = ",", col.names = F, row.names = F, quote = F)
##Ｒで次のスクリプトを実行し、Ｇ－ＢＬＵＰの予測精度を評価する。その中で、抗病性表現型変数の名前はsur.statusであり、Ori.Locationとb.weightは固定項目とし、それぞれ採集場所と体重である。
library(data.table)
library(pROC)
library(asreml)s
No.subset <- 10
ref.pheno <- asreml.read.table(“ref.pheno4cv.csv”, sep = “,”, header = T)
cv.coor <- as.matrix(fread(“coor4cv.ref.csv”))
cv.res <- matrix(NA, nrow = (ncol(cv.coor) + 1), ncol = No.subset)
colnames(cv.res) <- paste(‘s’, 1: ncol(cv.coor), sep = “”)
rownames(cv.res) <- c(paste(“CV”, 1: ncol(cv.coor), sep = “”), “Average”)
for (i in 1: No.subset) {
load(paste(“Ginv.ref.s”, i, “.bin”, sep = “”))
for (j in 1: ncol(cv.coor)) {
pheno4cv <- ref.pheno[-na.omit(cv.coor[, j]), ]
gblup <- asreml(sur.status ~ Ori.Location + b.weight,
random = ~vm(IID, Ginv),
dense = ~vm(IID, Ginv),
family = asr_binomial(“probit”),
data = pheno4cv,
maxiter = 100)
y4val <- ref.pheno$sur.status[na.omit(cv.coor[, j])]
gebv4val <- coef(gblup)$random[na.omit(cv.coor[, j])]
cv.res[j, i] <- roc(y4val, gebv4val)$auc
}
}
cv.res [nrow(cv.res), ] <- colMeans(cv.res, na.rm = T)
write.table(cv.res, “Ref.cv.results.gblup.csv”, sep = “,”, row.names = F, quote = F)
実行後、ファイル「Ｒｅｆ．ｃｖ．ｒｅｓｕｌｔｓ．ｇｂｌｕｐ．ｃｓｖ」の最終行情報が、即ちSNP数の変化に応じてＧ－ＢＬＵＰが抗病性ＧＥＢＶを予測する精度である。
ＡＳＲｒｅｍｌ－Ｒで、次のモデルを使用して表現型データをフィッティングさせ、Ｇ－ＢＬＵＰの予測精度を評価する。
yijk=Φ(Ori.Locationi+b.weightj+ak)
上述のモデルで、yijkは表現型値であり、場所ｉから採集された個体ｋを体重ｊのときに感染させた場合の表現型を示す。Ori.Locationiは固定効果であり、被検個体の採集場所を示す。b.weightjは共変量であり、感染実験時の被検個体の体重を示す。akは個体のランダムな相加効果である。Φは標準正規分布の累積関数を示す。Ｌｉｎｕｘ環境で次のコマンドを実行する。
##Ｒで次のスクリプトを実行し、３列Ｇ行列の逆行列を生成する
library(data.table)
library(Rcpp)
sourceCpp(“matrix_multiplication.cpp”)
ginv <- function(invG) {
Ginv <- data.frame(row = rep(1 : nrow(invG), nrow(invG)),
column = rep(1 : nrow(invG), each = nrow(invG)),
value = as.numeric(invG),
lower.mat = as.logical(lower.tri(invG, diag = T)))
Ginv <- Ginv[Ginv$lower.mat == T, c("row", "column", "value")]
Ginv <- Ginv[order(Ginv$row, Ginv$column), ]
return(Ginv)
}
No.subsets <- 10
IID <- fread(“iid.ref.txt”, header = F)$V1
A22 <- diag(1, length(IID))
w <- 0.95
for (i in 1: No.subsets) {
load(paste(“ref.snp.s”, i, “.bin”, sep = “”))
pi <- colSums(M) / (2 * nrow(M))
P <- matrix((2 * pi), nrow = nrow(M), ncol = ncol(M), byrow = T)
Z <- as.matrix(M - P)
sum2pq <- 2 * (sum(pi * (1- pi)))
Zt <- t(Z)
ZZt <- eigenMatMult(Z, Zt)
G <- (ZZt / sum2pq)
Ga <- (w * G) + ((1 - w) * A22)
invG <- solve(Ga)
Ginv <- ginv(invG)
attr(Ginv, “rowNames”) <- IID
attr(Ginv, “INVERSE”) <- TRUE
save(Ginv, file = paste(“Ginv.ref.s”, i, “.bin”, sep = “”))
}
##実行前に、次のスクリプトをファイル「ｍａｔｒｉｘ＿ｍｕｌｔｉｐｌｉｃａｔｉｏｎ．ｃｐｐ」に保存し、作業パスに配置する。
// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]
#include <RcppArmadillo.h>
#include <RcppEigen.h>
// [[Rcpp::export]]
SEXP armaMatMult(arma::mat A, arma::mat B){
arma::mat C = A * B;
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMatMult(Eigen::MatrixXd A, Eigen::MatrixXd B){
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMapMatMult(const Eigen::Map<Eigen::MatrixXd> A, Eigen::Map<Eigen::MatrixXd> B){
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
##Ｒで次のスクリプトを実行し、５倍交差検証に必要なファイル「ｃｏｏｒ４ｃｖ．ｒｅｆ．ｃｓｖ」を生成する。
library(data.table)
fold <- 5
No.cv <- 5
No.indiv <- 500
res <- matrix(nrow = (No.indiv / fold), ncol = (fold * No.cv))
for (i in 1: fold) {
coor <- data.table(seq = c(1: No.indiv), or = rnorm(No.indiv))
setorder(coor, -or)
res[, (i * 5 - 4): (i * 5)] <- coor$seq
}
for (i in 1: ncol(res))
res[, i] <- sort(res[ ,i])
res[which(res == 499)] <- NA
res[which(res == 500)] <- NA
write.table("coor4cv.ref.csv", sep = ",", col.names = F, row.names = F, quote = F)
##Ｒで次のスクリプトを実行し、Ｇ－ＢＬＵＰの予測精度を評価する。その中で、抗病性表現型変数の名前はsur.statusであり、Ori.Locationとb.weightは固定項目とし、それぞれ採集場所と体重である。
library(data.table)
library(pROC)
library(asreml)s
No.subset <- 10
ref.pheno <- asreml.read.table(“ref.pheno4cv.csv”, sep = “,”, header = T)
cv.coor <- as.matrix(fread(“coor4cv.ref.csv”))
cv.res <- matrix(NA, nrow = (ncol(cv.coor) + 1), ncol = No.subset)
colnames(cv.res) <- paste(‘s’, 1: ncol(cv.coor), sep = “”)
rownames(cv.res) <- c(paste(“CV”, 1: ncol(cv.coor), sep = “”), “Average”)
for (i in 1: No.subset) {
load(paste(“Ginv.ref.s”, i, “.bin”, sep = “”))
for (j in 1: ncol(cv.coor)) {
pheno4cv <- ref.pheno[-na.omit(cv.coor[, j]), ]
gblup <- asreml(sur.status ~ Ori.Location + b.weight,
random = ~vm(IID, Ginv),
dense = ~vm(IID, Ginv),
family = asr_binomial(“probit”),
data = pheno4cv,
maxiter = 100)
y4val <- ref.pheno$sur.status[na.omit(cv.coor[, j])]
gebv4val <- coef(gblup)$random[na.omit(cv.coor[, j])]
cv.res[j, i] <- roc(y4val, gebv4val)$auc
}
}
cv.res [nrow(cv.res), ] <- colMeans(cv.res, na.rm = T)
write.table(cv.res, “Ref.cv.results.gblup.csv”, sep = “,”, row.names = F, quote = F)
実行後、ファイル「Ｒｅｆ．ｃｖ．ｒｅｓｕｌｔｓ．ｇｂｌｕｐ．ｃｓｖ」の最終行情報が、即ちSNP数の変化に応じてＧ－ＢＬＵＰが抗病性ＧＥＢＶを予測する精度である。
【0029】
３、ＢａｙｅｓＣπの予測精度を評価する
Ｇ－ＢＬＵＰと同じモデルを使用し、ＲでＢａｙｅｓＣπの予測精度を評価する。Ｂａｙｅｓ法は時間がかかるため、ＳＮＰサブセットごとに1つのコマンドを指定する。以下では、ＳＮＰサブセット「ｓ１」を例とし、ＢａｙｅｓＣπの予測精度を推定する方法を説明する。
Ｌｉｎｕｘ環境で次のコマンドを実行する：
#Ｒで次のスクリプトを実行し、ＳＮＰサブセット「ｓ１」を使用する時、ＢａｙｅｓＣπでの予測精度を評価する。
library(data.table)
library(BGLR)
library(pROC)
set.seed(123)
i <- 1 #別のＳＮＰサブセットを使用するためにはｉの値を変更する
load(paste(“ref.snp.s”, i, “.bin”, sep = “”))
ref.pheno <- fread (“ref.pheno4cv.csv”)
cv.coor <- as.matrix(fread(“coor4cv.ref.csv”))
No.CV <- ncol(cv.coor)
gebv <- matrix(nrow = nrow(ref.pheno), ncol = No.CV)
cv.res <- matrix(nrow = No.CV + 1)
colnames(cv.res) <- paste(“s”, i, sep = “”)
rownames(cv.res) <- c(paste(“CV”, 1: ncol(dr.coor), sep = “”), “Average”)
for (i in 1: No.CV) {
pheno4trn <- ref.pheno[-na.omit(cv.coor[, i]), ]
geno4trn <- M[-na.omit(cv.coor[, i]), ]
ETA4trn <- list(FIX = list(~factor(Ori.Location) + b.weight, data = pheno4trn, model = 'FIXED'), MAR = list(X = geno4trn, model = 'BayesC'))
bc <- BGLR(y = pheno4trn$sur.status,
ETA = ETA4trn,
response_type = “ordinal”,
nIter = 20000,
burnIn = 5000,
thin = 10,
saveAt = “bc_”,
verbose = F)
MarEff <- bc$ETA$MAR$b
gebv[, i] <- M %*% MarEff
y4val <- ref.pheno$sur.status[na.omit(cv.coor[, i])]
gebv4val <- gebv[na.omit(cv.coor[, i]), i]
cv.res[i, ] <- roc(y4val, gebv4val)$auc
}
cv.res[nrow(cv.res), ] <- colMeans(cv.res, na.rm = T)
write.table(cv.res, paste(“Ref.results.cv.BC.s”, i, “.csv”, sep = “”), sep = “,”, row.names = F, quote = F)
実行後、ファイル「Ｒｅｆ．ｒｅｓｕｌｔｓ．ｃｖ．ＢＣ．ｓ１．ｃｓｖ」の最終行情報は、即ちＳＮＰサブセット「ｓ１」を使用する時、ＢａｙｅｓＣπで予測した抗病性の精度である。上述のスクリプトを１０回実行し（ｉの範囲は１～１０である）、SNP数の変化に応じてＢａｙｅｓＣπが東星斑の抗ビブリオ・ハーヴェイＧＥＢＶを予測する精度を取得する（図２）。
Ｇ－ＢＬＵＰと同じモデルを使用し、ＲでＢａｙｅｓＣπの予測精度を評価する。Ｂａｙｅｓ法は時間がかかるため、ＳＮＰサブセットごとに1つのコマンドを指定する。以下では、ＳＮＰサブセット「ｓ１」を例とし、ＢａｙｅｓＣπの予測精度を推定する方法を説明する。
Ｌｉｎｕｘ環境で次のコマンドを実行する：
#Ｒで次のスクリプトを実行し、ＳＮＰサブセット「ｓ１」を使用する時、ＢａｙｅｓＣπでの予測精度を評価する。
library(data.table)
library(BGLR)
library(pROC)
set.seed(123)
i <- 1 #別のＳＮＰサブセットを使用するためにはｉの値を変更する
load(paste(“ref.snp.s”, i, “.bin”, sep = “”))
ref.pheno <- fread (“ref.pheno4cv.csv”)
cv.coor <- as.matrix(fread(“coor4cv.ref.csv”))
No.CV <- ncol(cv.coor)
gebv <- matrix(nrow = nrow(ref.pheno), ncol = No.CV)
cv.res <- matrix(nrow = No.CV + 1)
colnames(cv.res) <- paste(“s”, i, sep = “”)
rownames(cv.res) <- c(paste(“CV”, 1: ncol(dr.coor), sep = “”), “Average”)
for (i in 1: No.CV) {
pheno4trn <- ref.pheno[-na.omit(cv.coor[, i]), ]
geno4trn <- M[-na.omit(cv.coor[, i]), ]
ETA4trn <- list(FIX = list(~factor(Ori.Location) + b.weight, data = pheno4trn, model = 'FIXED'), MAR = list(X = geno4trn, model = 'BayesC'))
bc <- BGLR(y = pheno4trn$sur.status,
ETA = ETA4trn,
response_type = “ordinal”,
nIter = 20000,
burnIn = 5000,
thin = 10,
saveAt = “bc_”,
verbose = F)
MarEff <- bc$ETA$MAR$b
gebv[, i] <- M %*% MarEff
y4val <- ref.pheno$sur.status[na.omit(cv.coor[, i])]
gebv4val <- gebv[na.omit(cv.coor[, i]), i]
cv.res[i, ] <- roc(y4val, gebv4val)$auc
}
cv.res[nrow(cv.res), ] <- colMeans(cv.res, na.rm = T)
write.table(cv.res, paste(“Ref.results.cv.BC.s”, i, “.csv”, sep = “”), sep = “,”, row.names = F, quote = F)
実行後、ファイル「Ｒｅｆ．ｒｅｓｕｌｔｓ．ｃｖ．ＢＣ．ｓ１．ｃｓｖ」の最終行情報は、即ちＳＮＰサブセット「ｓ１」を使用する時、ＢａｙｅｓＣπで予測した抗病性の精度である。上述のスクリプトを１０回実行し（ｉの範囲は１～１０である）、SNP数の変化に応じてＢａｙｅｓＣπが東星斑の抗ビブリオ・ハーヴェイＧＥＢＶを予測する精度を取得する（図２）。
【0030】
４、異なるＳＮＰ数における２つのゲノム選択方法の予測精度を比較する
SNP数の変化に応じてＧ－ＢＬＵＰとＢａｙｅｓＣπが東星斑の抗ビブリオ・ハーヴェイＧＥＢＶを予測する精度の変化傾向を図２に示す。その中からわかるように、ＳＮＰ数の減少に応じて、Ｇ－ＢＬＵＰとＢａｙｅｓＣπによる予測精度は低下する傾向を示し、ＳＮＰ数が５０ｋの場合、Ｇ－ＢＬＵＰとＢａｙｅｓＣπによる予測精度が最も高く、それぞれ０．６７０と０．６８１であって、ＳＮＰマーカー数が２０ｋから３ｋに減少すると、Ｇ－ＢＬＵＰとＢａｙｅｓＣπによる予測精度は許容範囲内で減少し、マーカー数が２ｋ未満の場合、２つの方法による予測精度は０．６０６と０．６１１までに急速に減少する。要するに、ＢａｙｅｓＣπの予測精度はＧ－ＢＬＵＰよりわずかに優れ、ＢａｙｅｓＣπを抗病性ＧＥＢＶを推定する方法として利用することができる。ＧＥＢＶの推定に使用されるマーカー数が多いほど、遺伝子型のコストも高くなるため、本発明は、２０ｋのＳＮＰ部位を使用して抗病性ＧＥＢＶを推定することにより、ゲノム選択方法の予測精度を保証するだけでなく、遺伝子型のコストを適切に削減し、計算にかかる時間を節約することができる。
SNP数の変化に応じてＧ－ＢＬＵＰとＢａｙｅｓＣπが東星斑の抗ビブリオ・ハーヴェイＧＥＢＶを予測する精度の変化傾向を図２に示す。その中からわかるように、ＳＮＰ数の減少に応じて、Ｇ－ＢＬＵＰとＢａｙｅｓＣπによる予測精度は低下する傾向を示し、ＳＮＰ数が５０ｋの場合、Ｇ－ＢＬＵＰとＢａｙｅｓＣπによる予測精度が最も高く、それぞれ０．６７０と０．６８１であって、ＳＮＰマーカー数が２０ｋから３ｋに減少すると、Ｇ－ＢＬＵＰとＢａｙｅｓＣπによる予測精度は許容範囲内で減少し、マーカー数が２ｋ未満の場合、２つの方法による予測精度は０．６０６と０．６１１までに急速に減少する。要するに、ＢａｙｅｓＣπの予測精度はＧ－ＢＬＵＰよりわずかに優れ、ＢａｙｅｓＣπを抗病性ＧＥＢＶを推定する方法として利用することができる。ＧＥＢＶの推定に使用されるマーカー数が多いほど、遺伝子型のコストも高くなるため、本発明は、２０ｋのＳＮＰ部位を使用して抗病性ＧＥＢＶを推定することにより、ゲノム選択方法の予測精度を保証するだけでなく、遺伝子型のコストを適切に削減し、計算にかかる時間を節約することができる。
【実施例4】
【0031】
実施例４：東星斑候補グループＧＥＢＶを予測するためのＳＮＰマーカーを選び出す
実施例３の結果に基づいて、本発明は、候補グループの抗病性ＧＥＢＶを推定するために２０，０００個（２０ｋ）のＳＮＰ部位を選択するが、その中には、東星斑の抗ビブリオ・ハーヴェイに関連する３０部位と東星斑ゲノム上に均等に分布する１９，９７０個の部位とが含まれる。抗病性に関連する部位は、ＳＥＤ ＩＤ ＮＯ：１～ＳＥＤ ＩＤ ＮＯ：１５のいずれかの配列の第３６位と第１０７位に位置し、ゲノムワイド関連解析により取得される。他のＳＮＰ部位は東星斑ゲノム上に均等に分布され、部位はＳＥＤ ＩＤ ＮＯ：１６～ＳＥＤ ＩＤ ＮＯ：１００００のいずれかの配列の第３６位と第１０７位に位置する。
実施例３の結果に基づいて、本発明は、候補グループの抗病性ＧＥＢＶを推定するために２０，０００個（２０ｋ）のＳＮＰ部位を選択するが、その中には、東星斑の抗ビブリオ・ハーヴェイに関連する３０部位と東星斑ゲノム上に均等に分布する１９，９７０個の部位とが含まれる。抗病性に関連する部位は、ＳＥＤ ＩＤ ＮＯ：１～ＳＥＤ ＩＤ ＮＯ：１５のいずれかの配列の第３６位と第１０７位に位置し、ゲノムワイド関連解析により取得される。他のＳＮＰ部位は東星斑ゲノム上に均等に分布され、部位はＳＥＤ ＩＤ ＮＯ：１６～ＳＥＤ ＩＤ ＮＯ：１００００のいずれかの配列の第３６位と第１０７位に位置する。
【0032】
Ｌｉｎｕｘ環境で次のコマンドを実行する：
#ＰＬＩＮＫを使用して遺伝子型データを準備する
plink2 --bpfile snps.imputed --autosome-num 24 --keep indiv.ref.txt --make-bpgen --out geno4gwas
#ＧＣＴＡを使用してＧＲＭ行列を生成する
gcta --bpfile geno4gwas --autosome-num 24 --make-grm --out grm.gcta
#ＧＣＴＡを使用してＧＲＭ疎行列を生成する
gcta --grm grm.gcta --make-bK-sparse 0.05 --out sparse.grm.gcta
#抗病性表現型と共変量をそれぞれテキストファイル「ｐｈｅｎｏ４ｇｃｔａ．ｐｈｅ」と「ｑｃｏｖ４ｇｃｔａ．ｑｃｏｖ」へ整理する。共変量には体重とＰＣＡ分析の上位１０個の主成分ベクトルが含まれる
#ＧＣＴＡを使用してゲノムワイド関連解析を実行する
gcta --bpfile geno4gwas --grm-sparse sparse.grm.gcta --pheno pheno4gcta.phe-qcovar qcov4gcta.qcov --fastGWA-mlm --autosome-num 24 --out res.gwas
#Ｒを使用して解析結果を読み取り、抗病性性状に関連する３０個のＳＮＰ情報を出力する
library(data.table)
res <- fread(“res.gwas.fastGWA”)
setorder(res, P)
sel.loci <- res[1: 30, ]
setorder(sel.loci, CHR, POS)
write.table(sel.loci, “selected.loci.info.txt”, sep = “＼t”, quote = F)
実行後、抗病性に関連する３０個のＳＮＰ部位の番号、染色体番号、物理的位置などの情報をファイル「ｓｅｌｅｃｔｅｄ．ｌｏｃｉ．ｉｎｆｏ．ｔｘｔ」に保存する。その後、ＧＥＢＶの計算に使用される２０ｋのＳＮＰを選択し、これらの２０ｋのＳＮＰには上述３０個のＳＮＰが含まれる。
#ＰＬＩＮＫを使用して遺伝子型データを準備する
plink2 --bpfile snps.imputed --autosome-num 24 --keep indiv.ref.txt --make-bpgen --out geno4gwas
#ＧＣＴＡを使用してＧＲＭ行列を生成する
gcta --bpfile geno4gwas --autosome-num 24 --make-grm --out grm.gcta
#ＧＣＴＡを使用してＧＲＭ疎行列を生成する
gcta --grm grm.gcta --make-bK-sparse 0.05 --out sparse.grm.gcta
#抗病性表現型と共変量をそれぞれテキストファイル「ｐｈｅｎｏ４ｇｃｔａ．ｐｈｅ」と「ｑｃｏｖ４ｇｃｔａ．ｑｃｏｖ」へ整理する。共変量には体重とＰＣＡ分析の上位１０個の主成分ベクトルが含まれる
#ＧＣＴＡを使用してゲノムワイド関連解析を実行する
gcta --bpfile geno4gwas --grm-sparse sparse.grm.gcta --pheno pheno4gcta.phe-qcovar qcov4gcta.qcov --fastGWA-mlm --autosome-num 24 --out res.gwas
#Ｒを使用して解析結果を読み取り、抗病性性状に関連する３０個のＳＮＰ情報を出力する
library(data.table)
res <- fread(“res.gwas.fastGWA”)
setorder(res, P)
sel.loci <- res[1: 30, ]
setorder(sel.loci, CHR, POS)
write.table(sel.loci, “selected.loci.info.txt”, sep = “＼t”, quote = F)
実行後、抗病性に関連する３０個のＳＮＰ部位の番号、染色体番号、物理的位置などの情報をファイル「ｓｅｌｅｃｔｅｄ．ｌｏｃｉ．ｉｎｆｏ．ｔｘｔ」に保存する。その後、ＧＥＢＶの計算に使用される２０ｋのＳＮＰを選択し、これらの２０ｋのＳＮＰには上述３０個のＳＮＰが含まれる。
【0033】
Ｌｉｎｕｘ環境で次のコマンドを実行する：
#ゲノムを均等にカバーする２０ｋのＳＮＰを選び出す
##ＰＬＩＮＫ２を使用してゲノムを均等にカバーする１９，９７０個のＳＮＰを抽出し、－－ｂｐ－ｓｐａｃｅの値を増やすと取得されるＳＮＰ数が少なくなり、逆に－－ｂｐ－ｓｐａｃｅの値を減らすと取得されるＳＮＰ数が多くなる。
plink2 --bpfile snps.imputed --extract comm.loci.txt --bp-space 37267 --make-bpgen --out snps.20k.part1
##Ｒを使用し、最後２０ｋのＳＮＰを抽出するためのファイルを整理する
library(data.table)
bim.p1 <- fread(“snps.20k.part1.bim”)
loci.info <- fread(“trait.realted.loci.txt”)
loci.20k <- c(bim.p1$V2, loci.info$SNP)
write.table(loci.20k, “snps.20k.txt”, sep = “＼t”, col.names = F, row.names = F, quote = F)
実行後、ファイル「ｓｎｐｓ．２０ｋ．ｔｘｔ」には２０ｋのＳＮＰ番号が含まれており、これらの番号はＰＬＩＮＫ２に読み取られ、対応する部位を抽出することができる。ファイルにはすべて２０，０００個のＳＮＰ部位が含まれており、その中には東星斑の抗病性に関連する３０個のＳＮＰ部位と、東星斑の参照ゲノムを均一にカバーする１９，９７０個のＳＮＰ部位が含まれる。図３は、１Ｍ範囲ごとで東星斑の参照ゲノム上のこれら２０ｋのＳＮＰの分布を示しており、これらの部位は基本的に東星斑のゲノムを均等にカバーしていることがわかる。
#ゲノムを均等にカバーする２０ｋのＳＮＰを選び出す
##ＰＬＩＮＫ２を使用してゲノムを均等にカバーする１９，９７０個のＳＮＰを抽出し、－－ｂｐ－ｓｐａｃｅの値を増やすと取得されるＳＮＰ数が少なくなり、逆に－－ｂｐ－ｓｐａｃｅの値を減らすと取得されるＳＮＰ数が多くなる。
plink2 --bpfile snps.imputed --extract comm.loci.txt --bp-space 37267 --make-bpgen --out snps.20k.part1
##Ｒを使用し、最後２０ｋのＳＮＰを抽出するためのファイルを整理する
library(data.table)
bim.p1 <- fread(“snps.20k.part1.bim”)
loci.info <- fread(“trait.realted.loci.txt”)
loci.20k <- c(bim.p1$V2, loci.info$SNP)
write.table(loci.20k, “snps.20k.txt”, sep = “＼t”, col.names = F, row.names = F, quote = F)
実行後、ファイル「ｓｎｐｓ．２０ｋ．ｔｘｔ」には２０ｋのＳＮＰ番号が含まれており、これらの番号はＰＬＩＮＫ２に読み取られ、対応する部位を抽出することができる。ファイルにはすべて２０，０００個のＳＮＰ部位が含まれており、その中には東星斑の抗病性に関連する３０個のＳＮＰ部位と、東星斑の参照ゲノムを均一にカバーする１９，９７０個のＳＮＰ部位が含まれる。図３は、１Ｍ範囲ごとで東星斑の参照ゲノム上のこれら２０ｋのＳＮＰの分布を示しており、これらの部位は基本的に東星斑のゲノムを均等にカバーしていることがわかる。
【実施例5】
【0034】
実施例５：ＢａｙｅｓＣπと好ましい２０ｋのＳＮＰマーカーの予測効果を検証する
１、交差検証を利用して予測効果を評価する
実施例３で記載されたのと同じモデルを使用し、実施例４で取得された２０ｋのＳＮＰを検証する。実施例３の結果と比較すると、次のことがわかる：東星斑点ゲノムの均一に分布された２０ｋのＳＮＰ部位（実施例４で選択した３０個の抗病性に関連するＳＮＰ部位を含まない）を使用したＢａｙｅｓＣπの予測精度は０．６７０であり（図２）、これらの３０個の抗病性に関連する部位を含む２０ｋのＳＮＰ部位を使用した後、ＢａｙｅｓＣπの予測精度は０．６８２であり、０．０１２が増加され、実施例３で５０ｋのＳＮＰを使用したＢａｅｙｓＣπの予測精度（０．６８１）に近い。従って、本発明により提供される２０ｋのＳＮＰを使用すると、ゲノム選択を無事完了できるだけでなく、良好な予測効果も有する。
１、交差検証を利用して予測効果を評価する
実施例３で記載されたのと同じモデルを使用し、実施例４で取得された２０ｋのＳＮＰを検証する。実施例３の結果と比較すると、次のことがわかる：東星斑点ゲノムの均一に分布された２０ｋのＳＮＰ部位（実施例４で選択した３０個の抗病性に関連するＳＮＰ部位を含まない）を使用したＢａｙｅｓＣπの予測精度は０．６７０であり（図２）、これらの３０個の抗病性に関連する部位を含む２０ｋのＳＮＰ部位を使用した後、ＢａｙｅｓＣπの予測精度は０．６８２であり、０．０１２が増加され、実施例３で５０ｋのＳＮＰを使用したＢａｅｙｓＣπの予測精度（０．６８１）に近い。従って、本発明により提供される２０ｋのＳＮＰを使用すると、ゲノム選択を無事完了できるだけでなく、良好な予測効果も有する。
【0035】
２、検証グループの抗病性ＧＥＢＶを利用して予測効果を評価する
抗病性表現型を有するが参照グループには含まれていない個体を検証個体とし、ＢａｙｅｓＣπと実施例４で取得された２０ｋのＳＮＰを使用して検証グループの抗病性ＧＥＢＶを算出し、生存個体と死亡個体の平均ＧＥＢＶを比較し、本発明の上記方法の実現可能性を検証する。Ｌｉｎｕｘ環境で次のコマンドを実行する：
#候補となる個体の遺伝子型を抽出する
##ＰＬＩＮＫ２で次のコマンドを実行する
plink2 --bpfile snps.imputed --autosome-num 24 --extract snps.20k.txt --keep indiv.ref.txt --export-allele export.A.txt --export A --out geno.Ref.20k
plink2 --bpfile snps.imputed --autosome-num 24 --extract snps.20k.txt --keep indiv.val.txt --export-allele export.A.txt --export A --out geno.Val.20k
awk ‘{print $2}’ geno.Val.20k.raw | sed ‘1d’ > iid.val.20k.txt
#検証グループの抗病性ＧＥＢＶを計算する
library(data.table)
library(BGLR)
set.seed(123)
pheno.ref <- fread(“pheno.Ref.csv”)
pheno.val <- fread(“pheno.Val.csv”)
iid.val <- fread(“iid.val.20k.txt”, header = F)$V1
M.ref <- fread(“geno.Ref.20k.raw”)
M.val <- fread(“geno.Val.20k.raw”)
M.ref <- as.matrix(M.ref[, c(1: 6):= NULL])
M.val <- as.matrix(M.val[, c(1: 6):= NULL])
ETA4MarEst <- list(FIX = list(~factor(Ori.Location) + b.weight, data = pheno.ref, model = ‘FIXED’), MAR = list(X = M.ref, model = ‘BayesC’))
BC <- BGLR(y = pheno$sur.status, ETA = ETA4MarEst,
response_type = “ordinal”,
nIter = 20000, burnIn = 5000, thin = 10,
saveAt = “BC4Val_”, verbose = T)
MarEff <- BC$ETA$MAR$b
gebv <- data.frame(IID = iid.val, GEBV = M.val %*% MarEff)
gebv <- merge(gebv, pheno.val, sort = F)
mean(gebv$GEBV[gebv$sur.status == 0])
mean(gebv$GEBV[gebv$sur.status == 1])
length(gebv$GEBV[gebv$sur.status == 0])
length(gebv$GEBV[gebv$sur.status == 1])
実行後、２９８匹の検証グループの抗病性ＧＥＢＶを予測することに成功し、生存個体と死亡個体の平均ＧＥＢＶを図４に示した。その中で、生存個体は１１０匹で、平均ＧＥＢＶは０．４１３で、全体の平均GEBVより高く、死亡個体は１８８匹で、平均ＧＥＢＶは０．２５５で、全体の平均ＧＥＢＶよりも低かった。このことからわかるように、ＧＥＢＶは個体の抗病性と正の相関があり、ＧＥＢＶが高い個体は抗病性が強いため、ＧＥＢＶを使用して個体を選択し、その中から抗病性の強い個体の選び出して繁殖に利用できる。
抗病性表現型を有するが参照グループには含まれていない個体を検証個体とし、ＢａｙｅｓＣπと実施例４で取得された２０ｋのＳＮＰを使用して検証グループの抗病性ＧＥＢＶを算出し、生存個体と死亡個体の平均ＧＥＢＶを比較し、本発明の上記方法の実現可能性を検証する。Ｌｉｎｕｘ環境で次のコマンドを実行する：
#候補となる個体の遺伝子型を抽出する
##ＰＬＩＮＫ２で次のコマンドを実行する
plink2 --bpfile snps.imputed --autosome-num 24 --extract snps.20k.txt --keep indiv.ref.txt --export-allele export.A.txt --export A --out geno.Ref.20k
plink2 --bpfile snps.imputed --autosome-num 24 --extract snps.20k.txt --keep indiv.val.txt --export-allele export.A.txt --export A --out geno.Val.20k
awk ‘{print $2}’ geno.Val.20k.raw | sed ‘1d’ > iid.val.20k.txt
#検証グループの抗病性ＧＥＢＶを計算する
library(data.table)
library(BGLR)
set.seed(123)
pheno.ref <- fread(“pheno.Ref.csv”)
pheno.val <- fread(“pheno.Val.csv”)
iid.val <- fread(“iid.val.20k.txt”, header = F)$V1
M.ref <- fread(“geno.Ref.20k.raw”)
M.val <- fread(“geno.Val.20k.raw”)
M.ref <- as.matrix(M.ref[, c(1: 6):= NULL])
M.val <- as.matrix(M.val[, c(1: 6):= NULL])
ETA4MarEst <- list(FIX = list(~factor(Ori.Location) + b.weight, data = pheno.ref, model = ‘FIXED’), MAR = list(X = M.ref, model = ‘BayesC’))
BC <- BGLR(y = pheno$sur.status, ETA = ETA4MarEst,
response_type = “ordinal”,
nIter = 20000, burnIn = 5000, thin = 10,
saveAt = “BC4Val_”, verbose = T)
MarEff <- BC$ETA$MAR$b
gebv <- data.frame(IID = iid.val, GEBV = M.val %*% MarEff)
gebv <- merge(gebv, pheno.val, sort = F)
mean(gebv$GEBV[gebv$sur.status == 0])
mean(gebv$GEBV[gebv$sur.status == 1])
length(gebv$GEBV[gebv$sur.status == 0])
length(gebv$GEBV[gebv$sur.status == 1])
実行後、２９８匹の検証グループの抗病性ＧＥＢＶを予測することに成功し、生存個体と死亡個体の平均ＧＥＢＶを図４に示した。その中で、生存個体は１１０匹で、平均ＧＥＢＶは０．４１３で、全体の平均GEBVより高く、死亡個体は１８８匹で、平均ＧＥＢＶは０．２５５で、全体の平均ＧＥＢＶよりも低かった。このことからわかるように、ＧＥＢＶは個体の抗病性と正の相関があり、ＧＥＢＶが高い個体は抗病性が強いため、ＧＥＢＶを使用して個体を選択し、その中から抗病性の強い個体の選び出して繁殖に利用できる。
【実施例6】
【0036】
実施例６：東星斑候補グループの抗病性ＧＥＢＶを計算する
１００匹の東星斑の幼魚を選択して候補グループとし、全ゲノムリシーケンシングを使用して候補グループ遺伝子型を取得し、遺伝子型情報をファイル「ｃａｎｄｉｄａｔｅｓ．ｇｅｎｏ．ｖｃｆ．ｇｚ」に保存する。シーケンシングが完了した後、実施例４に記載されたのと同じ２０ｋのＳＮＰを抽出し、実施例３の参照グループと組み合わせ、ＢａｙｅｓＣπを利用して候補グループの抗病性ＧＥＢＶを計算する。高い抗病性ＧＥＢＶを有する個体の抗病性が強いことをさらに説明するために、本発明は感染実験を利用して候補グループの抗病性を測定し、生存個体を抗病性が高い個体として認定し、死亡個体を感染され易い個体として認定した。これら１００匹の東星斑の候補グループには、それぞれ５０匹の抗病性が高い個体と５０匹の感染され易い個体が含まれる。Ｌｉｎｕｘ環境で次のコマンドを実行する：
#候補グループ遺伝子型ファイルから実施例４で記載された２０ｋのＳＮＰを抽出する
##Ｒで次のスクリプトを実行する
library(data.table)
snps.20k <- fread(“snps.20k.txt”)
pos <- as.integer(matrix(unlist(strsplit(snps.20k$V1, ‘:’)), ncol = 2, byrow = T)[, 2])
chr <- gsub(“rs”, “chr”, matrix(unlist(strsplit(snps.20k$V1, ‘:’)), ncol = 2, byrow = T)[, 1]
info <- data.table(CHR = chr, POS = pos)
fwrite(info, “info.20k.txt”, sep = “＼t”, col.names = F, quote = F)
##ＶＣＦｔｏｏｌｓで次のコマンドを実行する
vcftools --gzvcf candidates.geno.vcf.gz --positions info.20k.txt --recode -- recode-INFO-all --stdout | bgzip -c -f -@ 10 > candidates.geno.20k.vcf.gz
##ＰＬＩＮＫ２で次のコマンドを実行する
plink2 --vcf candidates.geno.20k.vcf.gz --autosome-num 24 --make-bpgen --out Can.geno.20k
##Ｒで次のスクリプトを実行する
library(data.table)
bim <- fread(“Can.geno.20k.bim”)
bim$V2 <- paste(“rs”, bim$V1, ‘:’, bim$V4, sep = “”)
fwrite(bim, “Can.geno.20k.bim”, sep = “＼t”, col.names = F, quote = F)
##ＰＬＩＮＫ２で次のコマンドを実行する
plink2 --bpfile Can.geno.20k --autosome-num 24 --export-allele export.A.txt --export A --out geno.Can.20k
awk ‘{print $2}’ geno.Val.20k.raw | sed ‘1d’ > iid.can.20k.txt
#候補グループＧＥＢＶを計算する。最初の５０個の個体は感染され易い個体であり、最後の５０個の個体は抗病性が高い個体である。
##Ｒで次のスクリプトを実行する
library(data.table)
library(BGLR)
setseed(123)
pheno.ref <- fread(“pheno.Ref.csv”)
iid.can <- fread(“iid.can.20k.txt”, header = F)$V1
M.ref <- fread(“geno.Ref.20k.raw”)
M.can <- fread(“geno.Can.20k.raw”)
M.ref <- as.matrix(M.ref[, c(1: 6):= NULL])
M.can <- as.matrix(M.can[, c(1: 6):= NULL])
ETA4MarEst <- list(FIX = list(~factor(Ori.Location) + b.weight, data = pheno.ref, model = ‘FIXED’), MAR = list(X = M.ref, model = ‘BayesC’))
BC <- BGLR(y = pheno$sur.status, ETA = ETA4MarEst,
response_type = “ordinal”,
nIter = 20000, burnIn = 5000, thin = 10,
saveAt = “BC_”, verbose = F)
MarEff <- BC$ETA$MAR$b
gebv <- data.frame(IID = iid.can, GEBV = M.val %*% MarEff, pheno = c(rep(0, 50), rep(1, 50)))
fwrite(gebv, “GEBV.candidates.csv”, sep = “,”, quote = F)
実行後、１００匹の東星斑候補グループの抗病性ＧＥＢＶの計算に成功し、５０匹の抗病性が高い個体のＧＥＢＶ平均値は０．６４３で、感染され易い個体のＧＥＢＶ平均値（０．５２５）よりも高かった（図５）。また、ＧＥＢＶの上位２０％個体の中には、抗病性が高い個体が１５匹で、感染され易い個体（５匹）よりもはるかに多かった。これは、ＧＥＢＶが高い個体は抗病性が強いため、子孫を繁殖させる親として使用し、子孫の耐病性を向上させることができることをさらに示している。同時に、本発明の実施例３で選択されたＢａｙｅｓＣπ法と実施例４で好ましく選択された２０ｋのＳＮＰが良好な予測効果を有し、抗病性の強い個体のスクリーニングに使用することができることも示している。
１００匹の東星斑の幼魚を選択して候補グループとし、全ゲノムリシーケンシングを使用して候補グループ遺伝子型を取得し、遺伝子型情報をファイル「ｃａｎｄｉｄａｔｅｓ．ｇｅｎｏ．ｖｃｆ．ｇｚ」に保存する。シーケンシングが完了した後、実施例４に記載されたのと同じ２０ｋのＳＮＰを抽出し、実施例３の参照グループと組み合わせ、ＢａｙｅｓＣπを利用して候補グループの抗病性ＧＥＢＶを計算する。高い抗病性ＧＥＢＶを有する個体の抗病性が強いことをさらに説明するために、本発明は感染実験を利用して候補グループの抗病性を測定し、生存個体を抗病性が高い個体として認定し、死亡個体を感染され易い個体として認定した。これら１００匹の東星斑の候補グループには、それぞれ５０匹の抗病性が高い個体と５０匹の感染され易い個体が含まれる。Ｌｉｎｕｘ環境で次のコマンドを実行する：
#候補グループ遺伝子型ファイルから実施例４で記載された２０ｋのＳＮＰを抽出する
##Ｒで次のスクリプトを実行する
library(data.table)
snps.20k <- fread(“snps.20k.txt”)
pos <- as.integer(matrix(unlist(strsplit(snps.20k$V1, ‘:’)), ncol = 2, byrow = T)[, 2])
chr <- gsub(“rs”, “chr”, matrix(unlist(strsplit(snps.20k$V1, ‘:’)), ncol = 2, byrow = T)[, 1]
info <- data.table(CHR = chr, POS = pos)
fwrite(info, “info.20k.txt”, sep = “＼t”, col.names = F, quote = F)
##ＶＣＦｔｏｏｌｓで次のコマンドを実行する
vcftools --gzvcf candidates.geno.vcf.gz --positions info.20k.txt --recode -- recode-INFO-all --stdout | bgzip -c -f -@ 10 > candidates.geno.20k.vcf.gz
##ＰＬＩＮＫ２で次のコマンドを実行する
plink2 --vcf candidates.geno.20k.vcf.gz --autosome-num 24 --make-bpgen --out Can.geno.20k
##Ｒで次のスクリプトを実行する
library(data.table)
bim <- fread(“Can.geno.20k.bim”)
bim$V2 <- paste(“rs”, bim$V1, ‘:’, bim$V4, sep = “”)
fwrite(bim, “Can.geno.20k.bim”, sep = “＼t”, col.names = F, quote = F)
##ＰＬＩＮＫ２で次のコマンドを実行する
plink2 --bpfile Can.geno.20k --autosome-num 24 --export-allele export.A.txt --export A --out geno.Can.20k
awk ‘{print $2}’ geno.Val.20k.raw | sed ‘1d’ > iid.can.20k.txt
#候補グループＧＥＢＶを計算する。最初の５０個の個体は感染され易い個体であり、最後の５０個の個体は抗病性が高い個体である。
##Ｒで次のスクリプトを実行する
library(data.table)
library(BGLR)
setseed(123)
pheno.ref <- fread(“pheno.Ref.csv”)
iid.can <- fread(“iid.can.20k.txt”, header = F)$V1
M.ref <- fread(“geno.Ref.20k.raw”)
M.can <- fread(“geno.Can.20k.raw”)
M.ref <- as.matrix(M.ref[, c(1: 6):= NULL])
M.can <- as.matrix(M.can[, c(1: 6):= NULL])
ETA4MarEst <- list(FIX = list(~factor(Ori.Location) + b.weight, data = pheno.ref, model = ‘FIXED’), MAR = list(X = M.ref, model = ‘BayesC’))
BC <- BGLR(y = pheno$sur.status, ETA = ETA4MarEst,
response_type = “ordinal”,
nIter = 20000, burnIn = 5000, thin = 10,
saveAt = “BC_”, verbose = F)
MarEff <- BC$ETA$MAR$b
gebv <- data.frame(IID = iid.can, GEBV = M.val %*% MarEff, pheno = c(rep(0, 50), rep(1, 50)))
fwrite(gebv, “GEBV.candidates.csv”, sep = “,”, quote = F)
実行後、１００匹の東星斑候補グループの抗病性ＧＥＢＶの計算に成功し、５０匹の抗病性が高い個体のＧＥＢＶ平均値は０．６４３で、感染され易い個体のＧＥＢＶ平均値（０．５２５）よりも高かった（図５）。また、ＧＥＢＶの上位２０％個体の中には、抗病性が高い個体が１５匹で、感染され易い個体（５匹）よりもはるかに多かった。これは、ＧＥＢＶが高い個体は抗病性が強いため、子孫を繁殖させる親として使用し、子孫の耐病性を向上させることができることをさらに示している。同時に、本発明の実施例３で選択されたＢａｙｅｓＣπ法と実施例４で好ましく選択された２０ｋのＳＮＰが良好な予測効果を有し、抗病性の強い個体のスクリーニングに使用することができることも示している。
【図1】
【図2】
【図3】
【図4】
【図5】
【配列表】
【国際調査報告】