macでインフォマティクス

macでインフォマティクス

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

インタラクティブなヒートマップを簡単に作成できるwebツール shinyheatmap

 

 ヒートマップソフトウェアは、一般に、静的ヒートマップソフトウェア (static heatmap software) [論文より ref.1-9]とインタラクティブヒートマップソフトウェア (interactive heatmap software) [ref.10-20]の2つのカテゴリに分類することができる。静的ヒートマップは、元のデータからカラー画像として生成されたゲノムアクティビティのスナップショット画像である。インタラクティブヒートマップはダイナミックパレットで、ヒートマップの内容をズームインまたはズームアウトして、特定の領域、クラスタ、または単一の遺伝子を調べると同時に、特定の行にマウスポインタを置くことができ、個々の細胞の内容(例えば、遺伝子名、発現レベル、および列名)に関する情報を収集できる。インタラクティブなヒートマップは、大きな遺伝子発現データセットを視覚化するために特に重要である。個々の遺伝子標識は、大型の入力データマトリックスの静的ヒートマップに見られる共通の欠点であるテキストオーバーラップによって、最終的に解読できなくなる。そのため、インタラクティブなヒートマップは、大規模な遺伝子発現データセットの全体的なlandscapeを調べる際に人気がある。同時に、ユーザーがヒートマップの特定の領域を拡大して見ることも可能になる(すなわち、さまざまな解像度レベルで) 。現在、さまざまな解像度で何百万ものデータポイントを視覚的に拡大することができる最新のライブラリが急務となっている。一般に、ビッグデータビジュアライゼーションソフトウェア[ref.22]のフロントエンドアルゴリズムとバックエンドアルゴリズムのオンザフライ計算には、インタラクティブなナビゲーションとさまざまな解像度レベルでのスムーズなスケーリングを実現する新しいソフトウェアインフラストラクチャが必要である。

 静的ヒートマップは依然として多くの研究において好ましいタイプのpublication figureであるが、インタラクティブヒートマップは、個々の数値がユーザ指定の色としてレンダリングされるデータセットの特定のセクタを強調し視覚化するために、 PCA、差分表現、遺伝子オントロジー、ネットワーク解析などの統合された統計的およびゲノム解析スイートとインタラクティブなヒートマップソフトウェアを相乗させることなどにより、ヒートマップ視覚化分野をデータ分析分野に段階的に移行している[ref.18 、23]。しかしながら、現在存在するインタラクティブなヒートマップソフトウェアは入力ファイルサイズ上限により機能的に有用性の範囲が制限される。例えばヒートマップ生成にpheatmap Rパッケージ[ref.9]を採用しているClustviz [ref.23]では、パフォーマンス上の理由から1000行を超える入力データセットは推奨されていない[ref.24]。同様に、MicroScopeでは、入力データセットの差分解析を最初に行い、統計的に有意な遺伝子のみを包含しレンダリングされる行数を縮小するよう促される[ref.18]。標準的な考え方は、静的ヒートマップがズーム可能ではないため、読みにくさなどのさまざまな要因の組み合わせによる大きなヒートマップの生成を避けることだった。大規模なインタラクティブなヒートマップでは、スーパーコンピュータレベルのメモリリソースで効率的な遅延フリーのズーミングとパンニングを実行する必要があるため、計算上不可能である(ref.25-31)。大きなヒートマップには多くの情報が含まれているため、標準的な推奨アプローチでは、入力データ行列をプリエンプティブに小さなサイズにサブセット化することが行われており[ref.32]、不明瞭な解釈がある。

 しかしながら、多くの場合、NGS主導の研究では10^4程度のデータセット(例えば、個々のエキソンを表す400,000行までのHTA 2.0アレイ[ref.33]などのトランスクリプトーム研究)が生成される。同様に、シングルセルのRNA-seq研究では、数千から数十万の細胞のデータセットが作成されることがある(ref.34,35)。現在、このような大きなデータをインタラクティブに視覚化することは、この方向の既存の努力にもかかわらず、既存の最先端の方法論でも不可能である[ref.36,37]。そのような前例のないサイズスケールでインタラクティブなヒートマップを視覚化する計算能力をアンロックすることで、研究者は高次元の数値データを細胞の色付きのグリッドとして調べることができる。ますます洗練されているインタラクティブなヒートマップソフトウェアの登場と、ビッグデータの登場、インタラクティブな探索方法に対するコミュニティの関心が高まるにつれて、大規模でインタラクティブなヒートマップの作成を妨げる計算上の制限に対処する緊急の必要性が生じている。このようなヒートマップは、全体的な遺伝子発現パターンおよび個々の遺伝子の両方のlandscapeを視覚化するために有益であろう。これらの目的を達成するために、著者らは、Webブラウザで高度にカスタマイズ可能な静的なインタラクティブなヒートマップを効率的に作成することができる、超高速で低メモリの使いやすいヒートマップソフトウェアスイートを提案する。

 

Github

 

実行方法

http://shinyheatmap.comにアクセスする。

f:id:kazumaxneo:20180918204322p:plain

 

テストデータのCSVを表示してみる。

ログの下のミドルサイズデータセットをクリックしてダウンロードする。 f:id:kazumaxneo:20180919203702p:plain

 

ダウンロードしたファイルをBrowseボタンからアップロードする。

 

Static Heatmap 

f:id:kazumaxneo:20180920092355j:plain

デフォルトでは、低発現の遺伝子が緑、中発現の遺伝子が色なし、高発現の遺伝子が赤になっている。

f:id:kazumaxneo:20180920092411j:plain

中発現を黒に変え、Choose Y Font Size:をいじってY軸のフォントサイズを大きくした。

f:id:kazumaxneo:20180920092637j:plain

Apply Clustering:からクラスタリングも実行できる。クラスタリング手法(デフォルトではcomple linkage clusterling)と距離(デフォルトではユークリッド距離)はそれぞれLinkage Algorithm:Distance Metric:から変更可能。 f:id:kazumaxneo:20180920093935j:plain

 

上のタブからInteractive Heatmapに切り替える。

f:id:kazumaxneo:20180920094044j:plain

マウスオーバーで遺伝子名が表示される。

f:id:kazumaxneo:20180920094131j:plain

ドラッグして囲むことで、特定の領域だけ拡大できる。これは特に字が潰れるような大きなデータセットで役に立つ。

f:id:kazumaxneo:20180920094256j:plain

右上のアイコンからも拡大縮小が可能。Panボタンを使えば自由に移動することもできる。

f:id:kazumaxneo:20180920094343j:plain

 

左下からクラスタリングしたデータもダウンロードできる。 

 

データが重いと多少待ち時間が発生するようです。慌てずに操作してください。 巨大なデータなら fastheatmap(Super high performance interactive heatmap software)も試して見て下さい(論文 図3(リンク)より)。

http://fastheatmap.com

 

引用

shinyheatmap: Ultra fast low memory heatmap web interface for big data genomics
Khomtchouk BB, Hennessy JR, Wahlestedt C

PLoS One. 2017 May 11;12(5):e0176334

 

メタゲノムから16Sなどのターゲットアセンブリを行う MATAM

2022/06/24 追記

 

Preprintより

 ショットガンのメタゲノムシーケンシングは、未知の微生物の多様性が未知のまま残っている、ヒトの微生物から土壌や海洋のサンプルまで、さまざまな用途で、未培養の微生物サンプルを研究する未曾有の機会を提供する。
 メタゲノム研究の主な目標は、微生物の多様性と生態学的構造を特徴づけることである。これはいくつかの系統発生マーカー遺伝子の1つに焦点を当てることによって達成されることが多い[ref.12]、[ref.28]。細菌集団のゴールデンスタンダードのマーカーは、Silva [ref.24]、RDP [ref.2]またはGreenGenes [ref.3]のような、キュレーションされたリファレンスデータベースで何百万もの配列が利用可能な16SリボソームRNA(rRNA、~1500bp)である。アンプリコンシーケンシングのような伝統的なアプローチは、マーカー配列の小部分の分析に限られる。これは、属を超えて、十分に正確な分類学的レベルで生物を同定するための強い技術的限界につながる[ref.23]。種、または株にマーカー配列をアサインするためには、1-kbあたりのエラー数がほとんどない状態で全長rRNAを回収できる必要がある。メタゲノムアセンブラは、ゲノム全体を扱うように最適化されており、非常によく似た配列を区別するのに苦労しているため、この作業には適していない[ref.27]。この点に関して、メタゲノムリードのサブセットを全長16S rRNAコンティグにアセンブリするために、最近、EMIRGE [ref.19]およびREAGO [ref.33]のようなマーカー指向のアプローチが開発され、環境サンプルのtaxonomic assignmentの精度を改善することを目指している。 EMIRGEはベイズ法を用いて16S rRNA全長配列を反復的に再構成する。 REAGOは、Infernal [ref.20]を使用してrRNAリードを同定し、suffix/prefix arrayを使用してリード間の正確なオーバーラップを検索することによってoverlap graphを構築する。しかし、そのようなツールは、低い存在種にリカバリー率で対処することも含め、いくつかの制限を示している。
 本研究では、エラー率とキメラ形成のリスクを最小限に抑えるように慎重に設計されたオーバーラップグラフの構築と利用に基づく新しい手法であるMATAMを提示する。 MATAMは、シミュレートされたデータと実際のシーケンシングデータの両方で検証された。ほぼ全長の16S rRNAを再構築することができ、シークエンシングデプスのばらつきやコミュニティの複雑さに強くなる。 

 

インストール

配布されているdockerコンテナを使ってmac os10.13でテストした。

依存

  • gcc v4.9.0 or superior, (full C++11 support, <regex> included, and partial C++14 support)
  • C++ libraries: rt, pthread, zlib
  • Samtools v1.x or superior
  • automake, make, cmake v3.1 or superior
  • Python 3
  • pip
  • numpy
  • Apache Ant
  • Java SE 7 JDK. OpenJDK is ok (openjdk-7-jdk paquet on debian)
  • bzip2
  • google sparse hash library (libsparsehash-dev paquet on debian)

Hardware requirements

  • We recommand running MATAM with at least 10Go of free RAM. You can try running MATAM with less RAM if --max_memory is set to a lower value (eg. --max_memory 4000 for 4Go).

本体 GIthub

 

#conda
mamba create -n matam python=3.9 -y
conda activate matam
mamba install -c bioconda matam -y

#またはdockerイメージを使う。
docker pull bonsaiteam/matam

> docker run --rm bonsaiteam/matam matam_assembly.py

$ docker run --rm bonsaiteam/matam matam_assembly.py

CRITICAL - the following arguments are required: -i/--input_fastx

usage: matam_assembly.py [-h] -i FASTX [-d DBPATH] [-o OUTDIR] [-v] [--cpu CPU] [--max_memory MAXMEM] [--best INT]

                         [--evalue REAL] [--score_threshold REAL] [--straight_mode] [--coverage_threshold INT]

                         [--min_identity REAL] [--min_overlap_length INT] [-N INT] [-E INT] [--seed INT]

                         [--quorum FLOAT] [-a {SGA}] [--read_correction {no,yes,auto}] [--contigs_binning]

                         [--min_scaffold_length MIN_SCAFFOLD_LENGTH] [--perform_taxonomic_assignment]

                         [--training_model {16srrna,fungallsu,fungalits_warcup,fungalits_unite}]

                         [--rdp_cutoff RDP_CUTOFF] [--keep_tmp] [--debug] [--resume_from STEP] [--filter_only]

 

MATAM assembly

 

optional arguments:

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

 

Main parameters:

  -i FASTX, --input_fastx FASTX                                          Input reads file (fasta or fastq format)

  -d DBPATH, --ref_db DBPATH                                             MATAM ref db. Default is

                                                                         $MATAM_DIR/db/SILVA_128_SSURef_NR95

  -o OUTDIR, --out_dir OUTDIR                                            Output directory.Default will be

                                                                         "matam_assembly"

  -v, --verbose                                                          Increase verbosity

 

Performance:

  --cpu CPU                                                              Max number of CPU to use. Default is 1 cpu

  --max_memory MAXMEM                                                    Maximum memory to use (in MBi). Default is

                                                                         10000 MBi

 

Read mapping:

  --best INT                                                             Get up to --best good alignments per read.

                                                                         Default is 10

  --evalue REAL                                                          Max e-value to keep an alignment for. Default

                                                                         is 1e-05

 

Alignment filtering:

  --score_threshold REAL                                                 Score threshold (real between 0 and 1). Default

                                                                         is 0.9

  --straight_mode                                                        Use straight mode filtering. Default is

                                                                         geometric mode

  --coverage_threshold INT                                               Ref coverage threshold. By default set to 0 to

                                                                         desactivate filtering

 

Overlap-graph building:

  --min_identity REAL                                                    Minimum identity of an overlap between 2 reads.

                                                                         Default is 1.0

  --min_overlap_length INT                                               Minimum length of an overlap. Default is 50

 

Graph compaction & Components identification:

  -N INT, --min_read_node INT                                            Minimum number of read to keep a node. Default

                                                                         is 1

  -E INT, --min_overlap_edge INT                                         Minimum number of overlap to keep an edge.

                                                                         Default is 1

  --seed INT                                                             Seed to initialize random generator. Default is

                                                                         picking seed from system time

 

LCA labelling:

  --quorum FLOAT                                                         Quorum for LCA computing. Has to be between

                                                                         0.51 and 1. Default is 0.51

 

Contigs assembly:

  -a {SGA}, --assembler {SGA}                                            Select the assembler to be used. Default is SGA

  --read_correction {no,yes,auto}                                        Set the assembler read correction step. 'auto'

                                                                         means assemble the components with read

                                                                         correction enabled when the components coverage

                                                                         are sufficient (20X) and assemble the other

                                                                         components without read correction. Default is

                                                                         auto

 

Scaffolding:

  --contigs_binning                                                      Experimental. Perform contigs binning during

                                                                         scaffolding.

  --min_scaffold_length MIN_SCAFFOLD_LENGTH                              Filter out small scaffolds

 

Taxonomic assignment:

  --perform_taxonomic_assignment                                         Do the taxonomic assignment

  --training_model {16srrna,fungallsu,fungalits_warcup,fungalits_unite}  The training model used for taxonomic

                                                                         assignment. Default is 16srrna

  --rdp_cutoff RDP_CUTOFF                                                Sequences assigned (by RDP) with a confidence

                                                                         score < 0.8 (at genus level) will be tagged as

                                                                         an artificial "unclassified" taxon

 

Advanced parameters:

  --keep_tmp                                                             Do not remove tmp files

  --debug                                                                Output debug infos

  --resume_from STEP                                                     Try to resume from given step. Steps are:

                                                                         reads_mapping, alignments_filtering,

                                                                         overlap_graph_building, graph_compaction,

                                                                         contigs_assembly, scaffolding,

                                                                         abundance_calculation, taxonomic_assignment

  --filter_only                                                          Perform only the first step of MATAM (i.e Reads

                                                                         mapping against ref db with sortmerna to filter

                                                                         the reads). Relevant options for this step

                                                                         correspond to the "Read mapping" section.

 

ラン

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

#カレントディレクトリと/homeを共有してスタート
docker run -itv $PWD:/home bonsaiteam/matam

mkdir /home/datadir
python /matam/index_default_ssu_rrna_db.py -d /home/datadir/ --max_memory 10000

 

 2、SSU rRNAのリカバリー。データベースと解析するFASTAを指定する。

matam_db_preprocessing.py -i ref_db.fasta -d /home/datadir \
--cpu 4 --max_memory 10000 -v

分類はRDP Classifier(link)を元に行われている。 fastqには対応していない。

  

引用

MATAM: reconstruction of phylogenetic marker genes from short sequencing reads in metagenomes
Pericard P, Dufresne Y, Couderc L, Blanquart S, Touzet H

Bioinformatics. 2018 Feb 15;34(4):585-591

 

samのフィルタリングツール SAMsift

 

 

SAMsiftはKarel BřindaさんがGithubで公開されている、samを様々な条件でフィルタリングできるツール。

 

インストール

mac os10.13のPython 3.6.2 :: Anaconda 3-5.0.0 でテストした。

本体 GIthub

#Anaconda環境
conda install -c bioconda samsift

#Anaconda使ってないならpipで導入
pip install samsift

> samsift -h

$ samsift -h

 

Program: samsift (advanced filtering and tagging of SAM/BAM alignments using Python expressions)

Version: 0.2.5

Author:  Karel Brinda <kbrinda@hsph.harvard.edu>

 

Usage:   samsift [-i FILE] [-o FILE] [-f PY_EXPR] [-c PY_CODE] [-m STR] 

                    [-0 PY_CODE] [-d PY_EXPR] [-t PY_EXPR]

 

Basic options:

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

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

  -i FILE               input SAM/BAM file [-]

  -o FILE               output SAM/BAM file [-]

  -f PY_EXPR            filtering expression [True]

  -c PY_CODE            code to be executed (e.g., assigning new tags) [None]

  -m STR                mode: strict (stop on first error)

                              nonstop-keep (keep alignments causing errors)

                              nonstop-remove (remove alignments causing errors) [strict]

 

Advanced options:

  -0 PY_CODE            initialization [None]

  -d PY_EXPR            debugging expression to print [None]

  -t PY_EXPR            debugging trigger [True]

 

 

ラン

ここではSAMsiftのGithubレポジトリ:test.sam(samsift/tests/tests/)を使っている。

> head -n 5 samsift/tests/tests/test.sam 

$ head -n 5 samsift/tests/test.sam 

@SQ SN:gi|561108321|ref|NC_018143.2| LN:4411709

r00001 0 gi|561108321|ref|NC_018143.2| 2798553 60 100M * 0 0 AAAGGAATGCATCACCAATCTCGTCGGCGAGGCTTCGTGCCTCCTCGCCCGGCGCGCGGGTCGCGCCCGGGTCACTCTCGCCGGCGAACCCGACATAGGC 2/454483235524234201022/1252422523223533243212222421513222236220215432131222242315235121223130142647 NM:i:2 MD:Z:51C23C24 AS:i:90 XS:i:0

r00002 0 gi|561108321|ref|NC_018143.2| 445039 60 100M * 0 0 TCATCGGCAGATGCCCCCCGAGCGCAGCTTCCACGCGGCGCCGCGACGCCGTGTGCTGGTTCGTCACCGCCCGACCGACGCGGGCCCAGTGGTCGAGCTG 11/242522323122622254122142322321121236232321205072426233223223204322022022633343113212322232-462422 NM:i:2 MD:Z:17G45A36 AS:i:90 XS:i:0

r00003 0 gi|561108321|ref|NC_018143.2| 2014256 60 100M * 0 0 GGCAGCCCCGATGGGGGACCAGATCGGCCCCTGGGCCGCCTCGCTCACCGCGAGGCTCGCGAGGAGTCCGACCAGCGCGGCGCCCACGGCCACAATCACG 2220242203322212123344242421/2-0042322131152322232426231022262222244315/353/227312462213575444122223 NM:i:2 MD:Z:19G50G29 AS:i:90 XS:i:0

r00004 0 gi|561108321|ref|NC_018143.2| 3675128 60 100M * 0 0 GTTGGGCGGGGTGGAATTCATCCATTTCATTCAGTGCCCGTTGCGAATCCCCAAGCTACCCCGACGGCGACCAGAGGATGTCGATGGGGACGGCGGCGAG 28542253205043252212122151202.212332342342273312355423123122235524002520/2421/242130/623202/22221632 NM:i:0 MD:Z:100 AS:i:100 XS:i:0

 

1、アライメントクオリティ(AS)が94以上を抽出し、bam出力

filtered.bam
samsift -i tests/test.bam -f 'AS>94' -o filtered.bam 

ASはSAMの最後のOptional tags(Optional field)にある。

f:id:kazumaxneo:20180915174705p:plain

ASとMAPQの違いはSEQanswersの回答を参照(リンク)。

 

2、アンマップのリードだけ出力

samsift -i tests/test.bam -f 'FLAG & 0x04' -o filtered.bam 

f:id:kazumaxneo:20180915174309p:plain

SAMフォーマットPDFより抜粋(SAM format  v1.6 PDF)。FLAGの0x04(10進数の4)はアンマップを意味する。

 

3、"ACCAGAGGAT"を含むリードだけ抽出

samsift -i tests/test.bam -f 'SEQ.find("ACCAGAGGAT")!=-1'

 

4、AとTのみ含むリードだけ抽出(それ以外の文字が配列中に見つかれば除かれる)

samsift -i tests/test.bam -f 're.match(r"^[AT]*$", SEQ)'

pythonのre.モジュールでの正規表現に従う。

 

5、25% ランダムサンプリング(randomモジュールを使っている)

samsift -i tests/test.bam -f 'random.random()<0.25'

 

6、乱数生成器を固定して25% ランダムサンプリング

samsift -i tests/test.bam -f 'random.random()<0.25' -0 'random.seed(42)'

ペアエンドリードのそれぞれに対して同じ値を使ったりする。

 

リストに名前があるリードだけ出力

samsift -i tests/test.bam -0 'q=open("tests/qnames.txt").read().splitlines()' -f 'QNAME in q'

 

 

タグ追加。リード長としてlnタグを追加し、平均ベースクオリティとしてabタグを追加する。

samsift -i tests/test.bam -c 'ln=len(SEQ);ab=1.0*sum(QUALa)/ln'

オリジナル NM:i:2 MD:Z:51C23C24 AS:i:90 XS:i:0

 ↓

出力    NM:i:2 MD:Z:51C23C24 AS:i:90 XS:i:0 ln:i:100 ab:f:50.62

 

 

全リードをアンマップに書き換える。

samsift -i tests/test.bam -c 'a.flag|=0x4;a.reference_start=-1;a.cigarstring="";a.reference_id=-1;a.mapping_quality=0'

 

引用

GitHub - karel-brinda/samsift: SAMsift: advanced filtering and tagging of SAM/BAM alignments using Python expressions.

 

複雑なクエリ表現に対応し、BAMを様々な条件でフィルタリングできる BAMQL

 

 Binary Alignment / Map(BAM)は、リファレンスゲノムとのアラインメント後に大量のゲノムリードデータを保持するための共通フォーマットを提供している。リードには捕捉情報が追加されており、例えばFASTAやFASTQファイルには含まれていないターゲット位置や染色体などの情報、およびメイトペアに関する情報が表示されている。

 さらなる分析のために、よくリードのサブセットのみが必要とされる。 BAMファイルを操作するための標準ソフトウェアであるSAMtoolsには、BAMファイルのフィルタリングを実行し、一致するリードを収集するオプションが用意されている[ref.1]。しかしインタフェースは表現力もなくユーザフレンドリでもない: ユーザは名前でなくビットフラグを数値で指定する必要がある。さらに、選択条件はビットフラグがマッチするかどうかによるフィルタリングに制限され、最低クオリティスコアや一致するリードグループ(@RG)などのフィルタリングは行えない。これにより、クエリの表現力が制限される。たとえば、品質の悪いリードを選択することはできない。

 汎用プログラミング言語を使用してBAMファイル操作ライブラリでファイルをフィルタリングすることができるが、サブセット操作のほとんどはad hocな作業なので、手間がかかり定型的なコードがかなり必要になる。さらに、BAM APIによって必要とされるものと、言語によって行われるチェックとの間には相違がある。たとえば、indexを使用してファイルのサブセットのみを読み取ることは可能だが、indexの要求とフィルタリング式の間で正確に一致する染色体名が必要になる。不一致が発生すると、複雑な言語混在クエリのために検出が困難なデータが漏れてしまうことがある。サブセット化プロセスを簡素化しつつ、強力なクエリを実行できるようにするため、BAM Query Language(BAMQL)を開発した。BAMQLクエリは、述語と論理結合によって指定される。 述語には、BAM flags、マッピング情報、補助読リード情報、ポジションおよびシーケンスと一致するものが含まれる。 たとえば、1番染色体上にあるすべてのリードとBAMファイルを出力するには、次の操作を行う。

bamql -f input.bam -o output.bam ’chr(1) & paired?’

使用可能なすべての連結子と述語の概要を表1に示す。

f:id:kazumaxneo:20180916105914j:plain

 論文より転載。画像をクリックで論文にジャンプ。

 

マニュアル

http://artefacts.masella.name/bamql_queries.html

BAMQLに関するツイート

 

インストール

ubuntu14.04でテストした。

依存

  • LLVM 3.4 - 4.0, HTSlib, and libuuid(ビルドに必要)

本体 GIthub

#Ubuntu
sudo apt-add-repository ppa:boutroslab/ppa && sudo apt-get update && sudo apt-get install bamql

 > bamql -h

$ bamql -h

bamql [-b] [-I] [-o accepted_reads.bam] [-O rejected_reads.bam] [-v] -f input.bam {query | -q query.bamql}

Filter a BAM/SAM file based on the provided query. For details, see the man page.

-b The input file is binary (BAM) not text (SAM).

-f The input file to read.

-I Do not use the index, even if it exists.

-o The output file for reads that pass the query.

-O The output file for reads that fail the query.

-q A file containing the query, instead of providing it on the command line.

-v Print some information along the way.

 

ラン

1番染色体にマッピングされているリードを抽出 (論文表1の"chr")

bamql -I -f input.bam -o passed.bam -O failed.bam  'chr(1)'
  • -I   Do not use the index, even if it exists.
  • -o  The output file for reads that pass the query.
  • -O  The output file for reads that fail the query.
  • -f    The input file to read.
  1. chr(glob) Matches the chromosome name to which the read is mapped.

chr以外のヘッダー、例えばplasmidLとかなら'chr(plasmidL)'と記載する。

 

複数条件も、and、or、if else, elsifなど指定できる。chr1の100-10000に適切にマッピングされたペアエンドリードだけ抽出。duplicated readsは除く。

bamql -I -f input.bam -o passed.bam -O failed.bam 'proper_pair? & !duplicate? &!supplementary? & mapping_quality(0.95) & position(100,10000) & chr(1)' 
  1. ! expr Is satisfied if expr is not satisfied.
  2. expr & expr Is satisfied only if both operands are satisfied.
  3. paired? The read is paired in sequencing.

  4. duplicate? The read is either a PCR or an optical duplicate. That is, another read with the same sequence occurs at exactly the same position in the reference genome.

  5. mapping_quality(probability) Matches the read if the probability of error is less than probability. The mapping quality is approximated in the SAM format, so this will be imperfect. For clarity, setting the probability to 0 will be so stringent as to reject all reads, while setting it to 1 will be so liberal as to accept all reads.
  6. supplementary? The alignment is supplementary.

  7. position(start,end)
    Matches all sequences that cover the range of position from specified start to end, i.e. N matches any base.

 

詳細はオンラインマニュアルを読んでください。

http://artefacts.masella.name/bamql_queries.html 

引用

BAMQL: a query language for extracting reads from BAM files

Masella AP, Lalansingh CM, Sivasundaram P, Fraser M, Bristow RG, Boutros PC

BMC Bioinformatics. 2016 Aug 11;17:305

 

fastqの配列をランダムに変化させる fastq-anonymous

 

インストール

mac os10.13のPython 3.6.2 :: Anaconda 3-5.0.0 でテストした。

本体 GIthub

pip install fastq-anonymous


#Anaconda環境なら
conda install -c bioconda fastq-anonymous

fastq-anonymous -h

$ fastq-anonymous -h

usage: fastq-anonymous [-h] [-v] [-m]

 

Change the sequence of a fastq file to enable sharing of confidential

information, for troubleshooting of tools.

 

optional arguments:

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

  -v, --version  Print version and exit.

  -m, --mask     Mask all nucleotides using N

 

ラン

fastqを指定する。

cat reads.fq | fastq-anonymous > anonymous_reads.fq
  •  -m   Mask all nucleotides using N

ヘッダー名が匿名になり、配列もランダムに置き換えられる。クオリティ、リード長、リード数は維持される。

$ head anonymous_reads.fastq 

@dummy0

GACCAGCACCTCGACCGGTCGCTGTGCTCAGAGACTTTATGCGCAAGTTAGCCGTCTTCCCTCCTGTTTGGTTACGTCCGCAGGTCGGGATCACATTCAGATGGATCGGAGCCCGCCACGATGTGAATTAAGCCCGATTTCGCAGGTCGACATAAGTGTTGTTATCTGCAACCCGCTGTAATTGTGACGTAGGTCTCTCCCCAACCATGCACCACGAAATTGTAACAACGTCGGATCATACGCAACTTGT

+

#9GGEF=CG*GFG?G#=G1#GGGGGG#G#GGEG;#GGGFGGGCG?GCGGGEEGGGFGDGGGG8GFGGGFGGGGGGFFGGGCGG*GFGG@FAGGGGEFGGGFGGFGGGG?GGGGG?EGGFGGG+3GGGFGGGCFGGEGGGGE,FGCGGGDGGGGGCGGCGGGGGGGCGGFGGGGGGEGGCGFGGGGGGG,GGGGGGGF:GGGFGGGGFGGGGGGGGGEGGGFGGGFGCFGGGGGGGGFGGGGGGGGCCCCC

 

@dummy1

GATTTTAGTCGCCACGAGTTCCTTCCTCATGATATTTGGTCAGCTTTTCTGATATGCAACTGCCACTGCAACTCATAGCCAGCATCAGCCTCGCCACGGCCCCACATGTAGCACACTGTGCAGCGTGCTATTGGACTACGAGTCCCCCTCGAGTTGGACACCCCACGAGTTCAGGAGGCCAGGGCTACTATCCCAAACGGTGGTCGGATCCTGCACGTGGCCGGAAACAGAAATCAAACAATTTGCACAC

+

G27GGEGG<GGGGGGA+GGGG#GE2C##G*GAG+GFGGGGGGGFC>GGCGGGGGCGGGGG*GGGGCGBFGGGGGG1GGFGFGGGGGGGFGGGGGGGGGG<GFCFGGF,GGGGGGCGGGGGGGGGGGGGGGGFGGBFGGGGGCFB<=GGGGGFGGGGGGGG=GGGGGFFGGFGGGGGCGGGFGGGFGGGGGGGGCGGCGCGGGFABCGG@GEGGCGGGGGGGG,GGGGDGGGCGGGGCGGCGGCGGC@CCC

 

gzip圧縮のfastqを使う場合、解凍して読み込む。gzip圧縮して出力する。

gunzip -c reads.fastq.gz | fastq-anonymous | gzip > anonymous_reads.fastq.gz

 

引用

GitHub - wdecoster/fastq-anonymous: Change the sequence of a fastq file to enable sharing of confidential information for troubleshooting

 

細菌性髄膜の原因バクテリア種をFASTA配列から検出するBMScan

 

 毎年世界で120万件が発生すると推定される細菌性髄膜炎(bacterial meningitis)は、公衆衛生上の懸念事項として残る、生命を脅かす感染症である [論文より ref.1]。数多くの病原体が細菌性髄膜炎を引き起こす可能性があり、病原体あたりの致命率や病気の有病率は地域、国、年齢によって異なる[ref.2]。細菌性髄膜炎は、一般に髄膜炎菌、肺炎連鎖球菌、リステリア菌、ヘモフィルスインフルエンザ菌および大腸菌によって引き起こされる[ref.3]。髄膜炎を引き起こすバクテリア種の同定は、症例管理および疾病監視のための重要なステップである。

 種を同定する古くからの標準はDNA-DNAハイブリダイゼーション(DDH)であり、これは2つの生物の間の距離を計算するためにDNAプール間の配列類似性に依存する技術である。 DDHでは、同種比較の伝統的なカットオフ値は70%と決定された[ref.4 link]。この方法の複雑さのために、細菌の表現型の特徴をターゲットにした種の同定方法が研究室の方法として開発されている。特定の種を確認するためには、しばしば複数の表現型を捉える方法が必要とされる。

 ゲノムデータの生成が容易になるにつれて、全ゲノム配列決定(WGS)ベースのツールが開発されており、代表的なゲノムのコレクションとの比較で種同定を行う動きが加速している。ゲノムデータが手に入ることで、それを下流の分析に利用できるという追加の利点も出てくる。これらのWGSベースの比較ツールの1つは、2つのゲノム間の相同なヌクレオチド断片を比較することによってゲノム類似性を評価するAverage Nucleotide Identity (ANI)である[ref.5]。 ANIは、原核生物種同定のためのゴールデンスタンダードな全ゲノムメソッドであると考えられている[ref.6]。 95%ANIは70%DDHに匹敵すると報告されている[ref.7]。 ANIの2つの一般的な実装は、それぞれBLASTアルゴリズム[ref.8]とMUMmerメソッド[ref.9]を使用するANI BLAST(ANIb)とANI MUMmer(ANIm)である。これらのANI法はゲノム間の遺伝的類似性を評価するための高レベルの分解能を提供するが、トレードオフは長い計算時間であり、大きなリファレンスコレクションに対して複数のゲノムをスキャンすることは不可能である。

 この限界に対処するために、遺伝的距離を推定するためのk-merベースの比較を用いて実行時間を改善することに焦点を当てたツールが開発された。これらのツールの主な例は、2つのゲノム間の距離を推定するMinHashアルゴリズムを適用するMash [ref.10]である(紹介)。彼らの論文でOndov et al. は、Mashによって推定された遺伝的距離は、およそ1-ANIと強く相関することを示した。すなわち、0.05のMash distance が0.95のANIに対応する。

 著者らは、Mashが種の描写のためにANIによって提供されたものと同等であったかどうかを調べた。次に、Mashを用いて様々な種について種特異的な閾値を確立した。最後に、著者らは、著者らの焦点となっている種、ならびにいくつかのclosely related speciesおよびsister-speciesからなる新しいゲノムコレクションを用いて、これらの閾値の正確さを検証した。これらの結果を用いて、BMScanという、ゲノム全体の類似性が系統的な種特異的な閾値を上回っている場合に、characterizeされていないクローンを著者らのコレクションに迅速かつ正確にアサインするプログラムを開発した。

 

f:id:kazumaxneo:20180915112202p:plain

BMScan Workflow: The figure above shows the workflow for a query in BMScan. The input query is an assembly file in FASTA format. 論文より転載。

 

 

BMScanに関するツイート

 

BMScanのターゲットとなるバクテリア

#### 細菌性髄膜炎を起こすバクテリア
N. meningitidis
H. influenzae
S.pneumoniae
L.monocytogenes
E.coli

#### Other Neisseria sp. of interest
N. polysaccharea
N. bergeri
N. weaveri
N. subflava
N. mucosa
N. lactamica
N. gonorrhoeae
N. elongata
N. cinerea

#### Other Haemophilus sp. of interest
H. parainfluenzae
H. parahaemolyticus
H. haemolyticus

 

インストール

本体 bitbucket

git clone https://bitbucket.org/ntopaz/bmscan.git #700MBほどあるので注意
cd bmscan/bin/

python identify_species.py -h

$ python identify_species.py -h

usage: identify_species.py [-h] [-d INDIR] [-f FILE] [-v] -o OUT [-j]

                           [-t THREADS]

 

Script for quickly determining species

 

optional arguments:

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

  -d INDIR, --indir INDIR

                        Input Directory: Directory of FASTA files to analyze

  -f FILE, --file FILE  Tab-delimited input containing Name and File Path

                        (Name Filepath)

  -v, --verbose         verbose standard output (default = false)

  -o OUT, --out OUT     Output File name

  -j, --json            Only output json file (default = false)

  -t THREADS, --threads THREADS

                        Number of max threads to use (default=1)

 

ラン

ランするには、シーケンシングデータをアセンブリして得たFASTAを収納したディレクトリを指定する。

identify_species.py -d assenbly_dir -j -t 4 -o outpit_dir

 ジョブが終わると、指定ディレクトリに.jsonと.csvが出力される。

 

JSONビューアjqで整形表示する。

brew install jq #jqをインストール
cat output/species_analysis_1536983455.385941.json  | jq

$ cat species_analysis_1536983455.385941.json  | jq

{

  "Neisseria_meningitidis_38277.fasta": {

    "mash_results": {

      "notes": "Hit above threshold",

      "species": "Neisseria meningitidis",

      "top_hit": "M36993_HUY3701A58_cleaned.fasta",

      "score": 0.99659813,

      "source": "BML_collection",

      "mash_pval": "0",

      "mash_hash": "871/1000"

    }

  }

}

 

引用

BMScan: using whole genome similarity to rapidly and accurately identify bacterial meningitis causing species
Topaz N, Boxrud D, Retchless AC, Nichols M, Chang HY, Hu F, Wang X

BMC Infectious Diseases 201818:405

 

MinHashを使った高速なANI計算ツール fastANI

2019 1/09 cocndaインストール追記 ,2/12 不要な文を削除, 4/12 dockerリンク追加

2020 4/2 インストール手順修正2022 03/28 help更新

 

 さまざまな生態学的背景と進化の歴史を持つ原核生物ゲノムのコレクションが公開されている。このゲノムデータの大洪水は、微生物生態学と進化における重要な問題をより堅固に評価する機会を提供し、大きなゲノムデータの分析のための既存のバイオインフォマティクスアプローチを進歩させる必要性を強調している。そのような疑問の1つは、バクテリア(および他の微生物)が離散したクラスター(discrete clusters )(species)を形成するか、または水平遺伝子伝達(HGT)の頻度が高く、ゆっくりと減衰するキネティクスのために、連続的な遺伝的多様性が観察されるかどうかである。Closely relatedな少数のゲノムに基づく研究は、genetic continuum(連続性)が優先される可能性があることを示している(例えば、ref.1)。他方、HGTは種間の境界を歪めるほど頻繁ではないかもしれない、あるいは同じ種内の生物が種間の生物に比べてより頻繁に遺伝子交換を行い、それによって異なるクラスターを維持するとの研究がある(例えばref.2)。これら全ての研究に対する重要な批判は、それらが典型的に実験室の単離されたゲノムを用いて行なっているため、培養バイアスによって自然の多様性を適切に表していないかもしれないというものである。または、培養可能少数の系統由来の利用可能なゲノムしか使っていないため、堅牢な結論を出すことはできないという批判もある。(一文略)。種の定義は重要なアカデミックの練習であるだけでなく、大きな実践的な結果をもたらす。例えば、病原体の診断、国境を越えて生物を輸送することができる規制、隔離されるべき病原体、または人間や動物、植物に有益な有機物や有機物の混合物に関するコミュニケーションなど、すべて種がどのように定義されているかに根ざしている。

 種の境界を評価する際の1つの基本的な作業は、2つのゲノム間の遺伝的関連性の推定である。The whole-genome average nucleotide identity(ANI)は、同じ種に属する生物が典型的にはそれらの間に≥95%のANIを示すため(ref.3,4)、近年この作業のための堅牢な方法として浮上している。 ANIは、任意の2つのゲノム間で共有されるすべてのオーソロガス遺伝子の平均ヌクレオチド同一性を表し、同一またはclosely relatedな種の株間比較で強い分解能(すなわち、80-100%のANIを示す)を示す。オルソログ遺伝子はゲノム間で大きく異なる可能性があるので、ANI測定は厳密にコアゲノムの進化的関連性を表すものではない。それにもかかわらず、それはバクテリア遺伝子プールの液体的性質を考慮し、したがって暗黙に共有機能を考慮しているので、種(ref.3)を定義するための伝統的なDNA-DNAハイブリダイゼーションをよく反映している。種を定義し、それらの進化的一意性を評価するための16S rRNA遺伝子のシーケンシングと比較して、ANIは、closely relatedなゲノム間でも高分解能を有するなどのいくつかの重要な利点を持っている。最後に、16S rRNA遺伝子などの普遍的に保存された遺伝子をコードしない(e.g., due to mis-assembly)メタゲノムまたはシングルセル技術を使用して、環境から回収されたドラフト(不完全な)ゲノム配列中のANIを推定することができ、ユニバーサル遺伝子ベースのアプローチと比較して少なくとも数百の共有遺伝子を使うことで研究および分類できる配列の数が大幅に拡大する。したがって、ANIは、DNA-DNAハイブリダイゼーションを、ゲノムの関連性の標準的な尺度として置き換える可能性について、国際的に認められている。(一文略)その強みにもかかわらず、伝統的なANI計算はアライメントベースの検索[例えばBLAST(ref.7)]に基づいており、アライメントアルゴリズムの二次的な時間複雑さのために計算上高価なままである。

 ANI計算アルゴリズムはいくつかのバリエーションが提案されているが(ref.8-10)、これらは主に共有遺伝子を同定するための特定のアプローチを変更し、それらがすべてアライメントベースであるため実質的に計算を高速化しない。したがって、これらの手法および一般的に利用可能な計算資源に基づき、現在までに入手可能な数十万のオーダーの微生物ゲノムでANI値を計算することは、ほとんど不可能である。重要なことに、利用可能なゲノムデータは現存する原核生物多様性のわずかな部分であると推定され(ref.11)、決定された新しいゲノムの数は指数関数的に増加し続ける。したがって、利用可能なデータと将来のデータをスケールアップするためには、新しい計算ソリューションが必要である。
 このような解決策のいくつかは、他の科学分野での「ビッグデータ」分析から概念を借りて最近提案されている。 MinHashは、World Wide Web(ref.12)で検索エンジンでほぼ重複するWeb文書を検出するために開発された、2セット間の類似性を迅速に評価する手法である。最近、この技術は、ゲノムアセンブリ(ref.13,14)やロングリードマッピング問題(ref.15)など、バイオインフォマティクスにおける新しい高速アルゴリズムの設計に成功した。 Ondov et al(ref.16)は、このテクニックを使用したANIの高速推定のためのMASHと呼ばれる最初の概念実証実装を提供した(紹介)。ANI計算についてMASHはアライメントベースよりも数桁高速であると報告されているが、ANIの計算問題に対してのMinHashの直接的採用は、不完全なドラフトゲノムについては不正確であることが判明している(ref.17) 。また、中程度に異なるゲノム
について(例えば80~90%のANIを示す時)、MASHがANIとどれくらいうまく近づけることができるかには限界があり、それは、ANIが共有された領域のみ計算するのと異なり、MASHの相同性検索が共有されていない領域についても相同性検索を行うためである。
 本研究では、MinHashベースのアラインメントフリーのシーケンスマッピングエンジンとしてMashmap(15)を利用した新しいアルゴリズムFastANIを開発し、ANI計算の計算上のボトルネックを軽減する。 FastANIは、80%〜100%のヌクレオチド同一性範囲に関連する完全およびドラフトクオリティのゲノムの両方について、アライメントベースのANI値と本質的に同一のANI値を提供する。したがって、FastANIにより、大量のゲノムコホートのペアワイズANI値の正確な推定、またはクエリーのドラフトゲノムと利用可能な原核生物ゲノムの完全コレクションの比較が可能になる。

 

 

f:id:kazumaxneo:20180914125716j:plain

Fig. 5. a. Graphical illustration of FastANI’s work-flow for computing ANI between a query genome and a reference genome. 

(Preprintより転載。写真をクリックでPreprint PDFにジャンプします。)

 

図1にMASHとの比較あり。ANIが低い時のズレが問題だが、データセットD1の結果を見ると、fastANIの方がANIとよく相関している。データセットD3の結果もそうで、MASHはANIが高い時もズレがでる傾向にある。

  

 

インストール

ubuntu16.04でテストした。

依存

本体 GIthub リリースよりlinux向け実行ファイルのバイナリをダウンロードする。

#condaでも導入可能(maclinux
mamba create -n fastani -y
conda activate fastani
mamba install -c bioconda -y fastani

fastANI -h #version 1.33

# ./fastANI -h

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

fastANI is a fast alignment-free implementation for computing whole-genome Average Nucleotide Identity (ANI) between genomes

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

Example usage:

$ fastANI -q genome1.fa -r genome2.fa -o output.txt

$ fastANI -q genome1.fa --rl genome_list.txt -o output.txt

 

SYNOPSIS

--------

./fastANI [-h] [-r <value>] [--rl <value>] [-q <value>] [--ql <value>] [-k

          <value>] [-t <value>] [--fragLen <value>] [--minFraction <value>]

          [--visualize] [--matrix] [-o <value>] [-v]

 

OPTIONS

--------

-h, --help

     print this help page

 

-r, --ref <value>

     reference genome (fasta/fastq)[.gz]

 

--rl, --refList <value>

     a file containing list of reference genome files, one genome per line

 

-q, --query <value>

     query genome (fasta/fastq)[.gz]

 

--ql, --queryList <value>

     a file containing list of query genome files, one genome per line

 

-k, --kmer <value>

     kmer size <= 16 [default : 16]

 

-t, --threads <value>

     thread count for parallel execution [default : 1]

 

--fragLen <value>

     fragment length [default : 3,000]

 

--minFraction <value>

     minimum fraction of genome that must be shared for trusting ANI. If

     reference and query genome size differ, smaller one among the two is

     considered. [default : 0.2]

 

--visualize

     output mappings for visualization, can be enabled for single genome to

     single genome comparison only [disabled by default]

 

--matrix

     also output ANI values as lower triangular matrix (format inspired from

     phylip). If enabled, you should expect an output file with .matrix

     extension [disabled by default]

 

-o, --output <value>

     output file name

 

-v, --version

     show version

 

dockerイメージ

https://hub.docker.com/r/staphb/fastani

 

 ラン

1、シングルのクエリゲノムとリファレンスゲノム間の比較

fastANI -q query.fa -r ref.fa -o output

ジョブが終わると"output"ができる。

 

2、シングルのクエリゲノムと複数のリファレンスゲノム間の比較。リストファイルはゲノムのFASTAのパスを1行に1つずつ記載したものを指定する。

#リファレンスリストを作る。
cd genome/
ls *fasta > reference_list #注 実行時のカレントパスにfastaがないならフルパス指定する必要がある

fastANI -t 12 -q query.fa --rl reference_list -o output

 

3、複数クエリゲノムと複数のリファレンスゲノム間の比較

fastANI -t 12 --ql query_list --rl reference_list --matrix -o output

outputとoutput.matrixができる。Refseqのバクテリアゲノムを適当に100選んでAll vs All比較を行ったが、39秒しかかからなかった("-t 8"フラグをつけて実行)。

 

4、genoPlotRを使い2ゲノム間の比較結果を可視化する。レポジトリのRスクリプトを使う。

#genoPlotRがなければ入れておく(リンク)。R内で
> install.packages("genoPlotR", repos="http://R-Forge.R-project.org")

#ゲノム比較実行
fastANI -q query.fa -r ref.fa --visualize -o output
#outputとoutput.visualができる。できなければgenoPlotRが入っていないかもしれない。

Rscript scripts/visualize.R B_quintana.fna B_henselae.fna fastani.out.visual

出力。 

f:id:kazumaxneo:20180914123327j:plain

 

引用

High-throughput ANI Analysis of 90K Prokaryotic Genomes Reveals Clear Species Boundaries

Chirag Jain, Luis M. Rodriguez-R, Adam M. Phillippy, Konstantinos T. Konstantinidis, Srinivas Aluru

bioRxiv preprint first posted online Nov. 27, 2017;

doi: http://dx.doi.org/10.1101/225342.

 

Preprint紹介後、2018年11月にNature Communicationsにアクセプトされています。

High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries
Chirag Jain, Luis M. Rodriguez-R, Adam M. Phillippy, Konstantinos T. Konstantinidis & Srinivas Aluru
Nature Communications volume 9, Article number: 5114 (2018)

 

 

 

 Mashmap

MashmapはD-GENIESでも使われています。

fastANIを使ったクラスタリング