macでインフォマティクス

macでインフォマティクス

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

エラーの多いロングリードのエラー修正とアセンブリを行う NECAT

2020 2/7  パラメータエラー修正、2/8 わかりにくい表現を修正、3/12 わかりにくい説明を修正

2021 1/5 論文引用、9/8 関連研究のリンク追記

2023/04/12 コマンド微修正(.pl削除)

 

 

 ナノポアのロングリードはde novoゲノムアセンブリで有利だが、ゲノム研究への適用は、これらロングリードの複雑なエラーによって依然として妨げられている。 ここでは、ナノポアのリードの複雑なエラーを克服するために設計されたエラー修正およびde novoアセンブリツールであるNECATを開発した。 Nanoporeのリードを高精度にすばやく修正するために、 adaptive read selection and two-step progressive methodを提案した。 Nanoporeリードの全長を利用するために、2段階アセンブラを導入した。 NECATは、エラー修正とNanoporeリードのde novoアセンブリの両方で優れたパフォーマンスを実現する。 NECATは、35Xカバレッジのヒトゲノムをアセンブルするのに7,225 CPU時間しか必要とせず、NG50で2.28倍の改善を達成する。 さらに、ヒトWERI細胞株のアセンブリでは、29 MbpのNG50が示された。 Nanoporeリードの高品質なアセンブリにより、構造変化の検出における偽陽性を大幅に減らすことができる。

 

インストール

ubuntu19.04の計算機を使ってテストした。

Sucessfully tested  on

本体 Github

#linux向けバイナリ
wget https://github.com/xiaochuanle/NECAT/releases/download/SourceCodes20200119/necat_20200119_Linux-amd64.tar.gz
tar xzvf necat_20200119_Linux-amd64.tar.gz
cd NECAT/Linux-amd64/bin

#パスを通す
export PATH=$PATH:$(pwd)

#conda 2022/06/02
mamba create -n necat -y
conda activate necat
mamba install -c bioconda necat -y

 > necat

Usage: necat correct|assemble|bridge|config cfg_fname

    correct:     correct rawreads

    assemble:    generate contigs

    bridge:      bridge contigs

    config:      generate default config file

 

 

 

実行方法

1、デフォルトパラメータのコンフィグファイル作成

以下のコマンドを打つとデフォルトパラメータのコンフィグファイルが作成される。

necat config config.txt

> cat config.txt

PROJECT=

ONT_READ_LIST=

GENOME_SIZE=

THREADS=4

MIN_READ_LENGTH=3000

PREP_OUTPUT_COVERAGE=40

OVLP_FAST_OPTIONS=-n 500 -z 20 -b 2000 -e 0.5 -j 0 -u 1 -a 1000

OVLP_SENSITIVE_OPTIONS=-n 500 -z 10 -e 0.5 -j 0 -u 1 -a 1000

CNS_FAST_OPTIONS=-a 2000 -x 4 -y 12 -l 1000 -e 0.5 -p 0.8 -u 0

CNS_SENSITIVE_OPTIONS=-a 2000 -x 4 -y 12 -l 1000 -e 0.5 -p 0.8 -u 0

TRIM_OVLP_OPTIONS=-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 1 -a 400

ASM_OVLP_OPTIONS=-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 0 -a 400

NUM_ITER=2

CNS_OUTPUT_COVERAGE=30

CLEANUP=1

USE_GRID=false

GRID_NODE=0

GRID_OPTIONS=

SMALL_MEMORY=0

FSA_OL_FILTER_OPTIONS=

FSA_ASSEMBLE_OPTIONS=

FSA_CTG_BRIDGE_OPTIONS=

POLISH_CONTIGS=true

 

ランするために以下のように編集した(赤字部分)。"read_list.txt"はリードのパスを示したテキスト。下の説明参照。

> cat config.txt

PROJECT=Arabidopsis

ONT_READ_LIST=read_list.txt

GENOME_SIZE=130000000

THREADS=40

MIN_READ_LENGTH=3000

PREP_OUTPUT_COVERAGE=40

OVLP_FAST_OPTIONS=-n 500 -z 20 -b 2000 -e 0.5 -j 0 -u 1 -a 1000

OVLP_SENSITIVE_OPTIONS=-n 500 -z 10 -e 0.5 -j 0 -u 1 -a 1000

CNS_FAST_OPTIONS=-a 2000 -x 4 -y 12 -l 1000 -e 0.5 -p 0.8 -u 0

CNS_SENSITIVE_OPTIONS=-a 2000 -x 4 -y 12 -l 1000 -e 0.5 -p 0.8 -u 0

TRIM_OVLP_OPTIONS=-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 1 -a 400

ASM_OVLP_OPTIONS=-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 0 -a 400

NUM_ITER=2

CNS_OUTPUT_COVERAGE=30

CLEANUP=1

USE_GRID=false

GRID_NODE=0

GRID_OPTIONS=

SMALL_MEMORY=0

FSA_OL_FILTER_OPTIONS=

FSA_ASSEMBLE_OPTIONS=

FSA_CTG_BRIDGE_OPTIONS=

POLISH_CONTIGS=true

 

もう1つ、リードのパスを示したファイルread_list.txtを準備する。suffixをチェックされるので”.txt"をつける。findで作るなら以下のように

find ./ONT_fastq/nanopore*gz -type f > read_list.txt

相対パスになるが、カレントからランするならこれでO.K(違うならフルパスで)

補足;使えるシーケンスデータは、ロングリードの rawまたはgzipedのfastqとfast。また.fqだとエラーになる。.fastqにした。

 

2、error correction

作成したconfigファイルを指定する。

necat correct config.txt 

指定したproject name、ここではArabidopsis/に結果の1-consensus/cns_final.fasta.gzが出力される。 全てのリードがエラー修正されるわけではない点に注意(Github参照)。

 

3、Assembly

以下のコマンドを打つ。自動で1の出力を認識してアセンブリが開始される。

necat assemble config.txt

4-fsaに中間ファイルやcontigs.fasta、polished_contigs.fastaが出力される。

 

4、Bridging

necat bridge config.txt

6-bridge_contigs/bridged_contigs.fasta が出力される。

 

計算機クラスタでランするためのオプションも用意されています。 Githubで確認して下さい。

 

感想

リファレンスのあるナズナのONTリードを使ってテストしたところ、Flyeと同じくらいの連続性と精度を持った配列が出力されていました。

(データの質や生物種によって傾向は変わる可能性があります)

引用

Fast and accurate assembly of Nanopore reads via progressive error correction and adaptive read selection

Ying Chen, Fan Nie, Shang-Qian Xie, Ying-Feng Zheng, Thomas Bray, Qi Dai, Yao-Xin Wang, Jian-feng Xing, Zhi-Jian Huang, De-Peng Wang, Li-Juan He, Feng Luo, Jian-Xin Wang, Yi-Zhi Liu, Chuan-Le Xiao

bioRxiv preprint first posted online Feb. 2, 2020

 

2021 1/5

Efficient assembly of nanopore reads via highly accurate and intact error correction

Ying Chen, Fan Nie, Shang-Qian Xie, Ying-Feng Zheng, Qi Dai, Thomas Bray, Yao-Xin Wang, Jian-Feng Xing, Zhi-Jian Huang, De-Peng Wang, Li-Juan He, Feng Luo, Jian-Xin Wang, Yi-Zhi Liu & Chuan-Le Xiao
Nature Communications volume 12, Article number: 60 (2021)

 

関連


2021 9/8

バナナゲノムのT2Tアセンブリ(バージョン4アセンブリ)の報告。ONTのリードも使用されていて、いくつかのアセンブラを検討した結果、NECATが一次アセンブリに使用されています(Supplementary Table.8に比較結果)。

 

2022/06/05

PacbioとNanoporeのリードを使った比較論文。NECATとNextDenovoが良いパフォーマンスを示しています(Table.S2)。


 

KEGGの遺伝子アノテーション結果を要約する GAEV

2020 10/20 追記

 

 非モデルアセンブリのコンピュータ アノテーション付き遺伝子の生物学的機能と、これらの遺伝子の産物によって形成される分子パスウェイの説明は、種のさまざまな固有の生物学的属性(生理学、生活史、行動など)の遺伝的基盤を特定するために重要である。 DNA /タンパク質データベースに対する計算検索、例えばNCBI Blast(Boratyn et al、2013)、UniProt(Bateman et al、2017)、InterPro(Finn et al、2017)、Blast(Camacho et al、2009)、InterProScan(Jones et al、2014)、Hmmer(Mistry et al、2013)などのツールを使うと個々の遺伝子機能を予測できる。対照的に、単一の種の遺伝子スイート全体によってエンコードされた分子パスウェイを描くことは、特に非モデル種の場合、はるかに困難な作業になる。集中的に研究されたモデル生物に由来する分子パスウェイへの遺伝子のマッピングは、このニーズに対処するためのエントリーポイントを提供している。

 遺伝子を既知の分子パスウェイにマッピングするために、Kyoto Encyclopedia of Genes and Genomes (KEGG) は包括的なWebサービスを提供している(Kanehisa et al、2017; Kanehisa&Goto、2000; Kanehisa et al、2016a)。KEGGは、ゲノムシーケンスの生物学的解釈のための統合データベースである。遺伝子の分子機能は、オルソロググループ、つまりKEGG Orthology(KO)を使用して分類される。KEGGには、KEGGパスウェイ、BRITE階層、およびKEGGモジュールも含まれている。これらはすべてKOノードのネットワークである。 WebサーバーBlastKOALAおよびGhostKOALAを通じて提供されるKEGG自動アノテーションサービスを使用して、完全/部分ゲノムアセンブリまたはメタゲノムデータセットからの一連の遺伝子の分子機能およびそれらのエンコードされた分子パスウェイにアノテーションを付けることができる(Kanehisa et al、2016b)。モデルではない種については、KAAS(KEGG Automatic Annotation Server)Webサービスを使用して、遺伝子の完全なセットまたはランダムなセットにアノーテーションを付けて、それらの分子機能を記述し、特定された分子パスウェイにマッピングできる。アノテーション結果は、各遺伝子のKO番号、KEGGパスウェイデータベースにマップされた遺伝子、およびBRITEにマップされた遺伝子で構成される。それにもかかわらず、結果として得られるパスウェイとBRITE階層の完全なセットは、KEGGが提供する一時URLを介してのみ表示され、分析が完了してから数日間しか利用できない。これらの結果は、管理されたKEGGパスウェイまたはBRITE階層のいずれかで編成されているが、KAASは、遺伝子機能とパスウェイの統合的な遺伝子中心のビュー、つまり、遺伝子機能と各遺伝子のすべての関連分子パスウェイの完全な要約jは提供しない。

 想像できるように、KEGGオルソロジーKEGGパスウェイに基づく遺伝子の function annotationを統合すると、新しくアセンブリされたゲノムまたはメタゲノムデータセットの予測遺伝子と関連するパスウェイの両方を特徴付ける効率的な方法を提供できる。KEGGデータベースが提供するAPIインターフェースを使用してKEGGパスウェイを取得するための多数の計算パッケージ(例:Moutselos et al、2009; Wrzodek et al、2011)が存在するにもかかわらず、これらのパッケージはどれも、新しくアセンブリされたゲノムに含まれる分子パスウェイの完全なセットを再構築することはできない。KEGGの非常に有益なリソースを利用する手段を提供して非モデル種のゲノム配列と分子パスウェイにアノテーションを付けるため、KEGG APIを使用してKEGGオルソロジーアノテーションKEGGパスウェイマッピングの結果を統合するGene Annotation Easy Viewer(GAEV)を開発した。 GAEVは、遺伝子機能とパスウェイの遺伝子中心の視点、つまり、遺伝子機能の完全な要約と、各遺伝子に関連する可能性のあるすべての分子パスウェイを提供することを目的としている。これは、MEGAN(Huson et al、2016)やMinPath(Ye&Doak、2009)などの他のKEGG関連ソフトウェアとは異なる。 GAEVはPython 3で実装されており、独立したパッケージとして使用できる。

 

インストール

anaconda3.7環境でテストした(macos10.14使用)。

本体 Github

git clone https://github.com/UtaDaphniaLab/Gene_Annotation_Easy_Viewer.git
cd Gene_Annotation_Easy_Viewer/gene_annotation_easy_viewer/

  

実行方法

1、リストの準備

KAAS、BlastKOALA、GhostKOALA、またはKofamKOALA(paper)を使い、配列にKEGG  Identifier(KO)をアサインしたファイルを準備する。

Download KO listからダウンロードできる。

f:id:kazumaxneo:20200206154350p:plain

このようなファイルになる(Githubのexampleデータ)。

f:id:kazumaxneo:20200206152539p:plain

 

2、実行

GAEV を使うには本体を叩く。

python GAEV.py

対話式で進める。バッチモードもあるが1を選択した。

$ python GAEV.py 

 

Would you like to:

   1) Create a new data file and generate a table from a new dataset of KO numbers

   2) Create a new data file and generate a table from a new dataset of KO numbers (Batch)

   3) Generate a new table from an existing data file

   4) Generate a new table from an existing data file (Batch)

 

Input a digit for your choice: 

 

リストファイルを指定する(カレントにないならフルパスで)。

Input a digit for your choice: 1

 

Please enter either the relative or absolute path to the input file below

query.ko.txt

(ファイルのパスをコピペしたなら最後にスペースが入る。残っているとエラーになる。)

 

抽出開始

f:id:kazumaxneo:20200206153257p:plain

 

何らかのキーでフィルタリングするか聞かれる。ここではno

Extracting data from KEGG (complete)

 

Would you like to filter the data?

   1) Yes

   2) No

 

 

 

出力ファイル名を指定

Enter the output file name or press ENTER to use the default name [query.ko]:

 

output

 

出力内容を聞かれる。ここでは遺伝子名のテーブルと代謝マップにアサインされた結果が欲しいので2を選択。

Would you like to:

   1) Generate a table of genes (HTML + txt)

   2) Generate a table of genes and pathways (HTML + txt)

   3) Generate a table of pathways (HTML + txt)

   4) Generate a only a table of genes without links to pathway maps (txt, tab-delimited, small size)

 

   Note: The txt file is tab delimited and easily manipulated in text editors, but does not contain embedded links.

 

 

3、出力

IDに遺伝子名と注釈が付いた。

f:id:kazumaxneo:20200206154815p:plain

 

GAEVの出力内容選択時に2か3を選択していればhtmlレポートも出力される。隠しているが、正しく遺伝子名がアサインされている。

f:id:kazumaxneo:20200206153939p:plain

f:id:kazumaxneo:20200206154227p:plain

クリックするとジャンプする。

f:id:kazumaxneo:20200206154107p:plain

 

追記

KAASの結果と比較すると違いが見つかることがありました。便利なツールですが注意して使って下さい。

引用

Gene Annotation Easy Viewer (GAEV): Integrating KEGG's Gene Function Annotations and Associated Molecular Pathways

Huynh T, Xu S

Version 3. F1000Res. 2018 Mar 29 [revised 2019 May 9];7:416

 

KEGGのパスウェイアノテーションwebサービス KAAS

 2020 2/6  タイトル修正

 

 近年、完全(complete)なゲノムとドラフトゲノムの数は急速に増加しており、これらのゲノムの遺伝子の機能的特性と生物学的役割の特定を自動化することがますます重要になっている。 KEGGデータベースでは、Smith–Watermanスコアを使用したベストヒット情報と手動キュレーションに基づいて、完全なゲノムの遺伝子にKEGG orthology (KO) identifiersまたはK番号の注釈が付けられている。各K番号は遺伝子のオーソロググループを表し、KEGGパスウェイマップまたはBRITE機能階層内(link)のオブジェクトに直接リンクされている。ここでは、KAAS(KEGG Automatic Annotation Server:http://www.genome.jp/kegg/kaas/)と呼ばれるWebベースのサーバーを開発した。つまり、K番号をゲノムの遺伝子に自動的に割り当てる迅速な方法の実装で ある。これにより 、KEGGパスウェイとBRITE階層の再構築が可能になる。この方法は、配列の類似性、双方向のベストヒット情報、およびいくつかのヒューリスティックに基づいており、手動でキュレーションされたKEGG GENESデータベースと比較した場合に高い精度を達成している。

 

KAAS help

https://www.genome.jp/kegg/kaas/help.html

 

webサービス

https://www.genome.jp/kegg/kaas/にアクセスする。

f:id:kazumaxneo:20200206103932p:plain

 complete or draft genomeであれば1番上を選択。

f:id:kazumaxneo:20200206113812p:plain

アミノ酸配列をペーストするか、アミノ酸配列のmulti fastaをローカルからアップロードする。cDNA/EST などに由来する塩基配列も指定できる(blastのみ)。

f:id:kazumaxneo:20200206114024p:plain

 

GHOSTXよりblastの方が精度は高いが(論文の図6参照)、配列数が多いならGHOSTXに切り替えるとジョブが早く終わる(SBH (single-directional best hit)にするとさらに半分の時間になる)。

 

クエリ名、メールアドレス、KOアサインのための遺伝子セット(link)を指定する。

f:id:kazumaxneo:20200206115145p:plain

 

Computeボタンを押すと、指定したメールアドレスに自動でメールが届く。メール中のリンクをクリックするとジョブがサブミットされる。

f:id:kazumaxneo:20200206121215p:plain

 

出力

f:id:kazumaxneo:20200206205847p:plain

 

f:id:kazumaxneo:20200206205757p:plain


サーバーからデータが消える前に以下のツールを走らせておくと良いと思います。結果を要約してくれます。

 

 

追記 

この記事と同じタイミングでRyanさんがblastKOALAについて説明されています。

リンク先の説明のように、Procaryotesの自動アノテーション結果をクエリに使用して、パスウェイ解析やfunctional annotation を行えます。私はメタゲノムやDe novo transcriptome解析にも利用しています。 

  

引用
KAAS: an automatic genome annotation and pathway reconstruction server

Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M

Nucleic Acids Res. 2007 Jul;35(Web Server issue):W182-5. Epub 2007 May 25

 

参考

隠れマルコフモデルと適応的閾値を用いた KEGG オーソログ予測法の開発とウェブツールの構築

https://www.genome.jp/tools/kofamkoala/Takuya_Aramaki_Master_Thesis_2019.pdf

 

パスウェイデータベースの紹介とKEGG PATHWAYの使い方@AJACS安芸 

動画と下の方にあるKAASの説明を見てください。

 

SPAdesの出力をフィルタリングする CVLFilter

 

 イルミナのシーケンシングは酵母ゲノミクスに革命をもたらし、現在、市販のドラフトゲノムシーケンシングの価格は200ドル未満になった。人気のあるSPAdesアセンブラにより、あらゆる酵母種のde novoゲノムアセンブリを簡単に生成できる。ただし、ゲノムアセンブリを作成することは日常的になっているが、それらに含まれるものを理解することには依然として困難を伴う。ここでは、SPAdesが各scaffoldsの長さとカバレッジについて提供する情報をグラフ化して、アセンブリの性質を調査し、考えられる問題を診断する方法を示す。ミトコンドリアDNA、リボソームDNA、および酵母プラスミドに由来するscaffoldsは、その高いカバレッジによって識別できる。 multiplexシーケンシングでの他のサンプルからの相互汚染は、カバレッジの低さで識別できる。イルミナのプロトコルで分子標準として頻繁に使用されるバクテリオファージPhiX174およびLambda DNAに由来するscaffoldsも検出できる。Interspecies hybridsなどのヘテロ接合性の高い酵母ゲノムのアセンブリには、多くの場合、2種類のscaffoldsが含まれる。2つの対立遺伝子が2つの別々のscaffoldsにアセンブリし、それぞれがカバレッジレベルCを有するゲノムの領域と、 単一のscaffoldsにco-assembled され(collapsed)カバレッジレベル2Cを持つ領域に分かれる。 Microsoft ExcelまたはGoogleシートを使用して実行できるCoverage-vs.-Length(CVL)プロットでデータを視覚化すると、ゲノムアセンブリの構造を理解し、異常なscaffoldsまたはコンティグを検出する簡単な方法が提供される。 CVLプロットで特定された汚染配列を除去するためにアセンブリをフィルター処理できるPythonスクリプトを提供する。

 

 

f:id:kazumaxneo:20200206032016p:plain

CVL plot of the assembly from Torulaspora delbrueckii strain L17, which is highly heterozygous.  論文より転載

 

インストール

condaでpython2.7の仮想環境を作ってテストした(macos10.14、anaconda3.7使用)。

Github

Git clone https://github.com/APDLS/CVLFilter.git
cd CVLFilter/

#python2環境で実行
conda create -n python2 python=2.7
conda activate python2
python CVLFilter.py

 

実行方法

spades出力のscaffoldsを選抜する。対話形式でファイルを指定して進める。

python CVLFilter.py

$ python CVLFilter.py -h

Enter the name of the scaffolds/contig input file: 

 

 

フルパスでscaffolds.fastaを指定

Enter the name of the scaffolds/contig input file: /data/spades_output/scaffolds.fasta  


次に最低サイズ(bp)を 指定

Enter the minimum contig length to retain: 300

 

最後に最低カバレッジを 指定

Enter the minimum contig coverage to retain: 10

done

Working...

Done

Filtered scaffolds output to the file: <open file 'scaffolds_filtered.fasta', mode 'w' at 0x1057a38a0>

scaffolds_filtered.fastaが出力される。

 

coverageとlenghのplot:CVLplotを出力するRscriptも提供されています。Githubで確認して下さい。

引用

Coverage-Versus-Length Plots, a Simple Quality Control Step for de Novo Yeast Genome Sequence Assemblies

Douglass AP, O'Brien CE, Offei B, Coughlan AY, Ortiz-Merino RA, Butler G, Byrne KP, Wolfe KH

G3 (Bethesda). 2019 Mar 7;9(3):879-887

 

 

カバレッジトラックを視覚化する SparK

2020 3/1 コマンド修正

2024/04/17 help更新

 

Integrative Genomics Viewer(IGV)やUCSCゲノムブラウザなど、NGSデータの表示に利用できる洗練されたリソースが存在するが、領域のエクスポートとpublication用の図の組み立ては依然として困難である。特に、トラックの外観のカスタマイズとtrack replicatesのオーバーレイは、手動で時間のかかるプロセスである。ここでは、RNA-seq、ChIP-seq、ATAC-seqなど、NGSベースのトラックからpublication可能な高解像度のベクターグラフィックフィギュアを自動生成するツールであるSparKを紹介する。 SparKの新しい機能には、replicatesの平均化、標準偏差トラックのプロット、大幅に変更された領域の強調表示が含まれる。 SparKはPython 3で記述されており、主要なOSプラットフォームで実行可能である。コマンドラインプロンプトを使用して図を生成すると、後で変更するのが非常に簡単になる。たとえば、プロットのゲノム領域を変更する必要がある場合、またはトラックを追加/削除する必要がある場合、すべてを手動で再エクスポートおよび再組み立てするプロセスなしで、数秒以内に図を簡単に再生成できる。 SparKでプロットした後、出力SVGベクターグラフィックファイルの変更は、テキスト、線、色など、簡単に行える。 SparKはGitHubで公開されている:https://github.com/harbourlab/SparK

 

 

インストール

macos10.14のanaconda3.7環境でテストした。

依存

  • Python 3 (in theory Python 2 should work too, but I highly reccomend using Python 3)
  • numpy
#bedgraphを作るにはdeeptoolsも必要。
conda install -c bioconda deeptools

本体 Github

git clone https://github.com/harbourlab/SparK.git
cd SparK/

python SparK.py -h

$  python SparK.py -h

usage: SparK.py [-h] [-pt PLOT_TYPE] [-ps SHOW_PLOTS] -pr REGION -cf CONTROL_FILES [CONTROL_FILES ...] [-tf TREAT_FILES [TREAT_FILES ...]] [-cg CONTROL_GROUPS [CONTROL_GROUPS ...]] [-tg TREAT_GROUPS [TREAT_GROUPS ...]] [-gl GROUP_LABELS [GROUP_LABELS ...]]

                [-l LABELS [LABELS ...]] [-gs GROUP_AUTOSCALE] [-es EXCLUDE_FROM_GROUP_AUTOSCALE [EXCLUDE_FROM_GROUP_AUTOSCALE ...]] [-eg EXCLUDE_GROUPS [EXCLUDE_GROUPS ...]] [-f FILLS [FILLS ...]] [-cs CUSTOM_Y_AXIS_SCALES [CUSTOM_Y_AXIS_SCALES ...]]

                [-dc DISPLAY_CHROM_LOCATION] [-gtf GTFFILE] [-tss DRAWTSS] [-genestart DRAW_GENESTART] [-sp SPARK] [-sc SPARK_COLOR [SPARK_COLOR ...]] [-sm SMOOTHEN] [-o OUTPUT_NAME] [-bed BED_FILES [BED_FILES ...]] [-bedcol BED_COLOR [BED_COLOR ...]]

                [-bedlab BED_LABELS [BED_LABELS ...]] [-w TRACK_WIDTH] [-dg DISPLAY_GENES [DISPLAY_GENES ...]] [-dt DISPLAY_TRANSCRIPTS [DISPLAY_TRANSCRIPTS ...]] [-wg WRITE_GENENAMES] [-scale DISPLAY_SCALEBAR]

 

SparC_args

 

optional arguments:

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

  -pt PLOT_TYPE, --plot_type PLOT_TYPE

                        choices: standard, STD, sine

  -ps SHOW_PLOTS, --show_plots SHOW_PLOTS

                        choices: all, averages

  -pr REGION, --region REGION

                        example: chr1:1647389-272634

  -cf CONTROL_FILES [CONTROL_FILES ...], --control_files CONTROL_FILES [CONTROL_FILES ...]

                        separate by space

  -tf TREAT_FILES [TREAT_FILES ...], --treat_files TREAT_FILES [TREAT_FILES ...]

                        separate by space

  -cg CONTROL_GROUPS [CONTROL_GROUPS ...], --control_groups CONTROL_GROUPS [CONTROL_GROUPS ...]

                        group numbers separate by spacse

  -tg TREAT_GROUPS [TREAT_GROUPS ...], --treat_groups TREAT_GROUPS [TREAT_GROUPS ...]

                        group numbers separate by space

  -gl GROUP_LABELS [GROUP_LABELS ...], --group_labels GROUP_LABELS [GROUP_LABELS ...]

                        set group labels

  -l LABELS [LABELS ...], --labels LABELS [LABELS ...]

                        set labels for controls and treatment

  -gs GROUP_AUTOSCALE, --group_autoscale GROUP_AUTOSCALE

                        set to "yes" to autoscale all groups(plots), except the groups excluded in -eg

  -es EXCLUDE_FROM_GROUP_AUTOSCALE [EXCLUDE_FROM_GROUP_AUTOSCALE ...], --exclude_from_group_autoscale EXCLUDE_FROM_GROUP_AUTOSCALE [EXCLUDE_FROM_GROUP_AUTOSCALE ...]

                        group numbers of groups to be excluded from autoscale

  -eg EXCLUDE_GROUPS [EXCLUDE_GROUPS ...], --exclude_groups EXCLUDE_GROUPS [EXCLUDE_GROUPS ...]

                        Exclude groups from the analysis

  -f FILLS [FILLS ...], --fills FILLS [FILLS ...]

                        track colors. One or two colors in hex format for control and treatment tracks

  -cs CUSTOM_Y_AXIS_SCALES [CUSTOM_Y_AXIS_SCALES ...], --custom_y_axis_scales CUSTOM_Y_AXIS_SCALES [CUSTOM_Y_AXIS_SCALES ...]

                        Enter y-axis values for all groups(plots). All groups need a value. Enter "D" for each group no value should be assigned, e.g. to keep autoscaling functionality for some groups

  -dc DISPLAY_CHROM_LOCATION, --display_chrom_location DISPLAY_CHROM_LOCATION

                        set to "no" if you do not want to plot the chromosomal coordinates

  -gtf GTFFILE, --gtffile GTFFILE

                        link gtf file for drawing genes here

  -tss DRAWTSS, --drawtss DRAWTSS

                        set to "no" if TSS sites should not be indicated

  -genestart DRAW_GENESTART, --draw_genestart DRAW_GENESTART

                        set to "yes" if TSS sites should be indicated

  -sp SPARK, --spark SPARK

                        display significant change "yes"

  -sc SPARK_COLOR [SPARK_COLOR ...], --spark_color SPARK_COLOR [SPARK_COLOR ...]

                        spark color

  -sm SMOOTHEN, --smoothen SMOOTHEN

                        smoothen tracks, int

  -o OUTPUT_NAME, --output_name OUTPUT_NAME

                        output graph name, str

  -bed BED_FILES [BED_FILES ...], --bed_files BED_FILES [BED_FILES ...]

                        bed files to be plotted

  -bedcol BED_COLOR [BED_COLOR ...], --bed_color BED_COLOR [BED_COLOR ...]

                        colors of bed files in hex

  -bedlab BED_LABELS [BED_LABELS ...], --bed_labels BED_LABELS [BED_LABELS ...]

                        set labels for bed tracks

  -w TRACK_WIDTH, --track_width TRACK_WIDTH

                        width of the track, default = 150, int

  -dg DISPLAY_GENES [DISPLAY_GENES ...], --display_genes DISPLAY_GENES [DISPLAY_GENES ...]

                        genes to display from the gtf file

  -dt DISPLAY_TRANSCRIPTS [DISPLAY_TRANSCRIPTS ...], --display_transcripts DISPLAY_TRANSCRIPTS [DISPLAY_TRANSCRIPTS ...]

                        display custom transcripts. By default, all transcripts annotated in the gtf file will be merged and displayed as one gene. Alternatively all can be plotted seperatelly by setting this to "all". Further, Transcript IDs can be listed to

                        plot only certain transcripts

  -wg WRITE_GENENAMES, --write_genenames WRITE_GENENAMES

                        write genename instead of transcript ID when transcripts are plotted. Set to "yes".

  -scale DISPLAY_SCALEBAR, --display_scalebar DISPLAY_SCALEBAR

                        set to "no" to remove scalebar

(r4) kamisakakazumanoMac-Studio:SparK kamisakakazuma$ python SparK.py -pr Aggregatilineales_A4b_gen._1_1:1-10000000 -cf ../uesaka.bdg -w 300 -cd 300

usage: SparK.py [-h] [-pt PLOT_TYPE] [-ps SHOW_PLOTS] -pr REGION -cf CONTROL_FILES [CONTROL_FILES ...] [-tf TREAT_FILES [TREAT_FILES ...]] [-cg CONTROL_GROUPS [CONTROL_GROUPS ...]] [-tg TREAT_GROUPS [TREAT_GROUPS ...]] [-gl GROUP_LABELS [GROUP_LABELS ...]]

                [-l LABELS [LABELS ...]] [-gs GROUP_AUTOSCALE] [-es EXCLUDE_FROM_GROUP_AUTOSCALE [EXCLUDE_FROM_GROUP_AUTOSCALE ...]] [-eg EXCLUDE_GROUPS [EXCLUDE_GROUPS ...]] [-f FILLS [FILLS ...]] [-cs CUSTOM_Y_AXIS_SCALES [CUSTOM_Y_AXIS_SCALES ...]]

                [-dc DISPLAY_CHROM_LOCATION] [-gtf GTFFILE] [-tss DRAWTSS] [-genestart DRAW_GENESTART] [-sp SPARK] [-sc SPARK_COLOR [SPARK_COLOR ...]] [-sm SMOOTHEN] [-o OUTPUT_NAME] [-bed BED_FILES [BED_FILES ...]] [-bedcol BED_COLOR [BED_COLOR ...]]

                [-bedlab BED_LABELS [BED_LABELS ...]] [-w TRACK_WIDTH] [-dg DISPLAY_GENES [DISPLAY_GENES ...]] [-dt DISPLAY_TRANSCRIPTS [DISPLAY_TRANSCRIPTS ...]] [-wg WRITE_GENENAMES] [-scale DISPLAY_SCALEBAR]

SparK.py: error: unrecognized arguments: -cd 300

 

実行方法

1、使用にはbedgraphファイルが必要。deeptoolsのbamCoverageコマンドでつくる。 -bs 1を指定するが、chip-seqではこのコマンドは推奨されない(Githubより)。MACS2から出力する。

bamCoverage -b input.bam -o output.bdg -bs 1 -of bedgraph
  • --binSize < INT> bp, -bs     Size of the bins, in bases, for the output of the bigwig/bedgraph file. (Default: 50)
  • -of {bigwig, bedgraph}

 

2、SparKの実行。

bedgraphファイルと領域を指定して実行する。chr1:1-1000のカバレッジを視覚化する。

python SparK.py -pr chr1:1-10000 -cf control.bdg
  • -pr, --region    example: chr1:1647389-272634
  • -cf , --control_files   CONTROL_FILES

f:id:kazumaxneo:20200204231447p:plain


Githubに例があるように、複数ファイルを指定することもできる。

python SparK.py -pr chr:1-10000 -cf sample1.bdg sample2.bdg -gl sample1 sample2

-gl, --group_labels    set group labels

f:id:kazumaxneo:20200301225233p:plain

 

gtfを指定すると、図の下にgenomic featureが表示される。同時に-dgフラグも立てると指定した遺伝子のみ図の下に表示できる(表記が重なるのを防ぐ)。

python SparK.py \
-pr chr12:6520512-6640512 \
-cf HepG2_H3K27AC_1_ENCFF495QSO.bigWig.bdg HepG2_H3K27AC_2_ENCFF348RLL.bigWig.bdg HepG2_H3K4me3_1_ENCFF699DRO.bigWig.bdg HepG2_H3K4me3_2_ENCFF400FYO.bigWig.bdg \
-gff gencode.v24.primary_assembly.annotation.txt \
-gl H3K27AC H3K4me3 H3K27AC-2 H3K4me3-2 \
-dg GAPDH IFFO1 NOP2 CHD4 LPAR5
  • -dg    genes to display from the gff file
  • -gff    link gff file for drawing genes here

 

オーバーレイ表示もできる。どれとどれをオーバーレイ表示するか-tg <NUM>,  -cg <NUM>で指定する。幅を150から300に増やす。出力名はtest.svgにする。

python SparK.py -pr chr:1-10000 \
-cf sample1_cont.bdg -tf sample1_treat.bdg\
-tg 1 -cg 1 \
-gl sample1 \
-w 300 \
-o test.svg

f:id:kazumaxneo:20200301232028p:plain

カバレッジに差があるサンプルを使ってしまったので分かりにくいが、下の方の赤くなっている部分がカバレッジの少ない方のサンプル。そのほか、カバレッジ変動が激しい領域をだけ色を変えて強調表示するなどがある。

 

追記

y軸の値を指定する(新機能)。

 python SparK.py -pr chr1:1-10000 -cf sample1_cont.bdg -w 300 -cs 300
  • -cs     Enter y-axis values for all groups(plots). All groups need a value. Enter "D" for each group no value should be as

 

更新についてはGIthubを確認してください。

引用

SparK: A Publication-quality NGS Visualization Tool

Stefan Kurtenbach, J. William Harbour

bioRxiv preprint first posted online Nov. 16, 2019

 

関連


 

 

 

HMMを使用してメタゲノムの遺伝子を見つけるwebサービス MetaHMM

 

大規模な臨床および環境のメタゲノムデータセットの高速で手頃なシーケンスにより、医療およびバイオテクノロジーのアプリケーションにおける新しい視野が開かれている。今日、地球上の微生物の約1%しか記述できていないと考えられており、メタゲノム解析はサンプル中のほとんどが未知の種を対象としている。極限環境の微生物群集は、高いバイオテクノロジーの可能性を秘めた遺伝子を含んでいる可能性があり、疾患に関連する臨床メタゲノムは、未知の病原体および既知の疾患の病理学的メカニズムを明らかにする可能性がある。サンプル内の分類群の種レベルの識別は今日可能ではないようであるが、隠れマルコフモデル( HMM)のような、人工知能ツールを含む多数の手法を使用して、これらのサンプルで既知の機能を持つ新規遺伝子を検索できる。ここでは、検索対象の遺伝子の相同性に基づいた自動モデル構築が可能で、メタゲノムで最も近い一致を見つけることができる、使いやすいWebサーバーMetaHMMについて説明する。 Webサーバーは、すでに非常に成功しているビルディングブロックを使用する:Clustal Omegaを適用することにより複数のアライメントを実行し、hmmbuildのHMMERコンポーネントを使用して隠れマルコフモデルを構築し、メタゲノム内の指定されたモデルに類似するシーケンスを見つけるためにhmmsearchを使用する。ウェブサーバーはhttps://metahmm.pitgroup.orgで公開されている。

 

 

使い方

http://pitgroup.org/metahmm/ にアクセスする。

f:id:kazumaxneo:20200203211101p:plain

 

入力はUniProt のaccession numbersになる。MetaHMMは、Clustal Omegaを使って入力配列をアラインし、hmmbuildを使ってこれらのタンパク質のモデルを構築、相同なタンパク質をメタゲノムから探索する。

f:id:kazumaxneo:20200203212136p:plain

ここでは右上の"Use the example model"ボタンをクリック。

 

UniProt accession numbers がスペース区切りで入力された。

f:id:kazumaxneo:20200203211151p:plain

 

メタゲノムのデータセットを選択する。メタゲノムデータはiMicrobeポータルに由来する。

f:id:kazumaxneo:20200203220711p:plain

解析にはしばらく時間がかかる。

 

出力

f:id:kazumaxneo:20200203220929p:plain

output_all.csv

f:id:kazumaxneo:20200203221033p:plain

出力について(README)

f:id:kazumaxneo:20200203221133p:plain

 

ユーザーカスタムのメタゲノムを検索することも可能。独自のHTTPまたはFTPストレージを用意し(ファイルホスティングサイトへのリンクは不可)、そこにFASTAファイル(fastqは不可)としてアップロードし、直接リンクを提供する。最大サイズは非圧縮で1GBになっている。

引用
MetaHMM: A webserver for identifying novel genes with specified functions in metagenomic samples

Szalkai B, Grolmusz V.

Genomics. 2019 Jul;111(4):883-885

 

CAARS

 

 大規模なRNAシーケンス(RNA-Seq)は、ゲノムシーケンスの実用的な代替手段として、特に比較分析のために非モデル種でよく使用される(Ozsolak and Milos、2011; Todd et al、2016; Wang et al 、2009)。しかし、トランスクリプトームアッセイのショートリードからの完全長の転写配列へのアセンブリは、反復領域、発現レベル変動、選択的スプライシング、シーケンシングエラー、および組成バイアスに関連する困難な問題を提起する(Garber et al、2011)。さらに、これらの配列を遺伝子ファミリーにクラスター化すること、それらのアラインメント、および遺伝子ツリー再構築のステップはすべて、合意された基準が無い、比較ゲノミクス研究が直面する課題を表している。
 近縁種のゲノムデータの存在に応じて、転写産物のアセンブリにさまざまな戦略を使用できる(Conesa et al、2016; Ockendon et al、2016)。シーケンスされたゲノムを持つ姉妹種が利用できない場合、リードは重複するシーケンスに基づいてデノボでアセンブルされる(例えばTrinity、 Grabherr et al(2011))。それ以外の場合は、ゲノムガイドアセンブリが使用できる(例:Tophat、Trapnell et al (2009)、 Cufflinks、Trapnell et al(2010))。その場合、リードはこのガイドゲノムに合わせて行われ、ローカルの転写産物アセンブリに使用されるリードのクラスターが作成される。この戦略は明らかに、非常にclosely relatedな種、すなわち種を超えてリードのマッピングが可能な場合に限定される。より遠い関係にある種では、RNA-Seqアセンブリに対するアプローチは提案されていないが、開発が続いている。特に、Johnson et al(2013)によるTarget Restricted Assembly Method (TRAM)  、自動のaTRAM(Allen et al、2015)では、BLAST(Camacho et al、2009)を使用してリファレンスゲノムとの配列類似性によってリードが繰り返しコレクションされ、その後遺伝子配列が再構築、アセンブリされる。 k-mersアプローチに基づいてKollectorで異なる実装が提案された(Kucuk et al、2017)。これらの方法は有望な結果を示しているが、RNA-Seqデータおよび一度に数千の遺伝子に使用するようには設計されていない。
 アセンブリ後、転写産物には遺伝子名のアノテーションを付けるのが理想的である。一般に、トランスクリプトームアノテーションは、すでにアノテーションが付けられている種のトランスクリプトームとトランスクリプトームの間の配列類似性に基づいている。このステップは、通常、種固有のduplicationsを処理できないBLAST(Camacho et al、2009)を使用して、Reciprocal Best Hits(RBH)(Rivera et al、1998)によって処理されることがほとんどである(Altenhoff and Dessimoz、2009; Tekaia、2016)。多くの遺伝子が重複しているため、これは問題である。たとえば、Ensemblデータベース(Herrero et al、2016; Yates et al、2016)では、すべてのヒト遺伝子の10%がマウス遺伝子と1対1のorthology関係を持たない。
 原則として、アノテーションをRBHではなく遺伝子系統に依存することで、複雑な相同関係を処理することができる(Chen et ak、2007; Kristensen et al、2011; Kuzniarアノテーション、2008 et al; Tekaia、2016)。このようなアプローチを取ることを勧める。アノテーション付きの転写産物からの遺伝子は、de novo(Kristensen et al、2011)または既存のファミリー[EnsemblCompara(Herrero et al、2016)、TreeFam(Finn et 、2014)、Hogenom(Penel et al、2009)、PhylomeDB(Huerta-Cepas et al、2014)]から相同な遺伝子ファミリーにクラスター化される。次に、再構築された転写産物は、配列の類似性に基づいて遺伝子ファミリーに統合される。これらの拡大された遺伝子ファミリーのアライメントとツリーが再構築される。ツリーの品質は、種のツリーによって提供される情報を使用する再構築方法を使用することで改善できる(Boussau et al、2013; Ullah et al、2015)。最後に、遺伝子ツリーを種のツリーと調整して、種分化、重複、および遺伝子の欠失(loss)にアノテーションを付ける(Kristensen et al、2011)。遺伝子進化のシナリオに基づいて、オルソロジーとパラロジーの関係が導き出され、遺伝子名がアノテーション付きの配列から新しい転写産物に伝播される(Kristensen et al、2011)。このアプローチでは、正確なアノテーションは正確な遺伝子ツリーの結果である。
 ここでは、CAARSという名前の自動化ツールを紹介する。CAARSはよく似た生物、またはより距離のある生物の情報を元にガイドアセンブリしたりアノテーションすることで、非モデル生物のトランスクリプトーム全体をアセンブルし、アノテーションを付ける。CAARSはリファレンス遺伝子のアライメントに依存しており、下流の比較分析に直接使用できる高品質の系統樹とオルソロジー関係を持つ相同遺伝子セットを出力する。 CAARSは、確立されたパイプラインを使って、トランスクリプトームの完全性、トランスクリプトの正確性、およびアノテーションの正確性を改善する。その高品質の出力遺伝子ツリーのおかげで、CAARSは、回復できるオルソログの数の点でEnsemblComparaを改善している。

 

f:id:kazumaxneo:20200121222505p:plain

CAARS overview.   論文より転載

 

wiki

https://github.com/CarineRey/caars/wiki/Tutorial#2-running-caars-example-on-mouse-test-data

  

インストール

Github

#dockerイメージが提供されている
docker pull carinerey/caars

 

実行方法

ランにはRNA seqのリードのほか、サンプルシート、ツリーファイル、遺伝子ファミリーのアラインメントファイルを使う。

caars --outdir OUTPUT_DIR --sample-sheet sample_sheet.tsv --species-tree /home/user/data/species_tree.nw --alignment-dir GENE_FAMILIES_MSA_DIR --seq2sp-dir GENE_FAMILIES_SEQ2SP_DIR --np 2 --memory 5 

 

 作成中

 

 

引用
CAARS: comparative assembly and annotation of RNA-Seq data
Rey C, Veber P, Boussau B, Sémon M

Bioinformatics. 2019 Jul 1;35(13):2199-2207

 

関連