---
title: "[第24回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_24)のウェブ資料"
date: '最終更新: `r format(Sys.Date(), "%Y/%m/%d")`'
output: 
  html_document:
    toc: TRUE
    toc_float: TRUE
    code_download: TRUE
---
```{r klippy, echo=FALSE, include=TRUE}
klippy::klippy()
```

# はじめに
  - [乳酸菌NGS連載第22回のPDF](https://www.iu.a.u-tokyo.ac.jp/kadota/JSLAB_22_kadota.pdf)
  - [dada2](https://bioconductor.org/packages/release/bioc/html/dada2.html)：[Callahan et al., Nat Methods, 2016](https://pubmed.ncbi.nlm.nih.gov/27214047/)   
- [SRP125723](https://www.ncbi.nlm.nih.gov/sra/?term=SRP125723) ([Hallmaier-Wacker et al., Sci Rep., 2018](https://pubmed.ncbi.nlm.nih.gov/29674641/))  
今回用いるリードデータです。オリジナルは計64サンプル分のデータを取得していますが、以下の12サンプル分のデータのみ解析します。  

| サンプルID | DNA抽出キット | 保存バッファ | リード数 |
|:---:|:---:|:---:|:---:|
| [SRR6325813](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325813) | QMINI | Native | 282517 |
| [SRR6325815](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325815) | QMINI | LysisBuffer | 266259 |
| [SRR6325816](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325816) | QMINI | RNA-later | 359984 |
| [SRR6325817](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325817) | QMINI | PBS | 325660 |
| [SRR6325828](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325828) | MOBIO | RNA-later | 206939 |
| [SRR6325829](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325829) | MOBIO | PBS | 261767 |
| [SRR6325830](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325830) | MOBIO | Native | 238627 |
| [SRR6325831](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325831) | MOBIO | LysisBuffer | 225548 |
| [SRR6325859](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325859) | GENIAL | LysisBuffer | 290383 |
| [SRR6325860](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325860) | GENIAL | Native | 253908 |
| [SRR6325861](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325861) | GENIAL | PBS | 120789 |
| [SRR6325862](https://www.ncbi.nlm.nih.gov/sra/?term=SRR6325862) | GENIAL | RNA-later | 258792 |

# W1 クオリティコントロール
## W1.1 事前準備
　第23回のW6.2と基本的に同じスクリプトです。`fnFs`と`fnRs`がQC前のファイルの相対パス情報、`filtFs`と`filtRs`がQC後のファイルの相対パス情報になります。
```{r, message=FALSE, error=TRUE, warning=TRUE}
basho <- "."                              # 相対パスなので、現在の作業ディレクトリを"."で表現
# Forward側
fnFs <- list.files(pattern = "_1.fastq", path = basho, full.names = TRUE)
hoge <- strsplit(x = basename(fnFs), split = "_")# basename(fnFs)の文字列を"_"で分割
sampleID <- sapply(X = hoge, FUN = `[`, 1)# hogeの1番目の要素のみを抽出
filtFs <- file.path(basho, "processed", paste0(sampleID, "_F_filt.fastq"))# ファイルパス情報を作成
names(filtFs) <- sampleID                 # sampleIDをfiltFsベクトルの要素名にする
# Reverse側
fnRs <- list.files(pattern = "_2.fastq", path = basho, full.names = TRUE)
hoge <- strsplit(x = basename(fnRs), split = "_")# basename(fnRs)の文字列を"_"で分割
sampleID <- sapply(X = hoge, FUN = `[`, 1)# hogeの1番目の要素のみを抽出
filtRs <- file.path(basho, "processed", paste0(sampleID, "_R_filt.fastq"))# ファイルパス情報を作成
names(filtRs) <- sampleID                 # sampleIDをfiltRsベクトルの要素名にする
```

```{r, include=FALSE}
library(dada2)                            # dada2パッケージのロード
library(Biostrings)                       # Biostringsパッケージのロード
```

- [乳酸菌NGS連載第4回のPDF](https://www.iu.a.u-tokyo.ac.jp/kadota/JSLAB_4_kadota.pdf)
  
## W1.2 `plotQualityProfile`関数
　`plotQualityProfile`関数は、fastqファイルを入力として、リード中の塩基ポジションごとのクオリティスコアのプロットを出力します（約10分）。  

### W1.2.1 QC前のForward側
　まずは`fnFs`オブジェクトの中身を表示させています。これは、QC前のForward側のfastqファイルの相対パス情報になります。`plotQualityProfile`関数の入力として、`fnFs`オブジェクトがそのまま利用されていることがわかります。入出力の関係からも明らかですが、入力は文字列ベクトルですので、ベクトルの要素数がファイル数に相当します。これが図1と同じものになります。
```{r, message=FALSE, error=TRUE, warning=FALSE, fig.dim = c(8, 7)}
fnFs                                      # fnFsの中身を確認
plotQualityProfile(fnFs)                  # fnFsのクオリティスコアのプロット
```

　入力が12個のファイル名情報ですので、出力として12個のプロットが得られます。各プロットの横軸はリード中の塩基のポジション（5'端が左、3'端が右）、縦軸はクオリティスコアです。[連載第4回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_4)でも解説していますが、高ければ高いほどよいと解釈します。

### W1.2.2 QC前のReverse側
　次に、`fnRs`オブジェクトの中身を表示させています。これは、QC前のReverse側のfastqファイルの相対パス情報になります。`plotQualityProfile`関数の入力として、`fnRs`オブジェクトをそのまま与えています。
```{r, message=FALSE, error=TRUE, warning=FALSE, fig.dim = c(8, 7)}
fnRs                                      # fnRsの中身を確認
plotQualityProfile(fnRs)                  # fnRsクオリティスコアのプロット
```

## W1.3 `filterAndTrim`関数
　`filterAndTrim`関数は、リードのfastqファイルを入力として、クオリティの低い領域のトリム、Nを含むリードのフィルタリング、サンプル中に含まれるPhiXというポジティブコントロールサンプル由来リードのフィルタリングなどを実行します（約5分）。主な出力結果であるQC後のリードファイルは`filtFs`と`filtRs`ですが、QC前後のリード数の情報は、`out`オブジェクトに保存されます。なお、実行中に「Multithreading has been DISABLED, as forking is not supported on .Platform$OS.type 'windows'」のようなメッセージが出ても、並列計算ができないだけで計算自体は動いているので問題ありません。
```{r, message=FALSE, error=TRUE, warning=FALSE}
out <- filterAndTrim(                     # filterandTrim関数を実行し、結果をoutオブジェクトに格納
  fwd = fnFs,                             # QC前のForward(F)側のファイル名
  filt = filtFs,                          # QC後のForward(F)側のファイル名
  rev = fnRs,                             # QC前のReverse(R)側のファイル名
  filt.rev =filtRs,                       # QC後のReverse(R)側のファイル名
  maxN = c(0, 0),                         # リードに含まれるNの許容数はF側R側ともに0
  maxEE = c(2, 5),                        # リード中のエラー率の和の閾値はF側が2、R側が5
  truncLen = c(225, 160),                 # F側は225、R側は160 bpになるまで右側をトリム
  trimLeft = 20,                          # リードの左側を20 bpトリム
  rm.phix = TRUE,                         # PhiX配列由来のリードを除去
  compress = FALSE,                       # 結果を圧縮ファイルにしない
  multithread = TRUE                      # マルチスレッドで計算を高速化する
)
```

この関数で与えている情報を以下にまとめます：

- `fwd = fnFs`  
QC前のファイルの相対パス情報（Forward側）を指定する部分です。オプション名の`fwd`は、forwardの略だと解釈します。  
- `filt = filtFs`  
QC後のファイルの相対パス情報（Forward側）を指定する部分です。オプション名の`filt`は、フィルタリング（filtering）後のファイル名情報を指定せよという意味だと解釈します。  
- `rev = fnRs`  
QC前のファイルの相対パス情報（Reverse側）を指定する部分です。オプション名の`rev``は、reverseの略だと解釈します。  
- `filt.rev = filtRs`  
QC後のファイルの相対パス情報（Reverse側）を指定する部分です。オプション名の`filt.rev`は、フィルタリング（filtering）後のreverse側ファイル名情報を指定せよという意味だと解釈します。  
- `maxN = c(0, 0)`  
リードに含まれるNの許容数を指定する部分です。`c(0, 0)`は、0という要素が2つある数値ベクトルを意味し、F側R側ともにNを含むリードを1つも許容しないという意味です。要素が2つになっているのは、Forward側とReverse側の2種類のファイルを入力として与えているためです。例えば、もし`c(1, 3)`なら、Forward側は、塩基配列中にNが全く含まれないリードだけでなく、1個含むリードまで許容する（2個以上Nを含むリードを捨てる）ことを意味します。また、Reverse側は、0～3個Nを含むリードまで許容することを意味します。  
- `maxEE = c(2, 5)`  
リード中のエラー率の和の閾値を指定する部分です。`c(2, 5)`は、2つの要素からなる数値ベクトル（1番目の要素が2、2番目の要素が5）を意味し、F側が2、R側が5という意味です。この場合はリード長が250 bpありますが、塩基ごとにクオリティスコア（つまりエラー率）の情報が付与されています。それを全部足していったエラー率の合計が2未満のリードのみ残す（それ以上のリードは捨てる）のがForward側、5未満のリードのみ残す（それ以上のリードは捨てる）のがReverse側だと解釈すればよいです（参考論文は[Edgar and Flyvbjerg, 2015](https://pubmed.ncbi.nlm.nih.gov/26139637/)）。閾値が2というのは、2/250 = 0.008ということになるので、塩基当たりの平均で0.8%のエラー率まで許容するのがForward側、5/250 = 2%のエラー率まで許容するのがReverse側なのだと理解します。Reverse側の値が大きいのは、こちらのほうが全体的にクオリティスコアが低いためです。  
- `truncLen = c(225, 160)`  
Forward側は225 bp、Reverse側は160 bpになるまで右側（3'末端）をトリム  せよという指令です。このオプションは、Forward側とReverse側でトリムしたい長さが異なる場合に利用します。補足として、もしForward側Reverse側関係なく223 bpになるまで右側（3'末端）をトリムしたいだけであれば、このオプションは丸々削除し、以下の`trimLeft`と同じノリで、`trimRight = 27`とすればよいです（リード長は250 bpなので、250 - 223 = 27）。なお、 `truncLen = c(225, 160)`の場合はReverse側は250 - 160 = 90 bpだけ右側をトリムさせることに相当します。基本的に`trimRight`と共存させないことをおススメします。もし`trimRight`を90 bp以下の数値を指定すると（例：`trimRight = 27`）それは何の影響も与えない一方で、90より大きい数値を指定すると全リードが除去されてしまいますのでご注意ください（なぜそうなるかという理由云々というよりはそういう仕様だということです）。
- `trimLeft = 25`  
これはそのまんまなのですが、リードの左側（5'端）を25 bp分トリムせよという指令です。これは、Forward側Reverse側関係なくクオリティスコア分布のx軸の左側の25 bp分（つまり1～25番目の塩基）は、50～150塩基付近に比べて相対的にスコアが低く見えるからです。別に25じゃなくて20でも30とかでもよいのですが、なんとなく25 bpかなという感じです。  
- `rm.phix = TRUE`  
PhiXというIlluminaのコントロールライブラリ由来のリード配列を除去せよという指令です。PhiXは配列既知のバクテリオファージ（スモールゲノム）です。サンプル由来リードではありませんし同定が簡単ですので、このオプションは無条件でTRUEが基本です。
- `compress = FALSE`  
出力結果を圧縮ファイルにはしないという指令です。これはTRUEにしても特に問題はないと思います。
- `multithread = TRUE`  
複数のCPUを使って計算せよという指令です。他の計算も同時進行で走らせていて、そちらが遅くなるのはいやな場合はFALSEでもよいです。特にCPU数を把握していないユーザは、とりあえずTRUEのままでやってみて差し支えありません。

## W1.4 QC後のスコア分布を確認
　`plotQualityProfile`関数は、fastqファイルを入力として、リード中の塩基ポジションごとのクオリティスコアのプロットを出力します（約10分）。  

### W1.4.1 QC後のForward側
　まず、`filtFs`オブジェクトの中身を表示させています。これは、QC後のForward側のfastqファイルの相対パス情報になります。`plotQualityProfile`関数の入力として、`filtFs`オブジェクトをそのまま与えています。
```{r, message=FALSE, error=TRUE, warning=FALSE, fig.dim = c(8, 7)}
filtFs                                    # filtFsの中身を確認
plotQualityProfile(filtFs)                # filtFsクオリティスコアのプロット
```

### W1.4.2 QC後のReverse側
　次に、まず`filtRs`オブジェクトの中身を表示させています。これは、QC後のReverse側のfastqファイルの相対パス情報になります。`plotQualityProfile`関数の入力として、`filtRs`オブジェクトをそのまま与えています。
```{r, message=FALSE, error=TRUE, warning=FALSE, fig.dim = c(8, 7)}
filtRs                                    # filtRsの中身を確認
plotQualityProfile(filtRs)                # filtRsのクオリティスコアのプロット
```

## W1.5 QC前後のリード数を確認
　W1.3で得た`filterAndTrim`関数の実行結果が格納された`out`オブジェクトを用いて、QC前後でリード数がどのように変わったかを確認します。ここで得たものが、表1のリード数情報になります。

### W1.5.1 基本形
　`out`の中身を表示させているだけです。reads.in列がQC前、reads.out列がQC後のリード数です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
out                                       # outの中身を確認
```

### W1.5.2 割合の情報を追加
　先ほどの結果より、`out`の行列は1列目がQC前、2列目がQC後のリード数であることはすぐにわかります。ここではまず、2列目/1列目として残ったリードの割合を`ratio`というオブジェクトにまず格納しています。次に、元々2列の情報からなる`out`行列の右側に、先ほど得た`ratio`オブジェクトを`cbind`関数を用いて列方向で結合して出力しています。ここで得たものが、表1bの割合の情報になります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
ratio <- out[,2]/out[,1]                  # outの2列目/1列目を計算し、ratioに格納
cbind(out, ratio)                         # outの右側にratioの情報を列方向で連結
```
　トリムやフィルタリングで残ったリードの割合の最大値と最小値は、それぞれ以下の通りです：

- 最大値とそのサンプル名：`r max(ratio)`、`r rownames(out)[which.max(ratio)]`  
- 最小値とそのサンプル名：`r min(ratio)`、`r rownames(out)[which.min(ratio)]`  

## W1.6 ユニーク配列の数
 ユニーク配列（unique sequences）とは、長さはもちろんのこと、同一ポジションの塩基が完全に一致する配列のことです。ここでは、[dada2](https://bioconductor.org/packages/release/bioc/html/dada2.html)パッケージが提供する`getUniques`関数を用いてユニーク配列の数をカウントする例を示します。ここで得たものが、表1のユニーク配列の情報になります。
 
### W1.6.1 基本的な考え方
　いきなり250 bpの長いリードで結果を示しても理解しづらいので、ここではまず第23回の図2bの左下に見えている7本のリード配列をサンプルデータとして用います。そして、実質的に`getUniques`関数と同じ挙動を示す`table`関数を用いて感覚をつかみます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- c("ACCGAGGAGGTCGGCGCCTC",         # 第23回の図2bの左下の7配列をhogeに格納
          "ACGGAGGAGGTCGGCGCCTC",         # 第23回の図2bの左下の7配列をhogeに格納
          "ACGGAGGACGTCGGCGCCTA",         # 第23回の図2bの左下の7配列をhogeに格納
          "ACGGAGGAGGTCGGCGCCTC",         # 第23回の図2bの左下の7配列をhogeに格納
          "ACGGAGGAGGTCGGCGCCTC",         # 第23回の図2bの左下の7配列をhogeに格納
          "ATGGAGGAGGTCGGCGCCTC",         # 第23回の図2bの左下の7配列をhogeに格納
          "ACGGATGAGGTCGGCGACTC")         # 第23回の図2bの左下の7配列をhogeに格納
table(hoge)                               # hogeベクトル中の要素の種類ごとの出現頻度を確認
```
　`hoge`は文字列ベクトルであり、計`r length(hoge)`つのリード塩基配列が格納されていることがわかります。つまり、リード数は`r length(hoge)`です。そして、よく眺めると、2, 4, 5番目のリードのみが同一であり、それ以外はどこかのポジションの塩基が異なっていることがわかります。つまり、ユニーク配列の数は`r length(table(hoge))`です。この事実を認識したうえで`table(hoge)`実行結果を眺めると、`table`関数の挙動がよくわかると思います。つまり、`table`関数にベクトルを入力として与えると、ユニークな要素を同定するとともに、要素ごとの出現頻度を返してくれるのです。

　この場合の`table(hoge)`実行結果の要素数がユニークな配列数となりますので、ベクトルの要素数を調べる`length`関数を引き続いて実行す
れば、「ユニークな配列数は`r length(table(hoge))`個」という結果が得られるのです。
```{r, message=FALSE, error=TRUE, warning=FALSE}
length(table(hoge))                       # ベクトル中のユニークな要素数を返す
```

　同様に、`table(hoge)`実行結果として得られる出現頻度情報を`sum`関数を用いて足すと、元のリード数情報（つまり`r sum(table(hoge))`個）が得られます。

```{r, message=FALSE, error=TRUE, warning=FALSE}
sum(table(hoge))                          # 出現頻度の総和（つまりリード数）を計算
```

　なお、`table(hoge)`実行結果は、デフォルトでは要素（この場合は塩基配列）をアルファベット順にソートした状態で返されます。出現頻度を降順（多 $\rightarrow$ 少）にした状態にするには、以下のように`sort`関数を`decreasing = T`というオプションつきで引き続いて実行すればよいです。`decreasing = T`はソート結果を降順（多 $\rightarrow$ 少）にせよというオプションです。
　
```{r, message=FALSE, error=TRUE, warning=FALSE}
sort(table(hoge), decreasing = T)         # ベクトル中の要素ごとの出現頻度を降順で返す
```

　このあと解説する`getUniques`関数は、上記スクリプトのような塩基配列ごとの出現頻度をデフォルトで降順（多 $\rightarrow$ 少）の状態で返してくれます。

### W1.6.2 `getUniques`関数
　次に、`getUniques`関数の利用例を示します。この関数は、fastqファイルを読み込んで、ユニーク配列ごとの出現頻度を降順（多 $\rightarrow$ 少）で返してくれます。ここではまず、「QC前のForward側のfastqファイル群の名前」情報を格納した`fnFs`オブジェクトの1番目の要素（つまり"./SRR6325813_1.fastq"）のみを入力として与えて、`getUniques`関数実行結果として得られるユニーク配列ごとの出現頻度情報に、さらに`sum`関数を実行しています。つまり、`getUniques`関数を利用してリード数情報を得ようとしています。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- fnFs[1]                           # fnFsの1番目の要素をhogeに格納
hoge                                      # 中身を確認
sum(getUniques(hoge))                     # リード数情報
```
　得られた結果（つまり`r sum(getUniques(hoge))`個）は、確かに第23回の表2で示されているリード数と一致しています。次に、出現頻度の多い上位3つを`head`関数を用いて示します。1つ1つのリードが250 bp（250文字）あるので見づらいかと思いますが、1番多いのは`r getUniques(hoge)[1]`個、2番目に多いのは`r getUniques(hoge)[2]`個と解釈します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
head(getUniques(hoge), n = 3)             # hoge中のユニーク配列の最初の3つ（塩基配列と出現頻度）を確認
```

　次に、出現頻度の少ない上位3つを`tail`関数を用いて示します。予想通りですが、3つとも1個だと解釈します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
tail(getUniques(hoge), n = 3)             # hoge中のユニーク配列の最後の3つ（塩基配列と出現頻度）を確認
```

　**ユニーク配列の数は、以下のように`length`関数を用いて得ることができます。**
```{r, message=FALSE, error=TRUE, warning=FALSE}
length(getUniques(hoge))                  # ユニークな配列数
```

　このサンプル（`r hoge`）のリード数は`r sum(getUniques(hoge))`個ですので、ユニーク配列1本あたり`r sum(getUniques(hoge))` \/ `r length(getUniques(hoge))` = `r sum(getUniques(hoge))/length(getUniques(hoge))`個の同一塩基配列のリードが存在すると解釈できます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
sum(getUniques(hoge)) / length(getUniques(hoge)) # ユニークな配列1本あたりの出現頻度を計算
```

### W1.6.3 自作の`getUniNum`関数
　様々な書き方（定義の仕方）がありますが、例えば以下のようにすれば、ファイル中のユニーク配列の数を得る`getUniNum`という名前の関数を定義することができます。この関数の実体は、`getUniques`関数実行に引き続いて`length`関数を実行することです。なお、`x`は入力として与える情報に相当します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
getUniNum <- function(x) length(getUniques(x))   # ユニーク配列数を得るための関数を定義
```

　一度定義してしまえば、あとは通常の関数と同様の利用法ができます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- fnFs[1]                           # fastqファイル名をhogeに格納
getUniNum(hoge)                           # ユニーク配列数を確認
```

　`getUniNum`関数は、1つのfastqファイルの入力しか想定していません。そのため、下記のように「QC前のForward側のfastqファイル群の名前」情報を格納した`fnFs`オブジェクトをそのまま入力として与えるとエラーが出ます。エラーが出ることが分かっている状態でエラーメッセージを眺めると、確かに「...Unrecognized format...」などと書かれている意味がなんとなく分かると思います。
```{r, message=FALSE, error=TRUE, warning=FALSE}
fnFs                                      # 中身を確認
getUniNum(fnFs)                           # ユニーク配列数を確認
```

### W1.6.4 `lapply`関数を利用
　12個のfastqファイルそれぞれのリード数を一度に実行する1つのやり方は、`lapply`関数の利用です。これは、`fnFs`の要素ごとに`getUniNum`関数を実行してくれます。実行結果は[[1]]のような**リスト形式**で返されますが、表示されているユニーク配列数は正しいです。
```{r, message=FALSE, error=TRUE, warning=FALSE}
lapply(fnFs, getUniNum)                   # fnFs中のサンプルごとのユニーク配列数をリスト形式で確認
```

### W1.6.5 `unlist`関数も併用
　結果をリスト形式ではなくベクトル形式で得たい場合は、以下で示すように`unlist`関数を追加で実行すればよいです。
```{r, message=FALSE, error=TRUE, warning=FALSE}
unlist(lapply(fnFs, getUniNum))           # fnFs中のサンプルごとのユニーク配列数をベクトル形式で確認
```

### W1.6.6 `sapply`関数を利用
　以下で示すように`sapply`関数を利用してもかまいません。このあたりは好みですので、**多様なやり方を知り、見慣れない関数や出力結果に狼狽しないことのみが重要**です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
sapply(fnFs, getUniNum)                   # fnFs中のサンプルごとのユニーク配列数を確認
```

### W1.6.7 QC前のForward側
　実質的に直前の項と同じ内容ですが、`sapply`実行結果を`uniseq_fnFs`に格納しています。これが「表1、(a)QC前、ユニーク配列の数、Forward」の基礎情報です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
uniseq_fnFs <- sapply(fnFs, getUniNum)    # fnFs中のサンプルごとのユニーク配列数をuniseq_fnFsに格納
uniseq_fnFs                               # 中身を確認
```

### W1.6.8 QC前のReverse側
　「表1、(a)QC前、ユニーク配列の数、Reverse」の基礎情報です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
uniseq_fnRs <- sapply(fnRs, getUniNum)    # fnRs中のサンプルごとのユニーク配列数をuniseq_fnRsに格納
uniseq_fnRs                               # 中身を確認
```

### W1.6.9 QC後のForward側
　「表1、(b)QC後、ユニーク配列の数、Forward」の基礎情報です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
uniseq_filtFs <- sapply(filtFs, getUniNum)# filtFs中のサンプルごとのユニーク配列数をuniseq_filtFsに格納
uniseq_filtFs                             # 中身を確認
```

### W1.6.10 QC後のReverse側
　「表1、(b)QC後、ユニーク配列の数、Reverse」の基礎情報です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
uniseq_filtRs <- sapply(filtRs, getUniNum)# filtRs中のサンプルごとのユニーク配列数をuniseq_filtRsに格納
uniseq_filtRs                             # 中身を確認
```

# W2 ノイズ除去
　[dada2](https://bioconductor.org/packages/release/bioc/html/dada2.html)の主要部分は、ノイズ除去と重複除去（denoising and dereplication）の2つのステップから構成されます。ここでは、その最初のステップを行います。シークエンサには癖があり、例えば「本当はAという塩基なのに間違ってTを出力しがち」とかがあります。このとき、「AからTへの置換（A to T）」というのを、この業界では「A2T」のようにtoを2に置き換えて表現する場合が多いです（特にLinuxの変換系のプログラム名でそのような傾向があります）。このような組み合わせは、理論上$4^2 = 16$通りあります（A2A, A2C, A2G, A2T, C2A, ..., T2G, T2T）。また、塩基ごとに付与されているクオリティスコア$q$が低いほどシークエンスエラーが起こりやすい（エラー率$p$が高い）のですが、この$p$と$q$の関係は以下の数式で表現されます：
$$
q = - 10 \times \log_{10} (p)
$$

## W2.1 `learnErrors`関数
　`learnErrors`関数は、リードファイルを入力として、Forward側とReverse側独立にエラーモデルの推定を行います。通常はQC後のfastqファイルが入力になりますので、例えばForward側の場合は、それらのファイルの相対パス情報である`filtFs`を与えます。`fls`はただのオプション名です。ここでは、得られたエラーモデル情報をForward側は`errFs`に、そしてReverse側は`errRs`オブジェクトに格納しています。
```{r, message=FALSE, error=TRUE, warning=FALSE}
errFs <- learnErrors(fls = filtFs, multithread = TRUE) # filtFsの一部を用いたエラーモデル推定
errRs <- learnErrors(fls = filtRs, multithread = TRUE) # filtRsの一部を用いたエラーモデル推定
```

## W2.2 `plotErrors`関数
　`plotErrors`関数は、`learnErrors`関数の実行結果を入力として、16通りの塩基置換の種類ごとにクオリティスコアとエラー頻度（推定エラーモデルを含む）の関係をプロットした結果を出力します。

### W2.2.1 Forward側
　Forward側の`learnErrors`関数の実行結果である`errFs`を入力としてプロットしています。`nominalQ = TRUE`は、赤線をつけよという指令です。これが図2と同じものになります。
```{r, message=FALSE, error=TRUE, warning=FALSE, fig.dim = c(8, 7)}
plotErrors(errFs, nominalQ = TRUE)        # errFsの中身（エラー頻度やエラーモデル推定結果）をプロット
```

### W2.2.2 Reverse側
　Reverse側の`learnErrors`関数の実行結果である`errRs`を入力としてプロットしています。
```{r, message=FALSE, error=TRUE, warning=FALSE, fig.dim = c(8, 7)}
plotErrors(errRs, nominalQ = TRUE)        # errRsの中身（エラー頻度やエラーモデル推定結果）をプロット
```

## W2.3 `dada`関数
　`dada`関数は、前処理が終わったリードファイル（Forward側の場合は`filtFs`）と先ほど得たエラーモデル情報（Forward側の場合は`errFs`）を入力として、denoisingを実行します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
dadaFs <- dada(derep = filtFs, err = errFs, multithread = TRUE) # filtFsとerrFsを入力としてエラー補正を実行
dadaRs <- dada(derep = filtRs, err = errRs, multithread = TRUE) # filtRsとerrRsを入力としてエラー補正を実行
```

## W2.4 結果の確認
　`head`関数は、行列の場合は最初の6行分、ベクトルやリストの場合は、最初の6つ分の要素を出力します。つまりデフォルトは`n = 6`です。ここでは`n = 2`オプションをつけて先ほどの`dada`関数実行結果（`dadaFs`および`dadaRs`）の中身を最初の2サンプル分に限定して確認しています。
```{r, message=FALSE, error=TRUE, warning=FALSE}
head(dadaFs, n = 2)                       # dadaFs中の最初の2つを確認
head(dadaRs, n = 2)                       # dadaRs中の最初の2つを確認
```

## W2.5 ユニーク配列数を確認
　W1.6.3で定義した自作関数`getUniNum`をここで改めて再定義して、`dada`関数実行後にユニーク配列の数がどれだけ減っているかを確認します。ここの数値が表2のユニーク配列数の基礎情報となります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
getUniNum <- function(x) length(getUniques(x))# ユニーク配列数を得るための関数を定義
sapply(dadaFs, getUniNum)                 # dadaFs中のサンプルごとのユニーク配列数を確認
sapply(dadaRs, getUniNum)                 # dadaRs中のサンプルごとのユニーク配列数を確認
```

## W2.6 リード数を確認
　W1.6.3の自作関数`getUniNum`はユニーク配列数を得るためのものでしたが、`length`関数部分を`sum`関数に置き換えることで、リード数を得るための関数に書き替えることができます。ここでは、それを`getReadNum`という自作関数として定義しています。ここの数値が表2のリード数の基礎情報となります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
getReadNum <- function(x) sum(getUniques(x))# リード数を得るための関数を定義
sapply(dadaFs, getReadNum)                # dadaFs中のサンプルごとのリード数を確認
sapply(dadaRs, getReadNum)                # dadaRs中のサンプルごとのリード数を確認
```

## W2.7 ユニーク配列の出現頻度
　第23回の図2でも示されているように、dada2は出現頻度の少ないユニーク配列を解析の過程で除去します。それを確かめるべく、`dada`関数実行結果（`dadaFs`と`dadaRs`）に対して、ユニーク配列ごとの出現頻度の主に下位を確認します。以前行ったW1.6.2のスクリプトを参考にしています。

### W2.7.1 出現頻度の上位
　まずは1番目のサンプル（サンプルID：SRR6325813）に対して、上位4つを眺めます。1つ1つのリードが長いので見づらいかと思いますが、1番多いのは`r getUniques(dadaFs[[1]])[1]`個、2番目に多いのは`r getUniques(dadaFs[[1]])[2]`個と解釈します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
head(getUniques(dadaFs[[1]]), n = 4)      # dadaFs中の1番目の要素に対して、ユニーク配列の最初の4つ（塩基配列と出現頻度）を確認
```

### W2.7.2 出現頻度の下位
　次に、下位4つを確認します。さきほどの`head`関数の結果から、`getUniques`関数実行結果は出現頻度でソートされているものだと思っていましたが、そうはなっていないことがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
tail(getUniques(dadaFs[[1]]), n = 4)      # dadaFs中の1番目の要素に対して、ユニーク配列の最後の4つ（塩基配列と出現頻度）を確認
```

　このような場合は、例えばW1.6.1などを参考にして、自分で降順（多 $\rightarrow$ 少）にソートしてから、改めて`tail`関数を実行するなどして対処します。W1.6.2の時点では、このサンプル（`dadaFs[[1]]`）は出現頻度が1回のものが存在しましたが、`dada`関数実行後はなくなっていることがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- sort(getUniques(dadaFs[[1]]), decreasing = T) # dadaFs中の1番目の要素に対して、ユニーク配列の出現頻度が多い順にソートした結果をhogeに格納
tail(hoge, n = 4)                         # hoge中の最後の4つのを確認
```

　ここで知りたいことは、主に出現頻度が最小のものが1ではなくなっているかどうかです。それゆえ、例えば以下のような感じで、`min`関数を用いて出現頻度の最小値のみを得るための`getMinNum`関数を自作します。W1.6.3、W2.5、およびW2.6を参考にしています。
```{r, message=FALSE, error=TRUE, warning=FALSE}
getMinNum <- function(x) min(getUniques(x))# 出現頻度の最小値を得るための関数を定義
```

　あとは、W2.5やW2.6と同じようにして`dada`関数実行結果（`dadaFs`と`dadaRs`）を入力として実行すればよいです。Forward側Reverse側含め全サンプルで2になっていることから、確かにdada2は出現頻度が2以上あるリードのみに限定していることがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
sapply(dadaFs, getMinNum)                 # dadaFs中のサンプルごとのリード数の最小値を取得
sapply(dadaRs, getMinNum)                 # dadaRs中のサンプルごとのリード数の最小値を取得
```

# W3 重複除去
　ここではまず、`mergePairs`関数を用いてForwardとReverseのオーバーラップ領域のマージ（重複領域の除去）を行います。この作業はサンプルごとに行われるため、得られたユニーク配列（ASV）の一部はサンプル間で重複します。そのため、`makeSequenceTable`関数を用いて解析対象の計12サンプル全体でASVの重複除去を行います。これは、和集合を得るイメージになります。

## W3.1 `mergePairs`関数
　`mergePairs`関数は、`dada`関数実行結果の2つのオブジェクト（`dadaFs`と`dadaRs`）、および`filterAndTrim`関数実行結果の2つのオブジェクト（`filtFs`と`filtRs`）を入力として、ペア―ドエンドリードのマージを実行します。ここではデフォルトのオプションを明示していますが、Forward側とReverse側で12塩基以上のオーバーラップがあり（`minOverlap = 12`）、かつオーバーラップ領域のミスマッチがない（`maxMismatch = 0`）という条件でマージが実行されます。
```{r, message=TRUE, error=TRUE, warning=FALSE}
mergers <- mergePairs(                    # ペア―ドエンドのマージを実行し、結果をmergersに格納
  dadaF = dadaFs,                         # エラー補正後のForward側はdadaFs
  derepF = filtFs,                        # QC後のForward側はfiltFs
  dadaR = dadaRs,                         # エラー補正後のReverse側はdadaRs
  derepR = filtRs,                        # QC後のReverse側はfiltRs
  minOverlap = 12,                        # オーバーラップの閾値は12塩基以上
  maxMismatch = 0,                        # オーバーラップ領域の許容ミスマッチ数は0
  verbose = TRUE                          # 結果を冗長に表示する（要は画面にも出力せよという指令）
)
```

## W3.2 サンプルごとの結果
　ここでは、`mergePairs`関数実行結果として得られた`mergers`オブジェクトの中身を確認します。`mergers[[1]]`が1番目のサンプルの結果、`mergers[[2]]`が2番目のサンプルの結果、...、`mergers[[12]]`が12番目のサンプルの結果になります。

### W3.2.1 結果の概要を確認
　1番目のサンプル（SRR6325813）のForwardとReverseをマージした結果について、`head`関数で最初の6配列分の結果を示します。前半部分で見えているのは、マージ後の塩基配列（コンティグ）です。後半部分で行列っぽく見えているのは、abundance列がコンティグの出現頻度、forward列が`dadaFs`オブジェクト中の配列のインデックス（何番目の配列かという情報）、reverse列が`dadaRs`オブジェクト中の配列のインデックス、nmatch列が一致塩基数です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
head(mergers[[1]])                        # mergers中の1番目の要素(サンプル)について、最初の6つを確認
```

### W3.2.2 手作業でマージ
　どのようなマージが実際に行われているのかを確認すべく、上記の結果でabundanceが`r mergers[[1]]$abundance[5]`個の5番目の配列を例として示します。まず、入力配列のForward側は`dadaFs`オブジェクト中の`r mergers[[1]]$forward[5]`番目の配列、Reverse側は`dadaRs`オブジェクト中の`r mergers[[1]]$reverse[5]`番目の配列だと解釈します。つまり、以下の2つがマージ前の配列です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
dadaFs[[1]]$sequence[5]                   # dadaFs中の1番目のサンプル中の5番目の塩基配列を確認
dadaRs[[1]]$sequence[6]                   # dadaRs中の1番目のサンプル中の6番目の塩基配列を確認
```

　これらの配列長は、W1.3でも指定していますが、Forward側が`r nchar(dadaFs[[1]]$sequence[5])` bp、Reverse側が`r nchar(dadaRs[[1]]$sequence[6])` bpです。W3.2のnmatch列の5番目を眺めると、オーバーラップ（一致）塩基数が`r mergers[[1]]$nmatch[5]`個だということが分かります。このことから、ペア―ドエンドのマージによって、`r nchar(dadaFs[[1]]$sequence[5])` bpあるForward側のリードの右側に、Reverse側の5'端の(`r nchar(dadaRs[[1]]$sequence[6])` - `r mergers[[1]]$nmatch[5]`) = `r nchar(dadaRs[[1]]$sequence[6]) - mergers[[1]]$nmatch[5]` bp分だけ追加されたのだろうと解釈します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
nchar(dadaFs[[1]]$sequence[5])            # dadaFs中の1番目のサンプル中の5番目の塩基配列の文字数(配列長)を確認
nchar(dadaRs[[1]]$sequence[6])            # dadaRs中の1番目のサンプル中の6番目の塩基配列の文字数(配列長)を確認
```

　ただし、実際にマージされるのは、Reverse側の配列の逆相補鎖（reverse complement）の3'端側の`r nchar(dadaRs[[1]]$sequence[6]) - mergers[[1]]$nmatch[5]` bp分だという点に注意が必要です。`rc`関数を利用すれば、逆相補鎖を得ることができます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
rc(dadaRs[[1]]$sequence[6])               # dadaRs中の1番目のサンプル中の6番目の塩基配列の逆相補鎖を確認
```

### W3.2.3 リード数を確認
　W2.6で自作したリード数を得るための関数`getReadNum`をここで改めて定義し、W3.1で得た`mergePairs`関数実行結果である`mergers`オブジェクトを入力として実行します。ここの数値が表3のリード数の基礎情報となります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
getReadNum <- function(x) sum(getUniques(x))# リード数を得るための関数を定義
sapply(mergers, getReadNum)               # mergers中のサンプルごとのリード数を確認
```

### W3.2.4 ユニーク配列数を確認
　W2.5で自作したリード数を得るための関数`getUniNum`をここで改めて定義し、W3.1で得た`mergePairs`関数実行結果である`mergers`オブジェクトを入力として実行します。ここの数値が表3のユニーク配列数の基礎情報となります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
getUniNum <- function(x) length(getUniques(x))# ユニーク配列数を得るための関数を定義
sapply(mergers, getUniNum)                    # サンプルごとのユニーク配列数
```


## W3.3 サンプル全体の結果
　さきほど得たユニーク配列の数は、サンプルごとに得たものです。それゆえ、計12サンプル全体でみると、まだ重複が含まれます。ここでは、サンプルごとではなく、サンプル全体でのユニーク配列数の情報を得ることを目的とします。

### W3.3.1 `getSequences`関数
　`getSequences`関数は、入力オブジェクトから、（塩基）配列をベクトルとして抽出した結果を返します。W3.2.1で概要を確認しましたが、W3.1で得た`mergers`オブジェクトには、計12サンプルのペア―ドエンドのマージ後のユニーク配列だけでなく、その出現頻度など様々な情報が格納されています。それゆえ、まずは`mergers`オブジェクト内の塩基配列情報のみを`getSequences`関数を用いて抽出し、その結果を`hoge`オブジェクトに格納します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- sapply(mergers, getSequences)     # mergers中のサンプルごとの塩基配列をhogeに格納
```

### W3.3.2 `class`関数
　得られた`hoge`オブジェクトは、サンプルごとに切り分けられたリスト形式のオブジェクトであることが、それを確認するための`class`関数実行結果からわかります。つまり、表示結果が"list"になっていれば、それはリスト形式のオブジェクトだということです。
```{r, message=FALSE, error=TRUE, warning=FALSE}
class(hoge)                               # hogeオブジェクトのクラスを確認
```

### W3.3.3 `unlist`と`is.vector`関数
　W1.6.5でも使いましたが、リスト形式をベクトル形式に変換したい場合には、`unlist`関数を利用します。得られた`hoge2`オブジェクトがベクトルか否かをTRUE or FALSEで返す`is.vector`関数実行結果からも、確かにベクトル形式に変換できていることがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge2 <- unlist(hoge)                     # リスト形式をベクトル形式に変換
is.vector(hoge2)                          # ベクトルか否かを表示
```

### W3.3.4 `sapply`と自作関数
　`hoge2`ベクトルの要素数は、W3.2.3で得たサンプルごとのユニーク配列数の総和と一致するはずです。ここでは確かにそうであることを確認しているだけです。
```{r, message=FALSE, error=TRUE, warning=FALSE}
length(hoge2)                             # hoge2ベクトルの要素数を確認
sum(sapply(mergers, getUniNum))           # mergers中のサンプルごとのユニーク配列数の総和(つまり総リード数)を確認
```

### W3.3.5 ユニーク配列の和集合
　サンプル全体のユニーク配列数は、以下のように`unique`関数を実行した上で、`length`関数を実行すればよいです。この`r length(unique(hoge2))`という数字は、ユニーク配列の和集合（union）ともいえます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
length(unique(hoge2))                     # hoge2ベクトル中のユニークな要素数を確認
```

　上記のやり方以外に、`table`関数を用いてユニーク配列ごとの出現頻度を調べると、`table`関数実行結果の要素数が「サンプル全体のユニーク配列数」に相当することを利用することもできます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
length(table(hoge2))                      # hoge2ベクトル中のユニークな要素数を確認
```

### W3.3.6 ユニーク配列の積集合
　ほぼ余談ですが、ユニーク配列の和集合（union）が分かると、計12サンプル全てに存在するユニーク配列の数（積集合；intersection）も知りたくなると思います。ここではユニーク配列ごとの出現頻度を調べる`table`関数実行結果から、出現頻度が12回（サンプル数の合計が12個だからです）のもののみを抽出するやり方で示します。実行結果は`r sum(table(hoge2) == 12)`個です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
sum(table(hoge2) == 12)                   # hoge2ベクトル中の出現頻度が12回のユニーク配列数を確認
```

### W3.3.7 スクリプトを一まとめ
　ユニーク配列の和集合（`r length(unique(hoge2))`個）と積集合（`r sum(table(hoge2) == 12)`個）を得るためのスクリプトをひとまとめにしたのが以下になります。この一連の流れを理解しておくと、関数が無数に入れ子になっていても、ステップバイステップで理解していけばいいのだと安心できると思います。（本当は最終行の12という数値も汎用性の高いものにするとよいのですが、ここでは見た目の分かりやすいさを優先しています。）
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- sapply(mergers, getSequences)     # mergers中のサンプルごとの塩基配列をhogeに格納
hoge2 <- unlist(hoge)                     # hogeをベクトル形式に変換してhoge2に格納
length(unique(hoge2))                     # hoge2ベクトル中のユニークな要素数(ユニーク配列の和集合)を確認
sum(table(hoge2) == 12)                   # hoge2ベクトル中の出現頻度が12回のユニーク配列数(ユニーク配列の積集合)を確認
```

## W3.4 `makeSequenceTable`関数
　`makeSequenceTable`関数は、`mergePairs`関数実行結果を入力として、解析サンプル中のユニーク配列の和集合を得るとともに、どのサンプルにどのユニーク配列がいくつ存在するかという行列を結果として返します。行数はサンプル数、列数はサンプル全体でのユニーク配列の数になります。

### W3.4.1 まずは実行
　先ほど得た`mergePairs`関数実行結果である`mergers`オブジェクトを入力として実行し、結果を`seqtab`オブジェクトに格納するところまで行います。
```{r, message=FALSE, error=TRUE, warning=FALSE}
seqtab <- makeSequenceTable(mergers)      # mergers中のユニーク配列の和集合情報を取得し、サンプルごとのユニーク配列の出現頻度をseqtabに格納
```

### W3.4.2 結果の概要を確認
　結果は行列ですので、まずは`dim`関数で行数と列数を表示させます。行数は`r nrow(seqtab)`、列数は`r ncol(seqtab)`であることがわかります。このことから、今回解析する計`r nrow(seqtab)`サンプル全体としては、`r ncol(seqtab)`個のユニークな配列が存在するのだと解釈します。W3.3.5と同じ結果が得られていることがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
dim(seqtab)                               # seqtabの行数と列数を確認
```

### W3.4.3 `rowSums`と`sum`関数
　`rowSums`関数は、数値行列を入力として、行（row）ごとの和（sum）を返します。得られる結果はW3.2.3と同じです。同じノリで行列全体の総和を`sum`関数で計算すると、サンプル全体の総リード数が得られます。得られる結果（`r sum(seqtab)`個）は表2bのマージ後のリード数の総和と一致します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
rowSums(seqtab)                           # seqtabの行ごと(サンプルごと)の和を確認
sum(seqtab)                               # seqtab中の総配列数を確認
```

### W3.4.4 `getSequences`関数
　`getSequences`関数は、入力オブジェクトから、（塩基）配列をベクトルとして抽出した結果を返します（W3.3.1）。`r ncol(seqtab)`個のユニークな塩基配列を得たい場合は、`seqtab`オブジェクトを入力として`getSequences`関数を実行します。直後の`length(hoge)`実行結果から、確かに`hoge`の中身はr ncol(seqtab)`個の要素からなるベクトルだということがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- getSequences(seqtab)              # seqtab中の塩基配列をhogeに格納
length(hoge)                              # hoge中の要素数を確認
```

### W3.4.5 `nchar`関数
　`nchar`関数は、文字列ベクトルを入力として、要素ごとの文字数を返します（W3.2.2）。したがって、この場合の`hoge2`オブジェクトの要素数は、入力である文字列ベクトルの要素数と同じ（`r length(hoge)`）になります。`head`関数で最初の6つ分の中身を眺めると、図3で得た213 bpというコンティグ長と酷似した結果ばかり表示されていることがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge2 <- nchar(hoge)                      # hoge中の要素ごとの文字数をhoge2に格納
length(hoge2)                             # hoge2中の要素数を確認
head(hoge2)                               # hoge2中の最初の6つを確認
```

### W3.4.6 `table`関数
　実質的に先ほどまでの結果に`table`関数を実行しているだけです。`table`関数は、入力ベクトル中の要素の種類ごとの出現頻度を返します。したがって、この結果は、`r ncol(seqtab)`個のユニークな塩基配列の長さ分布情報を得ていることと同義です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- getSequences(seqtab)              # seqtab中の塩基配列をhogeに格納
hoge2 <- nchar(hoge)                      # hoge中の要素ごとの文字数をhoge2に格納
table(hoge2)                              # hoge2ベクトル中の要素の種類ごとに出現頻度をカウント
```

### W3.4.7 一気に書く
　一気に書くとこんな感じです。
```{r, message=FALSE, error=TRUE, warning=FALSE}
table(nchar(getSequences(seqtab)))        # seqtab中の塩基配列ごとの文字数をカウントし、その文字数の種類ごとの出現頻度を表示
```

　W3.3.7に似せて、以下のような書き方でも構いません。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- sapply(mergers, getSequences)     # mergers中のサンプルごとの塩基配列をhogeに格納
hoge2 <- unlist(hoge)                     # hogeをベクトル形式に変換してhoge2に格納
hoge3 <- nchar(unique(hoge2))             # hoge2中のユニークな要素ごとの文字数をhoge3に格納
table(hoge3)                              # hoge3ベクトル中の要素の種類ごとの出現頻度を確認
```

## W3.5 `removeBimeraDenovo`関数
　[dada2](https://bioconductor.org/packages/release/bioc/html/dada2.html)パッケージのメインであるdada関数は、エラー補正のみを行いキメラ配列の除去は行ってくれません。つまり、まだキメラ配列を含んでいます。`removeBimeraDenovo`関数は、`makeSequenceTable`関数の実行結果である`seqtab`オブジェクトを入力として実行し、PCRキメラを除去した結果を返します。

### W3.5.1 まずは実行
　さきほど得た`makeSequenceTable`関数の実行結果である`seqtab`オブジェクトを入力として、PCRキメラを除去した結果を`seqtab.nochim`オブジェクトに格納するところまで行います。
```{r, message=FALSE, error=TRUE, warning=FALSE}
seqtab.nochim <- removeBimeraDenovo(seqtab,              # seqtab中のキメラ配列除去した結果をseqtab.nochimに格納
                                    method = "consensus",# seqtab中のキメラ配列除去した結果をseqtab.nochimに格納
                                    multithread = TRUE,  # seqtab中のキメラ配列除去した結果をseqtab.nochimに格納
                                    verbose = TRUE)      # seqtab中のキメラ配列除去した結果をseqtab.nochimに格納
```

### W3.5.2 結果の概要を確認
　結果は行列ですので、まずは`dim`関数で行数と列数を表示させます。行数は`r nrow(seqtab.nochim)`、列数は`r ncol(seqtab.nochim)`であることがわかります。このことから、今回解析する計`r nrow(seqtab.nochim)`サンプル全体としては、`r ncol(seqtab.nochim)`個のユニークな配列が存在するのだと解釈します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
dim(seqtab.nochim)                        # seqtab.nochimの行数と列数を確認
```

### W3.5.3 `rowSums`と`sum`関数
　`rowSums`関数は、数値行列を入力として、行（row）ごとの和（sum）を返します。得られる結果はサンプルごとのリード数になります。これが同じノリで行列全体の総和を`sum`関数で計算すると、サンプル全体の総リード数が得られます。表2bのペア―ドエンドのマージ後のリード数からほとんど減っていないこともよくわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
rowSums(seqtab.nochim)                    # seqtab.nochimの行ごとの総和(除去後のサンプルごとのリード数)を表示
sum(seqtab.nochim)                        # seqtab.nochim中の全要素の総和(除去後の総リード数)を表示  
sum(seqtab)                               # seqtab中の全要素の総和(除去前の総リード数)を表示
sum(seqtab.nochim)/sum(seqtab)            # 除去後の総リード数/除去前の総リード数の割合
```

### W3.5.4 配列長分布と出現頻度
　得られた計`r ncol(seqtab.nochim)`個の配列長分布と出現頻度情報を得ます。これはW3.4.7と同じ書き方です。キメラ配列除去前と酷似した傾向であることがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
table(nchar(getSequences(seqtab.nochim))) # seqtab.nochim中の塩基配列ごとの文字数をカウントし、その文字数の種類ごとの出現頻度を表示
```

## W3.6 図4のスクリプト
　コンパクトにまとめるとこんな感じで書けます。
```{r, message=FALSE, error=TRUE, warning=FALSE, eval=FALSE}
mergers <- mergePairs(dadaF = dadaFs, derepF = filtFs,# ペア―ドエンドのマージを実行し、結果をmergersに格納
                      dadaR = dadaRs, derepR = filtRs)# ペア―ドエンドのマージを実行し、結果をmergersに格納

seqtab <- makeSequenceTable(mergers)      # mergers中のユニーク配列の和集合情報を取得し、サンプルごとのユニーク配列の出現頻度をseqtabに格納
dim(seqtab)                               # seqtabの行数と列数を確認
sum(seqtab)                               # seqtab中の総配列数を確認
table(nchar(getSequences(seqtab)))        # seqtab中の塩基配列ごとの文字数をカウントし、その文字数の種類ごとの出現頻度を表示

seqtab.nochim <- removeBimeraDenovo(seqtab, method = "consensus")# seqtab中のキメラ配列除去した結果をseqtab.nochimに格納
dim(seqtab.nochim)                        # seqtab.nochimの行数と列数を確認
sum(seqtab.nochim)                        # seqtab.nochim中の全要素の総和(除去後の総リード数)を表示
table(nchar(getSequences(seqtab.nochim))) # seqtab.nochim中の塩基配列ごとの文字数をカウントし、その文字数の種類ごとの出現頻度を表示
```

# W4 系統割り当て
　最後のステップでは、まず、系統推定のためのデータベースに相当するリファレンス配列セットを取得します。次に、`assignTaxonomy`関数を用いて属（genus）の推定を行います。最後に`addSpecies`関数を用いて種（species）の推定を行います。

## W4.1 リファレンス配列取得
　本稿は、[DADA2 Pipeline Tutorial (1.16)](https://benjjneb.github.io/dada2/tutorial.html)の流れに沿って説明しています。この中のAssign taxonomyという部分で、既知の16S rRNAのリファレンス配列を取得するためのページとして[https://benjjneb.github.io/dada2/training.html](https://benjjneb.github.io/dada2/training.html)へのリンクが示されています。リンク先で見られる[Silva 138.1 prokaryotic SSU taxonomic training data formatted for DADA2](https://doi.org/10.5281/zenodo.4587954)というページから取得可能な以下の2つのファイルを利用します。

- 属レベルの系統推定用：[silva_nr99_v138.1_train_set.fa.gz](https://zenodo.org/records/4587955/files/silva_nr99_v138.1_train_set.fa.gz)
- 種レベルの系統推定用：[silva_species_assignment_v138.1.fa.gz](https://zenodo.org/records/4587955/files/silva_species_assignment_v138.1.fa.gz)

## W4.2 `assignTaxonomy`関数
　`assignTaxonomy`関数は、系統推定したい塩基配列の文字列ベクトルと属レベルまでの情報を含むリファレンス配列ファイルを入力として、系統推定結果を返します。

### W4.2.1 まずは実行
　属レベルの系統推定を行います。入力のクエリ（問い合わせ配列）側は、W3.5.1でキメラ配列除去まで終わった後の`seqtab.nochim`オブジェクトです。入力のデータベース側は、"./silva_nr99_v138.1_train_set.fa.gz"のようにして与えています。これは、作業ディレクトリ（デスクトップのdatasetフォルダ）直下にある、silva_nr99_v138.1_train_set.fa.gzというファイルだという意味です。ダウンロードしたファイルを2行目のように指定してもよいですが、（1番左側に#を入れてコメントアウトしているので機能はしませんが...）1行目で示すように、オリジナルのURLを直接指定することもできます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
#fastafile1 <- "https://zenodo.org/records/4587955/files/silva_nr99_v138.1_train_set.fa.gz"
fastafile1 <- "./silva_nr99_v138.1_train_set.fa.gz" # リファレンス配列のパス情報をfastafile1に格納
taxa <- assignTaxonomy(seqs = seqtab.nochim, # seqtab.nochimをクエリ配列、fastafile1をリファレンス配列として系統割り当てを実行し、結果をtaxaに格納
                       refFasta = fastafile1)# seqtab.nochimをクエリ配列、fastafile1をリファレンス配列として系統割り当てを実行し、結果をtaxaに格納
```

### W4.2.2 結果の概要を確認
　結果は行列ですので、まずは`dim`関数で行数と列数を表示させます。行数は`r nrow(taxa)`、列数は`r ncol(taxa)`であることがわかります。次は、`colnames`関数で`taxa`行列の列名を表示させています。この結果を眺めて、生物の階級を「界（Kingdom）・門（Phylum）・綱（Class）・目（Order）・科（Family）・属（Genus）・種（Species）」と表すと、属までに6つの階級が存在するからこのようになっているのだと納得できます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
dim(taxa)                                 # taxaの行数と列数を表示
colnames(taxa)                            # taxaの列名を表示
```

## W4.3 `addSpecies`関数
　`addSpecies`関数は、さきほどの属レベルの系統推定結果である`taxa`オブジェクト、および種レベルの系統推定用のリファレンス配列ファイル（[silva_species_assignment_v138.1.fa.gz](https://zenodo.org/records/4587955/files/silva_species_assignment_v138.1.fa.gz)）を入力として与えて、種レベルの系統推定結果情報を足したものを返します。

### W4.3.1 まずは実行
　以下のような感じで実行します。W4.2.1と同様、（1番左側に#を入れてコメントアウトしているので機能はしませんが…）1行目で示すように、オリジナルのURLを直接指定することもできます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
# fastafile2 <- "https://zenodo.org/records/4587955/files/silva_species_assignment_v138.1.fa.gz"
fastafile2 <- "./silva_species_assignment_v138.1.fa.gz" # リファレンス配列のパス情報をfastafile12に格納
taxa2 <- addSpecies(taxtab = taxa,        # taxaをクエリ配列、fastafile2をリファレンス配列として系統割り当てを実行し、結果をtaxa2に格納
                    refFasta = fastafile2)# taxaをクエリ配列、fastafile2をリファレンス配列として系統割り当てを実行し、結果をtaxa2に格納
```

### W4.3.2 結果の概要を確認
　結果は行列ですので、まずは`dim`関数で行数と列数を表示させます。行数は`r nrow(taxa2)`、列数は`r ncol(taxa2)`であることがわかります。次は、`colnames`関数で`taxa2`行列の列名を表示させています。生物の階級である「界（Kingdom）・門（Phylum）・綱（Class）・目（Order）・科（Family）・属（Genus）・種（Species）」の種の結果が足されていることが分かります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
dim(taxa2)                                # taxa2の行数と列数を表示
colnames(taxa2)                           # taxa2の列名を表示
```

## W4.4 ファイル出力
　ここでは、W3.5.1で得た`seqtab.nochim`オブジェクトとW4.3.1で得た`taxa2`オブジェクトの内容をうまく連結して、タブ区切りテキストファイルとして保存する作業を行います。

### W4.4.1 行列の転置
　`seqtab.nochim`の行数と列数はそれぞれ`r nrow(seqtab.nochim)`と`r ncol(seqtab.nochim)`です。その一方で連結対象である`taxa2`の行数と列数は、それぞれ`r nrow(taxa2)`と`r ncol(taxa2)`です。従って、うまく連結すべく、`seqtab.nochim`行列のほうを`t`関数を用いて転置して、同一行が同一のASVの情報になるようにできるかどうか確認しています。
```{r, message=FALSE, error=TRUE, warning=FALSE}
dim(seqtab.nochim)                        # seqtab.nochimの行数と列数を表示
dim(taxa2)                                # taxa2の行数と列数を表示
dim(t(seqtab.nochim))                     # seqtab.nochim行列を転置して、行数と列数を表示
```

### W4.4.2 ASV配列情報の取得
　ここでは、最終的に`t(seqtab.nochim)`と`taxa2`を列方向で連結させるのですが、ちょっと小細工をしてASVの塩基配列情報のみを文字列ベクトルとして得た結果を`ASV`というオブジェクトとして得ようとしています。ここでは、`taxa2`オブジェクトの行名情報を`ASV`に格納するほうを有効化していますが、#でコメントアウトしている`seqtab.nochim`オブジェクトの列名情報を`ASV`に格納するのでも構いません。そのあとは、`ASV`ベクトルの要素数を`length`関数で確認したり、最初の3つの要素（つまりASVの塩基配列情報）を表示して確認しているだけです。
```{r, message=FALSE, error=TRUE, warning=FALSE}
#ASV <- colnames(seqtab.nochim)           # seqtab.nochimの列名をASVに格納(するのでもよい)
ASV <- rownames(taxa2)                    # taxa2の行名をASVに格納
length(ASV)                               # ASVの要素数を確認
head(ASV, n = 3)                          # ASV中の最初の6つを表示
```

### W4.4.3 `cbind`と`write.table`関数
　`cbind`関数を用いてファイル出力させたい情報を列方向で連結し、それを`write.table`関数でファイルに保存しています。図5は、得られたファイルをExcelで開いてスクリーンショットをとったものになります。念のため、ここでもタブ区切りテキストファイル版（[output_JSLAB24.txt](https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.txt)）とxlsx版（[output_JSLAB24.xlsx](https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.xlsx)）をダウンロードできるようにしておきます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
tmp <- cbind(ASV, t(seqtab.nochim), taxa2)# ASV、seqtab.nochim行列の転置、taxa2を列方向で連結してtmpに格納
dim(tmp)                                  # tmpの行数を列数を表示
write.table(x = tmp, file = "output_JSLAB24.txt", sep = "\t", row.names = F)# tmpの中身を指定したファイル名でタブ区切りテキストファイルとして、行名をつけずに保存
```

### W4.4.4 `identical`関数
　さきほどの`cbind`の部分で`t(seqtab.nochim)`と`taxa2`が本当に同一のASV配列になっているか不安な方は、`seqtab.nochim`行列の列名情報である`colnames(seqtab.nochim)`と`taxa2`の行名情報が順番も含めて同一であるかを以下のコマンドで確認しておくとよいです。`identical`関数は、TRUEなら同じ、FALSEなら違っていると解釈します。先ほどのW4.4.3でいきなり`cbind`関数で列同士を連結している印象をもったかもしれませんが、実際には行名部分の一部が同じであることを目視確認したり、以下で示すように`identical`関数などを用いてTRUEという結果を得てからにしたり、あるいはもっと高度な方法で連結させます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
identical(colnames(seqtab.nochim), rownames(taxa2))# seqtab.nochimの列名とtaxa2の行名の2つのベクトルが同一かどうかを確認
```

　いきなり上記のやつだと挙動がよくわからないかもしれませんので、念のため以下に簡単な例を示しておきます。`hoge2`と`hoge3`の比較時のみ`r identical(hoge2, hoge3)`になっていることがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge1 <- c("G", "H", "U")                 # 3つの要素からなる文字列ベクトルを作成してhoge1に格納
hoge2 <- c("G", "H", "U", "K")            # 4つの要素からなる文字列ベクトルを作成してhoge2に格納
hoge3 <- c("G", "H", "U", "K")            # 4つの要素からなる文字列ベクトルを作成してhoge3に格納
hoge4 <- c("H", "U", "G")                 # 3つの要素からなる文字列ベクトルを作成してhoge4に格納
identical(hoge1, hoge2)                   # hoge1とhoge2の2つのベクトルが同一かどうかを確認
identical(hoge1, hoge3)                   # hoge1とhoge3の2つのベクトルが同一かどうかを確認
identical(hoge1, hoge4)                   # hoge1とhoge4の2つのベクトルが同一かどうかを確認
identical(hoge2, hoge3)                   # hoge2とhoge3の2つのベクトルが同一かどうかを確認
identical(hoge2, hoge4)                   # hoge2とhoge4の2つのベクトルが同一かどうかを確認
identical(hoge3, hoge4)                   # hoge3とhoge4の2つのベクトルが同一かどうかを確認
```

# おわりに

## W5.1 参考にしたページ
- [DADA2 Pipeline Tutorial (1.16)](https://benjjneb.github.io/dada2/tutorial.html)
- [dada2](https://bioconductor.org/packages/dada2/)の[Introduction to dada2](https://www.bioconductor.org/packages/release/bioc/vignettes/dada2/inst/doc/dada2-intro.html)ページ

## W5.2 解析環境の確認
　`packageVersion`関数は、個別のインストール済みパッケージのバージョン情報を知りたい場合に利用します。今回動作確認したdada2のバージョンは`r packageVersion("dada2")`だということがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
packageVersion("dada2")                   # dada2のバージョンを表示
```

　もっと広い意味での解析環境情報を知りたい場合は、`sessionInfo()`の出力結果に基づいて記載します。このPCの場合は、R本体のバージョンが4.3.2であることや、OSがWindows 11であることなどがわかります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
sessionInfo()                             # 解析環境を表示
```