macでインフォマティクス

macでインフォマティクス

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

tRNAをゲノムやraw fastqから探す tRNAscan-SE 2.0

 2019 5/7 タイトル修正

 

 トランスファーRNA(tRNA)は、遺伝暗号のタンパク質への翻訳者として働き全ての生物に存在している。 tRNA scan-SE(ref.1)は、ゲノム中のtRNA遺伝子を同定およびアノテーションするための最も広く使用されているツールである。 8000以上の引用があり、そのユーザーはRNA生物学者、シークエンシングセンター、データベースアノテーター、そして他の基礎研究者を含む。 UNIXベースのソフトウェアを扱う専門知識を持たない可能性がある研究者の使いやすさを向上させるために、tRNA scan-SE On-line Webサイト(ref.2、3)には、迅速で詳細なtRNA分析が掲載されている。tRNA scan-SEを使用して予測されたtRNAは、何千ものゲノムについてGenomic tRNAデータベース(GtRNA db)(ref.4、5)で利用可能であり、生命の3つのドメインすべてにわたって高品質のtRNAコレクションを閲覧できる。
 オリジナルのtRNA scan-SEの実装は、非常に貴重なRfamデータベース(ref.7)よりも前に、ゲノム内のRNA遺伝子にアノテーションを付けるための共分散モデル(CM)(ref.6)の大規模な使用を開拓した。同じRNAファミリーの構造的にアラインしたメンバーをトレーニングすることにより、共分散モデルは、一次配列と二次構造の両方の情報を統合できるstochastic context-free grammarsを介してRNA保存を捉える。集中的な計算要件にもかかわらず、共分散モデルは、tRNAおよび他の多くの構造化RNAを見出すことにおいて比類のない感度および特異性をもたらす。tRNA共分散モデルとのアラインメントによって、任意の所与の配列からtRNAを検索できる。共分散モデルを構築するために使用されるtRNAのトレーニングセットに応じて、検索は任意のtRNA配列の一般的な検出、またはクレード特異的な特徴を有するtRNAを検出するためのより専門的な検索(例えば真核細胞質型tRNA)または特異的tRNAの種類(例:開始メチオニンtRNA)に利用できる。したがって、tRNA scan-SEは、高品質のトレーニングセットアライメントによってのみ制限され、さまざまな種類のtRNAに対して簡単に「調整」できる。この柔軟性により、以前のバージョンのtRNA scanSで真核生物、バクテリアアーキア、または細胞小器官のtRNAのクレード固有の検索モードが可能になった。同じ柔軟性は、tRNAシーケンストレーニングセットによってのみ制限される、強力な新しいクレードとアイソタイプ固有の検索モードのためのフレームワークを提供する。全ゲノムにわたるCM検索とアライメントに必要な計算負荷を減らすために、オリジナルのtRNA scan-SEは、推定tRNAを同定するための初回通過スクリーニングとして2つの高速で高感度なアルゴリズムを使用した(ref.8,9)。次にプログラムは推定tRNAのみをCMにアラインメントし、アイソタイプやイントロン境界などの重要な情報を同定した。初期の一般的およびドメイン特異的なCMは、ゴールドスタンダードのSprinzlデータベースから抽出された1,415のtRNAのアライメントを用いて構築され(ref.10、11)、CMの構築およびアライメントはCOVEを用いて行われた(ref.6)。
 その最初の実施以来、tRNA遺伝子およびtmRNA遺伝子の両方を検出するARAGORN(ref.12)を含む、他の多数の tRNA検出および分類方法が開発されてきた。 DOGMA(ref.13)、ARWEN(ref.14)、MITOS(ref.15)、およびtRNAファインダー(ref.17)は、さまざまなタイプのオルガネラゲノムのtRNAにアノテーションを付けるために設計されている。 SPLITS(ref.18)は、微生物ゲノム中のスプリットおよびイントロンを含むtRNAを見つけるために設計されている。これらの方法はすべてtRNA scan-SEを改善または補完するように設計されており、特にtRNA scan-SEのコア検出ソフトウェアに依存しているものが多くある。
 tRNA scan-SEは、過去20年間tRNA検出のための信頼性が高くアクセスしやすいツールであり続けているが、新しいアルゴリズム、新しいデータ、および性能と機能予測精度の向上のための新しい戦略を組み込むための大きな修正が必要だった。ここでは、次のような機能強化が行われた最新バージョンのtRNA scan-SEについて説明する。(1)Infernal 1.1共分散モデル検索ソフトウェアの統合という形での改良型共分散モデル検索技術(ref.19)。 (2)新たにシーケンシングされた何千ものゲノムからのより広く代表的な多様性のtRNA遺伝子を利用する最新の検索モデル(表1)。 (3)完全な一連のアイソタイプ特異的tRNA共分散モデルからの比較情報に基づくtRNAのより良い機能的分類、および(4)タンパク質で最も使用される可能性が高い真核生物のtRNAを同定するための新しい「高信頼」フィルター翻訳。

 

インストール

 condaで導入できる。macのminiconda3-4.3.30環境でテストした。

mamba install -c bioconda trnascan-se

> tRNAscan-SE -h

$ tRNAscan-SE -h

 

tRNAscan-SE 2.0 (December 2017)

Copyright (C) 2017 Patricia Chan and Todd Lowe

                   University of California Santa Cruz

Freely distributed under the GNU General Public License (GPLv3)

 

 

Usage: tRNAscan-SE [-options] <FASTA file(s)>

 

  Scan a sequence file for tRNAs

   -- default: use Infernal & tRNA covariance models

      with eukaryotic sequences

      (use 'Search Mode Options' below to scan other types of sequences)

 

Search Mode Options:

 

  -E                          : search for eukaryotic tRNAs (default)

  -B                          : search for bacterial tRNAs

  -A                          : search for archaeal tRNAs

  -M <model>                  : search for mitochondrial tRNAs

                                  options: mammal, vert

  -O                          : search for other organellar tRNAs

  -G                          : use general tRNA model (cytoslic tRNAs from all 3 domains included)

  --mt <model>                : use mito tRNA models for cytosolic/mito detemination

                                  (if not specified, only cytosolic isotype-specific model scan will be performed)

  -I                          : search using Infernal

                                  default use with -E, -B, -A, or -G; optional for -O

      --max                   : maximum sensitivity mode - search using Infernal without hmm filter (very slow)

  -L                          : search using the legacy method (tRNAscan, EufindtRNA, and COVE)

                                  use with -E, -B, -A or -G

  -C  --cove                  : search using COVE analysis only (legacy, extremely slow)

                                  default use with -O

  -H  --breakdown             : show breakdown of primary and secondary structure components to

                                  covariance model bit scores

  -D  --nopseudo              : disable pseudogene checking

 

Output options:

 

  -o  --output <file>         : save final results in <file>

  -f  --struct <file>         : save tRNA secondary structures to <file>

  -s  --isospecific <file>    : save results using isotype-specific models in <file>

  -m  --stats <file>          : save statistics summary for run in <file>

                                  (speed, # tRNAs found in each part of search, etc)

  -b  --bed <file>            : save results in BED file format of <file>

  -a  --fasta <file>          : save predicted tRNA sequences in FASTA file format of <file>

  -l  --log <file>            : save log of program progress in <file>

  --detail                    : display prediction outputs in detailed view

  --brief                     : brief output format (no column headers)

 

  -? #                       : '#' in place of <file> chooses default name for output files

  -p  --prefix <label>        : use <label> prefix for all default output file names

 

  -d  --progress              : display program progress messages

  -q  --quiet                 : quiet mode (credits & run option selections suppressed)

  -y  --hitsrc                : show origin of hits (Ts=tRNAscan 1.4, Eu=EufindtRNA,

                                  Bo=Both Ts and Eu, Inf=Infernal)

 

Specify Alternate Cutoffs / Data Files:

 

  -X  --score <score>         : set cutoff score (in bits) for reporting tRNAs (default=20)

  -g  --gencode <file>        : use alternate genetic codes specified in <file> for

                                  determining tRNA type

  -z  --pad <number>          : use <number> nucleotides padding when passing first-pass

                                  tRNA bounds predictions to CM analysis (default=8)

  --len <length>              : set max length of tRNA intron+variable region for legacy search mode

                                  (default=116bp)

Misc Options:

 

  -h  --help                  : print this help message

  -c  --conf <file>           : tRNAscan-SE configuration file (default: tRNAscan-SE.conf)

  -Q  --forceow               : do not prompt user before overwriting pre-existing

                                  result files  (for batch processing)

 

  --match <EXPR>              : search only sequences with names matching <EXPR> string

                                  (<EXPR> may contain * or ? wildcard chars)

  --search <EXPR>             : start search at sequence with name matching <EXPR> string

                                  and continue to end of input sequence file(s)

Special Advanced Options (for testing & special purposes)

 

  -U                          : search for tRNAs with alternate models defined in configuration file

 

  -t  --tscan                 : search using tRNAscan only (defaults to strict params)

  --tmode <mode>              : explicitly set tRNAscan params, where <mode>=R or S

                                  (R=relaxed, S=strict tRNAscan v1.3 params)

 

  -v  --verbose <file>        : save verbose tRNAscan 1.3 output to <file>

  --nomerge                   : Keep redundant tRNAscan 1.3 hits (don't filter out multiple

                                  predictions per tRNA identification)

  -e  --eufind                : search using Eukaryotic tRNA finder (EufindtRNA) only

                                  (defaults to Normal seach parameters when run alone,

                                  or to Relaxed search params when run with Cove)

  --emode <mode>              : explicitly set EufindtRNA params, where <mode>=R, N, or S

                                  (relaxed, normal, or strict)

 

  --iscore <score>            : manually set "intermediate" cutoff score for EufindtRNA

  -r  --fsres <file>          : save first-pass scan results from EufindtRNA, tRNAscan, or

                                  Infernal hmm in <file> in tabular results format

  --mid                       : fast scan mode - search using Infernal with mid-level strictness of hmm filter

  -F  --falsepos <file>       : save first-pass candidate tRNAs in <file> that were then

                                  found to be false positives by second-pass analysis

  --missed <file>             : save all seqs that do NOT have at least one

                                  tRNA prediction in them (aka "missed" seqs)

  --thread <number>           : number of threads used for running infernal (default is to use available threads)

 

 

 

実行方法

genomeのfastaかraw fastqを選択する。

tRNAscan-SE input_genome.fa -o output
  • -E     search for eukaryotic tRNAs (default)
  • -B     search for bacterial tRNAs
  • -A     search for archaeal tRNAs
  • -M <model> : search for mitochondrial tRNAs
  • -O     search for other organellar tRNAs
  • -o     save final results in <file>

 

web版の使い方

http://trna.ucsc.edu/tRNAscan-SE/にアクセスする。

f:id:kazumaxneo:20190502125433p:plain

Downloadからはソースコードをダウンロードできる。

 

Sequence sourceから種を選択する。各検索モードは選択された種または系統発生グループからのtRNAトレーニングデータに基づいている。

f:id:kazumaxneo:20190503193342p:plain

 

Search modeから使用する確率モデルを選択する。標準の検索モードはほとんどの場合高速で十分にsensitiveだが、極端に異なるtRNAを含むショートリードの場合、「Infernal without HMM filter」モードにするとわずかに感度が向上する。

f:id:kazumaxneo:20190503193452p:plain

 

ファイルを選択ボタンから、ゲノム配列のfasta、またはraw fastqを指定する。

f:id:kazumaxneo:20190503194113p:plain

必要に応じてadvanced user向けパラメータを設定する。

f:id:kazumaxneo:20190503194739p:plain

Run tRNAscan-SEボタンを押して検索を実行する。

 

検出されたtRNAが表示される。

f:id:kazumaxneo:20190503194245p:plain

 

予測された各tRNAは、アイソタイプ特異的モデルでさらに評価される。

f:id:kazumaxneo:20190503195120p:plain

 

最後にRun statisticsがまとめられる。

f:id:kazumaxneo:20190503195133p:plain

f:id:kazumaxneo:20190503195143p:plain

f:id:kazumaxneo:20190503195149p:plain

f:id:kazumaxneo:20190503195151p:plain


引用

tRNAscan-SE 2.0: Improved Detection and Functional Classification of Transfer RNA Genes
Patricia P. Chan, Brian Y. Lin, Allysia J. Mak,  Todd M. Lowe

bioRxiv preprint first posted online Apr. 30, 2019

 

関連


複数データベースを統合した包括的な薬剤耐性遺伝子データベース ARGminer

2019 5/4 タイトル修正

 

 薬剤耐性(AMR)は、世界保健機関(WHO)によって世界規模の主要な健康上の脅威として認識されている。 AMRは2050年までに指数関数的に増加し、実質的なヒトの罹患率と死亡率をかなり増やすと予測されている。したがって、モニタリングを強化してAMRの拡散に取り組むためには、迅速な行動が必要である。これには以下のものが含まれる。環境源および経路を介して抗生物質耐性遺伝子(ARG)の蔓延を制御するメカニズムを理解する。 診療所で問題となる前に新しいARGを発見し、ARGアノテーションの新しい計算戦略を開発し、現在のARGリポジトリを拡張する。
 メタゲノムシーケンシングは、様々な環境の多様なARG、すなわち「resistomes」にアクセスする強力な手段を提供して新しいARGの発見とその相互作用をサポートしている。既存のメタゲノムアプローチは、配列類似性計算を通して抗生物質耐性属性を予測する方法に大きく依存しており、これは大きな制限を受ける。第一に、そのような類似性計算は、一貫した正確なARG識別を可能にする高品質で最新のARGリファレンス/アノテーションデータベースを必要とする。第二に、そのような分析の範囲は、シーケンスアラインメントにおいて使用されるパラメータカットオフの厳密さまたはアラインメントのための包括的な標的遺伝子の欠如のために、以前に特徴付けられたARGに限定される。

 所与のサンプル中に存在する全範囲のARGを広く正確に検出するためのメタゲノムベースのアプローチの能力を向上させるには、対応するデータベースのキュレーションを継続的に拡張し向上させることが必要である。しかしながら、偽陽性、すなわち、必ずしもAMR表現型を誘発しないARG様遺伝子を組み込むリスクは、拡大されたキュレーション努力に対する大きな障害となっている。したがって、手動による検査と潜在的なARGエントリの検証は、AMRデータベースとそのアプリケーションの有効性を保証するための重要な側面である。
 ARGの手動キュレーションは、一般に、公共データベースの維持管理に取り組む研究グループに関連する少数の専門家によって行われる。このプロセスは複雑で面倒で時間がかかる。例えば、抗生物質耐性データベース(ARDB)の最後の更新は2009年に行われ、したがって、bla NDM-1やmcrのような新しく発見されたARGは含まれていない。 ARGアノテーションの構成を簡素化するために設計されたMEGAResデータベースは、2016年12月から更新されていない。Mobile Genetic elements (MGEs)を介して伝達された証拠がある遺伝子を含むresquデータベースは、2013年以降更新されていない。包括的な抗生物質耐性データベースCARDは、最新のARGリソースであると広く考えられている。 2016年に最初に導入されたCARDは、2018年10月の最新リリース以来、合計21回更新されており、ARG配列およびメタデータ(例えば、抗生物質クラス、遺伝子名、およびメカニズム)に対応する変更が加えられている。これは、その分野の専門家であっても、ARGデータベースのキュレーションがいかに複雑で時間がかかるかを示している。
 命名法の矛盾や一塩基多型、ハウスキーピング遺伝子の排除などの現在利用可能なデータベースの制限に対処するため、強力なstructured antibiotic resistance database (SARG)などの試みがなされている。本著者ら自身の研究グループでは、以前にディープラーニングを使用してARGを予測するための計算アプローチDeepARGを導入した。機械学習モデルとともに、DeepARG-DBという名前のキュレーションデータベースもリリースした。このデータベースは、手動キュレーション、ARGの文献レビュー、およびシーケンスアラインメントを用いたARGのアノテーションを採用している。 DeepARG-DBは、2017年7月に最初にリリースされ、2018年8月に最後に更新された。しかし、deepARGデータベースは複数のリソースからのアノテーションに依存しているため、他のデータベースからのエラーの伝播に敏感である。これは、新しいARGを簡単に統合したり、現在のARGのアノテーションを検証したりするために、さまざまなリソースからすべてのARG情報を取り込むことができる特別なツールを有効にする必要性を強調している。

 膨大な数のARGのキュレーションおよび手動検証における困難を克服するために、この複雑なタスクをより単純でより小型のマイクロタスクに分解する新規なアプローチが提案されている。この方法論の核心は、AMRリソースの大要を集約すること、および非専門家、すなわち一般大衆およびドメイン専門家が集合的にARGデータベースのキュレーションを実行できるようにARG情報を単純化するクラウドソーシング戦略を展開することからなる。生物学におけるクラウドソーシングの適用、特にデータキュレーションのための適用は、新しいものではなく、患者のオンライン投稿からの医学的に関連するtermの特定、PubMedで説明されている疾患へのアノテーション、データベースおよび体系的な薬物適応症の調査、生物医学オントロジー、および遺伝子 - 疾患相互作用、など様々な分野を含む。興味深いことに、ほとんどの研究で、クラウドソーシングはエキスパートキュレーションと同じくらい効果的であることが証明されている。
 すべてのARGリソースを網羅する大きな問題は、標準化された遺伝子命名法の欠如である。特に、ARGの命名は細菌遺伝子の命名の一般的な命名法には従っていない。例えば、マクロライド耐性遺伝子は、クラスが括弧で示されるように構造化されている (e.g.,​ole​(B),​srm(​B),​vga​(B)or​ere​(B))。テトラサイクリン遺伝子と比較すると、この遺伝子命名法は根本的に異なる。なぜなら、テトラサイクリン遺伝子では、決定基が遺伝子名の後に大文字として置かれているからである (e.g., tetA, tetB, tetC) 。同時に、それらの命名法は、ベータラクタマーゼ遺伝子にアノテーションをつけるために提案された遺伝子規約とは異なる。(一部略)したがって、ARGの命名規則および命名規則の多様性とバリエーションは、一貫性のあるARGキュレーションを複雑にし、大きく妨げる。

 ここではARGの手動キュレーションを強化するオンラインプラットフォームであるARGminerを紹介する。 ARGminerを使用すると、CARD、DeepARG-DB、ARDB、MEGARes、UniProt、NDARO、SARG、ResFinder、およびARGANNOTなど、いくつかのARGリソースから入手可能なすべての情報を整理および取得できる。手動のクラウドソースベースのキュレーションは、検証を助け、一貫性を達成するために自然言語処理NLP)で広く使用されている技術であるword embeddingsに基づく機械学習モデルによって強化されている。一方、プラスミドファージやウイルスのようなMGEは、ARGの普及に重要な役割を果たしている。したがって、ARGminerはまた、それぞれ病原体またはMGEによるARGの潜在的な伝播についての情報を提供するPATRICおよびACLAMEデータベースともインターフェースする。 ARGminerプラットフォームは、AMRの拡大と闘うことに貢献するという共通の欲求に動機付けられた研究者と市民の幅広いコミュニティによるARGアノテーション標準化のための協調的かつ統合的なアプローチを促進するオープンソースプロジェクトとして設計、構築、実装されている。 ARGminerには、サイエンスコミュニティがARGデータベースの最新の更新と開発に積極的に取り組んでいくことを目的として、ユーザーが質問を投稿したり、薬剤耐性に関する解決策や議論を共有したりするためのコミュニティブログもある(補足図1を参照)。 ARGminerに関連するすべてのデータは、ソースコードと同様に、http://bench.cs.vt.edu/argminerの公開リポジトリから無料で入手できる。

 

Instructions

https://bench.cs.vt.edu/argminer/#/forum/selected_question;id=1541191047351

Blog

https://bench.cs.vt.edu/argminer/#/forum

 

Github


https://bench.cs.vt.edu/argminer/#/classify;gene_id=A0A0Q7HDK9 にアクセスする。

f:id:kazumaxneo:20190502203239p:plain

薬剤耐性遺伝子の検索

topから検索できる。catIと検索してみた。

f:id:kazumaxneo:20190502220556p:plain

 

該当するhitが表示される。

f:id:kazumaxneo:20190502220655p:plain

一番上のARDB (5) CARD (1)ヒットのchloramphenicolのcatIをクリックしてみる。

 

ARDBのhitはクラスchloramphenicolでSimilarity(sequence similarity)100、coverage(alignment coverage)100のcata1となっていた。

f:id:kazumaxneo:20190502220912p:plain

 

その下にはSuggested Gene Nomenclatureとの合致率も表示される。

f:id:kazumaxneo:20190502221829p:plain

 

さらに下には各データベースでの該当遺伝子の説明が記載されている。

f:id:kazumaxneo:20190502221301p:plain

 

その下にはアミノ酸配列も表示される。

f:id:kazumaxneo:20190502222023p:plain

 

さらに下にはPATRICからの、病原性バクテリアゲノムでの存在数が表示される。

f:id:kazumaxneo:20190502222050p:plain

1828 genomes contain this particular gene (P58777). From those, 894 ( 48.9% ) genomes are labeled as pathogens.

 894ゲノム ( 48.9% ) 保有となっている。

 

一番下にはACLAMEからのMGEs(wiki)とともに存在しているエビデンスがまとめられる。

f:id:kazumaxneo:20190502222427p:plain

 

 

データベースのダウンロード

現在の最新版は2019-04-24公開のARGminer-v1.1.1となっている。

f:id:kazumaxneo:20190502213943p:plain

 

ARGminerデータベースは2つのfastaファイルとしてリリースされている。

ARGminer-vxx.A.fasta:ARGminerによって検査されたか、他のキュレーションデータベース(CARD、Resfinder、ARG-ANNOT、SARG)で報告されたARGが含まれる。
ARGminer-vxx.B.fasta:このデータベースには、UniProtによる検査が必要な可能性のあるARGを含む、ARGminerからのすべてのARGが含まれる。

多様性を考量すると、Bタイプ利用が推奨されている。

 

リンクを右クリックしてダウンロードする。

f:id:kazumaxneo:20190502214706p:plain

fastaがダウンロードされる。 

 

引用

ARGminer: A web platform for crowdsourcing-based curation of antibiotic resistance genes
G. A. Arango-Argoty​,​ G. K. P. Guron​,​ E. Garner​,​ M.V. Riquelme,​ L. S. Heath​,​ ​ ​A. Pruden,​ P. J. Vikesland,​  L. Zhang

bioRxiv preprint first posted online Mar. 1, 2018

 

効率的なk-merカウンタ kmcEx

 

K-merは、それらの頻度と共に、エラー訂正、リピート検出、マルチプルシーケンスアラインメント、ゲノム構築などの基本的なビルディングブロックとして役立ち、k-merカウントにおける集中的な研究を引き付けた。ただし、k-merカウンタの出力自体は大きい。非常に多くの場合、メインメモリに収まるには大きすぎるため、ユーザビリティが非常に狭くなる。
 本著者らは、k-merならびにそれらの頻度を符号化し、良好なメモリ節約および検索効率を達成する新規なアイデアを紹介する。具体的には、カウントされたk-merを結合ビットアレイ(1つはk-mer表現用、もう1つは頻度符号化用)によって符号化するブルームフィルタのようなデータ構造を提案する。 5つのリアルデータセットでテストした結果、31-merのメモリ節約率の平均は、7つのハッシュ関数を持つ生の入力と比較して13.81と高いことがわかった。同時に、検索時間の複雑さはうまく制御され(事実上一定)、そして偽陽性率は2桁減少する。

 


インストール

ubuntu16.04でbuildしてテストした(docker使用、ホストOS macos10.14)。

依存

kmcEx is based on C++11.

Github

git clone https://github.com/lzhLab/kmcEx.git
cd kmcEx/
make

./kmcEx

# ./kmcEx 

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

           kmcEx: counted k-mer encoding & decoding                   

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

VERSION: 1.3

DATE   : Apr 2nd, 2019

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

 

1. USAGE

     kmcEx [options] <input_file_name> <output_file_name> <working_directory>

     kmcEx [options] <@input_file_names> <output_file_name> <working_directory>

2. OPTIONS

     1) REQUIRED

        input_file_name    - single file in FASTQ format (gziped or not)

        @input_file_names  - file name with list of input files in FASTQ format (gziped or not)

        working_directory  - save temporary files

     2) OPTIONAL

        -k<len>            - k-mer length (default: 31)

        -t<value>          - total number of threads (default: 4)

        -ci<value>         - exclude k-mers occurring less than <value> times (default: 1)

        -cs<value>         - maximal value of a counter (default: 1023)

        -nh<value>         - number of hash (default: 7)

        -nb<value>         - number of bit array (default: 5)

3. EXAMPLES

     kmcEx -k31 -nh7 -nb5  rs.fastq rs.res /tmp

     kmcEx -k31 -nh7 -nb5  @rs.lst rs.res /tmp

 

 

 

実行方法

fastq、出力、作業ディレクトリを指定して実行する。

kmcEx -k 31 -nh 7 -nb 5 -t 8 rs.fastq output /tmp
  • -t      total number of threads (default: 4)
  • -k     k-mer length (default: 31)
  • -nh   number of hash (default: 7)
  • -nb   number of bit array (default: 5)

***********

Stage 1: 100%

Stage 2: 100%

1st stage: 1.59445s

2nd stage: 0.901676s

Total    : 2.49613s

Tmp size : 61MB

 

Stats:

   No. of k-mers below min. threshold :            0

   No. of k-mers above max. threshold :            0

   No. of unique k-mers               :      4523834

   No. of unique counted k-mers       :      4523834

   Total no. of k-mers                :     62516374

   Total no. of reads                 :       254371

   Total no. of super-k-mers          :      5097219

 

kmcEx status:

   Number of hash functions              :  7

   Number of coupled-bit arrays          :  5

   Number of k-mers having count=1       :  669035

   Number of k-mers having count>1       :  3854799

   kmcEx model construction time         :  1.92668

   Memory usage for coupled-bit arrays   :  16MB

   Memory usage for hash_map             :  9MB

   (k-2)-mer BF size for count>1 k-mers  :  1MB

   k-mer BF size for count=1 k-mers      :  0MB

   (k-2)-mer BF size for count=1 k-mers  :  0MB

   Total memory usage                    :  27MB

   No.1 coupled-bit array occupancy rate :  0.337

   No.2 coupled-bit array occupancy rate :  0.339

   No.3 coupled-bit array occupancy rate :  0.340

   No.4 coupled-bit array occupancy rate :  0.341

   No.5 coupled-bit array occupancy rate :  0.333

 

-------

   kmcEx model is successfully saved in  :  /tmp/rs.res

カレントにoutput.kmc_preとoutput.kmc_sufができる。/tmp/outputにモデルファイルができる。

ls /tmp/output/

bit1.bin  bit2.bin  bloom.bin  bloom2.bin  hash.bin  km_back.bin  last_map.bin  param.conf

 

あとはkmcで処理できる。

全k-merを書き出す。

kmc_dump output count
head count #開く

# head count  

AAAAAAAATGGCGATCGCCGCTGTGAGCCAA 15

AAAAAAAATTACCGGGGGGCGATCGCCATTT 23

AAAAAAAATTTTCGGGCGATCGCCATTGATT 18

AAAAAAACTGCCATGGTTGCCCCCAATGAAG 16

AAAAAAAGCTTGATCTCCCACAGCCATGGTT 19

AAAAAAAGGCGATCGCAAAGGCTTTTCCCGC 1

AAAAAAAGGCGATCGCCACAGTCAATGACGA 11

AAAAAAATCCAAGGCGATCGCCTCTGCTACC 22

AAAAAAATGGCGATCGCCGCTGTGAGCCAAT 16

AAAAAAATTACCGGGGGGCGATCGCCATTTT 22

他の例

 

引用

kmcEx: memory-frugal and retrieval-efficient encoding of counted k-mers
Peng Jiang Jie Luo Yiqi Wang Pingji Deng Bertil Schmidt Xiangjun Tang Ningjiang Chen Limsoon Wong Liang Zhao
Bioinformatics, Published: 30 April 2019

 

 

SRAのRNA seqデータを素早く比較・分析する Digital expression explorer 2(手持ちのデータにも対応)

2021 1/9 ツイート追記

 

 10年前の最初の記述以来、RNAシーケンス(RNA-seq)はトランスクリプトームにおける強力な方法となり、非常に正確な遺伝子発現の定量を可能にした[ref.1]。シークエンシングのコストが下がるにつれて、RNA seqのデータは科学文献でより一般的になりつつある。再利用と透明性の向上を目的として、rawデータおよび処理済みファイルの形でこれらのデータをGene Expression Omnibus(GEO)およびSequence Read Archive(SRA)[ref.2、3]に寄託することは当分野における標準的な慣行であり、ジャーナルにとって必須の要件である。しかし実際には、生物学者による広範な再利用を妨げるいくつかの障害がある。まず、SRAのrawシーケンスデータを処理するには、かなりの計算資源とコマンドラインの専門知識が必要になる。第二に、GEOがホストする処理済みRNA-seqデータは、さまざまなソフトウェアツールとゲノムアノテーションセットを利用するさまざまなフォーマットで作成されているため、メタアナリシスが複雑になる。科学界にとってこれらのデータの価値とそれらを生成するための多大なコストにもかかわらず、RNA-seqデータ集約の取り組みはヒトとマウスに大部分限定されているか[ref.4,5]またはクローズドソース/subscriptionサービス[ref.6]になっている。 BgeeDBは、多くの動物種に関するさまざまなライフステージ段階 (excluding disease, treatments, or genetic perturbations)  でのベースラインサンプルの高品質測定に重点を置いて、アレイおよびシーケンスベースの発現データを提供している。 Expression Atlasは、有益なグラフィカルインターフェイスを備えた処理済みの発現マイクロアレイデータの最も包括的なリポジトリの1つだが、現在含まれているのは比較的少数のRNA-seqデータセットだけである[ref.8]。公共のトランスクリプトームデータの再利用を促進するために、本著者らは、多くの種類の下流分析と互換性があり、いくつかの主要な生物の一様に処理されたRNA-seqデジタル遺伝子レベルおよび転写レベルの発現データのオープンアクセスウェブベースリポジトリであるDigital Expression Explorer 2(DEE2)を開発した。

 DEE2は3つの部分からなる。(i)SRAからrawデータセットをダウンロードし処理するパイプライン。 (ii)処理されたファイルが収集され、フィルタリングされ、そして編成/格納され、そしてジョブ待ち行列が生成されるデータレポジトリ。 (iii)ユーザーがメタデータを検索し、興味のあるデータセットを取得できるWebサーバー。 DEE2の編成の概略図が論文図1に提供されている。データ処理ノードはウェブサーバからSRAランアクセッション番号を要求し、SRAからrawデータを取得する。 処理されたデータはWebサーバーに送信され、検証され、DEE2リポジトリサーバーに中継される。 リポジトリサーバはさらなる検証チェックを実行し、新しいデータセットリポジトリに組み込み、SRAdbV2 [ref.9]から対応するメタデータを収集し、未処理のジョブをキューに入れる。 リポジトリサーバーは、更新されたメタデータとジョブキュー情報をWebサーバーに送信する。 エンドユーザーは、Webブラウザコマンドライン、またはバルクダンプからデータを取得する。

 DEE2パイプラインはコンテナ化を使用して迅速なアプリケーション展開を可能にし、さまざまなコンピュータシステムにわたる分析の再現性を保証する。 エンドユーザーは自分のハードウェア上でDockerイメージ[ref.10]を実行して、種名とSRA実行登録で指定された目的のSRAデータセットを処理することができる。 処理が完了すると、ユーザーはすぐに出力にアクセスできるようになり、DEE2リポジトリサーバーによる検証後、データセットは一般に公開されるようになる。 このように、パワーユーザーは、確立された分析パイプラインを使用することによって利益を得ると同時に、公共リソースの拡大に貢献する。 Dockerイメージに関する1つの懸念は、それらが管理者「ルート」許可なしに、例えば共有高性能コンピューティングシステムのユーザによって実行されることができないということである。 この制限に対処するために、イメージをroot権限なしで利用できるSingularity [ref.11]またはUDocker [ref.12]で使用するため変換できる。

 

 

f:id:kazumaxneo:20190429165806p:plain

Overview of RNA-seq data processing, storage, and provision. Githubより転載

 


 

Digital Expression Explorer: RNA-seq analysis in minutes

 

DEE2 data quality metrics

dee2/qc_metrics.md at master · markziemann/dee2 · GitHub

 

GIthub

 

 

1、web版

http://dee2.io にアクセスする。

f:id:kazumaxneo:20190429114714p:plain

 

シロイヌナズナを指定し、GSE63462,SRP070529,SRR401430とタイプした。

f:id:kazumaxneo:20190430173230j:plain

 

19データセット見つかった。IDをクリックするとNCBI SRAにジャンプする。

f:id:kazumaxneo:20190430173331j:plain

上のSelect allをクリックして19ヒット全てにチェックを入れ、下のGet Countsボタンをクリック。

 

QC結果、Gene counts(STAR)およびTranscripts counts(Kallisto)結果が瞬時にダウンロードされる(解析フローは論文の図2参照)。

f:id:kazumaxneo:20190430173529j:plain

 

QC_Matrix.tsvをexcelで開いた。

f:id:kazumaxneo:20190430174019j:plain

一番下の行が評価基準に従って決定されたデータの質。“pass”, “warn” 、 そして“fail”がアサインされている。 評価は複数項目にわたり、このリンク先の基準に従っている。

dee2/qc_metrics.md at master · markziemann/dee2 · GitHub

論文表2も参照。 

 

geneレベルでの定量: GeneCountMatrix.tsv(STAR)

f:id:kazumaxneo:20190430173748j:plain

あっとゆう間ににリードカウントデータが得られた!

 

transcriptsレベルでの定量: TxCountMatrix.tsv(Kallisto)

(省略)

結果の違いは論文の引用ref.22参照(pubmed)。

 

あとは正規化して自由に比較できる。上に貼ったYou tube動画にあるように、オーサーらは取得したリードカウントの分析をDegustで行なうことを提案している(webベースで行うならiDEPも便利です)。

 

 

TopのSearch for keywordsからはキーワード検索もできる。

f:id:kazumaxneo:20190430174435j:plain

hit項目が多すぎる場合(>500)はバルクダウンロードのみ対応する。


 

2、コマンドライン

コマンドライン(コンソール)でも使用できる。できることはweb版と同じだが、コマンドライン版はDEE2に登録されていないデータの解析、手持ちのデータの解析に役に立つ(*1)。web版と同じ以下のモデル生物種に対応している。

  • A. thaliana
  • C. elegans
  • D. melanogaster
  • D. rerio
  • E. coli
  • H. sapiens
  • M. musculus
  • R. norvegicus
  • S. cerevisiae

特徴(Githubより)。

  • Intelligent adapter detection and clipping
  • Clipping of non-reference 5' bases (eg UMIs)
  • Strandedness detection
  • Parallel assignment of reads to genes and transcripts with STAR and Kallisto
  • Thorough quality control logs
  • Open source pipeline
  • Distributed approach using containers link
  • Ability to process own fastq files as well as those from SRA

 

DEE2はdockerイメージとして配布されているので、まずDEE2イメージを取得する。

docker pull mziemann/tallyup

 

例1、SRAのIDを指定してDEE2をランする。ここでは5つ指定。

docker run mziemann/tallyup \
ecoli SRR2637695,SRR2637696,SRR2637697,SRR2637698

 

ジョブが終わったらホストのカレントにコピーする("-alq"は停止したコンテナも含め最新起動のコンテナのIDのみ表示)。

docker cp $(docker ps -alq):/dee2/data/ .

 

data/ecoli/ディレクトリができる。出力を確認。

> ls -lh data/ecoli/

$ ls -lh data/ecoli/

total 744

drwxr-xr-x@ 9 user  staff   306B 10  5  2018 SRR057750

drwxr-xr-x@ 9 user  staff   306B  4 30 15:29 SRR2637695

-rw-r--r--@ 1 user  staff    88K  4 30 15:29 SRR2637695.ecoli.zip

drwxr-xr-x@ 9 user  staff   306B  4 30 15:34 SRR2637696

-rw-r--r--@ 1 user  staff    89K  4 30 15:34 SRR2637696.ecoli.zip

drwxr-xr-x@ 9 user  staff   306B  4 30 15:41 SRR2637697

-rw-r--r--@ 1 user  staff    89K  4 30 15:41 SRR2637697.ecoli.zip

drwxr-xr-x@ 9 user  staff   306B  4 30 15:50 SRR2637698

-rw-r--r--@ 1 user  staff    90K  4 30 15:50 SRR2637698.ecoli.zip

drwxrwxrwx@ 8 user  staff   272B 10  5  2018 SRR5985593_1

-rw-r--r--@ 1 user  staff    29B 10  5  2018 date.txt

サブディレクトリの1つSRR2637698/の中身を確認

> ls -lh data/ecoli/SRR2637698/

$ ls -lh data/ecoli/SRR2637698

total 544

-rw-r--r--@ 1 user  staff    86B  4 30 15:41 SRR2637698.attempts.txt

-rw-r--r--@ 1 user  staff     0B  4 30 15:50 SRR2637698.finished

-rw-r--r--@ 1 user  staff   122K  4 30 15:50 SRR2637698.ke.tsv

-rw-r--r--@ 1 user  staff    30K  4 30 15:50 SRR2637698.log

-rw-r--r--@ 1 user  staff   704B  4 30 15:50 SRR2637698.qc

-rw-r--r--@ 1 user  staff    42K  4 30 15:50 SRR2637698.se.tsv

-rwxr-xr-x@ 1 user  staff    61K  4 30 15:41 volunteer_pipeline.sh

TSVファイル

column -t data/ecoli/SRR2637698/SRR2637698.ke.tsv |head

$ column -t data/ecoli/SRR2637698/SRR2637698.ke.tsv |head

SRR2637698_target_id  SRR2637698_length  SRR2637698_eff_length  SRR2637698_est_counts  SRR2637698_tpm

AAC73112              66                 8.78165                19                     267.332

AAC73113              2463               2364                   503                    26.2902

AAC73114              933                834                    453                    67.1128

AAC73115              1287               1188                   584                    60.7393

AAC73116              297                198                    12                     7.48841

AAC73117              777                678                    839                    152.899

AAC73118              1431               1332                   558                    51.7611

AAC73119              954                855                    9361                   1352.79

AAC73120              588                489                    1423                   359.558

transcripts id、raw read count、transcripts length、TPM補正値(参考)などが出力される。 

 

例2、手持ちのシーケンシングデータを解析する。例えばヒトゲノムのデータ。

docker run -v $(pwd):/dee2/mnt mziemann/tallyup hsapiens -f sample1_R1.fq.gz,sample2_R1.fq sample1_R2.fq.gz,sample2_R2.fq 

シロイヌナズナならmziemann/tallyup athalianaとする。

 

感想

SRAのRNA seqデータのクオリティを分析し、比較できる形で瞬時にカウントデータを取得できるデータベースです。モデル生物の研究者の方々にはたいへん実用的なデータベースだと思います。

引用

Digital expression explorer 2: a repository of uniformly processed RNA sequencing data
Mark Ziemann, Antony Kaspi, Assam El-Osta
GigaScience, Volume 8, Issue 4, April 2019

 

関連


 

*1

dataが自動でアップロードされるとあるので注意してください 

genome trackを可視化する svist4get

 

 次世代シークエンシングは、生命科学の複数のハイスループットな方法を生み出した。その多くは、既存のゲノムアセンブリへのショートリードのマッピングに基づいている。マッピングされたリードの密度および計算により得られたゲノムシグナルトラックの可視化は,バイオインフォマティクスシーケンシングデータ解析における最も日常的なタスクの1つである。1つのアプローチは、専用のゲノムブラウザの使用である。UCSC Genome Browser [ref.1] やZenbu [ref.2] のような最も一般的な汎用ツールはウェブベースであり,アップロードされたユーザデータと共に既存のゲノムアノテーションインタラクティブな探索を可能にする。しかし、場合によっては、ユーザデータをリモートサーバにアップロードするのは都合が悪く、統合ゲノムビューア [ref.3] や統合ゲノムブラウザ [ref.4] のようなグラフィカルユーザインタフェースを備えたスタンドアロンアプリケーションの助けを借りて、あるいはコンソール環境 [ref.5] で直接、データを可視化し、探索することができる。最後に、ハイブリッドアプローチがあり、例えば、BioUMLバイオインフォマティクスプラットフォーム [ref.6] は、ウェブベース版とスタンドアロン版の両方でゲノムブラウジング機能を提供する。 Webベースのゲノムブラウザは、処理された公開データの探索的なデータ分析には適しているが、カスタムのペーパー品質のグラフィック生成や、進行中または非公開の実験からのデータの探索には、スタンドアロンのツールの方が適している。さらに、ゲノムのいくつかの遺伝子座について複数の画像生成プログラム方法を有することはしばしば便利である。スクリプトコマンドラインインターフェイス(論文表1)を使用してこのタスクを解決することを目的としたツールがいくつか存在する。最も進んだツールの中には、Gviz [ref.7] とggbio [ref.8] がある。これらは、ゲノムトラックとアノテーションのペアーパー品質の図を作成するためのBioconductor Rパッケージである。コマンドラインユーティリティを好むユーザはfluff [ref.9] とngs plot [ref.10] を使うことができる。これらのツールはデータ解析のためのいくつかの高度な機能を提供するが、特定のゲノムにおけるゲノムトラックセグメントの可視化への最小限のアプローチしか可能にしない。 svist4getはPython3で実装されており、複数のpypiパッケージを使用している。

 SVist4getはLinux環境で開発されテストされている。(一部略)入力データとして、svist4getはゲノムシグナル用のbedGraph形式とゲノムアノテーションのgtf形式をサポートしている。出力データとして、svist4getはベクトル・グラフィックをpdfで生成し、ラスターグラフィックをpngでエクスポートできまる。ImageMagickはラスター(ポイント)出力を提供するために使用される。 特定のゲノムウィンドウとゲノムシグナルトラックのセットが与えられると、svist4getは、必要に応じて、画像幅とゲノムウィンドウの可視長を考慮して、シグナルトラックのmoving-average smoothingを自動的に実行する。しかし、svist4getは純粋な視覚化ツールであるため、技術データの変換やリードデプス正規化などの前処理は、deeptools[ref.13]、bedtools[ref.14]、UCSCユーティリティ[ref.15]などの外部ツールを使用して実行する必要がある。 標準的なシナリオやデータ探索でのsvist4getの適用を容易にするために、コマンド行インターフェースでは、トランスクリプトーム研究で発生するいくつかの実用的な使用例を、ユーザー側のスクリプト作成の手間をかけずにカバーしている。さらに、svist4getにはPython APIが用意されており、Pythonプログラム内から追加のカスタマイズやプログラムによる使用が可能である。(以下略)

 

quick start guide

https://bitbucket.org/artegorov/svist4get/src/9a3dbaec1bd19c487d76efb5d02deddaf04b7ef8/docs/QSGUIDE.md?fileviewer=file-view-default

 

 

インストール

Python 3.6.3環境でテストした(docker使用、ホストOS mac10.12)。

依存

本体 Bitbucket

python3 -m pip install svist4get

svist4get -h

# svist4get -h 

 

svist4get (version 1.2.10): a simple visualization tool for genomic tracks from sequencing experiments. 

 

COMMAND-LINE PARAMETERS

_____________________________

[POST-INSTALL EXAMPLE DATA]

 

--sampledata

Creates the 'svist4get' folder in the current working directory.

The folder will contain sample data sets and adjustable

configuration file templates.

_____________________________

[MANDATORY ARGUMENTS]

-t <id>

Visualization of the genomic window of a particular

transcript (using GTF transcript_id)

  OR

-g <id>

Visualization of the genomic window of a particular

gene (using GTF gene_id)

  OR

-t <id> -w tss|tis <upstream> <downstream>

Visualization of the genomic window centered on a transcription start

site|translation initiation site of a particular transcript

        (using GFF transcript_id) with a fixed offsets upstream and downstream

  OR 

-w <contig> <start> <end>

Visualization an arbitrary genomic window using genomic coordinates

 

-fa <file.fa>

Path to the genome assembly multifasta for the sequence track

NOTE: the multifasta file must contain all contigs in a single file.

 

-gtf <file.gtf>

Path to the gtf file with the transcript annotation

_____________________________

[OPTIONAL ARGUMENTS]

-bg <file1.bedGraph> [<file2.bedGraph>...]

Paths to bedGraph files for genomic signal tracks

 

-pbg <file1A.bedGraph file1A.bedGraph> [ -pbg <file2A.bedGraph file2B.bedGrap>...]

Paths to 2 bedGraph files to be used for the paired genomic signal track

        (e.g. to plot signals for + and - DNA strands together). Positive and negative range

of the Y-axis will be used for the first and second file, respectively

(or vice versa if the reverse complementary transformation using -rc parameter is requested)

The track will use absolute values from each bedGraph.

Multiple paired tracks can be generated by using multiple '-pbg' parameters.

For this track, the signs of values in each bedGraph will be ignored

 

-c <path> |A4_p1|A4_p2|A4_l 

Path to a configuration file or name of the premade config file

(can be A4_p1, A4_p2, A4_l)

 

-o <name> 

Output file name. Generates a pdf (vector) or png (raster) file based on

the provided extension. Both png and pdf will be generated if the extension is not specified.

 

-hi

Hide intronic segments (default: false)

 

-rc

Reverse-complement transformation of the genomic window (default: off)

 

-bgb mean|median|max|min|none

This parameter sets the function to aggregate bedGraph signal values for a single bar of the bedGraph tracks.

  Default: 'mean'. Alternatives: 'median', 'max', 'min'. 'none' forces plotting the raw signal at 1nt resolution.

 

-bg_log 

Logarithmic (log2) scale for the each bedGraph (Y-axis).

 

-nts

Show reference genomic sequence with single-nucleotide resolution

(default: hide)

 

-aas auto|tics|codons

Style of the aminoacid sequence track: tics (suitable for long regions e.g. > 160 nts,

  shows only start and stop codons with colored tics) or codons (show each aminoacid as labeled block)

(default: auto)

 

-it <text>

Image title

 

-xs <N>

X axis tics step (nt) (default: auto)

 

-lb <N1> [<N2>..]

Y axis tics levels for a bedGraph track, a space-separated list

of values, a separate value for each bedGraph track (default: auto)

 

-ys <N1> [<N2>...]

Step between axis tics a the bedgraph track Y-axis, a space-separated

list of values, a separate value for each bedGraph track (default: auto)

 

-bul <N1> <N2>... | max

  Upper limit for a bedGraph track Y-axis, can be:

1. a space-separated list of values: a separate axis limit for each

  bedGraph track;

2. max: a maximum reachable value across all visible tracks.

Default: reachable maximum for each track separately.

 

-bll <N1> <N2>... | min

  Lower limit for a bedGraph track Y-axis, can be:

1. a space-separated list of values: a separate axis limit for each

  bedGraph track;

2. min: a minimum reachable value across all visible tracks.

Paired bedGraph uses absolute values for both signal tracks,

since the negative Y-axis range is used to display the 2nd track.

 

-bl <tex1t> <text2> 

Text labels for the bedGraph tracks, space-separated list of values,

a separate value required for each bedGraph track

(default: names of the bedGraph files)

 

-bgc <color1>[color2>..]

Space-separated list of colors for bedGraph tracks. The colors should be specified

  by names listed in the palette file (default: internal palette or palettes/palette.txt).

  The list will be cycled through if shorter than the number of bedGraph tracks.

  Paired and basic bedGraph tracks use joint color list.

 

-gi <s-e> [<s1-e1>...]

Additional track showing selected genomic intervals as segments,

space-separated list of start-end pairs

 

-gil <tex1t> <text2>...

Labels for genomic intervals, space-separated list of strings,

a separate value required for each gemomic intervals track

 

-gic <color1>[<color2>..]

Space-separated list of colors for genomic intervals tracks, uses the same logic as -bgc parameter

 

-st <transcript_id1> <transcript_id2>...

A complete list of transcript IDs to show (others will be hidden,

applicable to -w and -g modes) 

 

-tl <transcript_id> <label>

Allows to replace the default <transcript_id> label with an arbitrary text <label>.

For re-labeling several transcripts, several '-tl' keys can be used.

 

-hrf <N> 

Highlight a particular reading frame (+0, +1 or +2) 

_____________________________

[MISCELLANEOUS ARGUMENTS]

-h, --help

Show this help message and exit

 

-v, --version

Show program version

 

-hf <x1> <x2> <label>

Highlight and label a given segment of the genomic window

(global coordinates on the selected contig)

 

-hc <NtNtNt>

  Highlight selected codon on the aminoacid sequence track, e.g. -hc TGA

 

dockerイメージも置いておきます。

docker pull kazumax/svist4get 
#カレントと/dataをシェアしてラン
docker run -itv $PWD:/data/ kazumax/svist4get
> source ~/.profile
> svist4get -h #help

 

テストラン

テストデータ生成

svist4get --sampledata

ls -alh svist4get_data/

# ls -alh svist4get_data/

total 22M

drwxr-xr-x 4 root root 4.0K Apr 29 13:12 .

drwxr-xr-x 1 root root 4.0K Apr 29 13:25 ..

-rw-r--r-- 1 root root 2.7K Apr 29 13:11 A4_l.cfg

-rw-r--r-- 1 root root 2.8K Apr 29 13:11 A4_p1.cfg

-rw-r--r-- 1 root root 2.7K Apr 29 13:11 A4_p2.cfg

-rw-r--r-- 1 root root 2.6K Apr 29 13:11 MATa.fasta

-rw-r--r-- 1 root root   24 Apr 29 13:11 MATa.fasta.fai

-rw-r--r-- 1 root root 1.4K Apr 29 13:11 MATa.gtf

-rw-r--r-- 1 root root 183K Apr 29 13:11 RiboCov_cut.bedGraph

-rw-r--r-- 1 root root 114K Apr 29 13:11 RiboProElong_cut.bedGraph

-rw-r--r-- 1 root root  12M Apr 29 13:11 S.cerevisiae.dna.fa

-rw-r--r-- 1 root root  415 Apr 29 13:11 S.cerevisiae.dna.fa.fai

-rw-r--r-- 1 root root 9.1M Apr 29 13:11 S.cerevisiae.gtf

-rw-r--r-- 1 root root 2.8K Apr 29 13:11 alt_A4_p1.cfg

-rw-r--r-- 1 root root  260 Apr 29 13:11 alt_palette.txt

drwxr-xr-x 2 root root 4.0K Apr 29 13:12 fonts

drwxr-xr-x 2 root root 4.0K Apr 29 13:12 help

-rw-r--r-- 1 root root 216K Apr 29 13:11 mRNACov_cut.bedGraph

-rw-r--r-- 1 root root 4.5K Apr 29 13:11 mata_ribo1_minus.bedGraph

-rw-r--r-- 1 root root 5.3K Apr 29 13:11 mata_ribo1_plus.bedGraph

-rw-r--r-- 1 root root 3.2K Apr 29 13:11 mata_ribo2_minus.bedGraph

-rw-r--r-- 1 root root 4.7K Apr 29 13:11 mata_ribo2_plus.bedGraph

-rw-r--r-- 1 root root  15K Apr 29 13:11 mata_rna1_minus.bedGraph

-rw-r--r-- 1 root root  16K Apr 29 13:11 mata_rna1_plus.bedGraph

-rw-r--r-- 1 root root  13K Apr 29 13:11 mata_rna2_minus.bedGraph

-rw-r--r-- 1 root root  13K Apr 29 13:11 mata_rna2_plus.bedGraph

-rw-r--r-- 1 root root  217 Apr 29 13:11 palette.txt

-rw-r--r-- 1 root root  404 Apr 29 13:11 triplet_code.txt

 

sampleディレクトリができたらsvist4getを実行する。ゲノムのfasta、gtf、シーケンシングデータのbedgraphファイルを指定する。

視覚化する遺伝子はYFL031Wとする。このYFL031Wはgtfのtranscript IDから指定している。

svist4get \
-gtf svist4get_data/S.cerevisiae.gtf \
-fa svist4get_data/S.cerevisiae.dna.fa \
-bg svist4get_data/RiboCov_cut.bedGraph \
-t YFL031W
  • -gtf   Path to the gtf file with the transcript annotation
  • -fa    Path to the genome assembly multifasta for the sequence track
    NOTE: the multifasta file must contain all contigs in a single file.
  • -bg   Paths to bedGraph files for genomic signal tracks

     [MANDATORY ARGUMENTS]

  • -t <id>    Visualization of the genomic window of a particular transcript (using GTF transcript_id)   OR
  • -g <id>     Visualization of the genomic window of a particular gene (using GTF gene_id)     OR
  • -t <id> -w tss|tis <upstream> <downstream>    Visualization of the genomic window centered on a transcription start site|translation initiation site of a particular transcript (using GFF transcript_id) with a fixed offsets upstream and downstream  OR
  • -w <contig> <start> <end>   Visualization an arbitrary genomic window using genomic coordinates

出力。赤い波部分が指定したデータ; RiboCov_cut.bedGraph。

f:id:kazumaxneo:20190429134946j:plain


データのbedgraphは複数指定できる。下では3つ指定している。イントロンは非表示とする(-hi)。

svist4get \
-gtf svist4get_data/S.cerevisiae.gtf \
-fa svist4get_data/S.cerevisiae.dna.fa \
-bg svist4get_data/RiboProElong_cut.bedGraph svist4get_data/RiboCov_cut.bedGraph svist4get_data/mRNACov_cut.bedGraph \
-t YFL031W \
-it 'Gene HAC1' \
-bl 'A-Site Ribo-Seq' 'Coverage Ribo-Seq' 'RNA-Seq' \
-hi \
-hrf 0 \
-c A4_p1
  • -bl <tex1> <text2>    Text labels for the bedGraph tracks, space-separated list of values, a separate value required for each bedGraph track (default: names of the bedGraph files).
  • -hi     Hide intronic segments (default: false).
  • -hrf    Highlight a particular reading frame (+0, +1 or +2).
  • -it      Image title.
  • -c      Path to a configuration file or name of the premade config file (can be A4_p1, A4_p2, A4_l).

3つのtarackが視覚化された。タイトルは-it 'Gene HAC1'で指定し、それぞれのtrack名は-bl 'A-Site Ribo-Seq' 'Coverage Ribo-Seq' 'RNA-Seq'で指定している。

f:id:kazumaxneo:20190429135728j:plain


領域を指定して転写開始点付近を視覚化する(translation initiation site (tis)の-100から+30まで)。

svist4get \
-gtf svist4get_data/S.cerevisiae.gtf \
-fa svist4get_data/S.cerevisiae.dna.fa \
-bg svist4get_data/RiboProElong_cut.bedGraph svist4get_data/RiboCov_cut.bedGraph svist4get_data/mRNACov_cut.bedGraph \
-bl 'A-Site Ribo-Seq' 'Coverage Ribo-Seq' 'RNA-Seq' \
-t YOR030W \
-w tis 100 30 \
-it 'Upstream ORF of DFG16 gene' \
-gi 386730-386772 \
-gi 386824-386999 \
-gil 'Upstream ORF' 'CDS' \
-hf 386730 386772 ' '
  • -w <contig> <start> <end>   Visualization an arbitrary genomic window using genomic coordinates.
  • -gi <start-end>   Additional track showing selected genomic intervals as segments, space-separated list of start-end pairs.
  • -gil <tex1t> <text2>...   Labels for genomic intervals, space-separated list of strings, a separate value required for each gemomic intervals track.
  • -hf <x1> <x2> <label>
    Highlight and label a given segment of the genomic window (global coordinates on the selected contig).

f:id:kazumaxneo:20190429135808j:plain

 

出力名はHAC1とする。

svist4get \
-gtf svist4get_data/S.cerevisiae.gtf \
-fa svist4get_data/S.cerevisiae.dna.fa \
-bg svist4get_data/RiboProElong_cut.bedGraph svist4get_data/RiboCov_cut.bedGraph svist4get_data/mRNACov_cut.bedGraph \
-bl 'A-Site Ribo-Seq' 'Coverage Ribo-Seq' 'RNA-Seq' \
-g YFL031W \
-it 'Gene HAC1' \
-hf 75839 76091 intron \
-c svist4get_data/A4_p1.cfg \
-o HAC1
  • -g    Visualization of the genomic window of a particular gene (using GTF gene_id)

出力のHAC1.pdf

f:id:kazumaxneo:20190429135932j:plain

 

svist4get \
-gtf svist4get_data/S.cerevisiae.gtf \
-fa svist4get_data/S.cerevisiae.dna.fa \
-bg svist4get_data/RiboProElong_cut.bedGraph svist4get_data/RiboCov_cut.bedGraph svist4get_data/mRNACov_cut.bedGraph \
-bl 'A-Site Ribo-Seq' 'Coverage Ribo-Seq' 'RNA-Seq' \
-w I 108842 113506 \
-it 'CCR4 and FUN26 genes' \
-rc \
-c svist4get_data/A4_p1.cfg
  • -rc   Reverse-complement transformation of the genomic window (default: off)

f:id:kazumaxneo:20190429140040j:plain

 

 

svist4get \
-gtf svist4get_data/MATa.gtf \
-fa svist4get_data/MATa.fasta \
-pbg svist4get_data/mata_ribo1_plus.bedGraph svist4get_data/mata_ribo1_minus.bedGraph \
-pbg svist4get_data/mata_ribo2_plus.bedGraph svist4get_data/mata_ribo2_minus.bedGraph \
-pbg svist4get_data/mata_rna1_plus.bedGraph svist4get_data/mata_rna1_minus.bedGraph \
-pbg svist4get_data/mata_rna2_plus.bedGraph svist4get_data/mata_rna2_minus.bedGraph \
-w V01313.1 390 2430 \
-rc -it 'MAT locus of MATa yeast strain' \
-bl 'Ribo-Seq (sample1)' 'Ribo-Seq(sample2)' 'RNA-Seq (sample1)' 'RNA-Seq (sample2)' \
-bgc brightorange brightorange purple purple \
-c svist4get_data/alt_A4_p1.cfg \
-hf 1224 1305 ' ' \
-hf 1533 1638 'Translated open reading frames' \
-hf 1692 1932 ' '

f:id:kazumaxneo:20190429140132j:plain

 

bamからbedgraphの作成はbedtools genomecovなどで行える(*1)。

bedtools genomecovを"-bg"付きで実行する。

bedtools genomecov -ibam input.bam -bg > output.bg
  • -bg     Report depth in BedGraph format. For details, see: genome.ucsc.edu/goldenPath/help/bedgraph.html
  • -ibam    The input file is in BAM format. Note: BAM _must_ be sorted by position.

 

引用

svist4get: a simple visualization tool for genomic tracks from sequencing experiments

Artyom A. Egorov, Ekaterina A. Sakharova, Aleksandra S. Anisimova, Sergey E. Dmitriev, Vadim N. Gladyshev, Ivan V. Kulakovskiy
BMC Bioinformatics 2019 20:113

 

関連

 

*1

from bam to bedgraph

bam file to .bedgraph

 

(病原性の)大量のバクテリアゲノムの自動解析パイプライン TORMES

2019 12/20 インストール手順修正

2019 12/21, 12/22結果追記

 

連休中は不定期更新になります。よろしくお願いいたします。 

 

 ハイスループットシーケンシング(HTS)技術の進歩およびシーケンシングコストの削減は、全ゲノムシーケンシング(WGS)が多くの伝統的な実験室アッセイおよび手順に取って代わることができるようなものである。 HTSプラットフォームによって生成された大量のデータを活用するには、相当なコンピューティングスキルが必要である。これが、日常的なラボ手法としてのWGSの実装における主なボトルネックである。ゲノムシークエンシングの専門知識がなくても膨大な量の結果を研究者や臨床医に提示する方法も重要な問題である。
 この論文では、一連のバクテリアからIlluminaプラットフォームによって生成されたHTSデータのWGS分析を実行するための、オープンソースでユーザーフレンドリーなコマンドラインパイプラインであるTORMESを紹する。 TORMESは、バイオインフォマティクスバックグラウンドのではない研究者向けに設計されている。シーケンスクオリティフィルタリング、de novo assembly、レファレンスに対するドラフトゲノムの整列、ゲノムアノテーション、MLST、抗生物質耐性検索および病原性検索、遺伝子の比較やパンゲノムの比較などのバイオインフォマティクス分析のステップを自動化する。これは、インターネット接続を必要とせずに、非常に単純な指示に従ってrawシーケンシングデータから直接行うことができる。
 ゲノムシークエンシングの専門知識がなくても研究者や臨床医が理解できる形式で大量のデータをまとめることが、バクテリアゲノミクスの大きな問題として認識されている(Köseret al、2012)。 TORMESは処理中に生成されたすべてのファイルを保存し、分析が完了すると、結果は対話的なWebライクなレポートにまとめられ、修正、共有、そして人間工学的な比較が可能になる。
 レポートは、分析ごとに固有の、自動的に生成されたRMarkdownコードファイルを使用してR環境で生成される。また、より専門的なユーザーが分析を深め、ユーザー固有のレポートのコードを変更できるように、別のフォルダーにも保存されている。
TORMESは無制限数のサンプルに使用でき、さまざまなソース(臨床、糞便、動物および食品関連)からの多数の種(EscherichiaSalmonellaClostridium、およびKlebsiella を含む)の分離株における数百のバクテリアゲノムで試験され、概算として、これら数百のバクテリアゲノムの約50倍のシーケンシングデプスのTORMES分析には、124 GBのRAM と32コアのコンピューターで16時間かかった。

 

f:id:kazumaxneo:20190412103508j:plain

 

TORMESに関するツイート


インストール

 ubuntu18.04でテストした。

依存

Additional software when working with -g/--genera Escherichia.

Additional software when working with -g/--genera Salmonella.

本体 Github

#condaの仮想環境を作って導入するように設計されている。
wget https://anaconda.org/nmquijada/tormes-1.0/2019.04.25.180147/download/tormes-1.0.yml
conda env create -n tormes-1.0 --file tormes-1.0.yml
conda activate tormes-1.0

> tormes -h

# tormes -h

 

This is TORMES  version 1.0

Developed by Narciso M. Quijada <https://github.com/nmquijada/tormes>

 

usage: /root/.pyenv/versions/miniconda3-4.3.21/envs/tormes-1.0/bin/tormes <options>

 

OBLIGATORY OPTIONS:

-m/--metadata Path to the file with the metadata regarding the samples

The file must have an specific organization for the program to work

If you don't have any or you would like to have an example or extra information, 

please type: 

/root/.pyenv/versions/miniconda3-4.3.21/envs/tormes-1.0/bin/tormes example-metadata

-o/--output Path and name of the output directory

 

OTHER OPTIONS:

-a/--adapter Path to the adapters file

(default="/root/.pyenv/versions/miniconda3-4.3.21/envs/tormes-1.0/bin/../files/adapters.fa")

--assembler Select the assembler to use. Options available: 'spades', 'megahit'

(default='spades')

        -c/--config     Path to the configuration file with the location of all dependencies

                        (default="/root/.pyenv/versions/miniconda3-4.3.21/envs/tormes-1.0/bin/../files/config_file.txt")

--citation Show citation

--fast Faster analysis (default='0')

('megahit' is used as assembler and contig ordering and pangenome analysis are disabled)

--filtering Select the software for filtering the reads.

Options available: 'prinseq', 'sickle', 'trimmomatic'

(default="prinseq")

-g/--genera Type genera name to allow special analysis (default='none')

Options available: 'Escherichia', 'Salmonella'

-h/--help Show this help

--min_len Minimum length to the reads to survive after filtering (default=125) <integer>

--no_mlst Disable MLST analysis (default='0')

--no_pangenome Disable pangenome analysis (default='0')

-q/--quality Minimum mean phred score of the reads to survive after filtering (default=25) <integer>

-r/--reference Type path to reference genome (fasta, gbk) (default='none')

Reference will be used for contig ordering of the draft genome

-t/--threads Number of threads to use (default=1) <integer>

--title Path to a file containing the title in the project that will be used as title in the report

Avoid using special characters. TORMES will perform a default title if this option is not used

-v/--version Show version

 

 

For further explanation please visit: https://github.com/nmquijada/tormes

 

下で実行時にmauve.jarのパスが間違っていたので、エラーメッセージ通りmauveのディレクトリ名を修正した。

 

 

データベースの準備

tormes-setup

 

Setting up config_file.txt

 

TORMES is installed and ready to use. Enjoy!

完了。

minikrakenのデータベースなどもダウンロードされる。10GB近くある。また各ツールのパスを記載したconfigファイル(リンク)も自動作成される。

 

実行方法

1、metadata準備

ランにはfastqのパスや名前を記載したテキストファイル(メタデータ)を指定する。exampleメタデータを見てみる。

#exampleメタデータテキスト生成
tormes example-metadata

カレントにタブ区切りテキストsamples_metadata.txtができる。

以下のように1行目はコメント、2行目以降にサンプル情報を記載していく。1列目がサンプル名、2列目ペアエンドfastqのR1のフルパス、3列目にペアエンドfastqのR2のフルパス、4列目以降にコメントを記載する。コメント列は複数あっても良い。

f:id:kazumaxneo:20190412102548j:plain

上では6サンプル指定している。名前はsample1、sample2...にしている。

 

このような感じで用意する。コピペしてfastqのパス、Descriptionカラムと末尾のカラムを修正するだけで使えるはずです。

Samples Read1 Read2 Description Use_as_many_descrpition_colums_as_wanted
sample1 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 2018 AMR1
sample2 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 2012 AMR2
sample3 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 2006 AMR3
sample4 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 2000 AMR4
sample5 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 1994 AMR5
sample6 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 2019 AMR6
sample7 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 2019 AMR7
sample8 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 2019 AMR8
sample9 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 2019 AMR9
sample10 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 2019 AMR10
sample11 //tormes/SRRXXXXXXX_1.fq //tormes/SRRXXXXXXX_2.fq S.enterica isolated in 2019 AMR11

 

 

2、run

このメタデータを指定してランする。

tormes --metadata input_metadata.txt --output outdir --reference ref_genome.fasta --threads 32 --genera Salmonella

 EscherichiaSalmonellaのみ、追加分析が行える。上ではSalmonellaを指定している。たくさんのプロセスがあるため、threadは利用できるだけ指定する。

 

ジョブが終わると、出力ディレクトリに分析ごとのサブディレクトリができる。

f:id:kazumaxneo:20190412101240j:plain

 順調にランしているように見えたが、report.htmlの出力に失敗した。修正できたら追記します。

 

追記

ランできるようになった。

サマリーレポート

f:id:kazumaxneo:20191221110535p:plain

f:id:kazumaxneo:20191221110538p:plain

f:id:kazumaxneo:20191222160642p:plain

f:id:kazumaxneo:20191222160910p:plain

f:id:kazumaxneo:20191221110552p:plain

f:id:kazumaxneo:20191222160952p:plain

f:id:kazumaxneo:20191221110610p:plain

f:id:kazumaxneo:20191222160525p:plain

f:id:kazumaxneo:20191222160531p:plain

f:id:kazumaxneo:20191222160536p:plain

f:id:kazumaxneo:20191221110640p:plain

f:id:kazumaxneo:20191222160605p:plain

f:id:kazumaxneo:20191222161051p:plain

f:id:kazumaxneo:20191221113625p:plain

f:id:kazumaxneo:20191222160629p:plain


引用

TORMES: an automated pipeline for whole bacterial genome analysis

Quijada NM, Rodríguez-Lázaro D, Hernández M

Bioinformatics. 2019 Apr 8

 

関連


 

 

 

Bashの履歴検索(reverse search)を改善する McFly

 

McFlyはデフォルトのctrl-r によるBash履歴検索をインテリジェントなニューラルネットワークを使った検索エンジンで置き換える。エンジンは作業ディレクトリと最近実行されたコマンドのコンテキストを考慮に入れて優先順位を変更する。

 

インストール

homebrewを使ってmacos10.14に導入した。

依存

Github

 

brewで導入できる。

brew tap cantino/mcfly https://github.com/cantino/mcfly
brew install mcfly

#以下を~/.bashrcか~/.bash_profileに記載し、コンソールを再起動する(またはsource ~/.bashrc)。
if -r "$(brew --prefix)/opt/mcfly/mcfly.bash" ; then
source "$(brew --prefix)/opt/mcfly/mcfly.bash"
fi

iTerm2を使っている場合、optionの1つを外す必要がある。Githubを参照。

 

実行方法

Ctrl + rを押す。画面が一時的にclearされるので、コマンドを途中まで打つ。候補が表示されるので、十字キーで選択、Enterキーで実行する。

f:id:kazumaxneo:20190427001526p:plain

キャンセルするにはもう一度Ctrl + rを押す。

 

引用

GitHub - cantino/mcfly: Fly through your shell history. Great Scott!

 

他のTips


参考にしたページ