macでインフォマティクス

macでインフォマティクス

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

複数のアセンブラとk-merを使ったTranscriptome 自動アセンブリワークフロー Oyster River Protocol

 2018 11/2 コマンド追記 & 誤字修正

 2018 11/7 誤字修正

2019 4/6 docker追記

2019 6/17 追記、誤字修正

2019 6/21追記

2019 7/5 Step by step instructions link追記

 

 現代のシーケンシング技術は細胞内の代謝過程から人口変動パターンまで、非常に幅広い自然現象の基礎となるゲノムレベルのプロセスを深く理解する機会を提供してきた。トランスクリプトームシーケンシングは、特に functional genomics(論文より Lappalainen et al、2013; Cahoy et al、2008)に影響を与えており、何年か前には不可能だった発見も可能になってきている(Mortazavi et al、2008; Wang、Gerstein&Snyder、2009)。この理由の大部分は、これらの研究が実施されるスケールに起因している(Li et al、2017; Tan et al、2017)。 1つまたは少数の候補遺伝子に基づくadaptationの研究とは異なり(Fitzpatrick et al、2005; Panhuis、2006)、現代の研究は発現された転写産物を同時に全評価(トランスクリプトーム)することができる。ダイナミックレンジの高まりの直接的な結果として、スケール問題に加えて、低および高発現転写産物を同時に再構成および定量する能力が高まった(Wolf、2013; Vijay et al、2013)。最後に、遺伝子発現の変化を検出するためのより改善された方法が登場し(Robinson、McCarthy&Smyth、2010 (link); Love、Huber&Anders、2014 (link))、遺伝子発現変動の理解が目的の研究の解像度が向上した。

 広範囲にわたり人気になった直接的な結果として、トランスクリプトームのアセンブリのための多様なツールセットが存在しているが、潜在的には、再構成された転写産物のそれぞれが再構築に失敗する。最も初期の特別なトランスクリプトームアセンブラには、Trans-ABySS(Robertson et al、2010)、Oases(Schulz et al、2012)、SOAPdenovoTrans(Xie et al、2014)のパッケージがある。これらは、 de BruijnのgraphベースのゲノムアセンブラのABySS(Simpson et al、2009)、Velvet(Zerbino&Birney、2008)、SOAP(Li et al、2008)に機能が基づいている。これらの初期の努力から、トリニティ(Haas et al、2013)およびIDBA-Tran(Pengら、2013)のような一連のより専門化されたde novoトランスクリプトームアセンブラが生み出された。 De Bruijnのグラフアプローチは強力だが、より新しく開発されたソフトウェアはアルゴリズム的な視点で斬新な部分を探究しており、新しいメソッドがトランスクリプトームの再構成できなかった配列をアセンブリできると仮定すれば、大きな利点を提供できている。例えば、BinPacker(Liu et al、2016)(紹介)は、古典的なビン充填問題の後にアセンブリ問題をモデル化するde Bruijn graph手法を放棄し、一方、Shannon(Kannan et al、2016 preprint)は、ソフトウェアエンジニアが決定したヒューリスティックではなくinformation theoryを使う。これらの新しいアセンブラは、基本的に異なるアセンブリアルゴリズムを実装することにより、他のアセンブラが正確にアセンブルできないトランスクリプトームの部分を再構築する可能性がある。

 デノボアセンブリに利用可能な様々なツールに加えて、リードトリミングによる前処理(例えば、Skewer; Jiang et al、2014、Trimmomatic、Bolger、Lohse&Usadel、2014 、2011)、リードの正規化(khmer; Pell et al、2012)、リードのエラー訂正(SEECER; Le et al、2013、RCorrector; Song&Florea、2015、Reptile; Yang、Dorman&Aluru 、2010)がある。同様に、TransRate(Smith-Unna et al、2016)、BUSCO(Benchmarking Universal Single-Copy Orthologs、Simãoet al、2015)、Detonate(Li et al、2014)などのアセンブリされたトランスクリプトのクオリティを評価するベンチマークツールも同様に 開発されている。

(一部略)

 トランスクリプトーム再構成に関連する微妙な(そしてそれほど微妙ではない)方法論的課題は、非常に可変なアセンブリ品質をもたらす可能性がある。特に、ほとんどのツールはデフォルト設定を使用して実行されるが、これらのデフォルトは、特定の(多くの場合、不特定の)ユースケースまたはデータ型に対してのみ有効である。パラメータの最適化は本質的にデータセット依存型と階乗型であるため、特にパイプライン全体の網羅的な最適化は不可能である。このことを考慮すると、デノボトランスクリプトームアセンブリを実行するには、時間とリソース面で大きな投資が必要であり、各ステップを慎重に検討する必要がある。ここでは、様々な一般的な実験条件または分類群にわたって、高品質のトランスクリプトームアセンブリ作製をもたらす、エビデンスに基づくアセンブリプロトコルを提案する。

この論文では、トランスクリプトームアセンブリ用のオイスターリバープロトコル(ORP)1の開発について説明する。これは、Vijay et al (2013) に記載されている多くの欠点を明示的に考慮し、その改善を試みる。 本方法は複数のkmerとマルチアセンブラ戦略を活用している。この革新は非常に重要である。なぜなら、すべてのアセンブリソリューションがバイアスがある方法でシーケンシングリードデータを処理し転写産物をレカバリするからである。具体的には、アセンブリソフトウェアの開発では、アセンブリ問題自体の範囲を考慮して必要なヒューリスティックセットを使用する。各ソフトウェア開発チームが様々なアセンブリアルゴリズムを実装しながらこれらのヒューリスティックに関連するユニークなアイデアを持ち、個々のアセンブラはユニークなアセンブリ動作を示す。マルチアセンブラのアプローチを活用することで、あるアセンブラの強みが別のアセンブラの弱みを補うことができる。アセンブリヒューリスティックに関連するバイアスに加えて、アセンブリのkmer-lengthは転写再構成に重要な影響を及ぼす。より短いkmersは、abundantな転写産物より少ない転写産物を効率的に再構築することがよく知られている。この場合、複数の異なるkmer長でアセンブルし、結果として得られるアセンブリをマージすると、このタイプのバイアスを効果的に減らすことができる。これらの問題を認識して、私 (著者のこと) は、複数の異なるアセンブラの組み合わせと複数のアセンブリ - kmers長によって生じるアセンブリが、さまざまなメトリックにわたって、個々のアセンブリより優れていると仮定している。

強化されたパイプラインを開発することに加えて、この作業では、新しいアセンブリアルゴリズムやパイプラインを開発する際に他の研究者が使用できる、完全にベンチマークされたリファレンスアセンブリを利用できるようにしながら、他の多くの研究者がアセンブリ方法の比較を公表してきたが、これまでは単一のデータセットがいくつかの異なる方法でアセンブリされていた(Marchant et al、2016; Finseth&Harrison、2014)。

 

ソフトウエア

ORPはLinuxプラットフォームにインストールできる。Linuxbrew(Jackman&Birol、2016)のみ必要で、インストールにスーパーユーザー権限は必要ない。 このソフトウェアは、以下に説明するすべてのステップを調整するスタンドアローンmakefileとして実装されている。ソフトウェアはバージョン管理されており、共有と再利用を促進するためにオープンライセンスとなっている。 ユーザーガイドはhttp://oyster-river-protocol.rtfd.ioから入手できる。 

1、プリアセンブリ

Illuminaシーケンシングアダプターを、MacManes(2014 link)の推奨に従い、Trimmomaticバージョン0.36(Bolger、Lohse&Usadel、2014)プログラムを用いて除去した(リードの両端からPhred≦2を除去 *1)。トリミング後、MacManes&Eisen(2013)の推奨に従ってソフトウェアRCorrectorバージョン1.0.2(Song&Florea、2015)(紹介)を使用してエラー訂正した。

2、アセンブリ

トリムしてエラーコレクションされたデータセットを使い、3つの異なるトランスクリプトームアセンブラと3つの異なるkmer長で4つのユニークなアセンブリを作成した。まず、リードの正規化を行わずに、Trinityリリース2.4.0(Haas et al、2013)とデフォルト設定(k = 25)を使用してアセンブリを行なった。正規化しない決定は、正規化されたデータセットのパフォーマンスがわずかに悪化していることを示す以前の研究(MacManes、2015 preprint)に基づいている。次に、SPAdes RNAseqアセンブラ(バージョン3.10)(Chikhi&Medvedev、2014)(SPAdes紹介)によるアセンブリを、k-merサイズ55および75の別個の2つのランで行った。最後に、Shannonアセンブラのバージョン0.0.2(Kannan et al)によるアセンブリをk-merサイズ75で行った。これらのアセンブラは以下の事実に基づいて選択された:(1)オープンサイエンス開発モデルを使っており、エンドユーザがコードに貢献する可能性がある。(2)維持されており、積極的な開発が続いている。(3)異なるアルゴリズムを使っている。

3、OrthoFuseを使ったアセンブリのマージ

4つのアセンブリを結合するため、私(著者)はトランスクリプトームアセンブリを効果的にマージする新しいソフトウェアを開発した。簡単に説明すると、OrthoFuseはすべてのアセンブリをコンカテネートすることから始まり、ORPにパッケージ化されている、ヌクレオチド配列を受け入れるようにmodifyされてORPに組み込まれたOrthoFinder(Emms&Kelly、2015)を実行して転写産物グループを形成する。

一部略

Orthofinderの実行後、modifyされてORPに組み込まれたTransRateバージョン1.0.3(Smith-Unna et al、2016)(紹介)をマージされたアセンブリに使い、その後最高のコンティグスコアのcontigを各グループから選択し、新しいアセンブリファイルに配置する。得られたファイルは、すべての下流分析に使用できる。 

4、アセンブリの評価

すべてのアセンブリは、ORP-TransRate、Detonateバージョン1.11(Li et al、2014)(紹介)、shmlastバージョン1.2(Scott、2017, link)、およびBUSCOバージョン3.0.2(Simãoet al、2015)(紹介)を使用して評価された。 TransRateは、長さとマッピングベースのmetrcisによりスコアを生成し、トランスクリプトームアセンブリの連続性を評価する。Detonateはorthogonal analysisを実行し、アセンブリのスコアを生成する。 BUSCOは、すべての真核生物に見られる保存されたシングルコピーオーソログのアセンブリを検索することにより、アセンブリを評価する。具体的には、「完全なオーソログ」は、BUSCOグループ平均の長さの2標準偏差以内の転写産物として定義され、このメトリックを満たすと"complete orthologs"、満たさないと"fragmented"と判定される。

 

 

Oyster River Protocolに関するツイート

 

インストール

mac os10.12でテストした。

依存

  • Rcorrector
  • HMMER
  • Trimmomatic
  • Trinity
  • SPAdes
  • TransABySS
  • MCL
  • Metis
  • OrthoFuser
  • BLAST
  • seqtk
  • BUSCO (make sure to install databases)
  • TransRate (the ORP version packaged here)
  • Python modules numpy, scipy, biopython, cvxopt

本体 Github

 ツールの導入過程でAnacondaや様々な環境が導入される。ここでは仮想環境で実行した。

cd ~
git clone https://github.com/macmanes-lab/Oyster_River_Protocol.git
cd Oyster_River_Protocol
make
source ~/.profile

#configファイルを修正する
vi software/config.ini

> cat software/config.ini 

$ cat software/config.ini 

[busco]

lineage_path = /home/kazu/Oyster_River_Protocol/busco_dbs/eukaryota_odb9

force = True

quiet = False

 

[tblastn]

path = /home/kazu/Oyster_River_Protocol/software/anaconda/install/envs/orp_v2/bin/

 

[makeblastdb]

path = /home/kazu/Oyster_River_Protocol/software/anaconda/install/envs/orp_v2/bin/

 

[hmmsearch]

path = /home/kazu/Oyster_River_Protocol/software/anaconda/install/envs/orp_v2/bin/

パスを調べ、ここでは↑ のように修正した。

 

仮想環境"orp_v2"を立ち上げる。

source $HOME/Oyster_River_Protocol/software/anaconda/install/bin/activate orp_v2

ヘルプが表示されるか確認。フルパスで指定する。

>  ~/Oyster_River_Protocol/oyster.mk -h

$ ~/Oyster_River_Protocol/oyster.mk -h

 

 

*****  Welcome to the Oyster River Prptocol ***** 

*****  This is version 2.1.0 *****

 

Usage:

 

/path/to/Oyster_River/Protocol/oyster.mk main CPU=24 

MEM=128

STRAND=RF

READ1=1.subsamp_1.cor.fq

READ2=1.subsamp_2.cor.fq

RUNOUT=test

 

(orp_v2) kazu@aa796a2b1908:~/Oyster_River_Protocol$ /home/kazu/Oyster_River_Protocol/oyster.mk -h

Usage: make [options] [target] ...

Options:

  -b, -m                      Ignored for compatibility.

  -B, --always-make           Unconditionally make all targets.

  -C DIRECTORY, --directory=DIRECTORY

                              Change to DIRECTORY before doing anything.

  -d                          Print lots of debugging information.

  --debug[=FLAGS]             Print various types of debugging information.

  -e, --environment-overrides

                              Environment variables override makefiles.

  -f FILE, --file=FILE, --makefile=FILE

                              Read FILE as a makefile.

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

  -i, --ignore-errors         Ignore errors from commands.

  -I DIRECTORY, --include-dir=DIRECTORY

                              Search DIRECTORY for included makefiles.

  -j [N], --jobs[=N]          Allow N jobs at once; infinite jobs with no arg.

  -k, --keep-going            Keep going when some targets can't be made.

  -l [N], --load-average[=N], --max-load[=N]

                              Don't start multiple jobs unless load is below N.

  -L, --check-symlink-times   Use the latest mtime between symlinks and target.

  -n, --just-print, --dry-run, --recon

                              Don't actually run any commands; just print them.

  -o FILE, --old-file=FILE, --assume-old=FILE

                              Consider FILE to be very old and don't remake it.

  -p, --print-data-base       Print make's internal database.

  -q, --question              Run no commands; exit status says if up to date.

  -r, --no-builtin-rules      Disable the built-in implicit rules.

  -R, --no-builtin-variables  Disable the built-in variable settings.

  -s, --silent, --quiet       Don't echo commands.

  -S, --no-keep-going, --stop

                              Turns off -k.

  -t, --touch                 Touch targets instead of remaking them.

  -v, --version               Print the version number of make and exit.

  -w, --print-directory       Print the current directory.

  --no-print-directory        Turn off -w, even if it was turned on implicitly.

  -W FILE, --what-if=FILE, --new-file=FILE, --assume-new=FILE

                              Consider FILE to be infinitely new.

  --warn-undefined-variables  Warn when an undefined variable is referenced.

 

This program built for x86_64-pc-linux-gnu

Report bugs to <bug-make@gnu.org>

 

インストールはMakefileに従い自動に行われますが、ビルドが面倒で数回失敗したので、dockerイメージも上げておきます。testデータは捨ててあるので、このイメージで下記のテストランを行う場合はgithubから再取得して下さい。このイメージは16GBほどあります。

docker pull kazumax/oyster_river_protocol

#run
docker run -itv $PWD:/data kazumax/oyster_river_protocol
#user: kazu
#pass: 1898
su kazu
cd ~
source ~/.profiles
$HOME/Oyster_River_Protocol/software/anaconda/install/bin/activate orp_v2

#help
/home/kazu/Oyster_River_Protocol/oyster.mk -h

↑ 非推奨

2019 4/6 追記

オーサーがdockerイメージを準備してくれています。こちらを使って下さい。

How to install the ORP using Docker — Oyster River Protocol latest documentation

#docker hub version

#latest tagが付いてないのでversionチェックして入れてください
docker pull macmaneslab/orp:2.2.7

cd RNA_seq_raw_fastq
docker run -itv $PWD:/data -w /data macmaneslab/orp:2.2.7
> conda activate orp
> $HOME/Oyster_River_Protocol/oyster.mk \
TPM_FILT=1 \
MEM=50 \
CPU=40 \
READ1=pair_1.fq.gz \
READ2=pair_1.fq.gz \
RUNOUT=output

 

 

テストラン

Run test dataset 

cd sampledata
#実行(*2)
$HOME/Oyster_River_Protocol/oyster.mk main STRAND=RF \
MEM=15 \
CPU=8 \
READ1=test.1.fq.gz \
READ2=test.2.fq.gz \
RUNOUT=test

seqs] Processed 53898 reads in 1.990 CPU sec, 0.424 real sec

[main] Version: 0.7.17-r1188

[main] CMD: bwa mem -t 8 test /dev/fd/63 /dev/fd/62

[main] Real time: 0.821 sec; CPU: 2.110 sec

[bam_sort_core] merging from 0 files and 10 in-memory blocks...

-parsing file: test.sorted.bam

-done parsing file, examining orientations of reads.

 

 7| ##         

 6| ##         

 5| ##         

 4| ##         

 3| ##         

 2| ###       #

 1| ###  ##   #

   -----------

 

------------------------------------

|             Summary              |

------------------------------------

|         observations: 20         |

|       min value: -1.000000       |

|         mean : -0.988000         |

|       max value: -0.935000       |

------------------------------------

 

 

*****  See the following link for interpretation ***** 

*****  https://oyster-river-protocol.readthedocs.io/en/latest/strandexamine.html ***** 

 

 

 

*****  QUALITY REPORT FOR: test using the ORP version 2.1.0 ****

*****  THE ASSEMBLY CAN BE FOUND HERE: /home/kazu/Oyster_River_Protocol/sampledata/assemblies/test.ORP.fasta **** 

 

*****  BUSCO SCORE ~~~~~>           C:0.0%[S:0.0%,D:0.0%],F:0.3%,M:99.7%,n:303

*****  TRANSRATE SCORE ~~~~~>           0.37907

*****  TRANSRATE OPTIMAL SCORE ~~~~~>   0.53726

*****  UNIQUE GENES ORP ~~~~~>          37

*****  UNIQUE GENES TRINITY ~~~~~>      31

*****  UNIQUE GENES SPADES55 ~~~~~>     22

*****  UNIQUE GENES SPADES75 ~~~~~>     23

*****  UNIQUE GENES TRANSABYSS ~~~~~>   35

ジョブ終了

QCリード 、アセンブリ、レポート等が出力される。
出力

f:id:kazumaxneo:20181102102919j:plain

 

レポート解析だけ実行する。 

report.mk

$HOME/Oyster_River_Protocol/report.mk main \
ASSEMBLY=assembly.fasta \
CPU=24 \
LINEAGE=eukaryota_odb9 \
READ1=SRR2016923_1.fastq \
READ2=SRR2016923_2.fastq \
RUNOUT=SRR2016923

 

 他にもquant.mk、strandeval.mkなどのコマンドがあります(version2.1)。helpで使い方は確認してください。また、HPのversion logも定期的に確認してください。一部アセンブラやツールが、論文と変更されたりしています。今後もパフォーマンス改善のチューニングがされるかもしれません。

 

追記

Step by step instructions

orp_website/full_instructions.md at master · macmanes-lab/orp_website · GitHub

 

引用
The Oyster River Protocol: a multi-assembler and kmer approach for de novo transcriptome assembly
Matthew D. MacManes

PeerJ. 2018; 6: e5428

 

参考

*1

Frontiers | On the optimal trimming of high-throughput mRNA sequence data | Genetics

では、phread quality trimingなし、≤2、≤5、≤20のトリミング、がその後の処理に与える影響を調べている。

  

複数アセンブラを使ったde novo assemblyのワークフローの実際の流れに興味があれば、以下のprotocol.ioページも確認してみて下さい。丁寧に説明されています。

 

*2

FR、RF等はBiostarsの質問を参照

https://www.biostars.org/p/169942/

Paired reads:

  • RF: first read (/1) of fragment pair is sequenced as anti-sense (reverse(R)), and second read (/2) is in the sense strand (forward(F)); typical of the dUTP/UDG sequencing method.
  • FR: first read (/1) of fragment pair is sequenced as sense (forward), and second read (/2) is in the antisense strand (reverse)

Unpaired (single) reads:

  • F: the single read is in the sense (forward) orientation
  • R: the single read is in the antisense (reverse) orientation

 

追記

良い論文が出ています。

De novo transcriptome assembly: A comprehensive cross-species comparison of short-read RNA-Seq assemblers | GigaScience | Oxford Academic