macでインフォマティクス

macでインフォマティクス

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

スタンドアロンのGUI環境で動作する抗生物質耐性遺伝子検索ツール SSTAR

 

 Antimicrobial resistance(AR)は、近年急速に増加する現象であり、世界的な公衆衛生上の重大な脅威である(論文より ref.1,2)。 EnterobacteriaceaeやPseudomonasのような新興のカルバペネム抗生物質耐性バクテリアは、多くの治療法を無効にしているため特に懸念される。カルバペネム耐性の出現は、カルバペネムおよび他のβ-ラクタム系抗生物質を不活性化するタンパク質をコードする獲得性のカルバペネマーゼ遺伝子によって主に媒介される(ref.3)。拡張型β-ラクタマーゼ(ESBL)またはAmpC型β-ラクタマーゼの産生と組み合わせることができる、外膜ポーリンおよび排出メカニズムの改変または喪失を含む他のメカニズムもまた記載されている(ref.4)。 AR遺伝子は、水平遺伝子導入によって獲得でき、またバクテリアクロモソーム遺伝子座におけるspontaneousな突然変異のために生じ得る。新規の対立遺伝子の予測を含むAR決定基の迅速な同定は、AR病原体の分子疫学を理解し、感染制御介入を知らせるために重要であり、これは患者治療の指針として有用であり得る。

 近年、DNAシーケンシングのコストは急速に減少しており、この技術は世界中の大規模なコミュニティに利用可能になっている(ref.5)。アウトブレイク調査、疫学サーベイランス、およびARに関連する遺伝子の検出(ref.6-8)で、全ゲノムシーケンス(WGS)データが株の特徴付けに使用されている。 AR遺伝子を含む多くのオンラインデータベースは、現在一般に公開されている。最も完全で最新の2つのリポジトリは(論文執筆時点)、ResFinder(version3 database)(ref.9)(紹介)とComprehensive Antibiotic Resistance Database(CARD: CARD HP)(ref.10)である。これらのサービスでは、ユーザーはバクテリアゲノム配列をWebサービスにアップロードし、その後取得した抗菌抵抗性遺伝子のリストを取得できる。しかしこれらのWebサービスは比較的遅く、複数のゲノムファイルを照会したいユーザーのハードルになる可能性がある。

 SRST2(ref.11)、ARG-ANNOT(ref.12)、ABRicate (https://github.com/tseemann/ abricate)(紹介)などのオフラインソリューションも存在する。獲得された耐性遺伝子検索に加えて、ARG-ANNOTデータベースは、ARに関連することが知られているクロモソームターゲット遺伝子における点突然変異解析を含む。 SRST2、ARG-ANNOT、およびABRicateはローカルで実行され、ユーザーにWebベースのサービスよりも多くの制御を与える。 SRST2とABRicateはコマンドラインインターフェイスのみを提供し、ARG-ANNOTは古いソフトウェアパッケージ(ref.13)で動作する。 ARG-ANNOTおよびABRicateは、推定上の新しいAR遺伝子を検出する方法を提供しているが、上記のツールはいずれも簡単な検出方法を提供していない:既存のAR遺伝子、または耐性に関連する他の新規遺伝子またはバリアント。

 ここでは、修正されたARG-ANNOTデータベース、スタンドアロンのBLAST(ref.14)、および既知のAR遺伝子の検出を可能にする使いやすいグラフィカルユーザーインターフェイスを組み合わせたツールSSTAR(抗菌抵抗性のための配列検索ツール)推定された新しい変異体および切断された遺伝子を、迅速で、ローカルで、容易に更新可能なツールで予測する。

 

blastn解析をcontig配列と抗生物質耐性遺伝子間で実行し、耐性遺伝子を探すシンプルなツールです。

 

インストール

mac os10.13のjava1.8環境でテストした。

依存

本体 GIthub

本体をダウンロードして解凍する。

環境変数に以下のように定義しておく。(blastにパスが通ってなければフルパスで)。

export BLASTN="blastn"
export MAKEBLASTDB="makeblastdb"

 

ラン

1、contig配列の準備。

解析にはNGSデータをアセンブリして得られるcontig配列が必要になる。例えばspadesでd novoアセンブリを実行する(spades紹介)。オンラインで実行するなら、CyVerse(紹介)やGalaxy(training)を使う。

brew install spades
spades.py -1 pair1.fq.gz -2 pair2.fq.gz -k auto --careful -t 8 -o output

output/scaffolds.fastaができる。

 

2、実行

jarファイルを実行することでGUIパネルが出現する。

cd Sequence-Search-Tool-for-Antimicrobial-Resistance-SSTAR--master/
java -jar SSTAR.jar # *1

*1 windowsユーザーはSSTAR_windows.jarを実行する。

windowが出現する。

f:id:kazumaxneo:20180703102748p:plain

 

一番上のwindowにcontigを入力し、その下のwindowに検索した抗生物質耐性遺伝子のFASTAファイルを入力する。ダウンロードしたディレクトリに、ResfinderとARG-ANNOTから構築されたFASTAファイルがあるので、ここではそれを使う(Latest_AR_database/ResGANNOT_srst2.fasta)。3番目のウィンドウに検出下限閾値を記載し(0-1)、最後にIdentify resistance genes!をクリックして実行する。

 

blastn解析が終わると、検出された抗生物質耐性遺伝子と、耐性遺伝子が見つかったcontig名が表示される。ここではカナマイシン耐性カートリッジが入ったバクテリアのcontig配列を調べた。Aph3遺伝子がヒットしている。

f:id:kazumaxneo:20180703124123j:plain

 

引用

SSTAR, a Stand-Alone Easy-To-Use Antimicrobial Resistance Gene Predictor

Tom J. B. de Man, Brandi M. Limbago

mSphere. 2016 Jan-Feb; 1(1): e00050-15.

点変異も考慮して抗生物質耐性遺伝子を検出する PointFinder

 

 バクテリア間の遺伝子の水平伝播は、しばしば後天性抗菌剤耐性のメディエーターになると考えられている。しかし、抵抗性を付与するもう一つの重要な方法は突然変異抵抗性である。WGSは、水平伝播で獲得された耐性検出のための細菌分離株の表現型感受性試験に取って代わりうるものであることが以前に示されている(論文より ref.1.2)。バクテリアのクロモソーム突然変異へのマッピングのためのデータベースも結核菌に対して開発されている(ref.3)。足りないのは、ほとんどの細菌種で起こるクロモソームの点突然変異に起因した抗生物質耐性を容易に調べ解釈可能な結果を​​提供する適切なバイオインフォマティクス法である。

 この研究では、バクテリアのWGSデータにおいて、抗生物質抵抗性に関連するバクテリアのクロモソーム上の点突然変異検出のための新規ウェブツール、PointFinderが開発された。 PointFinderは、WGSデータを使い水平伝播で獲得された耐性遺伝子を検出する既存のWebサーバーツールResFinder、の拡張機能である。性能は、点突然変異を検出するためのマッピング法で比較され、両方の結果を表現型感受性試験の結果と比較して、これらの方法を標準的な表現型試験の代替として使用できる可能性を検証した。

 PointFinderは2つのデータベースで構成されている:fasta形式の全クロモソーム遺伝子データベース。およびコドン位置および置換に関する情報を含むクロモソーム突然変異データベース。 PointFinderは、BLASTnを使用して、クロモソーム遺伝子データベース内の各遺伝子のベストマッチを特定し、80%以上の同一性のヒットのみをさらに分析する。プログラムは、クエリの各位置(入力シーケンスで見つかったシーケンス)と、対象(データベースシーケンス)の対応する位置を比較する各アライメントを実行する。全てのミスマッチは保存され、クロモソーム突然変異データベースと比較される。ユーザーは、すべてのミスマッチを見たいのか、クロモソームデータベースのポジションにある既知のミスマッチだけを見たいのかを選択することができる。

 

 

インストール

依存

  • Biopython
  • Blast

Bitbucket

https://bitbucket.org/genomicepidemiology/pointfinder

git clone --recursive https://bitbucket.org/genomicepidemiology/pointfinder.git
cd pointfinder
./UPDATE_DB database
./VALIDATE_DB database

  

ラン

python pointfinder.py -i test.fsa -o OUTFOLDER -s e.coli 

 

 

webでは、ResFinderの機能の一部として利用可能。ResFinderは現在version3.0となっている。

https://cge.cbs.dtu.dk/services/ResFinder-3.0/

f:id:kazumaxneo:20180702210928p:plain

 アセンブリして作ったcontig配列かシーケンスデータなのかを指定する。ファイルを選択して、Uploadボタンをクリックしアップロードする。Upload完了後、自動で画面が切り替わるが、混雑しているときはそのまま停止してしまう場合がある。何度か繰り返せばうまく行く可能性もあるが、サーバーに負荷がかかっているのでしばらくしてからやり直すのが無難。

 

点変異を許容して検索するならChromosomal point mutationsにチェックを入れる。

f:id:kazumaxneo:20180702213138p:plain

 

f:id:kazumaxneo:20180702213140p:plain

 

出力例

CGE Server

f:id:kazumaxneo:20180702213307p:plain

ヒットした 抗生物質耐性遺伝子が表示される。

 

Chromosomal point mutationsにチェックが入っていれば、点変異により 抗生物質耐性遺伝子がコードするタンパク質のアミノ酸が変化するかどうかも調べられる。

f:id:kazumaxneo:20180703093308p:plain

f:id:kazumaxneo:20180703093221p:plain

 

検出される抗生物質耐性遺伝子一覧

https://cge.cbs.dtu.dk/services/ResFinder/database.php

 

引用

PointFinder: a novel web tool for WGS-based detection of antimicrobial resistance associated with chromosomal point mutations in bacterial pathogens
Ea Zankari Rosa Allesøe Katrine G Joensen Lina M Cavaco Ole Lund Frank M Aarestrup

Journal of Antimicrobial Chemotherapy, Volume 72, Issue 10, 1 October 2017, Pages 2764–2768

抗生物質耐性遺伝子を検出する KmerResistance

 

 抗生物質は、ヒトおよび家畜の両方で世界中で広く使用されており、疾患の治療または急速な成長を保証している。長年にわたり、これは抗生物質耐性菌の出現、選抜および普及のための好ましい条件を作り出してきた(ref.1)。

 バクテリアの耐性プロファイルを迅速かつ確実に決定することはサーベイランス(wiki)にとって重要であり、臨床的処置を導くためにも重要である。近年、次世代シーケンシング(NGS)技術は安価で迅速かつ正確になり、サーベイランスや迅速な臨床診断に日常的に利用されている(ref.2)。一つの懸念は、NGSが表現型感受性を予測する能力である。研究により、予測された感受性と測定された感受性が非常に高い一致を示すことが明らかにされている(ref.2)。

 WGSデータの遺伝子を同定するために、いくつかの異なる方法が開発されているが、どの方法が最適であるかについての合意はない。これらの方法は、参照データベースと比較する前にrawシーケンスリードからのコンティグのアセンブリすること、またはシーケンスリードをリファレンス配列に直接マッピングすること、の大きく2つのグループに分けることができる。Zankariら3)は、WGSデータから抗菌剤耐性を検出する第1のアプローチの1つを開発した。この方法は、WGSをアセンブリし、BLASTを用いて耐性遺伝子を同定することに基づいている。このアプローチの1つのリスクは、2つ以上のコンティグに分割された場合、遺伝子の同定が見落とされる可能性があることである。これは、データとアセンブリのクオリティが低い場合に発生する可能性がある。

 Inouyeら4)は、Bowtie2(ref.5)を用いてraw WGSデータを抵抗性遺伝子に直接マップして感度を上げ、それによってパフォーマンスを向上させるアプローチを提案した(ref.4 pubmed)。しかし、raw WGSデータを直接マッピングするアプローチでは、汚染物質などのWGSデータのノイズによる誤検出増加などの問題が発生する。

 感度を維持し、偽陽性率を低く抑えるために、ここでは新しいアプローチを提示する。 k-mers(長さkのDNA配列の断片)を使用して、raw WGSデータをリファレンスデータベースに対してマッピングし、耐性遺伝子を同定するだけでなく種を決定する。次いで、種のリファレンスゲノムに対しマッピングして、抗菌剤耐性予測を正規化する。

この新しいアプローチ(KmerResistance)と、SRST2およびResFinderを、ヒト(n = 143)およびブタ(n = 196)の合計339の細菌分離株の標準的な表現型感受性試験で比較した(これにはEscherichia coli、Salmonella Typhimurium、Enterococcus faecalisおよびEnterococcus faeciumおよび27種の抗生物質が含まれる)。

3つのツールについて

  • ResFinderは、シーケンスデータをアセンブリした配列と耐性遺伝子とのblast検索により耐性遺伝子を同定する。
  • SRST2はBowtie2を使いシーケンスデータをユーザー指定のデータベース(例えば耐性遺伝子)に対して直接マッピングして耐性遺伝子を同定する。いくつかの配列は高い類似性を共有するため、SRST2は、CD-hitを80%の同一性閾値で使い、クラスタ化する。 これにより、各クラスター/遺伝子の最もよく一致する対立遺伝子のみが報告されることが保証される。SRST2の方法は、同定された遺伝子のSNPs解析などさらなる分析も可能なものである。ResFinderと比較してSRST2はより感受性が高く、シーケンスエラーやシーケンスキャリーオーバーなどのコンタミネーションの可能性を考慮する必要がある。そのため、論文中では、SRST2を最小デプス5の閾値設定下でテストしている。
  • KmerResistanceは、raw WGSデータからバクテリアタイピングを実行するKmerFinder(10,11)上に構築された。 KmerFinderおよびKmerResistanceは、クエリーゲノムと耐性遺伝子データベースとの間に共存するk-merの数を調べる。正確な一致のみが報告され、これは、カバレッジおよびアイデンティティが同じ(つまりカバレッジ)として報告されることを意味する。データベース内の遺伝子間で同一のk-mersに起因する複数のヒットを避けるために、各k-merは、最初に、k-mer適合数が最も高い遺伝子にのみアサインされる。この後、最良のヒットにマッピングされたk-mersが取り除かれ、残りのリードで手順が繰り返される。その性質から「winner takes all strategy」と命名した。抗生物質耐性遺伝子の同定に加えて、KmerFinderと同じスキームに従って種の予測を行い、その結果を出力することもできる。

 

ラン

使い方

https://cge.cbs.dtu.dk/services/KmerResistance/instructions.php

webサーバー

シングルエンド、ペアエンドのfastqをアップロードする。

f:id:kazumaxneo:20180702201222p:plain

出力について

https://cge.cbs.dtu.dk/services/KmerResistance/output.php

 

O-157のNGSデータでテストした時の結果。サブサンプリングしてx10ほどにカバレッジを減らしているがうまく検出できている。

f:id:kazumaxneo:20180703094433p:plain

species results、抗生物質耐性遺伝子検索結果などがダウンロードできる。 

 

 

Center for Genomic Epidemiologyのサーバーはたくさんのツールがあり、常にジョブが走っているのかBusyなことが多いです。余裕を持って使ってください。

 

引用

Benchmarking of methods for identification of antimicrobial resistance genes in bacterial whole genome data

Clausen PT, Zankari E, Aarestrup FM, Lund O.

J Antimicrob Chemother. 2016 Sep;71(9):2484-8.

腸内細菌科(エンテロバクター科)のプラスミド同定ツール PlasmidFinder

 

 プラスミドは、自律複製が可能であり、異なるバクテリア種とクローンとの間で移動可能な二本鎖の環状または線状DNA分子である。既知のプラスミドのほとんどは、抗生物質耐性または病原性遺伝子のようなバクテリア宿主上で陽性選択される表現型を付与するため、同定されてきた。このような特徴は、異なる起源および地理的起源のバクテリア間で異なるプラスミド型の普及を助けている。抗微生物耐性遺伝子または病原性遺伝子を保有するプラスミドを獲得することによって、毒性または多剤耐性細菌クローンの蔓延が大幅に変わる可能性がある。したがって、異なるバクテリアクローンの分子疫学を研究するだけでなく、転移可能なプラスミドの分子疫学を研究し、理解することも重要である。この特定の目的のためには、プラスミドタイピングシステムが必要である。

 ほとんどのプラスミドは、複製をアクティブにし制御することができる機能をコードするレプリコンと呼ばれる特定の領域を含む(論文より ref.1)。 2005年以来、腸内細菌科のメンバー(ref.2)に存在する主要なプラスミドファミリーのレプリコンを、multiplex PCRでターゲットにするPCRベースのレプリコンタイピング(PBRT)スキームが利用可能である。この方法は当初、腸内細菌科(Enterobacteriaceae)の18の主要非適合性(Inc)群に属するプラスミドのレプリコンを検出するために開発された(ref.3)。より最近では、新規レプリコンおよびプラスミド型が全ゲノムおよびプラスミドハイスループットシーケンスによって同定され、PBRTを25の異なるレプリコン(ref.4、-9)の同定にまで拡張した。しかしながら、この方法はmultiplex PCRに基づいており、これはさらに多くの群をカバーするためには困難であり、特にこの変異がプライマー結合部位内にある場合にはレプリコン変異体の検出には必ずしも適していない。 (すなわち、耐性遺伝子含量、multilocus sequence typing [MLST]、系統樹、血清型などによって決定される配列タイプ[ST])、plasmid typing(PBRTまたは縮重プライマーMOBタイピング[ 10]による)は、疫学調査中に無関係な関連株の比較分析に使用されている。

 

(2段落省略)

 ここでは、腸内細菌科のプラスミドの迅速な同定に有用な2つの新しい使いやすいWebツールの設計について説明する。このツールは、プラスミドに関連する抗菌剤耐性の広がりの疫学的および臨床的微生物学的研究にとって重要である。現在、GenBankで利用可能なプラスミド配列の大部分は自動的にアノテーションされており、多くの場合、レプリコン配列の注釈はプラスミド分類スキーム(Incグループまたはrelaxaseグループ)を参照していない。従って、BLASTnによるNCBIへの直接的なsubmitは、研究中のプラスミドの系統を認識するために容易に使用することができない。    PlasmidFinder Webツールは、微生物学者に管理された、全ゲノムシーケンス(WGS)から腸内細菌科由来のプラスミドを同定するキュレートされたプラスミドレプリコンのデータベースであり、バイオインフォマティクススキルを要求するハイスループットのrawシーケンスリードアセンブリコンティグ、アセンブリされたサンガー配列の使用なしに利用することができる。PlasmidFinderは、WGS内のレプリコンの検出だけでなく、研究中のプラスミドをIncグループの既存情報系統にアサインし、可能な基準プラスミドを示唆する。pMLST Webツールは、現在利用可能なpMLSTスキームを持つ5つの非互換性グループについて、同じ種類のデータに対してpMLST解析を実行できる。

 

Bitbucket

https://bitbucket.org/genomicepidemiology/plasmidfinder

 

BinPacker に関するツイート

 

マニュアル

https://cge.cbs.dtu.dk/services/PlasmidFinder/instructions.php

 

webサーバ(使用法はIntroduction参照)

https://cge.cbs.dtu.dk/services/PlasmidFinder/

f:id:kazumaxneo:20180702141809j:plain

データをアップロードする。入力はcontigのFASTAか、NGSのリード(454、illumina、Ion torrent、SOLiD)。

 

指定した相同性の範囲内でplasmidとhitしたcontigが表示される。

f:id:kazumaxneo:20171218114715j:plain

 出力の詳細は公式の説明を参照。

https://cge.cbs.dtu.dk/services/PlasmidFinder/output.php

 

 

ローカルへの インストール

git clone https://bitbucket.org/genomicepidemiology/plasmidfinder.git 
cd plasmidfinder

cd /path/to/plasmidfinder
./INSTALL_DB test_database

# Check all DB scripts works, and validate the database is correct
./UPDATE_DB test_database
./VALIDATE_DB test_database

#Perlbrew is used to manage isolated perl environments. To install it run:
bash brew.sh

#At last PlasmidFinder has several Perl dependencies. To install them (this requires CPAN minus as package manager):
make install


perl plasmidfinder.pl #動作確認

 

plasmidseeker

引用

In silico detection and typing of plasmids using PlasmidFinder and plasmid multilocus sequence typing.

Carattoli A, Zankari E, García-Fernández A, Voldby Larsen M, Lund O, Villa L, Møller Aarestrup F, Hasman H.

Antimicrob Agents Chemother. 2014 Jul;58(7):3895-903. doi: 10.1128/AAC.02412-14. Epub 2014 Apr 28.

 

RNA seq用のde novoアセンブリツール BinPacker

 

RNA-seq法の出現によりmRNA発現レベルに関して前例のない正確さが提供されたため、転写、スプライシング変異および関連する機構の研究方法が大きく変わっている[論文より ref.1]。それらは、レアなスプライシングアイソフォームおよび低発現スプライシングアイソフォームを含む、すべてのスプライシングバリアントの正確な解明を可能にする。これは、ガンを含む異常なスプライシング[ref.1]に関連する様々なヒト疾患のメカニズムを研究するための多くの新しい扉を開いている。 RNA-seq法では、生成されたデータセットの解釈に関連する新たな課題が生じる。 PacBio RS IIシーケンサーからのシークエンシングリードは、複数のエクソンをカバーするのに十分な長さだが、より高いエラー率に悩まされているため、転写再構成の改善には一般的に使用されていなかった。したがって、ショートシーケンシングリードを使ったRNA-seq技術は依然として必要になっている。 1つの大きな課題は、ショートシーケンシングリードを、複数のスプライシングバリアントを含む可能性のある全長転写物に正確に組み立てる方法である。

 文献[ref.4-6]によれば、真核生物遺伝子は複数のアイソフォームを産生することができる様々な選択的スプライシングイベントが存在する。イベントタイプには、エキソンのスキップ、イントロン保持および相互に排他的なエキソンが含まれる。さらに複雑なことに、いくつかのエキソンは、選択的スプライシングプロセス中に転写産物に部分的に関与することがある。一見すると、トランスクリプトームアセンブリはゲノムアセンブリに似ているが、実際には異なり以下の事実により、トランスクリプトームアセンブリはより困難になる:(i)転写産物の発現レベルが非常に低かったり、非常に高かったりし得る[ref.8]。 (ii)各遺伝子座は、通常、種々の選択的スプライシングイベントに起因して複数の転写産物を産生する。(iii)低い発現レベルを有するいくつかの転写産物は、シーケンスエラーのために沈み得る[8,9]。したがって、トランスクリプトームアセンブラの成功は、これらの困難をすべて克服し、可変長、発現レベルおよびノイズの全長転写物を回収することができるもののはずである。

 トランスクリプトームアセンブリの計算ストラテジは、一般にab initioとde novoの2つのカテゴリーに分類される[ref,1,8]。リファレンスゲノムが利用可能である場合、Cufflinks [ref.10]やScripture [ref.11]のようなab initioアプローチは、通常RNA-Seqのリードをリファレンスゲノムにマッピングすることから始まり、オーバーラップアライメントを持つリードを結合グラフにマージする。よく研究された最小コスト最小経路カバーモデルは、RNA-seqデータセットを説明する最小セットの経路抽出にわずかに用いられる。非常に最近公開された(論文執筆時点)ab initioアセンブラであるStringTie [ref.12]は、RNA-Seqのりーどをリファレンスゲノムにマッピングし、最大スプライシンググラフを作成し、最大フローネットワークモデルを使用して転写産物を組み立てる。 ABySS [ref,13]、SOAPdenovo-Trans [ref.14]、Oases [ref.15]、IDBA-Tran [ref.16]などの新しいアプローチは、転写産物をリファレンスゲノムにマッピングすることなく直接使用する。これは、リファレンスゲノムが利用できないか、高度に断片化しているか、ガン組織のように実質的に変化している場合に重要になる。ゲノムアセンブリで使用される技術に基づいて開発されたこれらのデノボアプローチは、一般的なトランスクリプトームアセンブリ問題のすべてを解決するものではない[ref.7]。 de novoトランスクリプトームアセンブリ用に特別に設計されたTrinity [ref.8]は、最新のトランスクリプトームアセンブラの状態を大幅に改善した。これはショートリードをオーバーラップを通してコンティグへ拡張し、コンティグをgraphに接続し、このgraphからパスを抽出してbrute-force 戦略に基づいてスプライシングバリアントを構築することから始まる。Trinityは、ゲノムアセンブリ技術に根ざした以前のデノボアセンブラを改良しているが、そのソリューションを最適化するための適切なモデルを導入していないし、シーケンシングカバレッジのデプス情報もアセンブリ手順に組み込んでいない。しかしTrinitiyの著者らは転写産物中の異なるコード領域のカバレッジデプスの類似性の利用が有用であり得ることを指摘している。この目的のために、私たち(この論文の著者ら)は最近、Cufflinksで使用される技術を使用してトリニティの限界を克服できるように、CufflinksとTrinityの間を「橋渡しする」新たなトランスクリプトーム・アセンブラー、Bridger [17]を発表した。 Bridgerは、カバレッジ情報を適切なモデルを介して組み立て手順に組み込んでいるが、(1)Bridgerでは、weightがいくらか任意に(arbitrarilyに)定義されている:(ref.2) イン・エッジとアウト・エッジの両方を有するノードは、転写産物の終わりである可能性がない。したがって、改善の余地がまだ残っている。

 本論文では、Bridger [ref.17]で用いられている手法を用いて作成されたスプライシングgraph上のアイテムの軌道の集合を追跡することで問題を再構築することにより、フルサイズの転写産物を組み立てる新しいde novoアルゴリズムBinPackerを報告する。スプライシングgraph上のアイテムの軌跡の集合は、異なるサイズのアイテムの所定の数をパックするように定義された従来のビンパッキング問題とは異なるビンパッキング問題の一連の変形を解くことによって達成することができる可能な限り所定のサイズのビンに分割する。それぞれのビンは、ビンのサイズを超えないサイズの合計を持つアイテムのみを保持できる。Bridgerを実際のデータセットとシミュレートされたデータセットを使い、競合するアセンブリツールTrinity[ref.8]、ABySS [ref.13]、Trans-ABySS [ref.18]、SOAPdenovo-Trans [ref.14]、Oases [ref.15]、IDBA-Tran [ref.16]と比較した。シミュレーションデータセットは結果セクションで説明したように生成した。リアルデータセットは、2つの標準RNA-seqデータセット、1つのイヌと1つのヒト、および1つのstrand specificマウスRNA-seqデータセットを含む3つのデータセットが使用された。比較結果は、リアルデータとシミュレート両方で、BinPackerがトランスクリプトーム・アセンブラの評価に一般的に使用されているスタンダードな評価基準で比較されたアセンブラのほとんどを凌駕することを示していた。さらに驚くべきことに、犬のデータセットでは、最も最近にリリースされたab initioのアセンブラであるStringTie [ref.12]のよりも優れていた。

 

f:id:kazumaxneo:20180702125232j:plain

BinPackerアセンブリグローチャート。supplementaryより転載。

 

BinPacker に関するツイート


 インストール

ビルド済みのバイナリとソースコードがダウンロードできる。

SourceForge

https://sourceforge.net/projects/transcriptomeassembly/files/BinPacker_binary.tar.gz/download

解凍して中に入る。

cd BinPacker_binary/
./update
./BinPacker -h #ヘルプ

$ ./BinPacker 

    

** Error: data type is not specified! Please type -h option for help! **

    

[uesaka@cyano BinPacker_binary]$ ./BinPacker -h

    

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

    

BinPacker v1.0 usage:

    

** Required **

    

-s <string>: type of reads: ( fa or fq ).

    

-p <string>: type of sequencing: ( pair or single ).

    

If paired_end reads:

   -l <string>: left reads.

   -r <string>: right reads.

    

If single_end reads:

   -u <string>: single reads.

    

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

    

** Options **

    

-d: remove duplicated transcripts, at least one of whose ends has both in- and out- edges.

    

-o <string>: name of directory for output, default: ./BinPacker_Out_Dir/

    

-m <string>: strand-specific RNA-Seq reads orientation, default: double_stranded_mode.

             if paired_end: RF or FR;

             if single_end: F or R.

    

-k <int>: length of kmer, default:25.

    

-g <int>: gap length of paired reads, default: 200.

    

-S <int>: minimum coverage of kmer as a seed, default: 2.

    

-E <float>: minimum entropy of kmer as a seed, default: 1.5.

    

-C <int>: minimum coverage of kmer used to extend, default: 1.

    

-N <float>: minimum entroy of kmer used to extend, default: 0.0.

    

-J <int>: minimum of the coverage of a junction, default: 2.

    

-v: report the current version of BinPacker and exit.

    

** Note **

    

A typical command of BinPacker might be:

    

BinPacker -s fq -p pair -l reads.left.fq -r reads.right.fq

    

(If your data are strand-strand, it is recommended to set -m option.)

    

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

 BinPackerはバイナリファイルではなくシェルスクリプト

 

 dockerイメージも1つアップされている(リンク)。

docker pull ycogne/binpacker
#ホストの/Users/user/docker_shareとコンテナの/homeを共有ディレクトリとして起動。
docker run -i -t -v /Users/user/docker_share/:/home ycogne/binpacker
cd /usr/bin/BinPacker_1.0/

 

ラン

テストラン。

cd sample_test/
BinPacker -s fq -p pair -l reads.left.fq -r reads.right.fq
  • -s     type of reads: ( fa or fq ).
  • -p     type of sequencing: ( pair or single ).
  • -l      left reads.

  • -r      right reads.

テストの計算は数秒で終わる。出力 BInPacker_Out_Dir/にBinPacker.faとlogができる。

 

supplementaryのtable A-Dには、フルサイズまでアセンブリできたtranscriptsの数がまとめられています。 supplementaryには他のRNAアセンブリツールで使用したオプションなどもまとめられています(pubmed)。

 

引用

BinPacker: Packing-Based De Novo Transcriptome Assembly from RNA-seq Data.

Liu J, Li G, Chang Z, Yu T, Liu B, McMullen R, Chen P, Huang X.

PLoS Comput Biol. 2016 Feb 19;12(2):e1004772.

 

RNA seq用のターゲットアセンブリツール Kollector

 

 非モデル生物のための高品質のリファレンスゲノム配列の作製は、特に大きなゲノム(> 1Gbp)では依然として挑戦的な取り組みである。このようなプロジェクトでは、デノボでの全ゲノムアセンブリは、通常、数種の異なるタイプのDNAライブラリーの数十億のシーケンシングリードを必要とする。これらの大量のデータを処理し、それらを使用してゲノムを組み立てるには、通常、高性能コンピューティング環境、重要な専門知識、および特殊なソフトウェアへのアクセスが必要となる(論文より Nagarajan and Pop、2013)。リファレンス配列を生成するための魅力的な代替法は、目的の遺伝子/転写物配列のターゲットアセンブリによって達成され得る。たとえほとんどのトランスクリプトーム配列情報をもたない種であっても、関連する生物からの相同遺伝子配列のような、デノボアセンブリを助けるために使用され得る既存の配列が存在する可能性がある。これらのデータの利用はアセンブリ問題を局在化するのに役立ち、所望の配列(例えば、遺伝子領域)が完全に再構築されることを保証する。このメリット、全ゲノムアセンブリに比べて複雑さの低減および計算コストの削減である。しかしながら、実際には、ターゲット配列に関連するリードを同定するための計算コストは​​、ターゲット内に見出されない変異および新規配列のために依然として困難である。

 特定のターゲットの再構成のための最初のソリューションは、アライメントフリーのターゲットデノボアセンブリソフトウェアであるTASRで達成された(Warren and Holt、2011)。この方法に続いて、プロセスを導くためにリードアライメントを使用するMapembler(Peterlongo and Chikhi、2012)があり、よりメモリ効率が高く高速な代替手段を提示した。これらの先駆的なターゲットアセンブリテクニックは、もともと大きなショットガンデータから特定の転写バリアント、融合転写物または遺伝子を再構築するように設計されており、現在、ヒトの健康に関する研究に応用されている(Brown et al、2014; Warren et al、2012)(一部略)。

不完全な領域を再構築しギャップを埋めるするために、現代の方法のほとんどは、反復リードを呼び出すプロセスを採用している。 MITObim(Hahn et al、2013)、GRAbB(Brankovics et al、2016)およびaTRAM(Allen et al、2015)は、不完全な配列から新規領域を拡張するために初期の反復ステップに使われたリードを使う。 MITObimは、ミトコンドリアゲノムを組み立てるように設計されており、ターゲットと31-mer(長さ31の部分配列)の配列を共有するリードからターゲットが再構成されるまでシーケンスリードセットを何度も循環させる。 GRAbBは同様の方法で動作するが、一度に複数のターゲットのリードを使うように設計されているため、MITObimを計算上実行できなかった(Brankovics et al、2016)。最後に、aTRAMは関連するゲノムからオーソログを組み立てるために設計されており、BLAST(Altschul et al、1990)を利用して配列の一部を索引付けするため、複数回の反復では、より高いメモリ使用量になる。これらのツールのそれぞれは、複数のリード採用サイクルの後に、確立されたアセンブリツール(例えばVelvet; Zerbino and Birney、2008)を使用してアセンブリを実行し、その後の繰り返しにバイトシーケンスとして使われ、拡張される。

 RNA-Seq技術とde novo assembly tools(Grabherr et al、2011; Peng et al、2013; Robertson et al。、2010)の進歩により、非モデル生物からの高品質トランスクリプトームがますます利用可能になり、ターゲットアセンブリの貴重なリソースになっている。この論文では、全トランスクリプトームアセンブリを使用して全ゲノムショットガンシーケンシングリードをフィルタリングおよび分類し、対応する遺伝子座のデノボアセンブリローカライズすることができる、アライメントフリーのターゲットアセンブリパイプラインであるKollectorを紹介する。パイプラインは、BioBloom Tools(BBT)(Chu et al、2014)内に実装されたプログレッシブブルームフィルタと呼ばれる新しいデータ構造を使用して、ターゲット座に関連するシーケンスリードを収集し、de Bruijn  graph(Pevzner et al。、2001)アセンブラであるABySS(Simpson et al。、2009) でアセンブルする。 Kollectorはイントロン領域を反復的に拡張することができ、実際にはプログレッシブブルームフィルタを貪欲に配置することによって、以前の方法より反復回数は少ない。我々(著者ら)は、Caenorhabditis elegansとHomo sapiensの遺伝子を使い、Kollectorと公表された4つのアセンブリツールを比較による相対的効果を実証する。また、応用としてKollectorの比較ゲノムミクスおよびガンゲノミクスへの使用事例を示す。

 

 インストール

cent OSにインストールした。

 

依存

GMAP以外はbrewで導入できる 

brew tap brewsci/science
brew tap brewsci/bio
brew install abyss samtools bwa biobamtools

GMAPはgithubリンク)からダウンロードしてビルドする(./configure && make)。

 

本体 Github

https://github.com/bcgsc/kollector

ダウンロードすれば使える。

git clone https://github.com/bcgsc/kollector.git
cd Kollector/bin/
./kollector.sh #ヘルプ

$ ./kollector.sh 

Error: number of file args must be  3

Usage: kollector.sh [options] <seed> <pet_read1.fq> <pet_read2.fq> 

 

Description:

 

Do a targeted assembly  using ABySS. The input files are  

PET sequencing reads which must be a FASTA/FASTQ pair and a 

seed sequence FASTA file to recruit reads. The input files may be gzipped.

 

 

AbySS(1.5+),BioBloom Tools and GMAP should be in your path.

 

Options:

    

    -h        show this help message

    -j N      threads [1]

    -r N      min match length for tagging reads. Decimal value are

              the proportion of the valid k-mers and integer values

              will require that minimum number of bases to match [0.7]

    -s N      min match length for recruiting reads [0.50]

    -k N      k-mer size for ABySS contig assembly [32]

    -K N      k-mer size for read overlap detection [25]

    -n N      max k-mers to recruit in total [10000]

    -o FILE   output file prefix ['kollector']

    -p FILE   Bloom filter containing repeat k-mers for

              exclusion from scoring calculations; must match

              k-mer size selected with -K opt [disabled]

    -B N      pass bloom filter size to abyss 2.0.2 

              (B option, to be written: ex - 100M, optional)

kollector.sh、kollector-multiple.sh、kollector-extract.sh各々が使われるので、binにパスを通しておく。

export PATH="$HOME/kollector/bin:$PATH" 

 

ラン

 テストラン

cd kollector/test/

#ツールとC.elegansのテストデータのダウンロード、およびアセンブリzshとlinuxbrewが必要(apt installで導入可)。
make

 

 

kollector.sh <params> seed.fa read1.fa read2.fa
  •  <seed.fa> is the input transcript sequence in a form of FASTA file to recruit reads.
  • <read1.fa> and <read2.fa> are the PET sequencing reads and could be in a form of FASTA/FASTQ files. All the input files can be gzipped.

 

kollector.shを繰り返しランするには、kollector-multiple.shを使う。kollector-multiple.shは、-rを徐々に下げ、 kollector.shを繰り返し回す。アセンブリが成功したターゲットはその都度除外されていく。デフォルトでは試行回数は5回。

kollector-multiple.sh
  • -r <N>    min match length for tagging reads. Decimal value are the proportion of the valid k-mers and integer values will require that minimum number of bases to match [0.7]

    -max_iterations < N>  number of iterations to be performed [5]
  • -decrement <N>  decrement of the r parameter in each iteration [0.1] 

テスト時はエラーを起こした。

 

 

引用

Kollector: transcript-informed, targeted de novo assembly of gene loci.

Kucuk E, Chu J, Vandervalk BP, Hammond SA, Warren RL, Birol I.

Bioinformatics. 2017 Jun 15;33(12):1782-1788. 

 

SVシミュレーションや、SVのマージ、レポート生成ができる SURVIVOR

2019 10/29 インストール修正

2020 1/6 追記

 

 一塩基多型(SNP)、小さな挿入 - 欠失事象(indels)、トランスポゾン挿入および大きな構造変化(SV)を含む、様々な遺伝的変化が生物種に影響し得る。欠失、重複、挿入、逆位および転座を含むSVは、タイピングするのが最も困難であり、結果として最もよく説明されていない。

 それにもかかわらず、SVは様々な生物学的過程に強い影響を与えることは明らかである。特に、コピー数変化(CNV)は、農業上重要な形質、様々なヒト疾患、微生物、植物および動物の定量的形質に影響を及ぼす(論文より ref.1,2,3,4,5)。逆位は生殖的隔離(wiki)に影響することが知られている(ref.6,7,8,9,10,11,12,13)および他の進化的過程、例えば組換え(ref.8)や種間のハイブリダイゼーション(ref.14)。

我々(著者ら)は、最近、分裂酵母Schizosaccharomyces pombeをpopulation genomicsおよびquantitative trait analysis(QTL wiki)のモデルとして開発し始めた(ref.6,7,16,17,18)。このモデル生物は、小さく、しっかりアノテーションされた一倍体ゲノム(ref.19)を持ち、遺伝子操作およびハイスループット表現型解析のための豊富なツール(ref.20)があり、ゲノム規模および遺伝子中心のデータ(ref.21,22,23)の大量のリソースを持つという利点を兼ね備えている。

 分裂酵母のこれまでの分析では、自然発生的および人工的な逆位および相互転座の両方を記述し始めている(ref.6,7,18)。 SVの証拠とそのモデル種の影響を考慮すると、我々(著者ら)はSVの体系的調査が生物学的影響の理解を促進するであろうことを認識した。ここでは、161の分裂酵母ゲノムの最新の利用可能性と、定量的形質および生殖分離に関する広範なデータを利用して、S.pombeにおけるSVの性質および効果を説明する。

著者らは、SVが様々な量的形質および本質的な(もともと備わっている:内因性の)生殖隔離に強い影響を及ぼすことを示す。それらは、形質変動の平均11%に寄与する(より豊富なSNPが平均24%寄与する)が、CNVが最大の影響をもたらす。著者らは、CNVがクローン集団内で一過性であり、しばしばSNPによってうまくタグ付けされていないことを示す。また、再構成(逆位および転座)は生殖隔離に寄与するが、CNVはそうしないことを示す。

Methodより

 著者らは、ショートリードシーケンスデータのSVを評価するためいくつかのモジュールを含むSURVIVORツールキットを開発した。第一のモジュールは、リファレンスゲノムファイル(fasta)と各SVの数とサイズの範囲(挿入、削除、複製、逆位および転座)を与えられたSVをシミュレートする。リファレンスゲノムを読み取った後、SURVIVORは、提供されたパラメータに従ってSVの位置およびサイズをランダムに選択する。続いて、SURVIVORはリファレンスゲノムをそれに応じて変更してプリントする。さらに、SURVIVORは、シミュレートされたSVの位置を報告する拡張ベッドファイルを提供する。

 第2のモジュールは、Variant Call Format(VCF)ファイルおよび任意の既知SVリストに基づいてSVコールを評価する。 SVが正しいと同定されるのは、(i)同じSVサブタイプ(例えば、欠失)である、 (ii)同じ染色体上に報告される、(iii)シミュレートされ同定されたSVの開始および停止座標が1kb以内である(ユーザ定義可能)、を全て満たす場合である。

 SURVIVORの 第3モジュールは、3つのSVコーラーのVCFファイルからのコールをフィルタリングして結合するために使用された。著者らの場合、これらのファイルはDELLY、LUMPY、Pindelの結果であった。このモジュールには、メソッド固有の出力フォーマットをVCFフォーマットに変換するメソッドが含まれている。 SVは、3つのVCFファイルのいずれかのみに固有であれば除外された。 2つだった場合、同じ染色体上に存在し、その開始および停止座標は1kb以内であり、それらは同じ型の場合に重複すると定義された。最後に、SURVIVORはフィルタリングされたコールを1つのVCFファイル統合した。 SURVIVORはgithub.com/fritzsedlazeck/SURVIVORから入手できる。

 

SURVIVORは以下のような機能を持つ(Githubより)。

  • Simulate SVs and evaluate existing callers
  • Merge and compare SVs within a sample and among populations/samples.
  • Convert different formats to vcf files
  • Summarize the results within vcf files or results from SURVIVOR.

そのほか、リアルデータ(PacbioやONTのデータ)からエラープロファイルを取得し、ロングリードをシミュレーションする機能をもつ。

 

wiki

https://github.com/fritzsedlazeck/SURVIVOR/wiki

 

インストール

mac os10.13、python3.4.0環境でテストした。

本体 Github

git clone https://github.com/fritzsedlazeck/SURVIVOR.git
cd SURVIVOR/Debug
make

 Anaconda環境ならcondaでも導入可能。linuxのみ

#bioconda (link)
conda install -c bioconda survivor

> ./SURVIVOR

# SURVIVOR

Program: SURVIVOR (Tools for Structural Variations in the VCF format)

Version: 1.0.6

 

Usage: SURVIVOR <command> [options]

 

Commands:

-- Simulation/ Evaluation

simSV Simulates SVs and SNPs on a reference genome.

scanreads Obtain error profiles form mapped reads for simulation.

simreads Simulates long reads (Pacio or ONT).

eval Evaluates a VCF file after SV calling over simulated data.

 

-- Comparison/filtering

merge Compare or merge VCF files to generate a consensus or multi sample VCF files.

genComp Generates a pairwise comparison matrix based on any multi sample VCF file

filter Filter a vcf file based on size and/or regions to ignore

stats Report multipe stats over a VCF file

compMUMMer Annotates a VCF file with the breakpoints found with MUMMer (Show-diff).

 

-- Conversion

bincov Bins coverage vector to a bed file to filter SVs in low MQ regions

vcftobed Converts a VCF file to a bed file

bedtovcf Converts a bed file to a VCF file 

smaptovcf Converts the smap file to a VCF file (beta version)

bedpetovcf Converts a bedpe file ot a VCF file (beta version)

hapcuttovcf Converts the Hapcut2 final file to a VCF file using the original SNP file provided to Hapcut2

convertAssemblytics Converts Assemblytics to a VCF file

 

ラン

推奨パラメータ

https://github.com/fritzsedlazeck/SURVIVOR/wiki/Methods-and-Parameter

1、configファイルの作成。

SURVIVOR simSV parameter_file

parameter_fileファイルができる。中身を確認する。

> cat parameter_file

$ cat parameter_file 

PARAMETER FILE: DO JUST MODIFY THE VALUES AND KEEP THE SPACES!

DUPLICATION_minimum_length: 100

DUPLICATION_maximum_length: 10000

DUPLICATION_number: 3

INDEL_minimum_length: 20

INDEL_maximum_length: 500

INDEL_number: 1

TRANSLOCATION_minimum_length: 1000

TRANSLOCATION_maximum_length: 3000

TRANSLOCATION_number: 2

INVERSION_minimum_length: 600

INVERSION_maximum_length: 800

INVERSION_number: 4

INV_del_minimum_length: 600

INV_del_maximum_length: 800

INV_del_number: 2

INV_dup_minimum_length: 600

INV_dup_maximum_length: 800

INV_dup_number: 2

この定義に従いSVが導入される。変更が必要ならエディタで開いて修正、保存する(順番やスペースをいじってはならない)。

 

2、SVを含むゲノムfastaの作成。先ほど作ったconfigファイルを指定する。SNPs頻度(0-1の間)は0.1にした。

SURVIVOR simSV ref.fa parameter_file 0.1 0 simulated

上のコマンドは以下の順番に記載している。

1: Reference fasta file

2: Parameter file

3: SNP mutations frequency (0-1)

4: 0= simulated reads; 1= real reads

5: output prefix

4番目の数値(0 or 1)はシミュレーションのモードを表す。0はシミュレーションリードを発生させ、それをリファレンス(WTゲノム)にマッピングするモード。1はリアルシーケンスデータをSVを導入したリファレンス(MTゲノム)にマッピングするモード。いずれにしてもSVは検出されるが、挿入と欠失のコールなどは逆転する。

 

human chr20を使いテストしたが、simSVのランが2日経っても終わらなかった。CPU usageは100%のままになっているので、何か計算していると思われるが、logやエラーを立てるフラグがないためどうなっているかわからない。コマンドを勘違いしているのかもしれない。

 

3、コールされたSVの評価。

SURVIVOR eval caller.vcf simulated.bed 10 eval_res

 2のランが終われば追記します。

 

 

他にも以下のようなコマンドがある。

-- Comparison/filtering --

  • merge   Compare or merge VCF files to generate a consensus or multi sample vcf files.
  • filter   Filter a vcf file based on size and/or regions to ignore
  • stats   Report multipe stats over a VCF file
  • compMUMMer   Annotates a VCF file with the breakpoints found with MUMMer (Show-diff).

-- Conversion --

bincov   Bins coverage vector to a bed file to filter SVs in low MQ regions

  • vcftobed   Converts a VCF file to a bed file
  • bedtovcf   Converts a bed file to a VCF file 
  • smaptovcf   Converts the smap file to a VCF file (beta version)
  • bedpetovcf   Converts a bedpe file ot a VCF file (beta version)
  • hapcuttovcf   Converts the Hapcut2 final file to a VCF file using the original SNP file provided to Hapcut2
  • convertAssemblytics   Converts Assemblytics to a VCF file

 

ここではリアルデータのエラープロファイルを取り込んだシミュレーションロングリードを発生させる流れを記載する。

1、リアルデータを読み込みエラープロファイルファイルを作成。

SURVIVOR scanreadsにsamtools viewをパイプしてリードを読み込む。ターゲットは最小1000bpとする(ロングリードの場合、短かすぎるリードには異なるbiasが発生しているかもしれない。)。時間節約のため、viewに"-s 0.1"をつけて例えば10%だけサンプリングする。

samtools view -s 0.1 your_file.bam | ./SURVIVOR scanreads 1000 error_file

error.txtができる。ちなみにONT由来のbamでは動作したが、ショートリードのbamを使うと、パラメータを変えてもsegmentation faultを起こした。

 

2、エラープロファイルに従いリードを発生させる。

指定するパラメータはカバレッジだけになる。

SURVIVOR simreads ref.fa error_file 100 output.fa

ロングリードのFASTAファイルが出力される(fastqは出力不可)。

 

追記

SV callのマージ と視覚化(チュートリアルに載っています)

ls *vcf > sample_files
#1kb以内 のコールを統合、two SV caller以上がサポート、10bp以下除去
SURVIVOR merge sample_files 1000 2 10 0 0 10 merged.vcf

#視覚化
SURVIVOR genComp merged.vcf merged.mat.txt
#Rにて
t=read.table("merged.mat.txt",header=F)
dst <- data.matrix(t(t))
image(1:nrow(dst), 1:ncol(dst),log(dst), col=rev(heat.colors(10)),axes=F, xlab="", ylab="")
grid(col='black',nx=nrow(dst),ny=ncol(dst),lty=1)

f:id:kazumaxneo:20200106212437p:plain



 

引用

Transient structural variations have strong effects on quantitative traits and reproductive isolation in fission yeast.
Jeffares, Daniel C; Jolly, Clemency; Hoti, Mimoza; Speed, Doug; Shaw, Liam; Rallis, Charalampos; Balloux, Francois; Dessimoz, Christophe; Bähler, Jürg; Sedlazeck, Fritz J.
Nature communications, Vol. 8, 14061, 24.01.2017, p. 1-11. DOI:10.1038/NCOMMS14061