macでインフォマティクス

macでインフォマティクス

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

rDNA配列を解析するためのユーザーフレンドリーなバイオインフォマティクスパイプライン gDAT

 

 複数の生物を並行して行うハイスループットシーケンス(HTS)(メタバーコーディング)は、環境サンプル中の微生物群集を分析するための日常的で費用対効果の高い方法となっている。しかし、HTSデータの潜在的なエラーを特定するためには、慎重なデータ処理が必要であり、また、HTSによって生成された大量のデータは、下流の解析のためのコマンドラインツールを使いこなす必要がある。本論文では、最も一般的なコマンドラインツールを、使いやすいグラフィカルなインターフェースであるgDATに組み込んだパイプラインを紹介する。パイプラインは、Pythonスクリプト言語を使用することで、最新のWindowsmacOSLinuxオペレーティングシステムと互換性がある。パイプラインは、Sanger、454、IonTorrent、Illumina、PacBioの各配列の解析をサポートし、品質フィルタリングステップのカスタム変更が可能で、配列同定のためにオープンリファレンスとクローズドリファレンスの両方の運用分類学的ユニットピッキングを実装している。事前に設定したパラメータは、アーバスキュラー菌根菌由来のスモールサブユニット(SSU)rRNA遺伝子アンプリコンの解析に最適化されているが、このパイプラインは、幅広い生物を対象としたメタバーコーディング研究に広く適用できる。このパイプラインは、SSU遺伝子領域の一般的な真核生物用プライマーと、内部転写スペーサー(ITS)マーカー領域の真菌用プライマーを用いたデータでもテストされた。パイプラインの設計について説明し、Illuminaプラットフォームでシーケンスされた異なるマーカー領域を用いたサンプルデータセットの解析を行うことで、その性能と速度を評価した。グラフィカルなインターフェースは、必要に応じてコマンドラインを使用することもでき、再現性とログ機能を備えた迅速なデータ解析のためのアクセス可能なツールとなっている。ソフトウェアをオープンソースにすることで、コードへのアクセス性を最大限に高め、コミュニティによる精査やバグ修正を可能にしている。

 

マニュアル・チュートリアルは配布されているtar ballに含まれている。

f:id:kazumaxneo:20210615151756p:plain

 

マニュアルより

グラフィカルなインターフェースにより、ターミナルやコマンドライン・プロンプトを使わずにコマンドを作成することができます。これにより、パラメータの入力時のミスや、名前の入力ミスを防ぐことができます。MaarjAMデータベースを用いてSSU真菌データを解析するための一般的なワークフローを提供していますが、幅広い生物やアンプリコンに対応しています。解析は、de novo、closed-またはopen-OTU(operation taxonomic unit)ピッキングを使用して行うことができます。本パイプラインは、Sanger、454、IonTorrent、Illuminaのシングルエンドおよびペアエンドリード、PacBioのリードをサポートしています。パイプラインは、FASTA/QUALおよびFASTQファイルを含むパックされた*.gzファイルの読み込みもサポートしています。ファイルが*.zipまたは*.rarファイルでパックされている場合は、リードをクリーニングする前にファイルを解凍する必要があります。MaarjAMデータベースを用いたSSU真菌解析には、ほとんどの定義済みパラメータが最適です。

 

インストール 

リリースからmacos, windows, linuxバージョンをダウンロードできる。

https://github.com/ut-planteco/gDAT/releases

f:id:kazumaxneo:20210516210808p:plain

 

 

ダウンロードして解凍する。

f:id:kazumaxneo:20210516211024p:plain

Macos

f:id:kazumaxneo:20210615151202p:plain

 

 

 

実行方法

ターミナルで gdat.pyをランする。

python gdat.py

GUI版が立ち上がる。

f:id:kazumaxneo:20210615151625p:plain

 

 

テストラン

テストデータには、illuminaのWANDA-AML2 primersでシーケンスされたデータ;raw fastq/ (multiplexed)、それをdemultiplexingしたillumina_demultiplexed/、そして454のNS31フォワードプライマーでシーケンスされたデータ;454 (multiplexed)が含まれる。

f:id:kazumaxneo:20210616105049p:plain

(データは小さくなるよう処理されている)

 

1、シークエンシングリードのバーコードやプライマーの出現率や配列の品質分布の分析

Analyse FASTQ filesボタンを選択し、出てきたウィンドウのSelect FASTQ files separatelyボタンを選択し、個別にfastqを選択する。もしくはSelect folderでフォルダごとまとめて選択する。

f:id:kazumaxneo:20210616113016p:plain
データセットに応じて、バーコードの長さ、プライマーの長さ、またはそれらの組み合わせであるkmerの長さを定義する(パラメータについてはマニュアルを参照)。Runボタンで実行する。

 

出力

all.kmer20.stats.txt (R1だけ分析した)

最も豊富なkmerがトップに表示され、下の方に配列の品質分布が表示される。

f:id:kazumaxneo:20210616110115p:plain

もっとも豊富な20-merはPCRフォワードプライマーであるCAGCCGCGGTAATTCCAGCT配列になった。

 

2、リードのクリーニング

1. Clean and quality filter readsを選択する。

f:id:kazumaxneo:20210616111811p:plain

 

ウィンドウが出てくる。demultiplexするにはsingle/Paired-end fastq~を選択して、fastqを個別に選択する。demultiplex済みのfastqの場合(通常ユーザーが受け取るのはこちら)、Demultiplexed FASTA/FASTQ readsを選んで、フォルダを選択する。

f:id:kazumaxneo:20210616112607p:plain


single/Paired-end fastq~を選択した時

配列がdemultiplexされていない場合は、R1のフォワードリードとR2/R4(オリゴファイルが提供されているかどうかによる)のリバースリードのFASTQファイルを定義する必要がある。オリゴファイルが提供されていない場合(R2およびR3)、これらのファイル(R1およびR2)には、配列と一緒にバーコードがインラインで含まれている必要がある。FASTQのヘッダーに書かれているバーコードはサポートされていない。オリゴファイルR2およびR3が提供されている場合、これらのファイルにはバーコード配列が含まれており、R1のフォワードリードとR4のリバースリードを異なるサンプルに分離するために使用される。

f:id:kazumaxneo:20210616112843p:plain

 

配列がdemultiplexされていない場合、タブ区切りのサンプルシートファイルも定義する必要がある。

laelatu_samples.txt

f:id:kazumaxneo:20210616113939p:plain

このファイルには、サンプル名と使用したバーコード配列(使用するファイルに応じてフォワードまたはリバース)が含まれる。サンプル名、使用したバーコード配列(使用するファイルに応じてフォワードまたはリバース)、プライマー配列を含めることができるが、これらはprimerパラメータで上書きも可能。スプレッドシートファイルに情報がある場合は、データをテキストファイルにコピー&ペーストすることができる(各列にタブが自動的に生成される)。

f:id:kazumaxneo:20210616113856p:plain

パラメータはSSUアンプリコンを使用したIlluminaシーケンシング向けに最適化されている。パラメータは変更できるが、パラメータフィールドの最後にあるチェックボックスを有効にするのを忘れないようにする。

 

左下のforwardとreverseのprimerの選択をしてRunする。

f:id:kazumaxneo:20210616125844p:plain

demultiplexされて1つのfastq (R1とR2両方含むインターリーブfastq) として出力される。

 

demultiplexed FASTA/FASTQを選択した時

demultiplexされたfastqのフォルダを指定してRunする。

f:id:kazumaxneo:20210616130012p:plain

Run前にユーザーがプライマーを定義する必要がある(右下)。

=> demultiplexされて1つのfastq (R1とR2両方含むインターリーブfastq) として出力される。

 

3、ペアエンドリードのマージ (Combine paired-end reads)

FLAShソフトウェアを使用してフォワードリードとリバースリードをマージする。クリーニングステップで、インターリーブされたFASTQファイルを作成されており(奇数リードが順方向リード、偶数リードが逆方向リード)、これを指定する。配列は最後にマージされ、逆方向のリードは自動的に逆相補鎖に変換される。良好な結合率を得るためには、配列が十分にオーバーラップしている必要があり、低品質の領域を取り除くために品質に基づいて配列の末端をトリミングする必要がある。

f:id:kazumaxneo:20210616130830p:plain

demultiplexed.cleaned.combined.fastqとdemultiplexed.cleaned.combined.fastaが出力される。

 

4、(任意)配列のクラスタリング(Cluster reads)

計算量をさらに減らす方法の1つとして、配列同定の前に配列をクラスター化することがある。安全のためには、BLAST+による同定の際に使用した同一性(例:97%)よりも高い同一性(例:98%または99%でクラスタリング)を使用するのがベストである。クラスター化に成功するためには、配列が同じストランド上にあり、同じアンプリコン位置から始まる必要がある。キメラチェックはクラスタリングの前に行うこともできるが、処理を高速化する目的で、パイプラインはクラスタリング後のキメラチェックもサポートしている。入力には、FASTA ファイルを使用する(FASTQ は FASTA に変換する)。出力ファイル名は、入力名とクラスタリングIDに基づいて自動的に生成され、以下の拡張子を持つ。*.cluster98.uc (vsearchのクラスタセントロイドとクラスタエレメントを含む), *.cluster98.fasta (FASTA配列形式のクラスタセントロイド配列を含む), *.cluster97.cluster (従来のBLASTclustと同様に、各クラスタとそのエレメントを別の行に含む)。このクラスター・セントロイド・ファイルは、次のステップで結果を個々の配列にマッピングする際に使用する。ピボットテーブルを作成する際、クラスターファイル(*.uc)を定義すると、ピボットテーブルのセルのカウントにクラスター内のすべての配列を含めることができる。

f:id:kazumaxneo:20210616131327p:plain



5、キメラリードの除去(Remove chimeric reads)

PCRで増幅された配列データには、2つ以上の異なる生物から構築された人工的な配列であるキメラ配列が含まれることがある。キメラ検査には、データベース方式とde novo方式がある。de novo法では、入力配列を97%でクラスタリングし、クラスタリングされた配列を参照データベースとして用いて、任意の入力配列が単一生物のものであるか、異なる生物のものを組み合わせたものであるかの確率を算出する。De novo法、参照データベース法のいずれの場合も、配列をより小さな断片に分割し、各断片をデータベースから最も近い候補、またはクラスター化されたセントロイドの中の候補と照合することで行う。参照データベース方式では、データベース自体がキメラフリーであることを前提としている。

f:id:kazumaxneo:20210616131554p:plain

 

6、BLAST+を使ったTaxonomyアサイン(Identify reads with BLAST+)

FASTAから適切なリファレンスデータベース(MaarjAM、GB/NCBIなど)に対してBLAST+検索を行い、分類学上の割り当てを行う。

f:id:kazumaxneo:20210616131911p:plain



7、(任意)vsearchを使ったTaxonomyアサイン(Identify reads with vsearch)

FASTA形式の適切な参照データベース(MaarjAMなど)に対するvsesarch検索を用いて分類学上の割り当てを行う。出力はBLAST+の出力と一致するように変換されるため、異なる結果をマージすることが可能で、追加のステップはvsearchの結果と互換性がある。

f:id:kazumaxneo:20210616132133p:plain

 

8、BLAST+からピボットテーブルを生成(Generate pivot table from BLAST+)

サンプル対taxonのピボットテーブルを生成する。良好なヒットを得るためには、クエリ配列が、十分に完全なアラインメントと高い同一性を持つデータベースの対象配列と一致する必要がある。例えば、454配列の品質がリードの最後で低下する場合や、リファレンスセットの大部分と異なるプライマーを使用する場合など、より多くのヒットプールを取得するために、低いアライメント閾値を使用することが望ましい場合もあるし、trimmingパラメータを使って配列を短くカットし、再度クリーニングすることが望ましい場合もある。また、アンプリコン内の特定の位置にヒットした場合にフィルタリングするオプションもある。

パイプラインに含まれるMaarjAMデータベースの配列は、NS31プライマーとAML2プライマーの間でトリミングされている。この領域の中で最も変化しやすい部分は、NS31プライマーから70~300bpの間に位置している。この情報を利用してアライメント可変領域パラメータを指定すると、入力配列がアンプリコン内にランダムに配置されているイルミナのタグメンテーションを利用する際に便利である。このスクリプトを実行すると、2つのピボットテーブル(うち1つは転置(transpose)されている)、BLASTヒットがない配列を含むFASTAファイル、の合計3つのファイルができる。クラスタリングを使用した場合は、クラスタリングされたセントロイドを含む追加のnohit FASTAファイルが生成される。nohitsファイルを生成するには、ユーザーが入力FASTAファイルを定義する必要があり、そうでない場合はピボットテーブルのみが生成される。

f:id:kazumaxneo:20210616132914p:plain

 

9、クラスタから配列を選ぶ(Pick sequences from clusters)

nohit配列をさらに系統分析に利用したり、推定上の新種の分類を同定するためには、配列の多様性の中から代表的なサンプルを選択することが有用な場合がある。例えば、Glomeromycotina nohitsの中から新しい分類を探す場合、このスクリプトでは、nohit配列をクラスター化(例:97%)し、各クラスターに対していくつかの配列を返すことができる。nohitsをより深く分析することで、反復的な配列同定プロセスが可能になる。

 

全ての解析手順はマニュアルに書かれています。マニュアルを読みながら進めてください。

引用

User-friendly bioinformatics pipeline gDAT (graphical downstream analysis tool) for analysing rDNA sequences

Martti Vasar, John Davison, Lena Neuenkamp, Siim-Kaarel Sepp, J Peter W Young, Mari Moora, Maarja Öpik

Mol Ecol Resour. 2021 May;21(4):1380-1392