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

# はじめに
  - [乳酸菌NGS連載第22回のPDF](https://www.iu.a.u-tokyo.ac.jp/kadota/JSLAB_22_kadota.pdf)
  - [JSLAB23.Rmd](https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB23.Rmd)  

# PCRプライマー

# Operational taxonomic unit (OTU)
  - [Yarza et al., Nat Rev Microbiol
., 2014](https://pubmed.ncbi.nlm.nih.gov/25118885/)  
  - [MAPseq](https://github.com/jfmrod/mapseq)：[Matias Rodrigues et al., Bioinformatics, 2017](https://pubmed.ncbi.nlm.nih.gov/28961926/)  
クラスタリングを行わずに直接統割り当てを行うアプローチです。  
  - [VITCOMIC2](http://vitcomic.org/)：[Mori et al., BMC Syst Biol., 2018](https://pubmed.ncbi.nlm.nih.gov/29560821/)  
クラスタリングを行わずに直接統割り当てを行うアプローチです。 
  - [SILVA](http://www.arb-silva.de)：[Quast et al., Nucleic Acids Res., 2013](https://pubmed.ncbi.nlm.nih.gov/23193283/)  
16S rRNAのリファレンスデータベースです。  

# Amplicon sequence variant (ASV)
  - [DADA2](https://bioconductor.org/packages/release/bioc/html/dada2.html)：[Callahan et al., Nat Methods, 2016](https://pubmed.ncbi.nlm.nih.gov/27214047/)  
ノイズ除去と重複除去（denoising & dereplication）を行った後、系統割り当てを行うアプローチです。  
  - [Deblur](https://github.com/biocore/deblur)：[Amir et al., mSystems, 2017](https://pubmed.ncbi.nlm.nih.gov/28289731/)  
  ノイズ除去と重複除去（denoising & dereplication）を行った後、系統割り当てを行うアプローチです。 

# W1：解析データ
- [ヒトマイクロバイオーム計画(Human Microbiome Project; HMP)](https://www.hmpdacc.org/)  
[NIH](https://www.nih.gov/)主導のプロジェクトです。  
- W1.1：[HM-280](https://www.beiresources.org/Catalog/otherProducts/HM-280.aspx) ([Gevers et al., PLoS Comput Biol., 2012](https://pubmed.ncbi.nlm.nih.gov/23209389/))  
表1の内容は、[Hallmaier-Wacker et al., Sci Rep., 2018](https://pubmed.ncbi.nlm.nih.gov/29674641/)の[Table S1](https://static-content.springer.com/esm/art%3A10.1038%2Fs41598-018-24573-y/MediaObjects/41598_2018_24573_MOESM1_ESM.pdf)、および[HM-280](https://www.beiresources.org/Catalog/otherProducts/HM-280.aspx)のページ内にある[Product Information Sheet](https://www.beiresources.org/ProductInformationSheet/tabid/784/Default.aspx?doc=38929.pdf)から構成されています。  
- W1.2：[Product Information Sheet](https://www.beiresources.org/ProductInformationSheet/tabid/784/Default.aspx?doc=38929.pdf)  
[BEI Resources](https://www.beiresources.org/)から辿れる、38929.pdfというファイルです。  
- W1.3：[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サンプル分のデータを取得していますが、以下の表2としても示されている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.4：[For20200428LectureDADA2.tar.gz](http://palaeo.nig.ac.jp/Resources/TodaiAgri20/For20200428LectureDADA2.tar.gz)  
表2で示されている12サンプル分のリードファイルをまとめたものです。PCR増幅は515F-806Rプライマーで行い、増幅された300 bp弱のアンプリコンをIllumina MiSeqのペア―ドエンド(250 bp×2)で得たリードです。1サンプルにつき、Forward側のファイル(\*_1.fastq)とReverse側のファイル(\*_2.fastq)の2つのファイルが得られますので、合計24ファイルで構成されます。

# DADA2のインストール
- [DADA2](https://bioconductor.org/packages/release/bioc/html/dada2.html)：[Callahan et al., Nat Methods, 2016](https://pubmed.ncbi.nlm.nih.gov/27214047/)  

## W2：インストール実行
```{r, message=FALSE, error=TRUE, warning=FALSE, eval=FALSE}
if (!requireNamespace("BiocManager", quietly=T))
    install.packages("BiocManager")
BiocManager::install("dada2", update=F)
```

## W3：パッケージのロード
```{r, message=FALSE, error=TRUE, warning=FALSE, eval=FALSE}
library(dada2)
```

## W4：`list.files`関数
### W4.1：基本形と`path`オプション
`list.files`関数は、指定したディレクトリ内にあるファイル情報を得るためのものです。1つ目のコマンドは、オプションを全てデフォルトで実行する場合です。2つ目は、表示させたい場所を指定するための`path`オプションのデフォルト"."を明示しただけです。両者は同じ結果になります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
list.files()
list.files(path = ".")
```

### W4.2：`pattern`オプション(相対パス)
`pattern`オプションをつけることで、任意の文字列を含むファイルのみに限定することができます。以下は`"_1.fastq"`という文字を指定して、Forward側に限定するやり方です。1つ目のコマンドは、`path`オプションを明示せずデフォルトで実行した結果です。2つ目のコマンドは、`path`オプションのデフォルト"."を明示しただけです。両者は同じ結果になります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
list.files(pattern = "_1.fastq")
list.files(pattern = "_1.fastq", path = ".")
```

### W4.3：`pattern`オプション(絶対パス)
現在のディレクトリの場所を表示する`getwd`関数実行結果(つまり絶対パス情報)を`basho`というオブジェクトに格納し、それを`list.files`関数の`path`オプションに与えています。
```{r, message=FALSE, error=TRUE, warning=FALSE}
basho <- getwd()
basho
list.files(pattern = "_1.fastq", path = basho)
```

### W4.4：`full.names`オプション
`full.names`オプションをTRUEにすると、`path`オプションで指定したパス情報を含めたファイル名になります。1つ目のコマンドが相対パス、2つ目のコマンドが絶対パスになります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
list.files(pattern = "_1.fastq", path = ".", full.names = TRUE)
list.files(pattern = "_1.fastq", path = basho, full.names = TRUE)
```

## W5：文字列操作
### W5.1 `basename`関数(相対パスの場合)
`basename`関数は、パス情報を含むファイル名の文字列ベクトルを入力として、パス情報を含まないファイル名の文字列ベクトルを得る際に用いられます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
fnFs <- list.files(pattern = "_1.fastq", path = ".", full.names = TRUE)
fnFs
basename(fnFs)
```

### W5.2 `basename`関数(絶対パスの場合)
`basename`関数は、パス情報を含むファイル名の文字列ベクトルを入力として、パス情報を含まないファイル名の文字列ベクトルを得る際に用いられます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
basho <- getwd()
fnFs <- list.files(pattern = "_1.fastq", path = basho, full.names = TRUE)
fnFs
basename(fnFs)
```

### W5.3 `strsplit`関数
`basename`関数は、パス情報を含むファイル名の文字列ベクトルを入力として、パス情報を含まないファイル名の文字列ベクトルを得る際に用いられます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
basename(fnFs)
strsplit(basename(fnFs), "_")
```

### W5.4 `sapply`関数
`sapply`関数は、リスト形式のオブジェクトに対して逐次的な処理を行いたい場合に用いられます。ここでは、さきほどの`strsplit`実行結果を`hoge`というオブジェクトに一旦格納したのち、それを`sapply`関数の入力として用いています。出力結果は、表2のサンプルIDに相当します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- strsplit(basename(fnFs), "_")
sapply(X = hoge, FUN = `[`, 1)
```

### W5.5 `paste0`関数
`paste0`関数は、文字列同士の連結を行いたいときに用います。ここでは、さきほどの`sapply`実行結果を`sampleID`というオブジェクトに一旦格納したのち、それを`paste0`関数の入力として用いています。文字列ベクトルである`sampleID`の各要素の右側に"_F_filt.fastq"という文字列を連結させています。ここまでが、なんらかの処理実行結果を出力させるための任意のファイル名を生成する基本形です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
sampleID <- sapply(X = hoge, FUN = `[`, 1)
paste0(sampleID, "_F_filt.fastq")
```

### W5.6 `paste`関数
`paste`関数は、文字列同士の連結を行いたいときに用います。`paste0`関数と基本的に同じですが、`paste`関数のほうが間に挟む文字(区切り文字)を指定できる点が異なります。`paste`関数のデフォルトは、区切り文字を指定する`sep`オプションの値がスペース(`" "`)です。従って、`paste0`関数は、`paste`関数を`sep = ""`オプションをつけて実行することと同義です。
```{r, message=FALSE, error=TRUE, warning=FALSE}
sampleID <- sapply(X = hoge, FUN = `[`, 1)
paste(sampleID, "_F_filt.fastq")
paste(sampleID, "_F_filt.fastq", sep = "")
```

### W5.7 `file.path`関数
`file.path`関数は、ファイルのパス情報をうまく生成したいときに用います。ここでは、現在のディレクトリ直下に"processed"というディレクトリを作成し、その中に先ほどの`paste0`実行結果に相当する処理後のファイル名を連結させたものを作成します。
```{r, message=FALSE, error=TRUE, warning=FALSE}
basho <- getwd()
fnFs <- list.files(pattern = "_1.fastq", path = basho, full.names = TRUE)
hoge <- strsplit(basename(fnFs), "_")
sampleID <- sapply(X = hoge, FUN = `[`, 1)
filtFs <- file.path(basho, "processed", paste0(sampleID, "_F_filt.fastq"))
filtFs
```
一見すると`file.path`関数の有用性がわかりにくいヒトは、以下で示すような`file.path`の部分を`paste0`に置き換えた結果と比較すると納得できると思います。つまり、`file.path`のほうがスラッシュ(\/)をうまく挿入してくれます。
```{r, message=FALSE, error=TRUE, warning=FALSE}
paste0(basho, "processed", paste0(sampleID, "_F_filt.fastq"))
```
`file.path`関数の代わりに、`sep = "/"`オプションをつけて`paste`関数を実行しても構いません。
```{r, message=FALSE, error=TRUE, warning=FALSE}
paste(basho, "processed", paste0(sampleID, "_F_filt.fastq"), sep = "/")
```

### W5.8 `names`関数
これは少しわかりにくい概念ですが、`names`関数はベクトルの各要素に名前をつけたいときに用います。実際取り扱うのはForward側とReverse側の両方ですが、通常これらはペアとして取り扱います。それゆえ、対応する要素の中身の文字列（つまりファイル名）で異なっていたとしても、要素につけられた名前を一致させておくことでペアとして認識させることができるメリットがあります。ここでは、要素につける名前としてサンプルID（オブジェクト名としてはsampleID）を割り当てています。
```{r, message=FALSE, error=TRUE, warning=FALSE}
filtFs
names(filtFs) <- sampleID
filtFs
```

## W6 入出力ファイルの作成
### W6.1 絶対パスの場合
ここまではForward側のファイルのみでしたので、Reverse側も含めて全体をまとめて実行するスクリプトを示します。スクリプトに本当は不要な部分も含まれますが、Forward側とReverse側ともに同じように記載することを優先させています。`filtFs`と`filtRs`が前処理後のファイルの絶対パス情報になります。
```{r, message=FALSE, error=TRUE, warning=FALSE}
basho <- getwd()
# Forward側
fnFs <- list.files(pattern = "_1.fastq", path = basho, full.names = TRUE)
hoge <- strsplit(basename(fnFs), "_")
sampleID <- sapply(X = hoge, FUN = `[`, 1)
filtFs <- file.path(basho, "processed", paste0(sampleID, "_F_filt.fastq"))
names(filtFs) <- sampleID
filtFs
# Reverse側
fnRs <- list.files(pattern = "_2.fastq", path = basho, full.names = TRUE)
hoge <- strsplit(basename(fnRs), "_")
sampleID <- sapply(X = hoge, FUN = `[`, 1)
filtRs <- file.path(basho, "processed", paste0(sampleID, "_R_filt.fastq"))
names(filtRs) <- sampleID
filtRs
```

### W6.2 相対パスの場合
```{r, message=FALSE, error=TRUE, warning=FALSE}
basho <- "."
# Forward側
fnFs <- list.files(pattern = "_1.fastq", path = basho, full.names = TRUE)
hoge <- strsplit(basename(fnFs), "_")
sampleID <- sapply(X = hoge, FUN = `[`, 1)
filtFs <- file.path(basho, "processed", paste0(sampleID, "_F_filt.fastq"))
names(filtFs) <- sampleID
filtFs
# Reverse側
fnRs <- list.files(pattern = "_2.fastq", path = basho, full.names = TRUE)
hoge <- strsplit(basename(fnRs), "_")
sampleID <- sapply(X = hoge, FUN = `[`, 1)
filtRs <- file.path(basho, "processed", paste0(sampleID, "_R_filt.fastq"))
names(filtRs) <- sampleID
filtRs
```

