Trinotateを使ってde novo transcriptome のアセンブリ配列にアノテーションをつける

2020 9/16 pythonのバージョンを指定して導入, configファイルについて追記

2020 9/16 コメント追加, signalPとTMHMM、rnammerの初期設定追加

2020 9/27 わかりにくい表現を修正

Trinotarteは、Trinityののde novoトランスクリプトームアセンブリを機能的にアノテーションするためのプロトコルおよび支援ソフトウェア。

以前紹介したTrinotarteの説明が分かりづらかったので、簡潔にまとめ直します。

TrinotateはデフォルトではTrinityのアセンブリを使ってアノテーション付けを行います。そのため、別のde novoアセンブラを使った場合は少しだけ工夫が必要になります。下に書いた手順を確認して下さい、

wiki

TrinotateWeb: Graphical Interface for Navigating Trinotate Annotations and Expression Analyses

https://github.com/Trinotate/Trinotate.github.io/wiki/TrinotateWeb

インストール

macだと一部のバイナリ配布プログラムが動作しないのでlinuxを使う。ubuntu18.04LTSにて、condaでpython3.7の仮想環境を作ってテストした。

本体　Github 

#bioconda(link) 依存が多いツールは仮想環境を作って導入するのが無難
conda create -n trinotate python=3.7 -y
conda activate trinotate
conda install -c bioconda -y trinotate

#transdecoderは別途インストールが必要
conda install -c bioconda -y transdecoder

#最近のバージョンではhmmerとblastは入っているはずだが、無いなら導入
conda install -c bioconda -y hmmer
conda install -c bioconda -y blast

Trinotateのヘルプ

> autoTrinotate.pl

$ autoTrinotate.pl

##############################################################################

#

# Required:

#

#  --Trinotate_sqlite <string>                Trinotate.sqlite boilerplate database

#

#  --transcripts <string>                     transcripts.fasta

#

#  --gene_to_trans_map <string>               gene-to-transcript mapping file

#

#  --conf <string>                            config file

#

#  --CPU <int>                                number of threads to use.

#

#

##############################################################################

 RNAMMER、signalP、TMHMMはcondaでは入らない。それぞれダウンロードしてパスを通す。

まずHPからRNAMMER、signalP、TMHMMのプログラムをダウンロードする（ダウンロードを選ぶと記載したメールアドレスにダウンロードリンクが届く）。signalPはv５が出ているが、v4.1 の方をダウンロードする。

RNAMMER、signalP、TMHMM

https://services.healthtech.dtu.dk/software.php

signalPはダウンロード後初期設定が必要。signalp-4.1/signalPを開いてGENERAL SETTINGSのパスを修正する。ここでは/data/signalp-4.1/に実行スクリプトsignalpがあるので、以下のようにした。

FASTA.pmが見えないと言われたので
export PERL5LIB=/data/signalp-4.1/lib/で対処した。

TMHMMも初期設定が必要。中に入っているREADMEと TMHMM2.0.htmlに書いてあるが、tmhmmformat.plとtmhmmの１行目のperlパスを環境に合わせて修正する。また、decodeanhmm.Linux〜という実行形式ファイルが複数入っているので、動作するバージョンをdecodeanhmmにリネームする。ここではdecodeanhmm.Linux_x86_64をコピーしてリネーム。

> cp decodeanhmm.Linux_x86_64 decodeanhmm

rnammerも初期設定が必要。中に入っているrnammer-1.2.readmeに書いてあるが、$INSTALL_PATHのパスを環境に合わせて修正する。

他にも色々入れる必要がある。マニュアルの４ページ目によると、

http://arcsb2017.weebly.com/uploads/4/6/1/9/46199061/trinotate_install.pdf

hmmer2 を入れる（link）。

cd hmmer-2.3 
./configure 
make
make check
make install

rnammerのhmmsearchパスを修正する。導入したhmmsearch2のパスを記載する。

環境変数$PATHのhmmsearchもhmmsearch2シンボリックリンクに書き換えた。

mv /root/miniconda3/envs/trinotate/bin/hmmsearch /root/miniconda3/envs/trinotate/bin/hmmsearch3

ln -s /data/hmmer-2.3/src/hmmsearch /root/miniconda3/envs/trinotate/bin/hmmsearch

core-rnammerの”--cpu 1”を消す。２箇所ある。

rnammerのテストランはしたほうがいい。

cd rnammer-1.2.src/example/
rnammer -S bac -m lsu,ssu,tsu -gff - < ecoli.fsa 

出力

##gff-version2

##source-version RNAmmer-1.2

##date 2020-09-17

##Type DNA

# seqname           source                      feature     start      end   score   +/-  frame  attribute

# ---------------------------------------------------------------------------------------------------------

ecoli_section RNAmmer-1.2 rRNA 21067 21181 86.3 + . 5s_rRNA

ecoli_section RNAmmer-1.2 rRNA 16177 17706 1950.6 + . 16s_rRNA

ecoli_section RNAmmer-1.2 rRNA 18068 20969 3754.0 + . 23s_rRNA

# ---------------------------------------------------------------------------------------------------------

O.K！ 

もしXML::Simpleが無いなら、perl lib/にXML::Simpleをダウンロードして以下コマンドでビルドする。cpanmで入れると、gccのバージョン指定ミスなのかコンパイル時のチェックでこける。

https://metacpan.org/pod/XML::Simple

cd XML-Simple-2.25/
perl Makefile.PL
make
make install

 エラーがないか改めてテストする。

rnammer -S bac -m lsu,ssu,tsu -gff - < rnammer-1.2.src/example/

#自分のデータでもやったほうがいい。
rnammer -S bac -m lsu,ssu,tsu -gff - < tour_RNA_aasembly.fasta

 出力される.gffでrRNAが検出されているか確認する。ないとは思うが、TrinotateはrRNAがないデータだと、とrnammerのステップでエラー扱いになって止まる。

 データベースの準備

１、sqliteデータベースその他をダウンロードする。

Build_Trinotate_Boilerplate_SQLite_db.pl Trinotate

ダウンロードされた。

２、ダウンロードされたuniprot_sprot.pepを元に blastデータベースを作成する。

makeblastdb -in uniprot_sprot.pep -dbtype prot

 ３、続いてpfamデータベースを作成する。

#解凍
gunzip Pfam-A.hmm.gz

#index作成
hmmpress Pfam-A.hmm

レポジトリのスクリプトが必要。また、configファイル のテンプレートが置いてあるので、これをベースに修正していく。

git clone https://github.com/Trinotate/Trinotate.git
cp Trinotate/auto/conf.txt .

conf.txt

/dataにダウンロードしたプログラムsignalp-4.1、tmhmm-2.0c、rnammer-1.2.src、trinityrnaseqのレポジトリ、上のデータベースを配置して、以下のように修正した。condaで入れたツールはcondaのパスを記載。

実行方法

１、Trinityアセンブリ配列以外にgene-to-transcript mapping fileが必要。Trinityrnaseqのユーティリティで作成する（下記参照）。

git clone https://github.com/trinityrnaseq/trinityrnaseq.git

#trinityのアセンブリ結果のtrinity.fastaを指定する
perl trinityrnaseq/util/support_scripts/get_Trinity_gene_to_trans_map.pl trinity.fasta > gene_trans_map

> head gene_trans_map

$ head gene_trans_map

TRINITY_DN144_c0_g1 TRINITY_DN144_c0_g1_i1

TRINITY_DN179_c0_g1 TRINITY_DN179_c0_g1_i1

TRINITY_DN159_c0_g1 TRINITY_DN159_c0_g1_i1

TRINITY_DN159_c0_g2 TRINITY_DN159_c0_g2_i1

TRINITY_DN153_c0_g1 TRINITY_DN153_c0_g1_i1

TRINITY_DN153_c0_g2 TRINITY_DN153_c0_g2_i1

TRINITY_DN130_c0_g1 TRINITY_DN130_c0_g1_i1

TRINITY_DN106_c0_g1 TRINITY_DN106_c0_g1_i1

TRINITY_DN106_c0_g2 TRINITY_DN106_c0_g2_i1

TRINITY_DN154_c0_g1 TRINITY_DN154_c0_g1_i1

見ればわかるが、ヘッダーを少し変えただけ。このスクリプトに頼らずとも、配列のheaderを少し改良するだけで作成できる（*1）。

２、Trinotateを実行する。trinity.fastaと１で作成したgene_trans_map、configファイルを最低限指定する必要がある。

#trinityのアセンブリ結果のtrinity.fastaとgene_trans_mapを指定する
autoTrinotate.pl --Trinotate_sqlite Trinotate.sqlite --transcripts Trinity.fasta --gene_to_trans_map gene_trans_map --conf config.txt --CPU 12

ラン中にどこかのステップでエラーが起きたら、こちら

http://arcsb2017.weebly.com/uploads/4/6/1/9/46199061/trinotate_install.pdf

を確認してエラーが起きたステップのコマンドを実行し、エラーが起きるか調べてください。自分はsignalPのステップでエラーが起きた時、上記リンク先のコマンドを打ってオプションが認識されないことでsignalPのバージョンが違うことに気づきました。

TrinotateレポジトリのTrinotate/sample_data/data/Trinity.fastaはexampleデータですが、これを使ってうまくランするかテストすると、私の環境ではrnammerのランでエラーになるため使えませんでした（rnammerは正常に動いており、実際、rnammerだけ動かすと他のデータではrRNAを予測できても、このファイルではrRNAを予測できない）。注意して下さい。

追記

ラン最後のsqliteデータベースを作るところでTrinotateコマンドとextract_GO_assignments_from_Trinotate_xls.pl, annotation_importerがないとエラーが出てしまった。Trinotateのパスを調べてエラーの通りのパスにシンボリックリンクを置くことで対処した。

ln -s /root/miniconda3/envs/trinotate/bin/Trinotate /root/miniconda3/envs/trinotate/bin/..//

mkdir /root/miniconda3/envs/trinotate/bin/..//util
ln -s /root/miniconda3/envs/trinotate/bin/extract_GO_assignments_from_Trinotate_xls.pl /root/miniconda3/envs/trinotate/bin/..//util/

mkdir /root/miniconda3/envs/trinotate/bin/..//util/annotation_importer
ln -s /root/miniconda3/envs/trinotate/bin/import_transcript_names.pl /root/miniconda3/envs/trinotate/bin/..//util/annotation_importer/

Trinotateはチェックポイントを作っています。どこかのステップでランに失敗した場合、そのステップのコマンドが動作するようにして、失敗した同じコマンドを打てば、そのステップから再開でき、無駄な時間を取られずに済みます。

追記

Error, cannot find obo/go-basic.obo.gz

が出たら、ここを参照して修復する。

cd /root/miniconda3/pkgs/perl-extutils-makemaker-7.36-pl526_1/lib/5.26.2/x86_64-linux-thread-multi/
mkdir obo
cd obo/
cp <your_path>/Trinotate/PerlLib/obo/go-basic.obo.gz .
gzip -dv go-basic.obo.gz

追記

入力配列の数が多すぎると最後のデータベース作成でエラーになることがあるようです。50000配列（61MB）ではエラーを起こしましたが、半分ずつに分けてランすると最後までエラーを吐かずに動かすことができるようになりました。 

引用

The De Novo Transcriptome and Its Functional Annotation in the Seed Beetle Callosobruchus maculatus

Ahmed Sayadi, Elina Immonen, Helen Bayram, and Göran Arnqvist*

PLoS One. 2016; 11(7): e0158565. Published online 2016 Jul 21.

関連 

コメント

dockerイメージもdockerhubで探せば見つかります。導入に苦戦している方は探してみてください。

*1

https://github.com/Trinotate/Trinotate.github.io/issues/7

例えば

>TRINITY_DN144_c0_g1_i1 len=231 path=[227:0-162 228:163-230] [-1, 227, 228, -2]

というヘッダー名なら、

TRINITY_DN144_c0_g1 TRINITY_DN144_c0_g1_i1

にする。もう少しわかりやすく言うと、ヘッダーの赤字の部分

TRINITY_DN144_c0_g1_i1

をtab区切りで二回書く。２回目は"_i1"をつける。すなわち

TRINITY_DN144_c0_g1 TRINITY_DN144_c0_g1_i1

にする。これを全配列に対して実行するスクリプトを書けばよい。

シミュレーション精度と速度が改善された DeepSimulator1.5

2020 2/1 タイトル追加、文章追加、誤字修正

2020 2/2 誤字修正

2020 3/9 コマンド修正

　ナノポアシーケンスは、主要な第3世代シーケンステクノロジーの1つである。 Nanoporeデータの処理と分析を容易にするために、多くの計算ツールが開発された。以前、DeepSimulator1.0（DS1.0）を開発した。DeepSimulatorは、生の電気信号とリードの両方を生成するNanoporeシーケンスの最初のシミュレータである。ただし、DS1.0は高品質のリードを生成できるが、一部のシーケンスでは、シミュレートされた生の信号と実際の信号との相違が大きくなる場合がある。さらに、Nanoporeシーケンシング技術は、DS1.0がリリースされてから大きく進化した。したがって、これらの変更に対応するためにDS1.0を更新する必要がある。
ここではDeepSimulator1.5（DS1.5）を提案する。これらの3つのモジュールはすべて、DS1.0から大幅に更新されている。シーケンスジェネレーターについては、最新の実際のリードの機能を反映するようにサンプルリード長の分布を更新した。 DeepSimulatorのコアであるシグナルジェネレーターの観点から、もう1つのポアモデル、コンテキスト非依存ポアモデルを追加した。これは、以前のコンテキスト依存モデルよりもはるかに高速である。さらに、生成された信号を実際の信号により類似させるために、ポアモデル信号を後処理するローパスフィルターを追加した。 basecallerに関しては、GPUとCPUの両方をサポートできる最新の公式basecaller、Guppyのサポートを追加した。さらに、マルチプロセッシング制御、メモリ、およびストレージ管理に関連する複数の最適化が実装されており、DS1.5はDS1.0よりもはるかに使いやすく、軽量なシミュレータになっている。メインプログラムとデータはhttps://github.com/lykaust15/DeepSimulatorから入手できる。

wiki

https://github.com/lykaust15/DeepSimulator/wiki/Parameters-of-DS1.5

インストール

ubuntu16.04のMiniconda2.4.0.5環境でテストした（docker使用、ホストos macos10.14）。Anaconda2が使える古いubuntuマシン(12.04)でも試した。

依存

 Anaconda2 (https://www.anaconda.com/distribution/) or Minoconda2 (https://conda.io/miniconda.html). 

本体　Github

git clone https://github.com/lykaust15/DeepSimulator.git
cd ./DeepSimulator/
./install.sh

> ./deep_simulator.sh

$ ./deep_simulator.sh 

DeepSimulator v1.5 [Sep-26-2019] 

A Deep Learning based Nanopore simulator which can simulate the process of Nanopore sequencing.

USAGE: ./deep_simulator.sh <-i input_genome> [-o out_root] [-D multi_fasta] [-c CPU_num] [-S random_seed] [-B basecaller] 

[-n read_num] [-K coverage] [-l read_len_mean] [-C cirular_genome] [-m sample_mode] 

[-M simulator] [-e event_std] [-u tune_sampling] [-O out_align] [-G sig_out]

[-f filter_freq] [-s signal_std] [-P perfect] [-H home] 

Options:

***** required arguments *****

-i input_genome : input genome in FASTA format.

***** optional arguments (main) *****

-o out_root : Default output would the current directory. [default = './${input_name}_DeepSimu']

-c CPU_num : Number of processors. [default = 8]

-S random_seed : Random seed for controling the simulation process. [default = 0]

0 for a random seed. Use other number for a fixed seed for reproducibility.

-B basecaller : Choose from the following basecaller for the basecalling process. [default = 1] 

1: guppy_gpu, 2: guppy_cpu, 3: albacore.

***** optional arguments (read-level) *****

-n read_num : The number of reads need to be simulated. [default = 100] 

Set -1 to simulate the whole input sequence without cut (not suitable for genome-level).

-D multi_fasta : Whether the input fasta contains multi discontinuous sequences. [default = 1, separate different sequences] 

Set 0 to concatenate different sequences.

-K coverage : This parameter is converted to number of read in the program. [default = 0] 

If both K and n are given, we use the larger one.

-l read_len_mean : This parameter is used to control the read length mean. [default=8000]

-C cirular_genome : 0 for linear genome and 1 for circular genome. [default = 0]

-m sample_mode : Choose from the following distribution for the read length. [default = 3] 

1: beta_distribution, 2: alpha_distribution, 3: mixed_gamma_dis.

***** optional arguments (event-signal) *****

-M simulator : Choose context-dependent(0) or context-independent(1) simulator to generate event. [default = 1]

-e event_std : Set the standard deviation (std) of the random noise of the event. [default = 1.0]

-u tune_sampling : Tuning sampling rate to around eight for each event. [default = 1 to tune] 

Here eight is determined by 4000/450, where 4KHz is the signal sampling frequency, 

and 450 is the bases per second to pass the nanopore.

-O out_align : Output ground-truth warping path between simulated signal and event. [default = 0 NOT to output]

-G out_signal : Output simulated signal in txt format. [default = 0 NOT to output]

***** optional arguments (signal-signal) *****

-f filter_freq : Set the frequency for the low-pass filter. [default = 950] 

[hint]: a higher frequency value would result in better base-calling accuracy.

-s signal_std : Set the standard deviation (std) of the random noise of the signal. [default = 1.0] 

[hint]: tune event_std, filter_freq and signal_std to simulate different sequencing qualities.

-P perfect : 0 for normal mode (with length repeat and random noise). [default = 0]

1 for perfect pore model (without 'event length repeat' and 'signal random noise'). 

2 for generating almost perfect reads without any randomness in signals (equal to -e 0 -f 0 -s 0).

***** home directory *****

-H home : Home directory of DeepSimulator. [default = 'current directory']

Dockerhubにdcokerイメージも出ています。

#HP (link)
docker pull shkao/deepsimulator:1.5
docker run --rm -itv $PWD:/data shkao/deepsimulator:1.5
> /opt/DeepSimulator/deep_simulator.sh

実行方法

テストゲノムを鋳型に10000リード発生させる。CPU版のguppyでbasecallする。

deep_simulator.sh -i example/artificial_human_chr22.fasta -n 10000 -B 2

• -B basecaller : Choose from the following basecaller for the basecalling process. [default = 1] 1: guppy_gpu, 2: guppy_cpu, 3: albacore.
• -n read_num : The number of reads need to be simulated. [default = 100] 
  Set -1 to simulate the whole input sequence without cut (not suitable for genome-level).

テスト（presetモデル）

v1.0でシミュレート

v1.5

縦軸の単位が揃っていないので注意してください（*1）。

追記

実行する時はディスクの空き容量に注意してください。生のfast5 を出すため 、ゲノム全体を一定数カバーするリードを発生させるとかなりの容量になります。10万リード発生させた時は、fast5ディレクトリだけで80GBほどのファイルサイズになりました。

引用

DeepSimulator1.5: a more powerful, quicker and lighter simulator for Nanopore sequencing
Yu Li, Sheng Wang, Chongwei Bi, Zhaowen Qiu, Mo Li,  Xin Gao

Bioinformatics,  Published: 08 January 2020

v1.0

*1

v1.0と同じ条件で同じリード数だけシミューレートすると、（ラン後のfast5のディレクトリサイズは1.5倍以上になっていたのも関わらず）v1.5はv1.0の1/10くらいの時間で終わった。

ONT cDNA ロングリードのエラー修正を行うisONcorrect

2021 1/5  論文引用

　ロングリードを使用したトランスクリプトームシーケンスは、細胞の転写ランドスケープを理解するための強力な方法であることが証明されている（Wyman et al、n.d .; Bayega et al、2018; Byrne、Cole、et al、2019）。ロングリードテクノロジーにより、ほとんどの転写産物をエンドツーエンドでシーケンスできるため、ショートリードでは必要になる複雑なトランスクリプトームアセンブリ手順を克服できる（Gordon et al、2015; Liu et al、2017）。特に、オックスフォードナノポア（ONT）プラットフォームは、その移植性、低コスト、高スループットにより、ロングリードトランスクリプトームシーケンスの最先端のテクノロジーとなっている（Sessegolo et al、2019; Jenjaroenpun et al、2018）。これにより、選択的スプライシングパターン（Byrne et al、2017）、対立遺伝子特異的発現（Byrne et al、2017）、またはタイピング（Cole et al、2019）、RNA修飾（Leger et al、2019; Sessegolo et al、2019; Wongsurawat et al、nd）、新しいアイソフォームの発見（Workman et al、2019; Clark et al、2019; Sessegolo et al、2019）、およびメタトランスクリプトームのサンプル（Semmouri et al.2019）などの研究が可能になる。
　ただし、これまでのONTトランスクリプトーム研究の範囲は、エラー率が比較的高いために制限されていた。Direct RNAシーケンシングおよびcDNAシーケンシングの両方で約14％である（Workman et al、2019）。この制限を克服するための最も一般的なアプローチは、リードをリファレンストランスクリプトーム（たとえばヒトではGENCODE）に対してアラインさせることである（Wyman et al、n.d .; Workman et al、2019）。 これにより、高品質のリファレンスが利用できない場合、限られた技術となり、多くの非モデル生物が排除されてしまう。さらに、リファレンスが利用できる場合でも、通常は個人、細胞、または環境間の配列の違いをキャプチャしないため、欠落した遺伝子座または非常に可変性の高い遺伝子座からのリードの不整合が生じる。これは、リファレンスではを捕捉しない個人間で高い配列多様性を持つ複雑な遺伝子ファミリーで特に問題があることが示されている（Sahlin et al、2018）。エラー率を減らすためのいくつかの実験的な（ウェットの）アプローチが存在する（Lebrigand et al、2019; Cole et al、2019; Volden et al、2018）。しかし、これらは通常、スループットの低下と実験的なオーバーヘッドという犠牲を伴う。
　一方、計算によるエラー修正は、スループットや実験プロトコルをカスタマイズする必要性に影響を与えることなくエラー率を減らす非常に有望なアプローチである。ゲノムリードのエラーを修正するように設計されたツールがある（（Koren et al、nd）、（Tischler and Myers、nd）、（Salmela et al、2016）、（Xiao et al、2017） 、（Chin et al、2013））。しかし、トランスクリプトームのエラー修正は挑戦的であり、同じ遺伝子または遺伝子ファミリー遺伝子座からのリード内の構造的変動性、および、例えば、選択的スプライシング、可変転写によるリード内の非常に可変で領域特異的なカバレッジ、およびさまざまな転写産物の量のため、ゲノムの場合とは異なる。実際、最近の研究では、ゲノム用に設計されたエラー修正プログラムをONTトランスクリプトームデータに適用すると、過剰な修正によるexonの脱落または追加によるアイソフォームランドスケープの変更、またはカバレッジの低いサイトでのリードの分割など、望ましくないダウンストリーム効果があることがわかった（LIma et al、2019）。 ONTトランスクリプトームデータのエラー修正の可能性を実現するには、カスタムアルゴリズムを設計する必要がある。最近の論文はクラスタリング（Sahlin and Medvedev 2019; Marchet et al.2019）とこのデータの向きの問題（Ruiz-Reche et al。2019）に取り組んでいるが、現在、ONTトランスクリプトームリードのエラー修正に利用できるツールはない。
　このホワイトペーパーでは、エラー率を約1％に減らすエラー修正トランスクリプトームcDNA ONTデータの方法を提示する。これにより、リファレンスフリーのトランスクリプトーム解析のための費用対効果の高い完全なcDNA ONTトランスクリプトームシーケンスのためのメソッドを適用する可能性を示す。これらのエラー率は、isONcorrectと呼ばれる新しい計算エラー修正方法によって達成できる。これは、異なるアイソフォームからのリード間で共有されるシーケンス領域を活用する。 IsONcorrectはhttps://github.com/ksahlin/isONcorrectからダウンロードできる。ここでは、ショウジョウバエcDNAについて、変更されたストランドPCS109プロトコル、PCS109スパイクイン（SIRV）データ、およびin silicoデータを使用して生成されたデータを使用してメソッドを評価する。本方法は、ONTトランスクリプトームシーケンスのはるかに広範なアプリケーションへの扉を開く。

インストール

macos10.14のanaconda3.7環境にて、オーサーが推奨する以下の手順で仮想環境を作成してテストした。

依存

• python version >=3.
• spoa (1.1.5)
• edlib (1.1.2)
• NumPy (1.16.2)

本体　Github 

conda create -n isoncorrect python=3 pip -y
conda activate isoncorrect
pip install isONcorrect
conda install -c bioconda spoa -y

> isONcorrect --help

$ isONcorrect --help

usage: isONcorrect [-h] [--version] [--fastq FASTQ] [--k K] [--w W] [--xmin XMIN] [--xmax XMAX] [--T T] [--exact] [--disable_numpy] [--max_seqs_to_spoa MAX_SEQS_TO_SPOA] [--max_seqs MAX_SEQS]

                   [--exact_instance_limit EXACT_INSTANCE_LIMIT] [--set_w_dynamically] [--verbose] [--compression] [--outfolder OUTFOLDER]

De novo error correction of long-read transcriptome reads

optional arguments:

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

  --version             show program's version number and exit

  --fastq FASTQ         Path to input fastq file with reads (default: False)

  --k K                 Kmer size (default: 9)

  --w W                 Window size (default: 10)

  --xmin XMIN           Upper interval length (default: 14)

  --xmax XMAX           Lower interval length (default: 80)

  --T T                 Minimum fraction keeping substitution (default: 0.1)

  --exact               Get exact solution for WIS for evary read (recalculating weights for each read (much slower but slightly more accuracy, not to be used for clusters with over ~500 reads) (default: False)

  --disable_numpy       Do not require numpy to be installed, but this version is about 1.5x slower than with numpy. (default: False)

  --max_seqs_to_spoa MAX_SEQS_TO_SPOA

                        Maximum number of seqs to spoa (default: 200)

  --max_seqs MAX_SEQS   Maximum number of seqs to correct at a time (in case of large clusters). (default: 1000)

  --exact_instance_limit EXACT_INSTANCE_LIMIT

                        Activates slower exact mode for instance smaller than this limit (default: 0)

  --set_w_dynamically   Set w = k + max(2*k, floor(cluster_size/1000)). (default: False)

  --verbose             Print various developer stats. (default: False)

  --compression         Use homopolymenr compressed reads. (Deprecated, because we will have fewer minmimizer combinations to span regions in homopolymenr dense regions. Solution could be to adjust upper interval legnth

                        dynamically to guarantee a certain number of spanning intervals. (default: False)

  --outfolder OUTFOLDER

                        A fasta file with transcripts that are shared between samples and have perfect illumina support. (default: None)

テストラン

容量が大きいので注意。

git clone https://github.com/ksahlin/isONcorrect.git
cd isONcorrect/
isONcorrect --fastq test_data/isoncorrect/0.fastq \
 --outfolder output_dir

Githubには、実際のfull-length ONT cDNA sequencesの選抜からエラーコレクションまでの流れが説明されています。確認して下さい。

"Processing and error correction of full-length ONT cDNA reads is acheved by the pipeline of running pychopper --> isONclust --> isONcorrect"

引用

Error correction enables use of Oxford Nanopore technology for reference-free transcriptome analysis
Kristoffer Sahlin​, Botond Sipos, Phillip L James,​ Daniel J Turner,​ Paul Medvedev

bioRxiv preprint first posted online Jan. 8, 2020

2021 1/5

Error correction enables use of Oxford Nanopore technology for reference-free transcriptome analysis

Kristoffer Sahlin, Botond Sipos, Phillip L. James & Paul Medvedev
Nature Communications volume 12, Article number: 2 (2021)

関連

UniprotのID変換webサービスを使い、UniProt accessionsからタンパク質のアノテーションを得る (ID mapping)

2020 2/4 追記

UniProtのRetrieve/ID mappingサービスを使用すると、UniProt accessions IDからGenbankの配列、PDBのID、Entrez Gene ID、GI nnumber、タンパク質のアノテーションなどに変換できる。

Converting UniProt identifiers to external identifers (or vice versa)

https://www.ebi.ac.uk/training/online/course/uniprot-exploring-protein-sequence-and-functional/how-use-uniprot-tools/batch-retrieval-id--0

使い方

Retrieve/ID mapping

https://www.uniprot.org/uploadlists/ にアクセスする。

Retrieve/ID mappingの使い道はたくさんあるが、ここでは、UniProtのaccessions IDからタンパク質のアノテーションを得るために利用する。Uniprotをデータベースとして、Diamondによるblast検索を実行して得たtabularファイルになる。

IDをペーストする。

入力フォーマットを指定する。

次に出力フォーマットを指定する。ここではUniprotKBを指定（またはGene name）。

出力された。

結果は様々な形式でダウンロードできる。

excel形式でダウンロードしてexcelで開いた。

追記

Uniprot propteomeへのDIAMOND blastx検索結果をクエリにする

UniProtKB => Uniref100

出力

引用

Retrieve/ID mapping

関連

高速なヒトゲノムのアセンブラ Peregrine

　初期のヒトゲノムプロジェクトと安価なDNAシークエンシング技術の技術の開発は、学術研究とゲノム情報を使用して人間の健康を改善する産業の両方を進歩させた。それは、遺伝子型と表現型の関連と多くの重要かつ臨床関連のアプリケーションのための貴重な情報につながる。 Oxford Nanopore Technology（ONT）およびPacific Biosciences（PacBio）のシーケンサーは、第2世代のシーケンスよりも桁違いに長いDNA配列を読み取ることができる。リード長が長くなると、de novoアセンブリが比較的簡単になり、より連続したアセンブリを生成できるようになる。ゲノムのdenovo再構築により、リファレンスを事前情報として使用することへの依存が軽減される。リファレンスに依存するリシークエンシングアプローチは、リファレンスから大幅に逸脱したゲノム構造を探索するのに効果的ではない場合がある。直接的なヒトゲノムアセンブリの最近の研究により、現在のリファレンスゲノムにない新しい配列が発見された。構造変異を発見する体系的なアプローチは、多くの新しい構造変異と、ヒトのSNPおよびsmall indels以外のより大きな変異のより包括的なカタログを生成するために、より多くのサンプルを必要とする。

（一段落省略）
　ロングリードゲノムアセンブリの場合、現在のほとんどのロングリードアセンブラのにオーバーラップレイアウトコンセンサス（OLC）パラダイムが使用される。リードのquadraticな比較（おそらくリード数に応じて二次的に増加するという意味合い）が計算効率をさらに向上させるための主要なボトルネックのままである。たとえば、最初にノイジーなPacBioのリード用に設計されたhierarchical genome assembly process HGAPは、2つの重複するステップを取る、1つはエラー修正用で、もう1つはアセンブリグラフ生成用で、ノイズの多い配列からヒトゲノムをアセンブリするのに20,000から30,000 CPU時間必要である。最近開発されたほとんどの高性能アセンブラ、たとえばFlye、wtdbg、Shastaは、2つの完全なリードの間でこのような計算的に高価な明示的なオーバーラッピングステップを避けるため、新しい戦略を適応させる。このような最適化は、ノイズの多いリードを効率的にアセンブリするために必要になる可能性がある。一方、ゲノムアセンブリのオーバーラップしたリードの計算効率を向上させるために、より正確なコンセンサスリードから始めることができる。より良いリード精度を活用することで、オーバーラップ検出の計算の複雑さを軽減できることがわかった。新しいゲノムアセンブラPeregrine（https://github.com/cschin/peregrine）を開発した。コンセンサスリードから始まり、ヒトゲノムサイズのアセンブリを行うために、効率的なインデックススキームを利用して、計算時間を20 cpu時間まで短縮した。
ここで紹介するゲノムアセンブリアプローチは、クラスター計算用の特別なセットアップなしで、ゲノムアセンブリを日常作業にするのに効果的である。ゲノムプロジェクトのこのような基本的なプロセスを簡素化することは、リシークエンシングアプローチからの情報の欠落を回避するde novoゲノムアセンブリを定期的に生成するための鍵である。現在、de novoアプローチのコストは、リシークンシングよりも高価である。しかしながら、速いペースで進む計算方法とシークエンシング技術の向上により、de novoでヒトゲノムを生成するためのコストは、間も無く到来するオーダーメイド医療で受け入れられる価格まで下がるかもしれない。我々（本著者ら）の方法はまた、現在のヒトのリファレンスGRCh38では利用できない新規のヒトゲノム配列をキャプチャすることを可能にし、pan-human-genomics references の構築に役立つ。ヒトゲノムのより包括的な見解を提供する可能性を考え、我々（本著者ら）は、今回の全ゲノムアプローチがリシークンシングアプローチでは簡単に明らかにできない遺伝子疾患の重要な情報を提供することを望んでいる。（複数段落省略）

　単一のノードで、クラウドアクセス可能なハードウェアを使用して、約2時間でヒトゲノムアセンブリを実行できることを実証した。 これにより、ソフトウェアインフラストラクチャ要件とグリッドコンピューティングの専門スキルの両方が不要になる。
この普遍的なアクセスは、DNAシーケンス技術の迅速な進歩と相まって、リファレンスグレードのゲノムアセンブリの日常的な生成を可能にする。サンプル収集から結果まで1〜2日という速さになる。

インストール

オーサーらが準備したdockerイメージをpullしてテストした。

本体　Github

dockerhub (link)
docker pull cschin/peregrine:latest

Usage:

  pg_run.py asm <reads.lst> <index_nchunk> <index_nproc>

                            <ovlp_nchunk> <ovlp_nproc>

                            <mapping_nchunk> <mapping_nproc>

                            <cns_nchunk> <cns_nproc>

                            <sort_nproc>

                            [--with-consensus]

                            [--with-L0-index]

                            [--output <output>]

                            [--shimmer-k <shimmer_k>]

                            [--shimmer-w <shimmer_w>]

                            [--shimmer-r <shimmer_r>]

                            [--shimmer-l <shimmer_l>]

                            [--best_n_ovlp <n_ovlp>]

                            [--mc_lower <mc_lower>]

                            [--mc_upper <mc_upper>]

                            [--aln_bw <aln_bw>]

                            [--ovlp_upper <ovlp_upper>]

  pg_run.py (-h | --help)

  pg_run.py --verison

Options:

  -h --help                   Show this help

  --version                   Show version

  --with-consensus            Generate consensus after getting the draft contigs

  --with-L0-index             Keep level-0 index

  --output <output>           Set output directory (will be created if not exist) [default: ./wd]

  --shimmer-k <shimmer_k>     Level 0 k-mer size [default: 16]

  --shimmer-w <shimmer_w>     Level 0 window size [default: 80]

  --shimmer-r <shimmer_r>     Reduction factore for high level SHIMMER [default: 6]

  --shimmer-l <shimmer_l>     number of level of shimmer used, the value should be 1 or 2 [default: 2]

  --best_n_ovlp <n_ovlp>      Find best n_ovlp overlap [default: 4]

  --mc_lower <mc_lower>       Does not cosider SHIMMER with count less than mc_low [default: 2]

  --mc_upper <mc_upper>       Does not cosider SHIMMER with count greater than mc_upper [default: 240]

  --aln_bw <aln_bw>           Max off-diagonal gap allow during overlap confirmation [default: 100]

  --ovlp_upper <ovlp_upper>   Ignore cluster with overlap count greater ovlp_upper [default: 120]

テストラン

docker run -it --rm -v $PWD:/wd cschin/peregrine:latest test

 ラン途中で確認が出る、2019/10現在はまだプレリリースで、再配布は許可されていない。agreeならyes。

またランした時に追記します。

引用 

Human Genome Assembly in 100 Minutes

Chen-Shan Chin, Asif Khalak

bioRxiv preprint first posted online Jul. 17, 2019

WikiPathwaysのCytoscapeプラグイン WikiPathways App for Cytoscape

　このホワイトペーパーでは、Cytoscape用のオープンソースWikiPathwaysアプリ（http://apps.cytoscape.org/apps/wikipathways）を紹介する。WikiPathwaysアプリは、データの視覚化とネットワーク分析のためにバイオロジカルパスウェイをインポートするために使用できる。 WikiPathwaysは、手動でダウンロードしたりWebサービスを使用したりするための完全にアノテーションが付けられたパスウェイマップを提供する、オープンでcollaborativeなバイオロジカルパスウェイデータベースである。 WikiPathwaysアプリを使用すると、ユーザーは2つの異なるビューでパスウェイを読み込むことができる。

　2013年8月の最初のリリースから2014年5月の論文のサブミットまでの2000以上のダウンロードは、ネットワーク生物学分野でのアプリの重要性を強調している。

Tutorial

https://cytoscape.org/cytoscape-tutorials/protocols/wikipathways-app/#/title

Visualizing Expression Data

https://cytoscape.org/cytoscape-tutorials/protocols/mapping-data/index.html#/

Cytoscape User manual

WikiPathways

インストール

macos10.14にて、Cytoscape v3.7.2を使ってテストした。

Cytoscape App Store

１、リンク先の右端にあるdowloadから、ビルド済みの.jarファイルをダウンロードする。

２、cytoscape3.7（link）を立ち上げ、.jarファイルを読み込む。Apps => App Manager...

３、ウィンドウ左下のinstall from fileを選択。

 ４、プラグインとしてcytoscapeにインストールされた。ウィンドウを閉じる。

使い方

ここからはcytoscapeで作業する。Networkの左下にある▼マークをクリック、一覧からwikipathwaysを選択。

  ウィンドウ内にタイプして、pathwayを検索する。ここではstatin（スタチン）とタイプ。　ウィンドウが出現するので、読み込むパスウェイを選択する。

ここではヒトのStatin pathwayを選択した。左下のimport as Pathwayボタンをクリックしてロードする。

補足==============================================================

ウィンドウ右下のpreviewボタンを押すとパスウェイを確認できる。

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

読み込まれた。ロード時、ノードには対応するEnsemblの識別子がつけられる（下の表を見ると、右端のカラムにEnsembl identifiersがついている）。

Ensemblの識別子が付いているため、対応するEnsembl識別子が付いたCSVを読み込めば、その識別子に対応するnodeやedgeに、数値の大小に応じた色をアサインできる（nodeやedgeの色、太さなどを変えるのはcytoscape標準の機能）。

ここでは発現量情報のテーブルを読み込んで、nodeやedgeに発現量に応じた色を付ける。こちらのデモデータを使用する。1カラム目はEnsembl識別子、２カラム目はcase/controlのlog2 fold changeになっている。

読み込んでいく。File => import => Table from Fileを選択。

テーブル読み込みウインドウが出てきたらTo a Network Collectionを選択。

Key Column fot NetworkはEnsemblを選択 (判別に使う情報の列がどれかを教えてあげる)。

実際にノードに色をつけていく。Style タブ=>Fill colorを選択。

Fil Color左の＞をクリックすれば上のように展開される。

columnはhyper〜を選択。

Mapping typeはContinuous mappingを選択。

fold changeに応じて色がついた。

細かな色指定もできます。グラデーションのパレットをクリック。

最大と最小の値とその時の色、境界、など全てマニュアル指定できる。

パスウェイデータを読み込むときにimport as networkを選ぶと、エッジとノードからなる純粋なグラフ(ネットワーク)として読み込まれる。

ネットワークとして読み込み、ノードにfold changeによる色をつけた（手順は上の通り）。

引用

WikiPathways App for Cytoscape: Making biological pathways amenable to network analysis and visualization

Kutmon M, Lotia S, Evelo CT, Pico AR

Version 2. F1000Res. 2014 Jul 1 [revised 2014 Sep 11];3:152

参考

CYTOSCAPEを使ったデータの可視化

PDF版

http://motdb.dbcls.jp/?plugin=attach&pcmd=open&file=AJACS62_Cytoscape.pdf&refer=AJACS62

インタラクティブなdot plotビューア dot

ゲノムとゲノムのアライメントを視覚化する方法の１つに、両方のゲノムのアライメントの概要を提供するドットプロットがある。一般的なゲノムブラウザ（IGVやUCSCゲノムブラウザなど）が1つの次元でゲノムデータをプロットするのとは異なり、ドットプロットでは、一方の軸にリファレンスゲノムが配置され、もう一方の軸にクエリのゲノムが配置される。 2つのゲノム間のアライメントは、両方のゲノムの座標に従って配置される。

DNAnexusインタラクティブなドットプロットビューアDotは、ラージゲノムをすばやくプロットできる。ユーザーは関心ある領域にズームインして、２ゲノム間のアラインメントを調べることができる。

Dot: An Interactive Dot Plot Viewer for Comparative Genomics

https://blog.dnanexus.com/2018-01-11-dot-an-interactive-dot-plot-viewer-for-comparative-genomics/

Meet Dot: An Interactive Dot Plot Viewer for Comparative #Genomics https://t.co/sZJIJH9Snl Learn more at #PAG18

— DNAnexus, Inc. (@dnanexus) 2018年1月12日

Github

ヘルパースクリプト

> python DotPrep.py -h

$ python DotPrep.py -h

usage: DotPrep.py [-h] --delta DELTA [--out OUT]

                  [--unique-length UNIQUE_LENGTH] [--overview OVERVIEW]

Take a delta file, apply Assemblytics unique anchor filtering, and prepare

coordinates input files for Dot

optional arguments:

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

  --delta DELTA         delta file

  --out OUT             output file

  --unique-length UNIQUE_LENGTH

                        The total length of unique sequence an alignment must

                        have on the query side to be retained. Default: 10000

  --overview OVERVIEW   The number of alignments to include in the coords.idx

                        output file, which will be shown in the overview for

                        Dot. Default: 1000

on line

1, nucmer （オプションなしで実行する）

nucmer --prefix=nucmer_output ref_genome.fa query_genome.fa

2, helper scriptのラン（必要に応じてパラメータ設定する）

python DotPrep.py --delta nucmer_output.delta --out OUTPUT

出力

https://dnanexus.github.io/dot/にアクセスする。.coordsとindexファイルを読み込む。

視覚化された。直感的に操作できるので説明は省略する。

引用

GitHub - MariaNattestad/dot: Dot: An interactive dot plot viewer for comparative genomics

関連