macでインフォマティクス

macでインフォマティクス

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

たくさんのスモールゲノムを比較したり、複数メタゲノムアセンブリのde-replicationを行う dRep

2019 5/7 インストール追記

20196/16 パラメータ追記

20196/16  upしたdocker イメージのエラー修正

2019 6/18 link追加

 

 メタゲノム研究により、シーケンシングされ、ドラフト品質ゲノムが解読される微生物ゲノムの数は毎年急速に拡大している。大きなゲノムセットを包括的に比較するための迅速なアルゴリズムが開発されているが、ドラフト品質のゲノムでは正確ではない。ここでは、不正確だがゲノム距離の迅速かつ推定、および、正確だが遅いANIの計算を順番に適用することによってペアワイズゲノム比較の計算時間を短縮するプログラムdRepを提示する。 dRepは、以前に開発されたアルゴリズムに対してベンチマークされた場合、パーフェクトなrecallとprecisionを維持しながら28倍の速度向上を実現する。我々(著者ら)は、timeseriesデータセットからのゲノムリカバリーでdRepを実証する。時系列の各メタゲノムデータセットは別々にアセンブリされ、それらから、dRepを用いて同一のゲノムグループを同定した。この手順による結果は、時系列データのco-­assemblyでリカバリーされたゲノムセットと比べ、はるかにクオリティの高いゲノムリカバリーを達成していた。

 

drepマニュアル

https://drep.readthedocs.io/en/latest/overview.html

 

f:id:kazumaxneo:20180921221120p:plain

マニュアルより転載。 

 

論文のfull textにはアクセスできなかったので分からないが、Prepinrt*1では、よく使われる生後1ヶ月の未熟児の腸内細菌叢 (pubmed)の時系列メタゲノムシーケンシングデータセットSRA link)を用いて、co-aassemblyと、個別アセンブリ+ dRep(によるde-replication)のアセンブリを比較している。Prepirntの実験結果は、個別アセンブリ+ dRepの方がco-assemblyとその後のanvi’o を使ったマニュアルキュレーションより良好な結果を示している(下の図)。

f:id:kazumaxneo:20181028210958j:plain

Prepinrt(*1)のFigure2を転載。

 

dRepに関するツイート


インストール

依存

  • python3(drep -> python >=3.6,<3.7.0a0
  • Mash is used to rapidly compare all genomes in a pair-wise manner(紹介
  • MUMmer is used to perform more actuate comparisons between genomes which are shown to be similar with Mash(紹介

Optional

Accessory

  • Centrifuge can be used to perform rough taxonomic assignment of bins(紹介

参考

checkmはpython2.7系のコード(python >= 2.7 and < 3.0)、他はpython3に移行しているので、pyenvでpythonのバージョンを管理しているなら、pyenvで両方をglobalにして使う必要がある。

例えばminiconda2-4.0.5とminiconda3-4.0.5をpyenvで管理しているなら

> pyenv global miniconda2-4.0.5 miniconda3-4.0.5

#ANIcalculator
wget https://ani.jgi-psf.org/download_files/ANIcalculator_v1.tgz
tar -zxvf ANIcalculator_v1.tgz
cd ANIcalculator_v1
cp ANIcalculator /usr/local/bin/


git clone https://github.com/abadona/qsimscan.git
cd qsimscan/
make -j 8
#パスを通しておく

本体 Github

#Bioconda (link)
conda install -c bioconda drep

#pip
pip install drep

#インストール確認
dRep bonus output_directory --check_dependencies

dRep bonus output_directory --check_dependencies

Loading work directory

Checking dependencies

mash.................................... all good        (location = /usr/local/bin/mash)

nucmer.................................. all good        (location = /opt/conda/bin/nucmer)

checkm.................................. all good        (location = /usr/local/bin/checkm)

ANIcalculator........................... all good        (location = /gANI/current/ANIcalculator)

prodigal................................ all good        (location = /opt/conda/bin/prodigal)

centrifuge.............................. all good        (location = /usr/local/bin/centrifuge)

(

O.K (*バージョンアップに伴い現在は必要なものが変わってきています)

> dRep

               ...::: dRep v2.0.5 :::...

 

  Choose one of the operations below for more detailed help.

  Example: dRep dereplicate -h

 

  Workflows:

    compare         -> Perform rapid pair-wise comparison on a list of genomes

    dereplicate     -> De-replicate a list of genomes

 

  Single operations:

    filter          -> Filter a genome list based on size, completeness, and/or contamination

    cluster         -> Compare and cluster a genome list based on MASH and ANIn/gANI

    choose          -> Choose the best genome from each genome cluster

    evaluate        -> Evaluate genome de-replication

    bonus           -> Other random operations (determine taxonomy / check dependencies)

    analyze         -> Make figures

    

> dRep compare -h

# dRep compare -h

usage: dRep compare [-p PROCESSORS] [-d] [-h] [-ms MASH_SKETCH]

                    [--S_algorithm {gANI,ANIn,ANImf}]

                    [-n_PRESET {normal,tight}] [-pa P_ANI] [-sa S_ANI]

                    [--SkipMash] [--SkipSecondary] [-nc COV_THRESH]

                    [-cm {total,larger}] [--clusterAlg CLUSTERALG] [--run_tax]

                    [--tax_method {percent,max}] [-per PERCENT]

                    [--cent_index CENT_INDEX] [--warn_dist WARN_DIST]

                    [--warn_sim WARN_SIM] [--warn_aln WARN_ALN]

                    [-g [GENOMES [GENOMES ...]]]

                    work_directory

 

positional arguments:

  work_directory        Directory where data and output    

                        *** USE THE SAME WORK DIRECTORY FOR ALL DREP OPERATIONS ***

 

SYSTEM PARAMETERS:

  -p PROCESSORS, --processors PROCESSORS

                        threads (default: 6)

  -d, --debug           make extra debugging output (default: False)

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

 

GENOME COMPARISON PARAMETERS:

  -ms MASH_SKETCH, --MASH_sketch MASH_SKETCH

                        MASH sketch size (default: 1000)

  --S_algorithm {gANI,ANIn,ANImf}

                        Algorithm for secondary clustering comaprisons:

                        ANImf = (RECOMMENDED) Align whole genomes with nucmer; filter alignment; compare aligned regions

                        ANIn  = Align whole genomes with nucmer; compare aligned regions

                        gANI  = Identify and align ORFs; compare aligned ORFS

                         (default: ANImf)

  -n_PRESET {normal,tight}

                        Presets to pass to nucmer

                        tight   = only align highly conserved regions

                        normal  = default ANIn parameters (default: normal)

 

CLUSTERING PARAMETERS:

  -pa P_ANI, --P_ani P_ANI

                        ANI threshold to form primary (MASH) clusters

                        (default: 0.9)

  -sa S_ANI, --S_ani S_ANI

                        ANI threshold to form secondary clusters (default:

                        0.99)

  --SkipMash            Skip MASH clustering, just do secondary clustering on

                        all genomes (default: False)

  --SkipSecondary       Skip secondary clustering, just perform MASH

                        clustering (default: False)

  -nc COV_THRESH, --cov_thresh COV_THRESH

                        Minmum level of overlap between genomes when doing

                        secondary comparisons (default: 0.1)

  -cm {total,larger}, --coverage_method {total,larger}

                        Method to calculate coverage of an alignment

                        (for ANIn/ANImf only; gANI can only do larger method)

                        total   = 2*(aligned length) / (sum of total genome lengths)

                        larger  = max*1

                         (default: larger)

  --clusterAlg CLUSTERALG

                        Algorithm used to cluster genomes (passed to

                        scipy.cluster.hierarchy.linkage (default: average)

 

TAXONOMY:

  --run_tax             generate taxonomy information (Tdb) (default: False)

  --tax_method {percent,max}

                        Method of determining taxonomy

                        percent = The most descriptive taxonimic level with at least (per) hits

                        max     = The centrifuge taxonomic level with the most overall hits (default: percent)

  -per PERCENT, --percent PERCENT

                        minimum percent for percent method (default: 50)

  --cent_index CENT_INDEX

                        path to centrifuge index (for example,

                        /home/mattolm/download/centrifuge/indices/b+h+v

                        (default: None)

 

WARNINGS:

  --warn_dist WARN_DIST

                        How far from the threshold to throw cluster warnings

                        (default: 0.25)

  --warn_sim WARN_SIM   Similarity threshold for warnings between dereplicated

                        genomes (default: 0.98)

  --warn_aln WARN_ALN   Minimum aligned fraction for warnings between

                        dereplicated genomes (ANIn) (default: 0.25)

 

I/O PARAMETERS:

  -g [GENOMES [GENOMES ...]], --genomes [GENOMES [GENOMES ...]]

                        genomes to cluster in .fasta format (default: None)

 

Example: dRep compare output_dir/ -g /path/to/genomes/*.fasta

dRep dereplicate -h

# dRep dereplicate -h

usage: dRep dereplicate [-p PROCESSORS] [-d] [-h] [-l LENGTH]

                        [-comp COMPLETENESS] [-con CONTAMINATION]

                        [--noQualityFiltering] [-ms MASH_SKETCH]

                        [--S_algorithm {gANI,ANImf,ANIn}]

                        [-n_PRESET {normal,tight}] [-pa P_ANI] [-sa S_ANI]

                        [--SkipMash] [--SkipSecondary] [-nc COV_THRESH]

                        [-cm {total,larger}] [--clusterAlg CLUSTERALG]

                        [-comW COMPLETENESS_WEIGHT]

                        [-conW CONTAMINATION_WEIGHT]

                        [-strW STRAIN_HETEROGENEITY_WEIGHT] [-N50W N50_WEIGHT]

                        [-sizeW SIZE_WEIGHT] [--run_tax]

                        [--tax_method {percent,max}] [-per PERCENT]

                        [--cent_index CENT_INDEX] [--warn_dist WARN_DIST]

                        [--warn_sim WARN_SIM] [--warn_aln WARN_ALN]

                        [-g [GENOMES [GENOMES ...]]]

                        [--checkM_method {taxonomy_wf,lineage_wf}]

                        [--genomeInfo GENOMEINFO]

                        work_directory

 

positional arguments:

  work_directory        Directory where data and output    

                        *** USE THE SAME WORK DIRECTORY FOR ALL DREP OPERATIONS ***

 

SYSTEM PARAMETERS:

  -p PROCESSORS, --processors PROCESSORS

                        threads (default: 6)

  -d, --debug           make extra debugging output (default: False)

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

 

FILTERING OPTIONS:

  -l LENGTH, --length LENGTH

                        Minimum genome length (default: 50000)

  -comp COMPLETENESS, --completeness COMPLETENESS

                        Minumum genome completeness (default: 75)

  -con CONTAMINATION, --contamination CONTAMINATION

                        Maximum genome contamination (default: 25)

  --noQualityFiltering  Don't run checkM or do any quality filtering- will

                        ignore con and comp settings if genomeInfo.csv not

                        also provided. NOT RECOMMENDED! See docs for details

                        (default: False)

 

GENOME COMPARISON PARAMETERS:

  -ms MASH_SKETCH, --MASH_sketch MASH_SKETCH

                        MASH sketch size (default: 1000)

  --S_algorithm {gANI,ANImf,ANIn}

                        Algorithm for secondary clustering comaprisons:

                        ANImf = (RECOMMENDED) Align whole genomes with nucmer; filter alignment; compare aligned regions

                        ANIn  = Align whole genomes with nucmer; compare aligned regions

                        gANI  = Identify and align ORFs; compare aligned ORFS

                         (default: ANImf)

  -n_PRESET {normal,tight}

                        Presets to pass to nucmer

                        tight   = only align highly conserved regions

                        normal  = default ANIn parameters (default: normal)

 

CLUSTERING PARAMETERS:

  -pa P_ANI, --P_ani P_ANI

                        ANI threshold to form primary (MASH) clusters

                        (default: 0.9)

  -sa S_ANI, --S_ani S_ANI

                        ANI threshold to form secondary clusters (default:

                        0.99)

  --SkipMash            Skip MASH clustering, just do secondary clustering on

                        all genomes (default: False)

  --SkipSecondary       Skip secondary clustering, just perform MASH

                        clustering (default: False)

  -nc COV_THRESH, --cov_thresh COV_THRESH

                        Minmum level of overlap between genomes when doing

                        secondary comparisons (default: 0.1)

  -cm {total,larger}, --coverage_method {total,larger}

                        Method to calculate coverage of an alignment

                        (for ANIn/ANImf only; gANI can only do larger method)

                        total   = 2*(aligned length) / (sum of total genome lengths)

                        larger  = max*2

                         (default: larger)

  --clusterAlg CLUSTERALG

                        Algorithm used to cluster genomes (passed to

                        scipy.cluster.hierarchy.linkage (default: average)

 

SCORING CRITERIA

Based off of the formula: 

A*Completeness - B*Contamination + C*(Contamination * (strain_heterogeneity/100)) + D*log(N50) + E*log(size)

 

A = completeness_weight; B = contamination_weight; C = strain_heterogeneity_weight; D = N50_weight; E = size_weight:

  -comW COMPLETENESS_WEIGHT, --completeness_weight COMPLETENESS_WEIGHT

                        completeness weight (default: 1)

  -conW CONTAMINATION_WEIGHT, --contamination_weight CONTAMINATION_WEIGHT

                        contamination weight (default: 5)

  -strW STRAIN_HETEROGENEITY_WEIGHT, --strain_heterogeneity_weight STRAIN_HETEROGENEITY_WEIGHT

                        strain heterogeneity weight (default: 1)

  -N50W N50_WEIGHT, --N50_weight N50_WEIGHT

                        weight of log(genome N50) (default: 0.5)

  -sizeW SIZE_WEIGHT, --size_weight SIZE_WEIGHT

                        weight of log(genome size) (default: 0)

 

TAXONOMY:

  --run_tax             generate taxonomy information (Tdb) (default: False)

  --tax_method {percent,max}

                        Method of determining taxonomy

                        percent = The most descriptive taxonimic level with at least (per) hits

                        max     = The centrifuge taxonomic level with the most overall hits (default: percent)

  -per PERCENT, --percent PERCENT

                        minimum percent for percent method (default: 50)

  --cent_index CENT_INDEX

                        path to centrifuge index (for example,

                        /home/mattolm/download/centrifuge/indices/b+h+v

                        (default: None)

 

WARNINGS:

  --warn_dist WARN_DIST

                        How far from the threshold to throw cluster warnings

                        (default: 0.25)

  --warn_sim WARN_SIM   Similarity threshold for warnings between dereplicated

                        genomes (default: 0.98)

  --warn_aln WARN_ALN   Minimum aligned fraction for warnings between

                        dereplicated genomes (ANIn) (default: 0.25)

 

I/O PARAMETERS:

  -g [GENOMES [GENOMES ...]], --genomes [GENOMES [GENOMES ...]]

                        genomes to cluster in .fasta format (default: None)

  --checkM_method {taxonomy_wf,lineage_wf}

                        Either lineage_wf (more accurate) or taxonomy_wf

                        (faster) (default: lineage_wf)

  --genomeInfo GENOMEINFO

                        location of .csv file containing quality information

                        on the genomes. Must contain: ["genome"(basename of

                        .fasta file of that genome), "completeness"(0-100

                        value for completeness of the genome),

                        "contamination"(0-100 value of the contamination of

                        the genome)] (default: None)

 

Example: dRep dereplicate output_dir/ -g /path/to/genomes/*.fasta

 

導入に手間取ったので、環境構築した後のイメージをpdocker hubにpushしておきます。pullする場合、別のツールのテストもしていたので、少しイメージが大きくなっています。ご注意下さい。

追記

checkmのstepでエラーが出たのでBug fixしてpushし直しました。

docker pull kazumax/drep

#カレントパスと/shareを共有して立ち上げる。"--rm"をつけるとexitで廃棄
docker run --rm -itv $PWD:/data/ kazumax/drep

#usage:
source ${HOME}/.bash_profile
#check dependency
dRep bonus output_directory --check_dependencies

 

実行方法

1、ゲノムの比較

dRep compare output_directory -g RefSeq/*.fasta

結果は可視化される。

 

RefSeqのゲノム(ダウンロード)をいくつか選んで比較してみた。

出力ディレクト

f:id:kazumaxneo:20181028150243j:plain

f:id:kazumaxneo:20181028150232j:plain

f:id:kazumaxneo:20181028150247j:plain

f:id:kazumaxneo:20181028150242j:plain

出力の解説リンク

 

 

2、de-replication

dRep dereplicate outout_directory -g path/to/genomes/*.fasta \
-p 12 -l 50000
  • -p     threads (default: 6)
  • -l      Minimum genome length (default: 50000)
  • -g     genomes to cluster in .fasta format (default: None)
  • --checkM_method {taxonomy_wf, lineage_wf} Either lineage_wf (more accurate) or taxonomy_wf (faster) (default: lineage_wf)

結果は可視化される。

出力ディレクト

f:id:kazumaxneo:20181028150243j:plain

figures

f:id:kazumaxneo:20190617185427j:plain

f:id:kazumaxneo:20190617185452j:plain

f:id:kazumaxneo:20190617185459j:plain

f:id:kazumaxneo:20190617185513j:plain


出力の詳細はmanual参照。

de-replicationは例えば以下の有名な論文で使用されています。パラメータも記載されています。

https://www.nature.com/articles/s41467-018-03317-6

引用

dRep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replicatio
Olm MR, Brown CT, Brooks B, Banfield JF

ISME J. 2017 Dec;11(12):2864-2868

 

*1

Preprint

https://www.biorxiv.org/content/biorxiv/early/2017/02/13/108142.full.pdf

 

関連ツール

BMScan

Mash

Checkm

Centrifuge

 fastANI


 

 

*1:aligned length / genome 1), (aligned_length / genome2

*2:aligned length / genome 1), (aligned_length / genome2