Draft genome of leopard gecko, Eublepharis macularius

サンプル収集と配列決定

Eublepharis macularius: NCBI taxonomy ID 481883; specimen ID TG1477,雄の尾組織からゲノムを抽出した (fig. 1)． すべての組織は、University of Minnesota animal use protocols 0810A50001 および 1108A03545 に従って収集された。 この動物は、Woodland Park Zoo (シアトル) のインド由来の動物とNational Zoo (Washington, DC) のパキスタン由来の輸入動物に由来する系統の30数世代にわたる近親交配から飼育下で生まれたものであった。 170から20kbの勾配挿入サイズを持つ合計7つのペアエンドライブラリーを構築し、Illumina HiSeq 2000プラットフォームでメーカーの説明書に従って配列決定した（Illumina, San Diego, California, USA）。 ロングインサートサイズライブラリー（2、5、10、20kb）では、シーケンスのリード長は49bpであり、ショートインサートサイズライブラリー（170、500、800bp）では、シーケンスのリード長は100および150bpであった（表1）。 最終的に合計303Gb（136X）の生配列が得られた（表1）。 アセンブルの前に、SOAPdenovoパッケージのソフトウェアSOAPfilterを用いて、低品質リードやライブラリ構築時のPCR増幅に起因する重複リードの除去を含む厳密な品質管理を行った。 シーケンスエラーはSOAPec (version 2.02) のk-mer frequencyメソッドを用いて修正した。 この結果、187 Gb (84X)の高品質な配列が得られ、ゲノムアセンブリを行うことができた（Table 1）

Leopard gecko Eublepharis maculariusの例 (image from Tony Gamble)

ゲノムアセンブリ

まず、170および500bpインサートサイズライブラリーの54Gbクリーンシーケンスを用いて、ヒョウモントカゲモドキのゲノムサイズを推定するために17-mer分析を実施しました。 簡単に説明すると、リードを17bpのスライド式短鎖に分割し、最初の塩基対を除いて16bpずつオーバーラップさせた。 17-merのカウント分布はポアソン分布に従った（Additional file 1）。 表2 17-mer解析によるゲノムサイズ推定値の統計量。 ゲノムサイズは式にしたがって推定した。 ゲノムサイズ＝# Kmers/Peak of depth

次に、SOAPdenovo（バージョン2.0）を用いて、コンティグ構築、スキャフォールド、ギャップフィリングの3段階で高品質のヒョウモントカゲモドキのアセンブリを行いました。 コンティグ構築では、SOAPdenovoを用いて、ショートインサートライブラリーから得られた高品質リードをペアエンド情報を無視したkmerに分割してde Bruijnグラフを作成し、kmerのマージ、チップクリップ、バブルのマージ、低カバレッジリンクの削除を行った。 次に、de Bruijnグラフで曖昧さのない接続を示すコンティグを収集した。 一連のkmer長をテストし、最も長いN50値を持つコンティグアセンブリを生成するために33-merが選択された。 scaffoldingステップでは、大小両方の挿入ライブラリから得られたリードをコンティグ配列にマッピングし、リードペアからの距離情報を用いてスキャフォールドを構築した。 スキャフォールド内のギャップを埋めるため（ギャップフィリングステップ）、まず170 bp挿入ライブラリから重複するペアエンドリードをCOPEで接続し、次にKgfを用いて接続したリードと他の短い挿入サイズのライブラリからのリードを使用してギャップを閉じました。 さらに、GapCloser を用いて、リードペアの一端がコンティグに一意に整列し、他端がギャップ内に位置するリードの局所アセンブルを行いました。 その結果、総長 2.0 Gb、scaffold および contig の N50 がそれぞれ 664 および 20 kb の leopard gecko ゲノムアセンブリが得られ、これは以前に報告された Gekko japonicus ゲノムアセンブリと同程度であった (Table 3) 。 ヒョウモントカゲモドキゲノムと既報の11種の爬虫類ゲノム（Anolis carolinensis , Python molurus bivittatus , Ophiophagus hannah , Alligator sinensis , Alligator mississippiensis.）のN50を比較したところ、ヒョウモントカゲモドキのN50は664kb、コンティグのN50は20kbであり、既報のゲオミクロンゲノム（Anolis carolinensis）と同等であった。 Gavialis gangeticus, Crocodylus porosus, Chelonia mydas, Pelodiscus sinensis, Pogona vitticeps , and Chrysemys picta bellii ) さらに、我々の結果は同等以上の品質であることが確認されました（表4）。

ゲノム完全性の評価

CEGMAとBUSCO .を使用してアセンブリの完全性を評価しました。 これらは、進化的に予想される遺伝子の含有量を用いてゲノムの完全性を定量的に評価するものです。 CEGMAの評価では、248個の超保存コア真核生物遺伝子のうち225個（91％）が組み込まれ、そのうち210個（85％）が完全であることが示された。 BUSCO解析では、脊椎動物の遺伝子3023個のうち、58個が完全、18％が断片と同定され、24％はアセンブリに欠落しているとされた。 表5 CEGMAで評価したヤモリゲノムにおける真核生物中核遺伝子（CEG）の網羅率。 すべてのCEGは、タンパク質配列の保存の度合いに基づいて4つのグループに分けられた。 グループ1には保存度の低いCEGが含まれ、グループ4には保存度の高いものが含まれる

繰り返しアノテーション

我々は、相似型とde novo手法を組み合わせ、豹ヤモリのゲノムからトランスポゾーム（TE）とその他の繰り返し要素を同定した。 相同性に基づく方法では、Repbase TEライブラリ（RepBase21.01）に対してRepeatMaskerで、TEタンパク質データベースに対してRepeatProteinMaskで検索を行い、既知のTEを同定した。 de novo法では、まずRepeatModeler (http://www.repeatmasker.org/RepeatModeler.html, version 1.0.5) とPilerを用いてde novoヒョウモントカゲモドキリピートライブラリを構築し、その後RepeatMaskerを用いてヒョウモントカゲモドキのリピート配列を注釈付けを行った。 最後に、TRFを用いて、以下のパラメータでタンデムリピートを予測した。 Match = 2, Mismatch = 7, Delta = 7, PM = 80, PI = 10, Minscore = 50.です。 その結果、合計851 Mbの非冗長反復配列が同定され、ヒョウモントカゲモドキゲノムの42 %を占めた。 最も優勢な要素は、長い散在核要素（LINE）であり、全TE配列の30 %、ゲノムの13 %を占めた（表7）。

遺伝子予測

相同性ベース、de novo、トランスクリプトームベース手法を組み合わせて、ヒョウモントカゲモドキのタンパク質コード化遺伝子を予測しました。

相同性ベースの方法では、Ensemblデータベース（release-73）からTaeniopygia guttata、Homo sapiens、Anolis carolinensis、Pelodiscus sinensis、Xenopus tropicalisの遺伝子セットをダウンロードし、その遺伝子セットを用いて、ヒョウモントカゲモドキの遺伝子を予測した。 これらの相同タンパク質配列をヒョウモントカゲモドキゲノムアセンブリに TBLASTN で E-value cutoff 1e-5 でアラインメントし、BLAST ヒットを GenBlastA で候補遺伝子座にリンクさせた。 次に、GeneWiseを用いて候補遺伝子座のゲノム配列とその近傍3kbの配列を抽出し、遺伝子モデルを決定した。 最後に、フレームエラーがあるエクソンを1つだけ持つ偽遺伝子を、これらの遺伝子座はおそらくレトロトランスポジションに由来するため、フィルターにかけた。

de novo法では、相同性ベースの遺伝子セットからオープンリーディングフレーム（ORF）がそのままで、GeneWiseスコアが最も高いヒョウモントカゲモドキ遺伝子1000件をランダムに選び、デフォルトパラメータを用いてAugustus遺伝子予測ツールのトレーニングを実施した。 その後、Augustusを使用して、リピートマスクされたゲノム配列に対してde novo遺伝子予測を行った。 ORFが不完全な遺伝子モデルや、タンパク質コード長<150 bpの小さな遺伝子はフィルタリングで除外された。 最後に、予測された遺伝子をSwissProtデータベースに対してBLASTP検索を行った。 トランスポーズ、トランスポゾン、レトロトランスポゾン、レトロウイルス、レトロトランスポゾン、逆転写酵素、トランスポザーゼ、レトロウイルスのいずれかのキーワードを含むSwissProtタンパク質にマッチする遺伝子がフィルタリングされました。

次に、NCBIデータベースから取得した肝臓、唾液腺、嗅腺、皮膚組織のヒョウモントカゲモドキRNA-seqデータ（アクセッション番号SRR629643、ERR216315、ERR216316、ERR216322、ERR216325、ERR216304、ERR216306）によりトランスリプトームベース遺伝子予測を実施した． Tophat (v1.3.3) を使用して RNA-seq リードをヒョウモントカゲモドキのゲノムアセンブリと比較してスプライスジャンクションを特定し、cufflinks (v2.2.1) を使用してアラインメントした RNA-seq リードを使用して転写物をアセンブルしました。

最後に、相同性、デノボ、トランスクリプトームベースの解析結果を統合し、トランスクリプトームベースの証拠>相同性ベースの証拠>デノボベースの証拠の優先順位に基づいた非冗長参照遺伝子セットとしました。 (1)

GlimmerHMM に含まれる trainGlimmerHMM ツールを用いて、Augustus のトレーニングに使用した 1000 個の高品質遺伝子を用いてマルコフモデルを推定したところ、以下のようになりました。 トランスクリプトームデータから組み立てられた各転写物のコーディングポテンシャルを、マルコフモデルを用いて同定しました。 ORFが完全な転写産物は抽出され、同じ遺伝子座からの複数のアイソフォームは、最長のORFを保持することによって折りたたまれました。 優先度の高い遺伝子モデルが優先度の低いモデルと重複している場合（重複長>100 bp）、後者は削除された。 (3)

Homology-based gene models not supported by transcriptome-based evidence but supported by homologous evidence from least two species were added to the core gene set.転写産物に基づく証拠はないが、少なくとも2つの生物種からの相同性証拠によってサポートされる相同性に基づく遺伝子モデルが、コア遺伝子セットに追加された。

以上の結果、ヒョウモントカゲモドキのゲノムアセンブリにおいて、合計24,755個の非冗長タンパク質コーディング遺伝子をアノテーションした。

Functional annotation of protein-coding genes

BLASTPを用いて機能データベースTrEMBLおよびSwissProtを検索し、全ヒョウモントカゲモドキタンパク質コーディング遺伝子に対して名前を割り当てた（表8）。 次に、ヒョウモントカゲモドキのタンパク質配列をKyoto Encyclopaedia of Genes and Genomes（KEGG）データベースとBLASTPで検索し、その遺伝子が関与している可能性のある分子経路を同定した。 タンパク質のドメインやモチーフは、InterProScan (version 5.16) を用いて、7つの異なるモデル (Profilescan, blastprodom, HmmSmart, HmmPanther, HmmPfam, FPrintScan, PatternScan) でアノテーションされました。 その結果、予測されたヒョウモントカゲモドキのうち20,958個のタンパク質が保存された機能モチーフを持つことが明らかになりました。 また、対応するInterProエントリーから、15,873個のヒョウモントカゲモドキのタンパク質に割り当てられた1028個のGene Ontology (GO) termを得た。

利用可能性と要件

• プロジェクト名：Leopard gecko genome annotation scripts

• プロジェクトのホームページをご覧ください。 https://github.com/gigascience/paper-xiong2016

• オペレーティングシステム。 Linux

• プログラミング言語。 PERL

• その他の要件: なし

• ライセンス: MIT

• 非学会員による使用に対する制限: なし