macでインフォマティクス

macでインフォマティクス

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

KEGGパスウェイと発現プロファイルデータを視覚的に統合するウェブサーバー KeggExp

 

発現プロファイルデータを解析する際の知識発見には、効果的な可視化が重要である。しかし、発現プロファイルデータをKEGGパスウェイマップと視覚的に統合するための既存のツールは、広範なインタラクティブな可視化操作を欠いている。KeggExpは、1つのパスウェイのマップ、パスウェイ内の遺伝子のデンドログラムとヒートマップ、1つの遺伝子の散布図を同時に表示し、また、異なる発現遺伝子、ヒートマップから選択された共発現遺伝子、ユーザー入力遺伝子を含むパスウェイマップ上の特定の遺伝子を強調表示するためのインタラクティブな操作を提供する。KeggExpを使用すると、プログラミングスキルのない研究者を含む研究者は、発現プロファイルデータを分析する際に、キーとなる遺伝子やパスウェイを決定するために、そのインタラクティブな操作を活用することができる。

 

webサービス

ExpressVisのKeggExpタブ(http://www.fgvis.com/expressvis/KeggExp)にアクセスする。

f:id:kazumaxneo:20210302085421p:plain

 

exampleデータを読み込む。Try sample1をクリック。

f:id:kazumaxneo:20210302085527p:plain

 

パスウェイ図とヒートマップ図はそれぞれ拡大縮小とスクロールができる。マウスホイールでパスウェイを拡大した。

f:id:kazumaxneo:20210302090750p:plain

 

パスウェイのアサインされたタンパク質名をクリックすると、その遺伝子の発現がハイライト表示される。

f:id:kazumaxneo:20210302090601p:plain

 


ヒートマップは階層的クラスタリングされた上て表示されている。ヒートマップ横のデンドログラムをクリックすると該当するタンパク質がハイライト表示される。

f:id:kazumaxneo:20210302092035p:plain

 

ヒートマップ図の外観は歯車マークからカスタマイズできる。

f:id:kazumaxneo:20210302092226p:plain

 

高さを変更した。

f:id:kazumaxneo:20210302092331p:plain

 

デンドログラムを高くした。

f:id:kazumaxneo:20210302092414p:plain

 

sample2は3つのマップを切り替え表示できるようになっている。切り替えするには左上の矢印マーク(<-  ->)をクリックする。

f:id:kazumaxneo:20210302092710p:plain

パスウェイ図に合わせてヒートマップ図も切り替わる。

 

ユーザーのデータを読み込むにはimport dataをクリックする。

f:id:kazumaxneo:20210302092904p:plain

 

まず種を選択。

f:id:kazumaxneo:20210302093022p:plain

現在は10 species対応している。

 

タブ区切り、またはexcel形式の発現行列ファイルをアップロードする。1列目は Entrez IDまたはEnsembl IDを記載し、2列目は遺伝子シンボルを記載する。3列目以降にそれぞれのサンプルのlog変換していない発現量値を記載する。

f:id:kazumaxneo:20210302093104p:plain

 

どのカラムがどのような情報を持っているのか指示する。

f:id:kazumaxneo:20210302093509p:plain

 

引用
KeggExp: a web server for visual integration of KEGG pathways and expression profile data
Xian Liu, Mingfei Han, Chen Zhao, Cheng Chang, Yunping Zhu, Changhui Ge, Ronghua Yin, Yiqun Zhan, Changyan Li, Miao Yu, Fuchu He, Xiaoming Yang

Bioinformatics. 2019 Apr 15;35(8):1430-1432

 

BreakID

 

標的としたシーケンシングデータから融合イベントを同定するためのBreakpoint Identification(BreakID)と呼ばれるツールを開発した。BreakIDは、不一致リードペアやスプリットリードを支持証拠とし、1ヌクレオチド分解能で遺伝子融合のブレークポイントを同定することができる。ガン細胞株で確認された融合イベントを用いて検証を行った結果、BreakIDは500×のシークエンシング深度で90.63%の高感度、100%のPPVを達成し、他の融合検出ツールよりも優れた性能を発揮することが証明された。BreakIDが、臨床や研究のシーケンシングシナリオに関わる融合イベントの検出と解析に広く普及することを期待している。ソースコードhttps://github.com/SinOncology/BreakID で自由に入手できる。

 

インストール 

ubuntu18.04LTSでテストした。リリースからビルド済みのプログラムを入手できると説明にあるが、2021年3月現在公開されていない。自分でビルドする必要がある。

Github

mkdir nib_files 
cd nib_files
rsync -avzP rsync://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/ .
gunzip *.gz
wget https://users.soe.ucsc.edu/~kent/src/blatSrc35.zip
unzip blatSrc35.zip

git clone https://github.com/SinOncology/BreakID.git
cd BreakID/
sh install.sh

ビルドに失敗する。

 

 

引用

BreakID: genomics breakpoints identification to detect gene fusion events using discordant pairs and split reads
Linfang Jin, Jinhuo Lai, Yang Zhang, Ying Fu, Shuhang Wang, Heng Dai, Bingding Huang
Bioinformatics, Volume 35, Issue 16, 15 August 2019, Pages 2859–2861

ユーザーが関心のある遺伝子の周辺配列の相同性を調べる Flanker

2021 3/1  renameコマンド修正

 

 関心のある遺伝子の周囲の配列を解析することは、特に抗菌薬耐性遺伝子などの水平遺伝子導入におけるmobile genetic elements (MGEs) の役割を理解する上で非常に重要であることが多い。ここでは、一貫した形式で遺伝子のフランキング配列のアラインメントフリーなクラスタリングを行うPythonパッケージFlankerを紹介する。FlankerはMash距離に基づいてフランキング配列をクラスタリングし、配列間の類似度や類似度の程度を簡単に比較できるようにしている。さらに、Flankerは柔軟にパラメータを設定して、上流と下流の領域を別々に特徴付けしたり、可変長のフランキング配列を調査したりして、出力を微調整することができる。重要な carbapenemase genes(blaOXA-48とblaKPC-2/3)のプラスミド関連キャリッジを記述した最近の2つのデータセットにFlankerを適用し、既知のものと以前には明らかにされていなかった構造変異を含む、異なるフランキング配列のクラスター(フランキングパターン)を識別することに成功したことを示した。その結果、フランキングパターンが地理的地域やカルバペネム表現型と関連していることを示し、それらが疫学的マーカーとして有用である可能性を示唆している。FlankerはMITライセンスのもと、https://github.com/wtmatlock/flanker で自由に利用できる。

 本原稿で行った解析は、Flanker Githubページ(https://github.com/wtmatlock/flanker)で提供されているBinder環境で再現することができる。

 

Documentation

https://flanker.readthedocs.io/en/latest/

 

 

f:id:kazumaxneo:20210228200918p:plain

Flanker. Githubより転載

 

インストール

Github

#conda + pip(ここでは高速なmambaを使う)
mamba create -n flanker -c bioconda python=3 abricate=1.0.1 mash -y
conda activate flanker
pip install git+https://github.com/wtmatlock/flanker

> flanker

$ flanker 

usage: flanker [-h] -i FASTA_FILE (-g GENE [GENE ...] | -log LIST_OF_GENES) [-cm] [-f FLANK] [-m MODE] [-circ] [-inc] [-db DATABASE] [-v [VERBOSE]] [-w WINDOW] [-wstop WINDOW_STOP] [-wstep WINDOW_STEP] [-cl] [-o OUTFILE] [-tr THRESHOLD]

               [-p THREADS] [-k KMER_LENGTH] [-s SKETCH_SIZE]

 

Flanker (version 0.1.5). If you use Flanker in your work, please cite us: Matlock W, Lipworth S, Constantinides B, Peto TEA, Walker AS, Crook D, Hopkins S, Shaw LP, Stoesser N.Flanker: a tool for comparative genomics of gene flanking

regions.BioRxiv. 2021. doi: https://doi.org/10.1101/2021.02.22.432255

 

optional arguments:

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

  -g GENE [GENE ...], --gene GENE [GENE ...]

                        Gene(s) of interest (escape any special characters). Use space seperation for multipe genes (default: None)

  -log LIST_OF_GENES, --list_of_genes LIST_OF_GENES

                        Line separated file containing genes of interest (default: False)

  -cm, --closest_match  Find closest match to query (default: False)

  -f FLANK, --flank FLANK

                        Choose which side(s) of the gene to extract (upstream/downstream/both) (default: both)

  -m MODE, --mode MODE  One of "default" - normal mode, "mm" - multi-allelic cluster, or "sm" - salami-mode (default: default)

  -circ, --circ         Is sequence circularised (default: False)

  -inc, --include_gene  Include the gene of interest (default: False)

  -db DATABASE, --database DATABASE

                        Choose Abricate database e.g. NCBI, resfinder (default: ncbi)

  -v [VERBOSE], --verbose [VERBOSE]

                        Increase verbosity: 0 = only warnings, 1 = info, 2 = debug. No number means info. Default is no verbosity. (default: 0)

 

required arguments:

  -i FASTA_FILE, --fasta_file FASTA_FILE

                        Input fasta file (default: None)

 

window options:

  -w WINDOW, --window WINDOW

                        Length of flanking sequence/first window length (default: 1000)

  -wstop WINDOW_STOP, --window_stop WINDOW_STOP

                        Final window length (default: None)

  -wstep WINDOW_STEP, --window_step WINDOW_STEP

                        Step in window sequence (default: None)

 

clustering options:

  -cl, --cluster        Turn on clustering mode? (default: False)

  -o OUTFILE, --outfile OUTFILE

                        Prefix for the clustering file (default: out)

  -tr THRESHOLD, --threshold THRESHOLD

                        mash distance threshold for clustering (default: 0.001)

  -p THREADS, --threads THREADS

                        threads for mash to use (default: 1)

  -k KMER_LENGTH, --kmer_length KMER_LENGTH

                        kmer length for Mash (default: 21)

  -s SKETCH_SIZE, --sketch_size SKETCH_SIZE

                        sketch size for mash (default: 1000)

 

 

実行方法

1、使用する配列はあらかじめ連結しておく。

cat *fsa > seq.fasta

 

任意のステップ

FASTA ヘッダの名前を変更して、元のファイルと一致するようにする (まだ一致していない場合)。これを行うスクリプトが提供されている。実行すると、fastaのヘッダがファイル名と同じになる。multi-fastaは想定されていない。1つのファイルに2つ以上の配列を含むファイルを使わないこと。

ls *fsa | sed 's/[.]fsa//' > input_files
git clone https://github.com/wtmatlock/flanker.git
python flanker/scripts/multi_fa_rename.py seq.fasta input_files david_plasmids_renamed.fasta
mv david_plasmids_renamed.fasta david_plasmids.fasta

 

2、 周囲の配列の関心のある遺伝子と検出サイズ、そして検出の向き(上流か下流か両方か)を指定してflankerを実行する。例では、blaKPC-2遺伝子の上流 (-f upstream) 方向について、 100bp単位 (-wstep) で0bp (-w) から 5000bp (-wstop) まで抽出して比較している。また、blaKPC-2遺伝子自身も含んでいる(--include_gene)。

flanker -f upstream -w 0 -wstop 5000 -wstep 100 -g blaKPC-2 -i seq.fasta --include_gene
  •  -f   Choose which side(s) of the gene to extract (upstream/downstream/both) (default: both)
  • -w   Length of flanking sequence/first window length (default: 1000)
  • -wstop   Final window length (default: None)
  • -wstep   Step in window sequence (default: None)
  • -g    Gene(s) of interest (escape any special characters). Use space seperation for multipe genes (default: None)
  • -i    Input fasta file (default: None)
  • --include_gene   Include the gene of interest (default: False)
  • -log    Line separated file containing genes of interest (default: False)

出力例

f:id:kazumaxneo:20210301122517p:plain

 

Documentではより発展的な使い方についても解説されています。確認してください。

https://flanker.readthedocs.io/en/latest/

 

引用

Flanker: a tool for comparative genomics of gene flanking regions

William Matlock, Samuel Lipworth, Bede Constantinides, Timothy E.A. Peto, A. Sarah Walker, Derrick Crook, Susan Hopkins, Liam P. Shaw, Nicole Stoesser

bioRxiv, Posted February 22, 2021

 

関連


原核生物の遺伝子配列をGenBankファイルから簡単に検索するためのウェブツール BAGET 2.0

 

 細菌ゲノムや古細菌ゲノムの完全シークエンスから単一の遺伝子配列とコンテキストを検索することは、ウェットベンチ生物学者にとっては気が遠くなるような作業である。既存のウェブベースのゲノムブラウザは、日常的に使用するには複雑すぎるか、原核生物ゲノムのサブセットしか提供していない。
 菌類や古細菌の遺伝子の塩基配列シンテニーをマウスを3回クリックするだけで閲覧できるWebサービス「BAGET 2.0 (Bacterial and Archaeal Gene Exploration Tool)」を開発した。ユーザーが提供したアノテーションされたゲノムも処理することができる。BAGET 2.0は毎日更新されるローカルデータベースに依存している。
 BAGET 2.0は、ChromeFirefox、Edge、OperaSafariなどの現行ブラウザに対応しているほか、Internet Explorer 11をサポートしている。BAGET 2.0は https://archaea.i2bc.paris-saclay.fr/baget/ から自由にアクセスできる。

 

help

https://archaea.i2bc.paris-saclay.fr/baget/Home/Help#genome_selection

 

webサービス

 

https://archaea.i2bc.paris-saclay.fr/baget/にアクセスする。f:id:kazumaxneo:20210226235629p:plain

 

ゲノムを指定するか、ファイルをアップロードする。

f:id:kazumaxneo:20210227103759p:plain

 

ファイルをアップロードする場合、GenBankファイル(full)を指定する(DNA配列が含まれているGenBnak(ジェンバンク)ファイル)。NCBIからファイルをダウンロードする際には、アノテーションの下にDNA配列が追加されるように、GenBank (full)オプションにチェックをつける。

f:id:kazumaxneo:20210227105917p:plain

 

FTPサーバーからダウンロードするなら、gbff.gzファイルを選ぶ。f:id:kazumaxneo:20210227110038p:plain

アップロードが完了するか、リストからゲノムが選ばれると、右のGene listウィンドウに認識された全ての遺伝子が表示される。

f:id:kazumaxneo:20210227110250p:plain

 

何かの遺伝子の配列を取り出す。右のウィンドウの遺伝子をクリックするか、左のウィンドウで検索する。ここではDnaAと検索。

f:id:kazumaxneo:20210227110539p:plain

 

 

片方はDnaA遺伝子ではなく、DnaAに関連するタンパク質をコードする遺伝子だった。DnaA遺伝子の方をクリックする。下のウィンドウにDnaA遺伝子の配列が表示された。

f:id:kazumaxneo:20210227110845p:plain

赤いのがDnaA遺伝子の塩基配列。周辺にコードされている遺伝子も表示されている。

下にはタンパク質配列も表示されている。

f:id:kazumaxneo:20210227111040p:plain

 

他のORFをダブルクリックすると、そのORFの配列に切り替わる。

f:id:kazumaxneo:20210227111210p:plain

 

サイズは変更できる。30kbに変更。

f:id:kazumaxneo:20210227113829p:plain

この30kbの配列をすべてダウンロードするならFull sequenceをクリックする。DNA配列がfasta 形式でダウンロードされる。

 

  • 提供されたゲノム情報は、サーバーによってロードされ、RAMで完全に処理される。
  • ユーザーデータは単一のユーザーセッションに制限されており、同時セッション間で共有されることはない。
  • すべてのユーザーデータは、セッションの終了時にメモリから完全に削除される。処理中に管理者がアクセスできることもない。
  • BAGETは、選択された遺伝子のアクセッション/バージョンを用いて、NCBIリポジトリへの直接外部リンクを提供している。
  • BAGETは、選択された遺伝子のアクセッション/バージョンを用いて、Uniprotのウェブサイトへの直接外部リンクを提供している。

 

PDFレポート機能もあるが、試した時は動作しなかった。

引用

BAGET 2.0: an updated web tool for the effortless retrieval of prokaryotic gene context and sequence
Benjamin Hepp, Violette Da Cunha, Florence Lorieux, Jacques Oberto
Bioinformatics, Published: 03 February 2021

 

BURST

 

 次世代のDNAシーケンシングデータが計算能力が追いつかないほどの速さで出現しているため、基本的なDNAアライメント/マッピングの問題に対する近似ヒューリスティックな解法がますます使われるようになってきている。逆説的なことに、データが増えれば増えるほど、それを解析するために使用されるアライメントアルゴリズムの精度が低下する。ミスマッチがあるという制約の下で、完全な感度と特異性を持つアルゴリズムはより高速なアライメントを約束する技術のために無視されてきた。BURSTは、証明された最適なアライメントアルゴリズムのルーツに戻り、デフォルトの動作モードでアライメント品質を犠牲にすることなく、数百万倍のスピードアップで働く。
 BURST は、数学的に最適化されたハイスループットなエンドツーエンドのショートリード DNA アライナーである。BURSTは、292-bpの1,200万微生物アンプリコン配列をGreengenes 参照データベースに対して、E7-4850v2サーバーでは10分以内、2013年製デュアルコアMacbook Airでは数時間でアラインメントすることができる。32コアのIvy Bridgeサーバー1台で、RefSeqの31.5GBの完全ゲノムサブセットに対して、毎秒10,000以上の100 bpリードを98%のアラインメント同一性でアラインメントできる。オプションのヒューリスティックモードにより、より低い類似度スコアでより高速なアラインメントを可能にする。

 

インストール

Github

リリースより使っているOS (WindowsLinuxMac)に適した実行形式ファイルをダウンロードする。

https://github.com/knights-lab/burst/releases

./burst_linux_DB15

# ./burst_linux_DB15

This is BURST [v1.0 DB 15]

 

BURST aligner (v1.0; DB15 version)

Compiled with AVX-128 and multithreading

 

Basic parameters:

--references (-r) <name>: FASTA/edx DB of reference sequences [required]

--accelerator (-a) <name>: Creates/uses a helper DB (acc/acx) [optional]

--queries (-q) <name>: FASTA file of queries to search [required if aligning]

--output (-o) <name>: Blast6/edb file for output alignments/database [required]

 

Behavior parameters:

--forwardreverse (-fr): also search the reverse complement of queries

--whitespace (-w): write full query names in output (include whitespace)

--xalphabet (-x): Allow any alphabet and disable ambiguity matching

--nwildcard (-y): Allow N,X to match anything (in query and reference)

--taxonomy (-b) <name>: taxonomy map (to interpolate, use -m CAPITALIST)

--mode (-m) <name>: Pick an alignment reporting mode by name. Available modes:

  BEST (report first best match by hybrid BLAST id)

  ALLPATHS (report all ties with same error profile)

  CAPITALIST (minimize set of references AND interpolate taxonomy) [default]

  FORAGE (report all matches above specified threshold)

  ANY (report any valid hit above specified threshold)

--makedb (-d) [name qLen]: Create a database from input references

  [name]: Optional. Can be DNA, RNA, or QUICK [QUICK]

  [qLen]: Optional. Max query length to search in DB [500]

 

Performance parameters:

--dbpartition (-dp) <int>: Split DB making into <int> chunks (lossy) [1]

--taxacut (-bc) <num>: allow 1/<int> rank discord OR % conf; 1/[10]

--taxa_ncbi (-bn): Assume NCBI header format '>xxx|accsn...' for taxonomy

--skipambig (-sa): Do not consider highly ambiguous queries (5+ ambigs)

--taxasuppress (-bs) [STRICT]: Suppress taxonomic specificity by %ID

--id (-i) <decimal>: target minimum similarity (range 0-1) [0.97]

--threads (-t) <int>: How many logical processors to use [4]

--shear (-s) [len]: Shear references longer than [len] bases [500]

--fingerprint (-f): Use sketch fingerprints to precheck matches (or cluster db)

--prepass (-p) [speed]: use ultra-heuristic pre-matching

  [speed]: Optional. Integer, maximum search effort [16]

--heuristic (-hr): allow relaxed comparison of low-id matches

--noprogress: suppress progress indicator

 

--help (-h): Shows this help screen with version info

 

Example: burst -r myRefs.fasta -q myQs.fasta -o outputs.txt -i 0.98

 

Licensed under the GNU Affero General Public License v3.0

 

 

実行方法

1、データベースの作成

burst -r MyDB.fasta -a MyDB.acx -o MyDB.edx -d DNA -s

 

2、Search

burst -q myQueries.fasta -a MyDB.acx -r MyDB.edx -o output.txt

 

 

引用

BURST enables mathematically optimal short-read alignment for big data
Gabriel Al-Ghalith,Dan Knights

bioRxiv, Posted September 08, 2020

 

関連


ディープラーニングを用いた微生物ゲノムのビンニングツール Vamb

 

 メタゲノミクスワイドゲノム配列データからの微生物種の同定と再構築は、重要かつ挑戦的な課題である。現在の既存のアプローチは、複数のサンプルにわたる遺伝子またはコンティグの共分散情報と、配列中のk-mer組成情報に依存している。ここでは、最近のディープラーニングの進歩を利用して、変分自動エンコーダーを用いて、クラスタリングに先立って共分散情報と組成情報を符号化するアルゴリズムを開発した。本研究では、ディープネットワークがこれら2つの異種データセットを事前知識なしに統合できること、また、3つの異なるベンチマークデータセットから1.8倍から8倍の精度で完全なゲノムビンを再構成することで、既存の最先端技術を凌駕することを示している。さらに、ヒト腸内マイクロバイオームの約1,000万の遺伝子と1,270サンプルの遺伝子カタログにも適用した。結果、130万から180万の遺伝子をクラスタリングし、117から246個の高精度で完全な遺伝子群を再構成することができ、そのうち70個の遺伝子群は従来の方法と比較して全く新しいものであった。本手法 Variational Autoencoders for Metagenomic Binning (VAMB) は、https://github.com/jakobnissen/vamb で自由に利用できる。

 

インストール

Github

#conda、ここでは高速なmambaを使う
mamba create -n Vamb -y
conda activate Vamb
mamba install -c pytorch pytorch torchvision cudatoolkit=10.2
mamba install -c bioconda vamb -y

#pip
pip install https://github.com/RasmussenLab/vamb/archive/3.0.2.zip

> vamb

$ vamb 

usage: vamb outdir tnf_input rpkm_input [options]

 

Vamb: Variational autoencoders for metagenomic binning.

 

    Default use, good for most datasets:

    vamb --outdir out --fasta my_contigs.fna --bamfiles *.bam

 

    For advanced use and extensions of Vamb, check documentation of the package

    at https://github.com/RasmussenLab/vamb.

 

Help:

  -h, --help          print help and exit

 

Output (required):

  --outdir            output directory to create

 

TNF input (either fasta or all .npz files required):

  --fasta             path to fasta file

  --tnfs              path to .npz of TNF

  --names             path to .npz of names of sequences

  --lengths           path to .npz of seq lengths

 

RPKM input (either BAMs, JGI or .npz required):

  --bamfiles  [ ...]  paths to (multiple) BAM files

  --rpkm              path to .npz of RPKM

  --jgi               path to output of jgi_summarize_bam_contig_depths

 

IO options:

  -m                  ignore contigs shorter than this [100]

  -s                  ignore reads with alignment score below this [None]

  -z                  ignore reads with nucleotide identity below this [None]

  -p                  number of subprocesses to spawn [min(8, nbamfiles)]

  --norefcheck        skip reference name hashing check [False]

  --minfasta          minimum bin size to output as fasta [None = no files]

 

VAE options:

  -n  [ ...]          hidden neurons [Auto]

  -l                  latent neurons [32]

  -a                  alpha, weight of TNF versus depth loss [Auto]

  -b                  beta, capacity to learn [200.0]

  -d                  dropout [Auto]

  --cuda              use GPU to train & cluster [False]

 

Training options:

  -e                  epochs [500]

  -t                  starting batch size [256]

  -q [ [ ...]]        double batch size at epochs [25 75 150 300]

  -r                  learning rate [0.001]

 

Clustering options:

  -w                  size of window to count successes [200]

  -u                  minimum success in window [20]

  -i                  minimum cluster size [1]

  -c                  stop after c clusters [None = infinite]

  -o                  binsplit separator [None = no split]

 

 

実行方法

1、アセンブリ配列の結合(複数ある場合)

 concatenate.py /path/to/catalogue.fna.gz sample1/contigs.fasta
sample2/contigs.fasta

 

2、リードのアセンブリへのマッピング

minimap2 -d catalogue.mmi catalogue.fna.gz
minimap2 -t 8 -N 50 -ax sr catalogue.mmi sample1_R1.fq.gz sample1_R2.fq.gz | samtools view -F 3584 -b --threads 8 > sample1.bam
minimap2 -t 8 -N 50 -ax sr catalogue.mmi sample2_R1.fq.gz sample2_R2.fq.gz | samtools view -F 3584 -b --threads 8 > sample2.bam

 

3、ビニング。デフォルトではfasta配列を出力しないので、--minfastaオプションをつける。

vamb --outdir outdir --fasta catalogue.fna.gz  --bamfiles BAM1 BAM2  -o C --minfasta 200000

 

 

 

 Recommended workflow (Githubから)

  • 1) リードの前処理と品質チェック

FastQCと組み合わせたAdapterRemovalを使用している。

  • 2) 各サンプルを個別にアセンブルし、コンティグを取り出す。

各サンプルで個別に metaSPAdes を使用することを勧める。シングルサンプルをアセンブルした後、サンプルごとにビンスプリットすると最良の結果が得られるためです。

  • 3) コンティグヘッダがすべてユニークであることを確認しながら、fastaファイルを連結し、小さなコンティグをフィルタリングする。

重要: Vamb はシーケンスのエンコードニューラルネットワークを使用しており、ニューラルネットワークは小さなデータセットでオーバーフィットする。コンティグ数が 50,000 以下のデータセットではテストしていない。また、非常に短い配列をビニングしようとするべきではない。入力配列の長さのカットオフを決める際には、カットオフを低くしすぎると、ビンに入れにくいコンティグが残ってしまい、高くしすぎると、良いデータを捨ててしまうことになる。デフォルトでは2000bpの長さのカットオフを使用している。

コンティグヘッダはユニークでなければならない。さらに、ビンスプリットを使用したい場合、コンティグヘッダは {Samplename}{Separator}{X} という形式でなければならない。

Vambはメモリ効率が良いが、データセットが大きすぎてRAMに収まりきらない場合は、MASHのようなツールを使用して、類似のサンプルをより小さなバッチにまとめ、それらのバッチを個別にビニングすることもできる。

リードを複数のコンティグにマップできるようにすると、より正確なアバンダンスの推定が可能になる。ただし、1つのリードに対するすべてのBAMレコードはBAMファイル内で連続していなければならない。これは、ほとんどすべてのアライナーの出力におけるデフォルトの順序だが、アライメント位置でソートされた BAM ファイルを使用している場合や、マルチマッピングリードを使用している場合は、最初にリード名でソートする必要がある。

一般的に、コンティグAからのリードがコンティグBにアラインメントされている場合、VambはAとBを一緒に結合する。一般的に、コンティグ A からのリードをコンティグ B に合わせると、Vamb は A と B を一緒に結合する。minimap2を使用することを勧める。

MAPQの低いアラインメントをフィルタリングしてしまうと、相同配列へのアラインメントが除去されてしまい、深度推定にバイアスがかかる。

MetaBAT2のjgi_summarize_bam_contig_depthsプログラムは、Vambのparsebamモジュールよりも正確にBAMの深さを推定することがわかった。最良の結果を得るためには、MetaBAT2をダウンロードし、jgi_summarize_bam_contig_depthsを使って深さを推定し、Vambを--bamfilesの代わりに--jgiで実行することを勧める。また、snakemakeワークフローを使用することも検討してください。

  • 5) ランバンブ

デフォルトでは、Vamb はビンの FASTA ファイルを出力しない。レポジトリの例では、--minfasta 200000 オプションが設定されており、これはサイズが 200kbp 以上のすべてのビンが FASTA ファイルとして出力されることを意味する。

 

引用

Improved metagenome binning and assembly using deep variational autoencoders

Jakob Nybo Nissen, Joachim Johansen, Rosa Lundbye Allesøe, Casper Kaae Sønderby, Jose Juan Almagro Armenteros, Christopher Heje Grønbech, Lars Juhl Jensen, Henrik Bjørn Nielsen, Thomas Nordahl Petersen, Ole Winther & Simon Rasmussen
Nature Biotechnology (2021)

s

Binning microbial genomes using deep learning

Jakob Nybo Nissen, Casper Kaae Sønderby, Jose Juan Almagro Armenteros, Christopher Heje Grønbech, Henrik Bjørn Nielsen, Thomas Nordahl Petersen, Ole Winther, Simon Rasmussen

bioRxiv, Posted December 19, 2018

 

 

 

 

GO エンリッチメント解析を実行し、バックグラウンドセットと比較して過剰に存在する語彙を調べる FunSet

 

 遺伝子オントロジーエンリッチメント解析は、複雑な生物学的データセットから意味のある情報を抽出する効果的な方法を提供する。遺伝子セットの中で有意に過剰発現している語彙を特定することで、研究者は遺伝子が共有する生物学的特徴を明らかにすることができる。濃縮された語彙を抽出することに加えて、結果を生物学的解釈に資する方法で可視化することも重要である。
 ここでは、エンリッチメント解析を実行し、可視化するための新しいウェブサーバーであるFunSetを紹介する。このウェブサーバーは、バックグラウンドセットと比較してターゲットセットに統計的に過剰に存在する遺伝子オントロジーの語彙を特定する。エンリッチメントされた語彙は、語彙間の意味的類似性をキャプチャする2Dプロットに表示される。またスペクトル・クラスタリングによって語彙はクラスタリングされ、各クラスタの代表的な語彙を識別するオプションもある。FunSetは対話的にもプログラム的にも使用でき、ユーザーはエンリッチメント結果を表形式、SVGファイルとしてのグラフ形式、JSONcsvとしてのデータ形式の両方でダウンロードすることができる。分析の再現性を高めるために、ユーザーはオントロジーアノテーションの過去のデータにアクセスすることができる。スタンドアロンプログラムとWebサーバーのソースコードオープンソースライセンスで提供されている。

 

 

webサービス

FunSetにアクセスする。

f:id:kazumaxneo:20210224041326p:plain

 

遺伝子名をカンマ区切りで記載する。ここではexampleデータを入力。

f:id:kazumaxneo:20210223203142p:plain

(遺伝子を指定するためのフォーマットは、ヒトの場合はHGNCシンボル、ウシとイヌの場合はVGNCシンボル、モデル生物の場合はMODシンボル。)

 

エンリッチメント解析には、ターゲットセット(関心のある特性を持つ遺伝子)とバックグラウンドセットが必要。オプションで、バックグラウンドの遺伝子セットをアップロードする。アップロードしなければ、デフォルトでは、FunSetは選択した生物のすべてのアノテーションされた遺伝子をバックグラウンドとして選択する。

f:id:kazumaxneo:20210224041933p:plain

 

Ontologyのバージョンを指定する。

f:id:kazumaxneo:20210224041513p:plain

 

種名を選択する。デフォルトはhuman。

f:id:kazumaxneo:20210224041526p:plain

 

GOのカテゴリを指定する。

f:id:kazumaxneo:20210224041649p:plain

 

FDRを指定してサブミットする。

f:id:kazumaxneo:20210224041708p:plain

 

 

結果が表示されるまで数分かかった。

出力。

f:id:kazumaxneo:20210223203416p:plain

GOエンリッチメント解析結果の視覚化は、類似度インデックスから得られる距離行列上に多次元スケーリング(MDS)を使用して配置される。

 

右上にサマリーが表示される。自動で検出されたクラスタ数は2つだった(2色のプロットがそれぞれのクラスタを表す)。

f:id:kazumaxneo:20210224042319p:plain

 

図はインタラクティブに操作でき、マウスホイールで拡大縮小、左クリックしてドラッグすることでパンできる。語彙のノードを左クリックしてドラッグすれば、ノードを移動できる。

f:id:kazumaxneo:20210223204848p:plain

 

 

クラスタのノードをクリックすると、FDRやエンリッチメントサイズ(ES)などの関連データとともにエンリッチメントされた語彙が表示される。このパネルでは、特定の語彙をクリックすると、SVGグラフ内で赤くハイライト表示される。また、このパネルで語彙をクリックすると、語彙の説明も表示される。

f:id:kazumaxneo:20210223205148p:plain

GO termはAmigoにリンクされている。

 

右上のバーからクラスタ数は任意で変更できる。論文と同じ11に変更した。バーの操作では細かい指定が難しいので、数値を直接タイプする。

f:id:kazumaxneo:20210224042608p:plain

クラスタのサイズの上限は入力された語彙の数になる。ここでは、理屈上は473まで上がる(クラスタリングの意味はなくなる)。

 

結果が反映されるまで10秒程度かかった。各クラスタは自動で色分けされる。基本的にMDSの距離に応じてクラスタリングされているのが分かる。

f:id:kazumaxneo:20210224042733p:plain

 

クラスタのチェックマークを外すと、一時的に透明度を上げることができる。

f:id:kazumaxneo:20210224042937p:plain

 

繰り返しになるが、各クラスタのノードをクリックすると、クラスタを構成している語彙をリストで確認できる(下の写真の右側)。リストはFDRの低い順でソートされており、ノードのサイズはこのFDRを反映している(=> エンリッチ度合いを反映している)。

f:id:kazumaxneo:20210224043905p:plain

 

リスト左端のボタンをクリックすると青くなり、該当するノードが赤枠でハイライト表示される。注釈も付く。

f:id:kazumaxneo:20210224044423p:plain


図ははSVGファイルとしてダウンロードでき、濃縮された語彙と対応する遺伝子は表形式でダウンロードできます。図は現在表示されている領域がダウンロードされることに注意して下さい。

f:id:kazumaxneo:20210224044818p:plain

引用

FunSet: an open-source software and web server for performing and displaying Gene Ontology enrichment analysis
Matthew L. Hale, Ishwor Thapa & Dario Ghersi
BMC Bioinformatics volume 20, Article number: 359 (2019)

 

関連