BRAKER 3.0.0をインストールしたので、それぞれの遺伝子予測方法の結果を比較してみた。
Brakerによる遺伝子予測
GeneMark-EXおよびAUGUSTUSのtrainingに使用する配列情報により、3つの方法が提案されている。
テストデータ
- ドジョウ（中国系統）のドラフトゲノム
- 遺伝子予測に用いたRNA-seqデータ（Illumina short reads)
いくつかの異なる組織由来のRNA-seqデータを使用 これらをHisat2でリファレンスゲノムにマッピングし、BAM形式で読み込ませた。
DRR328565, SRR5997795, SRR7586852, SRR7586854,SRR14160170, SRR14160176
- タンパク質のアミノ酸配列データ
遺伝子構造予測のtraining（Braker2, Braker3）に使用するタンパク質のアミノ酸配列は、OrthoDBのVertebrataをダウンロードして使用した（proteins.fasta）。
詳細は以下のサイトから。
コマンド
Braker1 (RNA-seq only)
braker.pl --threads 20 --species=Mang_china \ --genome=JALDXJ000000000_Mang_assembly.sm.rename.fasta \ --bam=Mang_china_RNAseq/DRR328565.bam,\ Mang_china_RNAseq/SRR5997795.bam,\ Mang_china_RNAseq/SRR7586852.bam,\ Mang_china_RNAseq/SRR7586854.bam,\ Mang_china_RNAseq/SRR14160170.bam,\ Mang_china_RNAseq/SRR14160176.bam
Braker2 (protein sequences only)
braker.pl --threads 20 --species=Mang_china \ --genome=JALDXJ000000000_Mang_assembly.sm.rename.fasta \ --prot_seq=proteins.fasta
Braker3 (RNA-seq & protein sequences)
braker.pl --threads 20 --species=Mang_china \ --genome=JALDXJ000000000_Mang_assembly.sm.rename.fasta \ --prot_seq=proteins.fasta \ --bam=Mang_china_RNAseq/DRR328565.bam,\ Mang_china_RNAseq/SRR5997795.bam,\ Mang_china_RNAseq/SRR7586852.bam,\ Mang_china_RNAseq/SRR7586854.bam,\ Mang_china_RNAseq/SRR14160170.bam,\ Mang_china_RNAseq/SRR14160176.bam
解析結果の比較
Number of genes/transcripts
| algorithm | N_genes | N_transcripts |
|---|---|---|
| Braker1 | 45159 | 54210 |
| Braker2 | 48408 | 52575 |
| Braker3 | 30882 | 35050 |
Braker1, Braker2と比較して、Braker3ではgenes, transcriptsの数がかなり少なかった。単一の遺伝子を複数の異なる遺伝子として予測してしまうケースが減少しているのかもしれない。
Completeness (BUSCO5) %, n=3,640
| Complete (single + duplicated) | Complete (Single-copy) | Complete (Duplicated) | Fragmented | Missing | |
|---|---|---|---|---|---|
| Braker1 | 93.9 | 72.6 | 21.3 | 0.7 | 5.4 |
| Braker2 | 90.9 | 82.2 | 8.7 | 3.2 | 5.9 |
| Braker3 | 93.8 | 70.4 | 23.4 | 2.0 | 4.2 |
Braker1とBraker3のcompletenessの数値は大きく違わなかったが、Braker1の方がDuplicated BUSCOsの割合はやや少なかった。Braker2は、Braker1, Braker3と比較すると、completenessがやや低かったが、duplicated BUSCOsの割合はBraker1, Braker3よりもかなり低かった。また、missingの割合はBraker3がもっとも低かった。
Output GTF fileの形式の違い
今回の解析で気付いたことの1つだが、出力されるGTFファイルの形式が、Braker3のみ少し違っていた。Braker3では、geneおよびtranscriptを示す行がなく、各transcriptが"###"で区切られて出力される形式になっていた（以下の例）。
# Braker3 ### CM050755.1 AUGUSTUS start_codon 220137 220139 . + 0 transcript_id "g5.t1"; gene_id "g5"; supported "False"; CM050755.1 AUGUSTUS CDS 220137 224024 0.72 + 0 transcript_id "g5.t1"; gene_id "g5"; cds_type "single"; CM050755.1 AUGUSTUS stop_codon 224022 224024 . + 0 transcript_id "g5.t1"; gene_id "g5"; supported "True"; ### CM050755.1 AUGUSTUS start_codon 239666 239668 . + 0 transcript_id "g6.t1"; gene_id "g6"; supported "False"; CM050755.1 AUGUSTUS CDS 239666 239817 0.92 + 0 transcript_id "g6.t1"; gene_id "g6"; cds_type "initial"; CM050755.1 AUGUSTUS CDS 240465 240543 0.91 + 1 transcript_id "g6.t1"; gene_id "g6"; cds_type "terminal"; CM050755.1 AUGUSTUS intron 239818 240464 0.92 + 0 transcript_id "g6.t1"; gene_id "g6"; supported "True"; CM050755.1 AUGUSTUS stop_codon 240541 240543 . + 0 transcript_id "g6.t1"; gene_id "g6"; supported "False"; # Braker1 & Braker2 CM050755.1 AUGUSTUS gene 34204 47808 . + . g1 CM050755.1 AUGUSTUS transcript 34204 47808 0.99 + . g1.t1 CM050755.1 AUGUSTUS start_codon 34204 34206 . + 0 transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS CDS 34204 34495 1 + 0 transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS exon 34204 34495 . + . transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS intron 34496 47032 1 + . transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS CDS 47033 47140 1 + 2 transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS exon 47033 47140 . + . transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS intron 47141 47278 1 + . transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS CDS 47279 47416 1 + 2 transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS exon 47279 47416 . + . transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS intron 47417 47587 1 + . transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS CDS 47588 47808 0.99 + 2 transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS exon 47588 47808 . + . transcript_id "g1.t1"; gene_id "g1"; CM050755.1 AUGUSTUS stop_codon 47806 47808 . + 0 transcript_id "g1.t1"; gene_id "g1";
