Nanoporeのシーケンシングリードをマッピングして分析できるwebサーバー nanopipe

　近年、技術が絶えず改良され科学者や開業医を含む広範囲の顧客にとって利用しやすいためにDNAシーケンシングブームが起きている。メタゲノミクスから植物生理学、医学まで、多くの分野の研究者が、彼らの研究にシーケンシング実験を実施してきた。 Oxford Nanopore Technologies（ONT）は、技術的スキルとバイオインフォマティクスに関する知識を最低限しか必要としない携帯型マシンであるMinIONシーケンサーを導入することで、このプロセスを大幅にスピードアップした。したがって、DNAシーケンシング実験は、野外研究でも、小さな実験室でも、そして診療所での医学的応用のためにも実行可能になってきている[ref.1]。

　NCBI PubMedデータベースには、2009年から2018年の間（2018年9月1日まで）にpublishされた「Oxford Nanopore」という語句を含む261の科学論文が含まれている（その大部分は過去3年間にpublishされた）。これは、ONTシーケンサーの人気の高まりと、過去3年間におけるテクノロジとシーケンス品質の大幅な向上の両方を指している。例えば、MiniIONフローセルのR10バージョンが最近発表され、ホモポリマーストレッチについてもシーケンス品質の向上を約束している。ハイスループットおよびロングリードは、MinIONの多様な応用を可能にする。ウイルス学[ref.2 link、3 pubmed]、植物病理学および農業[ref.4、5]、結核研究[ref.6]、メタゲノミクスおよび食べ物[ref.7 link]、ならびに獣医学研究[ref.8 pubmed]。ポータブルプラットフォームとしては、野外生物多様性研究[ref.9]、現場での患者におけるエボラウイルスの検出[ref.10、11]、および宇宙でのシーケンシング[ref.12, pubmed]。 de novoゲノムアセンブリ、既存のゲノムアセンブリの改良、構造変異および長いリピート配列の発見など、ロングリードシーケンス研究の基本的な用途も忘れない[ref.13–16]。小型のバクテリアゲノムは1回のMinIONリードでカバーできる[ref.17, 18]、高解像度なゲノムアーキテクチャを提供できるため、ONTシーケンスは微生物学研究に適している。

　ONTはユーザーにMinIONの電気信号を一連のヌクレオチドに変換するbase callingを行うために必要なソフトウェアを提供している。それはMinNNOWとオフラインのAlbacoreである。どちらのアプリケーションも複雑で再帰的なニューラルネットワークアルゴリズムを利用している。これは最近コンピュータ科学で非常に人気がある。それはソフトウェアが既存のデータから学び、そのパフォーマンスを改善することを可能にする。base callingプロセスがONTシーケンステクノロジの精度を向上させるための中心であり、そのアルゴリズムが絶えず改善および更新されていることは注目に値する。出力は、FAST5および/またはFASTQファイルのコレクションになる。これらは、バイオインフォマティクスのシーケンス解析に使用されるファイルである。したがって、base calling元は、MinIONからデータ解釈への「ゲート」と呼ぶことができる。それにもかかわらず、ONTが提供する分析ツールの範囲は限られており、ユーザーに任されている一般的な処理を除いて、特定のアプリケーションだけに関係している。例えば、MetrichorプラットフォームをベースとしたEPI2MEソフトウェアスイートには、バーコード分析、メタゲノミクス、抗菌剤耐性分析、およびいくつかの技術的試験のためのアプリケーションが含まれている[ref.19]。最近、いくつかの研究グループがMinION特有のバイオインフォマティクスツールの開発に焦点を当てている[ref.20-22]。これらのほとんどはかなりのバイオインフォマティクス知識を必要とする。これらの条件は、MinIONベースのDNAシーケンスが、情報技術（IT）の経験が少ない医師や研究者にもたらすことができる利益を妨げる。このギャップを埋めるために、MinIONによって生成されたデータを迅速かつ容易に処理し、必要に応じてさらなるバイオインフォマティクス分析に必要なファイルを提供することができるWeb駆動の自動パイプラインでNanoPipeを開発した。

A schematic representation of the NanoPipe workflow. 論文より転載。

usage

Institute of Bioinformatics WWU Münster

使い方

http://bioinformatics.uni-muenster.de/tools/nanopipe2/generate/index.pl?lang=en

にアクセスする。

Target Fileにリファレンスゲノム、Quey Fileにシーケンスデータを選択する。複数シーケンスデータがある場合は圧縮して送る（zip archive.zip file1 file2 ...）。ファイルサイズは最大3GB。

いくつかのモデル生物のリファレンス配列はすでにサーバーに準備されている。Upload Fileをクリックして一覧から選択する。

あとはマッピングに使われる LASTのパラメータを選択してランする。テスト時は数時間で結果が得られた。

 
デモラン

View Testcaseをクリックすると、デモランの結果を表示できる。

Overview

マッピング率、LASTのパラメータなどが表示される。

Mapping distribution  

どのchrmosomeにどれだけマッピングされたか確認できる。

Alignment distribution

リードのアライメントされた長さの分布がまとめられる。

BAM FIles

ダウンロードして使う。bamとbam.bai、リファレンスfastaをダウンロードできる。

Target ID

Target IDのタブで表示するリファレンスを選ぶ（これより下は１つのリファレンスのみ分析対象になるため）

Nucleotide Plots

各ポジションの塩基が表示される。下のボタンで上流/下流に移動する。

 Consensus

コンセンサス配列が表示される。コンセンサス配列はDownloadボタンからダウンロードできる。

コンセンサス配列は、そのポジションでもっと割合が多い塩基がアサインされる（割合が80%以上になること）。ナノポアはエラー率が高いので、カバレッジが10以下は一律ギャップ（-）になる。

Polymorphisms

推定SNPが表示される。

以下の３つの条件で残ったものにウエイトがつけてまとめられる。

1. The coverage of the target nucleotide is lower than 80 % of the total coverage at that position.
2. For non-target nucleotides at that position: The coverage must be greater than 20% of the total coverage.
3. A position must have a coverage of more than 30% of the general target's maximum coverage. E.g.: The maximum read coverage of a target position is 5000, thus, SNPs with a coverage below 1500 are discarded.

ヒトゲノムとPlasmodium falciparum（熱帯熱マラリア原虫）は、dbSNP / PlasmoDBのデータベースとマッチするかも表示される。

Alignments

リファレンスとコンセンサス配列のpairwiseアライメント結果が表示される。

引用

NanoPipe—a web server for nanopore MinION sequencing data analysis
Victoria Shabardina, Tabea Kischka, Felix Manske, Norbert Grundmann, Martin C Frith, Yutaka Suzuki, Wojciech Makałowski
GigaScience, Volume 8, Issue 2, 1 February 2019

関連

スモールゲノムを可視化したり、複数ゲノムを比較して似た領域、異なる領域を可視化できる Gview

　グラフィックなゲノムマップは、ゲノムの特徴および配列の特徴を評価するために広く使用されている。 CGView（Circular Genome Viewer）ソフトウェアファミリーは、バクテリア、オルガネラ、ウイルスのゲノムマップを生成するためのツールの人気のあるコレクションである。 このレビューでは、オリジナルのCGViewプログラムの機能と、それ以降のコンパニオンアプリケーション（CGViewサーバやCGView Comparison Tool（紹介）など）の機能について説明する。 また、GView、CGViewのグラフィカルユーザーインターフェイス対応書き換え、および比較ゲノムのコレクションに対して共有または固有のゲノム領域を識別するためのいくつかの統合分析を提供するGViewサーバーについても説明する。 最後に、使いやすさを増しながら新しい機能を追加することを目的とした、CGViewに関連した現在の開発努力についていくつか述べておく。

Gview HP

https://www.gview.ca/wiki/GView/WebHome

Image Gallery

GView

インストール

https://www.gview.ca/wiki/GViewDownload/WebHome

このリンクからGView 1.7をダウンロードする。

java -jar gview.jar -h

$ java -jar gview.jar -h

usage: java -jar gview.jar [OPTION]...

-h, --help                  Print the usage.

-H, --height <integer>      Specifies the height in pixels for the image.

-c, --centerBase <integer>  Specifies the base pair value to center on.

-i, --inputFile <file>      The input file to parse.

-l, --layout <layout>       The layout of the input, 'circular' or 'linear'.

-f, --imageFormat <format>  The format of the output: [png, jpg, svg, svgz] (default:png)

-o, --outputFile <file>     The image file to create.

-v, --viewer                Loads up the input file in the genome viewer.

-s, --style <file>          The file to read the style information from.

-b, --birdsEyeView          Displays a birds eye view of the map.

                              This is ignored if -v is not used.

-g, --gffFile <file>        Specifies a .gff file with additional annotations to include.

                              Can be used multiple times for multiple files.

-W, --width <integer>       Specifies the width in pixels for the image.

-z, --zoomAmount <real>     The factor to zoom in by.

-t, --threads <integer>     The number of threads to use.

                              The default is the maximum threads the JVM can acquire.

-q, --quality <quality>     The initial rendering quality of the display, 'low' or 'high'.

                              The default quality is high.

--version                   Print version information.

Example:  java -jar gview.jar -i example_data/NC_007622.gbk -s example_styles/basicStyle.gss -l linear -W 1200 -H 900 -f png -o image.png

Example:  java -jar gview.jar -i example_data/NC_007622.gbk -s example_styles/gssExample.gss -l circular -v -b

またはjava webスタート版（リンク）を起動する。

実行方法

方法１

起動し、javaのGUIウィンドウ（ガワ）で操作する。

java -jar gview.jar 

起動直後

Filesにgenbankファイルを指定する。拡張子は".gbk"でないといけない。circularかlinearかを選択してBuildMapをクリックする。

可視化された。

linear表示

上のメニューで拡大縮小したりスクロールできる。左端の二つのボタンが拡大縮小になる。目のアイコンはバードビューで自由に拡大縮小出来る。目のアイコンの左隣のFit map to screenボタンをクリックすると、ウィンドウサイズに合わせ画像サイズがめいいっぱいの大きさに調節される。右端のボタンで画面スクロールできるが、あまり調整できないので、Fit map to screenボタンの方が使いやすい。

方法２

コマンドでgenbankファイルを指定する。exampleのgenbankファイルを使ってみる。

java -jar gview.jar -i08-5578.gbk -l linear -W 1200 -H 900 -f png -o image.png

exampleのgenbankファイル（ダイレクトダウンロードリンク）

GView Style Sheet (GSS)のファイルを指定すると、細かな設定ができる。exampleのgssファイル（ダイレクトダウンロードリンク）をダウンロードして指定してみる。

wget https://www.gview.ca/pub/GViewQuickStart/GViewWebStart/basicstyle.gss
java -jar gview.jar -i08-5578.gbk -f basicstyle.gss -l circular -W 1200 -H 900 -f png -o image1.png

wget https://www.gview.ca/pub/GViewQuickStart/GViewWebStart/gssexample.gss
java -jar gview.jar -i08-5578.gbk -f gssexample.gss -l circular -W 1200 -H 900 -f png -o image2.png

 image1.png

image2.png

GSSファイル指定なし

方法３

オンラインで使用する。https://server.gview.caにアクセスする。

ガイド

GView Server

タイプを選ぶ。

え

以下の９つのタイプがある。

1. Display genome features    Upload a rich sequence file (GenBank or EMBL) to display features in an interactive GView map. 
2. BLAST atlas    Create a BLAST atlas by uploading a reference genome and one or more related genomes. Regions will be displayed where there is similarity between the reference genome and one of the related genomes.
3. Core genomes   View the core regions of a genome by uploading a reference genome and a query file containing multiple related genomes. This will highlight the regions in the reference genome that also exist in all of the query genomes. 
4. Accesory genome   View the accessory regions of a genome by uploading a reference genome and a query file containing multiple related genomes. The resulting regions displayed will be those that had BLAST hits to some sequences in the query file, but not to every sequence in the query file.
5. Unique genome   View the unique regions of a genome by uploading a reference genome and a query file containing multiple related genomes. The resulting regions displayed will be those that had no BLAST hits to any sequence in the query file.
6. Signature genome   Find the signature regions of a genome by uploading a reference genome and a set of inclusion and exclusion genomes. The inclusion genomes should be very similar to the reference genome. The exclusion genomes should be similar to the reference genome, but lacking signature characteristics. The resulting regions displayed will be those that are found in the inclusion genomes but not in the exclusion genomes. 
7. Pangenome analysis   Create a pangenome by uploading a collection of similar query genomes. A seed genome is selected and other query genomes are compared to the seed to locate the unique regions. 
8. Reciprocal blast   Upload a reference sequence file and a query sequence file to perform a reciprocal BLAST and display features in an interactive GView map.
9. Custom analysis   Create a fully customized GView analysis.

詳細はリンク先のAnalysis typeを参照して下さい（リンク）。

ここでは"Display genome features"と"Pan genome"を試してみる。

まずDisplay genome featuresの流れを確認していく。TypeはDisplay genome featuresを選択、リファレンスのgenbankファイルを指定し、continueボタンをクリック（注 いずれのプログラムのランにもgenbankファイルが必要）。

メールアドレスを記入するとjob完了時に通知してくれる（未テスト）。

表示する項目を指定したり、色を変更する。

GC content表示にチェックをつけた。

featureの色を変えるには、上の画像の下の方にあるカラーボックスの左隣のテレビのようなアイコンをダブルクリックする。↓のウィンドウが出現するので、色を指定する。

okを押して閉じる。

最後にcompleteを押すとジョブが開始される。

テスト時は10分ほどで終わった。

画像をダウンロードしたり、ここからjava webスタートを立ち上げて、ローカルのGview.jarで編集することもできる。

次に"Pan genome"を試してみる。

Pan genomeを 選択。Pan genomeモードでは親ゲノムというものはないので、そのままContinueをクリック。

step2ではgenbankファイルを追加していく。

緑のボタンを押し、第２、第３..と、全ゲノムのgenbankを追加していく。最後にContinueをクリック。

step3でパラメータを設定してジョブを開始する。

終わると、比較した図が出力される。

java webスタートを立ち上げて、開き直した。表示項目を変更したい時にも便利です（Style -> Style Editor）。

exampleには様々な例が載ってますね。

https://server.gview.ca/examples

引用

Visualizing and comparing circular genomes using the CGView family of tools
Paul Stothard, Jason R. Grant, Gary Van Domselaar
Briefings in Bioinformatics, bbx081, Published: 26 July 2017

関連

editorial要約 Improving the usability and archival stability of bioinformatics software

2019 3/2 文章修正

　ゲノミクスおよびシーケンシング技術の急速な進歩は、圧倒的な量と多様性の新しいソフトウェアツールとしてパッケージされた分析アルゴリズムをもたらした[ref.1]。 このような計算ツールは、ライフサイエンスや医学研究者がますます複雑化するデータを分析し、困難な生物学的問題を解決し、そして新しい臨床翻訳に不可欠な基礎を築くのに役立ってきた。 確かに、ヒトゲノムの初期シーケンスから最新のハイスループットシーケンスデータの解析までシーケンスデータ解析のすべてのフェーズは、バイオインフォマティクスツール[ref.2]に依存している。

• 開発されたソフトウェアを使用可能にし[ref.3]、ソフトウェアツールにアクセスするためのURLを統一することがますます重要になってきている。
• 一貫して使用可能でアクセス可能なソフトウェアとは、同じ計算ツールを実行することで研究によって生成されたデータを再生成できる能力を持つものとして定義される。
•  さらに、ソフトウェアツールの使いやすさ、つまり「使いやすさ」が重要であり、アカデミアで開発されたソフトウェアツールは、産業界で開発されたツールと比較して、ユーザーにとって使い勝手がよくないことが度々ある。
•  毎年多数のツールがリリースされているため、これらの問題はソフトウェアツールの科学的有用性を制限する可能性がある。
• 計算生物学におけるこの問題の規模はまだ見積もられていないが、バイオインフォマティクスコミュニティは、不十分に維持または実行されたツールは「ビッグデータ」ドリブンの分野における進歩を最終的に妨げるであろうと警告している[ref.6 link]。
• サイエンスソフトウェアをうまく​​実装して幅広く配布するためには、数多くの独自の課題がある。 第一に、アカデミアでは、ソフトウェアツールは大学院生またはポスドク研究者からなる小グループによって開発されている[ref.7]。 これらのグループは、ソフトウェア開発グループと比較すると、短期間であまり包括的なトレーニングを受けずにソフトウエアを実装している。
• さらに、安定したサポートがあっても、ソフトウェアツールのユーザビリティは保証されない。これは、ユーザがツールをインストールして実行する能力として定義される。
• ソフトウェアの使用可能性および計算ツールの保存安定性の制限により、最終的にオリジナルpublishの結果を再現する能力が損なわれる。
•  アカデミックの世界では、開発者はソフトウェアエンジニアリング、特に専門的なユーザーエクスペリエンスとクロスプラットフォーム設計に関する正式なトレーニングを受けていないことがよくある。
• 多くの計算生物学ソフトウェア開発者は、ツールをインストールして実行するためのユーザーフレンドリーなインターフェースを提供するためのリソースを欠いている。
• 使いやすいツールの開発は、「依存関係」と呼ばれる、事前にインストールする必要がある多くのツールがサードパーティ製ソフトウェアに依存しているため、さらに複雑である。
• 現在、計算生物学の分野では、エンドユーザーが簡単にインストールできる十分に標準化されたアプローチが欠けている。
• ソフトウェア開発に対する制度的支援の欠如はこれらの課題を悪化させる。 政府機関は新しい計算方法に興味を持っているが、既存のツールの継続的な開発と維持のための資金は不十分である[ref.8]。
•  広く使用されているツールであっても、ツールの開発や提供を中止したりするなど、長期間のメンテナンスに必要なfundを突然失う可能性がある。
• さらに、ソフトウェア開発は学術的な採用や昇進にインセンティブを受けていない。むしろpublicationと資金調達に焦点を当てている。

•  2000年から2017年までにpublishされた24,490のオミックスソフトウェアリソースの実証分析では、全オミックスソフトウェアリソースの26％が、原著論文に掲載されているURLを通じて現在アクセスできないことを示している。
•  2000年から2017年にかけてアーカイブ不安定なリソースの割合が大幅に減少したが、そもかかわらず、研究期間中、毎年200のアーカイブ不安定なリソースが公開されている。
• 計算ツールの使いやすさとアーカイブの安定性を向上させるには、アカデミアのソフトウェア開発者間での協調的な努力が必要になる。 ソースコードをホストするように設計されたGitHubやSourceForgeなどのWebサイトでソフトウェアツールをホストすることを勧める。
•  著者らの研究では、そのようなオンラインリポジトリに読者を誘導するソフトウェアツールのURLは、高いアクセシビリティを持っている。 GitHubへのリンクの99％およびSourceForgeへのリンクの96％がアクセス可能だが、他の場所でホストされているリンクは72％のみがアクセス可能である。
• バイオインフォマティクスコミュニティは2001年以来これらのWebベースのサービスを使用しており、これらのリポジトリでホストされているソフトウェアツールの割合は、2012年の5％から2017年の20％に大幅に増加している[ref.9]。
• さらに、公開されているオミックスソフトウェアツールのインストール可能性についての体系的な評価は、いくつかの改善分野を示唆している。 分析に含まれた99のランダムに選択されたツールのうち、49％が「インストールが困難」と判断され、インストールに15分以上かかり、28％のツールが2時間以内にインストールされなかった。
• さらに、インストール可能性はソフトウェアツールの人気に影響を与えることがわかった。 正常にインストールされたツールは、2時間以内にインストールできなかったツールと比較して、引用数が大幅に増えた。
•  調査対象のツールをインストールするには、平均8つのコマンドが必要だったが、ユーザーマニュアルの上では平均3.9のコマンドしか示されていなかった。
•  調査されたソフトウェアツールのいくつかはパッケージマネージャを通して利用可能だった。それはユーザがインストール、アップグレードと構成を自動化することを可能にする。
• よく管理されたパッケージマネージャ（例：Bioconda [ref.10]）を通じて入手可能なツールは常にインストール可能である一方、パッケージマネージャなしのツールは、調査したケースの32％で問題を起こしやすいことがわった。
• バイオインフォマティクス研究者はより複雑化するデータセットと問題に取り組んでおり、我々のコミュニティは開発、ピアレビュー、そして公開されたソフトウェアパッケージに対して厳格で標準化されたアプローチを採用する必要がある。
•  アーカイブの不安定性に対する多くの解決策がすでに利用可能である。 たとえば、GitHubやSourceForgeなどのアーカイブ上安定したサービスでバイオインフォマティクスソフトウェアパッケージをホスティングすると、omicsツールの長期的なアクセス性が大幅に向上する。
• 開発者は、依存関係のサードパーティソフトウェアパッケージをダウンロードしてインストールできる使いやすいインストールインターフェイスを作成して提供できる。 あるいは、開発者はBioconda[ref.10]のようなパッケージマネージャでツールをラップすることができる。
• 期待される結果のデータセットが含まれていると、ツールを実験データに対して実行する前に、ツールが正常にインストールされ、正しく機能し得るか確認することができる。

• ジャーナルは、査読プロセス中にソフトウェアのユーザビリティとアクセシビリティに対する厳密で標準化されたアプローチを奨励する必要があるかもしれない。
•  レビュー担当者は、ソフトウェアツールを説明する文書に、インストールスクリプト、テストデータ、およびツールのインストールと実行の妥当性を自動的に確認できる機能などの関連項目を含めることを要求する場合がある。
•  ジャーナルは、私たちの分野で増大している別の問題に対処するための補完的なトップダウン戦略を提供することができる。バージョン管理、または多くのバージョンと構成からなる安定したソフトウェアシステムの維持など。

一部略

　バイオインフォマティクス研究およびハイスループットコンピューティング能力の劇的な変化は、いくつかのツールを無関係にする可能性があるが、それでも優れた精度で新しいツールを開発するためのバックグラウンドを提供している。 我々の研究は、使用可能でアーカイブ上安定なバイオインフォマティクスソフトウェアを作り上げることの課題を強調している。 計算生物学ソフトウェア開発の現在のモデルは、研究者が新しいツールを開発し公開することを奨励するが、既存のツールの維持を奨励するものではない。 さらにアカデミアでは、インストールや使用が簡単な計算ツールを開発する動機はほとんどない。 それにもかかわらず、我々（著者ら）の研究結果は、計算生物学の成長にとってソフトウェアの安定性とユーザビリティが重要性であることを示している。 この結果は、ソフトウェアの検証とアーカイブのための標準化されたアプローチに向けた協調的な取り組みを採用するための説得力のあるエビデンスを提供し、生物医学研究コミュニティにおけるソフトウェアツールの開発と保守のための資金とリソースの必要性を強調している。

感想

　私の主観では、ソフトウエアを簡単に導入できるかどうかが、ソフトウエアを自分の研究に取り込む上で最もクリティカルなステップになっていると感じます。例えば、素晴らしいソフトウエアを見つけたとして、その方法を使うと自分のデータはどうなるのか、（胸の高鳴りを抑えつつ）ソースコードのダウンロード&ビルドを始めても、端末の無味乾燥な黒画面にerrorが何度も表示されていると、次第に興味を失っていくからです（*1）。私だけかもしれませんが、すぐにインストールできないとモチベーションは大きく低下します。失敗した時は1行1行エラーメッセージを調べたりマニュアルを読み返して対処する訳ですが、数回トライアンドエラーしても解決出来ず、グーグリングしてあちこちたらい回しになっている内に頭の中に浮かんでくるのは、そもそも私の OSで動くのか？、本当にテストしてる？、初期化するか....、などの疑心疑惑、雑念です。パッケージマネージャや仮想環境による導入がサポートされていれば、このようなインストール時の時間の浪費は大幅に低減でき（*2）、純粋な研究外のタスクに多量の時間を消費しないで済みます（*3）。いったんソフトウエアが導入できれば、"tiny dataset"を使い最小限のオプションでジョブを開始してすぐに結果を得ることができ、十分な思考能力を残したまま結果を解釈できる余裕が生まれます（*4）。よって、導入ステップこそがソフトウエアを自分の研究に取り込む上で一番の障壁になっていると思っています。

この方の発言はツールの使い方やアルゴリズムについての言及かもしれませんが、導入時のトラブルの対処も基本は同じですね

really, the daily life of a bioinformatician is all about googling around and reading manuals and documentation.

— Ming Tang (@tangming2005) March 1, 2019

　このブログを始める前の話になりますが、私の場合、ソフトウエアの導入や実行手順でよくコケていたので、講習会で配布されたマニュアルを使ったり、インストール方法やコマンドについてまとめたwordの書類を使って対処療法的な対応を取っていました。しかしそのやり方では数が増えるにつれて管理が面倒になるのは目に見えていました。実際、徐々に目的の物を探す時間が増え、面倒さを感じる事が増えていました（ローカルマシン内のデータ検索能力はGoogleより大きく劣り、あのHDDにあるのは分かっているに見つけることが出来ない）。また計算マシンが変わると書類にアクセスできなくなる事や、人に教える時に非常に面倒な事もネックでした。こういった経験から、オンラインでアクセスできるindex付きメモがあれば自分の研究や知人のサポートに役立つだろう、加えて、もしかして会ったことのない他の人にとっても有用ではないか？、という気づきを得てこのブログを始めました。

　ブログでツールを管理し始めてからは、ソフトウエアの実行や、他のソフトウエアとの連携が飛躍的に楽になりました。思った以上にアクセスもいただいているようで、初期の動機はそれほど間違いではなかったかなと考えています。今後もどうぞよろしくお願いいたします。

引用

Improving the usability and archival stability of bioinformatics software
Serghei Mangul, Lana S. Martin, Eleazar Eskin, Ran Blekhman
Genome Biology 2019 20:47

PDF link

https://genomebiology.biomedcentral.com/track/pdf/10.1186/s13059-019-1649-8

関連

*1

私の場合、このブログで紹介しようと思って導入を試みたツールも合計200くらい導入に失敗して下書きのままです。失敗の原因は、論文の該当URLが消えている、ソースコードすらダウンロードできない、ソースコードをダウンロードできてもインストールに失敗する、ビルドは成功するがラン中にsegmentation error等を吐きexitしてしまう、エラーlogなくランは終わるが出力がゼロ、など様々です。オンラインで利用できるツールも、既にアクセスできなくなっている、トップページにはアクセスできるがジョブがスタートしない、ジョブはスタートするが終わらない（一番厄介）、などトラブルは尽きません。

*2

依存が増えてくるとパッケージマネージャを使ってもインストールに失敗する事はありますが、エラーメッセージを見ればどのツールとどのツールがコンフリクトしているか分かります（例えばconda）。原因がわかれば対処はずっと簡単です。

*3

失敗を重ねる事で勉強になる点はメリットかもしれません。

*4

正しくインストールできれば、publish済みのデータの再現性について調べることもできるようになります。反面、正しく導入できたか不明なまま再現性云々を言う人はいないでしょう。

IRLとIRRに挟まれたトランスポゾンのab initio挿入を 高感度に検出する panISa

20210910 誤字修正

2021 12/27 追記

2022 1/4 インストール手順変更

　panISaソフトウェアは、ショートリードデータから、最初に（すなわち、データベースを含まないアプローチで）NGSデータ上の挿入配列を検索する。 手短に言えば、ソフトウェアは、潜在的なISの開始位置および終了位置上のクリップされたリードを数えることによってアラインメント中の挿入サインを同定する。 これらのクリップされたリードは、IS挿入部位のダイレクトリピートとオーバーラップしている。 最後に、panISaはISの両側（IRLとIRR）の先頭を再構成して、inverted repeat領域を検索することによってISをvalidateする。

Githubより転載

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

ubuntu16.04のminiconda3.4.0.5環境でテストした。

依存

The program used the python library pysam (>=0.9) and request (>=2.12)

sudo apt update && apt install python-pysam python-requests emboss

本体　GIthub

mamba create -n panISA python=3.8 -y
conda activate panISA
#pipで導入できる
pip install panisa

 >  panISa.py -h

$ panISa.py -h

usage: panISa.py [options] bam

Search integrative element (IS) insertion on a genome using BAM alignment

positional arguments:

  bam                   Alignment on BAM/SAM format

optional arguments:

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

  -o [OUTPUT], --output [OUTPUT]

                        Return list of IS insertion by alignment,

                        default=stdout

  -q [QUALITY], --quality [QUALITY]

                        Min alignment quality value to conserve a clipped

                        read, default=20

  -m [MINIMUN], --minimun [MINIMUN]

                        Min number of clipped reads to look at IS on a

                        position, default=10

  -s [SIZE], --size [SIZE]

                        Maximun size of direct repeat region, default=20pb

  -p [PERCENTAGE], --percentage [PERCENTAGE]

                        Minimum percentage of same base to create consensus,

                        default=0.8

  -v, --version         show program's version number and exit

テストラン

bamを指定する（bwa推奨）。感度は関連するリードのリードデプスでのみ調整可能になっている。感度が高すぎるなら”-m <INT>”でリード数の閾値を増やす。

git clone https://github.com/bvalot/panISa.git
cd panISa/
panISa.py test/test.bam > output

• -m    Min number of clipped reads to look at IS on a position, default=10

出力

Chromosome End position End clipped reads Direct repeats Start position Start clipped reads Inverted repeats Left sequence Right sequence

contig-2000003 334451 59 GTCCTGGAGC 334442 40 No IR CAATGTCATCAACTTTGGAAATTATCCATAAATATCATATAATTAGCGCTCAAATCAGTGCATGGGAGNNGNC NNNNNGGCCATGGCGGCTGGCTGCTTCGGGGGGCTTGCCTTGGGCAGGGCTGCAGCTTAGGTTGATGACATTG

ランタイムは非常に短い。細菌のbamだと10秒程度で結果が得られる。深く読まれたデータなら、不完全（=集団の全体には起こっていないISイベント）なISの挿入も検出できる。

2021 12/27

結果をISFinderデータベースに問い合わせる。ネットに繋がっている必要がある。

ISFinder_search.py panISa_results　 > results

ISFinderのISとアサインされたクエリには右端にそのIS名が付く。

引用

GitHub - bvalot/panISa: panISa is a software to search insertion sequence (IS) on resequencing data (bam file)

panISa: ab initio detection of insertion sequences in bacterial genomes from short read sequence data

Panisa Treepong, Christophe Guyeux, Alexandre Meunier , Charlotte Couchoud , Didier Hocquet, Benoit Valot

Bioinformatics. 2018 Nov 15;34(22):3795-3800

2021 Deciphering the role of insertion sequences in the evolution of bacterial epidemic pathogens with panISa software

Charlotte Couchoud , Xavier Bertrand , Benoit Valot , Didier Hocquet

Microb Genom. 2020 Jun;6(6):e000356. 

関連

ゲノム中のISをスキャンする（アノテーションをつける）。ISの配列も出力される。

ナノポアのロングリードの長さやクオリティを分析する nanoQC

この論文ではOxford Nanopore TechnologiesとPacific Biosciencesのロングリードシーケンスデータの可視化と処理のために開発されたツールセット、NanoPackについて説明する。NanoPackツールはPython 3で書かれており、GNU GPL3.0ライセンスの下でリリースされている。 ソースコードはhttps://github.com/wdecoster/nanopackにあり、別々のスクリプトへのリンクとそのドキュメントもある。スクリプトは、Linux、Mac OS、およびLinux用のMS Windows 10サブシステムと互換性があり、グラフィカルユーザーインターフェース、http://nanoplot.bioinf.beのWebサービス、およびコマンドラインツールとして利用できる。

Base callingのレポートファイルを持っているなら、オンラインでも利用できる。

http://nanoplot.bioinf.be

リンク先のウィンドウ内にAlbacore/Guppyが出力するsequencing_summary.txt.ファイルをドラッグ&ドロップする（上限100MB）。

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

mac os10.14のanaconda3-5.1.0環境でテストした。

本体　GIthub

#anaconda環境ならcondaで導入できる
conda install -y -c bioconda nanoQC

#またはpipを使う
pip install nanoQC

 > nanoQC -h

$ nanoQC -h

usage: nanoQC [-h] [-v] [-o OUTDIR] [-l MINLEN] fastq

Investigate nucleotide composition and base quality.

positional arguments:

  fastq                 Reads data in fastq.gz format.

optional arguments:

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

  -v, --version         Print version and exit.

  -o OUTDIR, --outdir OUTDIR

                        Specify directory in which output has to be created.

  -l MINLEN, --minlen MINLEN

                        Filters the reads on a minimal length of the given

                        range. Also plots the given length/2 of the begin and

                        end of the reads.

実行方法

nanoQC input.fastq -o OUTDIR

出力

結果はBokehを使いインタラクティブなグラフとして可視化される。


引用

NanoPack: visualizing and processing long-read sequencing data
Wouter De Coster, Svenn D’Hert, Darrin T Schultz, Marc Cruts, Christine Van Broeckhoven
Bioinformatics, Volume 34, Issue 15, 1 August 2018, Pages 2666–2669

関連

ゲノムを比較してstructural rearrangementsを検出する SyRI

2019 12/17 論文引用追加 

　同じ種の半数体ゲノムは、典型的にはそれらのゲノム構造において高い類似性を示す広範囲のco-linear（シンテニー）領域を含む。しかし、これらのシンテニー領域は
異なるハプロタイプにおける異なる方向および／または位置によって特徴付けられるstructural rearrangements（SR）によって中断される。 SRは、逆位、転座、重複に分類することができる。（一部略）

　既存のＳＲ予測方法の多くは、リファレンス配列へのショートリードまたはロングリードのアライメントを利用する。SNPやsmall indelなどの局所的な違いは高精度で検出できるが、リードアラインメントだけでは複雑なSRを正確に予測することは困難である。対照的に、高品質のゲノムアセンブリによる比較は、通常、raw シーケンシングリードと比較してはるかに長く高品質であるため、正確なSR同定にとってより強力である［ref.6］。しかしながら、近年のde-novo全ゲノムアセンブリ生成を支える重要な技術的改良にもかかわらず[ref.7]、全ゲノムアラインメント（WGA）をゲノム差異の同定の基礎として使用するツールはほんの少ししかない [ref.8, pubmed]。たとえば、利用可能なツールには、de novoアセンブリの個々のscaffoldsをリファレンス配列と比較してアライメントのブレイクポイントを解析して逆位や転座を特定するAsmVar [ref.9]や、リファレンス配列にユニークにアライメントされたコンティグを利用してLarge indelや局所的なリピートの違いなど、さまざまなゲノムの違いを識別するAssemblyticsがある[ref.10]（紹介）。

　ここでは、ペアワイズWGAから生成されたgenome graphsを使用して、2つの関連ゲノム（通常は同一種由来）間のゲノム構造を同定するSyRI（Synteny and Rearrangement Identifier）を紹介する。 SyRIは、２つのゲノムの相同染色体間の全てのシンテニー領域を同定することから始まる。他のすべての領域は定義によりSRである（あるいはそうでなければそれらはシンテニック領域の一部である）ので、これはSR識別の問題をSR分類に変換する。 SyRIは、非シンテニー領域を逆位、転座および重複に分類する。 SyRIは、ゲノム全体にわたってリアレンジメントされた領域の分析を行い、ゲノム差異を最適化するためにSRにグローバルにアノテーションを付ける。転座と転移を区別するのが一般的だが、ここでは両方のタイプを転座と呼ぶ。さらに、転座と重複はまとめてTDと呼ぶ。最後に、SyRIは、リアレンジメント領域および非リアレンジメント領域を含むゲノム全体にわたる局所的変異を同定する。ローカルな変動および構造的なリアレンジメントは、それらのサイズまたは複雑さによって区別されないことに留意することが重要である。ローカルな変動には、large deletionまたはlarge insertionなどの大きな構造的変動も含まれ得るからである。その代わりにローカルな変動は、シンテニック内および構造的にリアレンジメントされた領域内に見出すことができる。これにより、変動が導入され、例えば、リアレンジメント領域のSNPと比較して、シンテニー領域のSNPを区別することが可能になる。リアレンジメントされた領域（およびその中のローカルな変動）は、それぞれの生物の子孫におけるメンデルの分離パターンに従わないであろうが、それらのコピー数の変化にさえつながる可能性があるので、この区別は重要である。著者らはSyRIを用いて、5つのモデル種の多様なゲノムを分析し、2つのA. thaliana株に見られる転座を遺伝的に検証し、50のF2組換えゲノムのIllumina全ゲノムシーケンシングデータを分析した。

Documentation

Synteny and Rearrangement Identifier (SyRI) | syri

インストール

ubuntu16.04のminiconda3.4.0.5環境でテストした。

依存

Python3

• Cython
• numpy
• scipy
• pandas
• python-igraph
• biopython
• psutil
• MUMmer3

conda install -y cython numpy scipy pandas biopython psutil
conda install -y -c conda-forge python-igraph
conda install -c bioconda mummer

本体　Github

git clone https://github.com/schneebergerlab/syri.git
cd syri/
python3 setup.py install
cd syri/bin/

> python syri

$ python syri 

usage: syri [-h] -c INFILE [-r REF] [-q QRY] [-d DELTA] [-o FOUT] [-k]

            [--log {DEBUG,INFO,WARN}] [--lf LOG_FIN] [--dir DIR]

            [--prefix PREFIX] [--seed SEED] [--nc NCORES] [--novcf] [--nosr]

            [-b BRUTERUNTIME] [--unic TRANSUNICOUNT] [--unip TRANSUNIPERCENT]

            [--inc INCREASEBY] [--no-chrmatch] [--nosv] [--nosnp] [--all]

            [--allow-offset OFFSET] [--cigar] [-s SSPATH]

syri: error: the following arguments are required: -c

kazu@edb2e2639563:~/syri/syri/bin$ ./syri -h

usage: syri [-h] -c INFILE [-r REF] [-q QRY] [-d DELTA] [-o FOUT] [-k]

            [--log {DEBUG,INFO,WARN}] [--lf LOG_FIN] [--dir DIR]

            [--prefix PREFIX] [--seed SEED] [--nc NCORES] [--novcf] [--nosr]

            [-b BRUTERUNTIME] [--unic TRANSUNICOUNT] [--unip TRANSUNIPERCENT]

            [--inc INCREASEBY] [--no-chrmatch] [--nosv] [--nosnp] [--all]

            [--allow-offset OFFSET] [--cigar] [-s SSPATH]

Input Files:

  -c INFILE             File containing alignment coordinates in a tsv format

                        (default: None)

  -r REF                Genome A (which is considered as reference for the

                        alignments). Required for local variation (large

                        indels, CNVs) identification. (default: None)

  -q QRY                Genome B (which is considered as query for the

                        alignments). Required for local variation (large

                        indels, CNVs) identification. (default: None)

  -d DELTA              .delta file from mummer. Required for short variation

                        (SNPs/indels) identification when CIGAR string is not

                        available (default: None)

optional arguments:

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

  -o FOUT               Output file name (default: syri)

  -k                    Keep internediate output files (default: False)

  --log {DEBUG,INFO,WARN}

                        log level (default: INFO)

  --lf LOG_FIN          Name of log file (default: syri.log)

  --dir DIR             path to working directory (if not current directory)

                        (default: None)

  --prefix PREFIX       Prefix to add before the output file Names (default: )

  --seed SEED           seed for generating random numbers (default: 1)

  --nc NCORES           number of cores to use in parallel (max is number of

                        chromosomes) (default: 1)

  --novcf               Do not combine all files into one output file

                        (default: False)

SR identification:

  --nosr                Set to skip structural rearrangement identification

                        (default: False)

  -b BRUTERUNTIME       Cutoff to restrict brute force methods to take too

                        much time (in seconds). Smaller values would make

                        algorithm faster, but could have marginal effects on

                        accuracy. In general case, would not be required.

                        (default: 60)

  --unic TRANSUNICOUNT  Number of uniques bps for selecting translocation.

                        Smaller values would select smaller TLs better, but

                        may increase time and decrease accuracy. (default:

                        1000)

  --unip TRANSUNIPERCENT

                        Percent of unique region requried to select

                        translocation. Value should be in range (0,1]. Smaller

                        values would selection of translocation which are more

                        overlapped with other regions. (default: 0.5)

  --inc INCREASEBY      Minimum score increase required to add another

                        alignment to translocation cluster solution (default:

                        1000)

  --no-chrmatch         Do not allow SyRI to automatically match chromosome

                        ids between the two genomes if they are not equal

                        (default: False)

ShV identification:

  --nosv                Set to skip structural variation identification

                        (default: False)

  --nosnp               Set to skip SNP/Indel (within alignment)

                        identification (default: False)

  --all                 Use duplications too for variant identification

                        (default: False)

  --allow-offset OFFSET

                        BPs allowed to overlap (default: 0)

  --cigar               Find SNPs/indels using CIGAR string. Necessary for

                        alignment generated using aligners other than nucmers

                        (default: False)

  -s SSPATH             path to show-snps from mummer (default: show-snps)

> ./chroder

$ ./chroder 

usage: chroder [-h] [-n NCOUNT] [-o OUT] [-noref] coords ref qry

chroder: error: the following arguments are required: coords, ref, qry

kazu@edb2e2639563:~/syri/syri/bin$ ./chroder -h

usage: chroder [-h] [-n NCOUNT] [-o OUT] [-noref] coords ref qry

positional arguments:

  coords      Alignment coordinates in a tsv format

  ref         Assembly of genome A in multi-fasta format

  qry         Assembly of genome B in multi-fasta format

optional arguments:

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

  -n NCOUNT   number of N's to be inserted

  -o OUT      output file prefix

  -noref      Use this parameter when no assembly is at chromosome level

実行方法

 SyRIはSRを同定するためにクロモソームレベルのアセンブリ配列を必要とする。ない場合、chroderを使いpseudo-chromosome配列を作成してから実行する（説明HP）。

#chromosomeレベルアセンブリでないならまずchroderを走らせる

#step1 nucmerのラン (HP)
nucmer --maxmatch -c 500 -b 500 -l 100 refgenome query_genome
delta-filter -m -i 90 -l 100 out.delta > out_m_i90_l100.delta 
show-coords -THrd out_m_i90_l100.delta > out_m_i90_l100.coords

#step2 chroder
chroder -o output out_m_i90_l100.coords ref.fa scaffolds.fa

#step3 syri
syri -c out_m_i90_l100.coords -d out_m_i90_l100.delta -r ref.fa -q output.fa

exampleランの説明がわかりやすいです。

https://schneebergerlab.github.io/syri/pipeline.html

plot例

Plotting genomic structure using plotsr | syri

引用

SyRI: identification of syntenic and rearranged regions from whole genome assemblies

Manish Goel, Hequan Sun, Wen-Biao Jiao, Korbinian Schneeberger

bioRxiv preprint first posted online Feb. 11, 2019

追記

SyRI: finding genomic rearrangements and local sequence differences from whole-genome assemblies

Manish Goel, Hequan Sun, Wen-Biao Jiao & Korbinian Schneeberger
Genome Biology volume 20, Article number: 277 (2019)

関連

効率的にペアエンドfastqを同期する Fastq-pair

2019 2/26 テストラン追加

2019 7/10 コメント追加

Fastqフォーマットのファイルは、シーケンスと品質の両方の情報を1つのファイルにまとめて含むため、DNAシーケンスを共有するための主要なファイルフォーマットとなっている（ref.1）。さらに、オーバーラップするペアを結合することによって獲得することができる潜在的により長い配列から得られるさらなる情報のために、ペアエンドシーケンシングがシーケンシングアプローチで支配的になった。たとえば、シーケンスリードアーカイブには、シングルリードライブラリの2倍のペアエンドライブラリが含まれている（2018年2月1日現在、SRAには2,233,015のペアエンドライブラリと1,110,884のシングルリードライブラリがあった）。
　ペアエンドリードデータの結合（pear（ref.2）など）、DNA配列のアセンブリ（spadeやmeta-spades（ref.3）など）、リードのリファレンス配列へのマッピング（例えばbowtie2（ref.4））に使用される多くの下流ツール））は、ペアになったシーケンスが同期していることを要求する。特に、これらのツールでは、1つのペアエンドシーケンスランの2ファイルが（i）各ファイルで同数のリードを持ち、（ii）左右（好みの用語に応じて前後）で実行される必要がある。シーケンスは各ファイルで同じ順序で現れる。対照的に、いくつかのアップストリームアプリケーションは同期してペアエンドシーケンスを提供しない。例えば、シーケンスリードアーカイブ（SRA）からシーケンスを検索するために広く使用されているfastq-dumpは、ファイル内のシーケンスの順序を自動的には同期せず、いくつかのトリミングプログラムは順序付けられていないペアエンドリードをもたらす（fastq-dumpコマンドのドキュメント化されていないオプション--split-3を使い、おそらくシーケンスを順番に分割できるが（ref.5））。
　シーケンステクノロジが向上するにつれて、fastqファイルのサイズは大きくなり、数十ギガバイトのファイルも珍しくない。これは、ペアエンドリードが同じ順序で行われるように、ペアリードの再同期をとるための課題をもたらす。著者らは、メモリと時間効率の良い方法で大きなfastqファイルを扱うため、fastq-pairを開発した。

　スタンダード（ref.1 link）によると、fastqファイルのシーケンス識別子は、先頭の@記号で示される識別子行の空白以外のすべての文字で構成されている。順方向および逆方向のリードは、通常、シーケンスIDの末尾の/ fまたは/ rで識別され、左および右は通常、シーケンスIDの末尾の/ 1および/ 2で識別される。他のほとんどのソフトウェアがそうであるように、我々は現在すべてのfastqファイルがシーケンスもクオリティ情報もラップしないと仮定する、従って、各シーケンスは正確に4行（@記号で始まる識別子行、DNAシーケンス行、スペーサーという +で始まりDNA配列の終わりを示す行、そして品質スコア行（ref.1）で表される。
　Fastq-pairは、2つのファイルのファイル名を指定することによってインスタンス化される。アルゴリズムは、最初のファイルのすべての識別子について、シーケンス識別子とファイル内の位置を含むオブジェクトのハッシュを作成する。これにより、シーケンス識別子、シーケンス、およびクオリティスコアを格納する簡単なソリューションと比較して、fastq-pairのメモリ要件が大幅に削減される。最初のファイルへのファイルポインタは、2番目のファイルが読み込まれている間はアクセス可能なままになる。2番目のファイルの各シーケンスについて、識別子がハッシュに存在する場合、2つのシーケンスとクオリティスコアは適切なペアのファイルに書き込まれる。識別子がハッシュに存在しない場合、シーケンスは適切なシングルトンファイルに書き込まれる。 2番目のファイルを反復処理した後、見られなかったシーケンスは適切なシングルトンファイルに書き込まれる。
（２段落省略）
　注意点が1つある。fastq-pairはファイルストリームへのランダムアクセスに大きく依存しているため、圧縮ファイルを処理できない。少なくとも1つのファイルはメモリか一時ファイルのどちらかに解凍する必要し、fastq-pairをメモリ効率良く扱えるようにする必要がある。
結論：fastq-pairは、fastqファイルが多くの下流の処理ステップで要求されるのと同じ順序でリードを含むことを保証する迅速でメモリ効率の良いアプローチを提供する。

インストール

mac os 10.14でテストした（cmake version 3.13.2）。

依存

We use cmake (version 3 or higher) for fastq_pair.

本体　Github

git clone https://github.com/linsalrob/fastq-pair.git
cd fastq-pair/
mkdir build && cd build
cmake ..
make && sudo make install

> fastq_pair

$ fastq_pair 

fastq_pair [options] [fastq file 1] [fastq file 2]

OPTIONS

-t table size (default 100003)

-p print the number of elements in each bucket in the table

-v verbose output. This is mainly for debugging

テストラン

ペアエンドfastqを指定するだけでランできる。 

fastq_pair -t 1000 test/left.fastq test/right.fastq

Writing the paired reads to test/left.fastq.paired.fq and test/right.fastq.paired.fq.

Writing the single reads to test/left.fastq.single.fq and test/right.fastq.single.fq

Left paired: 50 Right paired: 50

Left single: 200 Right single: 25

ジョブが終わると、元のfastqと同じディレクトリに４つのfastqが出力される（同期したペア2fastq、シングルトン2fastq）。 動作は非常に高速。200MBx2のペアエンドfastqを使うと、2.7秒でジョブは完了した（*1）。 

データによってはかなりリード数が減ってしまいます。注意して使って下さい。

引用

Fastq-pair: efficient synchronization of paired-end fastq files

John A. Edwards1​ ​ and Robert A. Edwards

bioRxiv preprint first posted online Feb. 19, 2019

*1

macbook pro 2015 15インチ上位モデル使用

まとめ

Splitting and pairing fastq files

https://edwards.sdsu.edu/research/sorting-and-paring-fastq-files/