macでインフォマティクス

macでインフォマティクス

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

アラインメントフリーの配列比較GUIツール CAFE

 

 配列比較は、分子配列間の関係を研究するために広く使用されている。配列比較のための主なツールは、global(ref.1)およびlocal(ref.2)シーケンスアラインメントを含むアラインメントベースの方法である。 BLAST(ref.3)のようなアラインメントベースのツールおよびRefSeq(ref.4)のようなシーケンスデータベースの出現により、アラインメントベースの方法は広範囲の用途で広く使用されている。それらの広範な用途にもかかわらず、位置合わせに基づく方法はいくつかの状況において適切ではない。第一に、遺伝子調節領域は一般に高度に保存されておらず、アラインメントに基づくアプローチは類似の転写因子によって結合されている関連調節領域を同定することを困難にする(ref.5)。第二に、次世代シーケンシング(NGS)技術は大量のショートリードを生成するため、ゲノム研究とメタゲノム研究の両方にそれらを組み合わせることは困難である。多くのサンプルにわたって長くアセンブリされたコンティグがなければ、ゲノムとメタゲノムを比較するためのアラインメントベースの方法では困難である(ref.6,7)。第三に、ウイルスは同じようなワードパターンの用法を持つ細菌宿主に感染する可能性が高く(ref.8、9)、したがって、ウイルスの宿主はそれらのワードパターンの用法に基づいて推定される可能性がある。しかしながら、アラインメントに基づく方法は、通常、ウイルス - 宿主感染性関連性の研究には適用できない。 

 アラインメントフリーのシーケンス比較法は、アラインメントに基づく方法が適切でないかまたは実際に実施するには時間がかかりすぎる場合に、配列間の関係を研究するための魅力的な代替法である(ref.10、11)。Briefing in Bioinformaticsの特別号で最近レビューされた(ref.12)、k-mer、最長共通部分列、最短欠如パターンなどに基づくものを含む、いくつかのタイプのアラインメントフリーアプローチが利用可能である。ここでは、k-mer数を使用したアライメントフリー統計に焦点を当てる。これらの手法は各シーケンスを k-mer(または同等に k-tuple、k-gram)カウント特徴空間に射影し、ここでシーケンス情報はk-mer頻度などの数値に変換される。追加の計算上の複雑さがある spaced k-mers を使用する相違測定は考慮しない。(一部略)

 CAFEは、CVTree(ref.15)、d2(ref.13)、d2S(ref.13)など、バックグラウンド調整されたk-merカウントに基づく最近開発されたメジャーの計算を大幅にスピードアップし、同時にメモリ使用量を削減する。さらに、CAFEは、Chebyshev(Ch)、Euclidean(Eu)、Manhattan(Ma)、d2 dissimilarity(ref.16)、Jensen-Shannon divergence(JS)(ref.17)、feature frequency profiles (FFP)(ref.18)、そしてCo-phylog(ref.19)などのk-merカウントに基づく10の従来の尺度を統合する。 CAFEはまた、JaccardやHammingのdistanceなど、k-merの有無に基づいて15の基準を提供する。さらに、霊長類、脊椎動物、微生物のゲノム配列からメタゲノムシーケンシングリードにに至るまで、実際のデータセットでCAFEを使用したアライメントフリーの非類似度測定の価値を示す。

 CAFEは、NGS技術からアセンブリされたゲノム配列とアセンブリされていないショットガンシーケンシングの両方のシーケンスデータを扱い、高速でメモリ効率の良いk-merカウントツールであるJELLYFISH(ref.20)によってk-merをカウントする。 JELLYFISHは、クエリシーケンスを並列に与えられたすべてのk-merカウントを含む圧縮データベースを生成する。その後、CAFEはデータベースをロードし、さまざまな非類似度に関して必要な変換情報を生成する。例えば、k-merの存在/不在に基づく尺度は、k-merカウントを存在/不在指標に二値化する。ほとんどの従来の尺度は、k-merカウントをk-merスペクトラムに正規化する。その上、予想されるk -merカウントは、CVTree、d2およびD2Sのようなバックグラウンド調整されたk-merカウントに基づく最近開発された尺度に関与している。そのような場合、配列のマルコフモデルは、それに応じて配列データから推定されるパラメータを用いて、基礎となる生成モデルとして想定される。マルコフ次数は、手動で設定することも、ベイズ情報量基準(BIC)を使用して自動的に選択することもできる(ref.21)。

 結果として得られるシーケンス間のペアワイズの非類似性は、対称行列を形成する。 CAFEは、非類似度行列を標準のPHYLIP形式で直接出力できる。あるいは、CAFEは、UPGMAアルゴリズムを使用して配列を樹状図にクラスター化すること、マトリックスのヒートマップ可視化、マトリックスtを射影することを含む、4種類の組み込み下流可視化分析を提供する。

 

 

インストール

本体 Github

Githubにあるリンクからダウンロードする。mac版とwindows版が用意されている。

 

 

使い方

解凍して、中にあるCAFEGUIをダブルクリックして起動。

f:id:kazumaxneo:20190719225545p:plain

またはターミナルから叩いて起動。

f:id:kazumaxneo:20190719225801p:plain

 

インターフェイスGithubより)

f:id:kazumaxneo:20190719234645p:plain

 

左上の一番左にあるLoadボタンをクリック

f:id:kazumaxneo:20190719225839p:plain


ここではテストデータの Phylip formatファイルを選択。

f:id:kazumaxneo:20190719225945p:plain

 

Dendrogram

f:id:kazumaxneo:20190719230538p:plain

 

右上のボタンから拡大縮小と図のダウンロードができる(.PNG)。

f:id:kazumaxneo:20190720001351p:plain

 

タブを切り替えることで視覚化方法を比較方法と視覚化方法を切り替えできる

f:id:kazumaxneo:20190720001416p:plain

 

principal coordinate analysis (PCoA)

f:id:kazumaxneo:20190719230605p:plain

 

Heatmap

f:id:kazumaxneo:20190719230623p:plain

 

Network

f:id:kazumaxneo:20190719230643p:plain

 

 

FASTAファイルを読み込んで実行するには、FASTAファイル( .faか .fna)のディレクトリを指定する。左から3番目のAdd all genomeを選択。

f:id:kazumaxneo:20190720000713p:plain

ディレクトリを指定

f:id:kazumaxneo:20190719234957p:plain

またはFASTAファイルを個別に読み込む。左から2番目のAdd one  genomeを選択。

f:id:kazumaxneo:20190720000804p:plain

 

読み込まれたファイルが表示される。

f:id:kazumaxneo:20190720000957p:plain

 

全てのFASTAを読み込んだら、パラメータを決める。

f:id:kazumaxneo:20190720001135p:plain

defaultではマンハッタン距離(Ma)、k-mer=8、k-mer配列の逆相補鎖は考慮しない。

 

右上のRunボタンを押して比較を実行

f:id:kazumaxneo:20190720000427p:plain

 

計算が終わると結果が表示される。

f:id:kazumaxneo:20190720000514p:plain

CUI環境で使うならcafe_mac を叩く。

> CAFE_mac/cafe_mac 

$ CAFE_mac/cafe_mac 

Start parsing the arguments... 

CAFE: aCcelerated Alignment-FrEe sequence analysis

Description: The program provides 29 alignment-free sequence distance measures.

Authors: Yang Lu and Prof. Fengzhu Sun, Computational and Molecular Biology, University of Southern California.

 

usage:

./cafe [options]* -D <dist> -I <fa_files> -K <intK>

 

Main arguments

-D <dist> Comma-separated list of distance measurements, E.g. -D D2star,Ma,CVtree. The options include: 

Conventional measures based on k-mer counts : 

Ch: Chebyshev distance 

Canberra: Canberra distance 

Chisq: Chi-Square distance 

Cosine: Cosine distance 

Co-phylog: Co-phylog distance with the seed C_{(k-1)/2,(k-1)/2}O_{1} when k is odd or C_{k/2-1,k/2}O_{1} when k is even 

D2: D2 distance 

Eu: Euclidean distance 

FFP: Feature frequency profiles (FFP) 

JS: Jensen-Shannon divergence 

Ma: Manhattan distance 

Pearson: Pearson distance 

Newly developed measures based on background adjusted k-mer counts: 

CVtree: CVtree distance 

D2shepp: D2shepp distance 

D2star: D2star distance 

Measures based on presence/absence of k-mers: 

Anderberg: Anderberg distance 

Antidice: anti-Dice distance 

Dice: Dice distance 

Gower: Gower distance 

Hamman: Hamman distance 

Hamming: Hamming distance 

Jaccard: Jaccard distance 

Kulczynski: Kulczynski distance 

Matching: Matching distance 

Ochiai: Ochiai distance 

Phi: Pearson Phi distance 

Russel: Russel-Rao distance 

Sneath: Sneath-Sokal distance 

Tanimoto: Rogers-Tanimoto distance 

Yule: Yule distance 

-I <fa_files> Comma-separated list of sequence fasta files, e.g. -I speciesA.fa,speciesB.fa,speciesC.fa. Pairwise similarity is calculated based upon the sequences specified with this option. 

-K <intK> Kmer Length

 

Options

-J <jfexe_path> Use jellyfish to accelerate kmer counting. <jfexe_path> denotes the file path of jellyfish executable file, e.g. jellyfish-2.2.4/bin/./jellyfish 

-L <lower> Only consider k-mer with occurrence >= <lower>. The default value is 0. 

-M <order> Markov Order involved in D2star and D2shepp. There are two possible options. The first option is one single value indicating that all the sequences use the same order. The second option is comma-separated list of orders. Notice that the length of the list should match the number of fasta files. The order value could be non-negative integer but less than Kmer length or "-1" with the special intention to automatically infer the suitable order (not suitable for JS). The default Markov Order is -1 (Automaticcaly determine by BIC).

-R Consider Reverse Complement in kmer counting. 

-S <dir> Save/Load calculated k-mer count binary files to the folder <dir>. Each input fasta file corresponds to particular model. 

-O <path> Output results to file at <path> 

-T <type> The output type as the input to downstream analysis, including: plain, phylip (as hierarchical clustering), cytoscape (as network analysis) and mds (Multidimensional Scaling as 2D plotting). E.g. -T mds. The default type is plain. 

 

Examples:

./cafe -M 0 -O output_path -S model_dir -T plain -I speciesA.fa,speciesB.fa -J jellyfish-2.2.4/bin/./jellyfish -K 10 -D D2star,Ma

./cafe -M 0 -S model_dir -I speciesA.fa,speciesB.fa -J jellyfish-2.2.4/bin/./jellyfish -K 10 -D D2star,Ma

./cafe -M 0 -L 2 -I speciesA.fa,speciesB.fa -J jellyfish-2.2.4/bin/./jellyfish -K 10 -D D2star,Ma -R


引用
CAFE: aCcelerated Alignment-FrEe sequence analysis
Yang Young Lu, Kujin Tang, Jie Ren, Jed A. Fuhrman, Michael S. Waterman,  Fengzhu Sun

Nucleic Acids Res. 2017 Jul 3; 45(Web Server issue): W554–W559

 

de novo transcriptomeのアセンブリツール TransLiG

 

 オルタナティブスプライシングは真核生物遺伝子における遺伝子調節の重要な形態であり、遺伝子機能の多様性ならびに疾患のリスクを増大させる[ref.1、2、3]。報告されているように[ref.4]、[ref.5]、ヒト遺伝子を含む真核生物遺伝子のほとんどはオルタナティブスプライシングの過程を経るため、異なる細胞1つの遺伝子が数十または数百のスプライシングアイソフォームを作り出せる。したがって、特定の条件下での全長転写物の同定は、その後の多くの生物学的研究において重要な役割を果たす。しかしながら、我々はまだ完全な人間の写しの風景からは程遠い、そしてその状況は非ヒト真核生物種にとってははるかに明らかではない[ref.6]。

 RNA-seqは、これまでにない正確さで、トランスクリプトームレベル全体での発現遺伝子の同定および存在量の測定を可能にする強力な技術である[ref.7、8、9、10]。 RNA-seqプロトコールはサンプリングされた発現トランスクリプトを入力として取り、ランのために2億以上のショートリードを生成し、そして各シーケンシングリードは一般に50-150塩基対の長さであり、完全長トランスクリプトを再構築することは大きな課題である。 第一に、異なる転写物は非常に異なる発現量を有する可能性があり、それは構築されたシーケンスグラフ(スプライシンググラフ、De brujinグラフなど)をかなり不均一な範囲にする。第二に、同じ遺伝子からの異なる転写物が選択的スプライシングのためにエキソン配列を共有し、スプライシンググラフをさらに複雑にする。第三に、大量のRNA-seqリードにはシークエンシングエラーが含まれているため、RNA-seqデータから低発現の転写産物を集めることはさらに困難である。上記のすべてが、トランスクリプトームアセンブリの問題を非常に困難にしている。

 近年、転写産物アセンブリの問題を解決するために開発された方法が増えてきており、それらの大部分は2つのアプローチに分類することができる:リファレンスベース(またはゲノムガイド)およびde novo [ref.11、12]。 Scallop [ref.13]、TransComb [ref.14]、StringTie [ref.6]、Cufflinks [ref.15]、およびScripture [ref.16]などのリファレンスベースのアプローチでは、通常、Hisa[ref.17]、Star [ref.18]、Tophat [ref.19]、SpliceMap [ref.20]、MapSplice[ref.21]、またはGSNAP [ref.22]などのアライメントツールを使用してRNAシーケンスのリードをリファレンスゲノムにマッピングする。 ならびに同じ遺伝子座からのリードは、スプライシンググラフを形成するためにクラスタに分類される。全ての発現された転写物はグラフを横断することによって集めることができ、この戦略によってアセンブリされた転写物は、リファレンスゲノムから恩恵を受けるので、一般にデノボの戦略によるものと比較してより高い精度を有する。しかしそのような高品質のリファレンスゲノムは現在ほとんどの種で利用できないため、実際には著しく制限される。

 リファレンスゲノムが利用できない、不完全である、高度に断片化されている、または癌組織におけるように実質的に改変されている場合、新規アセンブリは望ましいアプローチである。 BinPacker [ref.23]、Bridger [ref.24]、Trinity [ref.12]、IDBA-Tran [ref.25]、SOAPdenovo-trans [ref.26]、ABySS [ref.27]、Oases [ref.28]など、数多くのde novoアセンブラがある。この戦略は通常、それらの配列の重複に基づいてRNA-seqリードからスプライシンググラフを直接構築し、そして次に異なるアルゴリズムを使用してグラフをトラバースすることによって転写物を集める。 IDBA-Tran、SOAPdenovo-trans、ABySS、およびOasesなどのアセンブラは、ゲノム構築における重要な技法に基づいて開発されているため、一般に、トランスクリプトーム構築ではうまく機能しない。 Trinityは、トランスクリプトームのde novoアセンブリを処理するための方法を設計するための扉を開く。最初にk-mer拡張ストラテジーによってシーケンシングリードを長いコンティグに拡張し、次にそれらのコンティグをde bruijnグラフに接続し、最後にde bruijnグラフをトラバースすることによってすべての表現された転写産物を推論する。 Trinityの論文に見られるように、そのアプリケーションを妨げるいくつかの制限がある。アセンブリ手順において有用であるであろうデプス情報は適切に使用されず、そしてブルートフォース戦略がde bruijnグラフにおける転写産物を表す経路を検索するために適用され、それが偽陽性率としてユーザーをひどく悩ます。 Bridgerは、リファレンスベースのアセンブラCufflinksからde novoアセンブリへの最小パスカバーモデルの移植に成功し、徹底的な列挙を効果的に回避し、誤検知を大幅に減少させる。ただし、Trinityの論文に記載されているように、アセンブリ手順の開発に役立つはずのシーケンスデプス情報を十分に活用していない。続いて、最小パス数を制限することなく、ビンパッキングモデルによってシーケンスデプス情報を完全に使用するために、新しいアセンブラBinPackerが開発された。 BinPackerは同種の他のものよりも優れたパフォーマンスを発揮するが、ペアエンド情報アセンブリ手順に統合していないため、改善の余地がある。

 本稿では、フェージングパスを用いて開発された新しいde novoアセンブラTransLiGと、スプライシンググラフから繰り返し重み付き折れ線グラフを構築する方法を紹介する。 TransLiGでのフェージングパスのアイデアは、Scallop [ref.13](リファレンスベースのトランスクリプトームアセンブラ)から動機付けられた。(一部略)

 Transity、Bridger、およびBInPackerなどの同種のすべての顕著なツールよりも実質的に優れたものにするために、TransPathおよび反復折れ線グラフを構築することによってシーケンスデプスおよびペアエンド情報をアセンブリ手順に統合するTransLiGを開発した。人工データと実データの両方でテストした場合、テストされたマウスのデータにそれぞれで、精度はBinPackerとBridgerよりも6%高く、人工データではTrinityよりも15%近く、BinPacker、BridgerとBridgerよりも7%、14%、および21%高かった。 TransLiGは最高の精度を達成するだけでなく、テストしたすべてのデータセットで最高の感度を達成した。さらに、TransLiGはさまざまな評価パラメータで安定して最高のパフォーマンスを維持する。(以下略)

 

 

インストール

依存

  • boost_1_47_0 
#boost導入
wget http://sourceforge.net/projects/boost/files/boost/1.47.0/boost_1_47_0.tar.gz
tar zxvf boost_1_47_0.tar.gz
cd boost_1_47_0/

#ビルド 指定なしだと/usr/local/に導入されるが、ディストリビューション提供のboostが入っていると優先順でおかしくなる。
#できるだけ別の場所に入れる。例えば"$HOME/local"等
./bootstrap.sh
./b2 install --prefix=</BOOST/install_dir>

本体  SourceForge

cd TransLiG_1.0
./configure --with-boost=</BOOST/install_dir>
make -j 12

> ./TransLiG 

# ./TransLiG 

    

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

    

root@d6b44964ef68:/data/TransLiG_1.1# ./TransLiG -h

    

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

    

TransLiG v1.0 usage:

    

** Required **

    

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

    

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

    

If paired_end reads:

   -l <string>: left reads.

   -r <string>: right reads.

    

If single_end reads:

   -u <string>: single reads.

    

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

    

** Options **

    

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

    

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

             if paired_end: RF or FR;

             if single_end: F or R.

    

-t <int>: minimum length of transcripts, default: 200.

    

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

    

-K <int>: minimum length of kmer used to connect fragmented graphs, default: 21.

    

-c <int>: minimum coverage of nodes used to connect fragmented graphs, default: 30.

    

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

    

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

    

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

    

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

    

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

    

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

    

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

    

** Note **

    

A typical command of TransLiG might be:

    

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

    

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

    

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

root@d6b44964ef68:/data/TransLiG_1.1# 

 

 plugin/fastool/のビルドができておらずエラーが出た。別途ビルド。

cd plugins/fastool/ 
make
./fastool

  

実行方法

ペアエンドfastqを指定する。

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

ランの途中、fastoolによるfastq処理後にエラーが出る。

 

引用

TransLiG: a de novo transcriptome assembler that uses line graph iteration

Juntao Liu, Ting Yu, Zengchao Mu, Guojun Li
Genome Biology 2019 20:81

 

 

参考


 

アンプリコンシーケンスのペアエンドリードマージツール MeFiT

 

 

 次世代シークエンシング技術は,その開始以来,研究者が複雑なシステムから多面的な生物学的情報を抽出する方法を変え、ヒト疾患,環境科学、進化科学などの分野における研究を促進してきた。16S rRNA小サブユニット遺伝子、またはより一般的にはその一部のシークエンシングによる細菌生態学の探究は、様々なヒト/環境マイクロバイオーム研究[ref.2、3、4、5、6]のゴールドスタンダード技術として使用されている。このように、原核生物の16S rRNA遺伝子の異なる部分の予測可能な保存と変動性は、細菌ソースの高分解能同定と定量情報を提供するために利用されてきた [ref.6]。系統発生分析および定量のためのこの戦略は、従来のクローニングおよびシーケンシングまたはRT-PCRに基づくアプローチよりも有意に効率的であることが証明されている。これらの大規模な並列処理が可能で、費用対効果が高く、処理能力の高いシークエンシング技術は、現在では1回の実行で最大15Gbのゲノムデータを生成することができる [ref.7]。これらのNGSシステムによって生成されるデータの膨大な量および複雑さは、下流の分析を合理化するバイオインフォマティクスツールの開発を必要とする。

(一段落省略)

 ペアードエンドシーケンシングでは、正方向および逆方向リードの正確なマージは、特に微生物分類学的プロファイリングを含むが、これに限定されず多数の下流分析の結果に影響を及ぼす重要な第一段階である。各種ツール;たとえば、SHERA [ref.8]、FLASH [ref.9]、PANDAseq[ref.10]、COPE[ref.11]などは、ペアエンドデータをマージするために開発された。これらのツールは一般に、ペアのリードの末端のリード間の最良のオーバーラップを識別するためにシークエンスアラインメントを行う。クオリテイィスコアを考慮することでミスマッチ基準を解決し、単にクオリティの高い方にクオリティの低い方を置き換える。最近(論文執筆時点)、新しい方法CASPERが提案された[ref.12]。CASPERはクオリティスコアの差が有意でない場合を除いて,ミスマッチを解決するために伝統的なクオリティスコアに基づく方法を使用する。(一部省略)

 16S rRNA微生物プロファイリングの精度はシーケンシングデータの前処理精度に大きく依存する。そこで,Illumina MiSeqプラットフォームからのオーバーラップするペアエンドリードを効率的にマージし、それらをクオリティフィルタリングするMeFiTと名付けたツールを開発した。MeFiTは、ペアエンドリードをマージするためにCASPER[ref.12]を呼び出し、注意深いクオリティフィルタリングを含め、拡張する。

(以下略) 

 

f:id:kazumaxneo:20190523235426p:plain

The MeFiT pipeline. 論文より転載

 

インストール 

ubuntu16.0.4のminiconda2.4.0.5環境でテストした。

依存

  • Python version 2.7 (version 3.0 or greater may not work)
  • numpy 
  • HTSeq
  • CASPER (Context-Aware Scheme for Paired-End Read)
  • Jellyfish mer counter 
conda install -y numpy
conda install -c bioconda -y HTSeq
conda install -c bioconda -y bioconductor-casper
conda install -c bioconda -y jellyfish #(version2系にする)

casperが入らなければ直接入れてください(一番下のリンク参照)。

本体 Github

git clone https://github.com/nisheth/MeFiT.git
cd MeFiT/

>./mefit -h

$ ./mefit -h

usage: mefit [-h] -s S -r1 R1 -r2 R2 [-p P] [-nonovlp] [-n N]

             (-avgq AVGQ | -meep MEEP)

 

MeFiT - developed by Hardik I. Parikh, PhD

MeFiT - Merging and Filtering Tool for paired-end reads

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

 

For detailed information about the command - 

mefit -h

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

 

optional arguments:

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

  -s S        Sample Name

  -r1 R1      Forward Read - Sample R1 fastq

  -r2 R2      Reverse Read - Sample R2 fastq

  -p P        CASPER parameter file, tab-delimited

  -nonovlp    Save non-overlapping reads, default=False

  -n N        Length of Ns to insert between non-overlapping reads for

              merging, default=15

  -avgq AVGQ  Average-Q threshold, default=20

  -meep MEEP  meep-score threshold, default=1.0

 

実行方法

 ペアエンドfastqを指定する。

mefit -s sample -r1 sample_R1.fastq -r2 sample_R2.fastq -avgq 20
  • -avgq    Average-Q threshold, default=20
  • -r1   Forward Read - Sample R1 fastq
  • -r2   Reverse Read - Sample R2 fastq
  • -meep    meep-score threshold, default=1.0

 

オーバーラップしていないペアエンドもシュードマージする。N=15-bpで連結。

mefit -s sample -r1 sample_R1.fastq -r2 sample_R2.fastq \
-avgq 20 -nonovlp -n 15

 

引用
MeFiT: merging and filtering tool for illumina paired-end reads for 16S rRNA amplicon sequencing
Parikh HI, Koparde VN, Bradley SP, Buck GA, Sheth NU

BMC Bioinformatics. 2016 Dec 1;17(1):491

 

関連


アセンブリ過程でロングリードをフィルタリングする fpa

 

以前、ロングリードのアセンブリ前処理ツール yacrdを紹介した。

今回はアセンブリ過程でフィルタリングして出力を調節するfpaを紹介する。

 

以下のフィルタリングが行える (Githubより)。

  • internal match
  • containment
  • dovetails
  • self matching
  • read name match against regex
  • length of overlap
  • length of read in overlap

 

Pierre MarijonさんのBlog (More details and some experiment are present)

https://blog.pierre.marijon.fr/binary-mapping-format/

 

インストール

ubuntu16.0.4のminicona3.4.0.5環境でテストした。

ビルド依存

  • Rust >= 1.32
  • libgz
  • libbzip2
  • liblzma

その他 (option)

  • minimap2
  • miniasm

本体 Github

#bioconda (link) (*1)
conda install -c bioconda -y fpa

> fpa -h

# fpa -h

fpa 0.5 Krabby

Pierre Marijon <pierre.marijon@inria.fr>

fpa take long read mapping information and filter them

 

USAGE:

    fpa [OPTIONS] [SUBCOMMAND]

 

FLAGS:

    -h, --help       Prints help information

    -V, --version    Prints version information

 

OPTIONS:

    -z, --compression-out <compression-out>

            Output compression format, the input compression format is chosen by default [possible values: gzip, bzip2,

            lzma, no]

    -F, --format <format>                                  Force the format used [possible values: paf, m4]

    -i, --input <input>                                    Path to input file, use '-' for stdin [default: -]

        --internal-threshold <internal-match-threshold>

            A match is internal match if overhang length > match length * internal threshold this option set internal

            match [default: 0.8]

    -o, --output <output>                                  Path to output file, use '-' for stdout [default: -]

 

SUBCOMMANDS:

    drop      fpa drop mapping match this constraints

    gfa       fpa generate a overlap graph in gfa1 format with mapping passing filter

    help      Prints this message or the help of the given subcommand(s)

    index     fpa generate a index of mapping passing filter

    keep      fpa keep only mapping match this constraints

    rename    fpa rename reads with name you chose or with incremental counter

> fpa keep -h

# fpa keep -h

fpa-keep

fpa keep only mapping match this constraints

 

USAGE:

    fpa keep [FLAGS] [OPTIONS] [SUBCOMMAND]

 

FLAGS:

    -c, --containment      Keep only containment mapping

    -d, --dovetail         Keep only dovetail mapping

    -h, --help             Prints help information

    -i, --internalmatch    Keep only internal mapping

    -m, --same-name        Keep only mapping where reads have same name

    -V, --version          Prints version information

 

OPTIONS:

    -l, --length-lower <length_lower>                      Keep only mapping with length lower than value

    -L, --length-upper <length_upper>                      Keep only mapping with length upper than value

    -n, --name-match <name_match>                          Keep only mapping where one reads match with regex

    -s, --sequence-length-lower <sequence_length_lower>

            Keep only mapping where one reads have length lower than value

 

    -S, --sequence-length-upper <sequence_length_upper>

            Keep only mapping where one reads have length upper than value

fpa drop -h

# fpa drop -h

fpa-drop 

fpa drop mapping match this constraints

 

USAGE:

    fpa drop [FLAGS] [OPTIONS] [SUBCOMMAND]

 

FLAGS:

    -c, --containment      Drop containment mapping

    -d, --dovetail         Drop dovetail mapping

    -h, --help             Prints help information

    -i, --internalmatch    Drop internal mapping

    -m, --same-name        Drop mapping where reads have same name

    -V, --version          Prints version information

 

OPTIONS:

    -l, --length-lower <length_lower>                      Drop mapping with length lower than value

    -L, --length-upper <length_upper>                      Drop mapping with length upper than value

    -n, --name-match <name_match>                          Drop mapping where one reads match with regex

    -s, --sequence-length-lower <sequence_length_lower>    Drop mapping where one reads have length lower than value

    -S, --sequence-length-upper <sequence_length_upper>    Drop mapping where one reads have length upper than value

fpa gfa -h

# fpa gfa -h

fpa-gfa 

fpa generate a overlap graph in gfa1 format with mapping passing filter

 

USAGE:

    fpa gfa [FLAGS] [OPTIONS] [SUBCOMMAND]

 

FLAGS:

    -c, --containment      Keep containment overlap

    -h, --help             Prints help information

    -i, --internalmatch    Keep internal match overlap

    -V, --version          Prints version information

 

OPTIONS:

    -o, --output <output>    Write mapping passing filter in gfa1 graph format in path passed as parameter

fpa index -h

# fpa index -h

fpa-index 

fpa generate a index of mapping passing filter

 

USAGE:

    fpa index [OPTIONS] [SUBCOMMAND]

 

FLAGS:

    -h, --help       Prints help information

    -V, --version    Prints version information

 

OPTIONS:

    -f, --filename <filename>    Write index of mapping passing filter in path passed as parameter

    -t, --type <type>            Type of index, only reference read when it's query, target or both of them [default:

                                 both]  [possible values: query, target, both]

> fpa rename -h

# fpa rename -h

fpa-rename 

fpa rename reads with name you chose or with incremental counter

 

USAGE:

    fpa rename [OPTIONS] [SUBCOMMAND]

 

FLAGS:

    -h, --help       Prints help information

    -V, --version    Prints version information

 

OPTIONS:

    -i, --input <input>      Rename reads with value in path passed as parameter

    -o, --output <output>    Write rename table in path passed as parameter

 

実行方法

1、 keep - 制約に合致するマッピングを出力

2、 drop - 制約に合致するマッピングを除去

minimap/miniasmのロングリードアセンブリの途中にfpaのフィルタリングを加えるように設計されている。そのため、fpaはpafファイルを入力とし、pafファイルを出力する。miniasmによるRaw assemblyの第一ステップはminimapまたはminimap2によるリード同士のオーバーラップ(all-vs-all read self-mappings)なので、この第一ステップの出力(minimapの出力はdefaultではPAF)を受け取り、fpaでフィルタリングしてminiasmでアセンブリ(コンカテネートしてunitigをGFA出力(エラーレートはinputのリードと同じ))を行う。

例えばONTのrawアセンブリを以下の流れで実行するとする。
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | gzip -1 > reads.paf.gz
miniasm -f nanopore.fq.gz reads.paf.gz > raw_assebly.gfa


#この流れにfpaを組み込む。"fpa drop -l <NUM>"で指定以下のリードを捨てて出力。
#1000bp以上のリードだけpaf出力。
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | fpa drop -l 1000 | gzip - > filtered.paf.gz
#それからminiasmを実行する

#fpa dropなら逆の条件でフィルタリングできる。fpa keep -l <NUM>で指定以下のリードのみキープして出力
#10000bp以下のリードだけpaf出力。
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | fpa keep -l 10000 | gzip - > filtered.paf.gz

#つまりkeepとdropは逆の挙動をする。1000bp以上のリードだけキープして出力するなら"keep -L <NUM>"か"drop -l <NUM>"。keepを使うなら
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | fpa keep -L 1000 | gzip - > filtered.paf.gz

 

サイズセレクト

 1000bp以上10000-bp以下のリードのみにフィルタリング。

#keepを使うなら
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | fpa keep -L 1000 -l 10000| gzip - > filtered.paf.gz

#dropを使うなら
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | fpa drop -l 1000 -L 10000| gzip - > filtered.paf.gz

 

dovetail(蟻継ぎ)のオーバーラップを持つリードのみにフィルタリング。

minimap2 long_read.fasta long_read.fasta | fpa keep -d | gzip - > filtered.paf.gz
  • -d     Keep only dovetail mapping
       

インターナルのオーバーラップを持つリードのみにフィルタリング。

minimap2 long_read.fasta long_read.fasta | fpa keep -i | gzip - > filtered.paf.gz
  • -i     Keep only internal mapping 

 

汚染配列を除去。

minimap2 long_read.fasta long_read.fasta | fpa keep -c | gzip - > filtered.paf.gz
  • -c   Keep only containment mapping

 

汚染配列除去とサイズセレクト。

minimap2 long_read.fasta long_read.fasta | fpa keep -c keep -L 1000 | gzip - > filtered.paf.gz
  • -c   Keep only containment mapping

 

 

他の例はGithubを参考にしてください。 

 

fpa gfaコマンドも組み込めば、アセンブリグラフのGFA(version1)で出力できる。bandage等を使えばグラフを可視化できる。

minimap2 long_read.fasta long_read.fasta | fpa keep -d gfa -o dovetail.gfa > dovetail.paf

あくまでアセンブリ前提に設計されているので、GFAの出力(-o)以外にリダイレクトでPAFも出力する。fpa gfa -oで出力したGFAに、配列情報は含まれないので注意する。

 

 

fpaはロングリードのアセンブリを前提に設計されているので、アセンブリ前にロングリードを前処理したいならyacrdを使う。

引用

yacrd and fpa: upstream tools for long-read genome assembly

Pierre Marijon, Rayan Chikhi, Jean-Stéphane Varré 

bioRxiv preprint first posted online Jun. 18, 2019

 


 

*1

2019 7/4現在、macosだと0.4までしか入らない。v0.4だとkeepなどの一部コマンドが存在しない。 keepなどを使いたければソースからビルドする。

 

 

PGAPとPGAP-Xを組み込んだバクテリアのパンゲノム解析webサーバー PGAweb

2019 7/21追記

 

"PGAP-X: extension on pan-genome analysis pipeline"より

 パンゲノムの概念は2005年に提案されて以来[ref.12]、過去10年間でバクテリアゲノムの進化と動態を調査するために急速に採用されてきた[ref.3、6]。最近では、ウイルス[7]、真菌[ref.8]、植物[ref.9]の比較ゲノム解析にも広く利用されている。バクテリアゲノムのパンゲノム解析をより簡単かつ効率的にするために、Panseq [ref.10]、PGAT [ref.11]、PanCGHweb [ref.12]、PanGP [ref.13]、ITEP [ref.14]、PGAP [ref.15]を含むいくつかのプログラムとデータベースが開発されている。初期のプログラムやデータベースは主に限られた機能解析に焦点を合わせていたが、PGAPは機能遺伝子のクラスター分析、パンゲノムプロファイル分析、機能遺伝子の変異分析、種進化分析および遺伝子クラスターの機能強化分析を含む5つの共通分析モジュールを統合する。公開後、PGAPは60カ国以上から4000回以上ダウンロードされ、Mycobacterium [ref.16]、Bifidobacterium [ref.17]、Lactococcus [ref.18]などのさまざまな細菌のパンゲノム解析に広く使用されてきた。著者らの知る限り、PGAPはPanseqと共に、2014年末に最も人気のあるパッケージと報告された[ref.19]。

 しかし、パンゲノムツールの終わりのない改善は、データの解釈と視覚化である。これは、より良いデータマイニングの結果と、研究とpublicationのための質の高いグラフィックを提供する。過去数年間で、パンゲノムサイトからのデータを視覚化するために、いくつかのスタンドアロンプ​​ログラムとWebベースのサーバーが開発されてきた。しかし、これらのプログラムとサーバーは非常に限られた機能しか提供していない。さらに、それらはゲノム構造の観点からゲノムおよび遺伝子の両方の内部でオルソログ関係および遺伝的変異を示すことができない。この問題に取り組むために、著者らはゲノム指向のソフトウェア、PGAP-Xを開発した。これはゲノム構造の視界か​​らパンゲノム分析を実行する。 PGAP-Xは独立してデータ分析を実行するだけでなく、結果データを直接視覚化して解釈する。同様の機能を持つ他のプログラムによって生成された結果データは、互換性のあるデータ形式に変換された後、さらに分析および視覚化するためにPGAP-Xにインポートすることもできる。 PGAP-Xは、ゲノム構造において高い類似性を有する、同じ種またはclosely relatedな種からのそれらの株についてのゲノム構造および遺伝子含有量の多様性を十分に分析および提示するために使用され得る。

 PGAP-Xでは、分析プロセスは論理的に3つの層に分割されている(論文図1)。データインタフェース層、データ分析層、およびデータ視覚化層である。ユーザーは、データインターフェース層を使用して、パラメータを簡単にカスタマイズし、入力および出力データを管理できる。計算はデータ分析層を介して実行され、結果データはデータ視覚化層によって視覚化される。データ分析および可視化のためのすべてのモジュールは、それらの機能として4つの部分に編成することができる。1)全ゲノム配列アラインメントおよびゲノム構造の可視化。 2)オルソロガス遺伝子クラスタリングおよび保存による遺伝子分布の可視化。 3)全ゲノムプロファイル分析および全ゲノムプロファイル曲線の可視化。 4)遺伝子およびゲノム規模の両方からの変異分析および可視化。 

( 以下略)

 

”PGAweb: A Web Server for Bacterial Pan-Genome Analysis”より

 近年の微生物ゲノムデータの天文学的な増加は、種内および種間の全ゲノム分析のためのバイオインフォマティックツールに対する強い需要をもたらしている。この論文では、2つの主要なパンゲノム解析モジュール、PGAPとPGAP-Xを組み込んだ、バクテリアのパンゲノム解析のためのユーザーフレンドリーなWebベースのツールPGAwebを紹介する。 PGAwebは、オルソロガスクラスタリング、パンゲノムプロファイリング、変異と進化解析、および機能分類を含む、インタラクティブでカスタマイズ可能な主要機能を提供する。 PGAwebは、バクテリアゲノムのダイナミクスと進化を直感的に理解するのに役立つ、さまざまな視覚化手法を用いたゲノム構造ダイナミクスと配列多様性の特徴を提示する。 PGAwebはパラメータをワンクリックで設定できる直感的なインターフェースを持ち、http://PGAweb.vlcc.cn/で無料で利用できる。

 

PGAP-X HP

f:id:kazumaxneo:20190715210947p:plain

ソースコードコンパイルされたプログラムをダウンロードできる。

 

PGAwebマニュアル

Document

 


 

パンゲノムとは何かについては2005年のペーパーを読んでください。言葉の定義から説明されています。

https://www.sciencedirect.com/science/article/pii/S0959437X05001759?via%3Dihub

 

Pan-genomicなアプローチでバクテリアの種について研究した報告はいくつもあります。ここでは以下の論文リンクを載せておきます。

https://www.ncbi.nlm.nih.gov/pubmed/29802996

 

植物の例: イネ

https://www.nature.com/articles/s41586-018-0063-9

 

 

入力ファイル 

exampleデータリンク

1、PGAP-Xの場合リンク

入力ファイル形式はゲノム全長のFASTAファイルとそのテーブルファイルになる。テーブルファイルはNCBI genebankから変換するスクリプトが公開されているので、それを使えば良い。入力ファイルの拡張子以外は同じ名前で統一、suffixはFASTAが".fna"、テーブルファイルは".ptt"にする。すなわち、

GCF_000007045.1_ASM704v1_cds_from_genomic.fna

GCF_000007045.1_ASM704v1_cds_from_genomic.ptt

のようなファイルを準備する。

genebank (new format)からpttに変換するスクリプトがあるので(*1より)、それを使って変換する。BioperlのSeqIOが必要(参考 )。

  • create ptt from genebank format file

 

#!/usr/bin/env perl
use strict;
use Bio::SeqIO;

# This script takes a GenBank file as input, and produces a
# NCBI PTT file (protein table) as output. A PTT file is
# a line based, tab separated format with fixed column types.
#
# Written by Torsten Seemann
# 18 September 2006

my $gbk = Bio::SeqIO->new(-fh=>\*STDIN, -format=>'genbank');
my $seq = $gbk->next_seq;
my @cds = grep { $_->primary_tag eq 'CDS' } $seq->get_SeqFeatures;

print $seq->description, " - 0..",$seq->length,"\r\n";
print scalar(@cds)," proteins\r\n";
print join("\t", qw(Location Strand Length PID Gene Synonym Code COG 
Product)),"\r\n";

for my $f (@cds) {
   my $gi = '-';
   $gi = $1 if tag($f, 'db_xref') =~ m/\bGI:(\d+)\b/;
   my $cog = '-';
   $cog = $1 if tag($f, 'product') =~ m/^(COG\S+)/;
   my @col = (
     $f->start.'..'.$f->end,
     $f->strand >= 0 ? '+' : '-',
     ($f->length/3)-1,
     $gi,
     tag($f, 'gene'),
     tag($f, 'locus_tag'),
     $cog,
     tag($f, 'product'),
   );
   print join("\t", @col), "\r\n";
}

sub tag {
   my($f, $tag) = @_;
   return '-' unless $f->has_tag($tag);
   return join(' ', $f->get_tag_values($tag)) 

 

 ここではgenebank2ptt.plとして保存。

genebankからpttに変換。

#GCF_000007465.2_ASM746v2_genomic.gbffから変換
perl genebank2ptt.pl <GCF_000007465.2_ASM746v2_genomic.gbff > GCF_000007465.2_ASM746v2_genomic.ptt

たくさんあるならシェルを使いループで回す。 

 

genebak => ptt変換にはこちらのコードも利用できる(こちらもperlの標準ライブラリ以外にSeqIOを使用)。

git clone https://github.com/sgivan/gb2ptt.git
perl gb2ptt/bin/convert.pl -h

$ perl gb2ptt.pl -h

 

--debug

--verbose

--help

--infile

--rast (use with GenBank files downloaded from RAST)

 

genebakからpttに変換。

perl gb2ptt.pl --infile input.gbff

 input.gbff.pttが出力される。

  

2、PGAPの場合  (リンク

fullのgenebankファイル(.gb)か、ゲノム全長のFASTAファイル(.fna)とテーブルファイル(.ptt)、プロテインファイル(.faa)を3つセットでアップロードする。後者はgzip圧縮されていても良い。他の形式にも対応する。詳しくはマニュアル参照。

 

 

 PGAwebの使い方

http://pgaweb.vlcc.cn にアクセスする。 PGAPかPGAP-Xを選択する。 

f:id:kazumaxneo:20190715171623j:plain

PGAPはドラフトゲノムにも完全長ゲノムにも対応するが、 PGAP-Xは完全長ゲノムでみ正しく動作する。 PGAPよりPGAP-Xの方がゲノム比較と視覚化機能が強化されているが、オルソロガスクラスター探索やSNPsツリーなど、PGAPにしかない機能もある(わかりやすく言えば、分析結果の図が欲しいならPGAP-X、分析結果の表が欲しいならPGAP)。

 

ここではPGAP-Xを選択した。

f:id:kazumaxneo:20190715171708j:plain

 

Browseボタンを押してファイルをアップロードする。f:id:kazumaxneo:20190715172119j:plain

 

.fnaファイルと.pttファイルを両方選択してアップロードする(.fnaファイルだけアップロードすれば.同じディレクトリにあるpttを自動認識する訳ではない)。

f:id:kazumaxneo:20190715171804j:plain

拡張子以外のファイル名が同じなら同じゲノム由来と認識するので、同時にアップロードせず、fnaファイルと.pttファイルを順次アップロードしてもO.K。

 

ファイル選択後、左のUploadボタンを押してアップロード開始。

f:id:kazumaxneo:20190715172119j:plain

 

アップロードが終わると右側にアップロードされたファイルが並ぶ。

f:id:kazumaxneo:20190715171806j:plain

fnaとptt両方問題ないゲノムは背景が緑色になる。3ゲノム以上になると、次に進めるというメッセージがでる。

 

解析内容とパラメータを決定したらRUNする。

f:id:kazumaxneo:20190715172443j:plain

 

データ数に応じ、解析にはしばらく時間がかかる。

 

f:id:kazumaxneo:20190715172313j:plain

  

1、Genome Alignment
ゲノムアライメントの結果に基づいて、ゲノム構造が可視化される。 相同なDNA領域は同じ色でマークされる。 f:id:kazumaxneo:20190715221957j:plain

 

2、Orthologs analysis
オルソログ解析結果に基づいて、ゲノム上の遺伝子分布が視覚化される。 同じオルソログ遺伝子は同じ色でマークされる。 

f:id:kazumaxneo:20190715221959j:plain

 

3、Genetic Variation Analysis
全ゲノムアラインメント結果に基づいて各ゲノムの全変異部位が検出され視覚化される。 置換頻度がフィルタ条件以上である1つの領域または領域内のm個(mは変異数を表す)の置換部位が、高置換領域として識別される。 ペアワイズゲノム間のすべての変異部位が検出され出力テキストファイルに報告されるが、選択された系統の高置換領域のみが表示される。

f:id:kazumaxneo:20190715222047j:plain

 

4、Pan-genome Analysis

ゲノム解析モジュールでは、全ゲノムサイズとコア遺伝子サイズの曲線が同じグラフに表示される。

f:id:kazumaxneo:20190715222053j:plain 

 

引用

PGAP-X: extension on pan-genome analysis pipeline
Zhao Y, Sun C, Zhao D, Zhang Y, You Y, Jia X, Yang J, Wang L, Wang J, Fu H, Kang Y, Chen F, Yu J, Wu J, Xiao J

BMC Genomics. 2018 Jan 19;19(Suppl 1):36

 

PGAweb: A Web Server for Bacterial Pan-Genome Analysis
Xinyu Chen, Yadong Zhang, Zhewen Zhang, Yongbing Zhao, Chen Sun, Ming Yang, Jinyue Wang, Qian Liu, Baohua Zhang, Meili Chen, Jun Yu, Jiayan Wu, Zhong Jin, Jingfa Xiao

Front Microbiol. 2018; 9: 1910.Published online 2018 Aug 21

 

*1

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

 

関連


 

マルチプルアラインメント結果からコンセンサス配列を出力するEMBOSSのconsコマンド

 

タイトルの通りのコマンド。

 

HP

EMBOSS: cons

 

インストール

macos10.12の miniconda3-4.3.21環境でテストした。

condaやbrewで導入できる。

#bioconda (link) 
conda install -c bioconda -y emboss

#homebrew
brew install emboss

seqret -h

$ seqret -h

Read and write (return) sequences

Version: EMBOSS:6.6.0.0

 

   Standard (Mandatory) qualifiers:

  [-sequence]          seqall     (Gapped) sequence(s) filename and optional

                                  format, or reference (input USA)

  [-outseq]            seqoutall  [<sequence>.<format>] Sequence set(s)

                                  filename and optional format (output USA)

 

   Additional (Optional) qualifiers: (none)

   Advanced (Unprompted) qualifiers:

   -feature            boolean    Use feature information

   -firstonly          boolean    [N] Read one sequence and stop

 

   General qualifiers:

   -help               boolean    Report command line options and exit. More

                                  information on associated and general

                                  qualifiers can be found with -help -verbose

 

 

実行方法

たえばmafft紹介)でマルチプルシーケンスアラインメントを行い、出力からコンセンサス配列を得る。consと打って実行。

cons

入力のFASTAと出力のFASTA名を順番に入力する。

$ cons

Create a consensus sequence from a multiple alignment

Input (aligned) sequence set: 

 

またはinputとoutputのfasta名を指定して実行。

cons inut.maf output.fasta

 

引用

EMBOSS: the European Molecular Biology Open Software Suite.
Rice P, Longden I, Bleasby A

Trends Genet. 2000 Jun;16(6):276-7.

 

関連


マッピングベースのメタゲノム存在量プロファイリングを行う MiCoP

 

 微生物は、土壌、海水、人体など、地球上のほとんどすべての生態系に遍在している。単細胞生物はこれらの環境のそれぞれにおいて多くの重要な役割を果たしている[ref.1、2]。サンプル中に存在する微生物を特定することは、これらの生物によってどのような機能が実行されているのかを理解し、微生物群集の乱れがどのように様々な疾患につながるかを特徴付けるために重要である。伝統的に、微生物は培養ベースの技術によって研究されてきた。そこで、微生物は単離され、そして実験室環境で個々に研究されてきた。しかし、多くの微生物は培養できないことはよく知られている。それゆえ、それらは実験室環境で研究することはできない[ref.3]。さらに、実験室環境で微生物を研究する技術は、それらの自然の生息地における数百から数千の異なる微生物種間の複雑な関係を捉えることができない[ref.1、2]。ハイスループットシークエンシングは、微生物環境研究に革命をもたらし、宿主環境から直接何千もの微生物ゲノムの研究を可能にし、メタゲノミクスの分野を形成した。このアプローチは、伝統的な培養に依存するバイアスを回避し、人のさまざまな組織や自然の生息地における微生物群集の構成の研究を可能にする[ref.1、2]。メタゲノムプロファイリングは、16SリボソームRNA遺伝子シーケンシングを用いて、以前はバイアスのない方法で研究が不可能であった、真核生物およびウイルス病原体を含む様々な微生物を研究するために有用であることが証明された[ref.4、5、6]。

 マイクロバイオームと人の健康に影響を与える「 virome」と「eukaryome」の決定的な重要性にもかかわらず、ほとんどのメタゲノムプロファイリング方法は主にバクテリア古細菌の同定に焦点を合わせている[ref.7、8]。メタゲノムプロファイリングのためのいくつかの既存の方法は、特定の種に由来するものとしてリードをユニークに識別する「マーカー遺伝子」を使用することを提案した。この方法は、サンプル中のバクテリア古細菌の存在と相対的存在量を推定するのに効率的で正確であることが示されている[ref.9、10、11]。しかしながら、マーカー遺伝子に基づくアプローチは、ウイルスおよび真核生物ゲノムの同定に関していくつかの制限を有する。一つのアプローチは、「普遍的」と見なされるが種間で異なる遺伝子の違いを比較することを含む。このアプローチでは、そのマーカー遺伝子の特定の配列を示し、したがって種を一意に識別するリードを使用する[ref.9]。しかしながら、これはウイルスにとって問題が多く、それらは大部分新規の配列からなりそして単一の遺伝子を共有してはいない[ref.2、12、13]。別のアプローチでは、特定のクレードを一意に識別する配列[ref.10、11]を使用するが、ゲノムのこれら特定の領域にマッピングされる比較的少数のリードしか使用できず[ref.14]、感度が低下する[ref.15]。これは真核生物のゲノムにとって特に問題であり、それは通常長く、そして大部分が非コード領域から成り、不十分なリード利用をもたらす[ref.5、16、17]。 k-merに基づく最近のアプローチはこれらの問題を克服し、実行時間を劇的に改善した[ref.14、18]。しかしながら、これらのアプローチは、k-merの完全一致を必要とするために感度の低下を示す[ref.15、19]。さらに、それらはしばしば100GBを超える大量のメモリ使用量を要求し、多くのユーザはこれを利用できない[ref.20]。最後に、k-merベースの方法は、サンプルに実際には存在しない多数の低存在量種を予測してしまうことが観察されている[ref.15]。

 本稿では、MiCoP(Microbiome Community Profiling)を紹介する。これは、ウイルスと真核生物を高い精度と感度でプロファイリングできる計算方法である。これは高速マッピングベースのアプローチを利用することによって上記の問題を克服する。すなわち、リード使用率が高く、ウイルスや真核生物に対するバイアスを避け、そして存在量の少ない種に敏感である。マッピングベースのアプローチは、一般的なゴールドスタンダードの方法であるMegablast [ref.21]よりもさらに高い感度を持つことが観察されているが、マッピングベースのアプローチの速度とメモリ使用量に関する懸念があった[ref.22]。我々(著者ら)は、より小さなウイルスおよび真核生物の参照データベースを使用するとき、マッピングベースのアプローチが実行可能かつ好ましいことを実証する。

 まずBWA-MEMを使用してリードをリファレンスゲノムにマッピングし、次に2段階リードアサインメントプロセスを適用する。 BWA-MEMマッピングの結果から、一意にマッピングされたリードが直ちに割り当てられ、続いて、一意にマッピングされたリードの分布に基づいてマルチマッピングされたリードが確率的にゲノムに割り当てられる。各ゲノムにマッピングされたリードの数を計算し、ゲノムの長さで正規化することによって、種の存在量の推定を行う。 MiCoPはリードをリファレンスデータベースの特定のゲノムにマッピングするので、例えば、種からの異なる株またはクロモソームがリファレンスデータベースに別々にリストされている場合、それは種レベルよりも細かい粒度で微生物を検出することができる。本著者らは、その存在量推定パフォーマンスを2つの最も一般的な方法MetaPhlAn2 [ref.11]とKraken [ref.14]と比較することによって検証する。ウイルスおよび真核生物ゲノムからの模擬リードを使い改善された結果を実証し、MiCoPが以前に使用された方法よりもHuman Microbiome Projectデータでより多くのウイルスおよび真核生物を識別できることを示す。


f:id:kazumaxneo:20190714235132p:plain

MiCoP workflow.  論文より転載。

 


インストール

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

本体 GIthub

git clone https://github.com/smangul1/MiCoP.git
cd MiCoP/

#データベースのダウンロードとビルド
# Setup will download about 8.4 GB of data
./setup.sh

python run-bwa.py -h

$ python run-bwa.py -h

python: can't open file 'run-bwa.py-h': [Errno 2] No such file or directory

(ppr_meta) kazu@kazu:/media/kazu/1/MiCoP$ python run-bwa.py -h

usage: run-bwa.py [-h] [--fungi] [--virus] [--output OUTPUT] [--paired]

                  reads [reads ...]

 

Run BWA mem for viruses or fungi.

 

positional arguments:

  reads            Reads to run BWA on. One or two files, depending on single

                   or paired end. Required.

 

optional arguments:

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

  --fungi          Run BWA on fungi. This or --virus required.

  --virus          Run BWA on viruses. This or --fungi required.

  --output OUTPUT  Where to output BWA results to.

  --paired         Treat input reads file as interleaved paired end reads.

> python compute-abundances.py -h

$ python compute-abundances.py -h

usage: compute-abundances.py [-h] [--abundance_cutoff ABUNDANCE_CUTOFF]

                             [--fungi] [--min_map MIN_MAP] [--max_ed MAX_ED]

                             [--normalize {True,False}] [--output OUTPUT]

                             [--pct_id PCT_ID] [--read_cutoff READ_CUTOFF]

                             [--virus]

                             sam

 

Compute abundance estimations for species in a sample.

 

positional arguments:

  sam                   .sam file or file with list of SAM files to process.

                        Required.

 

optional arguments:

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

  --abundance_cutoff ABUNDANCE_CUTOFF

                        Organism abundance to count it as present.

  --fungi               Profile fungi.

  --min_map MIN_MAP     Minimum bases mapped to count a hit.

  --max_ed MAX_ED       Maximum edit distance from a reference to count a hit.

  --normalize {True,False}

                        Normalize species abundance by genome length or not.

                        Default: True

  --output OUTPUT       Output abundances file. Default: abundances.txt

  --pct_id PCT_ID       Minimum percent identity from reference to count a

                        hit.

  --read_cutoff READ_CUTOFF

                        Number of reads to count an organism as present.

  --virus               Profile viruses.

 

実行方法

1、データベースへのマッピング。virusかfiungiかを選択する。

#single-end
python run-bwa.py single.fq.gz [--virus OR --fungi]

#paired-end
python run-bwa.py pair_1.fq.gz pair_2.fq.gz [--virus OR --fungi]

#interleave
python run-bwa.py interleave.fq.gz --paired [--virus OR --fungi]

alignments.samが出力される。 

 

2、abundance profiling。--virusか--fungiを指定する(step1と同じ方を選ぶ)。

python compute-abundances.py alignments.sam [--virus OR --fungi]

出力

f:id:kazumaxneo:20190715035000p:plain


 

引用

MiCoP: microbial community profiling method for detecting viral and fungal organisms in metagenomic samples

Nathan LaPierre, Serghei Mangul, Mohammed Alser, Igor Mandric, Nicholas C. Wu, David Koslicki, Eleazar Eskin
BMC Genomics 2019 20 (Suppl 5) :423