メタゲノム・eDNA解析

先週の実験で、三四郎池の水や様々な食品からDNAを抽出し、PCR増幅しました。本日は、Oxford Nanoporeという次世代シーケンサーで増幅したPCR産物をシーケンスしたデータを使って解析を行います。

PCR増幅したサンプルリスト

サンプルプライマー領域対象生物PCR断片長
三四郎池の水ミトコンドリア12S rRNA魚160-190 bp
発酵食品バクテリア16S rRNAバクテリア1,500 bp程度
加工食品ミトコンドリア16S rRNA魚600 bp程度

本日のデータ解析のワークフローを以下に示します。

大雑把なデータ解析の流れとしては、 まずシーケンスデータの品質チェックをFastQCというツールを用いて行います。それから、コンピュータにインストールしたBLASTを使って、すべてのシーケンスデータをデータベースと照合させ、各リードがどの種と相同性があるかを調べます。それからMEGANを使ってBLASTのヒットを集計し、各サンプルに含まれていたバクテリアの種類を網羅的に解析します。

環境DNAのサンプルは、MEGANの結果で検出された種と、実際に池にいそうな魚なのかの考察や、使用したプライマーの種類(12S全長 or 12S MiFish)の違いが結果に与える影響などを考えてみてください。メタゲノムのほうでは、食品のデータはシャッフルしてお渡ししますので、含まれるバクテリアの種類から、食品サンプルが上記4つの食品のどれであるかを推定します。また、池の水と水槽の水で微生物叢にそれぞれ特徴があるかを考察してください。

サンプルの説明

三四郎池で採水したときの様子

1班：img_5208.jpg

2班：img_8531.jpeg

3班：img_8527.jpeg

4班：img_8528.jpeg

発酵食品

2班：チーズ

3班：ヨーグルト（R1)

加工食品

6: 蒲焼さん太郎

8: ちくわ

PCR結果まとめ

ナノポアシーケンス結果

サンプルリード数

A. 準備（ ツール・データのダウンロード ）

1. ツールのインストール

これからダウンロードするファイルを入れるディレクトリを作成し、その中にすべてのファイルをダウンロードする。とりあえずここではWindowsでダウンロードフォルダの中に「2022jissyu」というフォルダを作ったとする。 

#ダウンロードしたディレクトリに移動する。
cd ~/Downloads/2022jissyu/

#ダウンロードしたファイルを解凍する。
tar vxf seqkit_windows_amd64.exe.tar.gz
#ダウンロードしたファイルを解凍する。
tar vxf ncbi-blast-2.13.0+-x64-win64.tar.gz

2. データベースのダウンロード

B. クオリティチェック

FastQCは主にIllumina用のクオリティチェックツールなので、Nanoporeのデータに対しては適切な評価ができておらず、評価値の〇×は気にしなくてよいです。下記はシーケンスデータのクオリティスコアに関する、平均値等の情報。

Nanoporeのデータだとクオリティスコア10強(精度90%強)となるはずである。

クオリティスコアQは、エラーの生じる確率 perror から下記のように計算されます。 (出典： https://bi.biopapyrus.jp/rnaseq/qc/fastq-quality-score.html )

リード長の分布が想定される長さになっているかなども確認すること。 (16S全長は1.5 kbp程度、12S全長は1 kbp程度、12S MiFishは250 bp程度)

FASTQCを実行するには…

cd ~/Downloads/FastQC/
chmod 755 fastqc
./fastqc

FastQCでシーケンスファイルを開く

Ctrlキーを押しながらクリックすれば、複数のファイルを一度に選択できるので、必要なファイルを選択して開く。

結果のタブを開いてみると、次のようなQC結果が見られる。

C. BLAST

1. FASTQ→FASTA変換

FASTQファイルをメモ帳などで開いてみると次のように表示されます。

FASTQ形式は 塩基配列とクオリティが両方含まれるファイル形式になりますが、BLASTのプログラムはFASTQ形式では入力ファイルとして受け付けてくれません。BLASTが対応しているのは、塩基配列だけの情報であるFASTA形式になります。

FASTA形式は1つのシーケンスにつき、「>」 で始まる1行のヘッダ行と、2行目以降のACGTのシーケンス文字列で構成されます。

そこで、まずFASTQのクオリティを削除して、FASTA形式に変換します。 

./seqkit fq2fa ./2022nanopore/input_file.fastq -o output_file.fasta

input_file.fastq, output_file.fastaは適当に変更すること。

コマンド入力時の便利キー

カーソルの上・下キー前に入力したコマンドを呼び出す。
Ctrl＋Shift＋V貼り付け
tabキーファイル名・コマンドの自動補完
Ctrl+Cコマンド強制終了(blastを実行中に止めたい場合など)

2. BLASTデータベース作成

MitoFishデータベースのFASTAファイルから、makeblastdbコマンドによってBLASTのデータベースを作成します。ターミナルで下記のように入力します。 

# MitoFish (魚類用)
# MitoFishにヒトミトコンドリア配列を追加したファイルを作成する。
## Windowsの場合
cmd /C "type .\complete_partial_mitogenomes\mito-all sequence.fasta > mitofish-human.fasta"

## Mac/Linuxの場合
cat mito-all sequence.fasta > mitofish-human.fasta

#Blastデータベースを作成
./ncbi-blast-2.13.0+/bin/makeblastdb -in mitofish-human.fasta -dbtype nucl

-in には入力となるFASTAファイルを指定します。

-dbtype はFASTAファイルがDNA配列であれば「nucl」、アミノ酸配列であれば「prot」を指定します。

3. BLAST検索

FASTA形式に変換したシーケンスデータをクエリーとして、バクテリア16SデータベースもしくはMitoFishデータベースで相同性検索を行います。ターミナルで次のように入力してください。 

#バクテリア用
./ncbi-blast-2.13.0+/bin/blastn -db 16S_ribosomal_RNA -query input.fasta -num_threads 16 -out input.fasta.blastn

#魚類用
./ncbi-blast-2.13.0+/bin/blastn -db mitofish-human.fasta -query input.fasta -num_threads 16 -out input.fasta.blastn

-db には、BLASTデータベースファイルのprefixを指定します。今回の場合は元となったSILVAもしくはMitoFishのFASTAファイルを指定すれば良いです。

-query にはデータベースに対して相同性検索を行いたい問い合わせ配列が含まれるFASTAファイルを指定します。ここではFASTA形式に変換したシーケンスデータを指定します。(ファイル名は適宜変更すること。)

-num_threads には並列計算時に使用するCPUの数を指定します。演習で使用している研究室のLinuxはCPUが16スレッドあるので、ここでは例として16を指定しています。自分のPCで計算するときは、各自のPCに合わせて設定すること。

-out には出力ファイルの名前を書く。(拡張子はMEGANに読み込ませるために「.blastn」とし、ファイル名は出来れば「入力FASTAファイル名」＋「.blastn」としておくとMEGANがFASTAファイルを認識してくれて情報がリッチになる。)

4. BLAST結果例

下記のように入力すると、結果ファイルの中身を先頭30行だけ見ることが出来ます。 

## Windows (PowerShell)の場合
Get-Content -head 30 input.fasta.blastn

## Mac/Linuxの場合
head -n 30 input.fasta.blastn

# input.fasta.blastn は適当なファイル名に変更します。

D. MEGANによるBLAST結果集計

1．MEGANではLowest Common Ancestor (LCA)法によってtaxon (分類群)を推測している。LCA法とは、1つのリードがデータベース中の複数の配列にヒットした場合、ヒットした配列に共通の分類群を求める方法です。　http://bio4j.com/blog/2012/02/finding-the-lowest-common-ancestor-of-a-set-of-ncbi-taxonomy-nodes-with-bio4j/ ただ、MEGANのデフォルト値はIlluminaシーケンサー用(精度99.9%)に設定されているので、ナノポアのように80～90%程度の精度では変更する必要があります。

2．インストールしたMEGANを起動します。

3．MEGANが起動したら、「File」→「Import From BLAST」からBLAST結果を開きます。

4．ダイアログ右上のボタンをクリックして、BLASTの結果ファイルを選択し、「Open」をクリックします。(複数ファイルを選択できますが、結果がマージされてしまうので、ファイルを一つだけ選択します。}

5．「LCA Params」タブを開いて、Top Percent: の値を0.5に変更しておきます。このパラメータは、BLASTの結果の中で最もスコアの高いトップヒットからどの程度離れたヒットまで使用するかの閾値になります。ナノポアではシーケンス精度が悪く、無関係な生物も似たようなスコアでヒットしてしまうため、ほぼトップヒットしか使わないように厳しめに閾値を設定しておきます。それから、Min Score: の値をバクテリア16Sではリード長が1500bp程度なので1000、魚類16Sではリード長が600bp程度なので300、魚類12Sではリード長が200bp程度なので100などと指定し、スコアの低いリードをトリミングします。

6．結果が表示されたら、 種名まで表示するために、「Rank」→「Species」を選択します。

　　　　　　　　　　　　　　　　　　↓

7．マウスホイールで拡大・縮小できます。もしくはメニューバーの下記のボタンをクリックしても拡大できます。

8．各生物種に割り当てられたBLASTの生データを見たい場合、ノードを右クリックし、「Inspector」を開けば見れます。本当にその生物はいるのか確認したくなった時などに、アライメント長や塩基配列そのものを確認するのに便利です。

9．BLASTファイルを読み込むときに設定したLCAのパラメータを後から変更したい場合は、「Options」→「Change LCA Parameters」から可能です。

10．MEGANで複数サンプルを比較する場合は、blastnファイルを複数回開いたら、開いたウィンドウをそのままにした状態で、「File」→「Compare」とクリックして比較したいサンプルを選択し、「Apply」をクリックします。

11．複数結果の表示を変更する場合、円グラフのボタンなどを押すと良いです。

12．ExcelやR等で解析する場合は、データをタブ区切りテキスト形式でExportする必要があります。Exportしたいノードをクリックして(全部Exportする場合は全部選択(WidnowsならCtrl+a)して)、「File」→「Export」→「CSV Format」をクリックします。

13．Exportするデータを「taxonPathtocount」に変更します。(ほかのデータ形式でも勿論可)

14．適当な名前を付けて保存します。

E. データの転送、Excelでの解析

リモートデスクトップ先で「ファイル管理」を開いて、保存したファイルをコピーしてから、手元のWindowsでファイルエクスプローラーを開いて貼り付けるとファイルを簡単に転送できる。もしくは、メールやGoogle Driveなどで転送するなどしてもよい。

ExportしたファイルをExcelで開くには、Excelを起動しておき、ExportしたファイルをExcel上にドラッグアンドドロップすれば良いです。

Excelでデータの概要を把握するのに役立つテクニックとして、条件付き書式を設定することで、データの大小を一目でわかるようにできたりします。

そのほか、「データ」→「フィルター」を使ってみたり、グラフを描いてみたりするのが通常の解析の流れになるかと思います。

F. 明日の内容

各班次の内容について「目的」、「方法」、「結果」、「考察」の4つのパートを明確に区別してプレゼンテーションを作成する。班ごとに発表し、発表時間は質疑応答を入れて30分。 

1班．食品の品種判別 by サンガー
2班．三四郎池のeDNA
3班．発酵食品のメタゲノム、加工食品の品種判別 by ナノポア
4班．三四郎池のメタゲノム

各テーマごとに例えば下記のような項目について考察をすること。インターネットを積極的に使用して調べることを推奨します。また、ある程度調べてもわからないことがあればTA・スタッフに聞いてみてください。

・1班．食品の品種判別 by サンガー

NCBIのデータベースとMitoFishのデータベースを比べて、ヒットした種が同じかどうか調べ、どちらのデータベースのほうが良さそうか考えてみる。

ヒットした近縁種の配列をGenbankからダウンロードして加え、系統樹を描いてみる。

手法で詳しく説明して欲しい箇所：「DNA抽出」（使用したキットはDNeasy Blood & Tissue Kitsです。）

・2班．三四郎池のeDNA

検出された魚は三四郎池に棲息していそうな魚かどうか。

去年の三四郎池のデータとも比較してみる。http://suikou.fs.a.u-tokyo.ac.jp/yosh_data/2021jissyu/2021nanopore.zip (group2-12S-Sanshiroike1.fq) 今年はブラックバスが密かに話題になっているらしいが、去年と比べてブラックバスが増えていそうか。

手法で詳しく説明して欲しい箇所：「電気泳動、DNA精製」（使用したキットはFastGene™ Gel/PCRExtractionキットです。）

・3班．発酵食品のメタゲノム、加工食品の品種判別 by ナノポア

今回発酵食品で検出されるバクテリアはほぼ1種類だと思うので、精度の悪いナノポアのリードの精度を向上させる方法を実践してみてください。具体的にはGeneiousでマルチプルアライメントを作成して、コンセンサス配列を作ることで、NCBIのBlastで一致率99%程度のヒットが得られるようになることを確認し、ナノポアのリードはどのような間違いが多いのか考察してみてください。

加工食品から検出された魚は妥当でしょうか。

手法で詳しく説明して欲しい箇所：「PCR」（使用したDNAポリメラーゼはrepliQa HiFi ToughMixです。AmpliTaq GoldやEx Taqといった他の酵素と比較して、どういった特徴があるでしょうか。）

・4班．三四郎池のメタゲノム

検出されたバクテリアは淡水環境で良く検出されているでしょうか？

去年の三四郎池のデータとも比較してみる。http://suikou.fs.a.u-tokyo.ac.jp/yosh_data/2021jissyu/2021nanopore.zip (group2-16S-Sanshiroike1.fq)

手法で詳しく説明して欲しい箇所：「ナノポアシーケンシング」（使用したライブラリー調整キットはSQK-LSK110です。公式マニュアル：http://suikou.fs.a.u-tokyo.ac.jp/yosh_data/2022jissyu/amplicons-by-ligation-sqk-lsk110-ACDE_9110_v110_revT_10Nov2020-minion.pdf )

G. 課題

明日のプレゼン資料の完成版をファイルに保存して、11月15日（火）23:55までにITC-LMSの課題に提出すること。