第24回のウェブ資料

はじめに

• 乳酸菌NGS連載第22回のPDF
• dada2：Callahan et al., Nat Methods, 2016
  
• SRP125723 (Hallmaier-Wacker et al., Sci Rep., 2018)
  今回用いるリードデータです。オリジナルは計64サンプル分のデータを取得していますが、以下の12サンプル分のデータのみ解析します。

サンプルID	DNA抽出キット	保存バッファ	リード数
SRR6325813	QMINI	Native	282517
SRR6325815	QMINI	LysisBuffer	266259
SRR6325816	QMINI	RNA-later	359984
SRR6325817	QMINI	PBS	325660
SRR6325828	MOBIO	RNA-later	206939
SRR6325829	MOBIO	PBS	261767
SRR6325830	MOBIO	Native	238627
SRR6325831	MOBIO	LysisBuffer	225548
SRR6325859	GENIAL	LysisBuffer	290383
SRR6325860	GENIAL	Native	253908
SRR6325861	GENIAL	PBS	120789
SRR6325862	GENIAL	RNA-later	258792

W1 クオリティコントロール

W1.1 事前準備

　第23回のW6.2と基本的に同じスクリプトです。fnFsとfnRsがQC前のファイルの相対パス情報、filtFsとfiltRsがQC後のファイルの相対パス情報になります。

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ベクトルの要素名にする

• 乳酸菌NGS連載第4回のPDF

W1.2 plotQualityProfile関数

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

W1.2.1 QC前のForward側

　まずはfnFsオブジェクトの中身を表示させています。これは、QC前のForward側のfastqファイルの相対パス情報になります。plotQualityProfile関数の入力として、fnFsオブジェクトがそのまま利用されていることがわかります。入出力の関係からも明らかですが、入力は文字列ベクトルですので、ベクトルの要素数がファイル数に相当します。これが図1と同じものになります。

fnFs                                      # fnFsの中身を確認

##  [1] "./SRR6325813_1.fastq" "./SRR6325815_1.fastq" "./SRR6325816_1.fastq"
##  [4] "./SRR6325817_1.fastq" "./SRR6325828_1.fastq" "./SRR6325829_1.fastq"
##  [7] "./SRR6325830_1.fastq" "./SRR6325831_1.fastq" "./SRR6325859_1.fastq"
## [10] "./SRR6325860_1.fastq" "./SRR6325861_1.fastq" "./SRR6325862_1.fastq"

plotQualityProfile(fnFs)                  # fnFsのクオリティスコアのプロット

　入力が12個のファイル名情報ですので、出力として12個のプロットが得られます。各プロットの横軸はリード中の塩基のポジション（5’端が左、3’端が右）、縦軸はクオリティスコアです。連載第4回でも解説していますが、高ければ高いほどよいと解釈します。

W1.2.2 QC前のReverse側

　次に、fnRsオブジェクトの中身を表示させています。これは、QC前のReverse側のfastqファイルの相対パス情報になります。plotQualityProfile関数の入力として、fnRsオブジェクトをそのまま与えています。

fnRs                                      # fnRsの中身を確認

##  [1] "./SRR6325813_2.fastq" "./SRR6325815_2.fastq" "./SRR6325816_2.fastq"
##  [4] "./SRR6325817_2.fastq" "./SRR6325828_2.fastq" "./SRR6325829_2.fastq"
##  [7] "./SRR6325830_2.fastq" "./SRR6325831_2.fastq" "./SRR6325859_2.fastq"
## [10] "./SRR6325860_2.fastq" "./SRR6325861_2.fastq" "./SRR6325862_2.fastq"

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’」のようなメッセージが出ても、並列計算ができないだけで計算自体は動いているので問題ありません。

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）。閾値が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オブジェクトをそのまま与えています。

filtFs                                    # filtFsの中身を確認

##                            SRR6325813                            SRR6325815 
## "./processed/SRR6325813_F_filt.fastq" "./processed/SRR6325815_F_filt.fastq" 
##                            SRR6325816                            SRR6325817 
## "./processed/SRR6325816_F_filt.fastq" "./processed/SRR6325817_F_filt.fastq" 
##                            SRR6325828                            SRR6325829 
## "./processed/SRR6325828_F_filt.fastq" "./processed/SRR6325829_F_filt.fastq" 
##                            SRR6325830                            SRR6325831 
## "./processed/SRR6325830_F_filt.fastq" "./processed/SRR6325831_F_filt.fastq" 
##                            SRR6325859                            SRR6325860 
## "./processed/SRR6325859_F_filt.fastq" "./processed/SRR6325860_F_filt.fastq" 
##                            SRR6325861                            SRR6325862 
## "./processed/SRR6325861_F_filt.fastq" "./processed/SRR6325862_F_filt.fastq"

plotQualityProfile(filtFs)                # filtFsクオリティスコアのプロット

W1.4.2 QC後のReverse側

　次に、まずfiltRsオブジェクトの中身を表示させています。これは、QC後のReverse側のfastqファイルの相対パス情報になります。plotQualityProfile関数の入力として、filtRsオブジェクトをそのまま与えています。

filtRs                                    # filtRsの中身を確認

##                            SRR6325813                            SRR6325815 
## "./processed/SRR6325813_R_filt.fastq" "./processed/SRR6325815_R_filt.fastq" 
##                            SRR6325816                            SRR6325817 
## "./processed/SRR6325816_R_filt.fastq" "./processed/SRR6325817_R_filt.fastq" 
##                            SRR6325828                            SRR6325829 
## "./processed/SRR6325828_R_filt.fastq" "./processed/SRR6325829_R_filt.fastq" 
##                            SRR6325830                            SRR6325831 
## "./processed/SRR6325830_R_filt.fastq" "./processed/SRR6325831_R_filt.fastq" 
##                            SRR6325859                            SRR6325860 
## "./processed/SRR6325859_R_filt.fastq" "./processed/SRR6325860_R_filt.fastq" 
##                            SRR6325861                            SRR6325862 
## "./processed/SRR6325861_R_filt.fastq" "./processed/SRR6325862_R_filt.fastq"

plotQualityProfile(filtRs)                # filtRsのクオリティスコアのプロット

W1.5 QC前後のリード数を確認

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

W1.5.1 基本形

　outの中身を表示させているだけです。reads.in列がQC前、reads.out列がQC後のリード数です。

out                                       # outの中身を確認

##                    reads.in reads.out
## SRR6325813_1.fastq   282517    274412
## SRR6325815_1.fastq   266259    259243
## SRR6325816_1.fastq   359984    346364
## SRR6325817_1.fastq   325660    312068
## SRR6325828_1.fastq   206939    196357
## SRR6325829_1.fastq   261767    247957
## SRR6325830_1.fastq   238627    229740
## SRR6325831_1.fastq   225548    218387
## SRR6325859_1.fastq   290383    282994
## SRR6325860_1.fastq   253908    247186
## SRR6325861_1.fastq   120789    115542
## SRR6325862_1.fastq   258792    248696

W1.5.2 割合の情報を追加

　先ほどの結果より、outの行列は1列目がQC前、2列目がQC後のリード数であることはすぐにわかります。ここではまず、2列目/1列目として残ったリードの割合をratioというオブジェクトにまず格納しています。次に、元々2列の情報からなるout行列の右側に、先ほど得たratioオブジェクトをcbind関数を用いて列方向で結合して出力しています。ここで得たものが、表1bの割合の情報になります。

ratio <- out[,2]/out[,1]                  # outの2列目/1列目を計算し、ratioに格納
cbind(out, ratio)                         # outの右側にratioの情報を列方向で連結

##                    reads.in reads.out     ratio
## SRR6325813_1.fastq   282517    274412 0.9713115
## SRR6325815_1.fastq   266259    259243 0.9736497
## SRR6325816_1.fastq   359984    346364 0.9621650
## SRR6325817_1.fastq   325660    312068 0.9582632
## SRR6325828_1.fastq   206939    196357 0.9488642
## SRR6325829_1.fastq   261767    247957 0.9472432
## SRR6325830_1.fastq   238627    229740 0.9627578
## SRR6325831_1.fastq   225548    218387 0.9682507
## SRR6325859_1.fastq   290383    282994 0.9745543
## SRR6325860_1.fastq   253908    247186 0.9735258
## SRR6325861_1.fastq   120789    115542 0.9565606
## SRR6325862_1.fastq   258792    248696 0.9609880

　トリムやフィルタリングで残ったリードの割合の最大値と最小値は、それぞれ以下の通りです：

• 最大値とそのサンプル名：0.9745543、SRR6325859_1.fastq
  
• 最小値とそのサンプル名：0.9472432、SRR6325829_1.fastq

W1.6 ユニーク配列の数

ユニーク配列（unique sequences）とは、長さはもちろんのこと、同一ポジションの塩基が完全に一致する配列のことです。ここでは、dada2パッケージが提供するgetUniques関数を用いてユニーク配列の数をカウントする例を示します。ここで得たものが、表1のユニーク配列の情報になります。

W1.6.1 基本的な考え方

　いきなり250 bpの長いリードで結果を示しても理解しづらいので、ここではまず第23回の図2bの左下に見えている7本のリード配列をサンプルデータとして用います。そして、実質的にgetUniques関数と同じ挙動を示すtable関数を用いて感覚をつかみます。

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
## ACCGAGGAGGTCGGCGCCTC ACGGAGGACGTCGGCGCCTA ACGGAGGAGGTCGGCGCCTC 
##                    1                    1                    3 
## ACGGATGAGGTCGGCGACTC ATGGAGGAGGTCGGCGCCTC 
##                    1                    1

　hogeは文字列ベクトルであり、計7つのリード塩基配列が格納されていることがわかります。つまり、リード数は7です。そして、よく眺めると、2, 4, 5番目のリードのみが同一であり、それ以外はどこかのポジションの塩基が異なっていることがわかります。つまり、ユニーク配列の数は5です。この事実を認識したうえでtable(hoge)実行結果を眺めると、table関数の挙動がよくわかると思います。つまり、table関数にベクトルを入力として与えると、ユニークな要素を同定するとともに、要素ごとの出現頻度を返してくれるのです。

　この場合のtable(hoge)実行結果の要素数がユニークな配列数となりますので、ベクトルの要素数を調べるlength関数を引き続いて実行す れば、「ユニークな配列数は5個」という結果が得られるのです。

length(table(hoge))                       # ベクトル中のユニークな要素数を返す

## [1] 5

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

sum(table(hoge))                          # 出現頻度の総和（つまりリード数）を計算

## [1] 7

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

sort(table(hoge), decreasing = T)         # ベクトル中の要素ごとの出現頻度を降順で返す

## hoge
## ACGGAGGAGGTCGGCGCCTC ACCGAGGAGGTCGGCGCCTC ACGGAGGACGTCGGCGCCTA 
##                    3                    1                    1 
## ACGGATGAGGTCGGCGACTC ATGGAGGAGGTCGGCGCCTC 
##                    1                    1

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

W1.6.2 getUniques関数

　次に、getUniques関数の利用例を示します。この関数は、fastqファイルを読み込んで、ユニーク配列ごとの出現頻度を降順（多 \(\rightarrow\) 少）で返してくれます。ここではまず、「QC前のForward側のfastqファイル群の名前」情報を格納したfnFsオブジェクトの1番目の要素（つまり”./SRR6325813_1.fastq”）のみを入力として与えて、getUniques関数実行結果として得られるユニーク配列ごとの出現頻度情報に、さらにsum関数を実行しています。つまり、getUniques関数を利用してリード数情報を得ようとしています。

hoge <- fnFs[1]                           # fnFsの1番目の要素をhogeに格納
hoge                                      # 中身を確認

## [1] "./SRR6325813_1.fastq"

sum(getUniques(hoge))                     # リード数情報

## [1] 282517

　得られた結果（つまり282517個）は、確かに第23回の表2で示されているリード数と一致しています。次に、出現頻度の多い上位3つをhead関数を用いて示します。1つ1つのリードが250 bp（250文字）あるので見づらいかと思いますが、1番多いのは60791個、2番目に多いのは36760個と解釈します。

head(getUniques(hoge), n = 3)             # hoge中のユニーク配列の最初の3つ（塩基配列と出現頻度）を確認

## TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCCGGGAACTGCGTTCTGAACTGGGTGACTCGAGTGTGTCAGAGGGAGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCTCCTGGGACAACACTGACGTTCATGCCCGAAAGCGTGGGTAGCAAACA 
##                                                                                                                                                                                                                                                       60791 
## TACGAAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCAGCAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCAAAACTACTGAGCTAGAGTACGGTAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGACTGATACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACA 
##                                                                                                                                                                                                                                                       36760 
## TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGCAGGTGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACACTGAGGCGCGAAAGCGTGGGGAGCAAACA 
##                                                                                                                                                                                                                                                       25789

　次に、出現頻度の少ない上位3つをtail関数を用いて示します。予想通りですが、3つとも1個だと解釈します。

tail(getUniques(hoge), n = 3)             # hoge中のユニーク配列の最後の3つ（塩基配列と出現頻度）を確認

## TTCGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGCGCGCAGGCGGTCTTTTAAGTCTGATGTGGAAGCCCCCGGCTTAACCGGGGAGGGTCATTGGAAACTGGAAGACTGGAGTGCAGAAGAGGAGAGTGGAATTCCACGTGTAGCGGTGAAATGCGTAGATATGTGGAGGAACACCAGTGGCGAAGGCGACTCTCTGGTCTGTAACTGACGCTGAGGCGCGAGAGCGTGGGGAGCAAACA 
##                                                                                                                                                                                                                                                           1 
## TTCGTAGGTGGCAAGCGTTGTCCGGCTTTATTTGGCGTAAAGCGAGCGCAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCCCGGCTCAACCGGGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAGAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGATATATGGAGGAACACCAGTGGCGAAGGCGGCTCTCTGGTCTGTAACTGACGCTGAGGCTCGAAAGCGTGGGGAGCAAACA 
##                                                                                                                                                                                                                                                           1 
## TTCGTAGGTGGCAAGCGTTGTCCGGTTTTATTTGGCGTAAAGCGCGCGCAGGCGGTCTTTTAAGTCTGATGTGAAAGCCCCCGGCTTAACCGGGGAGGGTCATTGGAAACTGGAAGACTGGAGTGCAGAAGAGGAGAGTGGAATTCCACGTGTAGCGGTGAAATGCGTAGATATGTGGAGGAACACCAGTGGCGAAGGCGACTCTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA 
##                                                                                                                                                                                                                                                           1

　ユニーク配列の数は、以下のようにlength関数を用いて得ることができます。

length(getUniques(hoge))                  # ユニークな配列数

## [1] 47712

　このサンプル（./SRR6325813_1.fastq）のリード数は282517個ですので、ユニーク配列1本あたり282517 / 47712 = 5.9212986個の同一塩基配列のリードが存在すると解釈できます。

sum(getUniques(hoge)) / length(getUniques(hoge)) # ユニークな配列1本あたりの出現頻度を計算

## [1] 5.921299

W1.6.3 自作のgetUniNum関数

　様々な書き方（定義の仕方）がありますが、例えば以下のようにすれば、ファイル中のユニーク配列の数を得るgetUniNumという名前の関数を定義することができます。この関数の実体は、getUniques関数実行に引き続いてlength関数を実行することです。なお、xは入力として与える情報に相当します。

getUniNum <- function(x) length(getUniques(x))   # ユニーク配列数を得るための関数を定義

　一度定義してしまえば、あとは通常の関数と同様の利用法ができます。

hoge <- fnFs[1]                           # fastqファイル名をhogeに格納
getUniNum(hoge)                           # ユニーク配列数を確認

## [1] 47712

　getUniNum関数は、1つのfastqファイルの入力しか想定していません。そのため、下記のように「QC前のForward側のfastqファイル群の名前」情報を格納したfnFsオブジェクトをそのまま入力として与えるとエラーが出ます。エラーが出ることが分かっている状態でエラーメッセージを眺めると、確かに「…Unrecognized format…」などと書かれている意味がなんとなく分かると思います。

fnFs                                      # 中身を確認

##  [1] "./SRR6325813_1.fastq" "./SRR6325815_1.fastq" "./SRR6325816_1.fastq"
##  [4] "./SRR6325817_1.fastq" "./SRR6325828_1.fastq" "./SRR6325829_1.fastq"
##  [7] "./SRR6325830_1.fastq" "./SRR6325831_1.fastq" "./SRR6325859_1.fastq"
## [10] "./SRR6325860_1.fastq" "./SRR6325861_1.fastq" "./SRR6325862_1.fastq"

getUniNum(fnFs)                           # ユニーク配列数を確認

## Error in getUniques(x): Unrecognized format: Requires named integer vector, fastq filename, dada-class, derep-class, sequence matrix, or a data.frame with $sequence and $abundance columns.

W1.6.4 lapply関数を利用

　12個のfastqファイルそれぞれのリード数を一度に実行する1つのやり方は、lapply関数の利用です。これは、fnFsの要素ごとにgetUniNum関数を実行してくれます。実行結果は[[1]]のようなリスト形式で返されますが、表示されているユニーク配列数は正しいです。

lapply(fnFs, getUniNum)                   # fnFs中のサンプルごとのユニーク配列数をリスト形式で確認

## [[1]]
## [1] 47712
## 
## [[2]]
## [1] 44341
## 
## [[3]]
## [1] 61184
## 
## [[4]]
## [1] 60466
## 
## [[5]]
## [1] 40611
## 
## [[6]]
## [1] 50769
## 
## [[7]]
## [1] 41645
## 
## [[8]]
## [1] 37710
## 
## [[9]]
## [1] 49343
## 
## [[10]]
## [1] 43452
## 
## [[11]]
## [1] 24214
## 
## [[12]]
## [1] 44494

W1.6.5 unlist関数も併用

　結果をリスト形式ではなくベクトル形式で得たい場合は、以下で示すようにunlist関数を追加で実行すればよいです。

unlist(lapply(fnFs, getUniNum))           # fnFs中のサンプルごとのユニーク配列数をベクトル形式で確認

##  [1] 47712 44341 61184 60466 40611 50769 41645 37710 49343 43452 24214 44494

W1.6.6 sapply関数を利用

　以下で示すようにsapply関数を利用してもかまいません。このあたりは好みですので、多様なやり方を知り、見慣れない関数や出力結果に狼狽しないことのみが重要です。

sapply(fnFs, getUniNum)                   # fnFs中のサンプルごとのユニーク配列数を確認

## ./SRR6325813_1.fastq ./SRR6325815_1.fastq ./SRR6325816_1.fastq 
##                47712                44341                61184 
## ./SRR6325817_1.fastq ./SRR6325828_1.fastq ./SRR6325829_1.fastq 
##                60466                40611                50769 
## ./SRR6325830_1.fastq ./SRR6325831_1.fastq ./SRR6325859_1.fastq 
##                41645                37710                49343 
## ./SRR6325860_1.fastq ./SRR6325861_1.fastq ./SRR6325862_1.fastq 
##                43452                24214                44494

W1.6.7 QC前のForward側

　実質的に直前の項と同じ内容ですが、sapply実行結果をuniseq_fnFsに格納しています。これが「表1、(a)QC前、ユニーク配列の数、Forward」の基礎情報です。

uniseq_fnFs <- sapply(fnFs, getUniNum)    # fnFs中のサンプルごとのユニーク配列数をuniseq_fnFsに格納
uniseq_fnFs                               # 中身を確認

## ./SRR6325813_1.fastq ./SRR6325815_1.fastq ./SRR6325816_1.fastq 
##                47712                44341                61184 
## ./SRR6325817_1.fastq ./SRR6325828_1.fastq ./SRR6325829_1.fastq 
##                60466                40611                50769 
## ./SRR6325830_1.fastq ./SRR6325831_1.fastq ./SRR6325859_1.fastq 
##                41645                37710                49343 
## ./SRR6325860_1.fastq ./SRR6325861_1.fastq ./SRR6325862_1.fastq 
##                43452                24214                44494

W1.6.8 QC前のReverse側

　「表1、(a)QC前、ユニーク配列の数、Reverse」の基礎情報です。

uniseq_fnRs <- sapply(fnRs, getUniNum)    # fnRs中のサンプルごとのユニーク配列数をuniseq_fnRsに格納
uniseq_fnRs                               # 中身を確認

## ./SRR6325813_2.fastq ./SRR6325815_2.fastq ./SRR6325816_2.fastq 
##                97726                90228               115402 
## ./SRR6325817_2.fastq ./SRR6325828_2.fastq ./SRR6325829_2.fastq 
##               110355                65478                81542 
## ./SRR6325830_2.fastq ./SRR6325831_2.fastq ./SRR6325859_2.fastq 
##                75736                70284               107800 
## ./SRR6325860_2.fastq ./SRR6325861_2.fastq ./SRR6325862_2.fastq 
##                95390                45201                89439

W1.6.9 QC後のForward側

　「表1、(b)QC後、ユニーク配列の数、Forward」の基礎情報です。

uniseq_filtFs <- sapply(filtFs, getUniNum)# filtFs中のサンプルごとのユニーク配列数をuniseq_filtFsに格納
uniseq_filtFs                             # 中身を確認

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##      26947      25669      31854      31377      20033      23878      21710 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##      20537      28213      24314      12800      22918

W1.6.10 QC後のReverse側

　「表1、(b)QC後、ユニーク配列の数、Reverse」の基礎情報です。

uniseq_filtRs <- sapply(filtRs, getUniNum)# filtRs中のサンプルごとのユニーク配列数をuniseq_filtRsに格納
uniseq_filtRs                             # 中身を確認

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##      34499      33033      36310      34583      18587      20516      23508 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##      23016      39027      33721      14794      28879

W2 ノイズ除去

　dada2の主要部分は、ノイズ除去と重複除去（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オブジェクトに格納しています。

errFs <- learnErrors(fls = filtFs, multithread = TRUE) # filtFsの一部を用いたエラーモデル推定

## 109399275 total bases in 533655 reads from 2 samples will be used for learning the error rates.

errRs <- learnErrors(fls = filtRs, multithread = TRUE) # filtRsの一部を用いたエラーモデル推定

## 123202660 total bases in 880019 reads from 3 samples will be used for learning the error rates.

W2.2 plotErrors関数

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

W2.2.1 Forward側

　Forward側のlearnErrors関数の実行結果であるerrFsを入力としてプロットしています。nominalQ = TRUEは、赤線をつけよという指令です。これが図2と同じものになります。

plotErrors(errFs, nominalQ = TRUE)        # errFsの中身（エラー頻度やエラーモデル推定結果）をプロット

W2.2.2 Reverse側

　Reverse側のlearnErrors関数の実行結果であるerrRsを入力としてプロットしています。

plotErrors(errRs, nominalQ = TRUE)        # errRsの中身（エラー頻度やエラーモデル推定結果）をプロット

W2.3 dada関数

　dada関数は、前処理が終わったリードファイル（Forward側の場合はfiltFs）と先ほど得たエラーモデル情報（Forward側の場合はerrFs）を入力として、denoisingを実行します。

dadaFs <- dada(derep = filtFs, err = errFs, multithread = TRUE) # filtFsとerrFsを入力としてエラー補正を実行

## Sample 1 - 274412 reads in 26947 unique sequences.
## Sample 2 - 259243 reads in 25669 unique sequences.
## Sample 3 - 346364 reads in 31854 unique sequences.
## Sample 4 - 312068 reads in 31377 unique sequences.
## Sample 5 - 196357 reads in 20033 unique sequences.
## Sample 6 - 247957 reads in 23878 unique sequences.
## Sample 7 - 229740 reads in 21710 unique sequences.
## Sample 8 - 218387 reads in 20537 unique sequences.
## Sample 9 - 282994 reads in 28213 unique sequences.
## Sample 10 - 247186 reads in 24314 unique sequences.
## Sample 11 - 115542 reads in 12800 unique sequences.
## Sample 12 - 248696 reads in 22918 unique sequences.

dadaRs <- dada(derep = filtRs, err = errRs, multithread = TRUE) # filtRsとerrRsを入力としてエラー補正を実行

## Sample 1 - 274412 reads in 34499 unique sequences.
## Sample 2 - 259243 reads in 33033 unique sequences.
## Sample 3 - 346364 reads in 36310 unique sequences.
## Sample 4 - 312068 reads in 34583 unique sequences.
## Sample 5 - 196357 reads in 18587 unique sequences.
## Sample 6 - 247957 reads in 20516 unique sequences.
## Sample 7 - 229740 reads in 23508 unique sequences.
## Sample 8 - 218387 reads in 23016 unique sequences.
## Sample 9 - 282994 reads in 39027 unique sequences.
## Sample 10 - 247186 reads in 33721 unique sequences.
## Sample 11 - 115542 reads in 14794 unique sequences.
## Sample 12 - 248696 reads in 28879 unique sequences.

W2.4 結果の確認

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

head(dadaFs, n = 2)                       # dadaFs中の最初の2つを確認

## $SRR6325813
## dada-class: object describing DADA2 denoising results
## 270 sequence variants were inferred from 26947 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16
## 
## $SRR6325815
## dada-class: object describing DADA2 denoising results
## 328 sequence variants were inferred from 25669 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16

head(dadaRs, n = 2)                       # dadaRs中の最初の2つを確認

## $SRR6325813
## dada-class: object describing DADA2 denoising results
## 190 sequence variants were inferred from 34499 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16
## 
## $SRR6325815
## dada-class: object describing DADA2 denoising results
## 224 sequence variants were inferred from 33033 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16

W2.5 ユニーク配列数を確認

　W1.6.3で定義した自作関数getUniNumをここで改めて再定義して、dada関数実行後にユニーク配列の数がどれだけ減っているかを確認します。ここの数値が表2のユニーク配列数の基礎情報となります。

getUniNum <- function(x) length(getUniques(x))# ユニーク配列数を得るための関数を定義
sapply(dadaFs, getUniNum)                 # dadaFs中のサンプルごとのユニーク配列数を確認

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##        270        328        369        317        205        225        198 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##        185        404        298         83        155

sapply(dadaRs, getUniNum)                 # dadaRs中のサンプルごとのユニーク配列数を確認

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##        190        224        273        239        159        171        141 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##        131        262        192         56        109

W2.6 リード数を確認

　W1.6.3の自作関数getUniNumはユニーク配列数を得るためのものでしたが、length関数部分をsum関数に置き換えることで、リード数を得るための関数に書き替えることができます。ここでは、それをgetReadNumという自作関数として定義しています。ここの数値が表2のリード数の基礎情報となります。

getReadNum <- function(x) sum(getUniques(x))# リード数を得るための関数を定義
sapply(dadaFs, getReadNum)                # dadaFs中のサンプルごとのリード数を確認

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##     273973     258614     345867     311660     196147     247769     229465 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##     218210     282429     246714     115462     248439

sapply(dadaRs, getReadNum)                # dadaRs中のサンプルごとのリード数を確認

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##     272894     257784     344679     310610     195695     247249     228684 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##     217404     281362     245839     115080     247742

W2.7 ユニーク配列の出現頻度

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

W2.7.1 出現頻度の上位

　まずは1番目のサンプル（サンプルID：SRR6325813）に対して、上位4つを眺めます。1つ1つのリードが長いので見づらいかと思いますが、1番多いのは82848個、2番目に多いのは52232個と解釈します。

head(getUniques(dadaFs[[1]]), n = 4)      # dadaFs中の1番目の要素に対して、ユニーク配列の最初の4つ（塩基配列と出現頻度）を確認

## ATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCCGGGAACTGCGTTCTGAACTGGGTGACTCGAGTGTGTCAGAGGGAGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCTCCTGGGACAACACTGACGTTC 
##                                                                                                                                                                                                         82848 
## ATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCAGCAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCAAAACTACTGAGCTAGAGTACGGTAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGACTGATACTGACACTG 
##                                                                                                                                                                                                         52232 
## ATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTC 
##                                                                                                                                                                                                         38747 
## TCCGGAATTATTGGGCGTAAAGCGCGCGCAGGTGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACACTG 
##                                                                                                                                                                                                         35026

W2.7.2 出現頻度の下位

　次に、下位4つを確認します。さきほどのhead関数の結果から、getUniques関数実行結果は出現頻度でソートされているものだと思っていましたが、そうはなっていないことがわかります。

tail(getUniques(dadaFs[[1]]), n = 4)      # dadaFs中の1番目の要素に対して、ユニーク配列の最後の4つ（塩基配列と出現頻度）を確認

## TCCGGATTTATTGGGTTTAAAGGGTGCGTAGGTTGTTCGGTAAGTCAGCGGTGAAACCTGAGCGCTCAACGTTCAGCCTGCCGTTGAAACTGCCGGGCTTGAGTTCAGCGGCGGCAGGCGGAATTCGTGGTGTAGCGGTGAAATGCATAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGACTGATACTGACACTG 
##                                                                                                                                                                                                             4 
## TCCGGATTTATTGGGCGTAAAGCGAGCGCAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCCCGGCTCAACCGGGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAGAGTGGAATTCCATGTGTAGCGGTGAAATGCGCAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACGCTG 
##                                                                                                                                                                                                            13 
## TCCGGAATTATTGGGCGTAAAGCGCGCGTAGGCGGTTTTTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGAAAACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGACTGATACTGACACTG 
##                                                                                                                                                                                                             8 
## ATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACACTG 
##                                                                                                                                                                                                            15

　このような場合は、例えばW1.6.1などを参考にして、自分で降順（多 \(\rightarrow\) 少）にソートしてから、改めてtail関数を実行するなどして対処します。W1.6.2の時点では、このサンプル（dadaFs[[1]]）は出現頻度が1回のものが存在しましたが、dada関数実行後はなくなっていることがわかります。

hoge <- sort(getUniques(dadaFs[[1]]), decreasing = T) # dadaFs中の1番目の要素に対して、ユニーク配列の出現頻度が多い順にソートした結果をhogeに格納
tail(hoge, n = 4)                         # hoge中の最後の4つのを確認

## ATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCAGCAAGTTGGATGTGAAAGCCTTCGGCTCAACCGGAGAATTGCATCAGAAACTGTTGAACTTGAGTGCAGAAGAGGAGAGTGGAACTCCATGTGTAGCGGTGGAATGCGTAGATATATGGAAGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACACTG 
##                                                                                                                                                                                                             2 
## ATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCTGGTGGAGCGGTGGAATGCGCAGATATCAGGAGGAACACCAGTGGCGAAGGCGGTTCTCTGGGCCTTTCCTGACGCTG 
##                                                                                                                                                                                                             2 
## TCCGGATTTATTGGGCGTAAAGCGAGCGCAGACGGTTATTTAAGTCTGAAGTGAAAGCCCTCAGCTCAACTGAGGAATTGCTTTGGAAACTGGATGACTTGAGTGCAGTAGAGGAAAGTGGAACTCCATGTGTAGCGGTGAAATGCGTAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACACTG 
##                                                                                                                                                                                                             2 
## TTCGGATTCACTGGGCGTAAAGGGTGCGTAGGCGGTCGGGTAAGTCTGACGTGAAATCTCCGGGCCTAACCCGGAAACTGCGTCGGATACTATCCGGCTAGAGGATTGGAAGGGAGACTGGAATACTTGGTGTAGCAGTGAAATGCGTAGATATCAAGTGGAACACCAGTGGCGAAGGCGAGTCTCTGGACAATTCCTGACGCTG 
##                                                                                                                                                                                                             2

　ここで知りたいことは、主に出現頻度が最小のものが1ではなくなっているかどうかです。それゆえ、例えば以下のような感じで、min関数を用いて出現頻度の最小値のみを得るためのgetMinNum関数を自作します。W1.6.3、W2.5、およびW2.6を参考にしています。

getMinNum <- function(x) min(getUniques(x))# 出現頻度の最小値を得るための関数を定義

　あとは、W2.5やW2.6と同じようにしてdada関数実行結果（dadaFsとdadaRs）を入力として実行すればよいです。Forward側Reverse側含め全サンプルで2になっていることから、確かにdada2は出現頻度が2以上あるリードのみに限定していることがわかります。

sapply(dadaFs, getMinNum)                 # dadaFs中のサンプルごとのリード数の最小値を取得

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##          2          2          2          2          2          2          2 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##          2          2          2          2          2

sapply(dadaRs, getMinNum)                 # dadaRs中のサンプルごとのリード数の最小値を取得

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##          2          2          2          2          2          2          2 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##          2          2          2          2          2

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）という条件でマージが実行されます。

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                          # 結果を冗長に表示する（要は画面にも出力せよという指令）
)

## 268804 paired-reads (in 236 unique pairings) successfully merged out of 272596 (in 905 pairings) input.

## 253588 paired-reads (in 282 unique pairings) successfully merged out of 257333 (in 1053 pairings) input.

## 341395 paired-reads (in 299 unique pairings) successfully merged out of 344403 (in 772 pairings) input.

## 306377 paired-reads (in 249 unique pairings) successfully merged out of 310392 (in 708 pairings) input.

## 194243 paired-reads (in 182 unique pairings) successfully merged out of 195517 (in 486 pairings) input.

## 245046 paired-reads (in 206 unique pairings) successfully merged out of 247099 (in 542 pairings) input.

## 226534 paired-reads (in 168 unique pairings) successfully merged out of 228465 (in 594 pairings) input.

## 215157 paired-reads (in 171 unique pairings) successfully merged out of 217291 (in 651 pairings) input.

## 276081 paired-reads (in 371 unique pairings) successfully merged out of 280943 (in 1466 pairings) input.

## 241843 paired-reads (in 273 unique pairings) successfully merged out of 245487 (in 1030 pairings) input.

## 114574 paired-reads (in 72 unique pairings) successfully merged out of 115021 (in 190 pairings) input.

## 244352 paired-reads (in 136 unique pairings) successfully merged out of 247530 (in 464 pairings) input.

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列が一致塩基数です。

head(mergers[[1]])                        # mergers中の1番目の要素(サンプル)について、最初の6つを確認

##                                                                                                                                                                                                                sequence
## 1 ATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCCGGGAACTGCGTTCTGAACTGGGTGACTCGAGTGTGTCAGAGGGAGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCTCCTGGGACAACACTGACGTTCATGCCCGA
## 2 ATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCAGCAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCAAAACTACTGAGCTAGAGTACGGTAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGACTGATACTGACACTGAGGTGCGA
## 3 ATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGA
## 4 TCCGGAATTATTGGGCGTAAAGCGCGCGCAGGTGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACACTGAGGCGCGA
## 5 TCCGGATTTATTGGGTTTAAAGGGTGCGTAGGTTGTTCGGTAAGTCAGCGGTGAAACCTGAGCGCTCAACGTTCAGCCTGCCGTTGAAACTGCCGGGCTTGAGTTCAGCGGCGGCAGGCGGAATTCGTGGTGTAGCGGTGAAATGCATAGATATCACGAGGAACTCCGATTGCGAAGGCAGCTTGCCATACTGCGACTGACACTGAAGCACGA
## 6 ATCGGATTTACTGGGCGTAAAGCGTGCGTAGGCGGCTTATTAAGTCGGATGTGAAATCCCCGAGCTTAACTTGGGAATTGCATTCGATACTGGTGAGCTAGAGTATGGGAGAGGATGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGATGGCGAAGGCAGCCATCTGGCCTAATACTGACGCTGAGGTACGA
##   abundance forward reverse nmatch nmismatch nindel prefer accept
## 1     81414       1       1    132         0      0      1   TRUE
## 2     51794       2       2    132         0      0      1   TRUE
## 3     38101       3       3    132         0      0      1   TRUE
## 4     34709       4       4    132         0      0      1   TRUE
## 5     13645       5       6    132         0      0      1   TRUE
## 6     10653       6       5    132         0      0      1   TRUE

W3.2.2 手作業でマージ

　どのようなマージが実際に行われているのかを確認すべく、上記の結果でabundanceが13645個の5番目の配列を例として示します。まず、入力配列のForward側はdadaFsオブジェクト中の5番目の配列、Reverse側はdadaRsオブジェクト中の6番目の配列だと解釈します。つまり、以下の2つがマージ前の配列です。

dadaFs[[1]]$sequence[5]                   # dadaFs中の1番目のサンプル中の5番目の塩基配列を確認

## [1] "TCCGGATTTATTGGGTTTAAAGGGTGCGTAGGTTGTTCGGTAAGTCAGCGGTGAAACCTGAGCGCTCAACGTTCAGCCTGCCGTTGAAACTGCCGGGCTTGAGTTCAGCGGCGGCAGGCGGAATTCGTGGTGTAGCGGTGAAATGCATAGATATCACGAGGAACTCCGATTGCGAAGGCAGCTTGCCATACTGCGACTGACACTG"

dadaRs[[1]]$sequence[6]                   # dadaRs中の1番目のサンプル中の6番目の塩基配列を確認

## [1] "TCGTGCTTCAGTGTCAGTCGCAGTATGGCAAGCTGCCTTCGCAATCGGAGTTCCTCGTGATATCTATGCATTTCACCGCTACACCACGAATTCCGCCTGCCGCCGCTGAACTCAAGCCCGGCAGTTTCAACGGCAGGCTG"

　これらの配列長は、W1.3でも指定していますが、Forward側が205 bp、Reverse側が140 bpです。W3.2のnmatch列の5番目を眺めると、オーバーラップ（一致）塩基数が132個だということが分かります。このことから、ペア―ドエンドのマージによって、205 bpあるForward側のリードの右側に、Reverse側の5’端の(140 - 132) = 8 bp分だけ追加されたのだろうと解釈します。

nchar(dadaFs[[1]]$sequence[5])            # dadaFs中の1番目のサンプル中の5番目の塩基配列の文字数(配列長)を確認

## [1] 205

nchar(dadaRs[[1]]$sequence[6])            # dadaRs中の1番目のサンプル中の6番目の塩基配列の文字数(配列長)を確認

## [1] 140

　ただし、実際にマージされるのは、Reverse側の配列の逆相補鎖（reverse complement）の3’端側の8 bp分だという点に注意が必要です。rc関数を利用すれば、逆相補鎖を得ることができます。

rc(dadaRs[[1]]$sequence[6])               # dadaRs中の1番目のサンプル中の6番目の塩基配列の逆相補鎖を確認

## [1] "CAGCCTGCCGTTGAAACTGCCGGGCTTGAGTTCAGCGGCGGCAGGCGGAATTCGTGGTGTAGCGGTGAAATGCATAGATATCACGAGGAACTCCGATTGCGAAGGCAGCTTGCCATACTGCGACTGACACTGAAGCACGA"

W3.2.3 リード数を確認

　W2.6で自作したリード数を得るための関数getReadNumをここで改めて定義し、W3.1で得たmergePairs関数実行結果であるmergersオブジェクトを入力として実行します。ここの数値が表3のリード数の基礎情報となります。

getReadNum <- function(x) sum(getUniques(x))# リード数を得るための関数を定義
sapply(mergers, getReadNum)               # mergers中のサンプルごとのリード数を確認

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##     268804     253588     341395     306377     194243     245046     226534 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##     215157     276081     241843     114574     244352

W3.2.4 ユニーク配列数を確認

　W2.5で自作したリード数を得るための関数getUniNumをここで改めて定義し、W3.1で得たmergePairs関数実行結果であるmergersオブジェクトを入力として実行します。ここの数値が表3のユニーク配列数の基礎情報となります。

getUniNum <- function(x) length(getUniques(x))# ユニーク配列数を得るための関数を定義
sapply(mergers, getUniNum)                    # サンプルごとのユニーク配列数

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##        236        282        299        249        182        206        168 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##        171        371        273         72        136

W3.3 サンプル全体の結果

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

W3.3.1 getSequences関数

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

hoge <- sapply(mergers, getSequences)     # mergers中のサンプルごとの塩基配列をhogeに格納

W3.3.2 class関数

　得られたhogeオブジェクトは、サンプルごとに切り分けられたリスト形式のオブジェクトであることが、それを確認するためのclass関数実行結果からわかります。つまり、表示結果が”list”になっていれば、それはリスト形式のオブジェクトだということです。

class(hoge)                               # hogeオブジェクトのクラスを確認

## [1] "list"

W3.3.3 unlistとis.vector関数

　W1.6.5でも使いましたが、リスト形式をベクトル形式に変換したい場合には、unlist関数を利用します。得られたhoge2オブジェクトがベクトルか否かをTRUE or FALSEで返すis.vector関数実行結果からも、確かにベクトル形式に変換できていることがわかります。

hoge2 <- unlist(hoge)                     # リスト形式をベクトル形式に変換
is.vector(hoge2)                          # ベクトルか否かを表示

## [1] TRUE

W3.3.4 sapplyと自作関数

　hoge2ベクトルの要素数は、W3.2.3で得たサンプルごとのユニーク配列数の総和と一致するはずです。ここでは確かにそうであることを確認しているだけです。

length(hoge2)                             # hoge2ベクトルの要素数を確認

## [1] 2645

sum(sapply(mergers, getUniNum))           # mergers中のサンプルごとのユニーク配列数の総和(つまり総リード数)を確認

## [1] 2645

W3.3.5 ユニーク配列の和集合

　サンプル全体のユニーク配列数は、以下のようにunique関数を実行した上で、length関数を実行すればよいです。この1085という数字は、ユニーク配列の和集合（union）ともいえます。

length(unique(hoge2))                     # hoge2ベクトル中のユニークな要素数を確認

## [1] 1085

　上記のやり方以外に、table関数を用いてユニーク配列ごとの出現頻度を調べると、table関数実行結果の要素数が「サンプル全体のユニーク配列数」に相当することを利用することもできます。

length(table(hoge2))                      # hoge2ベクトル中のユニークな要素数を確認

## [1] 1085

W3.3.6 ユニーク配列の積集合

　ほぼ余談ですが、ユニーク配列の和集合（union）が分かると、計12サンプル全てに存在するユニーク配列の数（積集合；intersection）も知りたくなると思います。ここではユニーク配列ごとの出現頻度を調べるtable関数実行結果から、出現頻度が12回（サンプル数の合計が12個だからです）のもののみを抽出するやり方で示します。実行結果は26個です。

sum(table(hoge2) == 12)                   # hoge2ベクトル中の出現頻度が12回のユニーク配列数を確認

## [1] 26

W3.3.7 スクリプトを一まとめ

　ユニーク配列の和集合（1085個）と積集合（26個）を得るためのスクリプトをひとまとめにしたのが以下になります。この一連の流れを理解しておくと、関数が無数に入れ子になっていても、ステップバイステップで理解していけばいいのだと安心できると思います。（本当は最終行の12という数値も汎用性の高いものにするとよいのですが、ここでは見た目の分かりやすいさを優先しています。）

hoge <- sapply(mergers, getSequences)     # mergers中のサンプルごとの塩基配列をhogeに格納
hoge2 <- unlist(hoge)                     # hogeをベクトル形式に変換してhoge2に格納
length(unique(hoge2))                     # hoge2ベクトル中のユニークな要素数(ユニーク配列の和集合)を確認

## [1] 1085

sum(table(hoge2) == 12)                   # hoge2ベクトル中の出現頻度が12回のユニーク配列数(ユニーク配列の積集合)を確認

## [1] 26

W3.4 makeSequenceTable関数

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

W3.4.1 まずは実行

　先ほど得たmergePairs関数実行結果であるmergersオブジェクトを入力として実行し、結果をseqtabオブジェクトに格納するところまで行います。

seqtab <- makeSequenceTable(mergers)      # mergers中のユニーク配列の和集合情報を取得し、サンプルごとのユニーク配列の出現頻度をseqtabに格納

W3.4.2 結果の概要を確認

　結果は行列ですので、まずはdim関数で行数と列数を表示させます。行数は12、列数は1085であることがわかります。このことから、今回解析する計12サンプル全体としては、1085個のユニークな配列が存在するのだと解釈します。W3.3.5と同じ結果が得られていることがわかります。

dim(seqtab)                               # seqtabの行数と列数を確認

## [1]   12 1085

W3.4.3 rowSumsとsum関数

　rowSums関数は、数値行列を入力として、行（row）ごとの和（sum）を返します。得られる結果はW3.2.3と同じです。同じノリで行列全体の総和をsum関数で計算すると、サンプル全体の総リード数が得られます。得られる結果（2927994個）は表2bのマージ後のリード数の総和と一致します。

rowSums(seqtab)                           # seqtabの行ごと(サンプルごと)の和を確認

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##     268804     253588     341395     306377     194243     245046     226534 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##     215157     276081     241843     114574     244352

sum(seqtab)                               # seqtab中の総配列数を確認

## [1] 2927994

W3.4.4 getSequences関数

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

hoge <- getSequences(seqtab)              # seqtab中の塩基配列をhogeに格納
length(hoge)                              # hoge中の要素数を確認

## [1] 1085

W3.4.5 nchar関数

　nchar関数は、文字列ベクトルを入力として、要素ごとの文字数を返します（W3.2.2）。したがって、この場合のhoge2オブジェクトの要素数は、入力である文字列ベクトルの要素数と同じ（1085）になります。head関数で最初の6つ分の中身を眺めると、図3で得た213 bpというコンティグ長と酷似した結果ばかり表示されていることがわかります。

hoge2 <- nchar(hoge)                      # hoge中の要素ごとの文字数をhoge2に格納
length(hoge2)                             # hoge2中の要素数を確認

## [1] 1085

head(hoge2)                               # hoge2中の最初の6つを確認

## [1] 213 213 213 213 213 213

W3.4.6 table関数

　実質的に先ほどまでの結果にtable関数を実行しているだけです。table関数は、入力ベクトル中の要素の種類ごとの出現頻度を返します。したがって、この結果は、1085個のユニークな塩基配列の長さ分布情報を得ていることと同義です。

hoge <- getSequences(seqtab)              # seqtab中の塩基配列をhogeに格納
hoge2 <- nchar(hoge)                      # hoge中の要素ごとの文字数をhoge2に格納
table(hoge2)                              # hoge2ベクトル中の要素の種類ごとに出現頻度をカウント

## hoge2
## 205 212 213 214 215 217 227 231 
##   6  44 987  42   2   2   1   1

W3.4.7 一気に書く

　一気に書くとこんな感じです。

table(nchar(getSequences(seqtab)))        # seqtab中の塩基配列ごとの文字数をカウントし、その文字数の種類ごとの出現頻度を表示

## 
## 205 212 213 214 215 217 227 231 
##   6  44 987  42   2   2   1   1

　W3.3.7に似せて、以下のような書き方でも構いません。

hoge <- sapply(mergers, getSequences)     # mergers中のサンプルごとの塩基配列をhogeに格納
hoge2 <- unlist(hoge)                     # hogeをベクトル形式に変換してhoge2に格納
hoge3 <- nchar(unique(hoge2))             # hoge2中のユニークな要素ごとの文字数をhoge3に格納
table(hoge3)                              # hoge3ベクトル中の要素の種類ごとの出現頻度を確認

## hoge3
## 205 212 213 214 215 217 227 231 
##   6  44 987  42   2   2   1   1

W3.5 removeBimeraDenovo関数

　dada2パッケージのメインであるdada関数は、エラー補正のみを行いキメラ配列の除去は行ってくれません。つまり、まだキメラ配列を含んでいます。removeBimeraDenovo関数は、makeSequenceTable関数の実行結果であるseqtabオブジェクトを入力として実行し、PCRキメラを除去した結果を返します。

W3.5.1 まずは実行

　さきほど得たmakeSequenceTable関数の実行結果であるseqtabオブジェクトを入力として、PCRキメラを除去した結果をseqtab.nochimオブジェクトに格納するところまで行います。

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関数で行数と列数を表示させます。行数は12、列数は580であることがわかります。このことから、今回解析する計12サンプル全体としては、580個のユニークな配列が存在するのだと解釈します。

dim(seqtab.nochim)                        # seqtab.nochimの行数と列数を確認

## [1]  12 580

W3.5.3 rowSumsとsum関数

　rowSums関数は、数値行列を入力として、行（row）ごとの和（sum）を返します。得られる結果はサンプルごとのリード数になります。これが同じノリで行列全体の総和をsum関数で計算すると、サンプル全体の総リード数が得られます。表2bのペア―ドエンドのマージ後のリード数からほとんど減っていないこともよくわかります。

rowSums(seqtab.nochim)                    # seqtab.nochimの行ごとの総和(除去後のサンプルごとのリード数)を表示

## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830 
##     258952     242714     336540     302984     191864     241745     223312 
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 
##     209875     253659     227111     114053     241966

sum(seqtab.nochim)                        # seqtab.nochim中の全要素の総和(除去後の総リード数)を表示  

## [1] 2844775

sum(seqtab)                               # seqtab中の全要素の総和(除去前の総リード数)を表示

## [1] 2927994

sum(seqtab.nochim)/sum(seqtab)            # 除去後の総リード数/除去前の総リード数の割合

## [1] 0.9715782

W3.5.4 配列長分布と出現頻度

　得られた計580個の配列長分布と出現頻度情報を得ます。これはW3.4.7と同じ書き方です。キメラ配列除去前と酷似した傾向であることがわかります。

table(nchar(getSequences(seqtab.nochim))) # seqtab.nochim中の塩基配列ごとの文字数をカウントし、その文字数の種類ごとの出現頻度を表示

## 
## 205 212 213 214 215 217 227 231 
##   6  34 495  39   2   2   1   1

W3.6 図4のスクリプト

　コンパクトにまとめるとこんな感じで書けます。

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)の流れに沿って説明しています。この中のAssign taxonomyという部分で、既知の16S rRNAのリファレンス配列を取得するためのページとしてhttps://benjjneb.github.io/dada2/training.htmlへのリンクが示されています。リンク先で見られるSilva 138.1 prokaryotic SSU taxonomic training data formatted for DADA2というページから取得可能な以下の2つのファイルを利用します。

• 属レベルの系統推定用：silva_nr99_v138.1_train_set.fa.gz
• 種レベルの系統推定用：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を直接指定することもできます。

#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関数で行数と列数を表示させます。行数は580、列数は6であることがわかります。次は、colnames関数でtaxa行列の列名を表示させています。この結果を眺めて、生物の階級を「界（Kingdom）・門（Phylum）・綱（Class）・目（Order）・科（Family）・属（Genus）・種（Species）」と表すと、属までに6つの階級が存在するからこのようになっているのだと納得できます。

dim(taxa)                                 # taxaの行数と列数を表示

## [1] 580   6

colnames(taxa)                            # taxaの列名を表示

## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"

W4.3 addSpecies関数

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

W4.3.1 まずは実行

　以下のような感じで実行します。W4.2.1と同様、（1番左側に#を入れてコメントアウトしているので機能はしませんが…）1行目で示すように、オリジナルのURLを直接指定することもできます。

# 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関数で行数と列数を表示させます。行数は580、列数は7であることがわかります。次は、colnames関数でtaxa2行列の列名を表示させています。生物の階級である「界（Kingdom）・門（Phylum）・綱（Class）・目（Order）・科（Family）・属（Genus）・種（Species）」の種の結果が足されていることが分かります。

dim(taxa2)                                # taxa2の行数と列数を表示

## [1] 580   7

colnames(taxa2)                           # taxa2の列名を表示

## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"

W4.4 ファイル出力

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

W4.4.1 行列の転置

　seqtab.nochimの行数と列数はそれぞれ12と580です。その一方で連結対象であるtaxa2の行数と列数は、それぞれ580と7です。従って、うまく連結すべく、seqtab.nochim行列のほうをt関数を用いて転置して、同一行が同一のASVの情報になるようにできるかどうか確認しています。

dim(seqtab.nochim)                        # seqtab.nochimの行数と列数を表示

## [1]  12 580

dim(taxa2)                                # taxa2の行数と列数を表示

## [1] 580   7

dim(t(seqtab.nochim))                     # seqtab.nochim行列を転置して、行数と列数を表示

## [1] 580  12

W4.4.2 ASV配列情報の取得

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

#ASV <- colnames(seqtab.nochim)           # seqtab.nochimの列名をASVに格納(するのでもよい)
ASV <- rownames(taxa2)                    # taxa2の行名をASVに格納
length(ASV)                               # ASVの要素数を確認

## [1] 580

head(ASV, n = 3)                          # ASV中の最初の6つを表示

## [1] "ATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCCGGGAACTGCGTTCTGAACTGGGTGACTCGAGTGTGTCAGAGGGAGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCTCCTGGGACAACACTGACGTTCATGCCCGA"
## [2] "TCCGGATTTATTGGGCGTAAAGCGAGCGCAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCCCGGCTCAACCGGGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAGAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGATATATGGAGGAACACCAGTGGCGAAGGCGGCTCTCTGGTCTGTAACTGACGCTGAGGCTCGA"
## [3] "ATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGA"

W4.4.3 cbindとwrite.table関数

　cbind関数を用いてファイル出力させたい情報を列方向で連結し、それをwrite.table関数でファイルに保存しています。図5は、得られたファイルをExcelで開いてスクリーンショットをとったものになります。念のため、ここでもタブ区切りテキストファイル版（output_JSLAB24.txt）とxlsx版（output_JSLAB24.xlsx）をダウンロードできるようにしておきます。

tmp <- cbind(ASV, t(seqtab.nochim), taxa2)# ASV、seqtab.nochim行列の転置、taxa2を列方向で連結してtmpに格納
dim(tmp)                                  # tmpの行数を列数を表示

## [1] 580  20

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という結果を得てからにしたり、あるいはもっと高度な方法で連結させます。

identical(colnames(seqtab.nochim), rownames(taxa2))# seqtab.nochimの列名とtaxa2の行名の2つのベクトルが同一かどうかを確認

## [1] TRUE

　いきなり上記のやつだと挙動がよくわからないかもしれませんので、念のため以下に簡単な例を示しておきます。hoge2とhoge3の比較時のみ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つのベクトルが同一かどうかを確認

## [1] FALSE

identical(hoge1, hoge3)                   # hoge1とhoge3の2つのベクトルが同一かどうかを確認

## [1] FALSE

identical(hoge1, hoge4)                   # hoge1とhoge4の2つのベクトルが同一かどうかを確認

## [1] FALSE

identical(hoge2, hoge3)                   # hoge2とhoge3の2つのベクトルが同一かどうかを確認

## [1] TRUE

identical(hoge2, hoge4)                   # hoge2とhoge4の2つのベクトルが同一かどうかを確認

## [1] FALSE

identical(hoge3, hoge4)                   # hoge3とhoge4の2つのベクトルが同一かどうかを確認

## [1] FALSE

おわりに

W5.1 参考にしたページ

• DADA2 Pipeline Tutorial (1.16)
• dada2のIntroduction to dada2ページ

W5.2 解析環境の確認

　packageVersion関数は、個別のインストール済みパッケージのバージョン情報を知りたい場合に利用します。今回動作確認したdada2のバージョンは1.30.0だということがわかります。

packageVersion("dada2")                   # dada2のバージョンを表示

## [1] '1.30.0'

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

sessionInfo()                             # 解析環境を表示

## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 22631)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Japanese_Japan.utf8  LC_CTYPE=Japanese_Japan.utf8   
## [3] LC_MONETARY=Japanese_Japan.utf8 LC_NUMERIC=C                   
## [5] LC_TIME=Japanese_Japan.utf8    
## 
## time zone: Asia/Tokyo
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] Biostrings_2.70.1   GenomeInfoDb_1.38.5 XVector_0.42.0     
## [4] IRanges_2.36.0      S4Vectors_0.40.2    BiocGenerics_0.48.1
## [7] dada2_1.30.0        Rcpp_1.0.12        
## 
## loaded via a namespace (and not attached):
##  [1] SummarizedExperiment_1.32.0 gtable_0.3.4               
##  [3] hwriter_1.3.2.1             xfun_0.41                  
##  [5] bslib_0.6.1                 ggplot2_3.4.4              
##  [7] latticeExtra_0.6-30         Biobase_2.62.0             
##  [9] lattice_0.21-9              vctrs_0.6.5                
## [11] tools_4.3.2                 bitops_1.0-7               
## [13] generics_0.1.3              parallel_4.3.2             
## [15] tibble_3.2.1                klippy_0.0.0.9500          
## [17] fansi_1.0.6                 highr_0.10                 
## [19] pkgconfig_2.0.3             Matrix_1.6-1.1             
## [21] RColorBrewer_1.1-3          RcppParallel_5.1.8         
## [23] assertthat_0.2.1            lifecycle_1.0.4            
## [25] GenomeInfoDbData_1.2.11     farver_2.1.1               
## [27] stringr_1.5.1               deldir_2.0-4               
## [29] compiler_4.3.2              Rsamtools_2.18.0           
## [31] munsell_0.5.0               codetools_0.2-19           
## [33] snow_0.4-4                  htmltools_0.5.7            
## [35] sass_0.4.8                  RCurl_1.98-1.14            
## [37] yaml_2.3.8                  pillar_1.9.0               
## [39] crayon_1.5.2                jquerylib_0.1.4            
## [41] BiocParallel_1.36.0         cachem_1.0.8               
## [43] DelayedArray_0.28.0         ShortRead_1.60.0           
## [45] abind_1.4-5                 tidyselect_1.2.0           
## [47] digest_0.6.34               stringi_1.8.3              
## [49] reshape2_1.4.4              dplyr_1.1.4                
## [51] labeling_0.4.3              fastmap_1.1.1              
## [53] grid_4.3.2                  colorspace_2.1-0           
## [55] cli_3.6.2                   SparseArray_1.2.3          
## [57] magrittr_2.0.3              S4Arrays_1.2.0             
## [59] utf8_1.2.4                  withr_3.0.0                
## [61] scales_1.3.0                rmarkdown_2.25             
## [63] jpeg_0.1-10                 matrixStats_1.2.0          
## [65] interp_1.1-6                png_0.1-8                  
## [67] evaluate_0.23               knitr_1.45                 
## [69] GenomicRanges_1.54.1        rlang_1.1.3                
## [71] glue_1.7.0                  jsonlite_1.8.8             
## [73] plyr_1.8.9                  R6_2.5.1                   
## [75] MatrixGenerics_1.14.0       GenomicAlignments_1.38.2   
## [77] zlibbioc_1.48.0