研究メモブログの目次

トップ > 【1】16S rRNA 菌叢解析やってみた（データのインポート～DADA2 によるデノイジング by Qiime2）
最終更新日

【1】16S rRNA 菌叢解析やってみた（データのインポート～DADA2 によるデノイジング by Qiime2）

こんにちは～。ガイです。生きてます。博士課程の荒波に揉まれて瀕死になっていました。現在D3なのですが、なんとこの春から研究テーマが変わりました。その辺りからも崖っぷち度合いはお察しいただけると思います。大人の事情でふたたびバイオインフォっぽいテーマに戻ってきました。腸内細菌叢解析です。お金がないので、とりあえず16S 解析やってみましょうということで、デビューしました。今後も何回か同じ解析をやるので備忘録を書いていきます。

腸内細菌叢は超新人なのでお手柔らかにお願いいたします。間違いや思うところがあったら何卒コメントをください。

解析を進めるにあたって、北大生命科学院 修士卒の太田 圭祐さんのQiita 記事をたくさん参考にいたしました。本当にありがとうございます。

qiita.com

Qiime2 のインストール

インストールにはconda を使いました。

参考：Natively installing QIIME 2 — QIIME 2 2023.2.0 documentation

wget https://data.qiime2.org/distro/core/qiime2-2023.2-py38-linux-conda.yml
conda env create -n qiime2-2023.2 --file qiime2-2023.2-py38-linux-conda.yml
rm qiime2-2023.2-py38-linux-conda.yml
conda activate qiime2-2023.2

Qiime2 環境に入ります。

conda activate qiime2-2023.2
# スパコンの場合はsource activate qiime2-2023.2 じゃないと動かないかも？

インストールが上手くいっているかを確認

qiime --help

上手くいっていれば

Usage: qiime [OPTIONS] COMMAND [ARGS]...

  QIIME 2 command-line interface (q2cli)
  --------------------------------------

  To get help with QIIME 2, visit https://qiime2.org.

  To enable tab completion in Bash, run the following command or add it to
  your .bashrc/.bash_profile:

      source tab-qiime

  To enable tab completion in ZSH, run the following commands or add them to
  your .zshrc:

      autoload -Uz compinit && compinit
      autoload bashcompinit && bashcompinit
      source tab-qiime

Options:
  --version   Show the version and exit.
  --help      Show this message and exit.

Commands:
  info                Display information about current deployment.
  tools               Tools for working with QIIME 2 files.
  dev                 Utilities for developers and advanced users.
  alignment           Plugin for generating and manipulating alignments.
  composition         Plugin for compositional data analysis.
  cutadapt            Plugin for removing adapter sequences, primers, and
                      other unwanted sequence from sequence data.
  dada2               Plugin for sequence quality control with DADA2.
  deblur              Plugin for sequence quality control with Deblur.
  demux               Plugin for demultiplexing & viewing sequence quality.
  diversity           Plugin for exploring community diversity.
  diversity-lib       Plugin for computing community diversity.
  emperor             Plugin for ordination plotting with Emperor.
  feature-classifier  Plugin for taxonomic classification.
  feature-table       Plugin for working with sample by feature tables.
  fragment-insertion  Plugin for extending phylogenies.
  gneiss              Plugin for building compositional models.
  longitudinal        Plugin for paired sample and time series analyses.
  metadata            Plugin for working with Metadata.
  phylogeny           Plugin for generating and manipulating phylogenies.
  picrust2            Predicts gene families and pathways from 16S sequences.
  quality-control     Plugin for quality control of feature and sequence data.
  quality-filter      Plugin for PHRED-based filtering and trimming.
  rescript            Pipeline for reference sequence annotation and curation.
  sample-classifier   Plugin for machine learning prediction of sample
                      metadata.
  taxa                Plugin for working with feature taxonomy annotations.
  vsearch             Plugin for clustering and dereplicating with vsearch.

と表示されます。

 

生データのインポート

つぎに、QIIME 2 artifact に生配列データをインポートします。

参考：Importing data — QIIME 2 2022.2.0 documentation

Manifest file の作成

複数サンプルある場合には、サンプルの識別子と、そのサンプルのシーケンスおよびクオリティデータを含む fastq.gz または fastq のabsolute filepath を対応付けるファイル（ManifestFile.csv）を用意すると一気にインポートができます。

参考：https://docs.qiime2.org/2022.2/tutorials/importing/#fastq-manifest-formats

Excel で自作してもいいのですが、太田さんのブログを参考にコマンドで作成することもできます。生データ配列の入ったフォルダの二個上のフォルダで以下のコマンドを実行するとManifestFile.csv ファイルが出力されます。

ls */* | \
     awk -F'.' -v P="$(pwd)" 'BEGIN { print "sample-id,absolute-filepath,direction" } \
     {if(NR%2==0){print $3","P"/"$0",reverse"}else{print $3","P"/"$0",forward"}}' > ManifestFile.csv

自作する場合にはNanaimo さんの以下のブログが参考になります。

note.com

 

ManifestFile.csv の中身はこんな感じになります。パスは絶対パスです。サンプルID が一番左に来るようにします。

sample-id,absolute-filepath,direction
HAPPY-1,your/file/path/rawdata/HAPPY-1.1.fastq,forward
HAPPY-1,your/file/path/rawdata/HAPPY-1.2.fastq,reverse
HAPPY-2,your/file/path/rawdata/HAPPY-2.1.fastq,forward
HAPPY-2,your/file/path/rawdata/HAPPY-2.2.fastq,reverse
HAPPY-3,your/file/path/rawdata/HAPPY-3.1.fastq,forward
HAPPY-3,your/file/path/rawdata/HAPPY-3.2.fastq,reverse
(中略)
Wild-1,your/file/path/rawdata/Wild-1.1.fastq,forward
Wild-1,your/file/path/rawdata/Wild-1.2.fastq,reverse
Wild-2,your/file/path/rawdata/Wild-2.1.fastq,forward
Wild-2,your/file/path/rawdata/Wild-2.2.fastq,reverse
Wild-3,your/file/path/rawdata/Wild-3.1.fastq,forward
Wild-3,your/file/path/rawdata/Wild-3.2.fastq,reverse

qiime import コマンドを実行

qiime tools import \
--type 'SampleData[PairedEndSequencesWithQuality]' \
--input-path ManifestFile.csv \
--output-path qza/demux.qza \
--input-format PairedEndFastqManifestPhred33
オプションの説明
Usage: qiime tools import [OPTIONS]

  Import data to create a new QIIME 2 Artifact. See https://docs.qiime2.org/
  for usage examples and details on the file types and associated semantic types that can be imported. 
Semantic types — QIIME 2 2023.5.1 documentation
Options: 
--type TEXT  
The semantic type of the artifact that will be created upon importing. Use --show-importable-types to see what importable semantic types are available in the current deployment.[required] 
SampleData[SequencesWithQuality]: シングルエンドデータはこっち
SampleData[PairedEndSequencesWithQuality]: ペアエンドデータはこっち

--input-path PATH      
Path to file or directory that should be imported. [required] 
Manifest file があるならManifest file へのパス

--output-path ARTIFACT  
Path where output artifact should be written. [required] 
保存するときのファイル名と、保存先のパス

--input-format TEXT   
The format of the data to be imported. If not provided, data must be in the format expected by the semantic type provided via --type.
SingleEndFastqManifestPhred33：シングルエンドかつ品質スコアに使用される PHRED オフセット※が33 の場合
PairedEndFastqManifestPhred33：ペアエンドかつ品質スコアに使用される PHRED オフセット※が33 の場合
※FASTQ ファイル内に英数字で書かれたクオリティスコアを数値に変換するには、その英数字の ASCII コードを調べ、その ASCII コードから 33 もしくは64 を引く計算を行うが、その際の数値がPHRED オフセット
※※現在多くのシーケンサーは33（参考：クオリティスコア | シークエンサーが出力するクオリティスコアとシーケンシングエラー率)

--show-importable-types 
Show the semantic types that can be supplied to --type to import data into an artifact. 

--show-importable-formats 
Show formats that can be supplied to --input-format to import data into an artifact.

配列の精度を確認し、DADA2 によるデノイジング時の除去配列長を決定

まず、プライマー配列のカット

qiime cutadapt trim-paired

公式：https://docs.qiime2.org/2023.2/plugins/available/cutadapt/trim-paired/

切るときに、長さで除去範囲を指定する方法もありますが、配列で指定するオプションを使いました。

--p-front-f TEXT...

forward 配列の5' 末端からTEXTに代入した文字列を切り取る。

Sequence of an adapter ligated to the 5' end. The List[Str] adapter and any preceding bases are trimmed. Partial matches at the 5' end are allowed. If a ^ character is prepended, the adapter is only found if it is at the beginning of the read. Search in forward read.

--p-front-r TEXT...

reverse 配列の5' 末端からTEXTに代入した文字列を切り取る。

Sequence of an adapter ligated to the 5' end. The List[Str] adapter and any preceding bases are trimmed. Partial matches at the 5' end are allowed. If a ^ character is prepended, the adapter is only found if it is at the beginning of the read. Search in reverse read.

実行したスクリプト

qiime cutadapt trim-paired \
--i-demultiplexed-sequences ~/path/to/importした配列/demux.qza \
--p-cores 5 \
--p-front-f TCG（中略）CAG \
--p-front-r GTCT（中略）ATCC \
--o-trimmed-sequences qza/primer-trimmed-demux-paired-end.qza

うまくいくと以下の文章が表示されます。

Saved SampleData[PairedEndSequencesWithQuality] to: qza/primer-trimmed-demux-paired-end.qza

疑問：今回複数種類のプライマーを使ったので、複数回qiime cutadapt trim-paired コマンドを使ったのですが、あれでよかったのでしょうか…。わかる人がいれば教えてください。

クオリティチェック

以下のコマンドでサンプルごとのリードカウントをまとめます。

qiime demux summarize \
--i-data qza/demux.qza \
--o-visualization qzv/demux.qzv

可視化ファイル [.qzv] を開く（Qiime2 View にドラック＆ドロップ）

参考：https://note.com/nanaimo_/n/n0a7234c90ee6#281f2705-07bf-4314-ab51-34130fb21c7f

Qiime2 View https://view.qiime2.org/

Visualization > Interactive Quality Plot

上手くいけば以下のようなプロットが表示されます。

Interactive Quality Plot

 

デノイジング時の除去配列長を決定

以下の動画を参考に行いました。

www.youtube.com

このチュートリアルでは25% のデータがクオリティスコア30 を下回る配列を切り取ると言っていますが、わたしのデータでそれをやると失うものが大きすぎたので、わたしはクオリティスコア20 で行いました。

Interactive Quality Plot

Interactive Quality Plot 上で拡大して見たい箇所を、タッチパッドを長押ししながら囲うと拡大できます。さらに、棒グラフにカーソルを乗せると、プロット下のテーブルにその長さのリードのクオリティ情報の詳細が表示されます。上の図では295（The plot at position 295）の位置のクオリティ情報の詳細が出ています。この位置では25th のクオリティスコアは22 なので、クオリティスコア20 まで残す場合にはこの位置は切らずに済みます。

こんな風にして、どの長さ以上から25 th のクオリティスコアが20 を切るのか探し、Foward とReverse についてそれぞれメモしておきます。

わたしの場合はFoward リードが277、Reverse リードが158 からクオリティ20 を切りました。

 

ここまで終わったら、いよいよデノイジングです。

DADA2 によるデノイジング

denoise-paired: Denoise and dereplicate paired-end sequences

参考：denoise-paired: Denoise and dereplicate paired-end sequences — QIIME 2 2023.5.1 documentation

qiime dada2 denoise-paired \
--i-demultiplexed-seqs qza/primer-trimmed-demux-paired-end.qza \
--p-trunc-len-f 277 \
--p-trunc-len-r 158 \
--p-trim-left-f 0 \
--p-trim-left-r 0 \
--p-trunc-q 2 \
--p-n-reads-learn 1000000 \
--p-max-ee-f 2.0 \
--p-max-ee-r 2.0 \
--p-n-threads 4 \
--o-table qza/table.qza \
--o-denoising-stats qza/stats.qza \
--o-representative-sequences qza/rep-seqs.qza
オプションの説明

オプションのより詳しい説明は以下のQiita 記事がおすすめです。

Qiime2 を用いた 16S rRNA 菌叢解析 - Qiita

--i-demultiplexed-seqs input_path/primer-trimmed-demux-paired-end.qza(input)
--p-trunc-len-f 277 # Foward 配列をどの長さで切り捨てるか（先ほどクオリティスコアで決めた位置）。 0 にすると働かない。
--p-trunc-len-r 158 # Reverse 配列をどの長さで切り捨てるか（先ほどクオリティスコアで決めた位置）。 0 にすると働かない。
--p-trim-left-f 0 # Foward 配列を5'側からどの程度切り捨てるか。今回はプライマー除去を先に行っているので不要。デフォルトは0。
--p-trim-left-r 0 \ Reverse 配列を5'側からどの程度切り捨てるか。今回はプライマー除去を先に行っているので不要。デフォルトは0。
--p-trunc-q 2 #リードの開始点からこのクオリティスコア以下の配列は捨てられる。デフォルトは2。
--p-n-reads-learn 1000000 # error model の学習に使われるリードの数。デフォルトは1000000。
--p-max-ee-f 2.0 # 予想されるエラーの数がこの値よりも多いFoward リードは破棄される。デフォルトは2.0。
--p-max-ee-r 2.0 # 予想されるエラーの数がこの値よりも多いReverse リードは破棄される。デフォルトは2.0。
--p-n-threads 4 # 解析に使用するスレッド数
--o-table qza/table.qza # FeatureTable[Frequency] アウトプットのパスとアウトプットファイル名の指定
--o-denoising-stats qza/stats.qza # SampleData[DADA2Stats] アウトプットのパスとアウトプットファイル名の指定
--o-representative-sequences qza/rep-seqs.qza # FeatureData[Sequence] アウトプットのパスとアウトプットファイル名の指定

上手くいくと以下の文章が表示されます。

Saved FeatureTable[Frequency] to: qza/table.qza Saved FeatureData[Sequence] to: qza/rep-seqs.qza Saved SampleData[DADA2Stats] to: qza/stats.qza

デノイジング結果の可視化

qiime metadata tabulate \
  --m-input-file qza/stats.qza \
  --o-visualization stats.dada.qzv


qzv ファイルをQiime2 View https://view.qiime2.org/に投げ込みます。

DADA2 結果の可視化（stats.dada.qzv）
結果の説明
  • input：開始リード数
  • numeric filtered：DADA2 の最初のフィルタリングを通過したリード数
  • numeric percentage of input passed filter：DADA2 の最初のフィルタリングを通過したリード数のinput に対するパーセンテージ
  • numeric denoised：DADA2 によるデノイジングを通過したリード
  • numeric merged：マージ可能だったリード数
  • numeric percentage of input merged：開始リードの中で、マージ可能だった割合
  • numeric non-chimeric：保持されたノン-キメラリード。PCR が上手くいっていないと低くなる。
  • percentage of input non-chimeric：保持されたノン-キメラリードのinput に対する割合

参考：Denoising sequence data with DADA2 - YouTube(13:00~)

疑問：percentage of input non-chimeric ってどの程度残っていればOK なんだろう？今回は41~67% だったので、やや不安です。

DADA2 結果のさらなるフィルタリング

Total-frequency-based filtering：https://docs.qiime2.org/2023.2/tutorials/filtering/

以下のコマンドではリード数が1500 を切るサンプルを取り除きます。

qiime feature-table filter-samples  \
--i-table qza/table.qza  \
--p-min-frequency 1500 \
--o-filtered-table selected_table.dada.qza