コード領域のアミノ酸配列を考えてマルチプルアライメントを行うMUCSE

塩基配列からコード領域のアミノ酸配列を予測してマルチプルアライメントを行う場合、従来はギャップやミスを補正せず全ての配列をアミノ酸に変換してアライメントを行なっていた。しかしこのような一義的に変換する方法だと、シーケンスエラーや擬遺伝子のstopコドンに引っ張られてアライメント精度がおかしくなる場合があった。

オーサーの発表したMUCSEは、stopコドンなどにペナルティを設定して、擬遺伝子を含むリストでも従来のNeedleman-Wunsch algorithmと同等の精度とスピードでマルチプルアライメントを行う方法論である。

公式サイト

http://bioweb.supagro.inra.fr/macse/index.php?menu=intro

チュートリアルファイル

http://bioweb.supagro.inra.fr/macse/index.php?menu=download

ロシア語のチュートリアルが多く分かりにくい。

インストール

brewで導入できる。

brew install MUCSE

実行方法

brewで導入したなら"macse -prog <commnad>"の形で使う。

alignSequences

元のfastaファイル（先頭5配列）。

このような短いコード領域の配列なら、デフォルトのパラメーターでランする。

macse -prog alignSequences -seq tmem184a.fasta

数分で解析は終わる。

塩基出力（tmem184a_macse_NT.fasta）。

 アミノ酸出力（tmem184a_macse_AA.fasta）。

アミノ酸ベースのマルチプルアライメントなので、塩基で確認するとアライメント精度が低い部位もある。

出力ファイル名は-out_NTと-out_AAで変更できる。

 genetic codeがスタンダードではない生物の場合、コードを明示して解析する。例えばクロロプラストの遺伝子だと

macse -prog alignSequences -seq poacea_matk.fasta -gc_def 11

 NCBI genetic code（http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi）

擬遺伝子も混じっている可能性がある（stopコドンが入る可能性がある）データセットで、擬遺伝子がわかっているなら、ファイルを別々に入力し、ペナルティスコアを個別に与える。この方が精度が高くなる。

macse -prog alignSequences -seq ENAM_genes.fasta -seq_lr ENAM_pseudos.fasta -fs_lr 10 -stop_lr 10

• -seq　"reliable" sequences
• -seq_lr　"less reliable" sequences

 上ではpseudo geneのペナルティスコアを-fs_lr 10 -stop_lr 10にし、信頼できる遺伝子のペナルティスコア（-fsと-stop）をデフォルトでランしている。

NGSリードやcontigを使う場合

NGSのデータを直接使うとシーケンスエラーの可能性があるので、公式サイトではペナルティの与え方をはじめは-fs_lr 10 -stop_lr 15にして、後でチューニングするよう指示している。

macse -prog alignSequences -seq TMEM214.fasta -seq_lr TMEM214.fasta -fs_lr 10 -stop_lr 15

マルチプルアライメントはかなり計算が重くなる。100以上の配列があったりすると、計算が終わるのに非常に時間がかかる。

 splitAlignment

アライメント結果から、リストで与えた配列を別ファイルに出力する。

macse -prog splitAlignment -align tmem184a_NT.aln -subset primate_seqId.list

入力したリストは以下のようなものである。

user$ cat primate_seqId.list 

>Pongo

>Homo

>Gorilla

>Macaca

>Pan

他にもいくつかのコマンドがある。詳細は公式ページで確認してください。

2016年にversion2にアップデートされ、2018年夏にv2の論文が出ました。

https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msy159/5079334

引用

MACSE: Multiple Alignment of Coding SEquences Accounting for Frameshifts and Stop Codons

Vincent Ranwez , Sébastien Harispe, Frédéric Delsuc, Emmanuel J. P. Douzery

PLoS One. 2011;6(9):e22594

Needleman-Wunsch algorithm

グローバルアライメント | グローバルアライメントを求める Needleman–Wunsch アルゴリズム

メタゲノムデータをbinningして種を予測するMBBC

MBBCはメタゲノムをbinningする方法論。リード中のk-mer頻度とk-merカバレッジから分類とabundanceの見積もりを行う。2015年に論文が発表された。

マニュアル

http://eecs.ucf.edu/~xiaoman/MBBC/man1V1.html

インストール

ダウンロード

実行方法

GUIバージョンとターミナルで動かすCUIバージョンが提供されている。ここではGUI版の流れを説明する。

exampleデータを解析する。GUI-jarを開く。

 左上のopenをクリック

example/paired-end_reads/ba3mg8.fnaを選択（ペアードエンド）。

他のパラメータはdefaultのままにした。

Runボタンをクリックすると解析がスタートする。

数分で解析は終了する。fasta入力ディレクトリにba3mg8.fna_filtered_predictionsが出力される。

xneo$ cat /Users/kazumaxneo/Downloads/MBBC/MBBC_examples/paired-end_reads/ba3mg8.fna_filtered_predictions 

Predictions after EM algorithms: 

Predicted genome size: 1151747.76 1162868.69  

Predicted alpha: 28.55% 71.45%

Predicted lambda: 2.57 6.38

There are 2 species predicted.

Prediction after reads binning: 

Predicted genome size: 1049007.03 1204309.40  

Predicted alpha: 26.00% 74.00%

Predicted lambda: 2.57 6.38

uesaka-no-MacBook-Air-2:condetri-master kazumaxneo$

ツールに対しての感想

GUI実装されているが、シンプルなインターフェイスなので、アルゴリズムさえ構築されていればjavaのガワの作成は容易に感じる。統合開発環境に慣れた人なら１、２時間くらいで作れるんじゃないだろうか？もう少しこうゆうツールが出てきてもいいように思う。

引用

MBBC: an efficient approach for metagenomic binning based on clustering

Ying Wang, Haiyan Hu, Xiaoman Li

BMC Bioinformatics. 2015; 16: 36.

クオリティトリミングを行う condetri

condetriはペアリードを考量してクオリティトリミングが行えるperlのツール。

公式サイト

https://code.google.com/archive/p/condetri/

マニュアル

ダウンロードしたディレクトリにPDFマニュアルあり。

インストール

本体はperlスクリプトである。

GIthub

ダウンロードしてパスを通す。

実行方法

perl condetri.pl -fastq1=R1.fq -fastq2=R2.fq -prefix=output -cutfirst=i -hq=i -lq=i -frac=[0,1] -minlen=i -mh=i -ml=i -sc=i -rmN

• -fastq1=file 　Fastq(.gz) file. If a second file is given, the files are trimmed-fastq1=file Fastq(.gz) file. If a second file is given, the files are trimmed
• -fastq2=file　as a pair. The reads must have the same order in both files.
• -prefix=string　Prefix for the output file(s). The filtered fastq file(s) will be named prefix_trim1.fastq (and prefix_trim2.fastq if present). For pairs,  a third file will be given with unpaired reads (reads from pairs where one   low quality read has been removed).
• -cutfirst=i　Remove i first bases from the 5'end before any trimming [0].-cutfirst=i Remove i first bases from the 5'end before any trimming [0].
• -cutlast=i　Remove i bases from the 3'end before any trimming [0].-cutlast=i Remove i bases from the 3'end before any trimming [0].
• -rmN　 Remove non-ATCG bases from 5'end before any trimming [no].
• -hq=i　Hiqh quality threshold [25].-hq=i Hiqh quality threshold [25].
• -lq=i　Low quality threshold [10].
• -frac=[0,1]　 Fraction of read that must exceed hq [0.8].
• -lfrac=[0,1]　 Maximum fraction of bases with qual<lq [0].
• -minlen=i　 Min allowed read length [50].

ペアリードの順番が壊れないように動作するトリミングツールは他にもいくつかあります。Trimmomaticなどがよく知られていますが、ここではsickleを紹介しています。

condetriのパッケージにはPCRのduplicationを除くスクリプトも同梱されている。

引用

ConDeTri - A Content Dependent Read Trimmer for Illumina Data.

Smeds L, Künstner A (2011)

PLoS ONE 6(10): e26314. doi:10.1371/journal.pone.0026314

マルチプルアライメントを行う T-Coffee

 T-Coffee（Tree based Consistency Objective Function For AlignmEnt Evaluation）はマルチプルアライメントを行うツールである。始めに２つずつ配列を比較し、それから全部の配列を使いマルチプルアライメントを実行する。従来のclustalより高速に動作する。より大規模な解析ではclustal omegaなどの選択肢があるが、普通の解析ではT-Coffee で十分である。

公式サイト

http://www.tcoffee.org/Projects/tcoffee/

T-Coffeeマニュアル

http://tcoffee.readthedocs.io/en/latest/

T-Coffee webサーバー（EMBL-EBI）

http://www.ebi.ac.uk/Tools/msa/tcoffee/

インストール

公式サイトより自己解凍形式の.dmgファイルをダウンロードする。インストーラーを解凍し、指示通り進める。

インストールが終わったら、sourceしてターミナルで導入を確認する。

source ~/.bash_profile 
t_coffee -version

ラン 

マルチファスタファイルを解析する。

t_coffee input.fasta

標準では全コアを使い解析が行われる。計算が終わると.aln、.dnd、htmlファイルが出力される。

.alnファイル。clustalXのAppend Sequencesで開くこともできる。

head -20 /Users/user/Desktop/ISy203alignment.aln 

CLUSTAL FORMAT for T-COFFEE Version_11.00.8cbe486 http://www.tcoffee.org [MODE:  ], CPU=0.01 sec, SCORE=999, Nseq=10, Len=1174 

ISY203a         CAGAAGTGTTGAACGATAGTTATAAAGAGAAA-AAAGCTCTTTAAAATGA

ISY203b         CAGAAGTGTTGAACGATAGTTATAAAGAGAAAAAAAGCCCTTTAAAATGA

ISY203c         CAGAAGTGTTGAACGATAGTTATAAAGAGAAAAAAAGCCCTTTAAAATGA

ISY203d         CAGAAGTGTTGAACGATAGTGATAAAGAGAAAAAAAGCTCTTTAAAATGA

ISY203e         CAGAAGTGTTGAACGATAGTTATAAAGAGAAAAAAAGCCCTTTAAAATGA

ISY203f         CAGAAGTGTTGAACGATAGTTATAAAGAGAAAAAAAGCCCTTTAAAATGA

ISY203g         CAGAAGTGTTGAACGATAGTTATAAAGAGAAAAAAAGCCCTTTAAAATGA

ISY203j         CAGAAGTGTTGAACGATAGTTATAAAGAGAAAAAAAGCCCTTTAAAATGA

ISY203k         CAGAAGTGTTGAACGATAGTTATAAAGAGAAAAAAAGCCCTTTAAAATGA

ISY203x         CAGAAGTGTTGAACGATAGTTATAAAGAGAAAAAAAGCCCTTTAAAATGA

                ******************** *********** ***** ***********

ISY203a         TAGAATAAAGGGCAGAAATATTAGAATAATTGAAGCGATGATATCAAATT

ISY203b         TAGAATAAAGGGCAGAAATATTAGAATAATTGAAGCGATGATATCAAATT

ISY203c         TAGAATAAAGGGCAGAAATATTAGAATAATTGAAGCGATGATATCAAATT

ISY203d         TAGAATAAAGGGCAGAAATATTAGAATAATTGAAGCGATGATATCAAATT

ISY203e         TAGAATAAAGGGCAGAAATATTAGAATAATTGAAGCGATGATATCAAATT

ISY203f         TAGAATAAAGGGCAGAAATATTAGAATAATTGAAGAGATGATATCAAATT

html出力。

タンパク質のaccurateモード（アミノ酸のみ対応）。

t_coffee input.fasta -mode accurate

RNAのaccurateモード。

t_coffee input.fasta -mode rcoffee

高速モード。

t_coffee input.fasta -mode quickaln

低メモリモード（RNAに対応）。

t_coffee input.fasta -mode memory

引用

T-Coffee: A Novel Method for Fast and Accurate Multiple Sequence Alignment

CeÂdric Notredame1,2,3*, Desmond G. Higgins4 and Jaap Heringa1

J. Mol. Biol. (2000) 302, 205±217

メタゲノムの自動解析パイプライン MyCC

2019 7/6 インストール、ラン追記

MyCCは全プロセスを自動化したメタゲノム解析ツール。contigのfastaファイルを入力すると、配列の特性に従って自動で分類し、binning向けに色がついた図を描画し、さらにクラスタリングされたfastaまで出力することができる。既存のカバレッジやペアリードのつながりなどで分類する手法と比較して、クラスタリング精度が高いと言われている。また、全プロセスは完全に自動化されている。dockerのイメージや仮想イメージからMyCCを起動すれば、fastaを入力するだけで全て自動でランして結果を返してくれる。

 解析の流れは以下のようになっている。

１、prodicalでコード領域を予測。

２、fetchMGを使い、シングルコピーでユニバーサルなオーソログ遺伝子40を抽出。

３、speciesレベルのマーカー遺伝子をUCLUSTで抽出。

４、4-mer頻度からgenomic signitureを分析し、Barnes-HUt-SNEで２次元に分類。

５、非階層的クラスタリング手法のaffinity propagationを使いplotをクラスタリング。

インストール

ランにはfetchMG、Prodigal、UCLUSTなどが必要でであるが、オーサーらによりovaファイルが配布されている。virttual PC内で仮想イメージを起動することで、すぐにランできるようになっている。

 ovfファイルのダウンロード

https://sourceforge.net/projects/sb2nhri/files/MyCC/

• Virtualboxダウンロード=> リンク(OSXを選択）
• VMware=>リンク

VMwareへのインポートはこちらでは検証していません。 

 他にもdockerイメージやソースファイルが提供されている。

sb2nhri - Browse /MyCC at SourceForge.net

dockerイメージ

docker pull 990210oliver/mycc.docker:v1

#run
docker run 990210oliver/mycc.docker:v1 MyCC.py

> docker run 990210oliver/mycc.docker:v1 MyCC.py -h

Usage:

MyCC.py [inputfile] [4mer/5mer/56mer default:4mer] [Options]

Options:

-a A file having coverage information. 

-t Minimum contig length. [defaults: 1000 bp]

-lt A fraction of contigs for the first-stage clustering. [default: 0.7]

-ct Minimum contig lengh for first stage clustering. [bp]

-meta Change to meta mode of Prodigal. [default: single]

-p Perplexity for BH-SNE. [default: 20, range between 5 and 50]

-pm To set preferences all equal to the median of the other similarities for affinity propagation.

-st Maximum distance for sparse format of affinity propagation. [default: 500]

-mask To mask repetitive sequences.

-keep To keep temporary files.

ラン

virtual PC内にfastaを読み込ませる。

ターミナルからMyCCをラン。

MyCC.py metagenomic_contigs.fasta

dockerイメージを使ってランするなら

#カレントにmetagenome_contigs.fastaがある状態で以下のように打つ
docker run --rm -itv $PWD:/data/ -w /data \
 990210oliver/mycc.docker:v1 MyCC.py metagenome_contigs.fasta

テストでIndex of /dna/RD/Metagenome_RD/MetaBAT/Software/MockupのMockメタゲノム（25 bacteriaを混ぜた擬似メタゲノムのシーケンスデータ）のアセンブルデータを解析してみた。100MBのfastaファイルだったが、解析は30分程度で終了した。

出力された図。

23クラスター認識されている。

Cluster.summaryを開く。

 user$ cat /Users/user/Downloads/Cluster.summary 

Cluster WholeGenome N50 NoOfCtg LongestCtgLen AvgLenOfCtg Cogs

Cluster.1.fasta 4070502 70156 97 268872 21423 35

Cluster.2.fasta 4256229 524907 13 896543 192492 35

Cluster.3.fasta 11381431 24260 981 325606 5820 21

Cluster.4.fasta 3613209 84946 92 176933 19893 35

Cluster.5.fasta 3944578 1281749 7 1811839 441941 35

Cluster.6.fasta 4833441 1145236 12 1295155 203365 42

Cluster.7.fasta 2161903 112927 36 253843 32019 36

Cluster.8.fasta 8313943 163761 96 478860 44284 35

Cluster.9.fasta 3259866 254968 28 556331 63405 36

Cluster.10.fasta 3624126 555508 11 777357 178603 37

Cluster.11.fasta 4329879 100800 68 289305 32544 39

Cluster.12.fasta 4735197 604725 11 1268240 226359 32

Cluster.13.fasta 5575100 131134 60 388470 47221 33

Cluster.14.fasta 3261204 901082 15 1093628 132980 35

Cluster.15.fasta 1862577 233546 14 499879 81212 36

Cluster.16.fasta 5338418 159508 67 481874 41626 34

Cluster.17.fasta 4012035 625944 20 1020404 125987 35

Cluster.18.fasta 3176251 454949 13 747195 135258 37

Cluster.19.fasta 3611880 32457 178 121052 10257 36

Cluster.20.fasta 1749350 134168 26 274886 38596 33

Cluster.21.fasta 4087425 443282 16 712360 147145 36

Cluster.22.fasta 5209287 387726 27 926513 99767 35

Cluster.23.fasta 2125683 1546362 5 1546362 309272 36

 そのほか、affinity propagationによってクラスタリングされたfastaも保存されている。

リアルデータ（Lake_Huron_Sinkhole_Photosynthetic_Microbial_Mats_Metagenome）でテストしてみた（リンク）。15 GBx2のペアードエンドfastqのアセンブルにはmetaspadesを用いた（-k auto）。78000 contigできた。

下図：MyCC解析結果。

20クラスター検出された。

引用

Accurate binning of metagenomic contigs via automated clustering sequences using information of genomic signatures and marker genes

Hsin-Hung Lin & Yu-Chieh Liao

Scientific Reports 6, Article number: 24175

https://www.nature.com/articles/srep24175

Alignment-free Visualization of Metagenomic Data by Nonlinear Dimension Reduction

https://www.nature.com/articles/srep04516

DNA解析ソフトに近い機能を提供するwebツール集 SMS

2019 8/7リンク追加

2021 10/3リンクエラー修正

SMSは、NGSの登場よりずっと以前から使われているDNA/プロテインの編集や変換ができるツール集である。昔からあるDNA解析ソフトの大半の機能をカバーしている。webサーバー版とオフラインで動くローカル版がある。いずれもhtmlベースで動作する（内部でjavaが動いている）。登場は古いが現在でも管理されている。ローカル版でもHTMLベースで解析できるので、コマンドラインに不慣れなユーザーも使いやすいと思われる。

トップページ

http://www.bioinformatics.org/sms2/index.html

できること一覧

http://www.bioinformatics.org/sms2/about.html

SMS version1

http://www.bioinformatics.org/sms/

プライマー結合サイトのツールを例に説明する。

Primer Map

exampleデータをランする。制限酵素サイトを示すチェックを入れ、１行に105行表示に変えてラン。

　↓

プライマー結合サイトと制限酵素サイトが表示された。

このようにHTMLベースでデータ処理することが可能である。このようなツールが50近く利用できる。よく使いそうな機能に限定して機能を見ていく。

 シーケンス解析 

DNA配列の指定した領域を抽出する。

http://www.bioinformatics.org/sms2/range_extract_dna.html

アミノ酸配列の指定した領域を抽出する。

http://www.bioinformatics.org/sms2/range_extract_protein.html

DNA配列の逆相補鎖（reverse complementary）を表示する。

http://www.bioinformatics.org/sms2/rev_comp.html

DNAのcodon usageを計算する。

http://www.bioinformatics.org/sms2/codon_usage.html

CpGアイランド探す。

CpG Islands

DNAの分子量（ダルトン）を求める。

http://www.bioinformatics.org/sms2/dna_mw.html

タンパク質の分子量（ダルトン）を求める。

Protein Molecular Weight

 あるパターンの配列を探す（e.g., CTTで次がC or Tを探す）。

http://www.bioinformatics.org/sms2/dna_pattern.html

 似た塩基配列を探す（search homology）。

Fuzzy Search DNA 

似たアミノ酸配列を探す。

http://www.bioinformatics.org/sms2/fuzzy_search_protein.html

逆翻訳してdegenerate DNAを作成。

Multi Rev Trans

少しの変異導入で制限酵素サイトになるDNA配列を探す。

Mutate for Digest

コード領域を探す。

ORF Finder

プライマーのTmやプライマーダイマーの可能性を調べる。

PCR Primer Stats

塩基配列とプライマー配列を入力することでPCR増幅後の塩基配列を返す。

http://www.bioinformatics.org/sms2/pcr_products.html

タンパク質の等電点（pH）を計算。

http://www.bioinformatics.org/sms2/protein_iep.html

制限酵素処理後の配列を表示する。

Restriction Digest

DNAをタンパク質に変換（翻訳）する。

http://www.bioinformatics.org/sms2/translate.html

 (* 6-frameで翻訳するならこちらも便利じゃないでしょうか。)

タンパク質をDNAに逆変換（逆翻訳）する。

http://www.bioinformatics.org/sms2/rev_trans.html

Sequence Figures

アライメント結果を入力してマッチする部位に色をつける。

Color Align Properties

指定文字数で改行やスペースを入れ、DNA配列のポジションをわかりやすくする。

http://www.bioinformatics.org/sms2/group_dna.html

指定文字数で改行やスペースを入れ、タンパク質のポジションをわかりやすくする。

Group Protein

制限酵素サイト一覧を表示。

http://www.bioinformatics.org/sms2/rest_map.html

3-frameで翻訳する。

Translation Map

ランダム変異

指定頻度でDNAに変異を導入する。

http://www.bioinformatics.org/sms2/mutate_dna.html

指定頻度でタンパク質のアミノ酸残基を変化させる。

http://www.bioinformatics.org/sms2/mutate_protein.html

コード領域全長サイズのランダムなDNA配列を発生させる（ネガコンに使える）。

Random Coding DNA

DNAの指定した領域をランダムな塩基配列に置換する。

http://www.bioinformatics.org/sms2/random_dna_regions.html

タンパク質の指定した領域をランダムなアミノ酸配列に置換する。

Random Protein Sequence

DNA配列をランダムにシャッフルする。

http://www.bioinformatics.org/sms2/shuffle_dna.html

アミノ酸配列をランダムにシャッフルする。

Shuffle Protein

genetic code一覧。

http://www.bioinformatics.org/sms2/genetic_code.html

２つの配列をマージする（500000文字まで）。

http://www.bioinformatics.org/sms2/combine_fasta.html

フォーマット変換

EMBLをfastaに変換する。

http://www.bioinformatics.org/sms2/embl_fasta.html

EMBLのDNA配列をアミノ酸配列に変換する。

http://www.bioinformatics.org/sms2/embl_trans.html

Genbankをfastaに変換する。

GenBank to FASTA

GenbankのDNA配列をアミノ酸配列に変換する。

http://www.bioinformatics.org/sms2/genbank_trans.html

DNAから数値などの余計な文字を消去する。

http://www.bioinformatics.org/sms2/filter_dna.html

アミノ酸配列から数値などの余計な文字を消去する。

Filter Protein

アミノ酸１文字表記を３文字表記に変換する。

http://www.bioinformatics.org/sms2/one_to_three.html

アミノ酸３文字表記を１文字表記に変換する。

Three to One

fastaを指定した長さに分割する。

http://www.bioinformatics.org/sms2/split_fasta.html

マルチプルアライメントを行い分子系統樹を作成したければ、以下のMAFFT onlineを利用して下さい。全てオンラインでできます。

引用

The Sequence Manipulation Suite: JavaScript programs for analyzing and formatting protein and DNA sequences

Stothard P

Biotechniques. 2000 Jun;28(6):1102, 1104.

Fri Jun 17 16:17:06 2016 Valid XHTML 1.0; Valid CSS

http://www.bioinformatics.org/sms2/acknowledgments.html

関連

Genbankから遺伝子配列を簡単に取り出す

メタゲノムの簡単なシミュレートを行う BBMap

　メタゲノムをシミュレートするには、ゲノムごとのインサートサイズや増幅biasなどを考慮する必要があり、厳密に行うと計算が複雑になる。また計算リソースも高度に要求される。そのためGPUを使ったシミュレーションツールなども登場している。それに加えて、野外の真のメタゲノムでは1000を超える種がシーケンスされるとも言われている。また環境によって違いが大きすぎるので（＝例外が多すぎる）、腸内細菌叢のような限定された環境でない限り、シミュレーション自体あまり意味がないとも言える。そのためか、メタゲノムをシミュレーション可能な方法論はまだほとんど報告されていない。

　ただしメタゲノムのツールの簡単な検証をするだけならば、菌のバイオマスの差を反映できるシミュレーションツールで十分とも言える。それを叶えるツールとして、BBMapの機能の１つrandomreads.shがある。randomreads.shは配列のシミュレーションを行うスクリプトであるが、snpsやindelをランダムに発生させるオプションを持つ。また、metagenomeフラグをつけて走らせることで菌ごとにカバレッジを対数ランダムにふってリードを発生させることもでき、メタゲノムシミュレーションに最低限必要な機能を満たしている。BBMapを使ってメタゲノムをシミュレートする流れをまとめておく。

インストール

以前紹介した時に書いています。

ラン

catなどでコンカテネートしたfastaを作る。

cat refA.fa refB.fa refC.fa > refABC.fa

150-bpのペアリードを一千万リード発生させる。

randomreads.sh ref=refABC.fa length=150 paired reads=10000000 metagenome

• out　Output file. If blank, filename(s) will be autogenerated.
• ref　Reference file. Not needed if the reference is already indexed.
• reads　Generate this many reads (or pairs).
• paired　Set to true for paired reads.(default f).
• length　Generate reads of exactly this length.(default 100).
• metagenome　Assign scaffolds a random exponential coverage level, to simulate a metagenomic or RNA coverage distribution.
• interleaved　Set to true for interleaved output (rather than in two files) (default f).
• coverage　If positive, generate enough reads to hit this coverage

対数分布なので、十分なカバレッジを与えたと思っていても、カバレッジがほとんどないような クロモソームもできる。

リアルデータとしては、metabatの擬似メタデータが利用できます。"既知"の25サンプルを混ぜ込んでランされた答え合わせが可能なデータです。metabatのマニュアルからダウンロードできるのはbamとアセンブルされたfasta配列です。２ライブラリ分あります。

=>　リンク

リアルデータのメタゲノムが必要なら、例えば

https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP013413g

が利用できます。すでに論文になったデータです。

http://madsalbertsen.github.io/multi-metagenome/

引用

Biostars

https://www.biostars.org/p/235160/

How to generate random short reads from a reference genome [Archive] - SEQanswers