先週の実験で、三四郎池、葉山の海水、葛西臨海水族園の水槽の水からDNAを抽出し、12S rRNAの全長(1,000 bp程度)、もしくは12S rRNAのMiFishプライマー領域(160-190 bp)をPCR増幅しました。また、環境水や発酵食品の「ヨーグルト」、「納豆」、「イカの塩辛」からから抽出したDNAを、16S rRNA全長(1,500 bp程度)をPCR増幅しました。本日は、Oxford Nanoporeという次世代シーケンサーで増幅したPCR産物をシーケンスしたデータを使って解析を行います。

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

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

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

サンプルの説明

2020年食品サンプル

2020年環境水サンプル

PCR結果まとめ

2020年PCR結果まとめ

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

ツールのインストール

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

SILVA database ・・・16S rRNA (バクテリア・植物の葉緑体、合計約60万件(真核生物のミトコンドリア上の16Sは含まれない))、18S rRNA (真核生物核ゲノム、約8万件)を網羅的に集めて整理したデータベース。ただし哺乳類の網羅率は高くはなく、哺乳類は407件登録されいてるが、そのうちヒトが162件。(2019年調査)
https://www.arb-silva.de/fileadmin/silva_databases/current/Exports/SILVA_138.1_SSURef_NR99_tax_silva_trunc.fasta.gz
MitoFish ・・・魚類のミトコンドリアのデータベース。約2万種が登録されている。魚類ミトコンドリア上の16S rRNA、12S rRNAを探す場合に使用する。
http://mitofish.aori.u-tokyo.ac.jp/files/complete_partial_mitogenomes.zip
Nanoporeのシーケンスデータ・・・メタゲノムとeDNAのデータが入っています。
http://suikou.fs.a.u-tokyo.ac.jp/yosh_data/2020jissyu/nanopore.zip
- eDNA解析: 12S-xxx.fqのファイルを使用
- メタゲノム解析（食品）: 16S-food-xxx.fqのファイルを使用
- メタゲノム解析（水）: 16S-eDNA-xxx.fqのファイルを使用

B. クオリティチェック

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

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

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

クオリティスコアが 10 ならば、シーケンシングエラーが生じる確率は 10.0% であるから、読み取られた塩基の信頼度は 90.0% 。
クオリティスコアが 20 ならば、シーケンシングエラーが生じる確率は 1.0% であるから、読み取られた塩基の信頼度は 99.0% 。
クオリティスコアが 30 ならば、シーケンシングエラーが生じる確率は 0.1% であるから、読み取られた塩基の信頼度は 99.9% 。

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

FASTQCを実行するには、シーケンスデータをダウンロードしたディレクトリが「Downloads」フォルダの場合、次のようにする。

cd Downloads           # Downloadsフォルダへ移動する
unzip nanopore.zip     # nanopore.zipファイルを解凍する
fastqc nanopore/*.fq   # nanoporeフォルダ以下のすべてのFASTQファイルに対してFASTQCを実行する

出来上がったhtmlファイルを開いてみると、次のようなQC結果が見られる。

C. BLAST

FASTQ→FASTA変換

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

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

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

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

awk 'NR%4==1{print ">"substr($1,2)} NR%4==2{print $0}' input_file.fastq > output_file.fasta

> output_file.fasta をつけると、出力をファイルに保存できます。 input_file.fastq, output_file.fastaは適当に変更すること。

BLASTデータベース作成

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

## SILVA
#gzipファイルを解凍する
gzip -d SILVA_138.1_SSURef_NR99_tax_silva_trunc.fasta.gz

#Blastデータベースを作成
makeblastdb -in SILVA_138.1_SSURef_NR99_tax_silva_trunc.fasta -dbtype nucl


## MitoFish
#ファイルを解凍しておく
unzip complete_partial_mitogenomes.zip   #complete_partial_mitogenomes.zipを解凍する

#都合により、ヒトのミトコンドリアも追加しておく
#ヒトのFASTAファイルダウンロードする。curlはファイルをダウンロードするためのコマンド。
curl -o NC_012920.fasta "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=NC_012920&rettype=fasta&retmode=text"

#MitoFishにヒトを追加する。catはファイルの中身を表示するためのコマンド。
cat NC_012920.fasta complete_partial_mitogenomes.fa > MitoFish-Human.fasta

#Blastデータベースを作成
makeblastdb -in MitoFish-Human.fasta -dbtype nucl

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

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

BLAST検索

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

blastn -db SILVA_138.1_SSURef_NR99_tax_silva_trunc.fasta -query input.fasta -num_threads 8 -out input.fasta.blastn

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

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

-num_threads には並列計算時に使用するCPUの数を指定します。演習で使用しているノートPCはCPUが8スレッドあるので、8を指定します。

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

BLAST結果例

D. MEGANによるBLAST結果集計

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%程度の精度では変更する必要があります。

まずはMEGANを起動します。ターミナルで下記を入力してMEGANを起動してください。

~/megan/MEGAN

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

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

「LCA Params」タブを開いて、Top Percent: の値を0.5に変更しておきます。このパラメータは、BLASTの結果の中で最もスコアの高いトップヒットからどの程度離れたヒットまで使用するかの閾値になります。ナノポアではシーケンス精度が悪く、無関係な生物も似たようなスコアでヒットしてしまうため、ほぼトップヒットしか使わないように厳しめに閾値を設定しておきます。

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

　　　　　　　　　　　　　　　　　　↓

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

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

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

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

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

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

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

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

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

メールなどで転送しても良いし、ファイル置き場を研究室のサーバのほうに作ったので、Linuxのブラウザ(Firefox)で下記のWEBサイトにアクセスして

http://133.11.222.89:7080/

ID: training
Password: programming

でログインし、自分の名前のフォルダを作ってその中に保存し、自分のPCのブラウザで開いてダウンロードしてもよい。

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

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

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

F. 課題

食品1～6がどの食品だったか推理し、Meganから出力されたCSVファイルを読み込んだExcelファイルに記入する。11月5日（木）23:55までに、食品メタゲノムについて作成したExcelファイルをITC-LMSの課題に提出すること。最後まで完走したかどうかの確認なので、本日行った結果をとりあえず提出してもらえれば良いです。

G. 明日の実習

各班次の内容についてプレゼンテーションを作成する。班ごとに発表し、発表時間は質疑応答を入れて20分。

1．制限酵素
2．池・海メタゲノム
3．プラスミド導入
4．魚肉
5．eDNA
6．食品メタゲノム

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

・プラスミド抽出～制限酵素地図作成

EcoRIは1か所切断、BamHIは2か所切断のはずですが、それ以上にバンドが見えた版画多かったと思います。理由を考察してください。

・魚肉からの魚種判定

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

・食品からのメタゲノム解析

「イカの塩辛」、「ヨーグルト」、「カスピ海ヨーグルト」、「チーズ」をメタゲノムの結果からどの食品だったかを推測してください。理由も説明してください。

・水からのメタゲノム解析

三四郎池（7月、10月）と水族館の水槽でそれぞれ特異的なバクテリアはいるでしょうか？

・eDNA解析

水槽、三四郎池、相模湾でeDNAから見つかる魚にはどんな種類があるでしょう？12S 全長と12S MiFish領域で差はみられるでしょうか。どちらが実際の生物の種類・比率に近いでしょうか。

・プラスミド導入

Zoomによる説明と見学のみで、実際に手を動かして実験したわけではないのでイメージするのは難しいと思いますが、今回大腸菌へのプラスミド導入や魚卵へのプラスミドのインジェクションについて実験のプロトコールをWEB等で調べて説明してください。考察するのは難しいと思うので、考察の代わりに水産対象種でのゲノム編集やトランスジェニックの現状と将来展望について調べ、まとめてください。

メタゲノム・eDNA解析