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

# 概要

本稿では、菌叢解析ソフト [**Qiime2（2019.7 ver.）**](https://docs.qiime2.org/2019.7/)を用いて、細菌の系統分類マーカーである [16S rRNA 遺伝子（16S rDNA）](https://bifidus-fund.jp/keyword/kw009.shtml)のアンプリコン（PCR増幅産物）から、微生物群集構造を解析する方法**（16S アンプリコン解析）**を紹介する。本稿で紹介する解析フローやコマンドは一部を除いて別バージョンの Qiime2 でも共通である。（[2020.2 版](https://docs.qiime2.org/2020.2/)、[2020.6 版](https://docs.qiime2.org/2020.6/)、[2020.8 版](https://docs.qiime2.org/2020.8/)、[2021.2 版](https://docs.qiime2.org/2021.2/)では確認済。）
なお、系統分類マーカーだけでなくゲノム全体を対象にした[ショットガンメタゲノム解析](https://qiita.com/keisuke-ota/items/4e34ceffb33f2e983f09)は本稿の対象ではない。

# 菌叢解析手順

以下に、16S アンプリコンを用いた菌叢解析手順を 4 段階でまとめた。アンプリコン解析で得られるデータセットは、試料から抽出されたサンプルの read 数であって、試料に含まれる細菌の絶対的な存在量ではない。また、シーケンサーの性能や PCR の精度によって、サンプルの合計 read 数は変わる。したがって、サンプルは細菌の**組成**に基づいて比較されるべきである。

1. 試料からサンプルとして DNA を抽出する。
2. DNA をシーケンシングする。
3. 「塩基配列」と「各塩基配列の read 数」を得る。
4. read 数から、各サンプルの細菌組成を比較する。

## Qiime2 による解析の工程

16S アンプリコン解析において、Qiime2 では以下の工程で解析を実施する。

<img width="808" alt="スクリーンショット 2020-06-11 15.11.41.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/05ffd8dd-c64e-579d-acb8-811850f59bc8.png">

まずアンプリコン配列は、シーケンシング時に PCR で生じたエラーを修正される。ここで「塩基配列」と「各塩基配列の read 数」が得られる。これを、塩基配列のまま比較することも可能だが、データベース上の細菌に照合してから比較したり、細菌がもつ酵素の量を予測して比較したりすることもできる。比較は、**組成データ解析（Compositional data analysis:CoDA）**によって行われる。

### 各解析の概要

* Taxonomy analysis
16S rDNA のアンプリコンシーケンスをデータベース上の細菌に照合させる分析。どのような細菌が存在するのかを明らかにする。

* PICRUSt2 analysis
群集中の細菌がコードしている酵素の存在量を予測する解析。専用ツールの PICRUSt2 はデフォルトの Qiime2 にはインストールされないため、 別個にプラグインする必要がある。ただし、最新版の Qiime2 には対応していないことがあるため、本稿では 2019.７ 版をインストールしている。

* Differential analysis
群間で、特徴量（細菌や酵素）の存在量を比較する多変量解析。ANCOM、Balance、[ALDEx2](https://qiita.com/keisuke-ota/items/b58f9a28a5ecd53b8ecb) などがある。存在量に差がある特徴量の索敵を目的とする。

* Diversity analysis
群間で、特徴量（細菌や酵素）の多様性を比較する多変量解析。[α-多様性とβ-多様性](https://qiita.com/keisuke-ota/items/9701c583e153df467c05)がある。

# テストデータ

本稿では [IBD multi'omics database ](https://ibdmdb.org/tunnel/public/HMP2/16S/1806/rawfiles)  (IBDMDB) に公開されている腸生検の 16S rDNA アンプリコンシーケンスをテストデータとして採用し、クローン病（CD）患者、潰瘍性大腸炎（UC）患者、非炎症性腸疾患（nonIBD）被験者における腸内細菌叢および菌叢代謝を比較した。このデータセットには合計  178 サンプルが登録されている。

各データは、File Name（206172.tarなど）をクリックするとダウンロードできる。この記事を参考に動作確認をしてみたい場合は、いくつかのサンプルをダウンロードしてみるとよい。

![スクリーンショット 2020-11-30 9.50.50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/8ba3349b-c409-8a54-1f79-0c528235abf5.png)

ダウンロードした tar ファイルを展開すると .bz2 に圧縮されている 2つの FASTQ ファイル（拡張子.fq）が得られる。解析のため、.bz2 も展開しておこう。
![スクリーンショット 2020-11-30 10.02.22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/8e85cc66-a9dc-ab80-676e-a89255c48ccb.png)

末尾が 「.1.fq」 のファイルは forward reads で、末尾が「.2.fq」のファイルは reverse reads である。
## FASTQ 形式

テストデータには、シーケンシングによって得られた read の塩基配列と品質情報が **FASTQ 形式** で記載されている。各 read の情報は以下の 4 行でまとめられている。

* 1行目：ID
* 2行目：塩基配列
* 4行目：品質情報

> @M70287:115:000000000-B59YV:1:1101:18877:1845 1:N:0:
TACGGAGGATCCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGAGCGTAGGTGGACTGGTAAGTCAGTTGTGAAAGTTTGCGGCTCAACCGTAAAATTGCAGTTGATACTGTCAGTCTTGAGTACAGTAGAGGTGGGCGGAATTCGTGGTGTAGCGGTGAAATGCTTAGATATCACGAAGAACTCCGATTGCGAAGGCAGCTCACTGGACTGCAACTGACACTGATGCTCGAAAGTGTGGGTATCAAAC
+
>1>111111B>BF1AEEGGGGGHGG0EEHHHHHHGGFGHHHHGG/EEEGGGGBGG/1BF/GGDFHFGHHHHHGFFHBGHEBBEE>EEGFFH?ECEGBFHHFB1BBBGGGHHFHBDGHFHGGGFHH1GHHDG2DHHHGGEGGGGFHGFGGGHGHHHHGGGGGHHHHHHHHHHCHHHG09:AGGGGBFGGGGGGGGGGGFFFFFFFFFB/-BFFEBBBFFFFFFEFFFBFFFBF?FFF/;BFBFFBFFFFF#

品質情報は、各塩基が PCR で生じたエラー率 $P$ から計算された Q 値を指す。Q 値は以下のように計算される。

$$
Q =-10 \times log_{10}P
$$

FASTQ 形式では、Q 値に対応する1つの文字が、各塩基に割り当てられて記録されている。

# Qiime2 の成果物（artifact）

Qiime2 には、生データからインポートされた中間成果物（qzaファイル）と、それをブラウザに表示できるように変換した可視化成果物（qzvファイル）がある。

<img width="796" alt="スクリーンショット 2020-06-11 15.58.17.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/25b71551-986d-8b76-5b37-1e7e43ce8de5.png">

本稿の解析過程で作られる qza ファイルの中で、重要なものを以下にまとめた。ただし、PICRUSt2 analysis による成果物は除いている。

<img width="787" alt="スクリーンショット 2020-06-11 16.05.28.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/5a885939-d520-2ca9-e136-9dec9f2a3c3b.png">

## 各ファイルの説明

* demux.qza
 各サンプルのFASTQ データをまとめたファイル。

* rep-seqs.qza
 データセット中に存在する read の塩基配列をまとめたファイル。識別用の ID も割り当てられている。

* table.qza
 データセット中に存在する塩基配列の read 数をまとめたファイル。各塩基配列は ID によって区別されており、塩基配列の中身は記録されていない。  

* silva-132-99-nb-classifier.qza
  16S rDNA の塩基配列をまとめたデータベース。

* taxonomy.qza
　rep-seq.qza の塩基配列が、データベースにおけるどの細菌に属するのかをまとめたファイル。



# Qiime2 のインストール

## 環境

macOS Catalina 10.15.1

## conda にインストール
 
既に anaconda もしくは miniconda をインストールしている場合、Qiime2 (2019.7.ver.)の仮想環境を以下のコマンドで conda に構築することができる。公式サイトは[こちら](https://docs.qiime2.org/2019.7/install/native/)。

```terminal:ターミナル
wget https://data.qiime2.org/distro/core/qiime2-2019.7-py36-osx-conda.yml
conda env create -n qiime2-2019.7 --file qiime2-2019.7-py36-osx-conda.yml
rm qiime2-2019.7-py36-osx-conda.yml
```

conda の仮想環境を activate する。

```terminal:ターミナル
source activate qiime2-2019.7
```

qzvファイルは次のコマンドでブラウザに表示できる。

```terminal:ターミナル
qiime tools view ファイル名.qzv
```

## docker にインストール

miniconda や anaconda をインストールしていない場合は、docker で Qiime2 の仮想環境を構築することもできる。まず、docker image をダウンロードする。

```terminal:ターミナル
docker pull qiime2/core:2019.7
```

次のコマンドで、docker image からコンテナを作成し、Qiime2 を起動する。

```terminal:ターミナル
docker run -t -i -v $(pwd):/data qiime2/core:2019.7 /bin/bash
```

ただし、docker では ``qiime tools view ファイル名.qzv`` というコマンドを使用できず、ローカルで qzv ファイルをブラウザに表示させる手段がない。そこで、ローカルホストに可視化用のサーバー（Qiime2 viewer）を構築する。
まず、GitHub から[Qiime2 viewer](https://github.com/qiime2/q2view)をダウンロードする。

```terminal:ターミナル
git clone https://github.com/qiime2/q2view.git
```

次に、ファイルのルートディレクトリで次のコマンドを実行すると、qiime2 viewer が build される。

```terminal:ターミナル
npm install
npm run build
```

そして、次のコマンドで qiime2 viewer を起動する。

```terminal:ターミナル
npm start
```

この状態で http://localhost:4242/ にアクセスすると、ブラウザに viewer が表示される。

# 生データのインポート

次のように、ダウンロードした生データを１つのフォルダにまとめておこう。
![スクリーンショット 2020-11-30 10.18.31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/b6f1b936-1bfb-ebcd-8929-a8c18a32ddb9.png)

今回は、hmp2 というフォルダに生データをまとめたので、hmp2 フォルダで以下のコマンドを実行すると、``Manifest.csv`` が得られる。

```terminal:ターミナル
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
```
```:Manifest.csv
sample-id,absolute-filepath,direction
206534,/Users/home/hmp2/206534/HMP.2.206534.1.fq,forward
206534,/Users/home/hmp2/206534/HMP.2.206534.2.fq,reverse
206536,/Users/home/hmp2/206536/HMP.2.206536.1.fq,forward
206536,/Users/home/hmp2/206536/HMP.2.206536.2.fq,reverse
206538,/Users/home/hmp2/206538/HMP.2.206538.1.fq,forward
206538,/Users/home/hmp2/206538/HMP.2.206538.2.fq,reverse
（後略）
```

``awk``コマンドは、取得したレコードの各行に対して指定の処理を実行できる。
詳細を以下にまとめた。

| オプションなど  | 機能  |
| --- | --- |
| -F | レコードの区切り記号を指定  |
| Begin| awk コマンドによる処理を実行する前の命令を指定 |
| if(NR%2==0)  | 行数が 2 で割り切れる行でのみ実行することを指定  |
| else | if( ) が成立しない行でのみ実行することを指定 |
| print|  Manifest.csv に出力するテキストを指定 |
| $0|  テキストファイル全体を表示  |
| $3|  取得したテキストファイルの 3 区画目を表示 |

``ManifestFile.csv ``をもとにインポートする。

```terminal:ターミナル（生データがペアエンドシーケンスの場合）
qiime tools import \
--type 'SampleData[PairedEndSequencesWithQuality]' \
--input-path ManifestFile.csv \
--output-path demux.qza \
--input-format PairedEndFastqManifestPhred33
```
各オプションに代入する値を以下にまとめた。PHRED は、エンコードの方式を指す。FASTQ データに記載されている品質情報が、Q 値に 33 を加えた値であれば 33、64 を加えた値であれば 64 を選択する。両方を試してみて、解析エラーが表示されないモードを利用するとよい。

| Read type   | PHRED  | --type に代入   | --input-format  に代入 |
| ----------- | ------------ | ------------- | ---------------------- |
| Single Ends |      33      | 'SampleData[SequencesWithQuality]'|SingleEndFastqManifestPhred33|
| Single Ends |      33      | 'SampleData[SequencesWithQuality]'|SingleEndFastqManifestPhred64|
| Pair Ends   |      64      | 'SampleData[PairedEndSequencesWithQuality]'| PairedEndFastqManifestPhred33|
| Pair Ends   |      64      | 'SampleData[PairedEndSequencesWithQuality]'| PairedEndFastqManifestPhred64|

[公式サイト](https://docs.qiime2.org/2019.7/tutorials/importing/)では、``--input-format`` に``PairedEndFastqManifestPhred33V2``を代入しているが、Manifest.csv の形式が異なるだけである。


# 生データの denoising

生データには、PCR 増幅時に生じたエラー配列が含まれている。まずは生データのクオリティスコア（Q値）を可視化して、配列の精度を確認してみる。

```terminal:ターミナル
qiime demux summarize \
--i-data demux.qza \
--o-visualization demux.qzv
```

ペアエンドシーケンスの場合、forward reads と reverse reads の Q 値がそれぞれ表示される。一般に Q 値は 3'末端から5'末端にかけて低下し、さらに reverse reads の方が低い。PCR 触媒は forward reads の 3'末端から順番に消耗していくため、徐々に精度は低下する。
![forward.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/f487998c-b137-1940-f8e8-b8ef7c1ad45e.png)

![reverse.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/371230d4-2598-036e-a520-d11ddf36b198.png)

アンプリコンには、真の配列以外に PCR で生じるエラーが含まれている。以下の概念図では、各配列の相同性と存在量を円の位置と大きさで表している。
<img width="543" alt="スクリーンショット 2020-04-05 10.57.17.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/fe9ca6bc-481c-7ec8-974b-becf419e15e6.png">

## OTU
16S rDNA のアンプリコン配列が一定の閾値以上で一致する集団を OTU という。PCR のエラー配列と真の配列を区別できないことがある。
<img width="674" alt="スクリーンショット 2020-04-05 11.00.12.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/c40d529b-eeb5-4cea-dfaa-a47065723b3b.png">

そこで、PCRのエラー配列を修正あるいは除去し、真の配列から得られる細菌の分類単位を sub-OTU といい、従来の OTU と区別する。sub-OTU は **DADA2** もしくは **Deblur** というパッケージで作成できる。特に DADA2 で得られる sub-OTU は、**ASV (amplicon sequence varient)** と呼ばれる。

## DADA2 

DADA2 では、エラー配列を真の配列に修正し、加算している。どのように修正するかをエラーモデルの学習によって推定している。

<img width="560" alt="スクリーンショット 2020-05-07 10.40.18.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/dfd9f347-a733-8da4-3559-bd612e37a081.png">


```terminal:ターミナル（ペアエンドシーケンスの場合）
qiime dada2 denoise-paired \
--i-demultiplexed-seqs demux.qza \
--p-trunc-len-f 230 \
--p-trunc-len-r 220 \
--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 table.qza \
--o-denoising-stats stats.qza \
--o-representative-sequences rep-seqs.qza \
--verbose
```

* ``--p-trunc-len-f`` ``--p-trunc-len-r``

マージに使用する read の長さを決める。demux.qzv を確認し、Q 値が低下する後半の配列を除去したい場合は短く設定すると良い。ただし forward reads と reverse reads がマージするには最低でも 20bp 以上のオーバーラップ領域が必要とされている。例えば 255bp の V4 領域は合計 275 bp なければマージされない。V3 領域も含めた生データには、最低でも合計 500 bp 程度が必要である。

* ``--p-trim-left-f`` ``--p-trim-left-r``

3'末端側で除去する read 長を指定する。アダプター配列などを除去する場合に用いる。

* ``--p-trunc-q ``

Read の始まりから 1 塩基ずつ Q 値を調べ、この値より Q 値が高い最初の塩基を見つける。その塩基より先の配列を全てトリミングする。したがって、クオリティスコアの低い read の場合、この値が**高い**と塩基長が足りずにマージしなくなることがある。

この値より Q 値が低い塩基は ``--p-trunc-len-f`` や ``--p-trunc-len-r`` で設定した塩基長に到達するより前に見つかることもある。その場合 ``--p-trunc-len-f`` や ``--p-trunc-len-r`` で設定した塩基より短くトリミングされる。

* ``--p-max-ee-f`` ``--p-max-ee-r``

Read の始まりから 1 塩基ずつ P 値（Q 値ではない）を加算し、この値を超えた塩基より先の配列をトリミングする。``--p-trunc-q `` と対照的に、クオリティスコアの低い read の場合、この値が**低い**と塩基長が足りずにマージしなくなることがある。

品質の低い塩基が続くと ``--p-trunc-len-f`` や ``--p-trunc-len-r`` で設定した塩基長に到達するより前に見つかることもある。その場合 ``--p-trunc-len-f`` や ``--p-trunc-len-r`` で設定した塩基より短くトリミングされる。


* ``--p-n-reads-learn``

エラーモデルの学習に要する read 数を指定する。小さいと処理が早く終わるがエラーモデルの精度が落ちる。

* ``--p-n-threads ``

スレッド数を設定する。高いと処理が早く終わるが、それまで他の動作が鈍くなる。

* `` --o-table ``

各配列の存在量データを出力する。塩基の並びは含まれない。

* `` --o-representative-sequences ``

各配列の塩基データを出力する。

* `` --o-denoising-stats ``

DADA2 の各工程で入力された read 数を出力する。

<img width="1041" alt="スクリーンショット 2020-04-05 11.46.41.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/5de48e71-eb81-b98c-eb1c-0fda1bccfdca.png">

| 項目 | 内容 |
| -------- | -------- |
| Input    | 入力された read 数|
| Filtered | トリミング後の read 数。Q 値の基準が厳しいほど減少する。|
| Denoised |Denoising 後の read 数。エラーモデルに不寛容であるほど減少する。 |
| Merged   |マージ後の read 数。オーバーラップ領域が 20 bp 以下だと極端に減少する。|
| Non-chimeric | 最終生成 read 数。PCR の精度が低いと減少する。|

## Deblur

Deblur にはペアエンドシーケンスのマージが実装されていないため、事前に Vsearch で merge させておく。
トリミング長は自動で設定される。

```terminal:ターミナル
qiime vsearch join-pairs   \
--i-demultiplexed-seqs demux.qza   --o-joined-sequences demux-joined.qza
```

Deblur では、エラー配列を修正せずに除去する。DADA2 より配列の精度は高いが得られる read 数は少ない。
<img width="603" alt="スクリーンショット 2020-05-07 10.40.32.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/554f63b7-7032-d6af-fb4e-724c0f9940d6.png">


```terminal:ターミナル
qiime deblur denoise-16S \
--i-demultiplexed-seqs demux-joined.qza \
--p-trim-length 250 \
--p-mean-error 0.005 \
--p-indel-prob 0.01 \
--p-indel-max 3.0 \
--p-jobs-to-start 36 \
--o-table table.qza \
--o-stats stats.qza\
--o-representative-sequences rep-seqs.qza \
--verbose
```

* ``--p-trim-length``

read 長を設定する。DADA2 と異なり、ここで指定した read 長しか出力されない。

* ``--p-mean-error``

1塩基あたりの平均塩基数を設定する。

* ``--p-indel-prob``

1塩基あたりの最大インデル（indel）数を設定する。エラー配列の挿入（insertion）と欠失（deletion）を略してインデルという。

* `` --o-table `` `` --o-representative-sequences `` `` --o-denoising-stats ``

出力ファイルはいずれも DADA2 と共通である。

## 中間成果物の確認

DADA2 または Deblur を通して、微生物の存在量データ（カウントデータ）は``table.qza``に出力され、塩基配列データは``rep-seqs.qza``に出力される。
両ファイルのデータは、``sampleid``によって紐づけられている。



# メタデータの作成

次のように、メタデータファイルを tsv 形式で作成する。Excel などの表計算ソフトで作成したデータシートを tsv  形式にエクスポートすると簡単に作れる。
``sampleid``は一番左のカラムに用意すること。詳しくは[こちら](https://docs.qiime2.org/2020.2/tutorials/metadata/)。

```:Metadata.tsv
sampleid	Participant ID	Age at diagnosis	Education Level	Occupation	biopsy_location	consent_age	diagnosis	race	sex	read
206534	M2008	30	Bachelor's degree	Paid	Ileum	30	CD	White	Female	33709.0
206536	M2008	30	Bachelor's degree	Paid	Transverse colon	30	CD	White	Female	28684.0
206538	M2008	30	Bachelor's degree	Paid	Rectum	30	CD	White	Female	15375.0
（後略）
```
![スクリーンショット 2020-10-31 14.25.20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/cb5295f4-feb6-d4b8-20ba-38439302d1b8.png)

# サンプルの抽出
多変量解析には、read 数が足りず環境中の菌叢を正しく反映しているとはいえないサンプルは適さないため、除外する必要がある。
例えば、read 数 1500 以上のサンプルの抽出は、次のコマンドで実施できる。

```terminal:ターミナル
qiime feature-table filter-samples  \
--i-table table.qza  \
--p-min-frequency 1500 \
--o-filtered-table selected_table.qza
```

さらに``selected_table.qza``は、``Metadata.tsv`` における特定の条件を満たすサンプルのみ抽出することができる。抽出された``table.qza``も、全体データから得られた塩基データ``rep-seqs.qza``や細菌分類データ``taxonomy.qza``（後述）と対応している。

```terminal:ターミナル（男性データのみ抽出したい場合）
qiime feature-table filter-samples \
--i-table selected_table.qza \
--m-metadata-file Metadata.tsv \
--p-where "[sex]='Male’” \
--o-filtered-table Male/table.qza \
--output-dir Male
```

```terminal:ターミナル（診断時の年齢が30歳未満の男性データを抽出したい場合）
qiime feature-table filter-samples \
--i-table selected_table.qza \
--m-metadata-file Metadata.tsv \
--p-where " [Age at diagnosis]<30 AND [sex]='Male'” \
--o-filtered-table DisagnosisU30Male/table.qza \
--output-dir DisagnosisU30Male
```

# Taxonomy analysis

[公式サイト](https://docs.qiime2.org/2020.2/data-resources/)から、以下のうちいずれかのデータベースをダウンロードする。

* Silva 132 99% OTUs full-length sequences
* Silva 132 99% OTUs from 515F/806R region of sequences
* Greengenes 13_8 99% OTUs full-length sequences
* Greengenes 13_8 99% OTUs from 515F/806R region of sequences

Silva は現在も更新し続けているデータベースである。生データが V3 または V4 領域の場合は ``515F/806R region of sequences`` を使用し、他の領域の配列であれば ``full-length sequences`` を使用する。
Greengenes は現在更新されていないデータベースであり、新種の細菌は登録されていない。OTU で細菌を分類していた頃に使われていた。
特別な理由がなければ Silva をダウンロードするとよい。

塩基データ``rep-seqs.qza``から細菌分類データ``taxonomy.qza``を作成する。

```terminal:ターミナル
qiime feature-classifier classify-sklearn  \
--i-classifier silva-132-99-nb-classifier.qza  \
--i-reads rep-seqs.qza  \
--o-classification taxonomy.qza 
```

```terminal:ターミナル
qiime metadata tabulate \
--m-input-file taxonomy.qza \
--o-visualization taxonomy.qzv
```
存在量データ``table.qza``と合わせて棒グラフを得る。

```terminal:ターミナル
qiime taxa barplot   \
--i-table table.qza   \
--i-taxonomy taxonomy.qza   \
--m-metadata-file Metadata.tsv   \
--o-visualization taxa-bar-plots.qzv
```

``taxa-bar-plots.qzv`` の上部で棒グラフの表示を設定できる。各内容について以下で説明する。

![スクリーンショット 2020-10-31 14.59.54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/4b77a7ee-de08-0edb-30f3-17bc6f01309b.png)


* Taxonomic Level

以下の表に対応している。なお、塩基データ``rep-seqs.qza``の分類は sub-OTU に基づくものであり、Species より細かい。

| Level1  | Level2 | Level3 | Level4 | Level5 | Level6 | Level7  |
| ------- | ------ | ------ | ------ | ------ | ------ | ------- |
| Kingdom:界 | Phylum:門|Class:綱|Order:目|Family:科|Genus:属| Species :種|

* Color Palette

棒グラフのカラーグラデーションを決める。

* Sort Samples By

棒グラフ整列時の優先順を決める。上記の例では ``sex``（性別）、``biopsy_location``（生検部位）、``diagnosis``（診断名） の順に優先されて整列される。

## 特定の細菌の構成率の分布を可視化する方法
以下のように、特定の細菌（例では Firmicutes 門）に着目して構成率の分布を表す方法は、[こちら](https://qiita.com/keisuke-ota/items/cfd5d0cf835b99d24edd)に紹介している。

![visualization-3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/556342ac-7f2a-580a-0061-590e6e2c6307.png)

# PICRUSt2 analysis

菌叢代謝予測ツール PICRUSt2 は、次のコマンドでプラグインできる。

```terminal:ターミナル
conda install q2-picrust2=2019.7 -c conda-forge -c bioconda -c gavinmdouglas 
```

PICRUSt2 では、以下のデータベースに基づいて酵素群または代謝経路の存在量を予測する。

* EC (Enzyme Commision):酵素群
* KO (KEGG Orthology):酵素群
* MetaCyc Pathway :代謝経路（それぞれ数個〜数十個の酵素群を含む区分）

以下のコマンドで実行する。

```terminal:ターミナル
qiime picrust2 full-pipeline \
--i-table table.qza \
--i-seq rep-seqs.qza \
--output-dir q2-picrust2_output \
--p-threads 4 \
--p-hsp-method pic \
--p-max-nsti 2 \
--verbose
```
EC や KO は区分が細かいため、対象の酵素がピンポイントで定まっている場合に利用するとよい。一方で MetaCyc Pathway は、関連する酵素群をまとめて分類しているため、系全体の現象を解釈するときに便利である。

# データのエクスポート
`` table.qza `` を以下のコマンドで tsv 形式に変換する。

```terminal:ターミナル
qiime tools export \
--input-path table.qza \
--output-path export

biom convert \
-i export/feature-table.biom \
-o export/feature-table.tsv \
--to-tsv
```

同様に、``q2-picrust2_output/pathway_abundance.qza``を tsv 形式に変換する。

```terminal:ターミナル
qiime tools export \
--input-path q2-picrust2_output/pathway_abundance.qza \
--output-path pathway_export

biom convert \
-i pathway_export/feature-table.biom \
-o pathway_export/feature-table.tsv \
--to-tsv
```

これにより、 qiime2 外部の R パッケージである[ALDEx2](https://qiita.com/keisuke-ota/items/b58f9a28a5ecd53b8ecb)などで多変量解析ができる。

# Differential analysis

アンプリコン解析で得られるデータセットは、試料から抽出されたサンプルの read 数であって、試料に含まれる細菌の絶対的な存在量ではない。また、シーケンサーの性能や PCR の精度によって、サンプルの合計 read 数は変わる。したがって read 数はサンプル全体に対して何 % 占めるかという**組成データ**に変換して比較されるべきである。組成データに対する多変量解析のツールとしてANCOM、Balance、[ALDEx2](https://qiita.com/keisuke-ota/items/b58f9a28a5ecd53b8ecb)などがある。

## ANCOM

[ANCOM (Analysis of composition of microbiomes) ](https://www.tandfonline.com/doi/full/10.3402/mehd.v26.27663)は、Qiime2 の[公式チュートリアル](https://docs.qiime2.org/2019.7/tutorials/moving-pictures/)でも紹介されている、組成データに対する多変量解析手法の 1 つである。

### ANCOM モデル

あるデータセット全体に、$m$ 種の特徴量が存在し、群 $X$ に属するサンプル $i$ のカウントデータを $$X_i=(X_{i,1} \cdots \cdots, X_{i,m})$$ とおく。ここで、ある特徴量 $k$ に着目する。そして、$k$ を除く$m-1$ 種の特徴量 $l$ に対して、次の ANOVA モデルを総当たりで当てはめていく。

```math
ln \frac{X_{i,k}+1}{X_{i,l}+1} =a_{k,l}+b_{k,l,X}+e_{k,l}
```

$a_{k,l}$ ：特徴量 $k$ 、 $l$ における全サンプルの平均
$b_{k,l,X}$ ：群 $X$ の効果を表す偏回帰係数
$e_{k,l}$ ：過分散

モデルには、$$X_{i,l}=0$$ の場合でも計算できるように、pseudocount （$+1$） が組み込まれている。

同様に、群 $Y$ に属するサンプル $j$ のカウントデータを 
$$Y_j=(Y_{j,1} \cdots \cdots, Y_{j,m})$$
においても、特徴量 $k$ に対して $m-1$ 通りの ANOVA モデルを当てはめていく。

```math
ln \frac{Y_{j,k}+1}{Y_{j,l}+1} =a_{k,l}+b_{k,l,Y}+e_{k,l}
```

ANCOM では、

```math
b_{k,l,X}=b_{k,l,Y}= 0
```
という帰無仮説について ANOVA を実施する。すなわち各特徴量は、自身以外の m-1 通りの特徴量を相手に検定される。群X、Yの間に定量的な差がある特徴量は、どんな特徴量を相手にしても帰無仮説が棄却されるため、最大で m-1 個の帰無仮説が棄却される。
ANCOM では、ある特徴量 k に対して棄却された帰無仮説の数を $W_k$ とする。この値が大きいほど、群間に定量的な差がある特徴量と判断される。

### CLR:Centered log-ratio transformation について

ANCOM は定量的な差の有無を判定できるが、差の大きさは判定できない。そこで Qiime2 の ANCOM では、CLR という対数変換を用いて、差の大きさを独自で評価している。なお、ANCOM 専用の [python パッケージ](http://scikit-bio.org/docs/0.4.2/generated/generated/skbio.stats.composition.ancom.html)に CLR は実装されていない。

さて m 個の特徴量に対する組成データの標本空間 $S^{m-1} $ は、

$$
S^{m-1}=\left\\{ (x_{1}, x_{2}, \cdots,x_{m} ) \ | \ x_{i} \geq 0 \ (i=1,2, \cdots,m) \sum_{i=1}^{m}x_{i}= 1　\right\\}
$$

である。位相幾何学では、この m-1 次元空間を**単体空間（Simplex space）** という。サンプル i のカウントデータ $$X_i=(X_{i,1} \cdots \cdots, X_{i,m})$$は、
$$
\frac{X_i}{\sum_{k=1}^{m}X_{i,k}}=(x_{1}, x_{2}, \cdots,x_{m})
$$

によって、単体空間 $S^{m-1} $ 上に変換される。この空間において、特徴量1と2の差は、

$$
x_{2}-x_{1}= \frac{X_{2}-X_{1}}{\sum_{i=1}^{m} X_{i}}
$$

となり、無関係の特徴量にも影響される。このような性質は、差の大きさを評価するのに適さないといえる。この問題は、$S^{m-1}$ を $m-1$ 次元のユークリッド空間 $R^{m-1} $ に変換（全単射）する、**有心対数比変換（Centered log-ratio transformation:CLR）**によって解決することが知られている。
組成データ $ x =(x_{1}, x_{2}, \cdots,x_{m} )$ に対する CLR とは

$$
clr(x)=\left\\{log_e \frac{x_{1}}{g(X)}, log_e \frac{x_{2}}{g(X)}, \cdots,log_e \frac{x_{m}}{g(X)} \right\\}  
$$

である。ただし、$g(X)$ を幾何平均を

$$
g(x)=\sqrt[n]{x_{1}x_{2}\cdots x_{m}}
$$

とする。このとき、特徴量1と2の差は

$$
log_e \frac{x_{1}}{g(x)}-log_e \frac{x_{2}}{g(x)} = log_e \frac{x_{1}}{x_{2}}
$$

となり、幾何平均 $g(x)$ に依存しない。このような性質から、CLR は組成データを比較する有効な手段と考えられている。

### コマンド

Qiime2 の解析パイプラインに沿って ANCOM を実施してみる。
今回は次のコマンドで CD 群と対象群のサンプルを抽出し、両軍の菌叢を比較する。

```terminal:ターミナル
qiime feature-table filter-samples   \
--i-table selected_table.qza   \
--m-metadata-file Metadata.tsv   \
--p-where "[diagnosis]='Control' OR [diagnosis]='CD'"   \
--o-filtered-table CD_table.qza
```
上述の Taxonomy analysis で ``taxonomy.qza`` を求めている場合、sub-OTU で集計されている ``table.qza`` を genus （属）や family （科）などの上位階層で集計し直すこともできる。 
次のコマンドでは、sub-OTU のカウントデータを genus レベルに合算している。

```terminal:ターミナル
qiime taxa collapse   \
--i-table All/CD_table.qza   \
--i-taxonomy taxonomy.qza   \
--p-level 6   \
--o-collapsed-table CD_L6_table.qza
```
ANCOM モデルにおける pseudocount を実施する。すなわち全てののカウントデータに対して $+1$ される。

```terminal:ターミナル
qiime composition add-pseudocount   \
--i-table CD_L6_table.qza   \
--o-composition-table CD_L6_comp_table.qza
```
次のコマンドで、ANCOM を実施する。

```terminal:ターミナル
qiime composition ancom   \
--i-table CD_L6_comp_table.qza   \
--m-metadata-file Metadata.tsv   \
--m-metadata-column diagnosis   \
--o-visualization CD_L6_ancom.qzv
```
### 出力結果の解釈

ANCOM の結果は、**Volcano plot** という形式で出力される。縦軸は、上で説明したとおり、棄却された帰無仮説の個数を表す。横軸は、対象群における CLR 値の平均値から、CD 群におけるそれを差し引いた値である。群間の差が小さい場合、横軸は 0 に近づく。

<img width="643" alt="スクリーンショット 2020-05-07 20.44.40.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/f0ca76a7-82d5-9900-d91a-e74b1c6444ee.png">

ちなみに 3 群以上の比較も実施できる。その場合、横軸は、ANOVA における CLR 値の F 値をとる。群間に差が存在しない場合、F 値は 0 に近づく。非負の値をとることに注意しよう。

各特徴量の詳細な W は、tsv 形式で表示される。このとき、全サンプルの W に対する累積分布関数から閾値が設定され、W が閾値以上の特徴量には統計的有意差があるとみなされ ``True`` と表示される。 

<img width="1028" alt="スクリーンショット 2020-05-07 21.04.01.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/3e5de25a-d164-d97b-9cf1-6f4f42c87d7c.png">

``True`` と判定された特徴量に対しては、ブラウザ上にカウントデータの四分位数が表示される。Pseudocounts の影響で、最低値は必ず 1.0 となる。

<img width="1437" alt="スクリーンショット 2020-05-07 21.08.11.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/e0510fd9-6339-d119-d7db-189d77ef26e4.png">

この例では、該当の特徴量の中央値に着目すると、CD 群では 1.0 だが、 Control 群では二桁である。いずれの genus も、CD 群で低下傾向にあることがわかる。

### 運用上の注意点

ANCOM は、多くの特徴量で群間に差がない場合に、群間に差がある一部の特徴量を索敵するのに適している。もし多くの特徴量で群間に差があると、帰無仮説

```math
b_{k,l,X}=b_{k,l,Y}= 0
```
が棄却される k、l の組み合わせが無くなり、$W_k$ に差が付かなくなるためである。[公式チュートリアル](https://docs.qiime2.org/2019.7/tutorials/moving-pictures/)にも、 「ANCOM は25% 未満の特徴量に差があるデータセットを想定している」と書かれている。

## ALDEx2

Qiime2 外部の R パッケージで解析する。背景やコマンドは[こちら](https://qiita.com/keisuke-ota/items/b58f9a28a5ecd53b8ecb)で紹介している。

# Diversity analysis

Diversity analysis では、サンプルの多様性を表す指標([α-多様性とβ-多様性](https://qiita.com/keisuke-ota/items/9701c583e153df467c05))をもとに、群間を比較する。

## 系統樹の作成

多様性指標の中でも、α 多様性指標である Faith’s Phylogenetic Diversity (Faith PD) と、β 多様性指標である Unifrac 距離は、微生物の系統樹に基づいて計算される。
まずは次のコマンドで、塩基配列データ``rep-seqs.qza``から系統樹``rooted-tree.qza``を作成しておく。上述以外の指標は、系統樹を作成せずとも計算できる。

```terminal:ターミナル
qiime phylogeny align-to-tree-mafft-fasttree   \
--i-sequences rep-seqs.qza   \
--o-alignment aligned-rep-seqs.qza   \
--o-masked-alignment masked-aligned-rep-seqs.qza   \
--o-tree unrooted-tree.qza   \
--o-rooted-tree rooted-tree.qza
```

## 多様性指数の計算

以下の多様性指数は、次のコマンドによって ``core-metrics-results`` にまとめて出力できる。

* α 多様性指数

 -  Shannon index
 -  観測された OTU の数 
 -  Faith’s Phylogenetic Diversity (Faith PD)
 -  Pielou’s Evenness （Pielou の均衡度指数）

* β 多様性指数

 - Jaccard 距離
 - Bray Curtis 距離
 - Unweighted Unifrac 距離
 - Weighted Unifrac 距離


```terminal
qiime diversity core-metrics-phylogenetic   \
--i-phylogeny rooted-tree.qza   \
--i-table selected_table.qza   \
--p-sampling-depth 1500   \
--m-metadata-file metadata.tsv   \
--output-dir core-metrics-results
```
### サンプル深度について
多様性指標は、サンプルの全 read 数に影響されやすい。例えば、1000 read のサンプルと 10000 read のサンプルは正しく比較できないことがある。そこで、サンプルからのランダム抽出で read 数を希薄化 (Rarefaction) することがある。希薄時の抽出 read 数を**サンプル深度**といい、上記コマンドでは``--p-sampling-depth `` で 1500 に設定している。read 数がサンプル深度以下のサンプルは解析から外される。サンプル深度を小さくすると希少種が抽出されないことがあるが、サンプル深度を大きくすると解析対象のサンプル数が減ってしまう。解析から、どのような主張を導きたいかによってサンプル深度は決定する必要がある。

### α 希薄化 (alpha rarefaction)

適切なサンプル深度を決定するために、様々なサンプル深度で α 多様性指数を計算する操作を**α 希薄化(alpha rarefaction)**という。
次のコマンドでは、最大深度 10000 read を 21 段階に区切り、500 read ごとに希薄化し α 多様性指数を求めている。

```terminal:ターミナル
qiime diversity alpha-rarefaction   \
--i-table table.qza   \
--i-phylogeny rooted-tree.qza   \
--p-max-depth 10000   \
--p-steps 21   \
--m-metadata-file Metadata.tsv \
--o-visualization alpha-rarefaction.qzv
```

Shannon index は、[こちら](https://qiita.com/keisuke-ota/items/9701c583e153df467c05)で説明したとおり主要な種に影響を受けやすいため、サンプル深度が小さい段階から閾値に近い。下図では、CD 群、Control 群、UC 群の腸内細菌叢における shannon index を表しているが、1500 read 以上でほぼ一定である。このことから、1500 read までの抽出で腸内細菌の主要な種は含まれていると考えられる。

<img width="1066" alt="スクリーンショット 2020-05-10 14.17.17.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/76a2879e-1bee-2647-fa6e-0435bc7329fa.png">

一方、下図の observed OTUs （一度でも観測された OTU の数）は、サンプル深度と共に増え続けている。これは、サンプル深度を大きくするほど希少種が抽出されるためである。この傾向は特に菌叢多様性の高い土壌サンプルなどで見られる。多様な種を抽出することを優先してサンプル深度を大きくすると、解析対象のサンプル数が減ってしまうので注意すること。

<img width="1076" alt="スクリーンショット 2020-05-10 14.25.07.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/ef092f77-d153-6847-449a-d4256ba22590.png">

### その他について

なお、[こちら](https://qiita.com/keisuke-ota/items/9701c583e153df467c05)で紹介しているように、多様性指数には他の種類もあるため、必要に応じて個別に計算してみよう。


*  Chao1

```terminal:ターミナル
qiime diversity alpha   \
--i-table selected_table.qza  \
--p-metric 'chao1' \
--o-alpha-diversity core-metrics-results/chao1_vector.qza
```

* Simpson

```terminal:ターミナル
qiime diversity alpha   \
--i-table selected_table.qza  \
--p-metric 'simpson' \
--o-alpha-diversity core-metrics-results/simpson_vector.qza
```

* Simpson Index of Evenness

```terminal:ターミナル
qiime diversity alpha   \
--i-table selected_table.qza  \
--p-metric 'simpson_e' \
--o-alpha-diversity core-metrics-results/simpson_e_vector.qza
```
## α 多様性指数の比較

* shannon index の比較

```terminal:ターミナル
qiime diversity alpha-group-significance   \
--i-alpha-diversity core-metrics-results/shannon_vector.qza   \
--m-metadata-file Metadata.tsv   \
--o-visualization core-metrics-results/shannon_significance.qzv
```

Shannon index を縦軸に、CD 群、Control 群、 UC 群の箱ひげ図が表示される。

<img width="935" alt="スクリーンショット 2020-05-10 14.50.50.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/eff8b5da-1ef9-53f2-17d3-b7611ad9a1b8.png">

さらに、[クラスカル＝ウォリス検定](https://bellcurve.jp/statistics/glossary/1252.html)の結果が表示される。今回は、CD 群、Control 群、UC 群の 3 群比較であり、pairwise テストは 3 回実施されている。したがって、多重比較による[第一種の過誤](https://bellcurve.jp/statistics/glossary/1766.html)を防止するために、[BH法](https://products.sint.co.jp/aisia/blog/vol1-11)で偽陽性率（FDR）を調整した ``q-value`` を参考にするのが良い。

<img width="654" alt="スクリーンショット 2020-05-10 14.53.30.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/be10cc52-f22c-bbe0-fa4a-b157b796b1d3.png">


* Observed OTUs の検定

```terminal:ターミナル
qiime diversity alpha-group-significance   \
--i-alpha-diversity core-metrics-results/observed_otus_vector.qza   \
--m-metadata-file Metadata.tsv   \
--o-visualization core-metrics-results/observed_otus_significance.qzv
```

* Faith’s Phylogenetic Diversity (Faith PD) の検定

```terminal:ターミナル
qiime diversity alpha-group-significance   \
--i-alpha-diversity core-metrics-results/faith_pd_vector.qza   \
--m-metadata-file Metadata.tsv   \
--o-visualization core-metrics-results/faith_pd_significance.qzv
```

* Pielou’s Evenness （Pielou の均衡度指数） の検定

```terminal:ターミナル
qiime diversity alpha-group-significance   \
--i-alpha-diversity core-metrics-results/evenness_vector.qza   \
--m-metadata-file Metadata.tsv   \
--o-visualization core-metrics-results/evenness_significance.qzv
```

## β 多様性指数の比較

β 多様性の有意差検定には、ノンパラメトリック分散分析として知られる PERMANOVA(PERmutational Multivariate ANalysis Of VAriance) を採用する。

### PERMANOVA について

PERMANOVAでは、通常の分散分析と同様に、平均に対する二乗和（Sum of Squares:SS） から F 統計量を求める。

$$
F= \frac{ \frac{SS_T-SS_w}{a-1}}{\frac{SS_w}{N-a}} = \frac{ \frac{SS_a}{a-1}}{\frac{SS_w}{N-a}}
$$

$SS_T$ ：全てのサンプル間の β 多様性指標の平均に対する二乗和
$SS_W$ ：同一群内のサンプル間の β 多様性指標の平均に対する二乗和
$SS_A$ ：異なる群に属するサンプル間の β 多様性指標の平均に対する二乗和
$N$ ：サンプルの数
$a$ ：群の数

群間の差が小さい場合、すなわち群間で菌叢に違いがない場合、$SS_A$ が小さくなり、F も小さくなる。一方で、菌叢の差が大きい群間ほど F は大きくなる。

ただし通常の F 検定とは異なり、F 分布を帰無分布と仮定していない。帰無分布には、サンプルをランダムに並べ替え （permutation） て計算された $F^H$ の分布が採用されている。Qiime2 ではデフォルトで 999 回の並べ替えシミュレーションを実施しており、その際の $F^H$ の分布をもとに検定している。

### コマンド

検定コマンドでは、メタデータのどの項目で菌叢を比較するかを設定する ``--m-metadata-column`` が必須である。

* Jaccard 距離の比較

```terminal:ターミナル
qiime diversity beta-group-significance   \
--i-distance-matrix core-metrics-results/jaccard_distance_matrix.qza   \
--m-metadata-file Metadata.tsv   \
--m-metadata-column diagnosis   \
--o-visualization core-metrics-results/jaccard_distance_significance.qzv   \
--p-pairwise
```

* Bray Curtis 距離の比較

```terminal:ターミナル
qiime diversity beta-group-significance   \
--i-distance-matrix core-metrics-results/bray_curtis_distance_matrix.qza   \
--m-metadata-file Metadata.tsv   \
--m-metadata-column diagnosis   \
--o-visualization core-metrics-results/bray_curtis_distance_significance.qzv   \
--p-pairwise
```

* Unweighted Unifrac 距離の比較

```terminal:ターミナル
qiime diversity beta-group-significance   \
--i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza   \
--m-metadata-file Metadata.tsv   \
--m-metadata-column diagnosis   \
--o-visualization core-metrics-results/unweighted_unifrac_significance.qzv   \
--p-pairwise
```

* Weighted Unifrac 距離の比較

```terminal:ターミナル
qiime diversity beta-group-significance   \
--i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza   \
--m-metadata-file Metadata.tsv   \
--m-metadata-column diagnosis   \
--o-visualization core-metrics-results/weighted_unifrac_significance.qzv   \
--p-pairwise
```

β 多様性指標の成果物は下図のとおり。

<img width="947" alt="スクリーンショット 2020-05-10 15.24.18.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/86488326-8dc9-ec1c-e975-0243092ee074.png">

Group significance plots では、「同一群内のサンプル間の β 多様性指標の箱ひげ図」が最も左に表示され、他は「異なる群に属するサンプル間の β 多様性指標の箱ひげ図」となる。
例では、CD 群（n=75）、Control 群（n=45）、UC 群（n=34） を比較している。CD 群に着目すると、「同一群内のサンプル間の β 多様性指標」は、75×74=2775 通り計算され、「Control 群に属するサンプル間の β 多様性指標」は、75×45=3375 通り計算され、「UC 群に属するサンプル間の β 多様性指標」は、75×34=2550 通り計算されている。

PERMANOVAでは、CD 群、Control 群、UC 群の 3 群を比較しており、pairwise テストは 3 回実施されている。したがって、多重比較による[第一種の過誤](https://bellcurve.jp/statistics/glossary/1766.html)を防止するために、[BH法](https://products.sint.co.jp/aisia/blog/vol1-11)で偽陽性率（FDR）を調整した ``q-value`` を参考にするのが良い。