2024/09/06 追記
ゲノム配列決定は生物学者にとって日常的な作業となったが、遺伝子構造アノテーションの課題は依然として残っており、正確なゲノム・遺伝子研究を妨げている。SynGAPは、遺伝子のシンテニー情報を利用して、ゲノムの遺伝子構造アノテーションを正確かつ自動的に研磨するバイオインフォマティクスツールキットである。SynGAPは、遺伝子構造アノテーションの質の向上と、生物種間の統合的な遺伝子シンテニーのプロファイリングにおいて卓越した機能を提供する。さらに、発現変動インデックスは、系統的に近縁な種で観察される異なる形質の発現に関与する候補遺伝子を探索するための比較トランスクリプトーム解析用に設計されている。
シーケンシング技術と計算技術の進歩は、コストの低下と相まって、研究者が日常的にゲノムの塩基配列を決定し、目的とする高品質のアセンブリを得ることを可能にした。しかし、ゲノムアノテーションは、通常、反復DNA配列のマスキング、遺伝子構造アノテーション(GSA)、遺伝子機能アノテーションという3つの主要なステップを伴うが、生物学者にとっては依然として困難な作業であり、中でも遺伝子構造アノテーションは最も重要かつ困難なステップである。遺伝子構造アノテーションとは、ゲノム配列中の遺伝子の位置を決定し、遺伝子のエクソンとイントロンを正確に定義することである。遺伝子の転写が空間的・時間的に依存していることを考えると、GSAは非常に複雑になる可能性がある。一つの遺伝子が、alternative splicingやalternative start and termination siteによって複数の転写産物に転写される可能性がある。正確なGSAはゲノム・遺伝学研究にとって不可欠である。なぜなら、標準以下のGSAは下流の研究に大きな支障をきたし、誤ったバイオインフォマティクス解析や誤った機能ゲノミクス研究につながるからである。現在、遺伝子構造アノテーションのための様々なパイプラインやワークフローが開発されており、通常、第一原理または相同性ベースの予測やトランスクリプトーム支援アノテーションと統合されている。このようなパイプラインの代表例としては、AUGUSTUS、miniprot、MAKERなどがある。これらのパイプラインのどれもが優れているわけではなく、その結果、異なるゲノムアセンブリ間でGSAの品質に大きなばらつきが生じている。これは、第3世代のロングリードシーケンス技術の利用によるゲノムアセンブリ品質の大幅な向上に追いついていないことに起因する。ApolloやIGV-GSAmanなどのツールキットを用いたGSAの手動補正は、GSAを改善する効果的なアプローチと思われるが、包括的なトランスクリプトームやプロテオームデータに依存し、特にパンゲノムプロジェクトに取り組む場合には時間がかかる。
進化の過程で、染色体上の遺伝子の順序は、共通の祖先種から派生した近縁種で維持される。gene synteny として知られる、異なる種の染色体上の遺伝子の保存された共局在は、染色体の種間進化関係と、ゲノムのシャッフリングイベントの量や位置のような種内ゲノムの変化の両方に対する洞察を提供する。一般に、2つの生物種が近縁であればあるほど、その遺伝子のシンテニーの程度は高くなる。そのため、遺伝子シンテニーはしばしば比較ゲノミクスやトランスクリプトミクス解析に用いられ、相同ゲノムブロックを同定し、種間のオルソログ遺伝子をマッピングする。
シンテニー関係はオルソログ遺伝子の保存的配置を反映するので、異なる生物種間で整列したゲノム領域の遺伝子の比較解析に非常に適している。MCScanX、JCVI、WGDIのようなツールを用いれば、シンテニー解析により、シンテニックゲノム領域内のオルソログ遺伝子ペアや対になっていない遺伝子を容易に同定することができる。これらの対になっていない遺伝子は、遺伝子の欠失や挿入のようなゲノム配列の変化に起因する場合もあれば、単に遺伝子構造のアノテーションが不完全または不正確であることに起因する場合もある。後者の場合、シンテニー解析は近縁種の遺伝子構造アノテーションの相互修正と補完に利用できる。このシナリオに基づき、遺伝子のシンテニーに基づいて近縁種の遺伝子構造アノテーションの欠落を補い、不正確な遺伝子モデルを修正するツールキットSynGAP(Synteny-based Gene Structure Annotation Polisher)を開発した。また、シンテニーに基づくオルソログ遺伝子ペアの同定を改善することで、比較トランスクリプトーム解析へのSynGAPの応用を実証した。
インストール
Github
#conda(bioconda)
mamba install -c conda-forge -c bioconda syngap
#docker(hub)
docker pull yanyew/syngap:1.2.5
docker run -it yanyew/syngap:1.2.5
conda activate syngap
> syngap -h
usage: syngap [-h] {initdb,master,dual,triple,custom,genepair,evi,eviplot} ...
Synteny-based Gene structure Annotation Polisher (SynGAP) https://github.com/yanyew/SynGAP
positional arguments:
{initdb,master,dual,triple,custom,genepair,evi,eviplot}
initdb Import the pre-downloaded compressed file of masterdb for SynGAP master
master Polishment module for one species with masterdb provided by SynGAP
dual Polishment module for two species
triple Polishment module for three species
custom Polishment module for two species with given synteny block
genepair Pair genes from two species
evi Calculate the EVI for each gene pair between two species
eviplot Plot the given EVI file and highlight specific gene pairs
options:
-h, --help show this help message and exit
> syngap dual -h
usage: syngap dual [-h] --sp1fa SP1FA --sp1gff SP1GFF --sp2fa SP2FA --sp2gff SP2GFF [--sp1 SP1] [--sp2 SP2] [--annoType1 ANNOTYPE1] [--annoKey1 ANNOKEY1] [--annoparentKey1 ANNOPARENTKEY1] [--annoType2 ANNOTYPE2] [--annoKey2 ANNOKEY2] [--annoparentKey2 ANNOPARENTKEY2] [--datatype DATATYPE]
[--cscore CSCORE] [--threads THREADS] [--process PROCESS] [--evalue EVALUE] [--rank RANK] [--coverage COVERAGE] [--kmer1 KMER1] [--kmer2 KMER2] [--outs OUTS] [--intron INTRON]
options:
-h, --help show this help message and exit
--sp1fa SP1FA The genome squence file (.fasta format) for species1 [required]
--sp1gff SP1GFF The genome annotation file (.gff format) for species1 [required]
--sp2fa SP2FA The genome squence file (.fasta format) for species2 [required]
--sp2gff SP2GFF The genome annotation file (.gff format) for species2 [required]
--sp1 SP1 The short name for species1, e.g. Ath [default: sp1]
--sp2 SP2 The short name for species2, e.g. Ath [default: sp2]
--annoType1 ANNOTYPE1
Feature type to extract for species1 [default: mRNA]
--annoKey1 ANNOKEY1 Key in the attributes to extract for species1 [default: ID]
--annoparentKey1 ANNOPARENTKEY1
Parent gene key to group with --primary_only in jcvi [default: Parent]
--annoType2 ANNOTYPE2
Feature type to extract for species2 [default: mRNA]
--annoKey2 ANNOKEY2 Key in the attributes to extract for species2 [default: ID]
--annoparentKey2 ANNOPARENTKEY2
Parent gene key to group with --primary_only in jcvi [default: Parent]
--datatype DATATYPE The type of squences for jcvi, nucl|prot [default: nucl]
--cscore CSCORE C-score cutoff for jcvi [default: 0.7]
--threads THREADS Number of threads to use [default: 8]
--process PROCESS Process for gapanno, genblastg|miniprot [default: genblastg]
--evalue EVALUE Threshold for evalue in genBlast [default: 1e-5]
--rank RANK The number of ranks in genBlast output [default: 5]
--coverage COVERAGE Minimum percentage of query gene coverage of the HSP group in the genBlast output [default: 0.5]
--kmer1 KMER1 K-mer size for Indexing in miniprot [default: 5]
--kmer2 KMER2 K-mer size for the second round of chaining in miniprot [default: 4]
--outs OUTS Threshold of Score for miniprot output [default: 0.95]
--intron INTRON Max intron size allowed for miniprot output [default: 40k]
> syngap initdb -h
# syngap initdb -h
usage: syngap initdb [-h] --sp SP --file FILE
options:
-h, --help show this help message and exit
--sp SP The species type of masterdb to be imported, plant|animal [required]
--file FILE The compressed file of masterdb (.tar.gz) to be imported [required]
> syngap master -h
# syngap master -h
usage: syngap master [-h] --sp SP --sp1fa SP1FA --sp1gff SP1GFF [--sp1 SP1] [--annoType1 ANNOTYPE1] [--annoKey1 ANNOKEY1] [--annoparentKey1 ANNOPARENTKEY1] [--datatype DATATYPE] [--cscore CSCORE] [--threads THREADS] [--process PROCESS] [--evalue EVALUE] [--rank RANK] [--coverage COVERAGE]
[--kmer1 KMER1] [--kmer2 KMER2] [--outs OUTS] [--intron INTRON]
options:
-h, --help show this help message and exit
--sp SP The species type of the polished object, plant|animal [required]
--sp1fa SP1FA The genome squence file (.fasta format) for species1 [required]
--sp1gff SP1GFF The genome annotation file (.gff format) for species1 [required]
--sp1 SP1 The short name for species1, e.g. Ath [default: sp1]
--annoType1 ANNOTYPE1
Feature type to extract for species1 [default: mRNA]
--annoKey1 ANNOKEY1 Key in the attributes to extract for species1 [default: ID]
--annoparentKey1 ANNOPARENTKEY1
Parent gene key to group with --primary_only in jcvi [default: Parent]
--datatype DATATYPE The type of squences for jcvi, nucl|prot [default: nucl]
--cscore CSCORE C-score cutoff for jcvi [default: 0.7]
--threads THREADS Number of threads to use [default: 8]
--process PROCESS Process for gapanno, genblastg|miniprot [default: genblastg]
--evalue EVALUE Threshold for evalue in genBlast [default: 1e-5]
--rank RANK The number of ranks in genBlast output [default: 5]
--coverage COVERAGE Minimum percentage of query gene coverage of the HSP group in the genBlast output [default: 0.5]
--kmer1 KMER1 K-mer size for Indexing in miniprot [default: 5]
--kmer2 KMER2 K-mer size for the second round of chaining in miniprot [default: 4]
--outs OUTS Threshold of Score for miniprot output [default: 0.95]
--intron INTRON Max intron size allowed for miniprot output [default: 40k]
> syngap triple -h
# syngap triple -h
usage: syngap triple [-h] --sp1fa SP1FA --sp1gff SP1GFF --sp2fa SP2FA --sp2gff SP2GFF --sp3fa SP3FA --sp3gff SP3GFF [--sp1 SP1] [--sp2 SP2] [--sp3 SP3] [--annoType1 ANNOTYPE1] [--annoKey1 ANNOKEY1] [--annoparentKey1 ANNOPARENTKEY1] [--annoType2 ANNOTYPE2] [--annoKey2 ANNOKEY2]
[--annoparentKey2 ANNOPARENTKEY2] [--annoType3 ANNOTYPE3] [--annoKey3 ANNOKEY3] [--annoparentKey3 ANNOPARENTKEY3] [--datatype DATATYPE] [--cscore CSCORE] [--threads THREADS] [--process PROCESS] [--evalue EVALUE] [--rank RANK] [--coverage COVERAGE] [--kmer1 KMER1]
[--kmer2 KMER2] [--outs OUTS] [--intron INTRON]
options:
-h, --help show this help message and exit
--sp1fa SP1FA The genome squence file (.fasta format) for species1 [required]
--sp1gff SP1GFF The genome annotation file (.gff format) for species1 [required]
--sp2fa SP2FA The genome squence file (.fasta format) for species2 [required]
--sp2gff SP2GFF The genome annotation file (.gff format) for species2 [required]
--sp3fa SP3FA The genome squence file (.fasta format) for species3 [required]
--sp3gff SP3GFF The genome annotation file (.gff format) for species3 [required]
--sp1 SP1 The short name for species1, e.g. Ath [default: sp1]
--sp2 SP2 The short name for species2, e.g. Ath [default: sp2]
--sp3 SP3 The short name for species2, e.g. Ath [default: sp3]
--annoType1 ANNOTYPE1
Feature type to extract for species1 [default: mRNA]
--annoKey1 ANNOKEY1 Key in the attributes to extract for species1 [default: ID]
--annoparentKey1 ANNOPARENTKEY1
Parent gene key to group with --primary_only in jcvi [default: Parent]
--annoType2 ANNOTYPE2
Feature type to extract for species2 [default: mRNA]
--annoKey2 ANNOKEY2 Key in the attributes to extract for species2 [default: ID]
--annoparentKey2 ANNOPARENTKEY2
Parent gene key to group with --primary_only in jcvi [default: Parent]
--annoType3 ANNOTYPE3
Feature type to extract for species3 [default: mRNA]
--annoKey3 ANNOKEY3 Key in the attributes to extract for species3 [default: ID]
--annoparentKey3 ANNOPARENTKEY3
Parent gene key to group with --primary_only in jcvi [default: Parent]
--datatype DATATYPE The type of squences for jcvi, nucl|prot [default: nucl]
--cscore CSCORE C-score cutoff for jcvi [default: 0.7]
--threads THREADS Number of threads to use [default: 8]
--process PROCESS Process for gapanno, genblastg|miniprot [default: genblastg]
--evalue EVALUE Threshold for evalue in genBlast [default: 1e-5]
--rank RANK The number of ranks in genBlast output [default: 5]
--coverage COVERAGE Minimum percentage of query gene coverage of the HSP group in the genBlast output [default: 0.5]
--kmer1 KMER1 K-mer size for Indexing in miniprot [default: 5]
--kmer2 KMER2 K-mer size for the second round of chaining in miniprot [default: 4]
--outs OUTS Threshold of Score for miniprot output [default: 0.95]
--intron INTRON Max intron size allowed for miniprot output [default: 40k]
> syngap genepair -h
# syngap genepair -h
usage: syngap genepair [-h] --sp1fa SP1FA --sp1gff SP1GFF --sp2fa SP2FA --sp2gff SP2GFF [--sp1 SP1] [--sp2 SP2] [--annoType1 ANNOTYPE1] [--annoKey1 ANNOKEY1] [--annoparentKey1 ANNOPARENTKEY1] [--annoType2 ANNOTYPE2] [--annoKey2 ANNOKEY2] [--annoparentKey2 ANNOPARENTKEY2]
[--datatype DATATYPE] [--cscore CSCORE] [--evalue EVALUE] [--iTAK ITAK] [--threads THREADS]
options:
-h, --help show this help message and exit
--sp1fa SP1FA The squence file (.fasta format) for species1 [required]
--sp1gff SP1GFF The annotation file (.gff format) for species1 [required]
--sp2fa SP2FA The squence file (.fasta format) for species2 [required]
--sp2gff SP2GFF The annotation file (.gff format) for species2 [required]
--sp1 SP1 The short name for species1, e.g. Ath [default: sp1]
--sp2 SP2 The short name for species2, e.g. Ath [default: sp2]
--annoType1 ANNOTYPE1
Feature type to extract for species1 [default: mRNA]
--annoKey1 ANNOKEY1 Key in the attributes to extract for species1 [default: ID]
--annoparentKey1 ANNOPARENTKEY1
Parent gene key to group with --primary_only in jcvi [default: Parent]
--annoType2 ANNOTYPE2
Feature type to extract for species2 [default: mRNA]
--annoKey2 ANNOKEY2 Key in the attributes to extract for species2 [default: ID]
--annoparentKey2 ANNOPARENTKEY2
Parent gene key to group with --primary_only in jcvi [default: Parent]
--datatype DATATYPE The type of squences for jcvi, nucl|prot [default: nucl]
--cscore CSCORE C-score cutoff for jcvi [default: 0.7]
--evalue EVALUE Threshold for evalue in two-way blast [default: 1e-2]
--iTAK ITAK Perform iTAK to identify TFs and kinases (only for plants), yse|no [default: no]
--threads THREADS Number of threads to use [default: 8]
実行方法
SynGAP dual - 2つの生物種間の遺伝子構造アノテーションの研磨
2つの生物種間の遺伝子構造アノテーションを相互に補正するために設計されたモジュール。補正対象のゲノム配列とゲノムアノテーションを入力とする。
#docker imageを使うなら
docker run --rm -itv $PWD:/data/ -w /data yanyew/syngap:1.2.5
conda activate syngap
syngap dual \
--sp1fa=Athaliana_167_TAIR9.fa \
--sp1gff=Athaliana_167_TAIR10.gene.gff3 \
--sp2fa=Arabidopsis_halleri.Ahal2.2.dna.toplevel.fa \
--sp2gff=Arabidopsis_halleri.Ahal2.2.52.gff3 \
--sp1=Ath \
--sp2=Aha
-
--sp1fa The genome squence file (.fasta format) for species1 [required]
-
--sp1gff The genome annotation file (.gff format) for species1 [required]
-
--sp2fa The genome squence file (.fasta format) for species2 [required]
-
--sp2gff The genome annotation file (.gff format) for species2 [required]
-
--sp1 The short name for species1, e.g. Ath [default: sp1]
-
--sp2 The short name for species2, e.g. Ath [default: sp2]
- --threads Number of threads to use [default: 8]
出力例
SynGAP.gff3がポリッシュ済みゲノムアノテーションファイル(オリジナル+ポリッシュ済み)となる。
レポジトリによると、ポリッシュ済みのみのアノテーション(*.SynGAP.clean.gff3)、オリジナルのアノテーションででミスアノテーションと判定されたオリジナルのアノテーション(*.SynGAP.clean.miss_annotated.gff3)なども出力されるとある(上の画像にはなし)。
syngap triple - 3種を組み合わせた研磨
syngap triple \
--sp1fa=Athaliana_167_TAIR9.fa \
--sp1gff=Athaliana_167_TAIR10.gene.gff3 \
--sp2fa=Arabidopsis_halleri.Ahal2.2.dna.toplevel.fa \
--sp2gff=Arabidopsis_halleri.Ahal2.2.52.gff3 \
--sp3fa=Brassica_rapa_ro18.SCU_BraROA_2.3.dna.toplevel.fa \
--sp3gff=Brassica_rapa_ro18.SCU_BraROA_2.3.53.chr.gff3 \
--sp1=Ath \
--sp2=Aha \
--sp3=Bra
syngap initdb & syngap master - 既知高品質アノテーションをリファレンスとした研磨
また、1つの生物種の遺伝子構造アノテーションを、著者らがピックアップしたCoreセットで磨くこともできる。Coreセットには、高品質のゲノムアノテーションが施された動植物種が含まれている。
(レポジトリより転載)
まずデータベースplant.tar.gzをダウンロードする。
ダウンロードリンク(animalとplantに分かれている。展開して選択する。ただしセキュアな環境だとダウンロードできない可能性がある)
https://tbtools.cowtransfer.com/s/85ed3920aa7f47
それからインポートする。
syngap initdb --sp=plant --file=plant.tar.gz
- --sp The species type of masterdb to be imported, plant|animal [required]
- --file The compressed file of masterdb (.tar.gz) to be imported [required]
syngap masterを実行する。
syngap master --sp=plant \
--sp1fa=Brassica_rapa_ro18.SCU_BraROA_2.3.dna.toplevel.fa \
--sp1gff=Brassica_rapa_ro18.SCU_BraROA_2.3.53.chr.gff3 \
--sp1=Bra
syngap custom - 特定のシンテニー・ブロックのアノテーション研磨にのみ焦点を当てる場合やjcviではなく他のソフトウェアのシンテニー結果を使用したい場合、そのブロックを含む*.anchorsファイルを提供し、SynGAPカスタムを使用する。
syngap genepair -信頼性の高い異種間相同遺伝子ペアの生成
(レポジトリより)SynGAPには機能モジュールgenepairが組み込まれている。これは、改良されたシンテニー(SynGAP dualまたはtripleから)とベストレシプロカルBLASTを組み合わせることにより、信頼性の高い異種間相同遺伝子ペアを生成する。
syngap genepair \
--sp1fa=Can.fa \
--sp1gff=Can.SynGAP.gff3 \
--sp2fa=Sly.fa \
--sp2gff=Sly.SynGAP.gff3 \
--sp1=Can \
--sp2=Sly
syngap evi - 近縁種の遺伝子差発現解析
(レポジトリより)時系列トランスクリプトームデータの遺伝子発現レベル、発現レベルの差、発現傾向の差に基づいて計算される発現変動指数(EVI)というパラメータを使う。入力発現ファイルはタブ区切りのテキストファイルで、FPKM、RPKM、TPMなどの正規化された発現値(その中でもTPMの使用を推奨する)を入力する。
syngap evi \
--genepair=Can.Sly.final.genepair \
--sp1exp=Can.S1_S7.transcript.TPM.xls \
--sp2exp=Sly.S1_S7.transcript.TPM.xls
遺伝子ペアをEVIでランク付けしたファイルや、EVIの閾値.EVIが閾値を超えた遺伝子ペ(顕著な発現差シグナルを示すと考えられる)などが出力される。また、特定の遺伝子ペアに興味がある場合はsyngap eviplotを使うことができる(レポジトリより)。
論文より
- A. thalianaの遺伝子構造アノテーションを用い、100個の遺伝子アノテーションをランダムに削除したA. thaliana GSAの修正版と、A. thalianaと進化的距離の異なる他の3つの生物種との比較を行った。ランダムに削除された遺伝子アノテーションはSynGAPによりほぼ完全に検索された。検索された遺伝子モデルはオリジナルのものと高いCDS構造の一致率を保ち、高い遺伝子構造精度(GT-AGイントロンスプライスで93%以上、開始コドンで89%以上、停止コドンでほぼ100%)を示した。
- 進化的距離の増加とともにアノテーションの性能は低下するが、SynGAPは依然としてかなりのアノテーションを高い品質で取得できる。
- さまざまな植物ゲノムの組み合わせでSynGAPをテストした。すべてのゲノムはSynGAP dualで磨いた後、高信頼性の遺伝子アノテーションの数が増加した。シンテニー遺伝子ペアの数が増加しただけでなく、BUSCO評価によるGSAの完全性もいくつかの種で大きく改善された(論文図2)。動物ゲノムについてもすべての異なる組み合わせで同様の研磨効果が得られた。
-
比較トランスクリプトーム解析とは、異なる生物種における相同遺伝子間の発現パターンを比較することである。このアプローチの主な課題の1つは、2つの生物種間の相同遺伝子の最適な1対1の関係を確立することである。SynGAPはもう一つの機能モジュールであるgenepairを組み込み、改良されたシンテニー(SynGAPデュアルまたはトリプルから)とベストレシプロカルBLASTを組み合わせることにより、信頼性の高い異種間相同遺伝子ペアを生成する(図4a)。SynGAPによって発見されたペア遺伝子は、各生物種における遺伝子総数の大部分を占め、種内および同属間の遺伝子ペアリングでは65%以上に達したが、進化的距離が長くなるにつれて減少した。~ SynGAP genepairで得られた信頼性の高い遺伝子ペアは、比較トランスクリプトーム解析のリファレンス遺伝子として使用するのに適している。
引用
SynGAP: a synteny-based toolkit for gene structure annotation polishing
Fengqi Wu, Yingxiao Mai, Chengjie Chen & Rui Xia
Genome Biology volume 25, Article number: 218 (2024)