macでインフォマティクス

macでインフォマティクス

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

Dockerイメージをリバースエンジニアリングする Whaler

 

ブログで開発の動機は説明されています。 

Reversing Docker Images into Dockerfiles

https://samaritan.ai/blog/reversing-docker-images-into-dockerfiles/

 

インストール

mac os10.13のgo version go1.11 darwin/amd64でビルドした。

Go言語がないなら先に導入しておく。

Github

go get -u github.com/P3GLEG/Whaler
cd ~/go/src/github.com/P3GLEG/Whaler
go build .

> Whaler -h

$ Whaler -h

Usage of Whaler:

  -f string

    File containing images to analyze seperated by line

  -filter

    Filters filenames that create noise such as node_modules. Check ignore.go file for more details (default true)

  -sV string

    Set the docker client ID to a specific version -sV=1.36

  -v Print all details about the image

  -x Save layers to current directory

 

実行方法

分析するにはdokcerコンテナイメージ名を指定する。

Whaler image名

 

ubuntuコンテナイメージを使い試してみる。まずイメージを検索する。

sudo docker search ubuntu

 

$ sudo docker search ubuntu

NAME                                                   DESCRIPTION                                     STARS               OFFICIAL            AUTOMATED

ubuntu                                                 Ubuntu is a Debian-based Linux operating sys…   8426                [OK]                

dorowu/ubuntu-desktop-lxde-vnc                         Ubuntu with openssh-server and NoVNC            224                                     [OK]

rastasheep/ubuntu-sshd                                 Dockerized SSH service, built on top of offi…   171                                     [OK]

consol/ubuntu-xfce-vnc                                 Ubuntu container with "headless" VNC session…   130                                     [OK]

ansible/ubuntu14.04-ansible                            Ubuntu 14.04 LTS with ansible                   95                                      [OK]

ubuntu-upstart                                         Upstart is an event-based replacement for th…   90                  [OK]                

neurodebian                                            NeuroDebian provides neuroscience research s…   54                  [OK]                

1and1internet/ubuntu-16-nginx-php-phpmyadmin-mysql-5   ubuntu-16-nginx-php-phpmyadmin-mysql-5          44                                      [OK]

ubuntu-debootstrap                                     debootstrap --variant=minbase --components=m…   39                  [OK]                

nuagebec/ubuntu                                        Simple always updated Ubuntu docker images w…   23                                      [OK]

tutum/ubuntu                                           Simple Ubuntu docker images with SSH access     18                                      

i386/ubuntu                                            Ubuntu is a Debian-based Linux operating sys…   14                                      

1and1internet/ubuntu-16-apache-php-7.0                 ubuntu-16-apache-php-7.0                        13                                      [OK]

ppc64le/ubuntu                                         Ubuntu is a Debian-based Linux operating sys…   12                                      

eclipse/ubuntu_jdk8                                    Ubuntu, JDK8, Maven 3, git, curl, nmap, mc, …   6                                       [OK]

codenvy/ubuntu_jdk8                                    Ubuntu, JDK8, Maven 3, git, curl, nmap, mc, …   4                                       [OK]

darksheer/ubuntu                                       Base Ubuntu Image -- Updated hourly             4                                       [OK]

1and1internet/ubuntu-16-nginx-php-5.6-wordpress-4      ubuntu-16-nginx-php-5.6-wordpress-4             3                                       [OK]

pivotaldata/ubuntu                                     A quick freshening-up of the base Ubuntu doc…   2                                       

1and1internet/ubuntu-16-sshd                           ubuntu-16-sshd                                  1                                       [OK]

ossobv/ubuntu                                          Custom ubuntu image from scratch (based on o…   0                                       

smartentry/ubuntu                                      ubuntu with smartentry                          0                                       [OK]

1and1internet/ubuntu-16-healthcheck                    ubuntu-16-healthcheck                           0                                       [OK]

paasmule/bosh-tools-ubuntu                             Ubuntu based bosh-cli                           0                                       [OK]

pivotaldata/ubuntu-gpdb-dev                            Ubuntu images for GPDB development              0                                       

 

docker社が用意したイメージ、STARSが8426のubuntuをダウンロードする。

sudo docker pull ubuntu:16.04

イメージ名とイメージIDを確認。

docker images

$ docker images

REPOSITORY               TAG                 IMAGE ID            CREATED             SIZE

ubuntu-kazumax           latest              77672f0a2762        35 hours ago        2.45GB

ubuntu                   16.04               b9e15a5d1e1a        2 weeks ago         115MB

jeltje/varscan2          latest              6df55adfc7d5        5 weeks ago         725MB

ubuntu                   14.04               971bb384a50a        2 months ago        188MB

centos                   latest              49f7960eb7e4        3 months ago        200MB

lavenderca/muver         latest              fd0c2415cdde        4 months ago        1.24GB

metabat/metabat          latest              d9259b96b365        5 months ago        60.1MB

nileshtawari/chronqc     chronqc_1.0.4       a2bb7bfc794c        9 months ago        890MB

eitanbanks/manta-1.0.3   latest              61cd9a5d5f6f        19 months ago       288MB

molecular/breakdancer    latest              28543e51bb72        23 months ago       658MB

 

ubuntu16.04コンテナイメージを分析する。

Whaler -sV=1.36 b9e15a5d1e1a

$ Whaler -sV=1.36 b9e15a5d1e1a

Analyzing b9e15a5d1e1a

Docker Version: 17.06.2-ce

GraphDriver: overlay2

Environment Variables

|PATH=/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin

 

Image user

|User is root

 

Potential secrets:

Dockerfile:

RUN set -xe  \

&& echo '#!/bin/sh' > /usr/sbin/policy-rc.d  \

&& echo 'exit 101' >> /usr/sbin/policy-rc.d  \

&& chmod +x /usr/sbin/policy-rc.d  \

&& dpkg-divert --local --rename --add /sbin/initctl  \

&& cp -a /usr/sbin/policy-rc.d /sbin/initctl  \

&& sed -i 's/^exit.*/exit 0/' /sbin/initctl  \

&& echo 'force-unsafe-io' > /etc/dpkg/dpkg.cfg.d/docker-apt-speedup  \

&& echo 'DPkg::Post-Invoke { "rm -f /var/cache/apt/archives/*.deb /var/cache/apt/archives/partial/*.deb /var/cache/apt/*.bin || true"; };' > /etc/apt/apt.conf.d/docker-clean  \

&& echo 'APT::Update::Post-Invoke { "rm -f /var/cache/apt/archives/*.deb /var/cache/apt/archives/partial/*.deb /var/cache/apt/*.bin || true"; };' >> /etc/apt/apt.conf.d/docker-clean  \

&& echo 'Dir::Cache::pkgcache ""; Dir::Cache::srcpkgcache "";' >> /etc/apt/apt.conf.d/docker-clean  \

&& echo 'Acquire::Languages "none";' > /etc/apt/apt.conf.d/docker-no-languages  \

&& echo 'Acquire::GzipIndexes "true"; Acquire::CompressionTypes::Order:: "gz";' > /etc/apt/apt.conf.d/docker-gzip-indexes  \

&& echo 'Apt::AutoRemove::SuggestsImportant "false";' > /etc/apt/apt.conf.d/docker-autoremove-suggests

RUN rm -rf /var/lib/apt/lists/*

RUN sed -i 's/^#\s*\(deb.*universe\)$/\1/g' /etc/apt/sources.list

RUN mkdir -p /run/systemd  \

&& echo 'docker' > /run/systemd/container

CMD ["/bin/bash"]

 

 

インフォマティクスのツールには、公共のDocker Hubにイメージは共有されていても、Dockerfile(イメージビルドの時に実行した内容が記載される)が別に公開されてないイメージもたくさんあります。特に、よく使われるツールで、開発者以外の方が登録したイメージに多い印象です。仮想環境といえど、そもそもそのようなイメージは使うべきではないという議論はさておき、例えば妙にファイルサイズが大きなイメージ、ツールがどこに入っているかわからないイメージなどは気になりますね。そういったイメージもビルド手順みれば一目瞭然です。例えば使わないツールが親切心でapt-getされていて、それで容量が肥大しているなら、除いてカスタムビルドするだけでかなりディスクスペースを節約できます(立ち上げて修正し、eixtしてコミットし直してもいいと思いますが)。 活用して下さい。

 

引用

GitHub - P3GLEG/Whaler: Program to reverse Docker images into Dockerfiles

インタラクティブなDNA配列の2次元プロットを作成する Squiggle

 

 次世代シークエンシング技術の登場により、DNA配列解析は、バイオインフォマティクスと生物学の両方でますます一般的なツールとなっている。この理由から、注釈されていないDNA配列を迅速に検査する能力は極めて重要である。しかし、FASTAファイルに含まれるようなrawシーケンシングデータは、視覚的に類似した文字の長いシーケンスであり、人間が探索するにはあまり適していない。例えば、複数の高度に保存されたコード配列を見る場合、肉眼で見るのが困難な異なる文字が比較的少なくてもよい。配列データをよりよく表示するために、DNAを構成するAs、Ts、Gs、およびCsの長い配列を2次元のグラフィカル表現に変えることができる多数の視覚化アルゴリズムがある(Gates、1986; Yauetal、2003; Randic et al 、2003b; Qi and Qi、2007)。
文献に記載された数多くのrawDNAシーケンシング可視化方法にもかかわらず、文献の総説はこれらの方法のオープンソースソフトウェア実装を特定しなかった。したがって、オープンソースコマンドラインツールであるSquiggleと、塩基配列を対話型のブラウザベースの2次元視覚化に変換するPythonパッケージを提案する。 

Squiggleは、raw DNA配列のインタラクティブなWebベースの2次元グラフィカル表現を自動的に生成するソフトウェアツールである。使いやすさを考慮して構築されたSquiggleは、いくつかの従来のシーケンス視覚化アルゴリズムを実装し、人間の使い勝手を最大限にするように設計された新しいvisualization手法を導入している。

 

マニュアル

https://squiggle.readthedocs.io/en/latest/

 

インストール

mac os10.13のpython3.4環境下でテストした。

Github

pip install squiggle

> squiggle -h

$ squiggle -h

Usage: squiggle [OPTIONS] [FASTA]...

 

Options:

  -w, --width FLOAT               The width of the line. Defaults to 1.

  -p, --palette TEXT              Which color palette to use. Choose from boke

                                  h.pydata.org/en/latest/docs/reference/palett

                                  es.html. Defaults to Category10.

  --color / --no-color            Whether to plot the visualizations in color.

  --hide / --no-hide              Whether to hide sequences when clicked in

                                  the legend. Defaults to false if plotting

                                  one sequence and true if plotting multiple.

  --bar / --no-bar                Whether to show a progress bar when

                                  processing multiple sequences. Defaults to

                                  true.

  -t, --title TEXT                A title to display when plotting sequences

                                  together.

  --separate                      Whether to plot the visualizations

                                  separately.

  -c, --cols INTEGER              The number of columns when plotting

                                  separately. Defaults to a the closest to

                                  square layout as possible.

  --link-x / --no-link-x          Whether to link the x axes for separate

                                  plotting. Defaults to true.

  --link-y / --no-link-y          Whether to link the y axes for separate

                                  plotting. Defaults to false.

  -o, --output PATH               The output file for the visualization. If

                                  not provided, will open visualization in

                                  browser. The filetype must be .html

  --offline                       Whether to include the resources needed to

                                  plot offline when outputting to file.

                                  Defaults to false.

  --method [squiggle|gates|yau|yau-bp|randic|qi]

                                  The visualization method.

  -d, --dimensions WIDTH HEIGHT   The width and height of the plot,

                                  respectively. If not provided, will default

                                  to 750x500.

  --skip / --no-skip              Whether to skip any warnings. Defaults to

                                  false.

  --mode [seq|file|auto]          Whether to treat each sequence or file as

                                  the individual object. Defaults to automatic

                                  selection.

  --legend-loc [top_left|top_center|top_right|center_right|bottom_right|bottom_center|bottom_left|center_left|center]

                                  Where to put the legend, if applicable.

                                  Defaults to top left.

  --output-backend [canvas|svg|webgl]

                                  Which output backend to use while plotting.

                                  Defaults to canvas.

  -s, --downsample INTEGER        The downsampling factor. Useful for handling

                                  large sequences. Default value is 1.

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

 

 

実行方法

分析対象のFASTAを指定してsquiggleを叩く。(オフライン環境なら"--offline"をつける必要あり)。

squiggle human_HBB.fasta -o output.html

"-o outout.html"を除くと、デフォルトブラウザに標準出力される。

f:id:kazumaxneo:20180924183028p:plain

squiggle/example_seqsにあるテストデータを使っている。グラフは、インタラクティブなグラフを描けるBokehを使って作られていて、右上のボタンからリアルタイムに移動や拡大縮小ができる。

 

複数ゲノム

squiggle human_HBB.fasta chimpanzee_HBB.fasta norway_rat_HBB.fasta rhesus_HBB.fasta

f:id:kazumaxneo:20180924183134p:plain

 

分割描画

squiggle human_HBB.fasta chimpanzee_HBB.fasta --separate

f:id:kazumaxneo:20180924183428p:plain

出力のボタン操作はグラフ間で連動します(目盛りのリンクは--no-link-xと--no-link-yで無効化)。 

 

2ゲノム間の高発現な遺伝子を、26遺伝子までランダムサンプリングして比較

squiggle human_HBB.fasta chimpanzee_HBB.fasta --separate

 

グラフの高さと幅、タイトル、色を変更

squiggle human_HBB.fasta -d 650 150 -t β-globin -p Accent
  • -d <WIDTH HEIGHT>   The width and height of the plot, respectively. If not provided, will default to 750x500.
  • -t  TEXT     A title to display when plotting sequences together.
  • -p TEXT     Which color palette to use. Choose from bokeh.pydata.org/en/latest/docs/reference/palettes.html. Defaults to Category10.

f:id:kazumaxneo:20180924184617p:plain

Visualisation methodsはオンラインマニュアルにまとめられている。defaultでは"--method squiggle"(the 2bit encoding scheme that maps T to 00, C to 01, A to 10, and G to 11)だが、他のmethodに変更することもできる(link1, link2, link3)。

https://squiggle.readthedocs.io/en/latest/methods.html

 

 

引用

Squiggle: a user-friendly two-dimensional DNA sequence visualization tool

Benjamin D Lee

Bioinformatics 2018 ACCEPTED MANUSCRIPT

 

 

既知変異情報を利用して精度を上げたバリアントコールを行う IVC

 

 ゲノムのバリアント検出は、ゲノミクス、バイオインフォマティクス、生物医学研究およびその応用(1000 Genomes Project Consortium、2012,2015; Pabinger et al、2014)において非常に重要な意味を持つ。次世代シークエンシング(NGS)技術の最近の進歩により、ヒトや他の種に対して費用効果が高く、ハイスループットで大規模なシーケンシングデータを提供することが可能になった(Auton et al、2012; Lek et al、2016)。これは、ゲノムの変異検出を含む幅広いゲノミクスおよびバイオインフォマティクス研究を容易にする。 SNPs、IndelsおよびStructural Variantsを含むゲノムのバリアントコールのための多くのアルゴリズムおよびツールが開発されている。それらの大部分はベイズ法、隠れマルコフモデル、ロジスティック回帰モデルなどの統計的方法に基づいている(Bansal et al、2010; Challis et al、2012; Chen et al、2009; DePristo et al、2011; Garrison Marth、2012; Li、2011; Mose et al、2014; Wang et al、2013; Ye et al、2009)。これらの方法は十分に確立されているが、ゲノムの反復的または高密度に変化した領域におけるIndelsおよび構造変異を検出することには依然として苦しんでいる。 Li(2014)は、お互いに近接したIndelsではアライナーが正しいアライメントを取得するのが困難で、その結果、バリアントコーラーが正しいバリアントを取得するのを困難にしていることを示した。最近、Jiangら(2015)は、Indelsが現在評価されているよりも豊富であることを示している。さらに、Chaissonら(2014)は、一分子リアルタイムDNAシーケンシングを用いて一倍体のヒトゲノム(CHM1)を分析することにより、以前に報告されたよりも数倍多いIndelsを検出した。これらバリアントが欠けていると、腫瘍進化の特徴付けや治療反応の予測などの不正確な下流解析につながる可能性がある。従って、Indel検出の感度を改善する必要性が依然として存在する。

 現在のバリアントコールメソッドのほとんどは、外部のアライナーによるアライメント情報に基づいてバリアントを検出する(Pabinger et al、2014)。アラインメントプロセスを高速化するため、ほとんどのアライナは同じ領域からのそれぞれのリード間の依存関係は考えない。これにより、一貫性のないアライメントが生成され、バリアントコールが複雑になる(Li、2014)。研究者は、通常、関心のある領域(Albers et al、2011; DePristo et al、2011)でのミスアライメントを回復するためにリアライメントを実行する。いくつかの方法は、ローカルアセンブリを実行しunitigを構築し、それから、unitigをリファレンスゲノムにマッピングしてバリアントコールする(Carnevali et al、2012; Mose et al、2014)。他のいくつかの方法は、long Indels検出を改善するためにアライメントグラフと呼ばれる特別なデータ構造を利用する(Marschall et al、2013)。研究者はまた、バリアントコールを支援するため複数サンプルを利用する(Bansal et al、2010; Li、2011; Wang et al、2013)。それにもかかわらず、それらのツールは依然として、さまざまな呼び出しプロセスをサポートするための入力データの活用が十分ではない。そのため、通常、バリアントを正確に検出できるようにするには、高いシーケンシングカバレッジが必要になる。

 現在、ヒトおよび他の種の集団における多数のゲノムバリアントが公的データベースに集められている。例えば、dbSNP(Wheeler et al、2007)、dbVarおよびDGVa(Lappalainen et al、2013)、1000 Genomes Project(1999)などの多くの努力を払って多数のヒト個体から多くの変異が検出されている(1000ゲノムプロジェクトコンソーシアム、2015)またはExAC(Lek et al、2016)。したがって、新しいゲノムのバリアントコールをサポートするためにそれらを効率的に利用できる方法を設計することが望ましいであろう。この目的のために、SOAPsnp(Li et al、2009b)は、ヒトのデータの分析においてdbSNPを使用することにより、事前genotype確率を計算する。 Atlas(Shen et al、2010)は、トレーニングデータセットの事前SNP確率と事前エラー確率を使用している。それにもかかわらず、これらの方法は、既知のバリアント情報を効率的に利用してバリアントコールをサポートしていない。特に、これらのツールは、既知のバリアントを使用してのアライメント分析をサポートしていない。最近では、BWBBLE(Huang et al、2013)がリードのアラインメントをサポートするために既知のバリアントを利用できる少ないツールの1つだが、このメソッドはバリアントコール問題に適応していない。

 本稿では、既知バリアント情報を利用して既知のバリアントと未知のバリアントの両方を、特にIndels近傍で、検出する新規の方法を紹介する。この方法は、リファレンスゲノムとそれに関連する既知バリアントプロファイルを適切な方法で組み合わせて、リードアライメントとバリアントコールを効率的かつ正確に実行する手助けをする。次に、既知のバリアント情報を考慮に入れるアルゴリズムを開発して、リードのリファレンスへのアライメントを行う。この戦略により、リードをより正確にアライメントさせることができ、従来の方法と比較して少ないリード数(低いカバレッジ)であっても、バリアントをより正確にコールすることができる。特に、著者らの方法は、高いカバレッジでもコールが難しい、Indels近傍のコールが改善されている。低カバレッジであっても高精度なこの方法により、変異検出の実験コストが削減できる。

 

 IVC(Integrated Variant Calling)と呼ばれる本発明者らの方法は、(i)リファレンスゲノムと既知のバリアントプロファイルを適切な方法で組み合わせるメタ-ゲノム表現方法、および(ii)既知のバリアント感受性のあるアラインメントアルゴリズム、これは、最適なアラインメントを決定する際に既知のバリアントを考慮に入れることができるもの、の2つを主要なcontributionsに含んでいる。両方とも、標準的なリファレンスゲノムおよび従来のアラインメントアルゴリズムでは困難な、リードアライメントおよびバリアント検出プロセスにおいて既知バリアントを効率的かつ効果的に利用するように設計されている。 (iii)効率的にシードを見つけることができる反復ランダム化アルゴリズム、および(iv)アライメントプロセス中にバリアントプロファイルをアップデートするストラテジーは、リードアライメントの精度を向上させるのを助けることができ、それゆえバリアントコール精度も上がる 。方法については、論文の以降のセクションで説明する。

 

インストール

ubuntu16.04でテストした。

本体 GIthub

GO言語で実装されている。

インストールの詳細はGIthub参照。指示に従えばできる。

> ./ivc -h

$ ./ivc -h

2018/09/23 14:11:09 IVC - Integrated Variant Caller using next-generation sequencing data.

2018/09/23 14:11:09 IVC-main: Calling variants based on alignment between reads and reference multi-genomes.

Usage of ./ivc:

  -1 string

    pairend read file, first end

  -2 string

    pairend read file, second end

  -I string

    index directory

  -O string

    variant call output file

  -R string

    reference genome file

  -V string

    variant profile file

  -d float

    threshold of alignment distances

  -debug

    turn on debug mode.

  -e float

    gap extension cost

  -lmax int

    maximum length of seeds

  -lmin int

    minimum length of seeds

  -maxp int

    maximum number of paired-seeds

  -maxs int

    maximum number of seeds

  -mode int

    searching mode for finding seeds (1: random (default), 2: deterministic)

  -o float

    gap open cost

  -r int

    maximum number of iterations

  -s float

    substitution cost

  -start int

    starting position on reads for finding seeds

  -step int

    step for searching in deterministic mode

  -t int

    maximum number of CPUs

./ivc-index -h

$ ./ivc-index -h

2018/09/23 14:12:09 IVC - Integrated Variant Caller using next-generation sequencing data.

2018/09/23 14:12:09 IVC-index: Indexing reference genomes and variant profiles.

Usage of ./ivc-index:

  -I string

    index directory

  -R string

    reference genome file

  -V string

    variant profile file

  -debug

    turn on debug mode.

 

 

 

実行方法

1、リファレンスゲノムのindex作成

cd $GOPATH/src/github.com/namsyvo/IVC

ivc-index -R test_data/refs/chr1_ref.fasta \
-V test_data/refs/chr1_variant_prof.vcf \
-I test_data/indexes
  • -R     reference genome (FASTA format).
  • -V     known variant profile (VCF format).
  • -I      directory for storing index.

既知バリアントのvcfを指定している。indexes/にindexファイルが出力される。

f:id:kazumaxneo:20180923142131j:plain

 

2、バリアントコール

ivc.go -R test_data/refs/chr1_ref.fasta \
-V test_data/refs/chr1_variant_prof.vcf \
-I test_data/indexes \
-1 test_data/reads/chr1_dwgsim_100_0.001-0.01.bwa.read1.fastq \
-2 test_data/reads/chr1_dwgsim_100_0.001-0.01.bwa.read2.fastq \
-O test_data/results/chr1_variant_calls.vcf
  • -R     reference genome (FASTA format).
  • -V     known variant profile (VCF format).
  • -I     directory for storing index.
  • -1     the read file (for single-end reads) (FASTQ format).
  • -2     the second end file (for pair-end reads) (FASTQ format).
  • -O     variant call result file (VCF format).

既知バリアントのvcf、作成したindexesディレクトリ、リードファイルを指定している。results/にバリアントコールの.vcfファイルが出力される。

f:id:kazumaxneo:20180923142717j:plain

既知バリアントと合致するポジションも含めてコールされる。

 

引用

Leveraging known genomic variants to improve detection of variants, especially close-by Indels
Vo NS, Phan V

Bioinformatics. 2018 Mar 24

 

関連ツール


バクテリアのシーケンスエンリッチメント解析ツール SEER

 

 細菌形質の遺伝的基盤を決定しようとする研究は、伝統的に特定の原因遺伝子エレメントを特定するのではなく、目的の表現型に関連する新生クローンを同定することに限定されてきた(ref.1)。これは、バクテリアがクローン的に複製するという事実に起因しており、これはゲノムの大部分が任意の形質 (ref.2) と連鎖不平衡(LD)にあることを意味している。形質に関連する変異のリストのうちのどれが真に形質と因果関係であるかを決定するには、その形質が単一のクローン系統に特有に関連していないことを必要とする。(一部略)病原性のような臨床的に関連する形質の場合のように、表現型に対して完全に浸透 (浸透率 penetration rate) していない変異を見出すには、多数の試料を必要とし、より一般的な関連検査が必要となる(ref.4)。

これらの理由から、バクテリアの表現型のゲノムワイド関連研究(GWAS)は近年になってわずかになってきている(ref.2,5,6,7,8)。もともとヒトの一塩基多型(SNP)データのために開発された標準的なGWAS法は、バクテリアのゲノム突然変異に首尾よく適用可能であることが示されている(ref.6,7)。しかしながら、既知のバクテリア種の多くのゲノムの可塑性が高いことから、そのような方法は表現型の変異の遺伝的決定因子を部分的にしか同定できないと予想される。例えば遺伝子内容に関連する機構の発見を可能にするために、アライメントフリーの代替法も導入されている(ref.5,8)。これらの方法は、観察された表現型分布における差異の推定的説明に、一般的にk-mer、すなわち長さkのDNAワードを使用する。 k-mersの主な利点は、突然変異、indel、組換え、可変プロモーター構造および遺伝子含量の差異を含むゲノムのコレクション全体に存在するいくつかの異なるタイプの変異を捕捉する能力である。

K-mersは、アセンブリ (ref.9)、SNPjコール (ref.10) および距離推定 (ref.11) のためのバクテリアゲノミクスに使用されている。これまでのk-mersを用いた研究は、母集団構造をコントロールし系統推定を行うためにワードの増加と消失のモンテカルロシミュレーションを用いていた(ref.5)。一方で、SNPベースの研究では、コアアライメントに関するクラスタリングアルゴリズムと、得られたサンプルグループについてのstratified association testが使用されている( ref.6,7)。前者は、より小さなエフェクトサイズ(効果量)の関連を見つける必要があり、数百の分離株まで計算的にスケーリングできない。後者は、コアアライメントを探す感度が不十分であり、また多数のサンプルがある場合には生成が困難であり、 そもそも多様性があるかもしれない。

 ここでは、新たにアセンブリされたコンティグまたはrawシーケンシングデータを入力として使用するスタンドアロンパイプライン実装の、計算的に数万のゲノムへスケーラブルできるsequence element enrichment analysisツール(SEER)を紹介する。著者らは、SEERを大規模で多様な集団のシミュレーテッドデータおよびリアルデータの両方に適用し、遺伝子の存在およびコード領域におけるSNPの両方によって引き起こされる抗生物質耐性との関連性を正確に検出することができることを示す。

 

Usage

https://github.com/johnlees/seer/wiki/Usage

 

インストール

ubuntu18.04でテストした。

依存

You will also require

  • gcc >4.9 or equivalent
  • gcc libstdc++ >4.9

スタティックとダイナミックのバイナリファイルをダウンロードできる。自分でビルドする場合はGithubを参照。

https://github.com/johnlees/seer/releases

> ./seer -h

$ ./seer -h

seer: sequence element enrichment analysis

 

Required options:

  -k [ --kmers ] arg       dsm kmer output file

  -p [ --pheno ] arg       .pheno metadata

 

Covariate options:

  --struct arg             mds values from kmds

  --covar_file arg         file containing covariates

  --covar_list arg         list of columns covariates to use. Format is 1,2q,3 

                           (use q for quantitative)

 

Performance options:

  --threads arg (=1)       number of threads. Suggested: 2

 

Filtering options:

  --no_filtering           turn off all filtering and perform tests on all 

                           kmers input

  --max_length arg (=100)  maximum kmer length

  --maf arg (=0.01)        minimum kmer frequency

  --min_words arg          minimum kmer occurences. Overrides --maf

  --chisq arg (=10e-5)     p-value threshold for initial chi squared test. Set 

                           to 1 to show all

  --pval arg (=10e-8)      p-value threshold for final logistic test. Set to 1 

                           to show all

 

Other options:

  --print_samples          print lists of samples significant kmers were found 

                           in

  --version                prints version and exits

  -h [ --help ]            full help message

 

またはオーサーらが準備してくれている仮想イメージを利用する。(VirtualBox

ftp://ftp.sanger.ac.uk/pub/pathogens/pathogens-vm/pathogens-vm.latest.ova

 

ラン

以下の手順で解析を進める。

  1. k-merカウント。
  2. kmdsを走らせて.dsmファイル作成。
  3. .dsmファイルを使い、もう一度kmdsを走らせる。distance matrixテーブルが出力される。
  4. seerで解析。

 

1、初めに、アセンブリして作成したcontig配列をターゲットkにしてk-merカウントを実施する。クラスター環境でランできるなら、推奨されているようにdsmリンク)を使用する。一般的なシングルサーバ環境ならfsm-lite(リンク)を使用する。

fsm-lite -v -s 5 -S 95 -l fsm_list.txt -t tmp_index > fsm_kmers.txt

The -s and -S options control the frequency of k-mers to report. In the above example there are 100 samples and only k-mers appearing in between 5 and 95 samples are output. This is equivalent to the recommended 5% minor allele frequency cutoff.

 

 使えないならdskを使う。k-merサイズは固定となる。

dsk -file sample1_contigs.fa -abundance-min 1 -out sample1_dsk
dsk2ascii -file sample1_dsk.h5 -out sample1_dsk.txt
combineKmers -r 21mer_samples.txt -o all_21mers_new --min_samples 2

 

2、kmdsによる全k-merのサブサンプリングして.dsmファイルを作成。

distance.matrixを2ステップで作成する。まず分割してgzipする。

split -d -n l/16 fsm_kmers.txt fsm_out
rm fsm_kmers.txt
gzip fsm_out*

それからkmdsを走らせる。

kmds -k fsm_out{1..16}.gz --pheno metadata.pheno --no_filtering --no_mds --size 10000

 {1..16}は1から16まで16回ランすることを意味する。

2の作業については、オーサーらはmash(紹介)を使った方が計算が早く終わり好ましいとしてmashを使った手順も紹介している(リンク)。

 

 3、 2で作った.dsmを使い、distance matrixを計算。

ls *.dsm > subsampled_matrices.txt
kmds -p metadata.pheno --mds_concat subsampled_matrices.txt -o all_structure --threads 16 --write_distances

 

4、seerで解析する。

seer -k fsm_out{1..16}.gz --pheno metadata.pheno --struct all_structure --threads 4 --print_samples > significant_kmers.txt

Various filtering options are available. Note the the default is a MAF of >1%, unadjusted p-value of < 10e-5 and adjusted p-value of < 10e-8. 

 

 

インストールが簡単になったpyseerも発表されています (link)。また紹介します。

引用

Sequence element enrichment analysis to determine the genetic basis of bacterial phenotypes

Lees JA, Vehkala M, Välimäki N, Harris SR, Chewapreecha C, Croucher NJ, Marttinen P, Davies MR, Steer AC, Tong SY, Honkela A, Parkhill J, Bentley SD, Corander J

Nat Commun. 2016 Sep 16;7:12797.

https://www.nature.com/articles/ncomms12797

 

高感度で高速なプロテイン検索を行う MMseqs2

 

 DNAシーケンシングのスループットは、過去10年間で計算速度よりもはるかに速くなってきており、感度の高いシーケンス検索は、ラージメタゲノムデータセットの分析における主要なボトルネックになっている。それゆえ、著者らは、速度と感度のトレードオフの全レンジににわたって現在の検索ツールを改良し、400以上の倍の速度で、PSI-BLASTよりも良好な感度を達成し400倍以上高速なMMseqs2を開発した。

 2007年以来、シーケンシングコストが4桁低下した結果、テラバイトの配列を生産する大規模なメタゲノムプロジェクトが、医療、生物工学、微生物学、農業研究に応用されて多数行われている(論文より ref.1,2,3,4) 。計算分析の中心的なステップは、データベース内の類似の配列を検索することによってオープンリーディングフレームの機能を推論するアノテーションである。メタゲノミクスでは、現在、計算コストがシーケンシングコストに対して支配的になっており(ref.5,6,7)、その中でタンパク質検索は、感度は高いが低速なBLASTがはるかに高速な検索ツールで置き換えられているにもかかわらず、一般に計算資源の90%を消費している(ref,9,10,11,12)。しかし、それらのツールのスピードは感度低下を犠牲にして行われている。メタゲノミクスとメタトランスクリプトミクス研究で見つかる多くの種は、よくアノテーションづけられたゲノムとはclosely relatedではないので、unannotatableシーケンスの割合は、多くの場合、65から90パーセントに昇る(ref.2,13)。シーケンシングと計算コストとのギャップ拡大により、すぐにこの問題は悪化する。

 この課題に対処するために、著者らはープンソースの並列化ソフトウエアスイートMMseqs2を開発した。MMseqs2は以前のMMseqs(ref.14 pubmed)と比較してはるかに高感度で、反復的なプロファイル間シーケンスとプロファイル間検索をサポートし、さらに強化された機能を提供する(論文のSupplementary Table 1)。

 MMseqs2のサーチは3つのステージ、短いワード( 'k-mer')マッチステージ、ベクトル化されたungappedアライメント、およびgapped(Smith-Waterman)アライメントという3つのステージから構成されている(論文 図1a)。第1段階は、改善されたパフォーマンスにとって重要である。与えられたクエリ配列に対して、同じ対角線上に2つの連続する類似のk-mer一致を有するすべての標的配列を見つける(論文 図1b)。連続的なk-merの一致は、相同配列の場合、同じ対角線上にあることがよくあるが(それらの間にアライメントギャップがない場合)、偶然に起こる可能性は低い。ほとんどの高速なツールは正確なk-merマッチのみを検出するが、MMseqs2、MMseqs、BLASTは類似k-mer間のrマッチも見つける。このような類似k-mer照合は、MMseqs2が、類似性設定に応じて、多数の同様のk-merを生成することによって(1クエリに対して600 から 60,000)感度を失うことなく大きな単語サイズk = 7を使用することを可能にする。 MMseqs2の速度については、最もの内側のループ4(Supplementary Information Figure S 1 のマゼンタ色のフレーム内のループ)の最後の行のランダムメモリアクセスを排除する方法を発見することが重要だった。

 

f:id:kazumaxneo:20180922182940p:plain

The average AUC sensitivity versus speed-up factor relative to BLAST for 637,000 test searches.  The sensitivity is called area under the curve  (AUC).

ユーザーガイドより転載。AUCの定義についてはSupplementary Information参照。

 

MMseqs2/Linclust updates by following Martin


ユーザーガイド

https://github.com/soedinglab/mmseqs2/wiki#set-sensitivity--s-parameter

 

インストール

mac os 10.13でテストした。

MMseqs can be installed by compiling the binary from source, download a statically compiled version, using Homebrew or Docker.

MMseqs2  Githubレポジトリ

ソースからのビルドに必要なもの

  • git
  • g++ (4.6 or higher)
  • cmake (3.0 or higher)
#ソースからビルドする 
git clone https://github.com/soedinglab/MMseqs2.git
cd MMseqs2
mkdir build
cd build
cmake -DCMAKE_BUILD_TYPE=RELEASE -DCMAKE_INSTALL_PREFIX=. ..
make
make install
export PATH=$(pwd)/bin/:$PATH

他に、homebrewによるインストール、static バイナリイメージダウンロード、dockerイメージが利用できる。手順はMMseqs2レポジトリのマークダウンファイル (README.md)参照。

mmseqs

$ mmseqs 

MMseqs2 (Many against Many sequence searching) is an open-source software suite for very fast, 

parallelized protein sequence searches and clustering of huge protein sequence data sets.

 

Please cite: M. Steinegger and J. Soding. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature Biotechnology, doi:10.1038/nbt.3988 (2017).

 

MMseqs2 Version: aa4e079915dcfaba34fa58206144180b64ef6258

© Martin Steinegger (martin.steinegger@mpibpc.mpg.de)

 

Easy workflows (for non-experts)

  easy-search       Search with a query fasta against target fasta (or database) and return a BLAST-compatible result in a single step

  easy-linclust     Compute clustering of a fasta database in linear time. The workflow outputs the representative sequences, a cluster tsv and a fasta-like format containing all sequences.

  easy-cluster      Compute clustering of a fasta database. The workflow outputs the representative sequences, a cluster tsv and a fasta-like format containing all sequences.

 

Main tools  (for non-experts)

  createdb          Convert protein sequence set in a FASTA file to MMseqs sequence DB format

  search            Search with query sequence or profile DB (iteratively) through target sequence DB

  map               Fast ungapped mapping of query sequences to target sequences.

  cluster           Compute clustering of a sequence DB (quadratic time)

  linclust          Cluster sequences of >30% sequence identity *in linear time*

  createindex       Precompute index table of sequence DB for faster searches

  clusterupdate     Update clustering of old sequence DB to clustering of new sequence DB

 

Utility tools for format conversions

  createtsv         Create tab-separated flat file from prefilter DB, alignment DB, cluster DB, or taxa DB

  convertalis       Convert alignment DB to BLAST-tab format, SAM flat file, or to raw pairwise alignments

  convertprofiledb  Convert ffindex DB of HMM/HMMER3/PSSM files to MMseqs profile DB

  convert2fasta     Convert sequence DB to FASTA format

  result2flat       Create a FASTA-like flat file from prefilter DB, alignment DB, or cluster DB

  createseqfiledb   Create DB of unaligned FASTA files (1 per cluster) from sequence DB and cluster DB

 

Taxonomy tools      

  taxonomy          Compute taxonomy and lowest common ancestor for each sequence.

  lca               Compute the lowest common ancestor from a set of taxa.

 

 

An extended list of all tools can be obtained by calling 'mmseqs -h'.

 

Bash completion for tools and parameters can be installed by adding "source MMSEQS_HOME/util/bash-completion.sh" to your "$HOME/.bash_profile".

Include the location of the MMseqs2 binary is in your "$PATH" environment variable.

> mmseqs createdb -h

mmseqs createdb -h

mmseqs createdb:

converts a protein sequence flat/gzipped FASTA or FASTQ file to the MMseqs sequence DB format. This format is needed as input to mmseqs search, cluster and many other tools.

 

Please cite:

Steinegger, M. & Soding, J. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature Biotechnology, doi:10.1038/nbt.3988 (2017)

 

© Martin Steinegger <martin.steinegger@mpibpc.mpg.de>

 

Usage: <i:fastaFile1[.gz]> ... <i:fastaFileN[.gz]> <o:sequenceDB> [options]

 

misc options             default   description [value range]

  --dont-split-seq-by-len true      Dont split sequences by --max-seq-len                       

  --dont-shuffle         true      Do not shuffle input database                               

  --id-offset            0         numeric ids in index file are offset by this value          

 

common options           default   description [value range]

  --max-seq-len          65535     Maximum sequence length [1,32768]                           

  -v                     3         verbosity level: 0=nothing, 1: +errors, 2: +warnings, 3: +info

mmseqs search -h

$ mmseqs search -h

mmseqs search:

Searches with the sequences or profiles query DB through the target sequence DB by running the prefilter tool and the align tool for Smith-Waterman alignment. For each query a results file with sequence matches is written as entry into a database of search results (alignmentDB).

In iterative profile search mode, the detected sequences satisfying user-specified criteria are aligned to the query MSA, and the resulting query profile is used for the next search iteration. Iterative profile searches are usually much more sensitive than (and at least as sensitive as) searches with single query sequences.

 

Please cite:

Steinegger, M. & Soding, J. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature Biotechnology, doi:10.1038/nbt.3988 (2017)

 

© Martin Steinegger <martin.steinegger@mpibpc.mpg.de>

 

Usage: <i:queryDB> <i:targetDB> <o:alignmentDB> <tmpDir> [options]

 

prefilter options       default   description [value range]

  --comp-bias-corr      1         correct for locally biased amino acid composition [0,1]     

  --add-self-matches    false     artificially add entries of queries with themselves (for clustering)

  -s                    5.700     sensitivity: 1.0 faster; 4.0 fast default; 7.5 sensitive [1.0,7.5]

  -k                    0         k-mer size in the range [6,7] (0: set automatically to optimum)

  --k-score             2147483647 k-mer threshold for generating similar-k-mer lists          

  --alph-size           21        alphabet size [2,21]                                        

  --offset-result       0         Offset result list                                          

  --split               0         Splits input sets into N equally distributed chunks. The default value sets the best split automatically. createindex can only be used with split 1.

  --split-mode          2         0: split target db; 1: split query db;  2: auto, depending on main memory

  --split-memory-limit  0         Maximum system memory in megabyte that one split may use. Defaults (0) to all available system memory.

  --diag-score          1         use diagonal score for sorting the prefilter results [0,1]  

  --exact-kmer-matching 0         only exact k-mer matching [0,1]                             

  --mask                1         0: w/o low complexity masking, 1: with low complexity masking

  --min-ungapped-score  15        accept only matches with ungapped alignment score above this threshold

  --spaced-kmer-mode    1         0: use consecutive positions a k-mers; 1: use spaced k-mers 

 

align options           default   description [value range]

  -a                    false     add backtrace string (convert to alignments with mmseqs convertalis utility)

  --alignment-mode      2         How to compute the alignment: 0: automatic; 1: only score and end_pos; 2: also start_pos and cov; 3: also seq.id; 4: only ungapped alignment

  -e                    0.001     list matches below this E-value [0.0, inf]                  

  --min-seq-id          0.000     list matches above this sequence identity (for clustering) [0.0,1.0]

  --seq-id-mode         0         0: alignment length 1: shorter, 2: longer sequence          

  --alt-ali             0         Show up to this many alternative alignments                 

  -c                    0.000     list matches above this fraction of aligned (covered) residues (see --cov-mode)

  --cov-mode            0         0: coverage of query and target, 1: coverage of target, 2: coverage of query 3: target seq. length needs be at least x% of query length

  --realign             false     compute more conservative, shorter alignments (scores and E-values not changed)

  --max-rejected        2147483647 maximum rejected alignments before alignment calculation for a query is aborted

  --max-accept          2147483647 maximum accepted alignments before alignment calculation for a query is stopped

  --score-bias          0.000     Score bias when computing the SW alignment (in bits)        

 

profile options         default   description [value range]

  --pca                 1.000     pseudo count admixture strength                             

  --pcb                 1.500     pseudo counts: Neff at half of maximum admixture (0.0,infinity)

  --mask-profile        1         mask query sequence of profile using tantan [0,1]           

  --e-profile           0.001     includes sequences matches with < e-value thr. into the profile [>=0.0]

  --wg                  false     use global sequence weighting for profile calculation       

  --filter-msa          1         filter msa: 0: do not filter, 1: filter                     

  --max-seq-id          0.900     reduce redundancy of output MSA using max. pairwise sequence identity [0.0,1.0]

  --qid                 0.000     reduce diversity of output MSAs using min.seq. identity with query sequences [0.0,1.0]

  --qsc                 -20.000   reduce diversity of output MSAs using min. score per aligned residue with query sequences [-50.0,100.0]

  --cov                 0.000     filter output MSAs using min. fraction of query residues covered by matched sequences [0.0,1.0]

  --diff                1000      filter MSAs by selecting most diverse set of sequences, keeping at least this many seqs in each MSA block of length 50

  --num-iterations      1         Search iterations                                           

 

misc options            default   description [value range]

  --no-preload          false     Do not preload database                                     

  --rescore-mode        0         Rescore diagonal with: 0: Hamming distance, 1: local alignment (score only) or 2: local alignment

  --min-length          30        minimum codon number in open reading frames                 

  --max-length          32734     maximum codon number in open reading frames                 

  --max-gaps            2147483647 maximum number of codons with gaps or unknown residues before an open reading frame is rejected

  --contig-start-mode   2         Contig start can be 0: incomplete, 1: complete, 2: both     

  --contig-end-mode     2         Contig end can be 0: incomplete, 1: complete, 2: both       

  --orf-start-mode      0         Orf fragment can be 0: from start to stop, 1: from any to stop, 2: from last encountered start to stop (no start in the middle)

  --forward-frames      1,2,3     comma-seperated list of ORF frames on the forward strand to be extracted

  --reverse-frames      1,2,3     comma-seperated list of ORF frames on the reverse strand to be extracted

  --translation-table   1         1) CANONICAL, 2) VERT_MITOCHONDRIAL, 3) YEAST_MITOCHONDRIAL, 4) MOLD_MITOCHONDRIAL, 5) INVERT_MITOCHONDRIAL, 6) CILIATE, 9) FLATWORM_MITOCHONDRIAL, 10) EUPLOTID, 11) PROKARYOTE, 12) ALT_YEAST, 13) ASCIDIAN_MITOCHONDRIAL, 14) ALT_FLATWORM_MITOCHONDRIAL, 15) BLEPHARISMA, 16) CHLOROPHYCEAN_MITOCHONDRIAL, 21) TREMATODE_MITOCHONDRIAL, 22) SCENEDESMUS_MITOCHONDRIAL, 23) THRAUSTOCHYTRIUM_MITOCHONDRIAL, 24) PTEROBRANCHIA_MITOCHONDRIAL, 25) GRACILIBACTERIA, 26) PACHYSOLEN, 27) KARYORELICT, 28) CONDYLOSTOMA, 29) MESODINIUM, 30) PERTRICH, 31) BLASTOCRITHIDIA

  --use-all-table-starts false     use all alteratives for a start codon in the genetic table, if false - only ATG (AUG)

  --id-offset           0         numeric ids in index file are offset by this value          

  --add-orf-stop        false     add * at complete start and end                             

  --start-sens          4.000     start sensitivity                                           

  --sens-steps          1         Search steps performed from --start-sense and -s.           

  --slice-search        false     For bigger profile DB, run iteratively the search by greedily swapping the search results.

 

common options          default   description [value range]

  --sub-mat             blosum62.out amino acid substitution matrix file                         

  --max-seq-len         65535     Maximum sequence length [1,32768]                           

  --max-seqs            300       maximum result sequences per query (this parameter affects the sensitivity)

  --threads             1         number of cores used for the computation (uses all cores by default)

  -v                    3         verbosity level: 0=nothing, 1: +errors, 2: +warnings, 3: +info

 

mmseqs convertalis -h

$ mmseqs convertalis -h

mmseqs convertalis:

Convert alignment DB to BLAST-tab format, SAM flat file, or to raw pairwise alignments

 

Please cite:

Steinegger, M. & Soding, J. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature Biotechnology, doi:10.1038/nbt.3988 (2017)

 

© Martin Steinegger <martin.steinegger@mpibpc.mpg.de>

 

Usage: <i:queryDb> <i:targetDb> <i:alignmentDB> <o:alignmentFile> [options]

 

misc options   default   description [value range]

  --format-mode 0         output format 0: BLAST-TAB, 1: PAIRWISE, 2: BLAST-TAB + query/db length

  --no-preload false     Do not preload database                                     

  --db-output  false     Output a result db instead of a text file                   

 

common options default   description [value range]

  --threads    1         number of cores used for the computation (uses all cores by default)

  -v           3         verbosity level: 0=nothing, 1: +errors, 2: +warnings, 3: +info

 

 

実行方法

元々、クラスター環境での大規模なデータベース構築を想定して設計されているが、このブログではローカル環境でのコマンド実行方法のみ紹介する。

 

検索するには、まずクエリとデータベースそれぞれの配列をmmseq2のformatに変換する必要がある。テストファイルに対して実行してみる。

cd MMseqs2/examples/

#クエリのアミノ酸配列fastaのフォーマット変換
mmseqs createdb QUERY.fasta queryDB

#データベースのアミノ酸配列fastaのフォーマット変換
mmseqs createdb DB.fasta targetDB

 

クエリとデータベースのファイルを指定してホモロジーサーチを実行する。tmpは作業ディレクトリ。巨大なデータベースを使う場合、十分な容量とI/Oの高速なストレージの利用が推奨されている(詳細はユーザーガイド参照)。

mmseqs search queryDB targetDB resultDB tmp
  • -s    sensitivity: 1.0 faster; 4.0 fast default; 7.5 sensitive (default 5.7)

感度は-sで1.0から7.5まで変更できる (マニュアルの説明)。目的に合わせて変更する。上の表にもあるように、感度により計算時間も変わる。

$ mmseqs createdb examples/DB.fasta targetDB

Program call:

createdb examples/DB.fasta targetDB 

 

MMseqs Version:              aa4e079915dcfaba34fa58206144180b64ef6258

Max. sequence length         65535

Split Seq. by len            true

Do not shuffle input database true

Offset of numeric ids        0

Verbosity                    3

 

.Time for merging files: 0h 0m 0s 3ms

Time for merging files: 0h 0m 0s 4ms

Time for merging files: 0h 0m 0s 4ms

Time for merging files: 0h 0m 0s 3ms

Time for processing: 0h 0m 0s 126ms

 

結果をblast様のタブ仕分けファイルに変換する。

mmseqs convertalis queryDB targetDB resultDB resultDB.m8 --format-mode 0
  • --format-mode 0   output format 0: BLAST-TAB, 1: PAIRWISE, 2: BLAST-TAB + query/db length 

 

引用

MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets
Steinegger M, Söding J

Nat Biotechnol. 2017 Nov;35(11):1026-1028

 

MMseqs software suite for fast and deep clustering and searching of large protein sequence sets.
Hauser M, Steinegger M, Söding J

Bioinformatics. 2016 May 1;32(9):1323-30.

 

Clustering huge protein sequence sets in linear time

Martin Steinegger & Johannes Söding

Nat Commun. 2018; 9: 2542.

 

Kostabl labのANIやAAIを計算するwebツール (enveomicsコレクションの一部)

 ゲノムおよびメタゲノム解析は、生物学的研究のいくつかの分野でますます一般的になってきているが、頻繁に繰り返される特殊な分析は、論文のpublish後にはほとんど利用できないin-houseスクリプトとして報告されている。著者らは、微生物ゲノミクスおよびメタゲノミクスにおけるいくつかの反復的タスクおよび専門のタスクに積極的に維持されるスクリプトセットであるenveomicsコレクションについて説明し、またグラフィカルユーザインタフェースおよびいくつかの事例研究を提示する。(一部略)

enveomicsコレクションは、Artistic License 2.0(https://github.com/lmrodriguezr/enveomics)およびオンライン分析(http://enve-omics.ce.gatech.edu)の条件で自由に利用できる。

 

Toolsにアクセスする。

http://enve-omics.ce.gatech.edu

10近いツールが利用できるが、ここではANIとAAIの計算ツールに限定して紹介する。

 

1、ANI: Average Nucleotide Identity calculator

2ゲノム間のANI計算ツール

http://enve-omics.ce.gatech.edu/ani/

f:id:kazumaxneo:20180921163516p:plain

比較対象の2ゲノムのFastAファイルを選択する。NCBI accession numberを指定することもできる。windowサイズは1000でMinimum identityは70%だが、変えることもできる。

f:id:kazumaxneo:20180921163748p:plain

 

デモファイルの結果を見てみる。ANIの分布が棒グラフで表示され、ANIの平均値と中央値がプリントされている。

f:id:kazumaxneo:20180921164014p:plain

それぞれの領域のアライメント結果もダウンロードできる。

f:id:kazumaxneo:20180921164007p:plain

 

 

2、AAI: Average Amino acid Identity calculator

http://enve-omics.ce.gatech.edu/aai/

f:id:kazumaxneo:20180921163542p:plain

AAIはアミノ酸配列同士の比較。実行するには対象生物のアミノ酸FastAファイルを準備する必要がある。AAI計算のMinimum identityは20%、Minimum alignmentsは50a.aになっている。

 

 

3、All-vs-all ANI/AAI matrix calculator

http://enve-omics.ce.gatech.edu/g-matrix/

ゲノム総当たりのANIまたはAAIの比較。

f:id:kazumaxneo:20180921163607p:plain

入力のFastAをzipやgzipなどでフォルダごと圧縮して指定する。

"対応フォーマット: .zip, .tar, .tar.gz, and .tar.bz2. Packages must include only FastA files (up to 50 genomes)"

 

結果は行列 (matrix)と系統樹で出力される。計算量が多いため、かなり時間がかかる(*1)。テストデータ(link

f:id:kazumaxneo:20180921201148p:plain

f:id:kazumaxneo:20180921201206p:plain

 

引用

The enveomics collection: a toolbox for specialized analyses of microbial genomes and metagenomes

Luis M Rodriguez-R​, Konstantinos T Konstantinidis
PeerJ Preprints March 27, 2016

 

参考記事

http://enve-omics.gatech.edu/sites/default/files/2014-Rodriguez_R-Konstantinidis_Microbe_Magazine.pdf

 

*1

知人が総当たりAAIを実行した時は、結果のメールが来るまで2日かかったそうです(10-20程度の生物のタンパク質FastA)。

 

MetaMeta

 

 現在、環境サンプルをcharacterizeすることを目指して、ますます多くのメタゲノム分析ツールが利用可能になっている[論文より ref.1,2,3,4]。Whole metagenome shotgun (WMS)シーケンシングテクニックから生成される大量のデータにより動機づけられたメタゲノムのプロファイリングは、実際のシナリオでよりアクセスしやすく、より速く適用可能となり、メタゲノミクス分析の標準的な方法となってきている[ref.5,6,7]。 WMSシーケンシングデータに基づいて配列の分類を実行するツールにさまざまな種類がある。 1つの基本的なアプローチは、リファレンスや先行知識なしにショートリードから完全またはほぼ完全なゲノムを再構築することを目的とするde novo sequence assembly [ref.8,9,10]である。これはコミュニティ構成を評価するための最高の解決策を提供する。しかし、リード長の短さ、不十分なカバレッジ、類似のDNA配列、およびlow abundant な株のために、メタゲノミクスデータから有意義なアセンブリを生成することは非常に困難である[ref.11 link]。

 より一般的には、アセンブリなしで直接WMSリードを使う、リファレンスベースの手法がある。つまり、以前に取得されたゲノム配列に依存して解析が行われる。このカテゴリのアプリケーションでは、2つのスタンダード:taxonomic profilingツールとビニングツールが採用されている。Profilersは、WMSシーケンス全体を分析し、所定のリファレンス配列に基づいて生物およびそれらの相対存在量を予測することを目指している。ビニングツールは、特定のサンプル中の各シーケンスを個別に分類し、それらのそれぞれをリファレンスセットの最も有望な生物に連結することを目的とする。これら概念の違いにかかわらず、両方のツール群を微生物群集の特徴づけに使用することが可能だった。しかし、ビニングツールは各配列の個々の分類を生成し、分類学的プロファイラとして使用するために変換し、正規化する必要がある。

 これらの2つのカテゴリの中で利用可能な方法は、いくつかの技術、例えば、リードマッピング、k-merアライメント、および組成分析が含まれる。データベース、例えばcomplete ゲノム配列、マーカー遺伝子、タンパク質配列の構築に関する変形もまた一般的である。これらの技術の多くは、最新のシークエンシング技術の高いスループットと、利用可能な多数のリファレンスゲノム配列を処理するための計算コストを克服するために開発された。

 ツール、パラメータ、データベース、およびテクニックのためのいくつかのオプションの利用可能性は、研究者にとってどのメソッドを使用するか決めるシナリオを複雑にする。さまざまなツールは、さまざまなシナリオで良好な結果を提供し、複数の構成で多かれ少なかれ正確であるか、感度が高い。研究やサンプルの変化ごとにそれらの出力に依存することは難しい。さらに、複数の方法が使用される場合、異なる基準セットを使用するツール間の矛盾した結果を統合することは困難である。さらに、インストール、パラメータ、データベース作成、標準出力の欠如は容易に解決できない課題である。

 メタゲノムシーケンス分類ツールの共同実行と統合のための新しいパイプラインであるMetaMetaを提案する。 MetaMetaにはいくつかの利点がある:簡単なインストールとセットアップ、複数のツールとサンプルとデータベースのサポート、複数の結果を組み合わせた最終プロファイルの改善、すぐに使用できる並列化と高性能コンピューティング(HPC)の統合、データベースの自動ダウンロードと設定、カスタムデータベースの作成、統合された前処理ステップ(リードのトリミング、エラー訂正、およびサブサンプリング) 。MetaMetaは、正しい識別情報をマージし、誤った識別情報を適切にフィルタリングすることにより、より慎重なプロファイリング結果を実現する。 MetaMetaはSnakeMake [ref.12 link]で構築され、オープンソースである。(以下略)

 

Tool selectionより

MetaMetaは、CLARK [ref.19]、DUDes [ref.20]、GOTTCHA [ref.21]、Kraken [ref.22]、Kaiju [ref.23]、およびmOTUs [ref.24]の6つのツールセットで評価された。 選んだ動機は、部分的には、そのようなツールの性能を比較している最近の論文による[ref.3、4、25]。 CLARK、GOTTCHA、Kraken、およびmOTUは[ref.4]に従い非常に低い偽陽性数を達成した。 DUDesは、[ref.25]に従って精度と感度との間の良好なトレードオフを達成するin houseツールであった。 Kaijuはtranslatedデータベースを使用し、全ゲノムベースの方法に多様性をもたらす。

 

f:id:kazumaxneo:20180920211149j:plain

MetaMeta Pipeline. 論文より転載(pubmed)。

 

MetaMetaに関するツイート


インストール

cent os6のanaconda3-4.0.0環境でテストした。

本体 Github

#Anaconda環境ならcondaで導入可能だが、linuxのみ。
conda install -c bioconda metameta
#他とバッティングしていたようなので、仮想環境にインストールし直した


#仮想環境でテストする。pyenvでpythonのバージョン管理しているものとする
conda config --add channels defaults
conda config --add channels conda-forge
conda config --add channels bioconda

conda create -n metametaenv metameta=1.2.0
source activate metametaenv

 >metameta --help

$ metameta --help

usage: snakemake [-h] [--profile PROFILE] [--snakefile FILE] [--gui [PORT]]

                 [--cores [N]] [--local-cores N]

                 [--resources [NAME=INT [NAME=INT ...]]]

                 [--config [KEY=VALUE [KEY=VALUE ...]]] [--configfile FILE]

                 [--list] [--list-target-rules] [--directory DIR] [--dryrun]

                 [--printshellcmds] [--debug-dag] [--dag]

                 [--force-use-threads] [--rulegraph] [--d3dag] [--summary]

                 [--detailed-summary] [--archive FILE] [--touch]

                 [--keep-going] [--force] [--forceall]

                 [--forcerun [TARGET [TARGET ...]]]

                 [--prioritize TARGET [TARGET ...]]

                 [--until TARGET [TARGET ...]]

                 [--omit-from TARGET [TARGET ...]] [--allow-ambiguity]

                 [--cluster CMD | --cluster-sync CMD | --drmaa [ARGS]]

                 [--drmaa-log-dir DIR] [--cluster-config FILE]

                 [--immediate-submit] [--jobscript SCRIPT] [--jobname NAME]

                 [--cluster-status CLUSTER_STATUS] [--kubernetes [NAMESPACE]]

                 [--kubernetes-env ENVVAR [ENVVAR ...]]

                 [--container-image IMAGE] [--reason] [--stats FILE]

                 [--nocolor] [--quiet] [--nolock] [--unlock]

                 [--cleanup-metadata FILE [FILE ...]] [--rerun-incomplete]

                 [--ignore-incomplete] [--list-version-changes]

                 [--list-code-changes] [--list-input-changes]

                 [--list-params-changes] [--latency-wait SECONDS]

                 [--wait-for-files [FILE [FILE ...]]] [--benchmark-repeats N]

                 [--notemp] [--keep-remote] [--keep-target-files]

                 [--keep-shadow]

                 [--allowed-rules ALLOWED_RULES [ALLOWED_RULES ...]]

                 [--max-jobs-per-second MAX_JOBS_PER_SECOND]

                 [--max-status-checks-per-second MAX_STATUS_CHECKS_PER_SECOND]

                 [--restart-times RESTART_TIMES] [--attempt ATTEMPT]

                 [--timestamp] [--greediness GREEDINESS] [--no-hooks]

                 [--print-compilation]

                 [--overwrite-shellcmd OVERWRITE_SHELLCMD] [--verbose]

                 [--debug] [--runtime-profile FILE] [--mode {0,1,2}]

                 [--bash-completion] [--use-conda] [--conda-prefix DIR]

                 [--create-envs-only] [--use-singularity]

                 [--singularity-prefix DIR] [--singularity-args ARGS]

                 [--wrapper-prefix WRAPPER_PREFIX]

                 [--default-remote-provider {S3,GS,FTP,SFTP,S3Mocked,gfal,gridftp}]

                 [--default-remote-prefix DEFAULT_REMOTE_PREFIX]

                 [--no-shared-fs] [--version]

                 [target [target ...]]

 

Snakemake is a Python based language and execution environment for GNU Make-

like workflows.

 

positional arguments:

  target                Targets to build. May be rules or files.

 

optional arguments:

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

  --profile PROFILE     Name of profile to use for configuring Snakemake.

                        Snakemake will search for a corresponding folder in

                        /etc/xdg/xdg-ubuntu/snakemake and

                        /home/kazu/.config/snakemake. Alternatively, this can

                        be an absolute or relative path. The profile folder

                        has to contain a file 'config.yaml'. This file can be

                        used to set default values for command line options in

                        YAML format. For example, '--cluster qsub' becomes

                        'cluster: qsub' in the YAML file. Profiles can be

                        obtained from https://github.com/snakemake-profiles.

  --snakefile FILE, -s FILE

                        The workflow definition in a snakefile.

  --gui [PORT]          Serve an HTML based user interface to the given

                        network and port e.g. 168.129.10.15:8000. By default

                        Snakemake is only available in the local network

                        (default port: 8000). To make Snakemake listen to all

                        ip addresses add the special host address 0.0.0.0 to

                        the url (0.0.0.0:8000). This is important if Snakemake

                        is used in a virtualised environment like Docker. If

                        possible, a browser window is opened.

  --cores [N], --jobs [N], -j [N]

                        Use at most N cores in parallel (default: 1). If N is

                        omitted, the limit is set to the number of available

                        cores.

  --local-cores N       In cluster mode, use at most N cores of the host

                        machine in parallel (default: number of CPU cores of

                        the host). The cores are used to execute local rules.

                        This option is ignored when not in cluster mode.

  --resources [NAME=INT [NAME=INT ...]], --res [NAME=INT [NAME=INT ...]]

                        Define additional resources that shall constrain the

                        scheduling analogously to threads (see above). A

                        resource is defined as a name and an integer value.

                        E.g. --resources gpu=1. Rules can use resources by

                        defining the resource keyword, e.g. resources: gpu=1.

                        If now two rules require 1 of the resource 'gpu' they

                        won't be run in parallel by the scheduler.

  --config [KEY=VALUE [KEY=VALUE ...]], -C [KEY=VALUE [KEY=VALUE ...]]

                        Set or overwrite values in the workflow config object.

                        The workflow config object is accessible as variable

                        config inside the workflow. Default values can be set

                        by providing a JSON file (see Documentation).

  --configfile FILE     Specify or overwrite the config file of the workflow

                        (see the docs). Values specified in JSON or YAML

                        format are available in the global config dictionary

                        inside the workflow.

  --list, -l            Show availiable rules in given Snakefile.

  --list-target-rules, --lt

                        Show available target rules in given Snakefile.

  --directory DIR, -d DIR

                        Specify working directory (relative paths in the

                        snakefile will use this as their origin).

  --dryrun, -n          Do not execute anything.

  --printshellcmds, -p  Print out the shell commands that will be executed.

  --debug-dag           Print candidate and selected jobs (including their

                        wildcards) while inferring DAG. This can help to debug

                        unexpected DAG topology or errors.

  --dag                 Do not execute anything and print the directed acyclic

                        graph of jobs in the dot language. Recommended use on

                        Unix systems: snakemake --dag | dot | display

  --force-use-threads   Force threads rather than processes. Helpful if shared

                        memory (/dev/shm) is full or unavailable.

  --rulegraph           Do not execute anything and print the dependency graph

                        of rules in the dot language. This will be less

                        crowded than above DAG of jobs, but also show less

                        information. Note that each rule is displayed once,

                        hence the displayed graph will be cyclic if a rule

                        appears in several steps of the workflow. Use this if

                        above option leads to a DAG that is too large.

                        Recommended use on Unix systems: snakemake --rulegraph

                        | dot | display

  --d3dag               Print the DAG in D3.js compatible JSON format.

  --summary, -S         Print a summary of all files created by the workflow.

                        The has the following columns: filename, modification

                        time, rule version, status, plan. Thereby rule version

                        contains the versionthe file was created with (see the

                        version keyword of rules), and status denotes whether

                        the file is missing, its input files are newer or if

                        version or implementation of the rule changed since

                        file creation. Finally the last column denotes whether

                        the file will be updated or created during the next

                        workflow execution.

  --detailed-summary, -D

                        Print a summary of all files created by the workflow.

                        The has the following columns: filename, modification

                        time, rule version, input file(s), shell command,

                        status, plan. Thereby rule version contains the

                        versionthe file was created with (see the version

                        keyword of rules), and status denotes whether the file

                        is missing, its input files are newer or if version or

                        implementation of the rule changed since file

                        creation. The input file and shell command columns are

                        selfexplanatory. Finally the last column denotes

                        whether the file will be updated or created during the

                        next workflow execution.

  --archive FILE        Archive the workflow into the given tar archive FILE.

                        The archive will be created such that the workflow can

                        be re-executed on a vanilla system. The function needs

                        conda and git to be installed. It will archive every

                        file that is under git version control. Note that it

                        is best practice to have the Snakefile, config files,

                        and scripts under version control. Hence, they will be

                        included in the archive. Further, it will add input

                        files that are not generated by by the workflow itself

                        and conda environments. Note that symlinks are

                        dereferenced. Supported formats are .tar, .tar.gz,

                        .tar.bz2 and .tar.xz.

  --touch, -t           Touch output files (mark them up to date without

                        really changing them) instead of running their

                        commands. This is used to pretend that the rules were

                        executed, in order to fool future invocations of

                        snakemake. Fails if a file does not yet exist.

  --keep-going, -k      Go on with independent jobs if a job fails.

  --force, -f           Force the execution of the selected target or the

                        first rule regardless of already created output.

  --forceall, -F        Force the execution of the selected (or the first)

                        rule and all rules it is dependent on regardless of

                        already created output.

  --forcerun [TARGET [TARGET ...]], -R [TARGET [TARGET ...]]

                        Force the re-execution or creation of the given rules

                        or files. Use this option if you changed a rule and

                        want to have all its output in your workflow updated.

  --prioritize TARGET [TARGET ...], -P TARGET [TARGET ...]

                        Tell the scheduler to assign creation of given targets

                        (and all their dependencies) highest priority.

                        (EXPERIMENTAL)

  --until TARGET [TARGET ...], -U TARGET [TARGET ...]

                        Runs the pipeline until it reaches the specified rules

                        or files. Only runs jobs that are dependencies of the

                        specified rule or files, does not run sibling DAGs.

  --omit-from TARGET [TARGET ...], -O TARGET [TARGET ...]

                        Prevent the execution or creation of the given rules

                        or files as well as any rules or files that are

                        downstream of these targets in the DAG. Also runs jobs

                        in sibling DAGs that are independent of the rules or

                        files specified here.

  --allow-ambiguity, -a

                        Don't check for ambiguous rules and simply use the

                        first if several can produce the same file. This

                        allows the user to prioritize rules by their order in

                        the snakefile.

  --cluster CMD, -c CMD

                        Execute snakemake rules with the given submit command,

                        e.g. qsub. Snakemake compiles jobs into scripts that

                        are submitted to the cluster with the given command,

                        once all input files for a particular job are present.

                        The submit command can be decorated to make it aware

                        of certain job properties (input, output, params,

                        wildcards, log, threads and dependencies (see the

                        argument below)), e.g.: $ snakemake --cluster 'qsub

                        -pe threaded {threads}'.

  --cluster-sync CMD    cluster submission command will block, returning the

                        remote exitstatus upon remote termination (for

                        example, this should be usedif the cluster command is

                        'qsub -sync y' (SGE)

  --drmaa [ARGS]        Execute snakemake on a cluster accessed via DRMAA,

                        Snakemake compiles jobs into scripts that are

                        submitted to the cluster with the given command, once

                        all input files for a particular job are present. ARGS

                        can be used to specify options of the underlying

                        cluster system, thereby using the job properties

                        input, output, params, wildcards, log, threads and

                        dependencies, e.g.: --drmaa ' -pe threaded {threads}'.

                        Note that ARGS must be given in quotes and with a

                        leading whitespace.

  --drmaa-log-dir DIR   Specify a directory in which stdout and stderr files

                        of DRMAA jobs will be written. The value may be given

                        as a relative path, in which case Snakemake will use

                        the current invocation directory as the origin. If

                        given, this will override any given '-o' and/or '-e'

                        native specification. If not given, all DRMAA stdout

                        and stderr files are written to the current working

                        directory.

  --cluster-config FILE, -u FILE

                        A JSON or YAML file that defines the wildcards used in

                        'cluster'for specific rules, instead of having them

                        specified in the Snakefile. For example, for rule

                        'job' you may define: { 'job' : { 'time' : '24:00:00'

                        } } to specify the time for rule 'job'. You can

                        specify more than one file. The configuration files

                        are merged with later values overriding earlier ones.

  --immediate-submit, --is

                        Immediately submit all jobs to the cluster instead of

                        waiting for present input files. This will fail,

                        unless you make the cluster aware of job dependencies,

                        e.g. via: $ snakemake --cluster 'sbatch --dependency

                        {dependencies}. Assuming that your submit script (here

                        sbatch) outputs the generated job id to the first

                        stdout line, {dependencies} will be filled with space

                        separated job ids this job depends on.

  --jobscript SCRIPT, --js SCRIPT

                        Provide a custom job script for submission to the

                        cluster. The default script resides as 'jobscript.sh'

                        in the installation directory.

  --jobname NAME, --jn NAME

                        Provide a custom name for the jobscript that is

                        submitted to the cluster (see --cluster). NAME is

                        "snakejob.{rulename}.{jobid}.sh" per default. The

                        wildcard {jobid} has to be present in the name.

  --cluster-status CLUSTER_STATUS

                        Status command for cluster execution. This is only

                        considered in combination with the --cluster flag. If

                        provided, Snakemake will use the status command to

                        determine if a job has finished successfully or

                        failed. For this it is necessary that the submit

                        command provided to --cluster returns the cluster job

                        id. Then, the status command will be invoked with the

                        job id. Snakemake expects it to return 'success' if

                        the job was successfull, 'failed' if the job failed

                        and 'running' if the job still runs.

  --kubernetes [NAMESPACE]

                        Execute workflow in a kubernetes cluster (in the

                        cloud). NAMESPACE is the namespace you want to use for

                        your job (if nothing specified: 'default'). Usually,

                        this requires --default-remote-provider and --default-

                        remote-prefix to be set to a S3 or GS bucket where

                        your . data shall be stored. It is further advisable

                        to activate conda integration via --use-conda.

  --kubernetes-env ENVVAR [ENVVAR ...]

                        Specify environment variables to pass to the

                        kubernetes job.

  --container-image IMAGE

                        Docker image to use, e.g., when submitting jobs to

                        kubernetes. By default, this is

                        'quay.io/snakemake/snakemake', tagged with the same

                        version as the currently running Snakemake instance.

                        Note that overwriting this value is up to your

                        responsibility. Any used image has to contain a

                        working snakemake installation that is compatible with

                        (or ideally the same as) the currently running

                        version.

  --reason, -r          Print the reason for each executed rule.

  --stats FILE          Write stats about Snakefile execution in JSON format

                        to the given file.

  --nocolor             Do not use a colored output.

  --quiet, -q           Do not output any progress or rule information.

  --nolock              Do not lock the working directory

  --unlock              Remove a lock on the working directory.

  --cleanup-metadata FILE [FILE ...], --cm FILE [FILE ...]

                        Cleanup the metadata of given files. That means that

                        snakemake removes any tracked version info, and any

                        marks that files are incomplete.

  --rerun-incomplete, --ri

                        Re-run all jobs the output of which is recognized as

                        incomplete.

  --ignore-incomplete, --ii

                        Do not check for incomplete output files.

  --list-version-changes, --lv

                        List all output files that have been created with a

                        different version (as determined by the version

                        keyword).

  --list-code-changes, --lc

                        List all output files for which the rule body (run or

                        shell) have changed in the Snakefile.

  --list-input-changes, --li

                        List all output files for which the defined input

                        files have changed in the Snakefile (e.g. new input

                        files were added in the rule definition or files were

                        renamed). For listing input file modification in the

                        filesystem, use --summary.

  --list-params-changes, --lp

                        List all output files for which the defined params

                        have changed in the Snakefile.

  --latency-wait SECONDS, --output-wait SECONDS, -w SECONDS

                        Wait given seconds if an output file of a job is not

                        present after the job finished. This helps if your

                        filesystem suffers from latency (default 5).

  --wait-for-files [FILE [FILE ...]]

                        Wait --latency-wait seconds for these files to be

                        present before executing the workflow. This option is

                        used internally to handle filesystem latency in

                        cluster environments.

  --benchmark-repeats N

                        Repeat a job N times if marked for benchmarking

                        (default 1).

  --notemp, --nt        Ignore temp() declarations. This is useful when

                        running only a part of the workflow, since temp()

                        would lead to deletion of probably needed files by

                        other parts of the workflow.

  --keep-remote         Keep local copies of remote input files.

  --keep-target-files   Do not adjust the paths of given target files relative

                        to the working directory.

  --keep-shadow         Do not delete the shadow directory on snakemake

                        startup.

  --allowed-rules ALLOWED_RULES [ALLOWED_RULES ...]

                        Only consider given rules. If omitted, all rules in

                        Snakefile are used. Note that this is intended

                        primarily for internal use and may lead to unexpected

                        results otherwise.

  --max-jobs-per-second MAX_JOBS_PER_SECOND

                        Maximal number of cluster/drmaa jobs per second,

                        default is 10, fractions allowed.

  --max-status-checks-per-second MAX_STATUS_CHECKS_PER_SECOND

                        Maximal number of job status checks per second,

                        default is 10, fractions allowed.

  --restart-times RESTART_TIMES

                        Number of times to restart failing jobs (defaults to

                        0).

  --attempt ATTEMPT     Internal use only: define the initial value of the

                        attempt parameter (default: 1).

  --timestamp, -T       Add a timestamp to all logging output

  --greediness GREEDINESS

                        Set the greediness of scheduling. This value between 0

                        and 1 determines how careful jobs are selected for

                        execution. The default value (1.0) provides the best

                        speed and still acceptable scheduling quality.

  --no-hooks            Do not invoke onstart, onsuccess or onerror hooks

                        after execution.

  --print-compilation   Print the python representation of the workflow.

  --overwrite-shellcmd OVERWRITE_SHELLCMD

                        Provide a shell command that shall be executed instead

                        of those given in the workflow. This is for debugging

                        purposes only.

  --verbose             Print debugging output.

  --debug               Allow to debug rules with e.g. PDB. This flag allows

                        to set breakpoints in run blocks.

  --runtime-profile FILE

                        Profile Snakemake and write the output to FILE. This

                        requires yappi to be installed.

  --mode {0,1,2}        Set execution mode of Snakemake (internal use only).

  --bash-completion     Output code to register bash completion for snakemake.

                        Put the following in your .bashrc (including the

                        accents): `snakemake --bash-completion` or issue it in

                        an open terminal session.

  --use-conda           If defined in the rule, run job in a conda

                        environment. If this flag is not set, the conda

                        directive is ignored.

  --conda-prefix DIR    Specify a directory in which the 'conda' and 'conda-

                        archive' directories are created. These are used to

                        store conda environments and their archives,

                        respectively. If not supplied, the value is set to the

                        '.snakemake' directory relative to the invocation

                        directory. If supplied, the `--use-conda` flag must

                        also be set. The value may be given as a relative

                        path, which will be extrapolated to the invocation

                        directory, or as an absolute path.

  --create-envs-only    If specified, only creates the job-specific conda

                        environments then exits. The `--use-conda` flag must

                        also be set.

  --use-singularity     If defined in the rule, run job within a singularity

                        container. If this flag is not set, the singularity

                        directive is ignored.

  --singularity-prefix DIR

                        Specify a directory in which singularity images will

                        be stored.If not supplied, the value is set to the

                        '.snakemake' directory relative to the invocation

                        directory. If supplied, the `--use-singularity` flag

                        must also be set. The value may be given as a relative

                        path, which will be extrapolated to the invocation

                        directory, or as an absolute path.

  --singularity-args ARGS

                        Pass additional args to singularity.

  --wrapper-prefix WRAPPER_PREFIX

                        Prefix for URL created from wrapper directive

                        (default: https://bitbucket.org/snakemake/snakemake-

                        wrappers/raw/). Set this to a different URL to use

                        your fork or a local clone of the repository.

  --default-remote-provider {S3,GS,FTP,SFTP,S3Mocked,gfal,gridftp}

                        Specify default remote provider to be used for all

                        input and output files that don't yet specify one.

  --default-remote-prefix DEFAULT_REMOTE_PREFIX

                        Specify prefix for default remote provider. E.g. a

                        bucket name.

  --no-shared-fs        Do not assume that jobs share a common file system.

                        When this flag is activated, Snakemake will assume

                        that the filesystem on a cluster node is not shared

                        with other nodes. For example, this will lead to

                        downloading remote files on each cluster node

                        separately. Further, it won't take special measures to

                        deal with filesystem latency issues. This option will

                        in most cases only make sense in combination with

                        --default-remote-provider. Further, when using

                        --cluster you will have to also provide --cluster-

                        status. Only activate this if you know what you are

                        doing.

  --version, -v         show program's version number and exit

 

実行方法

ランにはデータベースディレクトリやシーケンシングデータのパスを記載した、以下のようなconfigファイルが必要になる。

workdir: "/home/user/folder/results/"
dbdir: "/home/user/folder/databases/"
samples:
  sample_name_1:
     fq1: "/home/user/folder/reads/file.1.fq"
     fq2: "/home/user/folder/reads/file.2.fq"

データベースはmetametaのGithubレポジトリのlinkからダウンロードできる。

f:id:kazumaxneo:20180920212550j:plain

 

 configファイルを全て埋めたら実行する。

metameta --configfile yourconfig.yaml --use-conda --keep-going --cores 24

 

 テストラン

#metametaをインストールしたディレクトリに移動する。pyenv下でanacondaを管理していたので以下の場所に移動した。
cd /home/kazu/.pyenv/versions/anaconda3-4.0.0/opt/metameta/
metameta --configfile sampledata/sample_data_custom_viral.yaml --use-conda --keep-going --cores 6

 大量のエラーが出る。ツール間のバッティングがありそうなので、condaで仮想環境を作りやり直してみる。

 

 

 

 

 

複数結果の統合、さらなるツールの追加などもサポートしています。

引用

MetaMeta: integrating metagenome analysis tools to improve taxonomic profiling
Piro VC, Matschkowski M, Renard BY

Microbiome. 2017 Aug 14;5(1):101