微生物群集解析ツール QIIME 2 を使う．Part2.（QC，FeatureTable生成編）

# QIIME 2 を使う．
[QIIME2](https://docs.qiime2.org/) はDNAのシーケンスデータから微生物解析を行うためのオープンソースのパイプラインである．クオリティの高いグラフや統計の処理を行うことが可能である．

今回は，前回の[微生物群集解析ツール QIIME 2 を使う．Part1.（インストール，ファイルインポート編）] (https://qiita.com/drafts/5d29bd8bbe1d19fdeec5/edit)の続きで，fastq ファイルのクオリティコントロール（QC）とFeature Table の作成を行う．

## Fastq ファイルのクオリティコントロール（QC）を行う．

cutadapt にて、アダプター配列の除去。
341f (--p-front-f CCTAYGGGRBGCASCAG)と806r (--p-front-r GGACTACNNGGGTATCTAAT)を想定。

```bash
# cutadapt による341f (--p-front-f CCTAYGGGRBGCASCAG)と806r (--p-front-r GGACTACNNGGGTATCTAAT)のプライマー配列除去	
qiime cutadapt trim-paired --i-demultiplexed-sequences demux-paired-end.qza --p-cores 4 --p-front-f CCTAYGGGRBGCASCAG --p-front-r GGACTACNNGGGTATCTAAT --o-trimmed-sequences primer-trimmed-demux-paired-end.qza --verbose

# Summarise the reads
qiime demux summarize --i-data primer-trimmed-demux-paired-end.qza --o-visualization primer-trimmed-demux-paired-end.qzv
```

[DADA2](https://www.ncbi.nlm.nih.gov/pubmed/27214047)はIllumina のアンプリコン配列を検出し修正するパイプラインである．Phix配列やキメラ配列をフィルタリングしてくれる．

公式チュートリアルでは，`dada2 denoise-single`になっているが，今回のデータはpaired-end のデータなので`dada2 denoise-paired`で実行する．
[deblur error "Argument to parameter 'demux' is not a subtype of SampleData[SequencesWithQuality]."](https://forum.qiime2.org/t/deblur-error-argument-to-parameter-demux-is-not-a-subtype-of-sampledata-sequenceswithquality/1105)

`dada2 denoise-paired` メソッドでは，4つのパラメータを与える．

`--p-trim-left-f x`：各シーケンスのforward 配列の初めから塩基長xまでを削除する．
`--p-trim-left-r y`：各シーケンスのReverse 配列の初めから塩基長yまでを削除する．

出力されたクオリティのプロットを見て、クオリティが20以下の部分をカットしつつ。
**採用する最大配列長には結合したペアエンドリードが十分重なる長さを指定することが求められる．**

`--p-trunc-len-f m`：各シーケンスのforward 配列から採用する最大配列長mを決定する． mを超える塩基長は除去する．
`--p-trunc-len-r n`：各シーケンスのreverse 配列から採用する最大配列長nを決定する． nを超える塩基長は除去する．

**v3-v4領域だと460bpなので、最低でも 20bp 以上のオーバーラップ領域を加味して、270,210とか480を超える長さで実行する．**

`--p-n-threads k`：解析に使用するスレッド数を指定する．

**プライマー配列は取り除く必要があるので、先頭側のカット長はプライマー配列長を通常指定する．**

```bash
qiime dada2 denoise-paired \
  --i-demultiplexed-seqs primer-trimmed-demux-paired-end.qza \
  --p-trim-left-f 0 \
  --p-trim-left-r 0 \
  --p-trunc-len-f 280 \
  --p-trunc-len-r 220 \
  --o-representative-sequences rep-seqs-dada2.qza \
  --o-table table-dada2.qza \
  --o-denoising-stats stats-dada2.qza \
  --p-n-threads 30
# 名前を変更する（別にしなくても良いがqiime2のチュートリアルでの汎用性を考えてここでは変更しておく）
mv rep-seqs-dada2.qza rep-seqs.qza
mv table-dada2.qza table.qza

```

## Feature Table の生成とサマライズを行う．
クオリティのフィルタリングが完了したら，結果を調べるためにFeature Table を生成する．

`feature-table summarize`:各サンプルに関連付けられているシーケンスの数、各特徴、それらの分布のヒストグラムを出力する．
`feature-table tabulate-seqs`:Feature ID のシーケンスへのマッピングとNCBI database に対するBLAST検索のリンクを提供する．

```bash
# metadata ファイルをダウンロードする
wget -O "sample-metadata.tsv" "https://data.qiime2.org/2017.12/tutorials/moving-pictures/sample_metadata.tsv"
#or
curl -sL "https://data.qiime2.org/2017.12/tutorials/moving-pictures/sample_metadata.tsv" > "sample-metadata.tsv"

# サマライズの実行
qiime feature-table summarize \
  --i-table table.qza \
  --o-visualization table.qzv \
  --m-sample-metadata-file sample-metadata.tsv

qiime feature-table tabulate-seqs \
  --i-data rep-seqs.qza \
  --o-visualization rep-seqs.qzv
```

## 補足
DAD2の公式ドキュメントにおいて、「[How can I remove primers/adapters/etc from my amplicon reads?](https://benjjneb.github.io/dada2/faq.html)」という項目にて、重要なことが書かれている。

プライマー配列の除去について、より複雑な状況（例えば、不均一性スペーサー、ITS領域などの可変長アンプリコン）の場合は、cutadaptツールまたはtrimmomaticツールをお勧めします。
プライマー配列が残っていると後の解析に影響してしまうので、しっかり除去しよう。

> If primers are at the start of your reads and are a constant length (a common scenario) you can use the trimLeft = c(FWD_PRIMER_LEN, REV_PRIMER_LEN) argument of dada2’s filtering functions to remove the primers.

> **For more complex situations (e.g. heterogenity spacers, variable length amplicons such as the ITS region) we recommend the cutadapt tool or the trimmomatic tool.** You can see a worked example, including verification of primer removal, in the ITS-specific version of the DADA2 workflow.

> **Please double-check that your primers have been removed! The ambiguous nucleotides in unremoved primers will interfere with the dada2 pipeline.**