macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

アセンブル結果をCore gene setの検出数で評価する BUSCO

2018 1/9 version3のコマンドを最後に追記

2019 2/24 論文追記

2019 7/5  verrsion3向けに説明をアップデート

2019 11/24 論文追記

2019 12/26 v4インストール追記

2020 3/26 v4 追記

2020 6/15 構成を変更

2020 7/7 v4 のtrancrripts の説明が間違っていたのを修正

2020 7/22 誤字修正

2020 8/26 docker link追記2020 9/27 最新の情報を一番上に移動

2021 2/5 busco5のインストール追記,  dockerのコマンド修正

 

 ゲノムのアセンブルやde novo transcriptomeの評価手法の1つに、Core gene setがアセンブルされた配列の中にどれだけあるか調べる方法がある(core genesは構成的に発現していると考える)。そのようなツールとしてCEGMAがよく知られている。CEGMAはversion2.5までアップーデートされていたが、2015年にはその後継のBUSCOが発表され、CEGMAはサポートが中止された。BUSCOはCEGMA3.0とも言うべきツールで、CEGMAより随分高速化されている(CEGMAサポート停止の経緯 Goodbye CEGMA, hello BUSCO! — ACGT)。ここではBUSCOを紹介する。

 BUSCOは、90%以上の種でsingle copy orthologousな遺伝子をcore geneとして定義し、それをデータベースにしてクエリー配列を検索する。100%でなく90パーセントで線を引いたのは、一部の種ではcore geneが複製したり、アセンブルが不完全で抜けているため見つからなかったなどの理由で100%保存にできなかったためらしい。そのため90という数値に生物的な意味はないが、core geneヒット数が90パーセントを大きく割るようならおそらくアセンブルが不完全と判断できる。

 

公式ページ

https://busco.ezlab.org

ランマニュアル

https://busco.ezlab.org/busco_userguide.html

 

インストール

GitLab

#bioconda (link)バージョン指定なしなら最新の4が入る(テスト時はバージョン指定しないとうまく入らなかった)
mamba create -n busco4 python=3.7 -y
conda activate busco4
mamba install -c bioconda -c conda-forge busco==4.1.3 -y

#v5 (link)2021 5/19現在5.1.2が最新
mamba create -n busco5 python=3.7 -y
conda activate busco5
mamba install -c bioconda -c conda-forge busco==5.1.2 -y

#docker (dockerhub link)(v4 latest)
docker pull ezlabgva/busco:v4.1.2_cv1
#run
docker run -itv $PWD:/data -w /data --rm ezlabgva/busco:v4.1.2_cv1

> busco -h

$  busco -h

usage: busco -i [SEQUENCE_FILE] -l [LINEAGE] -o [OUTPUT_NAME] -m [MODE] [OTHER OPTIONS]

 

Welcome to BUSCO 4.0.5: the Benchmarking Universal Single-Copy Ortholog assessment tool.

For more detailed usage information, please review the README file provided with this distribution and the BUSCO user guide.

 

optional arguments:

  -i FASTA FILE, --in FASTA FILE

                        Input sequence file in FASTA format. Can be an assembled genome or transcriptome (DNA), or protein sequences from an annotated gene set.

  -c N, --cpu N         Specify the number (N=integer) of threads/cores to use.

  -o OUTPUT, --out OUTPUT

                        Give your analysis run a recognisable short name. Output folders and files will be labelled with this name. WARNING: do not provide a path

  --out_path OUTPUT_PATH

                        Optional location for results folder, excluding results folder name. Default is current working directory.

  -e N, --evalue N      E-value cutoff for BLAST searches. Allowed formats, 0.001 or 1e-03 (Default: 1e-03)

  -m MODE, --mode MODE  Specify which BUSCO analysis mode to run.

                        There are three valid modes:

                        - geno or genome, for genome assemblies (DNA)

                        - tran or transcriptome, for transcriptome assemblies (DNA)

                        - prot or proteins, for annotated gene sets (protein)

  -l LINEAGE, --lineage_dataset LINEAGE

                        Specify the name of the BUSCO lineage to be used.

  -f, --force           Force rewriting of existing files. Must be used when output files with the provided name already exist.

  --limit REGION_LIMIT  How many candidate regions (contig or transcript) to consider per BUSCO (default: 3)

  --long                Optimization mode Augustus self-training (Default: Off) adds considerably to the run time, but can improve results for some non-model organisms

  -q, --quiet           Disable the info logs, displays only errors

  --augustus_parameters AUGUSTUS_PARAMETERS

                        Pass additional arguments to Augustus. All arguments should be contained within a single pair of quotation marks, separated by commas. E.g. '--param1=1,--param2=2'

  --augustus_species AUGUSTUS_SPECIES

                        Specify a species for Augustus training.

  --auto-lineage        Run auto-lineage to find optimum lineage path

  --auto-lineage-prok   Run auto-lineage just on non-eukaryote trees to find optimum lineage path

  --auto-lineage-euk    Run auto-placement just on eukaryote tree to find optimum lineage path

  --update-data         Download and replace with last versions all lineages datasets and files necessary to their automated selection

  --offline             To indicate that BUSCO cannot attempt to download files

  --config CONFIG_FILE  Provide a config file

  -v, --version         Show this version and exit

  -h, --help            Show this help message and exit

  --list-datasets       Print the list of available BUSCO datasets

 

 

 

version4のラン方法

1、利用できるデータベースを確認

busco --list-datasets

 

 2、ラン

ここではeukaryota_odb10データベースを指定する。

#genome
busco -m geno -i genome.fasta -o out_dir -l eukaryota_odb10 -c 20

#transcriptome
busco -m tran -i transcriptome.fasta -o out_dir -l eukaryota_odb10 -c 20

#proteome
busco -m prot -i proteome.fasta -o out_dir -l eukaryota_odb10 -c 20
  • -m    Specify which BUSCO analysis mode to run.
    There are three valid modes:
    - geno or genome, for genome assemblies (DNA)
    - tran or transcriptome, for transcriptome assemblies (DNA)
     - prot or proteins, for annotated gene sets (protein)
  • -c   Specify the number (N=integer) of threads/cores to use.

 

2021 08/14

version5は改良点が多いので別にまとめています。更新されたデータベースを使っており、真菌などではかなり結果が変わる可能性があります。詳しくは論文を読んでください。


 

引用

BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs
Felipe A. Simão, Robert M. Waterhouse, Panagiotis Ioannidis, Evgenia V. Kriventseva, and Evgeny M. Zdobnov
Bioinformatics, published online June 9, 2015

 

BUSCO: Assessing Genome Assembly and Annotation Completeness
Seppey M, Manni M, Zdobnov EM

Methods Mol Biol. 2019;1962:227-245. doi: 10.1007/978-1-4939-9173-0_14

 

関連

アセンブリmetricsを調べるツールには、TransrateやDRAPなどもあります。DRAPにはBUSCOも取り込まれています。


gVolanteではBUSCOv4とv5が使えるようになっています。

web上で簡単に使うことができます。

 

 

Insect Genomics

Using BUSCO to Assess Insect Genomic Resources

https://link.springer.com/protocol/10.1007/978-1-4939-8775-7_6

 

追記

Niwa Ryoさんがv4の使い方を説明されています。

 

 

 

===========================================================

以下は古いバージョンの情報

公式ページ 解析対象の生物に応じてHPからデータベースをダウンロードする。

=> バージョン4では指定データベースが自動ダウンロードされるため、手動での出すんロードは不要になっています。一番下のversion4を確認して下さい。

https://busco.ezlab.org

f:id:kazumaxneo:20180410124024j:plain

HPからはマニュアルもダウンロードできるようになっている。 

 

インストール

v4を導入

#2019 12/26追記 v4
conda install -c bioconda -y busco


#homebrew(インストール未確認)
brew install busco

*CEGMAもbrewでインストール可能。

VMのイメージもあります(リンク)。 

 > run_BUSCO.py -h

# run_BUSCO.py -h

usage: python BUSCO.py -i [SEQUENCE_FILE] -l [LINEAGE] -o [OUTPUT_NAME] -m [MODE] [OTHER OPTIONS]

 

Welcome to BUSCO 3.0.2: the Benchmarking Universal Single-Copy Ortholog assessment tool.

For more detailed usage information, please review the README file provided with this distribution and the BUSCO user guide.

 

optional arguments:

  -i FASTA FILE, --in FASTA FILE

                        Input sequence file in FASTA format. Can be an assembled genome or transcriptome (DNA), or protein sequences from an annotated gene set.

  -c N, --cpu N         Specify the number (N=integer) of threads/cores to use.

  -o OUTPUT, --out OUTPUT

                        Give your analysis run a recognisable short name. Output folders and files will be labelled with this name. WARNING: do not provide a path

  -e N, --evalue N      E-value cutoff for BLAST searches. Allowed formats, 0.001 or 1e-03 (Default: 1e-03)

  -m MODE, --mode MODE  Specify which BUSCO analysis mode to run.

                        There are three valid modes:

                        - geno or genome, for genome assemblies (DNA)

                        - tran or transcriptome, for transcriptome assemblies (DNA)

                        - prot or proteins, for annotated gene sets (protein)

  -l LINEAGE, --lineage_path LINEAGE

                        Specify location of the BUSCO lineage data to be used.

                        Visit http://busco.ezlab.org for available lineages.

  -f, --force           Force rewriting of existing files. Must be used when output files with the provided name already exist.

  -r, --restart         Restart an uncompleted run. Not available for the protein mode

  -sp SPECIES, --species SPECIES

                        Name of existing Augustus species gene finding parameters. See Augustus documentation for available options.

  --augustus_parameters AUGUSTUS_PARAMETERS

                        Additional parameters for the fine-tuning of Augustus run. For the species, do not use this option.

                        Use single quotes as follow: '--param1=1 --param2=2', see Augustus documentation for available options.

  -t PATH, --tmp_path PATH

                        Where to store temporary files (Default: ./tmp/)

  --limit REGION_LIMIT  How many candidate regions (contig or transcript) to consider per BUSCO (default: 3)

  --long                Optimization mode Augustus self-training (Default: Off) adds considerably to the run time, but can improve results for some non-model organisms

  -q, --quiet           Disable the info logs, displays only errors

  -z, --tarzip          Tarzip the output folders likely to contain thousands of files

  --blast_single_core   Force tblastn to run on a single core and ignore the --cpu argument for this step only. Useful if inconsistencies when using multiple threads are noticed

  -v, --version         Show this version and exit

  -h, --help            Show this help message and exit

 

 

注意

以下は旧バージョンの解説です。最新のv4の使い方は上のversion4を見てください。データベースの手動ダウンロードは必要なくなっています。

データベース

解析には、Core geneのデータベースを用意する必要がある。データセットは以下のリンクからダウンロードする。バクテリアも用意されている。

 

シロイヌナズナのデータなら陸上植物のデータベースを使う。wgetでダウンロードする。

wget http://busco.ezlab.org/v2/datasets/embryophyta_odb9.tar.gz
tar -zxvf embryophyta_odb9.tar.gz

ダウンロードが終わったら解凍する。遺伝子リストは解凍したディレクトリのinfo/~info.txt.gzに保存されている。中身を見てみる。

user$ gzip -dc embryophyta_3193_OrthoDB9_orthogroup_info.txt.gz |head -10 |cut -f 1-2

OrthoGroupID Description

EOG09360002 "Zinc finger, UBR-type"

EOG0936000A "Zinc finger, RING/FYVE/PHD-type"

EOG0936001N "PIK-related kinase"

EOG0936001W "Uncharacterized protein"

EOG09360025 "SNF2-related, N-terminal domain"

EOG0936003Z "P-loop containing nucleoside triphosphate hydrolases superfamily protein"

EOG0936004R "Telomere length regulation protein, conserved domain"

EOG0936004W "DNA-directed DNA polymerase, family B"

EOG0936008E "Zinc finger, RING/FYVE/PHD-type"

左端1、2列目にIDとdescriptionがある。植物のcore geneデータベースは1440遺伝子ある。この解凍してできたembryophyta_odb9/をランする時に指定する ("-l embryophyta_odb9")。

 

実行方法

Genome assembly

run_BUSCO.py -i scaffolds.fasta -o OUTPUT -l embryophyta_odb9 -m genome -c 8
  •  -i   Input sequence file in FASTA format. Can be an assembled
  • -o   output
  • -l  lineage,lineage Which BUSCO lineage to be used.(ここではembryophyta_odb9を指定)
  • -m      Specify which BUSCO analysis mode to run.
                            There are three valid modes:
                            - geno or genome, for genome assemblies (DNA)
                            - tran or transcriptome, for transcriptome assemblies (DNA)
                            - prot or proteins, for annotated gene sets (protein)
  • -c Number of threads/cores to use.

 

Transcriptome 

run_BUSCO.py -i transcripts.fasta -o OUTPUT -l embryophyta_odb9 -m transcriptome -c 8

アラビのRNAをde novoでアセンブルした配列(fasta)をクエリーとしている。

  • -l ここではembryophyta_odb9を指定

アラビのtranscriptomeデータをde novo assemblyして作った配列(54000 transcrpts)を解析したところ、10分くらいかかった。

 

full_table_OUTPUT~が詳細なデータとなる。またshort_summary_OUTPUT_PLANTに結果がまとめられている。summaryには以下のような情報が載っている。

#Summarized BUSCO benchmarking for file: transcripts.fasta 

#BUSCO was run in mode: trans

Summarized benchmarks in BUSCO notation

 : C:0%[D:0%],F:0%,M:0%,n:1440

 1357 Complete BUSCOs

 896 Complete and single-copy BUSCOs

 461 Complete and duplicated BUSCOs

 17 Fragmented BUSCOs

 66 Missing BUSCOs

 1440 Total BUSCO groups searched

 

植物用のcore geneのカタログは合計1440遺伝子あるが、1357見つかり、そのうち461はduplicateしているとの結果が出た。また、1440遺伝子のうち66が見つからないと出ている。

missing_buscos_list~には見つからなかった66のIDが載っている。orthoDBを検索して、どんな遺伝子がsingle copy orthologousとして見つからなかったか確認する。

Ortholog Search | cegg.unige.ch Computational Evolutionary Genomics Group

 

 今回の検証ではシロイヌナズナのリファレンスcDNAを使っている。ハウスキーピング遺伝子の発現がゼロになるとは考えにくいので、ナズナのゲノムにはmissingの66遺伝子の大半が存在しないということだろうか。

  

ラン途中にエラーが起きたら、fastaのヘッダーが悪さをしている可能性もあります。以下のようなコマンドでシンプルな名前に修正してランし直してみて下さい。

perl -ane 'if(/\>/){$a++;print ">name_$a\n"}else{print;}' input.fa > rename.fa

  

追記 1

version3

version3.02は少しコマンドが変更されている。

ここではdockerhubのイメージを使う(link)。

docker pull comics/busco

例えばバクテリアゲノムのアセンブル結果を評価

#ここでは中に入ってランする
docker run --rm -itv $PWD:/data/ comics/busco
cd /data/

run_BUSCO.py -i assembly.fasta -o busco -l bacteria_odb9/ -m genome -c 1