macでインフォマティクス

macでインフォマティクス

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

de novo transcriptomeのアノテーションツール dammit

2020 1/31 誤字修正

2020 2/1 コマンド修正

 

dammitは、単純なde novo transcriptome annotatorである。 アノテーションのプロセスの個々の部分は全てすでに存在しているが、既存の解決策は過度に複雑であるか、または無駄な非フリーソフトウェアに依存しているという観測から生まれた。

dammitは無料でオープンソースであり、フリーでオープンソースのエコシステムを中心に構築されている。 したがって、著者が十分に考慮していないプログラムは、依存から避けている。 これは、非フリーライセンスのプログラム、またはインストールと構成が非常に難しいプログラムのいずれかを意味する可能性がある。著者らは、アクセスはオープン性の一部であると考えている。

 

dammitはde novo transcriptomeのfunctional annotationを自動で行うことができるツール。pythonの仮想環境で実行するように設計されており、ユーザーは最小限の手間でジョブを実行できる。以前から公開されていたが、Biocondaのサポート前は、導入手順が煩雑で、正常に動かすまで手間だった。最近になってbiocondaのサポートが入り、condaで簡単に導入できるようになった。

 

HP

http://www.camillescott.org/dammit/

 

Annotating de novo transcriptomes with dammit¶

https://angus.readthedocs.io/en/2017/dammit_annotation.html

 

 

dammit annotationに関するツイート


dammitは以下のデータベースを使う(HPより)。

Pfam-A

Pfam-A is a collection of protein domain profiles for use with profile hidden markov model programs like hmmer. These searches are moderately fast and very sensitive, and the Pfam database is very well curated. Pfam is used during TransDecoder’s ORF finding and for annotation assignment.

Rfam

Rfam is a collection of RNA covariance models for use with programs like Infernal. Covariance models describe RNA secondary structure, and Rfam is a curated database of non-coding RNAs.

OrthoDB

OrthoDB is a curated database of orthologous genes. It attempts to classify proteins from all major groups of eukaryotes and trace them back to their ancestral ortholog.

BUSCO

BUSCO databases are collections of “core” genes for major domains of life. They are used with an accompanying BUSCO program which assesses the completeness of a genome, transcriptome, or list of genes. There are multiple BUSCO databases, and which one you use depends on your particular organism (*1).

uniref90

uniref is a curated collection of most known proteins, clustered at a 90% similarity threshold. This database is comprehensive, and thus quite enormous. dammit does not include it by default due to its size, but it can be installed and used with the --full flag.

 

 

インストール

ubuntu16.04に導入した。

依存

dammit, for now, is officially supported on GNU/Linux systems via bioconda. macOS support will be available via bioconda soon.

本体 Github

#Anaconda環境で実行
conda create -y -n dammit python=3 #Python3本体を含めた環境を作成
conda activate dammit #仮想環境をActiveにする
conda install -y -c bioconda -c conda-forge dammit #インストール

#dockerhub(link)(試してません)
docker pull pypi/dammit

dammit -h

$ dammit -h

usage: dammit [-h] [--debug] [--version] {migrate,databases,annotate} ...

 

# dammit: a tool for easy de novo transcriptome annotation

 

optional arguments:

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

  --debug

  --version             show program's version number and exit

 

dammit subcommands:

  {migrate,databases,annotate}

    databases           Check for databases and optionally download and

                        prepare them for use. By default, only check their

                        status.

    annotate            The main annotation pipeline. Calculates assembly

                        stats; runs BUSCO; runs LAST against OrthoDB (and

                        optionally uniref90), HMMER against Pfam, Inferal

                        against Rfam, and Conditional Reciprocal Best-hit

                        Blast against user databases; and aggregates all

                        results in a properly formatted GFF3 file.

 

——

 

データベースの準備

#フルダウンロード(2018年6月現在、9データベースがout of dataしていた)
dammit databases --install

#軽量版(OrthoDB、uniref、Pfam、Rfamが入らない)
dammit databases --install --quick

#BUSCOデータベース(リンク)plantsを導入
dammit databases --install --busco-group plants

 Macではhomeの.dammit/databasesにデータベースが保存された。

 

実行方法

dammit annotate transcriptome.fa --n_threads 12 --busco-group \
eukaryota -o output_dir

 FASTAファイル名のディレクトリができ、そこに結果が出力される。

 

 

配列が多いと、"--full"フラグ付きのランにはかなりの時間がかかります。ご注意下さい。

引用

GitHub - camillescott/dammit: just annotate it, dammit!

 

参考 

*1

BUSCOは以下のデータベースが利用できる。

 dammit databases: error: argument --busco-group: invalid choice: 'plant' (choose from 'tenericutes', 'dikarya', 'pezizomycotina', 'enterobacteriales', 'euarchontoglires', 'basidiomycota', 'firmicutes', 'sordariomyceta', 'bacteria', 'saccharomycetales', 'metazoa', 'proteobacteria', 'eurotiomycetes', 'ascomycota', 'cyanobacteria', 'betaproteobacteria', 'vertebrata', 'gammaproteobacteria', 'lactobacillales', 'aves', 'mammalia', 'actinopterygii', 'insecta', 'microsporidia', 'laurasiatheria', 'deltaepsilonsub', 'saccharomyceta', 'fungi', 'hymenoptera', 'bacteroidetes', 'alveolata_stramenophiles', 'protists', 'actinobacteria', 'spirochaetes', 'clostridia', 'bacillales', 'tetrapoda', 'rhizobiales', 'arthropoda', 'endopterygota', 'eukaryota', 'embryophyta', 'diptera', 'nematoda')

利用の際は

"dammit databases --install --busco-group xxx"

を実行しておく。

 

 

 

関連ツール

 

 

よく似たゲノム情報を使い不完全なゲノム情報しか持たない種のRNA seq解析の精度を上げる自動化されたツール Necklace

2018 10/31、11/2  タイトル、コード等修正、docker追加

2021 3/9 出力例追記

 

 

 シーケンシングされた種の数が増加しているが、ゲノムの大部分は不完全である。それらにはギャップが含まれていても、配置されていない領域が残っていてもよく、アノテーションが不十分な場合もある。ゲノムを持つ種のRNAシーケンシング(RNA-seq)はモデル生物と同じ手順に従う。つまりリードをゲノムにアライメントし、注釈付き遺伝子とオーバーラップするリードをカウントし、リードカウントに基づいて発現変動があるかをテストする[論文より ref. 1]。しかし、このアプローチは、多くの生物にとって重要なbiologyを見逃す可能性がある。リファレンス配列のギャップまたはアノテーションの欠落のために、遺伝子のセグメントが見逃される可能性がある。下流のdifferential expression 解析は、遺伝子数が過小評価されているため、統計分析の能力を低下させる可能性が高い。著者らは、遺伝子が複数のアセンブリscaffoldsにまたがる場合に、異なるアノテーションを持つ遺伝子としてアノテートされる例を観察した。しかし、これは遺伝子が単一scaffolds内に収まる場合でも起こり得る。最悪の場合、遺伝子全体が見逃される可能性がある。

 RNA-seq解析は、発現された遺伝子に関する情報をデータ自体から抽出することで、利用可能な遺伝子モデルを修復するゲノムガイドの方法かde novo assemblyの方法を使うのが理想的である。しかし、アセンブリを含む解析は、複数のアセンブリを統合する必要がある場合に複雑になる(一部略)。

 Davidson et alの論文[ef. 4](ツール紹介)で、著者らはsuperTranscriptomeの概念を導入した。そこでは、各遺伝子はすべてのエクソンを含む1つの配列によって表される。スーパートランスクリプトは、アセンブリアノテーションの異なるソースからのトランスクリプトームをコンパクトで統一されたリファレンスに結合できる便利な手段を提供した。チキンRNA seqに適用すると、チキンリファレンスゲノムには存在しない何百もの遺伝子セグメントを回収できることが示された。

 ここでは、不完全なリファレンスゲノムを持つ種について、[ref.4]で説明したプロセスを自動化するNecklaceというソフトウェアを紹介する。 NecklaceはRNA-seqリードへのパス、リファレンスゲノム、および1つ以上のリファレンスゲノムアノテーションを含む設定ファイルを入力として受け取る。デノボアセンブリはエラーが起こりやすいため、関連する(よくアノテーションけされた)種のコード配列の中に新たにアセンブリされた遺伝子が発見される必要がある。したがって、関連する種のゲノムおよび注釈もNecklaceに提供されなければならない。Necklaceは、ゲノムガイドおよびデノボアセンブリに含まれるステップを実行し、アセンブリされたトランスクリプトームを、ユーザーが関心ある種のアノテーションと組み合わせる。 superTranscriptomeを構築した後、Necklaceは、リードをsuperTranscriptomeにアライメントしてリードをカウントし、edgeR [ref.5]、DEseq [ref.6]、またはDEXseq [ref.7]などのよく確立されたツールによる試験の準備を行う。

 新しいデータセットにNecklaceを適用することを実証するために、著者らはsheep milkの公開RNA-seqデータを分析した。ヒツジリファレンスゲノムのみを使用する場合と比較して、Necklaceを使った分析では、遺伝子にアサインされたリードが18%多く、発現変動遺伝子は19%多く検出された。

 

 

f:id:kazumaxneo:20180624133710j:plain

necklaceパイプライン。論文より転載。

Necklaceは、bpipeフレームワーク[ref.8]を使用して構築されたパイプラインである。 Necklaceはいくつかのサブステージで構成されている。初期ステージのゲノガイドアセンブリおよびde novoアセンブリ、遺伝子群への転写産物のクラスタリング、スーパートランスクリプトームの構築のための再構築、およびリードをマッピングしてカウントし発現変動遺伝子とスプライシングバリアントの分析に備える。アライナーやアセンブラなどの外部ソフトウェアと、C / C ++で書かれた独自のユーティリティセットを管理する。 入力として、Necklaceはraw RNA-seqリードと該当する種のリファレンスゲノムおよびアノテーションファイルをとる。それに加えて、Necklaceは、よく研究されている種であるヒト、ショウジョウバエ酵母などからリファレンスゲノムとアノテーションを取り込む。 

リードカウントまでやってくれます。正規化は自分で実行する必要があります。

 

wiki

https://github.com/Oshlack/necklace/wiki

 

関連ツイート

 

インストール

ubuntu16.04を使ってテストした(docker使用。ホストOS: mac os10.12)。

依存

  • bash command line
  • Java 1.8 - You can also use earlier versions of Java, but need to perform de novo transcriptome assembly outside Necklace. See Options.
  • wget - Only required for installing Necklace and not for running it. If the required tools are installed manually, this is not needed.
  • g++

本体 Github

下記インストールスクリプトを実行するとcondaコマンドを途中実行する。他の依存とのコンフリクトが起きないようにも、仮想環境で実行するのが無難。javaの最新版は2018年現在version10だが、本ツールは8を要求する。以下の方法で1.8を導入した(確認は"java version")

#JDK1.8導入
sudo apt-get install software-properties-common
sudo apt-get update
sudo add-apt-repository ppa:webupd8team/java
sudo apt-get update
sudo apt-get install oracle-java8-installer
java -version

#本体
git clone https://github.com/Oshlack/necklace.git
cd necklace/

#linuxの場合
./install_linux64.sh

#macの場合 (30分ほどかかる)
./install_mac.sh

#パスが表示されるかcheck
cat tools.groovy

$ cat tools.groovy 

// Path to tools used by the pipeline

bpipe="/home/kazu/necklace/tools/bin/bpipe"

hisat2="/home/kazu/necklace/tools/bin/hisat2"

stringtie="/home/kazu/necklace/tools/bin/stringtie"

gffread="/home/kazu/necklace/tools/bin/gffread"

blat="/home/kazu/necklace/tools/bin/blat"

cluster="/home/kazu/necklace/tools/bin/cluster"

python="/home/kazu/necklace/tools/bin/python"

lace="/home/kazu/necklace/tools/bin/lace"

gtf2flatgtf="/home/kazu/necklace/tools/bin/gtf2flatgtf"

samtools="/home/kazu/.pyenv/shims/samtools"

Trinity="/home/kazu/necklace/tools/bin/Trinity"

bowtie2="/home/kazu/necklace/tools/bin/bowtie2"

make_blocks="/home/kazu/necklace/tools/bin/make_blocks"

featureCounts="/home/kazu/necklace/tools/bin/featureCounts"

上では直しているが、samtoolsは導入に失敗していたので、別途導入して、tools.groovyのbowtie2="" 行にパスを記載した (*1)。

> tools/bin/bpipe -h

# tools/bin/bpipe

 

Bpipe Version 0.9.9.2   Built on Mon Jul 18 03:11:39 UTC 2016

 

usage: bpipe [run|test|debug|execute] [options] <pipeline> <in1> <in2>...

             retry [test]

             stop

             history

             log

             jobs

             checks

             status

             cleanup

             query

             preserve

             register <pipeline> <in1> <in2>...

             diagram <pipeline> <in1> <in2>...

             diagrameditor <pipeline> <in1> <in2>...

 -b,--branch <arg>                Comma separated list of branches to

                                  limit execution to

 -d,--dir <arg>                   output directory

 -f,--filename <arg>              output file name of report

 -h,--help                        usage information

 -l,--resource <resource=value>   place limit on named resource

 -L,--interval <arg>              the default genomic interval to execute

                                  pipeline for (samtools format)

 -m,--memory <arg>                maximum memory in MB, or specified as

                                  <n>GB or <n>MB

 -n,--threads <arg>               maximum threads

 -p,--param <param=value>         defines a pipeline parameter, or file of

                                  paramaters via @<file>

 -r,--report                      generate an HTML report / documentation

                                  for pipeline

 -R,--report <arg>                generate report using named template

 -t,--test                        test mode

 -u,--until <arg>                 run until stage given

 -v,--verbose                     print internal logging to standard error

 -y,--yes                         answer yes to any prompts or questions

 

#導入が少し手間だったので、dockerhubにイメージを上げておきます。mac osでも使えるはずです。イメージにはNecklaceの下流解析のツールは入れていません。下流解析もコンテナ中で行うなら、必要なものを入れてcommitし直して使って下さい、

docker pull kazumax/necklace

#ホストのカレントディレクトリとイメージの/dataをシェアして起動
docker run -itv $PWD:/data/ kazumax/necklace
cd /root/necklace
tools/bin/bpipe -h

上のdockerイメージを使っても テストランは正常に終わりますが、ラン中の進捗メッセージが非常に多く、何か見過ごしがあるかもしれません。注意して使って下さい。

 

 

テストラン

本ツールはBpipeというワークフローマネジメントシステムを使って構築されている。詳細はbpipeのペーパーとdocument参照。ランは、bin/のbpipeを叩くことで実行する。このテストランにはmac pro 2012で1時間ほどかかる。

#テストデータのダウンロードと解凍(羊のRNA seq。自身のゲノム以外に、ヒトゲノムをrelated speciesのゲノムとして使っている)
wget https://github.com/Oshlack/necklace/wiki/asserts/sheep_small_demo_data.tar.gz
tar -zxvf sheep_small_demo_data.tar.gz

#ラン
tools/bin/bpipe run necklace.groovy data/data.txt

情報はconfigファイルに記載する。

> cat data/data.txt

$ cat data/data.txt

// sequencing data

reads_R1="data/SRR2932539_small_1.fastq.gz,data/SRR2932540_small_1.fastq.gz,data/SRR2932541_small_1.fastq.gz,data/SRR2932542_small_1.fastq.gz,data/SRR2932561_small_1.fastq.gz,data/SRR2932562_small_1.fastq.gz,data/SRR2932563_small_1.fastq.gz,data/SRR2932564_small_1.fastq.gz"

reads_R2="data/SRR2932539_small_2.fastq.gz,data/SRR2932540_small_2.fastq.gz,data/SRR2932541_small_2.fastq.gz,data/SRR2932542_small_2.fastq.gz,data/SRR2932561_small_2.fastq.gz,data/SRR2932562_small_2.fastq.gz,data/SRR2932563_small_2.fastq.gz,data/SRR2932564_small_2.fastq.gz"

 

//The genome and annotation

annotation="data/Ovis_aries.Oar_v3.1.90.chr14.gtf"

genome="data/Ovis_aries.Oar_v3.1.dna.toplevel.chr14.fa"

 

//The genome and annotation of a related species

annotation_related_species="data/Homo_sapiens.GRCh38.90.CDS.chr19.gtf"

genome_related_species="data/Homo_sapiens.GRCh38.dna.toplevel.chr19.fa"

例えばfastqディレクトリに2サンプルのペアエンドfastq、refディレクトリにリファレンスとアノテーション情報があるなら以下のように記載する。

#===========================================================

 // sequencing data

reads_R1="fastq/sample1_1.fq,fastq/sample2_1.fq"

reads_R2="fastq/sample1_2.fq,fastq/sample2_2.fq"

//The genome and annotation
annotation="ref/Ovis_aries.Oar_v3.1.90.chr14.gtf"
genome="ref/Ovis_aries.Oar_v3.1.dna.toplevel.chr14.fa"

//The genome and annotation of a related species
annotation_related_species="ref/Homo_sapiens.GRCh38.90.CDS.chr19.gtf"
genome_related_species="ref/Homo_sapiens.GRCh38.dna.toplevel.chr19.fa"

#===========================================================

 

 ジョブの最後には以下のように表示される。

||    Running time : 0.00 minutes                                             ||

||                                                                            ||

||                         Read assignment finished.                          ||

||                                                                            ||

|| Summary of counting results can be found in file "counts/block.counts"     ||

||                                                                            ||

===================== http://subread.sourceforge.net/ ======================//

 

======================================= Stage get_gene_info ========================================

 

======================================== Stage get_bp_info =========================================

 

====================================== Stage get_mapping_info ======================================

 

======================================== Pipeline Succeeded ========================================

04:46:37 MSG:  Finished at Tue Oct 30 04:46:37 UTC 2018

04:46:37 MSG:  Outputs are: 

relatives_superTranscriptome/genome_superT.relative.fasta

stats/mapping_info.txt

stats/size_info.txt

genome_superTranscriptome/genome_superT.fasta

de_novo_assembly/de_novo_assembly.fasta

... 30 more ...

kazu@7225a400f1a8:~/necklace$

 出力はsupertranscriptsのfasta、リードカウントファイルなど多岐に渡ります。詳細はwikiで確認して下さい。

https://github.com/Oshlack/necklace/wiki/Output

 

2021 3/9  出力例

f:id:kazumaxneo:20210309154552p:plain

f:id:kazumaxneo:20210309154554p:plain

 

 

引用

Necklace: combining reference and assembled transcriptomes for more comprehensive RNA-Seq analysis

Davidson NM, Oshlack A

Gigascience. 2018 May 1;7(5).

 

*1

samtoolsのバージョンが間違っていたため、初めエラーを起こした。

 

関連ツール

  


複数メタゲノムアセンブリのアセンブリ精度を比較して、種レベルでユニークな配列セットを得る dRep

2019 5/7 インストール追記、6/16 パラメータ追記、6/16  upしたdocker イメージのエラー修正、6/18 link追加

2021 4/29 インストール追記、5/18 インストール追記 (condaによるpplacerの導入)、5/27 タイトル変更、5/29, 6/30 compareコマンド追記

2022/06/14 ツイート追記

 

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

 

drep documentation

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を転載。

 

Documentより

Dereplicationとは、ゲノムセットの中から「同じ」ゲノム群を特定し、それぞれのセットから「最適」なゲノムを特定するプロセス。どの程度の類似性があれば「同一」とみなされるか、またどのように「最適」なゲノムを選択するかは、研究に応じて調整することができる。

 

2024/03/01

2022/06/14

 

インストール

依存

  • 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にして使う必要がある。

checkmは、インストール後、データベースをダウンロードしてパスを指定しておく必要がある(checkmの導入)。

 

#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
#パスを通しておく

#centrifuge
git clone https://github.com/infphilo/centrifuge
cd centrifuge
make -j 8
sudo make install prefix=/usr/local

#fastANI(binary)
wget https://github.com/ParBLiSS/FastANI/releases/download/v1.33/fastANI-Linux64-v1.33.zip
unzip
cp fastANI /usr/local/bin/

本体 Github

#Bioconda (link)
conda install -c bioconda drep

#pip
pip install drep

#pip & conda (2021 4/29)
mamba create -n drep -y python=3.7
conda activate drep
pip install drep
pip install checkm-genome

mamba install -c bioconda -y mash
mamba install -c bioconda -y mummer
mamba install -c bioconda -y fastANI
mamba install -c bioconda -y prodigal
mamba install -c bioconda -y pplacer

#インストール確認
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

 

                ...::: dRep v3.2.0 :::...

 

  Matt Olm. MIT License. Banfield Lab, UC Berkeley. 2017 (last updated 2020)

 

  See https://drep.readthedocs.io/en/latest/index.html for documentation

  Choose one of the operations below for more detailed help. 

  

  Example: dRep dereplicate -h

 

  Commands:

    compare            -> Compare and cluster a set of genomes

    dereplicate        -> De-replicate a set of genomes

    check_dependencies -> Check which dependencies are properly installed

    

> dRep compare -h

$ dRep compare -h

 

 

usage: dRep compare [-p PROCESSORS] [-d] [-h] [-g [GENOMES [GENOMES ...]]]

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

                    [-ms MASH_SKETCH] [--SkipMash] [--SkipSecondary]

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

                    [-nc COV_THRESH] [-cm {total,larger}]

                    [--clusterAlg {ward,average,complete,median,weighted,centroid,single}]

                    [--multiround_primary_clustering]

                    [--primary_chunksize PRIMARY_CHUNKSIZE]

                    [--greedy_secondary_clustering]

                    [--run_tertiary_clustering] [--warn_dist WARN_DIST]

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

                    work_directory

 

positional arguments:

  work_directory        Directory where data and output are stored    

                        *** 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 INPUT:

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

                        genomes to filter in .fasta format. Not necessary if

                        Bdb or Wdb already exist. Can also input a text file

                        with paths to genomes, which results in fewer OS

                        issues than wildcard expansion (default: None)

 

GENOME COMPARISON OPTIONS:

  --S_algorithm {ANImf,ANIn,fastANI,goANI,gANI}

                        Algorithm for secondary clustering comaprisons:

                        fastANI = Kmer-based approach; very fast

                        ANImf   = (DEFAULT) 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

                        goANI   = Open source version of gANI; requires nsmimscan

                         (default: ANImf)

  -ms MASH_SKETCH, --MASH_sketch MASH_SKETCH

                        MASH sketch size (default: 1000)

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

                        all genomes (default: False)

  --SkipSecondary       Skip secondary clustering, just perform MASH

                        clustering (default: False)

  --n_PRESET {normal,tight}

                        Presets to pass to nucmer

                        tight   = only align highly conserved regions

                        normal  = default ANIn parameters (default: normal)

 

GENOME CLUSTERING OPTIONS:

  -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)

  -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 and fastANI can only do larger method)

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

                        larger  = max*1

                         (default: larger)

  --clusterAlg {ward,average,complete,median,weighted,centroid,single}

                        Algorithm used to cluster genomes (passed to

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

 

GREEDY CLUSTERING OPTIONS

These decrease RAM use and runtime at the expense of a minor loss in accuracy.

Recommended when clustering 5000+ genomes:

  --multiround_primary_clustering

                        Cluster each primary clunk separately and merge at the

                        end with single linkage. Decreases RAM usage and

                        increases speed, and the cost of a minor loss in

                        precision and the inability to plot

                        primary_clustering_dendrograms. Especially helpful

                        when clustering 5000+ genomes. Will be done with

                        single linkage clustering (default: False)

  --primary_chunksize PRIMARY_CHUNKSIZE

                        Impacts multiround_primary_clustering. If you have

                        more than this many genomes, process them in chunks of

                        this size. (default: 5000)

  --greedy_secondary_clustering

                        Use a heuristic to avoid pair-wise comparisons when

                        doing secondary clustering. Will be done with single

                        linkage clustering. Only works for fastANI S_algorithm

                        option at the moment (default: False)

  --run_tertiary_clustering

                        Run an additional round of clustering on the final

                        genome set. This is especially useful when greedy

                        clustering is performed and/or to handle cases where

                        similar genomes end up in different primary clusters.

                        Only works with dereplicate, not compare. (default:

                        False)

 

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)

 

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

dRep dereplicate -h

 

$ dRep dereplicate -h

usage: dRep dereplicate [-p PROCESSORS] [-d] [-h] [-g [GENOMES [GENOMES ...]]]

                        [-l LENGTH] [-comp COMPLETENESS] [-con CONTAMINATION]

                        [--ignoreGenomeQuality] [--genomeInfo GENOMEINFO]

                        [--checkM_method {lineage_wf,taxonomy_wf}]

                        [--set_recursion SET_RECURSION]

                        [--checkm_group_size CHECKM_GROUP_SIZE]

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

                        [-ms MASH_SKETCH] [--SkipMash] [--SkipSecondary]

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

                        [-nc COV_THRESH] [-cm {total,larger}]

                        [--clusterAlg {median,centroid,ward,complete,single,average,weighted}]

                        [--multiround_primary_clustering]

                        [--primary_chunksize PRIMARY_CHUNKSIZE]

                        [--greedy_secondary_clustering]

                        [--run_tertiary_clustering]

                        [-comW COMPLETENESS_WEIGHT]

                        [-conW CONTAMINATION_WEIGHT]

                        [-strW STRAIN_HETEROGENEITY_WEIGHT] [-N50W N50_WEIGHT]

                        [-sizeW SIZE_WEIGHT] [-centW CENTRALITY_WEIGHT]

                        [-extraW EXTRA_WEIGHT_TABLE] [--warn_dist WARN_DIST]

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

                        work_directory

 

positional arguments:

  work_directory        Directory where data and output are stored    

                        *** 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 INPUT:

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

                        genomes to filter in .fasta format. Not necessary if

                        Bdb or Wdb already exist. Can also input a text file

                        with paths to genomes, which results in fewer OS

                        issues than wildcard expansion (default: None)

 

GENOME 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)

 

GENOME QUALITY ASSESSMENT OPTIONS:

  --ignoreGenomeQuality

                        Don't run checkM or do any quality filtering. NOT

                        RECOMMENDED! This is useful for use with

                        bacteriophages or eukaryotes or things where checkM

                        scoring does not work. Will only choose genomes based

                        on length and N50 (default: False)

  --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)

  --checkM_method {lineage_wf,taxonomy_wf}

                        Either lineage_wf (more accurate) or taxonomy_wf

                        (faster) (default: lineage_wf)

  --set_recursion SET_RECURSION

                        Increases the python recursion limit. NOT RECOMMENDED

                        unless checkM is crashing due to recursion issues.

                        Recommended to set to 2000 if needed, but setting this

                        could crash python (default: 0)

  --checkm_group_size CHECKM_GROUP_SIZE

                        The number of genomes passed to checkM at a time.

                        Increasing this increases RAM but makes checkM faster

                        (default: 2000)

 

GENOME COMPARISON OPTIONS:

  --S_algorithm {goANI,gANI,ANIn,fastANI,ANImf}

                        Algorithm for secondary clustering comaprisons:

                        fastANI = Kmer-based approach; very fast

                        ANImf   = (DEFAULT) 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

                        goANI   = Open source version of gANI; requires nsmimscan

                         (default: ANImf)

  -ms MASH_SKETCH, --MASH_sketch MASH_SKETCH

                        MASH sketch size (default: 1000)

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

                        all genomes (default: False)

  --SkipSecondary       Skip secondary clustering, just perform MASH

                        clustering (default: False)

  --n_PRESET {normal,tight}

                        Presets to pass to nucmer

                        tight   = only align highly conserved regions

                        normal  = default ANIn parameters (default: normal)

 

GENOME CLUSTERING OPTIONS:

  -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)

  -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 and fastANI can only do larger method)

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

                        larger  = max*2

                         (default: larger)

  --clusterAlg {median,centroid,ward,complete,single,average,weighted}

                        Algorithm used to cluster genomes (passed to

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

 

GREEDY CLUSTERING OPTIONS

These decrease RAM use and runtime at the expense of a minor loss in accuracy.

Recommended when clustering 5000+ genomes:

  --multiround_primary_clustering

                        Cluster each primary clunk separately and merge at the

                        end with single linkage. Decreases RAM usage and

                        increases speed, and the cost of a minor loss in

                        precision and the inability to plot

                        primary_clustering_dendrograms. Especially helpful

                        when clustering 5000+ genomes. Will be done with

                        single linkage clustering (default: False)

  --primary_chunksize PRIMARY_CHUNKSIZE

                        Impacts multiround_primary_clustering. If you have

                        more than this many genomes, process them in chunks of

                        this size. (default: 5000)

  --greedy_secondary_clustering

                        Use a heuristic to avoid pair-wise comparisons when

                        doing secondary clustering. Will be done with single

                        linkage clustering. Only works for fastANI S_algorithm

                        option at the moment (default: False)

  --run_tertiary_clustering

                        Run an additional round of clustering on the final

                        genome set. This is especially useful when greedy

                        clustering is performed and/or to handle cases where

                        similar genomes end up in different primary clusters.

                        Only works with dereplicate, not compare. (default:

                        False)

 

SCORING CRITERIA

Based off of the formula: 

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

 

A = completeness_weight; B = contamination_weight; C = strain_heterogeneity_weight; D = N50_weight; E = size_weight; F = cent_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)

  -centW CENTRALITY_WEIGHT, --centrality_weight CENTRALITY_WEIGHT

                        Weight of (centrality - S_ani) (default: 1)

  -extraW EXTRA_WEIGHT_TABLE, --extra_weight_table EXTRA_WEIGHT_TABLE

                        Path to a tab-separated file with two-columns, no

                        headers, listing genome and extra score to apply to

                        that genome (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)

 

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

dRep compare -h

$ dRep compare -h

usage: dRep compare [-p PROCESSORS] [-d] [-h] [-g [GENOMES [GENOMES ...]]]

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

                    [-ms MASH_SKETCH] [--SkipMash] [--SkipSecondary]

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

                    [-nc COV_THRESH] [-cm {total,larger}]

                    [--clusterAlg {ward,median,average,single,weighted,complete,centroid}]

                    [--multiround_primary_clustering]

                    [--primary_chunksize PRIMARY_CHUNKSIZE]

                    [--greedy_secondary_clustering]

                    [--run_tertiary_clustering] [--warn_dist WARN_DIST]

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

                    work_directory

 

positional arguments:

  work_directory        Directory where data and output are stored    

                        *** 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 INPUT:

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

                        genomes to filter in .fasta format. Not necessary if

                        Bdb or Wdb already exist. Can also input a text file

                        with paths to genomes, which results in fewer OS

                        issues than wildcard expansion (default: None)

 

GENOME COMPARISON OPTIONS:

  --S_algorithm {gANI,fastANI,ANImf,ANIn,goANI}

                        Algorithm for secondary clustering comaprisons:

                        fastANI = Kmer-based approach; very fast

                        ANImf   = (DEFAULT) 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

                        goANI   = Open source version of gANI; requires nsmimscan

                         (default: ANImf)

  -ms MASH_SKETCH, --MASH_sketch MASH_SKETCH

                        MASH sketch size (default: 1000)

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

                        all genomes (default: False)

  --SkipSecondary       Skip secondary clustering, just perform MASH

                        clustering (default: False)

  --n_PRESET {normal,tight}

                        Presets to pass to nucmer

                        tight   = only align highly conserved regions

                        normal  = default ANIn parameters (default: normal)

 

GENOME CLUSTERING OPTIONS:

  -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)

  -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 and fastANI can only do larger method)

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

                        larger  = max*3

                         (default: larger)

  --clusterAlg {ward,median,average,single,weighted,complete,centroid}

                        Algorithm used to cluster genomes (passed to

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

 

GREEDY CLUSTERING OPTIONS

These decrease RAM use and runtime at the expense of a minor loss in accuracy.

Recommended when clustering 5000+ genomes:

  --multiround_primary_clustering

                        Cluster each primary clunk separately and merge at the

                        end with single linkage. Decreases RAM usage and

                        increases speed, and the cost of a minor loss in

                        precision and the inability to plot

                        primary_clustering_dendrograms. Especially helpful

                        when clustering 5000+ genomes. Will be done with

                        single linkage clustering (default: False)

  --primary_chunksize PRIMARY_CHUNKSIZE

                        Impacts multiround_primary_clustering. If you have

                        more than this many genomes, process them in chunks of

                        this size. (default: 5000)

  --greedy_secondary_clustering

                        Use a heuristic to avoid pair-wise comparisons when

                        doing secondary clustering. Will be done with single

                        linkage clustering. Only works for fastANI S_algorithm

                        option at the moment (default: False)

  --run_tertiary_clustering

                        Run an additional round of clustering on the final

                        genome set. This is especially useful when greedy

                        clustering is performed and/or to handle cases where

                        similar genomes end up in different primary clusters.

                        Only works with dereplicate, not compare. (default:

                        False)

 

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)

 

Example: dRep compare 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

#2022/03/06 v3.0
docker pull kazumax/drep:3.0
source ${HOME}/.bashrc
dRep 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参照。

2021 5/18

1次クラスタリングでは90%でラフにクラスタリング(図の黒い破線)、2次クラスタリングではgANIを使って95%閾値でde-replication(mOTUsの取得)。5000bp以上の配列を対象とする。

dRep dereplicate outdir -g maxbin2/*fa \
-p 40 -l 5000 -pa 0.90 -sa 0.95 -comp 75 -con 25 -nc 0.1
  • -pa   ANI threshold to form primary (MASH) clusters  (default: 0.9)
  • -sa    ANI threshold to form secondary clusters (default: 0.99)
  • -nc    Minmum level of overlap between genomes when doing secondary comparisons (default: 0.1)
  • -l     Minimum genome length (default: 50000)
  • -comp   Minumum genome completeness (default: 75)
  • -con   Maximum genome contamination (default: 25)

 

compareコマンドを使ってbinned.fastaと既知ゲノム配列を比較する。

dRep compare outdir -g genome/*fna

 

メモ

  • 例えば、AがBに似ていて、BがCに似ていて、AとCが似ていないこのようなケースが存在すると仮定すると、 ANIがしきい値より大きいゲノムペアが異なるクラスタに入ってしまう可能性があるとされる(link)。

   => シングルモード(-clusterAlg single)で実行する

  • バクテリア以外の存在(ファージなど)が予想されて、checkMで評価できない場合、 -ignoreGenomeQualityというフラグを立てて、品質フィルタリングやゲノム選択時の完全性・汚染性の使用をオフにすることができる。
  • Mashは不完全なゲノム間の距離を過小評価し 、同一種のゲノムを複数のゲノムビンに分割してしまうことがあるため、第一段階のall versus all比較ではMashを使い、それから、精度が高い方法で2回目のクラスタリングを行う。

 

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

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

 

こちらのプレプリントにもパラメータ例があります。

https://www.biorxiv.org/content/10.1101/2021.04.02.438222v1.full.pdf

 

Paleofecesをターゲットにした 

Reconstruction of ancient microbial genomes from the human gut

https://www.nature.com/articles/s41586-021-03532-0

の論文でも使用されています。

引用

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


 2022/02/28

checkM2のランタイムを入れてもこちらの方がdRepよりずっと高速。ANI99以上にも対応している点も良い。


 

 

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

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

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

共発現遺伝子の自動クラスタリングツール CLust

 

 遺伝子の転写は、すべての生物において動的かつ協調的に調節される。そのような協調的調節は、時系列およびperturbation-responseデータセット (参考HP) における転写産物の協調的変化として明白である。遺伝子の転写は、遺伝子のプロモーター領域またはエンハンサー領域に位置するDNA /クロマチンelementsへの転写因子の結合によって調節される。典型的には、転写因子は、ゲノム中の全遺伝子数の約10%を構成し、転写の複雑な時空間パターンは、制御ネットワークにおけるこれらの遺伝子のcombinatorialなアクションを介して達成される[論文より ref.1]。これらのネットワークのcombinatorial な性質は、それらの動作が生来conditionalであることを意味する。すなわち、ある条件の下で同時発現しているように見える遺伝子は、必ずしもすべての条件下で同時発現されるとは限らない。この結果は、いずれかの実験的状況(例えば、いくつかの生物学的プロセスにまたがる時系列実験またはperturbation-response実験)において、すべての遺伝子が協調的に行動するわけではないということである。代わりに、遺伝子のサブセットは、実験的な状況の中では特定のレギュレーターの組み合わせにより協調的に行動するが、他の遺伝子は、実験のデザインとは独立したパターンに従う。したがって、与えられたobservation window(すなわち実験内容)内では、すべての遺伝子が限られた協調的振る舞いにアサインされることは期待されない[ref.2,3, link2,  pubmed3]。

 遺伝子のサブセットのみが特定の状況において同時発現される可能性が高いことを考えると、これらのサブセットの同定はデータ抽出問題であり、データ分割問題ではない。すなわち、遺伝子発現クラスタリングの目的は、特定の状況において検出される完全な遺伝子セットから協調的に行動する遺伝子のコホート(集団)を同定および抽出することであり、全遺伝子セットを分類することではない。実際には、クラスタリングは、それらがコーディネートされた挙動(すなわち、同時発現遺伝子のクラスター)を有する遺伝子の個別コホートの完全なセットを同定することを期待して遺伝子発現データに広く適用されており、それらの行動を示す遺伝子は、正しいクラスターに割り当てられる[ref.4,5]。しかし、共発現遺伝子のコホートを同定することを目的とする大部分の方法は、データ分割(例えば、k-means [ref.6]、階層的クラスタリング[ref.7]、自己組織化マップ[ref.8])(参考HP)に基づいている。これらのアプローチは、データ分割メトリックの数値最適化によって決定されるクラスタの数を用いて、すべての遺伝子を有限集合のクラスタに割り当てようと試みる[ref.9]。従って、調査中の文脈において同時発現されていない遺伝子もまた、それらの「最良適合」クラスターに割り当てられ、その結果、クラスターの大部分は、共発現遺伝子および非共発現遺伝子の両方を含むであろう。この結果は、共発現された遺伝子クラスターの生物学的特性の期待に従わない。すなわち、各クラスターは、問題の実験的または生物学的状況において強調した挙動を示す遺伝子のみを含み、どのクラスターも同一のプロファイルを持たない。データ分割メソッドが最も一般的に使用されており、多くの部分クラスタリング方法も提案されており [ref.5,10,11,12]、それの方法では、完全なデータセットクラスタ間で分割される必要はない。代わりに、クラスタに簡単に割り当てることができるサブセットを特定することを目指している。

 ここでは、5つのモデル生物から100の実際の生物学的データセットの分析を通して、データ分割ベースおよび部分クラスタリングベースの方法を遺伝子発現データに適用すると、アサインに信頼性のない相当数の遺伝子が出ることを示す、すなわち排他的に適合しない遺伝子で、そのクラスターでは除外されていたはずの遺伝子である。そのような信頼できないコンテンツは、これらのクラスタの約50%までを構成する。この問題を解決するために、遺伝子発現データからクラスタ抽出のためのclustと呼ばれる新しい方法を提供する。 Clustは、同時発現遺伝子クラスターの生物学的期待を満たす遺伝子の共発現クラスターを抽出するように設計されている。Clustは、データ分割方法および部分的クラスタリング方法よりも低い分散レベルで共発現されたクラスターを抽出するこことを示す。著者らはまた、7つの評価metricsを使い、clustによって生成されるクラスタが他のテストされたどのツールより良好であることを示す。さらに、clustによって抽出されたクラスターのfunctional termsが、他の方法によって生成されたものよりも同等に、または大幅に豊富であることを示す。最後に、clustを使い複数のデータセットで一貫して同時発現する遺伝子のクラスタを抽出し、このclustを活用することで、研究者が複数のデータセットを活用して、高精度に co-expressedの遺伝子クラスターを同定可能になることを示す。

 

f:id:kazumaxneo:20181027194058p:plain

 

Githubより

 

Features (Githubより)

  1. No need to pre-process your data; clust automatically normalises the data.
  2. No need to preset the number of clusters; clust finds this number automatically.
  3. You can control the tightness of the clusters by varying a single parameter -t
  4. It is okay if the datasets:
  • Were generated by different technologies (e.g. RNA-seq or microarrays)
  • Are from different species
  • Have different numbers of conditions or time points
  • Have multiple replicates for the same condition
  • Require different types of normalisation
  • Were generated in different years and laboratories
  • Have some missing values
  • Do not include every single gene in every single dataset

 

インストール

ubuntu16.04のminiconda2-4.0.5環境にてテストした(docker使用。ホストOS: mac os10.12)。

依存

Clust is a Python package that requires Python 2.7 but is not compatible with Python 3 yet.

  • numpy
  • scipy
  • matplotlib
  • sklearn
  • sompy
  • joblib
  • portalocker

本体 Github

> clust

$ clust  

usage: clust [-h] [-n <file or int> [<file or int> ...]] [-r <file>]

             [-m <file>] [-o <directory>] [-t <real number>]

             [-basemethods <string> [<string> ...]]

             [-K <integer> [<integer> ...]] [-s <real number>] [-d <integer>]

             [-fil-v <real number>] [-fil-c <integer>] [-fil-d <integer>]

             [--fil-abs] [--fil-perc] [--fil-flat] [--no-fil-flat]

             [-cs <integer>] [-q3s <real number>] [--no-optimisation]

             [--deterministic] [-np <integer>]

             datapath

 

/===========================================================================\

|                                   Clust                                   |

|     Optimised consensus clustering of multiple heterogeneous datasets     |

|                              Version v1.8.4                               |

|                                                                           |

|                            By Basel Abu-Jamous                            |

|                       Department of Plant Sciences                        |

|                         The University of Oxford                          |

|                      basel.abujamous@plants.ox.ac.uk                      |

+---------------------------------------------------------------------------+

|                                 Citation                                  |

|                                 ~~~~~~~~                                  |

| When publishing work that uses Clust, please cite:                        |

| Basel Abu-Jamous and Steven Kelly (2018) Clust: automatic extraction of   |

| optimal co-expressed gene clusters from gene expression data. Genome      |

| Biology 19:172; doi: https://doi.org/10.1186/s13059-018-1536-8.           |

+---------------------------------------------------------------------------+

| Full description of usage can be found at:                                |

| https://github.com/BaselAbujamous/clust                                   |

\===========================================================================/

 

positional arguments:

  datapath              The directory that includes the data files.

 

optional arguments:

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

  -n <file or int> [<file or int> ...]

                        Normalisation file or list of codes (default: 1000)

  -r <file>             Replicates structure file

  -m <file>             OrthoGroups (OGs) mapping file

  -o <directory>        Output directory

  -t <real number>      Cluster tightness (default: 1.0).

  -basemethods <string> [<string> ...]

                        One or more base clustering methods (default: k-means)

  -K <integer> [<integer> ...]

                        K values, e.g. 2 4 6 10 ... (default: 4 to 20 (step=4))

  -s <real number>      Outlier standard deviations (default: 3.0)

  -d <integer>          Min datasets in which a gene must exist (default: 1)

  -fil-v <real number>  Filtering: gene expression threshold (default: -inf)

  -fil-c <integer>      Filtering: number of conditions (default: 0)

  -fil-d <integer>      Filtering: number of datasets (default: 0)

  --fil-abs             Filter using absolute values of expression

  --fil-perc            -fil-v is a percentile of genes rather than raw value

  --fil-flat            Filter out genes with flat expression profiles (default)

  --no-fil-flat         Cancels the default --fil-flat option

  -cs <integer>         Smallest cluster size (default: 11)

  -q3s <real number>    Q3s defining outliers (default: 2.0)

  --no-optimisation     Skip cluster optimsation & completion

  --deterministic       Obsolete as all steps are already deterministic (v1.7.4+)

  -np <integer>         Number of parallel processes (default: 1)

gene expressionのファイルをディレクトリに用意すればランできる。書式は、1列目にgene ID、1行目にユニークな名前(下の画像 (Githubより) )。タブ区切りにする。

f:id:kazumaxneo:20181027204431p:plain

#実行方法 - ファイルをinput/にいれているなら
clust input/

 

 

テストラン 

ここでは3つの異なる生物種から replicatesデータ数も異なるデータセット画像)のクラスタリングを実行する(three datasets from two yeast species, two datasets from fission yeast, and one from budding yeast)。

Git clone https://github.com/BaselAbujamous/clust.git
cd clust/ExampleData/1_RawData/

#Data/のX1.txt、X2.txt、X3.txtを入力としてクラスタリングを実行
clust Data/ -r Replicates.txt -n Normalisation.txt -m MapIDs.txt

1データセットに複数のreplicatesがある時は、 画像(リンク)のようなファイルを準備して、上のテストのように"-r replicates_file"で指定する。複数の生物種のデータセットからのクラスタリング時は"-m MapIDs.txt"で指定する(Data from multiple species)。複数のトランスクリプトームテクニックからのクラスタリングData from multiple technologies (e.g. mixing RNA-seq and microarrays)にも対応している。

解析

/===========================================================================\

|                              RESULTS SUMMARY                              |

+---------------------------------------------------------------------------+

| Clust received 3 datasets with 9332 unique genes. After filtering, 9330   |

| genes made it to the clustering step. Clust generated 14 clusters of      |

| genes, which in total include 6124 genes. The smallest cluster includes   |

| 122 genes, the largest cluster includes 1006 genes, and the average       |

| cluster size is 437.428571429 genes.                                      |

+---------------------------------------------------------------------------+

|                                 Citation                                  |

|                                 ~~~~~~~~                                  |

| When publishing work that uses Clust, please include this citation:       |

| Basel Abu-Jamous and Steven Kelly (2018) Clust: automatic extraction of   |

| optimal co-expressed gene clusters from gene expression data. Genome      |

| Biology 19:172; doi: https://doi.org/10.1186/s13059-018-1536-8.           |

+---------------------------------------------------------------------------+

| For enquiries contact:                                                    |

| Basel Abu-Jamous                                                          |

| Department of Plant Sciences, University of Oxford                        |

| basel.abujamous@plants.ox.ac.uk                                           |

| baselabujamous@gmail.com                                                  |

\===========================================================================/

kazu@becc7e791608:~/clust/ExampleData/1_RawData$ l

total 68K

正規化方法は自動で算出される(マニュアル指定は-nを使う)。数値それぞれ異なる意味がある。詳細はGithub参照

出力

f:id:kazumaxneo:20181027181900p:plain

Clust generates the following output files (Githubより)

  • A table of clustering statistics
  • A table listing genes included in each cluster
  • Pre-processed (normalised, summarised, and filtered) datasets' files
  • Plotted gene expression profiles of clusters (a PDF file)

> cat /Results_27_Oct_18/Normalisation_actual #正規化ファイル

$ cat /Results_27_Oct_18/Normalisation_actual 

X1.txt 101 4

X2.txt 101 4

X3.txt 101 3 4

>  cat Results_27_Oct_18/Summary.tsv 

$ cat Results_27_Oct_18/Summary.tsv 

Starting data and time Saturday 27 October 2018 (09:14:00)

Ending date and time Saturday 27 October 2018 (09:15:34)

Time consumed 0 hours, 1 minute, and 34 seconds

Number of datasets 3

Total number of input genes 9332

Genes included in the analysis 9330

Genes filtered out from the analysis 2

Number of clusters 14

Total number of genes in clusters 6124

Genes not included in any cluster 3206

Min cluster size 122

Max cluster size 1006

Average cluster size 437.428571429

Total number of samples (including replicates) 38

Average number of samples (including replicates) per dataset 12.6666666667

Total number of conditions 38

Average number of conditions per dataset 12.6666666667

**********************************

Genes included in all datasets 0

Genes missed from 1 dataset 4911

Genes missed from 2 datasets 4421

Genes missed from 3 datasets 0

Genes missed from more than 3 datasets 0

Clusters_Objects.tsv

f:id:kazumaxneo:20181027200217p:plain

クラスターの厳密さは-t で指定する。デフォルトは1.0(Are you obtaining noisy clusters?)。

 

すでにノーマライズ済みのデータも使うことができる。

cd ExampleData/2_Preprocessed/

#このデータにはreplicatesがない。ノーマライズの必要もないので、関連オプションは全くなしでランする。
clust Data/ -t 5 -o MyResultsDirectory/
  • -o   Output directory
  • -t    Cluster tightness (default: 1.0).

 

引用

Clust: automatic extraction of optimal co-expressed gene clusters from gene expression data

Basel Abu-Jamous, Steven Kelly

Genome Biology, 201819:172

 

ハイブリッドアセンブリのためのアライメントフリー scaffolding graph構築ツール Fast-SG

2018 10/26 タイトル修正 

 

 ゲノム全体のデノボアセンブリの主要な課題は、リピートを解決することである[論文より 1,2]。リピートは、ゲノムの複数の位置で生じるほぼ同一のゲノム配列に対応する。この課題に対処するために、主に2つのタイプのアプローチが提案されている。一つはのショート・リードのペアの情報[ref.3]を使ったアプローチと、それ以外のロング・リード[ref.4]を使ったアプローチである。後者は、ロングリード中にリピートを完全にキャプチャすることが目的である。そのようなロングリードの非反復サフィックスプレフィックスシーケンスを使用してユニークなオーバーラップを計算すると、元のリードをコンティグと呼ばれる大きな配列に曖昧さなく拡張でき、時には(必ずしもそうではないが)完全なゲノム配列を直接推測できる。前者のアプローチは、genome scaffoldingと呼ばれる操作に関連している必要がある。ショートリードは、オーバーラップの計算[rerf.5]またはde Bruijn graph[ref.6]のいずれかによって、上記のように最初にコンティグに組み立てられる。しかし、この場合得られたコンティグは全ゲノムにまたがっていない。ほとんどの場合、ずっと短くなる。それらは、第2のステップで結合される(すなわち、一緒にリンクされる)必要がある。リンキング情報は、一般に、ペアエンドまたはメイトペアシーケンスによって提供される。一般的に、両端がシーケンスされている1kbより大きいゲノムフラグメントはメイトペアライブラリと呼ばれ、そうでなければペアエンドライブラリと呼ばれる。ペアリードを使用するgenome scaffoldingは、コンティグ間にギャップ(すなわち、未知の配列)を導入し、それによってゲノム配列全体には至らず、いわゆるscaffold配列のセットにつながる。したがって、scaffoldは、順序と方向が決まったコンティグのセットを表す。

 Genome scaffoldingの問題は、最初にHusonらによって定式化された[ref.7]。著者らが提案した方法は、ノードがコンティグを表現し、エッジがメイトペア(重み)の数、方向、および2つの異なるコンティグ間の距離を符号化する、scaffold graphと呼ばれるものを構築することから始まった。greedyアルゴリズムを使用して、ヒューリスティックにscaffold配列に対応する最適な経路を得る。

 Husonらの定式化から開発されたスキャフォールディング法の大部分は、同種のグラフを使用しており、超高速のショートリードアライナー[ref.8-10]でscaffoldの基礎を構築している[ref.3]。この分野におけるアルゴリズムの革新は、主に最適経路(通常は最大の重み)を選択する方法に重点を置いており、したがって、精度の高いlarge scaffoldsを得ることができる。ダイナミックプログラミング[ref.11]、breadth-first search [ref.12]、maximum weight matching[ref.13]、branch and bound[14]などの様々なアプローチが提案されている。

 新しいロングリードシークエンシング技術(Pacific Biosciences、Oxford Nanopore)は、高レベルのエラー(現在の平均15%)を含む非常に長いリード(> 10kb)を生成することによって、急激にゲノムアセンブリを変更した。これらの新技術は解決可能なリピート配列の景観を拡大した[ref.15]。現在、このようなロングリード[ref.4,16]を使用するde novoアセンブラは、バクテリアゲノムをFinishさせ、ヒトゲノムの高度に連続した再構成を生成することができる[ref.4,17]。しかし、オーバーラップを計算すること[ref.5]に基づく大規模なゲノムのデノボアセンブリは計算的に激しく[ref.4]、エラーを含むロングリードシーケンスのエラー訂正のためにかなりのカバレッジ(50X)を必要とし、これらの方法の大きなゲノムのデノボアセンブリへの幅広い適用を妨げている[17]。

 しかしながら、新規ライブラリ作成技術の補完的な長距離情報と関連させると、ロングライブラリを用いたデノボアセンブリは染色体レベルにスケーラビリティがあることが証明されている[18,19] [20,21]。このような新規の実験ライブラリは、イルミナ(Illumina)のマシンでシーケンスされ、従来のペアエンドリードに至る。したがってDOVETAILゲノミクス[20](参考ページ)は1-200kbの範囲で有用なリンク情報を生成し、一方10Xゲノミクス[ref.22]は、バーコードを賢く活用することで、最大100kbのリンクリードを生成する。どちらの技術も、アセンブリパイプライン内で長距離情報を使用してscaffolding graphを構築し、scaffolding配列を得るための独自のアルゴリズム解を適用する。両方の技術は、高価で時間のかかる長距離のメイトペアライブラリを作製するために必要な実験プロトコルを、ショートリードシーケンスで置き換えることを目指して考案された[23,24]。

 原則的に、ロングレンジの情報は、後者の、ロングリードの実際のサイズの制限範囲内で直接抽出することができる。このような情報はハイブリッドアセンブリ法を考案するために使用することができ、ショートリードアセンブリから高品質のコンティグがscaffolding graphのノードとして使用され、エッジがロングリードからのリンク情報を使用して作成され、scaffoldはショートリードscaffolderで形成される。しかし、現在、ロングリードからscaffolding graphを構築するためのアルゴリズムが欠けている。そのようなアルゴリズムは、効率的な既存のショートリードアルゴリズムを新しいハイブリッドアセンブリパイプラインを構成するのに再利用することを可能にする。

 既存の効率的なショートリードスキャナーと互換性がある構造標準を維持しながら、適度な計算資源を用いて、ショートリードまたはロングリードから超短時間でこのようなgraphを構築できることは、ここで取り上げる主な課題である。我々(著者ら)が提案しているFAST-SGは、さまざまな配列ソース(イルミナ、パシフィックバイオサイエンス、オックスフォードナノポア)からの情報だけでなく、アライメントフリーのアルゴリズム[ref.25]ストラテジーを使用し、スケーラビリティ、スピード、およびモジュール性を最大化するように考えられた。後者の特徴は、特に、大きなゲノムの効率的なアセンブリを可能にする新規のハイブリッドアセンブリパイプラインを定義することを可能にする。

 FAST-SGは包括的な一連の標準データセット[ref.3、26]とベンチマークを使用して広範にテストされた。我々(著者ら)は、FAST-SGが大きなゲノムのハイブリッドアセンブリを可能にし、特にカバレッジの浅いロングリードカバレッジデータ(5X〜10X)で有効であることを示す。このハイブリッド戦略は、バクテリア人工染色体(BACs、180kb)までのインサートサイズを有することができるいくつかの合成メイトペアライブラリーの構築から成り、長距離のscaffoldを生成するためにショートリードscaffolderと組み合わせることができる。このような戦略は、適度な計算資源でヒトゲノムサイズまでスケーリングできる。さらに、FAST-SGは従来のショート・リード・アライナーよりも高速(7X-15X)で、ショートリードのメイトペアデータを使用したスキャフォールディングの強力な代替手段であることを示す。

 我々(著者ら)は、FAST-SGの効果的なハイブリッドアセンブリのための手順を提供することによって結論を出し、提案する戦略がどのようにロングリードを使用してギャップを埋めるように拡張され、そしてscaffoldsのエラーを訂正するか議論する。

 

Fast-SGに関すつツイート

 

インストール

mac os10.12でテストした。

依存

  •  KMC

本体 Github

git clone https://github.com/adigenova/fastsg
cd fastsg/
make all

#KMCのダウンロード(pre-compile) ./KMC/bin/にないと認識しないので解凍後、KMC/bin/にKMCバイナリを放り込む
#mac
wget https://github.com/refresh-bio/KMC/releases/download/v3.0.0/KMC3.mac.tar.gz
tar -zxvf KMC3.mac.tar.gz
mv KMC3.mac KMC
cd KMC
mkdir bin && mv kmc* bin/ && cd ../

make test

 test結果

ont.I1000.FastSG_K15.sam

ont.I2000.FastSG_K15.sam

ont.I3000.FastSG_K15.sam

End of mapping in thread

Query of 56850535 kmer-pairs  in 34.01s rate 1671387.23 per second

Query of 286428 pairs  in 34.01s

Number of total mapped kmer-pairs 4751139

Number of total mapped reads-pairs 144987

Number of total long reads 1298

Number of total mapped pairs: 144988

touch test.fast-sg.K15.done

make[1]: Leaving directory '/Volumes/4TB-gene_new2/Fast-SG/fastsg'

 

 

ENDING Fri Oct 26 20:06:55 2018 FAST-SG with K=15

 

all done

正常に動いている。手動で行うなら

cd fastsg/ #fastsgのrootに移動
perl FAST-SG.pl -k 15 -l example/ecoli-reads.txt -r example/ecoli-illumina.fa.gz -p test -t 20

ヘルプ

perl FAST-SG.pl   #"FastSG++"をperlのラッパーで呼び出している。

$ perl FAST-SG.pl 

 

 

ERROR: Mandatory options need to be specified -k -l -r -p

 

 

 

Usage: FAST-SG.pl  [<arguments>]

 

 

Basic options:

 

-k *K-mer size [12-256] could be a single value or a range [15-40:5]

-l *Read configuration file

-r *Contig fasta file

-p *Output prefix

-t Number of thread to be used [1]

 

 

Advanced options:

 

Illumina reads options:

-x Vector hit size [10]

-y Mimimum score to report reads [6]

-z Mimimum score for pair rescue [4]

 

The value of -x must be larger than -y and -y larger than -z [-x > -y > -z]. 

To modify these values take into account the k-mer size, the length and error rate of the reads.

Specified values must be larger than defaults values.

 

Long reads options:

-u Vector hit size [20]

-v Mimimum score to report a pair of reads [15]

-w Mimimum score for pair rescue [4]

-s Length of syntethic read [200]

 

The value of -u must be larger than -v and -v larger than -w [-u > -v > -w]. 

To modify these values take into account the k-mer size, the length and error rate of the reads.

Specified values must be larger than defaults values.

 

Other options:

-c Full path to KMC3 directory [FAST-SG-DIR/KMC]

-N Show pipeline commands [-N 1]

-R Delete results to restart FAST-SG [-R 1]

--help Show help message

--version Display current version of FAST-SG

 

Note: Options with description starting with (*) are mandatory.

 

 

テストラン

NA12878のハイブリッドアセンブリを例にフローが説明されている。この流れを確認する。

https://github.com/adigenova/fast-sg/wiki/Hybrid-scaffolding-of-NA12878

 

1、illuminaのショートリードのde novoアセンブリ結果をダウンロード。

cd fastsg/ #fastsgのrootに移動

#DiscoverDENOVOによるilluminaショートリードのアセンブリ結果をダウンロードする
wget ftp://ftp.broadinstitute.org/pub/crd/Discovar/assemblies/51400.newchem/a.lines.fasta -O NA12878.asm.fa

#fa.faiを利用して2kb以上のcontigだけ取り出す。
samtools faidx NA12878.asm.fa
sort -rn -k2,2 NA12878.asm.fa.fai | awk '{if($2 >=2000){print $1}}' > NA12878.asm.2kb.lst
samtools faidx NA12878.asm.fa $(cat NA12878.asm.2kb.lst | xargs) > NA12878.asm.2kb.fa

 

2、ONTのロングリードをダウンロードし、コンカテネートしてFASTAに変換。このリンク先のデータ(link)が使われている(preprint)。

cat > downlod.txt <<EOF
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-16056159-FAF15665.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-17958431-FAF13748.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-2901545329-FAF10039.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-3439856925-FAF09968.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-3709819546-FAF09277.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-3976726082-FAF14035.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4109802543-FAF15694.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4111860526-FAF09713.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4178920553-FAF18554.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4244782843-FAF15630.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4245291640-FAF09640.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4249180049-FAF09701.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-82266371-FAF15586.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-87644245-FAF05869.fastq.gz
EOF

#wgetにテキストで与えダウンロード。
mkdir fastq && cd fastq
wget -i downlod.txt

#concatenate
zcat fastq/*.fastq.gz | awk 'BEGIN{h=1;s=2;i=1}{if(NR==h){print ">ONT_"i; h+=4;i++;} if(NR == s){print $0; s+=4;}}' | fold | gzip > ULTRA-LONG-RENAME-FOLD.fa.gz

 

3、configファイルを準備する。Fast-SGを使い、ONTのロングリードを使ってilluminaショートリードアセンブリをscaffoldingして、2kbから180kbまでのsynthetic librariesを作成する場合、以下の書式になる。

cat > ultra-long-conf.txt <<EOF
long ont ULTRA-LONG-RENAME-FOLD.fa.gz 2000,4000,6000,8000,10000,12000,14000,16000,18000,20000,30000,40000,50000,60000,70000,80000,100000,120000,150000,180000 1
EOF

以下のような書式になっている。

 

#shortリードの場合
#type libID Path(fwd) Path(rev) SAM(1:single 2:paired)
short lib1 example/ecoli.ill-sim.fwd.fq.gz example/ecoli.ill-sim.rev.fq.gz 1

#longリードの場合
#type libID Path(long read) Insert-sizes SAM(1:single 2:paired)
long pac example/ECOLI-PACSEQ.subset.fasta.gz 1000,2000,3000,5000 1
long ont example/ECOLI-ONT-1D.subset.fasta.gz 1000,2000,3000,5000 1

 example  https://github.com/adigenova/fast-sg/blob/master/example/ecoli-reads.txt

 

4、実行

perl FAST-SG.pl -k 22 -l ultra-long-conf.txt –r NA12878.asm.2kb.fa.gz -p NA12878-ultra –t 20 

 -k、-l、-r、-pはx必須。複数k-merを指定する際は"-k 15-25:5"のように指定する。エラー訂正されていないロングリードを使う場合は、15-22程度の短いk-merを使うことが推奨されている。

ジョブが終わると、insert librariesで指定したサイズのsamファイルができる。

*私の環境では4日経ってもジョブが終わらなかったので、一旦テストを終了します。  

 

5、samをforwardとreverseに分割する。

awk -f fastSG2scaff.awk -v name=$(basename ultra_ont_raw.I2000.fast-sg_K22.sam) -v k=22 ultra_ont_raw.I2000.fast-sg_K22.sam

 

6、ScaffMatch(リンク)を使ってScaffoldingを実行する。

scaffmatch -m -w UONT-K22-7X -c NA12878.asm.2kb.fa \
-s 201,416,632,847,1063,1279,1495,1710,1925,2141,3216,4291,5364,6434,7495,8555,10675,12781,15946,19133 \
-i 2010,4166,6323,8478,10636,12793,14950,17106,19259,21414,32169,42912,53641,64340,74955,85557,106753,127813,159463,191339 \
-p fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr \
-1 ultra_ont_raw.I2000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I4000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I6000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I8000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I10000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I12000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I14000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I16000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I18000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I20000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I30000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I40000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I50000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I60000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I70000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I80000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I100000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I120000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I150000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I180000.Fast-SG_K22.fwd.sam \
-2 ultra_ont_raw.I2000.Fast-SG_K22.rev.sam,ultra_ont_raw.I4000.Fast-SG_K22.rev.sam,ultra_ont_raw.I6000.Fast-SG_K22.rev.sam,ultra_ont_raw.I8000.Fast-SG_K22.rev.sam,ultra_ont_raw.I10000.Fast-SG_K22.rev.sam,ultra_ont_raw.I12000.Fast-SG_K22.rev.sam,ultra_ont_raw.I14000.Fast-SG_K22.rev.sam,ultra_ont_raw.I16000.Fast-SG_K22.rev.sam,ultra_ont_raw.I18000.Fast-SG_K22.rev.sam,ultra_ont_raw.I20000.Fast-SG_K22.rev.sam,ultra_ont_raw.I30000.Fast-SG_K22.rev.sam,ultra_ont_raw.I40000.Fast-SG_K22.rev.sam,ultra_ont_raw.I50000.Fast-SG_K22.rev.sam,ultra_ont_raw.I60000.Fast-SG_K22.rev.sam,ultra_ont_raw.I70000.Fast-SG_K22.rev.sam,ultra_ont_raw.I80000.Fast-SG_K22.rev.sam,ultra_ont_raw.I100000.Fast-SG_K22.rev.sam,ultra_ont_raw.I120000.Fast-SG_K22.rev.sam,ultra_ont_raw.I150000.Fast-SG_K22.rev.sam,ultra_ont_raw.I180000.Fast-SG_K22.rev.sam

 UONT-K22-7X/Scaffolds.faができる。

テスト中 

 

引用

Fast-SG: an alignment-free algorithm for hybrid assembly

Di Genova A, Ruz GA, Sagot MF, Maass 

Gigascience. 2018 May 1;7(5).

 

複数のSVコーラーを動かし、結果を統合する Parliament2

2018 10/26 エラー修正

2018 10/28 エラー修正

2019 3/2 追記

2019 6/11 twitter追記

2019 7/1 dockerインストールをlatestタグに修正=>エラーがあったため0.1.17に戻した

2019 12/18 説明微修正

2020 12/22 論文引用

 

 構造変異(SV)は、ゲノムの大きな(50bp+)変異である[論文より ref.1,2]。これらのバリエーションは、個々のショート・リードのサイズに近いか、最も頻繁には、リードより大きいため、ショート・リード・データで直接検出することは難しく、スプリット・リード・マッピング、ソフト・クリッピング、リードペアの距離と方向の変化、カバレッジデプス変化、またはヘテロ接合性の変化[ref.3]によって推測される。 SNPおよびIndelの検出精度は99.9%を超えるが、、Genome in a Bottle [ref.5] (GIAB HP)の最近のベンチマークセットによると、SVパイプラインの精度は著しく低い。

 Breakdancer [ref.6]、CNVnator [ref.7]、Crest [ref.8]、Delly [ref.9] Lumpy [ref.10]、Manta [ref.11]、Pindel [ref.12]などのさまざまなSVコーラーは、様々なheuristicsを利用して、マップされたリードの信号を読み取り構造変化を推測する。これらツールのアプローチの多様性により、SVの検出感度がSVタイプのバリエーション(たとえば、 DEL(DEL)、INS(INS)、DUPDUP)、INV(INVA)、TRA(Translocation)やサイズにより変化する。
 アプローチとドメインパフォーマンスのこの多様性は、アンサンブル戦略を魅力的にする。 2つの方法:MetaSV [ref.13](紹介)とParliament1 [ref.14]は、3段階オーバーラップ - マージ -戦略を使用して複数のコーラー結果から高品質のコンセンサスセットに結合する。 MetaSVとParliament1の両方は、検証ステップにアセンブリベースの方法を使用している。これはParliament1では非常に正確だが、計算量が多く、イベントの最大サイズが制限されている。 MetaSVとParliament1は、すでに生成されたSVバリアントコールから開始するため、個々のSVコーラーをインストールして実行する負担をユーザに負わせている。
 ここでは、調整可能な感度と特異性を維持し1万人のWGSコホートスケールにSVコールを経済的に拡張するために設計されたParliament2(Parliament1の改良版)を提示する。 Parliament2を使用すると、Breakdancer、Breakseq、CNVnator、Delly、Lumpy、Manta(SV検出ツール紹介Manta紹介)の一部またはすべてを実行してコールを生成できる。 SURVIVOR [ref.15](紹介)を使用して、これらのオーバーラップがあるコーからコンセンサスコールを生成する。 SVTyper [ref.16](紹介)を使用してこれらのコールを検証する。
 Parliament2は、CPU、ディスクI / O、および各プログラムが所与の時間に持つRAMのさまざまな要件を利用して、コール元を同時に実行するためのさまざまな並列化戦略を使用する。これによりParliament2は、最も遅いコーラーであるCNVnatorとほぼ同じウォールクロックとCPU時間ですべてのコーラーの解析を完了させることができる。 16コアマシンならば、2〜5時間で35X Whole Genome Sequence(WGS)サンプルを処理できる。 Parliament2は、個々のコーラーのraw resultと、サポートされているすべてのコーラーの遺伝子型コンセンサスVCFの両方を生成する。
Parliament2はオープンソースであり、コードベースとしても、またはすべてのコーラーを単一、または複数で簡単に実行できるDockerイメージとしても利用できる。Parliament2は、DNAnexus上で公開されているアプリケーションである。これにより、必要に応じて個々のコンポーネントだけを実行することもできる。さらに、Parliament2は容易に拡張可能で、ユーザーはツールを実行してParliament2の高速実行を開始できる。オプションで、SVVIZ [ref.17](紹介)を使用してSVのPDFイメージを提供することもできる。


2018年10月現在、Preprintが発表されています。first authorとsecond authorはDNAnexus社 (wiki) の研究者の方です(*6)。

 

アナウンス

https://blog.dnanexus.com/2018-09-24-parliament2-fast-structural-variant-calling/

 

f:id:kazumaxneo:20181025203908j:plain

Parliament2のワークフロー。各ツールを並列で走らせた後、SURVIVORで統合し、SVtyperでgenotypingを行う(上は順番が逆?)。最後にsvvizを使ってSVイベント領域のリードを可視化する。アナウンスのHPより転載。

 

Parliament2に関するツイート

 

2019 6/11 twitter追記

 

インストール

dockerのubuntu16.04にてテストした(ホストOSはmac os10.12)。

本体 Github

dockerイメージをpullして使う(buildするとエラーが出た)。

git clone https://github.com/dnanexus/parliament2.git
cd parliament2/
dokcer build -t parliament2 .

#またはpullする(tag link)
#0.1.17
docker pull dnanexus/parliament2:0.1.7

#latest
docker pull dnanexus/parliament2:latest

> docker run -it dnanexus/parliament2:0.1.7 -h

$ docker run -it dnanexus/parliament2:0.1.8 -h

flag needs an argument: 'h' in -h

See 'docker run --help'.

——

user~/Desktop/6803T1_all_tools_referebce_GT-S: $ docker run -it ccb464d859a8 -h

usage: parliament2.py [-h] --bam BAM [--bai BAI] -r REF_GENOME [--fai FAI]

                      [--prefix PREFIX] [--filter_short_contigs]

                      [--breakdancer] [--breakseq] [--manta] [--cnvnator]

                      [--lumpy] [--delly_deletion] [--delly_insertion]

                      [--delly_inversion] [--delly_duplication] [--genotype]

                      [--svviz] [--svviz_only_validated_candidates]

 

Parliament2

 

optional arguments:

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

  --bam BAM             The name of the Illumina BAM file for which to call

                        structural variants containing mapped reads.

  --bai BAI             (Optional) The name of the corresponding index for the

                        Illumina BAM file.

  -r REF_GENOME, --ref_genome REF_GENOME

                        The name of the reference file that matches the

                        reference used to map the Illumina inputs.

  --fai FAI             (Optional) The name of the corresponding index for the

                        reference genome file.

  --prefix PREFIX       (Optional) If provided, all output files will start

                        with this. If absent, the base of the BAM file name

                        will be used.

  --filter_short_contigs

                        If selected, SV calls will not be generated on contigs

                        shorter than 1 MB.

  --breakdancer         If selected, the program Breakdancer will be one of

                        the SV callers run.

  --breakseq            If selected, the program BreakSeq2 will be one of the

                        SV callers run.

  --manta               If selected, the program Manta will be one of the SV

                        callers run.

  --cnvnator            If selected, the program CNVnator will be one of the

                        SV callers run.

  --lumpy               If selected, the program Lumpy will be one of the SV

                        callers run.

  --delly_deletion      If selected, the deletion module of the program Delly2

                        will be one of the SV callers run.

  --delly_insertion     If selected, the insertion module of the program

                        Delly2 will be one of the SV callers run.

  --delly_inversion     If selected, the inversion module of the program

                        Delly2 will be one of the SV callers run.

  --delly_duplication   If selected, the duplication module of the program

                        Delly2 will be one of the SV callers run.

  --genotype            If selected, candidate events determined from the

                        individual callers will be genotyped and merged to

                        create a consensus output.

  --svviz               If selected, visualizations of genotyped SV events

                        will be produced with SVVIZ, one screenshot of support

                        per event. For this option to take effect, Genotype

                        must be selected.

  --svviz_only_validated_candidates

                        Run SVVIZ only on validated candidates? For this

                        option to be relevant, SVVIZ must be selected. NOT

                        selecting this will make the SVVIZ component run

                        longer.

——

 

実行方法

ランするのに必要なものは、.bamと.bam.bai、そのリファレンスのfasta.gzまたはfa.gzとおよび.faiファイルになる。fastaファイルは非圧縮にも対応している(*3)。

 

breakdancer、cnvnator、mantaでSVコールする。さらに結果を統合し、フィルタリングしてSVコールとそのリード分布可視化ファイルを得る。

#出力ディレクトリを指定する
cd /path/to/bam_data/ && mkdir output

docker run -v ${PWD}:/home/dnanexus/in \
-v ${PWD}/output/:/home/dnanexus/out dnanexus/parliament2:0.1.7 \
--bam input.bam --bai input.bam.bai -r ref_genome.fa.gz \
--breakdancer --cnvnator --manta \
--genotype --svviz --svviz_only_validated_candidates

fasta.gzを指定するとラン中に解凍される。2回目のランの時は注意する。

 

ジョブが終わると、各SVコーラー出力のvcfを保存したディレクトリ、マージしてgenotypingしたvcf、コール部位を1箇所ずつ可視化してtarにまとめたファイルなどが出力される。

f:id:kazumaxneo:20181025203902j:plain

 

 

 個別のSVコーラーのランも簡単に実行できる。

Cnvnatorのランのみ実行

docker run 
-v /home/uesaka/bam_data/:/home/dnanexus/in \
-v /home/uesaka/bam_data/output/:/home/dnanexus/out \
dnanexus/parliament2:0.1.7 \
--bam input.bam --bai input.bam.bai -r ref_genome.fa.gz \
--cnvnator

 

Mantaのランのみ実行

docker run 
-v /home/uesaka/bam_data/:/home/dnanexus/in \
-v /home/uesaka/bam_data/output/:/home/dnanexus/out \
dnanexus/parliament2:0.1.7 \
--bam input.bam --bai input.bam.bai -r ref_genome.fa.gz \
--cnvnator

 

全解析

docker run 
-v /home/uesaka/bam_data/:/home/dnanexus/in \
-v /home/uesaka/bam_data/output/:/home/dnanexus/out \
dnanexus/parliament2:0.1.7 \
--bam input.bam --bai input.bam.bai -r ref_genome.fa.gz \
--breakdancer --cnvnator --manta --breakseq \
--delly_deletion --delly_insertion --delly_inversion --delly_duplication \
--genotype --svviz

breakseqは1000genomeのhs37d5をリファレンスにした時のみ、結果を得られる(*5)。

 

補足

DNAnexus社の研究者が発表しているので、当然DNAnexusのクラウド解析サーバーでも使えるようになっています(*4)。確認のため、DNAnexusに登録して実行までの流れを調べましたが、マウスでファイルを追加してパラメータを決め、Runボタンを押すだけ実行直前まで進めました。費用の検討は必要ですが、この方法が一番簡単に使えそうです (*2)。データは、AWSやAzureのサーバーに置いてあっても使えるようです。

f:id:kazumaxneo:20181025144644j:plain

 

DNAnexusには、Parliament2以外にも様々なワークフローがready to goの状態で準備されています。例えば、最近論文が出たGoogleのDeepVariant(pubmed) (下の画像)も使えます (HP)(Googleから多額の出資があることも関係していそうですね *1)

f:id:kazumaxneo:20181025145332j:plain

 

必要な費用ですが、環境によって変わってきます(例えばスペック、転送量、使用時間など)。このあたり、他のクラウドコンピューティング環境と変わりませんね。

引用

Parliament2: Fast Structural Variant Calling Using Optimized Combinations of Calle

Samantha Zarate, Andrew Carroll, Olga Krashenina, Fritz J Sedlazeck, Goo Jun, William Salerno, Eric Boerwinkle, Richard Gibbs

bioRxiv preprint first posted online Sep. 23, 2018; doi:

doi.org/10.1101/424267.

 

2020 12/22

Parliament2: Accurate structural variant calling at scale
Samantha Zarate, Andrew Carroll, Medhat Mahmoud, Olga Krasheninina, Goo Jun, William J Salerno, Michael C Schatz, Eric Boerwinkle, Richard A Gibbs, Fritz J Sedlazeck
GigaScience, Volume 9, Issue 12, December 2020

 

参考

*1

マイクロソフトやGoogle、遺伝子情報プラットフォームのDNAnexusに5800万ドル出資 | HealthTechNews

マイクロソフト、ゲノム研究用クラウドツールの一般提供開始を発表 - News Center Japan

 

*2

無料だとデータのtransfer機能のみ使用できます。log in後、右上のupgrade nowからスペックを指定し、支払いをすませると解析できるようになります。

 

*3

svtyperのマージ付近でエラーが出たので、fastaのヘッダーに特殊文字が含まれないよう修正したら最後までランできるようになった。

 

*4

DNAnexus社のサーバーはAWSを使っているらしい。

 

*5

Google genomicsからダウンロードできる。

https://googlegenomics.readthedocs.io/en/latest/use_cases/discover_public_data/reference_genomes.html

 

*6

既にオーサーの方の1人はGoogle genomicsに移られたようなことをツイートされていますね(未確認)。間違ってたらすみません。

 

関連

Dockerについては前に簡単に説明しています。


大きなk-merも使うde Bruijn graph のアセンブリツール SKESA

 2019 4/12 dockerリンク追加

 

 NGSデータを分析するためのシーケンスアライメント、アセンブリ、変異検出、またはそれらのいくつかの組み合わせは、通常、バイオインフォマティクスパイプラインの主要なモジュールである[論文より ref.1,2,3,4,5,6]。微生物ゲノムシーケンシングの重要な用途は、食物サプライチェーン[ref.7,8,9]および病院[ref.10,11,12,13]での病原菌の大発生を検出することである。NGSを使った食品媒介性病原菌のサーベイランスおよびアウトブレイク調査の利点およびバイオインフォマティクスの課題を、Listeria Monocytogenesを例として、またリアルタイムアウトブレイク[ref.15]例を引用して検討した。両方のレビューでは、情報を使用する上での重大な挑戦であるNGSの新しいアセンブリが確認された。

 米国の州、連邦機関、および国際パートナー間の協力で、食品由来病原菌のシーケンシングデータをsubmitするNCBIのPathogen Detection Project (PDP)によりNGSベースのアウトブレイクの調査が加速している。サルモネラ、リステリア、エシェリヒアおよびシゲラ、カンピロバクターの4つの主要な食品媒介性病原菌でsubmitされたシーケンシングデータセットの数は、2016年の初めにPDPによって発表された44017から、2018年の85823に急拡大している。

(一部略)

 シーケンシングリードのde-novoアセンブラがいくつか公開されている[ref.21,22,23,24,25,26,27,28]。ploidity[ref.29]やメタゲノム[ref.30,31,32,33,34,35]、シングルセル[ref.36]、シーケンシング技術[ref.37]、またはいくつかのアセンブリを1つに結合することに特化してたりと[ref.38]、それぞれ特徴がある。アセンブラは、一倍体ゲノムであっても「エラーのない」アセンブリを保証しない。微生物ゲノムに加えて、ハプロイドアセンブリは、ヒトの特別なケース: hydatidiform mole (胞状奇胎) のような場合に重要である[39]。 PDPのようなユースケースでは、微生物ゲノムのアプリケーションは、母系から娘細胞へのゲノムデータの垂直継承に依存し、4Mbゲノムではほとんど変化のない、典型的には変異が10未満ということに基づいてリードセットクローン性パターンを解決する[ref.15]。このようなアプリケーションでは、非常に高いシーケンス品質のアセンブリが必要となるが、真のバリエーションを確実に検出することができる。

 シーケンシングマシンによって生成されたリードセットを使用してアセンブリを評価するには、同じサンプルのリードとほぼ完全な高品質のドラフトアセンブリの両方を含む公的に利用可能なベンチマークセットが必要になる。 FDA-ARGOS(HP link)は、Food and Drug Administration(US FDA)[ref.19]によって開発された、微生物ゲノムのアセンブリ品質評価のベンチマーク要件を満たすデータベースであり、バクテリアのregulatory grade のシーケンスから構成されている。

 ここでは、イルミナシーケンシング技術を使用して生成された微生物ゲノムのリードを迅速かつ高品質にデノボアセンブリするSKESA [skee-sa](strategic k-mer extension for scrupulous assemblies))を紹介する。 SKESAで使用されるHeuristicsは、Illuminaシーケンシングの低レベルのコンタミネーションやストランド固有のエラーによるアセンブリ品質への影響が低減されるように設計されている。エラー率の高い他のシーケンシングテクノロジでは、SKESAで使用される控えめなヒューリスティックにより、他のアセンブラで生成されたものよりアセンブリのcontiguityは低くなる。 SKESAは微生物ゲノムより大きなゲノムもアセンブリできるが、他のアセンブラと比較して、そのようなゲノムのためのプロファイリングはされていない。例えば、SKILSAは約40Mbの長さであるMonilinia fructigenaのSRR7262862データ(49.8 million reads with total length 9.2 Gb)とMonilinia laxaのSRR6748693データ(73.2 million reads with total length 18.3 Gb)を、10コア使用で3時間以下でアセンブリし、100コア使用で30分でアセンブリする。

 多くのデノボアセンブラ、例えばSPAdes [ref.24]やMegaHit [ref.35]は、アセンブリ時に複数のk-merサイズを使用する。 ALLPATHS_LG [ref.40]は、100塩基対が40塩基オーバーラップしている特別短いインサートライブラリープロトコルを使用する。オーバーラップを利用して160 bpのマージされたリードが生成されたが、最大k-merサイズは96が使用された。 SKESAの顕著な特徴は、リードのサブセットのミニアセンブリを使ってペアエンドのインサートサイズまでの長いk-merを生成することである。リードより長いk-merを使用することで、SKESAは、インサートサイズよりも短いが、リードよりも長いリピート領域を正確にアセンブリできる。これに対して、現在のアセンブラはすべて、ペアエンドのmateリードサイズまでのk-mersだけ使用している。

 本稿では、SKESAとSPAdes、MegaHitを、5つのタイプの微生物検査セットに対して適用した。(一部略) MegaHitは非常に高速なアセンブラであるため、SPAdesは、さまざまなテクノロジやサンプルのオプションを提供する多目的で広く使用されているアセンブラのため、この2つのアセンブラを選択した。ラージゲノムとスモールゲノムの両方をアセンブリできるよう設計されたABySS [ref.28](バージョン2.0.2)と不均一なカバレッジを扱うように設計されたIDBA [36](バージョン1.1.1)もテストされた(Additional file 1)。

(一部略)

 著者らは、微生物ゲノムのアセンブリでは、SKESAとMegaHitは同じくらい速く、SPAdesよりもはるかに速いことを示す。 SKESAはSRAから直接アクセスすることもできる。これにより、ファイルからの入力を読み取るよりも高速に処理できる。 QUAST [ref.42]によって計算された100Kbあたりのミスマッチ数、N50統計量によって測定されたアセンブリの連続性、およびリファレンスアセンブリの長さからのずれによって測定されたアセンブリ品質は、SKESAアセンブリの品質がSPAdesとMegaHitの品質より良好であることを示している。同じ入力では、スレッド数、メモリ、または実行回数に関係なく、同じ結果が得られる。大量のデータを処理し、回帰テストを必要とする本番システムでは、これは非常に重要な要件である。著者らのテストでは、SPAdesとMegaHitの両方は、同じ数のスレッドとメモリの設定であっても、繰り返しにより異なるアセンブリを生成していた。したがって、SKESAは、高い基本レベルの品質と下流の分析に十分なcontiguityを備え、低レベルの汚染を扱い、迅速かつ高精度なアセンブリを作成するすべての要件を満たしている。すなわち生産環境でrobustである。 SKESAは現在、SRAの微生物ゲノムを組み立てるためにNCBIのproductionに使用されており、PDPのワークフローに組み込まれている。 SKESAのソフトウェアは自由に利用可能であり[ref.43、44](see “Availability and requirements”)、クラウドでも利用できるようにする予定である。

 

 

 

インストール

mac os10.12のanaconda2-4.2.0 環境でテストした。

Github

git clone https://github.com/ncbi/SKESA
cd SKESA/

#If you would like to build NGS library for accessing reads from SRA, then do
make
#Otherwise, if reading inputs only from files, do
make -f Makefile.nongs

#またはcondaで導入(sraの直接アセンブリ機能はなし)
conda install -c bioconda skesa

 > skesa

$ skesa

skesa 

 

Provide some input reads

 

General options:

  -h [ --help ]                 Produce help message

  -v [ --version ]              Print version

  --cores arg (=0)              Number of cores to use (default all) [integer]

  --memory arg (=32)            Memory available (GB, only for sorted counter) 

                                [integer]

  --hash_count                  Use hash counter [flag]

  --estimated_kmers arg (=100)  Estimated number of unique kmers for bloom 

                                filter (M, only for hash counter) [integer]

  --skip_bloom_filter           Don't do bloom filter; use --estimated_kmers as

                                the hash table size (only for hash counter) 

                                [flag]

 

Input/output options : at least one input providing reads for assembly must be specified:

  --fasta arg                   Input fasta file(s) (could be used multiple 

                                times for different runs) [string]

  --fastq arg                   Input fastq file(s) (could be used multiple 

                                times for different runs) [string]

  --use_paired_ends             Indicates that a single (not comma separated) 

                                fasta/fastq file contains paired reads [flag]

  --contigs_out arg             Output file for contigs (stdout if not 

                                specified) [string]

 

Assembly options:

  --kmer arg (=21)              Minimal kmer length for assembly [integer]

  --min_count arg               Minimal count for kmers retained for comparing 

                                alternate choices [integer]

  --max_kmer_count arg          Minimum acceptable average count for estimating

                                the maximal kmer length in reads [integer]

  --vector_percent arg (=0.05)  Count for  vectors as a fraction of the read 

                                number (1. disables) [float (0,1]]

  --insert_size arg             Expected insert size for paired reads (if not 

                                provided, it will be estimated) [integer]

  --steps arg (=11)             Number of assembly iterations from minimal to 

                                maximal kmer length in reads [integer]

  --fraction arg (=0.1)         Maximum noise to signal ratio acceptable for 

                                extension [float [0,1)]

  --max_snp_len arg (=150)      Maximal snp length [integer]

  --min_contig arg (=200)       Minimal contig length reported in output 

                                [integer]

  --allow_snps                  Allow additional step for snp discovery [flag]

 

Debugging options:

  --force_single_ends           Don't use paired-end information [flag]

  --seeds arg                   Input file with seeds [string]

  --all arg                     Output fasta for each iteration [string]

  --dbg_out arg                 Output kmer file [string]

  --hist arg                    File for histogram [string]

  --connected_reads arg         File for connected paired reads [string]

 

実行方法

ペアエンドのfastqを指定する。

skesa --fastq pair1.fq,pair2.fq --cores 8 --memory 48 > skesa_assembly.fa

 

SRAのシーケンシングデータを直接アセンブリする。

skesa --sra_run SRR1960353 --cores 4 --memory 48 > SRR1960353.skesa.fa

* --sra_runのフラグはビルド時に -f Makefile.nongsを指定していると存在しない。

 

2つ以上のSRAを指定

skesa --sra_run SRR1695624 --sra_run SRR1745628 --cores 4 --memory 48 > SAMN03218571.skesa.fa 

他の例はGIthubのREADMEを参考にして下さい。 

 

 

 

簡単なベンチマーク

他のスモールゲノム用アセンブラとごくごく簡単に比較してみる。シーケンシングデータはbacteriaのアセンブリコンペティション GAGE-B(リンク)で使われたRhodobacter sphaeroidesのMiSeqシーケンシングデータを使う( Miseqのペアエンドデータ: SRA link)。MiSeqのリードならオーバーラップが十分あり、SKESAのアルゴリズムが活かせると考えた。SKESAやMEGAHITの特徴であるランタイムの短さはここでは調べない(*1)。

 

1、fastqの準備。

#1 GATGE-Bよりダウンロード
mkdir test && cd test/
wget https://ccb.jhu.edu/gage_b/datasets/R_sphaeroides_MiSeq.tar.gz
tar -zxvf R_sphaeroides_MiSeq.tar.gz
mv R_sphaeroides_MiSeq/raw/insert_540* .

#2 reenameと圧縮
mv insert_540_1__cov100x.fastq pair1.fq && gzip pair1.fq
mv insert_540_2__cov100x.fastq pair2.fq && gzip pair2.fq

#リファレンス
wget ftp://ftp.ncbi.nih.gov/genomes/Bacteria/Rhodobacter_sphaeroides_2_4_1_uid57653
#ftpのlinkが消えていたのでNCBI Refseqからダウンロードした

#===================================================================
#GAGE-Bでなく、例えばSRAをデータ "SRR1703350" をダウンロードするなら
fastq-dump --split-spot --skip-technical --split-files SRR1703350

#SRAからダウンロードした場合、BBsketchを使いベストマッチするゲノムを推定(リンク
sendsketch.sh in=SRR1703350_1.fq reads=100k
#ANI近似値が最大のゲノムをダウンロードする

pair1.fq.gzとpair2.fq.gzができる。

 

2、de novo assembly

#1 skesa
skesa --fastq pair1.fq.gz,pair2.fq.gz --cores 12 --memory 32 > skesa.fa

#2 megahit(紹介)
megahit -1 pair1.fq.gz -2 pair2.fq.gz -o megahit_out --k-min 21 --k-max 141 --k-step 12 -t 12

#3 spades1(紹介)
spades.py -k 21,31,41,51,61,71,81,91,101,111,127 \
-1 pair1.fq.gz -2 pair2.fq.gz \
-t 12 --careful -o spades1

#4 spades2
#spades v3.12で追加された--mergedのリードを組み込む。
#merged.fqの作成
bbmerge.sh in1=pair1.fq.gz in2=pair2.fq.gz \
out=merged.fq outu1=unmerged1.fq outu2=unmerged2.fq \
ihist=ihist.txt
#圧縮
gzip merged.fq
#assembly
spades.py -k 21,31,41,51,61,71,81,91,101,111,127 \
-1 pair1.fq.gz -2 pair2.fq.gz --merged merged.fq.gz \
-t 12 --careful -o spades2

#5 unicycler(紹介)
unicycler -1 pair1.fq.gz -2 pair2.fq.gz -l merged.fq.gz \
-o unicycler -t 12

#6 stride(紹介)
stride all pair1.fq.gz pair2.fq.gz -r 250 -i 550 -t 12

#a7 A5miseq(紹介)
a5_pipeline.pl pair1.fq.gz pair2.fq.gz a5miseq_out_dir

spadesは、v3.12からmerged.fqを使いパフォーマンスを上げるオプションがついたので、"--merged merged.fq"の有無による差も比較した。

 

3、Quastによる評価(リンク

#gene情報もあるが、ここでは無しで実行する。

quast skesa.fa megahit.fasta spades1.fasta spades2.fasta unicycler.fasta StriDe.fasta a5output.fasta\
-R GCF_000012905.2_ASM1290v2_genomic.fna \
-1 pair1.fq.gz -2 pair2.fq.gz \
-o quast_output \
-t 12

 report.html

f:id:kazumaxneo:20181025175926j:plain

SKESAは比較したアセンブラの中で最もキメラアセンブリ、ミスマッチ、indelsが少なかった。

f:id:kazumaxneo:20181025180151j:plain

連続性に関しては、NG50がspades やunicyclerの約半分だった。

 

icarus.html

f:id:kazumaxneo:20181025180013j:plain

全体を眺めていくと、下の図のように、SPAdesでアセンブリできていない領域がSKESAではアセンブリできていることもあった。

f:id:kazumaxneo:20181025181939j:plain

 

引用
SKESA: strategic k-mer extension for scrupulous assemblies
Souvorov A, Agarwala R, Lipman DJ

Genome Biol. 2018 Oct 4;19(1):153

 

*1

主観では、spadesより4倍~5倍以上短い時間で終わる。