| サンプル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 |
ここでは、タブ区切りテキストファイルとxlsx形式ファイルの読込を行います。入力ファイルが同じでも、読み込み時のオプション次第で、読み込み後のオブジェクト（この場合はhoge）の行数や列数、行名や列名の有無なども変わりますので、ご注意ください。
第24回の図5で見えているファイルの読込を行うのが主な目的です。タブ区切りテキストファイル（output_JSLAB24.txt）の読込、xlsx形式ファイル（output_JSLAB24.xlsx）の読込の2パターンを示します。但し後者の場合は、xlsx形式ファイルを読み込むための関数を提供しているパッケージのロードからになります。ここではそのうちの1つであるopenxlsxパッケージのロードから行います。ヒトによってこのパッケージのインストールができているかどうかは異なります。それゆえ、若干挙動も異なりますが、W1.1を一通りやれば全員タブ区切りテキストファイルでもxlsx形式ファイルでも読み込めるようになると思います。
read.table関数 read.table関数は、様々なテキストファイルを読み込むための関数です。header=TRUEは、1行目を列名として取り扱うという宣言です。row.names=1は、1列目を行名として取り扱うという宣言です。sep="\t"は、区切り文字がタブ(\t)であるという宣言です。例えばもし半角スペースで中身を区切りたい場合は、sep=" "とします。また、コンマ区切りにして読み込みたい場合は、sep=","にします。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t")# in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
dim(hoge) # hogeの行数と列数## [1] 580 19読み込み後のオブジェクトの行数は580、列数は19であることがわかります。 一番右から2番目のGenus列が属名、一番右のSpecies列が種名ですので、この中にどれくらい計19属（22種）の既知微生物群集HM-280（Gevers et al., PLoS Comput Biol., 2012）が含まれているのかなどを評価していくことになります。
read.xlsx関数 read.xlsx関数は、openxlsxパッケージから提供する、xlsx形式ファイルを読み込むための関数です。colNames=Tオプションは、1行目を列名として取り扱うという宣言になります（1行目が列名情報ではない場合はcolNames=Fとします）。それゆえ、openxlsxパッケージがインストールされていないヒトは、このパッケージのロード部分以降でエラーになります。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.xlsx" # 読み込みたいファイル情報をin_fに格納
library(openxlsx) # openxlsxパッケージのロード
hoge <- read.xlsx(in_f, colNames=T, rowNames=T)# in_fを、1行目を列名として読み込んだ結果をhogeに格納
dim(hoge) # hogeの行数と列数## [1] 580 19連載第18回でも言及した、openxlsxパッケージのインストールを行います。このパッケージ自体はCRANというリポジトリから提供されていますので、本当はもう少しシンプルに書けます。しかし、dada2のようにBioconductorというリポジトリから提供されているパッケージに対しては、その書き方は通用しません。以下はどちらのリポジトリから提供されているパッケージに対しても適用可能なインストール手段になります。
if (!requireNamespace("BiocManager", quietly=T))# もしBiocManagerパッケージが入っていなければ...
install.packages("BiocManager") # BiocManagerパッケージをインストールせよ
BiocManager::install("openxlsx", update=F)# openxlsxパッケージをインストールせよ実行時に「警告: パッケージ ‘openxlsx’ はバージョン 4.x.x の R の下で造られました 」のような警告が出ることがありますが、これは気にする必要はありません。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.xlsx" # 読み込みたいファイル情報をin_fに格納
library(openxlsx) # openxlsxパッケージのロード
hoge <- read.xlsx(in_f, colNames=T, rowNames=T)# in_fを、1行目を列名として読み込んだ結果をhogeに格納
dim(hoge) # hogeの行数と列数## [1] 580 19なお、上記の2-3行目の部分については、以下のようにまとめて書くヒトもいます。見慣れないうちは相当困惑しますので、今のうちに慣れておきましょう。「パッケージ名::そのパッケージが提供する関数名」というフォーマットになります。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.xlsx" # 読み込みたいファイル情報をin_fに格納
hoge <- openxlsx::read.xlsx(in_f, colNames=T, rowNames=T)# in_fを、1行目を列名として読み込んだ結果をhogeに格納
dim(hoge) # hogeの行数と列数## [1] 580 19 連載第23回の表1のタブ区切りテキストファイル（JSLAB23_table1.txt）を読込みます。このファイルの中身は、計19属（22種）の既知微生物群集HM-280（Gevers et al., PLoS
Comput Biol.,
2012）の情報になります。header=TRUEは、1行目を列名として取り扱うという宣言です。row.names=1は、1列目を行名として取り扱うという宣言です。sep="\t"は、区切り文字がタブ(\t)であるという宣言です。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB23_table1.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t") #in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
dim(hoge) # hogeの行数と列数## [1] 22 3読み込み後の行数は22、列数は3であることがわかります。以下は、hogeの最初の4行分を表示して確認しているだけです。
head(hoge, n = 4) # hoge中の最初の4つを確認## Organism copy_num_16S.rRNA
## 1 Acinetobacter baumannii, strain 5377 5
## 2 Actinomyces odontolyticus, strain 1A.21 2
## 3 Bacillus cereus, strain NRS 248 12
## 4 Bacteroides vulgatus, strain NTC 11154 7
## RefSeq_ID
## 1 NC_009085
## 2 NZ_DS264585, NZ_DS264586
## 3 NC_003909
## 4 NC_009614ここでは、W1で読み込んだ2つのファイルの中から必要な部分だけを抽出する操作などを行います。W1.1で読み込んだdada2実行結果ファイル（output_JSLAB24.txt）とW1.2で読み込んだ計19属（22種）の既知微生物群集の情報を含んだファイル（JSLAB23_table1.txt）
W1.1で読み込んだdada2実行結果ファイル（output_JSLAB24.txt）を再度読み込むところからスタートし、計12サンプル分のASVごとの出現頻度情報に相当する列（“SRR…”という名前の12列）をdataというオブジェクト名で、Genus列をgenusというオブジェクト名で、そしてSpecies列をspeciesというオブジェクト名で得ます。
colnames関数 ここでは、W1.1.1を再度実行した後、colnames関数を用いて行列hogeオブジェクトの列名情報を表示させます。read.table関数実行時に、header=TRUEオプションをつけているので、行列hogeに列名が付与されています。もしheader=TRUEオプションをつけずに読み込んでいれば、（hogeに列名情報が付与されないので）colnames関数を実行しても何も表示されません。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t")# in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
dim(hoge) # hogeの行数と列数## [1] 580 19colnames(hoge) # hogeの列名情報## [1] "SRR6325813" "SRR6325815" "SRR6325816" "SRR6325817" "SRR6325828"
## [6] "SRR6325829" "SRR6325830" "SRR6325831" "SRR6325859" "SRR6325860"
## [11] "SRR6325861" "SRR6325862" "Kingdom" "Phylum" "Class"
## [16] "Order" "Family" "Genus" "Species"grep関数 次に、計12サンプル分のASVごとの出現頻度情報に相当する列（“SRR…”という名前の12列）をdataというオブジェクトに格納します。作業自体はdata <- hoge[, 1:12]でもよいのですが、ここではせっかくですのでgrep関数を利用します。grep関数は、入力が文字列ベクトルの場合、検索したい文字列を含む要素番号を返します。まずやりたいことは、colnames(hoge)という文字列ベクトルの中から、"SRR"という文字列を含む要素番号情報を得ることです。pattern = "SRR"が検索文字列情報、x = colnames(hoge)が被検索対象になります。
colnames(hoge) # hogeの列名情報## [1] "SRR6325813" "SRR6325815" "SRR6325816" "SRR6325817" "SRR6325828"
## [6] "SRR6325829" "SRR6325830" "SRR6325831" "SRR6325859" "SRR6325860"
## [11] "SRR6325861" "SRR6325862" "Kingdom" "Phylum" "Class"
## [16] "Order" "Family" "Genus" "Species"obj <- grep(pattern = "SRR", x = colnames(hoge)) # colnames(hoge)の中から"SRR"という文字列を含む要素番号情報をobjに格納
obj # objの中身を確認## [1] 1 2 3 4 5 6 7 8 9 10 11 12 objの中身は1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
12ですので、確かに検索文字列である"SRR"を含む列番号情報を抽出できていることがわかります。
data） 列数が19のhoge行列の中から、objで指定した列番号のみを抽出した結果をdataオブジェクトに格納します。dim(data)やcolnames(data)の結果からも、確かに意図通りの列を抽出できていることがわかります。
dim(hoge) # hogeの行数と列数## [1] 580 19data <- hoge[, obj] # hoge行列中のobjで指定した列のみ抽出した結果をdataに格納
dim(data) # dataの行数と列数## [1] 580 12colnames(data) # dataの列名情報## [1] "SRR6325813" "SRR6325815" "SRR6325816" "SRR6325817" "SRR6325828"
## [6] "SRR6325829" "SRR6325830" "SRR6325831" "SRR6325859" "SRR6325860"
## [11] "SRR6325861" "SRR6325862"genus） 次は、580行 \(\times\)
19列のhogeオブジェクトの中から、Genus列をgenusというオブジェクト名で得ます。hogeオブジェクトがデータフレーム形式であれば（class(hoge)の実行結果が”data.frame”になっていれば）、hoge$Genusのような感じで$を間に挟んで行列hogeの列名情報を指定すれば、任意の列情報を抽出することができます。これは、列名情報がわかっていて列数が多い行列データの場合に、列番号情報を覚えておく必要がない便利な書き方です。
class(hoge) # hogeのクラス属性## [1] "data.frame"genus <- hoge$Genus # hoge中のGenus列を抽出した結果をgenusに格納
head(genus) # genus中の最初の6つを表示## [1] "Neisseria" "Enterococcus" "Escherichia-Shigella"
## [4] "Pseudomonas" NA "Acinetobacter"species） 次は、580行 \(\times\)
19列のhogeオブジェクトの中から、Species列をspeciesというオブジェクト名で得ます。hogeオブジェクトがデータフレーム形式であれば（class(hoge)の実行結果が”data.frame”になっていれば）、hoge$Speciesのような感じで$を間に挟んで行列hogeの列名情報を指定すれば、任意の列情報を抽出することができます。
class(hoge) # hogeのクラス属性## [1] "data.frame"species <- hoge$Species # hoge中のSpecies列を抽出した結果をspeciesに格納
head(species) # species中の最初の6つを表示## [1] "meningitidis" NA NA NA NA
## [6] NAW2.1の必要最小限のスクリプトを以下に示します。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t")# in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
obj <- grep(pattern = "SRR", x = colnames(hoge)) # colnames(hoge)の中から"SRR"という文字を含む要素番号情報をobjに格納
data <- hoge[, obj] # hoge行列中のobjで指定した列のみ抽出した結果をdataに格納
genus <- hoge$Genus # hoge中のGenus列を抽出した結果をgenusに格納
species <- hoge$Species # hoge中のSpecies列を抽出した結果をspeciesに格納W2.1で得た結果をまとめると以下のようになります。
data：数値行列
nrow(data)：行数（ASVの総数）は580ncol(data)：列数（サンプル数）は12genus：文字列ベクトル
length(genus)：要素数（ASVの総数）は580species：文字列ベクトル
length(species)：要素数（ASVの総数）は580 W1.2で読み込んだ菌株情報ファイル（JSLAB23_table1.txt）を再度読み込むところからスタートします。この場合は1列目のOrganism列がとりあえず欲しい情報になりますので、抽出結果をorganismオブジェクトとして得ます。次に、organismオブジェクトの中から、コンマ(,)の右側の菌株名情報を除外し、属名（genus_truth）と種名（species_truth）の情報のみをうまく抽出するのがここでの目的になります。属の場合はgenus_truth（正解情報）とW2.1.6で得たgenusを、そして種の場合はspecies_truth
（正解情報）とW2.1.6で得たspeciesを比較して、dada2実行結果の中にどれくらい計19属（22種）の既知微生物群集HM-280（Gevers et al., PLoS
Comput Biol., 2012）が含まれているのかなどを評価します。
read.table関数 ここでは、W1.2を再度実行した後、colnames関数を用いて行列hogeオブジェクトの列名情報を表示させるところまで行います。read.table関数実行時にheader=TRUEオプションをつけているので、行列hogeに列名が付与されています。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB23_table1.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t") #in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
dim(hoge) # hogeの行数と列数## [1] 22 3colnames(hoge) # hogeの列名情報## [1] "Organism" "copy_num_16S.rRNA" "RefSeq_ID" read.table関数実行時にrow.names=1オプションをつけているので、行列hogeに行名が付与されています。そのため、元ファイルは4列から構成されますが、dim(hoge)で得られる列数は3になります。
organism） 22行 \(\times\)
3列のhogeオブジェクトの中から、Organism列をorganismというオブジェクト名で得ます。hogeオブジェクトがデータフレーム形式の場合（class(hoge)の実行結果が”data.frame”になっている場合）、hoge$Organismのような感じで$を間に挟んで行列hogeの列名情報を指定すれば、任意の列情報を抽出することができます。
class(hoge) # hogeのクラス属性## [1] "data.frame"organism <- hoge$Organism # hoge中のOrganism列を抽出した結果をorganismに格納
head(organism) # organism中の最初の6つを表示## [1] "Acinetobacter baumannii, strain 5377 "
## [2] "Actinomyces odontolyticus, strain 1A.21 "
## [3] "Bacillus cereus, strain NRS 248 "
## [4] "Bacteroides vulgatus, strain NTC 11154 "
## [5] "Bifidobacterium adolescentis, strain E194a "
## [6] "Clostridium beijerinckii, strain NCIMB 8052 " organismオブジェクトを眺めると、コンマ(,)の右側が菌株名情報、左側が属名と種名になっていることがわかります。例えば1番目の生物種（Acinetobacter
baumannii, strain 5377
）の場合、属名はAcinetobacter、種名はbaumannii、菌株名はstrain
5377だと読み解きます。従って、organismオブジェクトから、属名と種名の情報のみをうまく抽出すれば、W2.1.6で得たgenusやspeciesオブジェクトと文字列の完全一致での比較ができます。
strsplit関数 第23回のW5.3でも軽く触れていますが、strsplit関数は、文字列ベクトルを入力として、任意の区切り文字で文字列を分割した結果をリスト形式で返します。ここでは、organismオブジェクトを入力として、区切り文字をコンマ（,）にして文字列を分割します。
hoge <- strsplit(x = organism, split = ",")# organism中の文字列ベクトルをコンマで分割した結果をhogeに格納
head(hoge) # hoge中の最初の6つを表示## [[1]]
## [1] "Acinetobacter baumannii" " strain 5377 "
##
## [[2]]
## [1] "Actinomyces odontolyticus" " strain 1A.21 "
##
## [[3]]
## [1] "Bacillus cereus" " strain NRS 248 "
##
## [[4]]
## [1] "Bacteroides vulgatus" " strain NTC 11154 "
##
## [[5]]
## [1] "Bifidobacterium adolescentis" " strain E194a "
##
## [[6]]
## [1] "Clostridium beijerinckii" " strain NCIMB 8052 "分割後は、「属名 種名」と「菌株名」にうまく分かれています。
sapply関数 第23回のW5.4と全く同じ使い方ですが、sapply関数は、リスト形式のオブジェクトに対して逐次的な処理を行いたい場合に用いられます。ベクトル中の指定した位置の情報を抽出する[という関数を実行し、（1番目の要素を抽出したいので）1という情報をオプションとして与えています。
hoge2 <- sapply(X = hoge, FUN = `[`, 1) # hoge中のベクトル形式の各要素から、1番目の要素の抽出結果をhoge2に格納
head(hoge2) # hoge2中の最初の6つを表示## [1] "Acinetobacter baumannii" "Actinomyces odontolyticus"
## [3] "Bacillus cereus" "Bacteroides vulgatus"
## [5] "Bifidobacterium adolescentis" "Clostridium beijerinckii"strsplit関数 さきほど得たhoge2オブジェクトは「属名
種名」となってます。さらにこれを分割すべく、ここでは、hoge2オブジェクトを入力として、区切り文字を半角スペース（
）にして文字列を分割した結果をhoge3として得ます。
hoge3 <- strsplit(x = hoge2, split = " ")# hoge2中の文字列ベクトルを半角スペースで分割した結果をhoge3に格納
head(hoge3) # hoge3中の最初の6つを表示## [[1]]
## [1] "Acinetobacter" "baumannii"
##
## [[2]]
## [1] "Actinomyces" "odontolyticus"
##
## [[3]]
## [1] "Bacillus" "cereus"
##
## [[4]]
## [1] "Bacteroides" "vulgatus"
##
## [[5]]
## [1] "Bifidobacterium" "adolescentis"
##
## [[6]]
## [1] "Clostridium" "beijerinckii"sapply関数 W2.2.4と同じノリですが、hoge3を入力として、1番目の要素を属名、2番目の要素を種名として得ます。今取り扱っているのは、計19属（22種）の既知微生物群集HM-280（Gevers et al., PLoS
Comput Biol.,
2012）の情報です。つまり、ここで得る属名と種名が、dada2実行結果にどれだけ含まれているかを調べるための答え合わせ用の正解情報になります。従って、ここではそれぞれのオブジェクト名をgenus_truthとspecies_truthとします。こうすることで、dada2実行結果として得た580個の要素からなる属名情報のgenusおよび種名情報のspeciesと区別することができます。
genus_truth <- sapply(X = hoge3, FUN = `[`, 1) # hoge3中のベクトル形式の各要素から、1番目の要素の抽出結果をgenus_truthに格納
species_truth <- sapply(X = hoge3, FUN = `[`, 2)# hoge3中のベクトル形式の各要素から、2番目の要素の抽出結果をspecies_truthに格納
head(genus_truth) # genus_truth中の最初の6つを表示## [1] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides"
## [5] "Bifidobacterium" "Clostridium"head(species_truth) # species_truth中の最初の6つを表示## [1] "baumannii" "odontolyticus" "cereus" "vulgatus"
## [5] "adolescentis" "beijerinckii"unique関数 連載第23回の表1、および計19属（22種）の既知微生物群集HM-280（Gevers et al., PLoS
Comput Biol.,
2012）という言葉からも想像がつきますが、計22個の要素からなるspecies_truth中のユニークな要素数は22個です。その一方で、計22個の要素からなるgenus_truth中のユニークな要素数は19個です。この理由は、18～19番目が同じStaphylococcus属、20～22番目が同じStreptococcus属だからです。第24回のW3.3.5でも利用していますが、unique関数は、入力ベクトル中のユニークな要素を返します。ここではそれを利用して、unique関数をもう一度同じオブジェクト名で上書き保存するテクニックを使って（その分不可逆にはなりますが…）、任意のベクトルをユニークな要素のみにします。
length(genus_truth) # genus_truth中の要素数## [1] 22length(species_truth) # species_truth中の要素数## [1] 22genus_truth # genus_truthを確認## [1] "Acinetobacter" "Actinomyces" "Bacillus"
## [4] "Bacteroides" "Bifidobacterium" "Clostridium"
## [7] "Deinococcus" "Enterococcus" "Escherichia"
## [10] "Helicobacter" "Lactobacillus" "Listeria"
## [13] "Neisseria" "Porphyromonas" "Propionibacterium"
## [16] "Pseudomonas" "Rhodobacter" "Staphylococcus"
## [19] "Staphylococcus" "Streptococcus" "Streptococcus"
## [22] "Streptococcus"species_truth # species_truthを確認## [1] "baumannii" "odontolyticus" "cereus" "vulgatus"
## [5] "adolescentis" "beijerinckii" "radiodurans" "faecalis"
## [9] "coli" "pylori" "gasseri" "monocytogenes"
## [13] "meningitides" "gingivalis" "acnes" "aeruginosa"
## [17] "sphaeroides" "aureus" "epidermidis" "agalactiae"
## [21] "mutans" "pneumoniae"unique(genus_truth) # genus_truth中のユニークな要素## [1] "Acinetobacter" "Actinomyces" "Bacillus"
## [4] "Bacteroides" "Bifidobacterium" "Clostridium"
## [7] "Deinococcus" "Enterococcus" "Escherichia"
## [10] "Helicobacter" "Lactobacillus" "Listeria"
## [13] "Neisseria" "Porphyromonas" "Propionibacterium"
## [16] "Pseudomonas" "Rhodobacter" "Staphylococcus"
## [19] "Streptococcus"unique(species_truth) # species_truth中のユニークな要素## [1] "baumannii" "odontolyticus" "cereus" "vulgatus"
## [5] "adolescentis" "beijerinckii" "radiodurans" "faecalis"
## [9] "coli" "pylori" "gasseri" "monocytogenes"
## [13] "meningitides" "gingivalis" "acnes" "aeruginosa"
## [17] "sphaeroides" "aureus" "epidermidis" "agalactiae"
## [21] "mutans" "pneumoniae"genus_truth <- unique(genus_truth) # genus_truth中のユニークな要素をgenus_truthに格納
species_truth <- unique(species_truth) # species_truth中のユニークな要素をspecies_truthに格納
length(genus_truth) # genus_truth中の要素数## [1] 19length(species_truth) # species_truth中の要素数## [1] 22 上記では、unique関数実行結果が不変のspecies_truthも意図的に含めています。これは、全体を統一的な書き方にすることを意図しています。unique関数実行結果が不変な場合は、含めても特に悪さをしないので書いてもいいのです。このように、念のための確認作業はわりとよくなされます。
W2.2の必要最小限のスクリプトを以下に示します。2行目と4行目で同じオブジェクト名（hoge）を重複して利用していますが、一時利用目的であれば、このように上書きしても特に問題ありません。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB23_table1.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t") #in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
organism <- hoge$Organism # hoge中のOrganism列を抽出した結果をorganismに格納
hoge <- strsplit(x = organism, split = ",")# organism中の文字列ベクトルをコンマで分割した結果をhogeに格納
hoge2 <- sapply(X = hoge, FUN = `[`, 1) # hoge中のベクトル形式の各要素から、1番目の要素の抽出結果をhoge2に格納
hoge3 <- strsplit(x = hoge2, split = " ") # hoge2中の文字列ベクトルを半角スペースで分割した結果をhoge3に格納
genus_truth <- sapply(X = hoge3, FUN = `[`, 1) # hoge3中のベクトル形式の各要素から、1番目の要素の抽出結果をgenus_truthに格納
species_truth <- sapply(X = hoge3, FUN = `[`, 2)# hoge3中のベクトル形式の各要素から、2番目の要素の抽出結果をspecies_truthに格納
genus_truth <- unique(genus_truth) # genus_truth中のユニークな要素をgenus_truthに格納
species_truth <- unique(species_truth) # species_truth中のユニークな要素をspecies_truthに格納W2.2で得た結果をまとめると以下のようになります。
genus_truth：文字列ベクトル
length(genus_truth)：要素数（ユニークな属の数）は19species_truth：文字列ベクトル
length(species_truth)：要素数（ユニークな種の数）は22W2の必要最小限のスクリプトを以下に示します。これは、W2.1.6とW2.2.8をマージしただけです。
# W2.1 dada2実行結果ファイルのほう
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t")# in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
obj <- grep(pattern = "SRR", x = colnames(hoge)) # colnames(hoge)の中から"SRR"という文字を含む要素番号を抽出した結果をobjに格納
data <- hoge[, obj] # hoge行列中のobjで指定した列のみ抽出した結果をdataに格納
genus <- hoge$Genus # hoge中のGenus列を抽出した結果をgenusに格納
species <- hoge$Species # hoge中のSpecies列を抽出した結果をspeciesに格納
# W2.2 菌株情報ファイルのほう
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB23_table1.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t") #in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
organism <- hoge$Organism # hoge中のOrganism列を抽出した結果をorganismに格納
hoge <- strsplit(x = organism, split = ",")# organism中の文字列ベクトルをコンマで分割した結果をhogeに格納
hoge2 <- sapply(X = hoge, FUN = `[`, 1) # hoge中のベクトル形式の各要素から、1番目の要素の抽出結果をhoge2に格納
hoge3 <- strsplit(x = hoge2, split = " ") # hoge2中の文字列ベクトルを半角スペースで分割した結果をhoge3に格納
genus_truth <- sapply(X = hoge3, FUN = `[`, 1) # hoge3中のベクトル形式の各要素から、1番目の要素の抽出結果をgenus_truthに格納
species_truth <- sapply(X = hoge3, FUN = `[`, 2)# hoge3中のベクトル形式の各要素から、2番目の要素の抽出結果をspecies_truthに格納
genus_truth <- unique(genus_truth) # genus_truth中のユニークな要素をgenus_truthに格納
species_truth <- unique(species_truth) # species_truth中のユニークな要素をspecies_truthに格納 ここではまず、連載第24回の図5で見えているdada2実行結果ファイル（output_JSLAB24.txt）のGenus列およびSpecies列から得た「580個の要素からなるgenusおよびspeciesオブジェクト」、そして計12サンプル分のASVごとの出現頻度情報に相当する列（“SRR…”という名前の12列）を格納した「580行
\(\times\)
12列のdataオブジェクト」を主な入力として、dada2実行結果の要約情報を見ていきます（W3.1）。
連載第23回の表1に相当する菌株情報ファイル（JSLAB23_table1.txt）には、計19属（22種）の既知微生物群集HM-280（Gevers et al., PLoS
Comput Biol.,
2012）に関する情報が含まれています。この情報と突き合わせて、どれだけ本来検出されるべきものが正しく検出できているかを評価していきます。実際に使う情報は、19個の要素からなる文字列ベクトルのgenus_truthオブジェクトと22個の要素からなる文字列ベクトルのspecies_truth
オブジェクトになります。
まず、dada2実行結果ファイル（output_JSLAB24.txt）の要約情報を見ていきます。重要なのは、予めExcelなどで眺めて全体像を把握しておくことです。これは、どのような結果が得られれば正しく処理されていると判断してよいかを知っておくことと同義です。具体的には、例えば「Genus列よりもSpecies列のほうが圧倒的に”NA”の要素が多い」とか「“NA”というのは第24回のW4でリファレンス配列中にヒットしなかったものに相当する」といったことを理解しておくことです。
na.omit関数 na.omit関数は、ベクトルが入力の場合は”NA”ではない要素番号情報を返し、行列が入力の場合は”NA”を含まない行を返します。ここではまず、“NA”の要素数が非常に多いspeciesベクトルを入力として、na.omit関数を実行します。
na.omit(species) # speciesから"NA"の要素を除去## [1] "meningitidis" "gingivalis" "vulgatus"
## [4] "odontolyticus" "radiodurans" "pylori"
## [7] "vulgatus" "vulgatus" "mutans"
## [10] "aerosaccus" "clevelandensis" "beijerinckii"
## [13] "aureus" "durum" "mucilaginosa"
## [16] "maltophilia" "lwoffii" "hydrogenalis"
## [19] "dentocariosa" "aquatilis" "vaginae"
## [22] "melaninogenica" "mucilaginosa" "nagyae"
## [25] "elegans" "uniformis" "octavius"
## [28] "magna" "oligocarboniphilum" "timonensis"
## [31] "carboxydivorans" "bromii" "carnosum"
## [34] "psychrotolerans" "hominis" "massiliensis"
## [37] "salivarius" "nigrescens" "geothermalis"
## [40] "kroppenstedtii" "stomatis" "mucilaginosa"
## [43] "formosensis" "dankookensis" "saccharivorans"
## [46] "neonati" "profundi" "perfoetens"
## [49] "edouardi" "wisconsinensis" "caenicola"
## [52] "massiliensis" "splanchnicus" "acetatigenes"
## [55] "pigrum" "laidlawii" "shahii"
## [58] "gandavensis" "hartonius"
## attr(,"na.action")
## [1] 2 3 4 5 6 8 9 12 13 14 17 20 21 23 24 25 26 27
## [19] 28 29 30 31 32 33 35 36 37 38 39 40 41 42 43 44 47 48
## [37] 49 50 51 52 53 56 57 58 59 60 61 62 64 65 67 68 69 71
## [55] 72 73 74 75 77 78 79 80 81 82 84 85 86 87 88 89 90 91
## [73] 93 94 95 97 98 99 100 101 102 103 104 105 107 108 109 110 111 112
## [91] 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
## [109] 133 134 135 137 138 139 140 141 142 143 144 146 147 148 149 151 152 153
## [127] 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171
## [145] 172 173 174 175 176 177 178 180 181 183 184 185 186 187 188 189 190 191
## [163] 192 193 194 195 196 197 200 201 202 203 204 205 206 207 208 209 210 211
## [181] 212 213 214 215 216 217 219 220 221 222 223 224 225 226 227 228 230 231
## [199] 233 235 236 237 238 239 240 241 242 243 246 247 248 249 250 251 252 253
## [217] 254 255 256 257 258 259 260 261 262 263 264 265 266 267 269 270 271 272
## [235] 273 274 275 276 277 279 280 281 283 284 285 288 289 290 291 292 293 294
## [253] 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312
## [271] 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330
## [289] 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 347 348 349
## [307] 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367
## [325] 368 369 370 371 372 373 374 375 377 379 380 381 382 383 384 385 386 387
## [343] 388 389 390 391 392 393 394 395 396 397 399 400 401 402 403 405 406 407
## [361] 408 409 411 412 413 414 415 416 417 419 420 422 423 424 425 427 428 429
## [379] 430 431 432 433 434 436 438 439 440 441 442 443 444 445 446 447 448 449
## [397] 450 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468
## [415] 469 470 471 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487
## [433] 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505
## [451] 506 507 508 510 511 512 513 514 515 516 517 518 519 520 521 522 524 525
## [469] 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 543 544
## [487] 545 546 547 548 549 550 551 552 554 555 556 557 558 559 560 561 562 563
## [505] 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580
## attr(,"class")
## [1] "omit" 表示結果とdada2実行結果ファイル（output_JSLAB24.txt）を見比べることで、“NA”でない1番最初の要素がmeningitidis、2番目がgingivalisなのだと解釈します。そして、“NA”でない要素数は59なのだと読み解きます。ちなみに、この”NA”でない要素数は、na.omit(species)実行結果に対して、さらにベクトルの要素数を返すlength関数を実行すれば得られます。この59という数値が、種レベルで系統割り当てできたASV数に相当します。つまり、種レベルだと、580個のASVのうち、59個の系統割り当てができているということです。
length(na.omit(species)) # species中の"NA"以外の要素数## [1] 59 次に、属レベルでどれだけ系統割り当てができたかを調べます。dada2実行結果ファイル（output_JSLAB24.txt）中のGenus列情報に相当するgenusオブジェクトを入力として、同様の解析を行います。
length(na.omit(genus)) # genus中の"NA"以外の要素数## [1] 343属レベルだと、580個のASVのうち、343個の系統割り当てができていることがわかります。
getNonnaNum さきほどまでは1つのベクトルを入力として、na.omitからlengthという流れで情報を得ましたが、dada2実行結果ファイル（output_JSLAB24.txt）を直接読み込んだ結果を入力として、列ごとに同様の解析を行うこともできます。そうすることで、Genusよりも高い階級、つまり「界（Kingdom）・門（Phylum）・綱（Class）・目（Order）・科（Family）」）についても系統割り当て結果の全体像が見られます。ここでは、「NAではない（Non-NA）数の情報を得るための関数」という意味で、以下のgetNonnaNum関数を定義します。
getNonnaNum <- function(x) length(na.omit(x)) # "NA"以外の要素数を得るための関数を定義自作関数については、第24回ののW1.6.3、W2.5、W2.6などでも解説しています。
apply関数 apply関数は、行列を入力として、行ごとや列ごとに任意の処理を行うための関数です。ここでは、apply関数の中で、先ほど自作したgetNonnaNum関数を実行する例を示します。dada2実行結果ファイル（output_JSLAB24.txt）を読み込んだ直後のhogeオブジェクトを入力として（X = hoge）、列ごとに（MARGIN = 2）、getNonnaNum関数を適用する（FUN = getNonnaNum）という書き方にしています。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t")# in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
result <- apply(X = hoge, MARGIN = 2, FUN = getNonnaNum)# hogeの各列に対して、getNonnaNum関数を実行した結果をresultに格納
result # resultの中身## SRR6325813 SRR6325815 SRR6325816 SRR6325817 SRR6325828 SRR6325829 SRR6325830
## 580 580 580 580 580 580 580
## SRR6325831 SRR6325859 SRR6325860 SRR6325861 SRR6325862 Kingdom Phylum
## 580 580 580 580 580 580 568
## Class Order Family Genus Species
## 558 525 454 343 59 実行結果であるresultオブジェクトには本来不要な結果も含まれていますが、Genus列やSpecies列の結果も含まれますので、答え合わせになるメリットもあります。実際の作業としては、まずはこれをやって手順や結果に間違いがないことを確認します。例えば、“SRR”から始まる計12サンプルのASVの出現頻度情報部分は、全て580という数値になっていて妥当だと判断します。理由は、ここの元情報は出現頻度なので、0以上の整数なはずであり、原理的にNAの数が1つも含まれないはずだからです。また、界（Kingdom）が580になるのも妥当です。理由は、第24回のW4で行った系統割り当て時に利用したリファレンス配列は、全てバクテリア由来だからです。dada2実行結果ファイル（output_JSLAB24.txt）のKingdom列を実際に眺めてみると、全要素がBacteriaになっていることがわかります。こんな感じで確認しながら進めていく姿勢が大事です。
grep関数 先ほどの実行結果には、“SRR”から始まる計12サンプルのASVの出現頻度情報部分の結果も含まれていました。ここでは、“SRR”という文字を含まない列を除外したものに対して、先ほどと同様の作業を実行します。grep関数は、W2.1.2でも利用しましたが、「検索したい文字列を含む、文字列ベクトル中の要素番号」を返すのが基本形です。dada2実行結果ファイル（output_JSLAB24.txt）を読み込んだ直後のhogeオブジェクトの列名情報に相当するのがcolnames(hoge)になりますので、W2.1.2でも利用したgrep関数を利用して、“SRR”という文字を（pattern = "SRR"）を含まない（invert = TRUE）要素番号を取得するというテクニックを使っています。得られたobjが「今回抽出したい列番号情報」に相当しますので、apply関数の入力をX = hoge[, obj]とすれば目的を達成できます。
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t")# in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
obj <- grep(pattern = "SRR", x = colnames(hoge), invert = TRUE)# colnames(hoge)の中から、"SRR"という文字を含まない要素番号をobjに格納
obj # objの中身## [1] 13 14 15 16 17 18 19result <- apply(X = hoge[, obj], MARGIN = 2, FUN = getNonnaNum)# hoge[, obj]の各列に対して、getNonnaNum関数を実行した結果をresultに格納
result # resultの中身## Kingdom Phylum Class Order Family Genus Species
## 580 568 558 525 454 343 59ここまでの結果が得られるスクリプトが、一般に論文のSupplementなどで提供されるものになります。このスクリプトだけを最初に見せられると難解なイメージを持たれると思いますが、得たい出力結果と実際の入力情報との兼ね合いで、いろいろと小細工をしないといけないものです。
table関数 第24回のW1.6.1でも利用しましたが、table関数は、ベクトルを入力として与えると、要素の種類ごとの出現頻度を返します。それゆえ、table関数実行結果のベクトルの要素数は、入力ベクトル中のユニークな要素数（但しNAは除く）と等しくなります。ここではまず、580個の要素からなるspeciesオブジェクトを入力として実行します。
table(species) # species中の要素の種類ごとの出現頻度## species
## acetatigenes aerosaccus aquatilis aureus
## 1 1 1 1
## beijerinckii bromii caenicola carboxydivorans
## 1 1 1 1
## carnosum clevelandensis dankookensis dentocariosa
## 1 1 1 1
## durum edouardi elegans formosensis
## 1 1 1 1
## gandavensis geothermalis gingivalis hartonius
## 1 1 1 1
## hominis hydrogenalis kroppenstedtii laidlawii
## 1 1 1 1
## lwoffii magna maltophilia massiliensis
## 1 1 1 2
## melaninogenica meningitidis mucilaginosa mutans
## 1 1 3 1
## nagyae neonati nigrescens octavius
## 1 1 1 1
## odontolyticus oligocarboniphilum perfoetens pigrum
## 1 1 1 1
## profundi psychrotolerans pylori radiodurans
## 1 1 1 1
## saccharivorans salivarius shahii splanchnicus
## 1 1 1 1
## stomatis timonensis uniformis vaginae
## 1 1 1 1
## vulgatus wisconsinensis
## 3 1length(table(species)) # species中のユニークな要素数## [1] 54連載第23回の表1でも示されているとおり、今回のサンプルは計19属（22種）の既知微生物群集（HM-280）から構成されています。従って、この54個という結果を見ただけで、最低でも(54 - 22)個の本来は存在しないはずの生物種が検出されたと解釈します。「最低でも」とわざわざ書いたのは、現時点ではこの22種全てを正しく検出できているかどうかが不明だからです。
なお、この54という数値は、以下でも得られます。na.omit関数を挟んでいるのは、まずは”NA”の結果を除外しておきたいからです。もしlength(unique(species))と書くと、“NA”の出現回数の結果も含まれてしまうため、55という数値になって、上記と合わなくなるからです。このように、関数によってデフォルトの”NA”への対処方針は異なりますのでご注意ください。
length(unique(na.omit(species))) # speciesベクトル中のユニークな要素数## [1] 54 次に、580個の要素からなるgenusオブジェクトを入力として実行します。この222個という結果を見ただけで、最低でも(222
- 19)個の本来は存在しないはずの属が検出されたと解釈します。
length(table(genus)) # genusベクトル中のユニークな要素数## [1] 222W3.1で得た結果をまとめると以下のようになります。「“NA”を除く要素数」よりも「“NA”を除くユニークな要素数」が少なくなるのは、ASV配列が微妙に異なっていても、同じ属または種に割り当てられることもあるからだと解釈すればよいです。
genus：文字列ベクトル
length(genus)：要素数（ASVの総数）は580length(na.omit(genus))：“NA”を除く要素数（系統割り当てできたASV数）は343length(table(genus))：“NA”を除くユニークな要素数（系統割り当てできたユニークな属の数）は222species：文字列ベクトル
length(species)：要素数（ASVの総数）は580length(na.omit(species))：“NA”を除く要素数（系統割り当てできたASV数）は59length(table(species))：“NA”を除くユニークな要素数（系統割り当てできたユニークな種の数）は54W3.1では、dada2実行結果として得られた580個のASV配列のうち、属レベルでは343個（ユニークな属の数は222個）、種レベルでは59個（ユニークな種の数は54個）について系統割り当てが行われたことを主に確認しました。ここでは、連載第23回の表1の、計19属（22種）の既知微生物群集（HM-280）をどれだけ検出できているかを調べます。例えば、属レベルの系統割り当て結果である222個のユニークな属数の中に、HM-280に属する計19属が全て含まれていれば、全て正しく検出できた（感度100%）と評価します。また、dada2結果の中にHM-280に属する計19属以外の余計なものが含まれていないかどうかも重要な指標といえます。もしdada2実行結果のASV配列全てがこの計19属のみから構成されていれば、適合率100%ということになります。ここでは、具体的な計算方法を含めて解説していきます。なお、用いる既知微生物群集情報は、W2.2で取得した以下の2つの文字列ベクトルです。
genus_truth：文字列ベクトル
length(genus_truth)：要素数（ユニークな属の数）は19species_truth：文字列ベクトル
length(species_truth)：要素数（ユニークな種の数）は22上記2つの正解情報が格納されたベクトルに対して、W2.1で得た以下の2つの予測結果（dada2実行結果）が格納されたベクトルと一致するものがあるかを調べるのがここでの作業になります。
genus：文字列ベクトル
length(genus)：要素数（ASVの総数）は580length(na.omit(genus))：“NA”を除く要素数（系統割り当てできたASV数）は343length(table(genus))：“NA”を除くユニークな要素数（系統割り当てできたユニークな属の数）は222species：文字列ベクトル
length(species)：要素数（ASVの総数）は580length(na.omit(species))：“NA”を除く要素数（系統割り当てできたASV数）は59length(table(species))：“NA”を除くユニークな要素数（系統割り当てできたユニークな種の数）は54is.element関数 is.element関数は、2つのベクトル(例：xとy)を入力として与えると、x中の要素ごとに、y中の要素と一致するものがあるかどうかを論理値（TRUE
or
FALSE）で返します。ここではまず、以下に示す2つの文字列ベクトルで考えます。
x <- c("A", "C", "D", "F", "J", "J") # 文字列ベクトルを作成してxに格納
y <- c("J", "C", "C", "G", "B") # 文字列ベクトルを作成してyに格納 is.element関数は、与えるベクトルの順番を気にする必要があります。例えば以下のようにis.element(x, y)として実行すると、実行結果の要素数は6になります。つまり、xの要素数と同じということです。なお、is.element(x, y)とx %in% yは同義です。出力結果の1番目の要素がFALSEになっているのは、x中の1番目の要素である”A”と同じものがyには存在しないことを意味します。同様に、出力結果の2番目の要素がTRUEになっているのは、x中の2番目の要素である”C”と同じものがyにも存在することを意味します。
is.element(x, y) # x中の各要素がyに含まれるか判定## [1] FALSE TRUE FALSE FALSE TRUE TRUEx %in% y # x中の各要素がyに含まれるか判定## [1] FALSE TRUE FALSE FALSE TRUE TRUE xとyを入れ替えてis.element(y, x)として実行すると、実行結果の要素数は5になります。なお、is.element(y, x)とy %in% xは同義です。出力結果の1番目の要素がTRUEになっているのは、y中の1番目の要素である”J”と同じものがxにも存在することを意味します。同様に、出力結果の5番目の要素がFALSEになっているのは、y中の5番目の要素である”B”と同じものがxには存在しないことを意味します。
is.element(y, x) # y中の各要素がxに含まれるか判定## [1] TRUE TRUE TRUE FALSE FALSEy %in% x # y中の各要素がxに含まれるか判定## [1] TRUE TRUE TRUE FALSE FALSE is.element関数を以下のように利用すれば、計19属のうち、どれを正しく検出できて、どれを検出できなかったかを調べることができます。genus_truthは、属の正解情報に相当する文字列ベクトルであり、要素数は19です。genusは、属の予測結果に相当する文字列ベクトルであり、要素数は580です。
obj <- is.element(genus_truth, genus) # genus_truth中の各要素がgenusに含まれるか判定した結果をobjに格納
obj # objの中身## [1] TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE
## [13] TRUE TRUE FALSE TRUE TRUE TRUE TRUE 出力結果はTRUE or
FALSEから構成される論理値ベクトルであり、要素数はlength(genus_truth)と同じ19です。このベクトル中のTRUEの要素数が正しく検出できた数に相当します。この程度であれば目でもカウントできますが、以下のようにsum関数で15だと分かります。
sum(obj) # obj中のTRUEの要素数## [1] 15 この時点で、計19属のうち15属を正しく検出できたことが分かったので、属レベルの検出感度は15
/ 19 = 0.7894737になります。感度（sensitivityまたはrecall）は、一般に「陽性と判定されるべきものを正しく陽性と判定する確率」と定義されます。この場合の「陽性（Positive）と判定されるべきもの」がgenus_truth中の19属になります。ここで陽性はdada2で検出された（系統割り当てできた）ユニークな属のことを指し、陰性（Negative）は検出されなかった（NAとなった）属のことを指します。また、「陽性と判定されるべきもののうち、正しく陽性と判定できたもの」のことをTrue
Positive (TP)といい
ます。この場合、TPの数は15個ということになります。一般に、「陽性と判定されるべきもののうち、正しく陽性と判定できなかった（間違って陰性と判定してしまった）もの」のことをFalse
Negative
(FN)といいます。従って、残りの4個はFNのカテゴリに含まれることになります。若干省略してTP
= 15、FN = 4のように表現すると、感度はTP / (TP +
FN)と定義できます。感度の値の取りうる範囲は0～1であり、1に近いほどよいです。
TP <- sum(obj) # obj中のTRUEの要素数をTPに格納
FN <- length(genus_truth) - TP # genus_truthの要素数 - TPをFNに格納
sensitivity <- TP / (TP + FN) # 感度を計算した結果をsensitivityに格納
sensitivity # sensitivityの中身## [1] 0.7894737陽性はdada2で検出された（系統割り当てできた）ユニークな属のことを指しますので、「陽性の数」は「dada2で検出された（系統割り当てできた）ユニークな属の数」、つまり222個になります。この222は、TPの15個とそれ以外（本当は検出されるべきではないにもかかわらず、検出してしまったもの）の207個に分割可能です。一般に、この「陽性でない（陰性）と判定すべきにもかかわらず、間違って陽性と判定したもの」のことをFalse Positive (FP)といいます。ここまでで、もう1つの評価指標である、TP / (TP + FP)と定義される適合率（Precision）を算出することができます。TP = 15、FP = 207ですので、属レベルの適合率は15 / (15 + 207) = 0.0675676になります。適合率の値の取りうる範囲は0～1であり、1に近いほどよいです。
FP <- length(table(genus)) - TP # genusのユニークな要素数 - TPをFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
precision # precisionの中身## [1] 0.06756757一般に、感度と適合率はトレードオフの関係にあります。感度を100%に近づけるには、（FPがどれだけ増えてもいいので）とにかくなんでも検出できるようにしてしまえばよいです。具体的には、ASVの絶対数を増やし、系統割り当て時の判定基準を緩めればよいのです。しかしそれはFPの増加を許容することと同義ですので、結果として適合率の低下に寄与します。
適合率を100%に近づけるには、FPが増えないように「陽性の数（系統割り当てできたユニークな属の数）」を絞ればよいです。具体的にはASVの配列数を減らし、系統割り当て時の判定基準を厳しくすればよいです。しかしそれはTPの減少を許容することと同義ですので、結果的に感度の低下に寄与します。この点を正しく認識して、1つの評価指標だけにこだわるのではなく、（感度と適合率のように）複数の指標を用いて多面的に評価することが重要です。
names属性は、ベクトルなどの各要素に名前を付与することです。以下にW3.2.2を再掲しますが、得られるobjという論理値ベクトルは、要素数が19個あります。
obj <- is.element(genus_truth, genus) # genus_truth中の各要素がgenusに含まれるか判定した結果をobjに格納
obj # objの中身## [1] TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE
## [13] TRUE TRUE FALSE TRUE TRUE TRUE TRUE この中には計4個のFALSEの要素がありますが、具体的な属の名前を知るには、「FALSEは(3,
6, 9, 15)番目の要素だからgenus_truth中の(3, 6, 9,
15)番目の要素と突き合わせて…」といった作業が必要になります。
genus_truth # genus_truthの中身## [1] "Acinetobacter" "Actinomyces" "Bacillus"
## [4] "Bacteroides" "Bifidobacterium" "Clostridium"
## [7] "Deinococcus" "Enterococcus" "Escherichia"
## [10] "Helicobacter" "Lactobacillus" "Listeria"
## [13] "Neisseria" "Porphyromonas" "Propionibacterium"
## [16] "Pseudomonas" "Rhodobacter" "Staphylococcus"
## [19] "Streptococcus" names属性の付与は、objというベクトルの何番目の要素がどの属名由来なのかがぱっと見でわかるように、要素に名前という属性を付与する作業になります。そうすることで、要素番号情報をイチイチ覚えなくてもよくなり、要素名（この場合は属の名前）で情報を抽出できるようになります。このあたりはread.table関数などを用いてファイルを読み込む際に、header=TRUEとかrow.names=1をつけるのと似たような思想になります。以下は、genus_truthをobjのnames属性として付与する例です。このような、ベクトルに対してnames属性を付与するテクニックはよく使われます。
names(obj) <- genus_truth # objにgenus_truthをnames属性として付与
obj # objの中身## Acinetobacter Actinomyces Bacillus Bacteroides
## TRUE TRUE FALSE TRUE
## Bifidobacterium Clostridium Deinococcus Enterococcus
## TRUE FALSE TRUE TRUE
## Escherichia Helicobacter Lactobacillus Listeria
## FALSE TRUE TRUE TRUE
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## TRUE TRUE FALSE TRUE
## Rhodobacter Staphylococcus Streptococcus
## TRUE TRUE TRUEこのような一目で分かる情報が得られると、例えば以下のような感じで、dada2実行結果ファイル（output_JSLAB24.txt）のGenus列に、“Actinomyces”が文字列検索で引っ掛かってくるであろう行番号の目安が得られたり、“Bacillus”が本当に文字列検索で引っ掛かってこない（integer(0)という結果に相当）ことの確認が可能になります。話の本筋からはそれていますが、このように得られた結果が本当に正しいかを都度検証する姿勢は、ミスを防ぐという意味で非常に重要です。
grep(pattern = "Actinomyces", x = genus) # genus中の"Actinomyces"のインデックス情報## [1] 11 90 204grep(pattern = "Bacillus", x = genus) # genus中の"Bacillus"のインデックス情報## integer(0) 種レベルの解析も同様であり、計22種のうち、どれを正しく検出できてどれを検出できなかったかを調べることができます。species_truthは、種の正解情報に相当する文字列ベクトルであり、要素数は22です。speciesは、種の予測結果に相当する文字列ベクトルであり、要素数は580です。
obj <- is.element(species_truth, species) # species_truth中の各要素がspeciesに含まれるか判定した結果をobjに格納
obj # objの中身## [1] FALSE TRUE FALSE TRUE FALSE TRUE TRUE FALSE FALSE TRUE FALSE FALSE
## [13] FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSEsum(obj) # obj中のTRUEの要素数## [1] 8この時点で、計22種のうち8種を正しく検出できたので（TP = 8）、種レベルの感度は8 / 22 = 0.3636364になります。感度はTP / (TP + FN)と定義でき、TP = 8、FN = 14になります。なお、分母の(TP + FN)は、「陽性と判定されるべき数（真の陽性数）」ですので、種の正解情報に相当する「ユニークな種の数」になります。
TP <- sum(obj) # obj中のTRUEの要素数をTPに格納
FN <- length(species_truth) - TP # species_truthの要素数 - TPをFNに格納
sensitivity <- TP / (TP + FN) # 感度を計算した結果をsensitivityに格納
sensitivity # sensitivityの中身## [1] 0.3636364陽性はdada2で検出された（系統割り当てできた）ユニークな種のことを指しますので、「陽性の数」は「dada2で検出された（系統割り当てできた）ユニークな種の数」、つまり54個になります。この54は、TPの8個とそれ以外（本当は検出されるべきではないにもかかわらず、検出してしまったもの）の46個（つまりFP = 46）に分割可能です。ここまでで、TP / (TP + FP)と定義される適合率（Precision）を算出することができます。種レベルの適合率は、8 / (8 + 46) = 0.1481481になります。
FP <- length(table(species)) - TP # speciesのユニークな要素数 - TPをFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
precision # precisionの中身## [1] 0.1481481 W3.2のスクリプトを以下にまとめます。ただし、属レベルの感度と種レベルの感度などを区別する必要があるので、最終的なオブジェクト名は変えています。これを実行するために必要な4つのオブジェクト（genus_truth,
genus, species_truth, and
species）は、W2.3を実行すれば得られます。
# 属レベルの感度
obj <- is.element(genus_truth, genus) # genus_truth中の各要素がgenusに含まれるか判定した結果をobjに格納
TP <- sum(obj) # obj中のTRUEの要素数をTPに格納
FN <- length(genus_truth) - TP # genus_truthの要素数 - TPをFNに格納
sensitivity <- TP / (TP + FN) # 感度を計算した結果をsensitivityに格納
# 属レベルの適合率
FP <- length(table(genus)) - TP # genusのユニークな要素数 - TPをFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
# オブジェクト名の変更
TP_genus <- TP # TPをTP_genusにコピー
FN_genus <- FN # FNをFN_genusにコピー
FP_genus <- FP # FPをFP_genusにコピー
sensitivity_genus <- sensitivity # sensitivityをsensitivity_genusにコピー
precision_genus <- precision # precisionをprecision_genusにコピー
# 種レベルの感度
obj <- is.element(species_truth, species) # species_truth中の各要素がspeciesに含まれるか判定した結果をobjに格納
TP <- sum(obj) # obj中のTRUEの要素数をTPに格納
FN <- length(species_truth) - TP # species_truthの要素数 - TPをFNに格納
sensitivity <- TP / (TP + FN) # 感度を計算した結果をsensitivityに格納
# 種レベルの適合率
FP <- length(table(species)) - TP # speciesのユニークな要素数 - TPをFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
# オブジェクト名の変更
TP_species <- TP # TPをTP_speciesにコピー
FN_species <- FN # FNをFN_speciesにコピー
FP_species <- FP # FPをFP_speciesにコピー
sensitivity_species <- sensitivity # sensitivityをsensitivity_speciesにコピー
precision_species <- precision # precisionをprecision_speciesにコピーここでは、重複なしの性能評価として、属レベルの感度と適合率、種レベルの感度と適合率を算出しました。ここでの「重複なし」は「ユニークな」と同義です。dada2実行結果ファイル（output_JSLAB24.txt）では、ASVの配列が違っていても同じ属や種に割り当てられることがあります。それらを重複カウントせずに、ユニークな属や種の数を調べた上で性能評価を行っています。性能評価指標を算出する際に何をどのように定義したかを正しく把握しておく必要があります。
W3.2では、dada2実行結果として得られた580個のASV配列のうち、222個の系統割り当てできたユニークな属の数と54個の系統割り当てできたユニークな種の数をベースとして性能評価を行いました。ここでは、343個の系統割り当てできた（重複を含む）属の数と59個の系統割り当てできた（重複を含む）種の数をベースとして性能評価を行います。なお、ここでは感度は算出せず、適合率のみ算出します。理由は、適合率はTP / (TP + FP) として計算されますが、その基礎情報であるTPとFPはともに重複ありで等価にカウント可能だからです。その一方で、感度はTP / (TP + FN) として計算されますが、FNは重複の有無にかかわらず属レベルでは常に4、種レベルでは常に14になり（W3.2.9）、原理的に重複（異なる複数のASV配列が同一の属に割り当てられうること）の影響を考慮できないからです。入力として用いる情報はW3.2と同じですが、念のため以下でも再掲します。
genus_truth：文字列ベクトル
length(genus_truth)：要素数（ユニークな属の数）は19species_truth：文字列ベクトル
length(species_truth)：要素数（ユニークな種の数）は22genus：文字列ベクトル
length(genus)：要素数（ASVの総数）は580length(na.omit(genus))：“NA”を除く要素数（系統割り当てできたASV数）は343length(table(genus))：“NA”を除くユニークな要素数（系統割り当てできたユニークな属の数）は222species：文字列ベクトル
length(species)：要素数（ASVの総数）は580length(na.omit(species))：“NA”を除く要素数（系統割り当てできたASV数）は59length(table(species))：“NA”を除くユニークな要素数（系統割り当てできたユニークな種の数）は54重複ありの場合は、以下のように書けば属レベルの適合率の値までを一気に得ることができます。
obj <- is.element(genus, genus_truth) # genus中の各要素がgenus_truthに含まれるか判定した結果をobjに格納
TP <- sum(obj) # obj中のTRUEの要素数をTPに格納
FP <- length(na.omit(genus)) - TP # genusのNAを除く要素数 - TPをFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
precision # precisionの中身## [1] 0.148688 W3.2.3との主な違いは、is.element関数の入力として与えている2つのベクトルの順序です。ここでは、580個の要素からなるgenus中の要素1つ1つに対して、それが19個の要素からなるgenus_truth中に存在するかどうかをTRUE
or
FALSEで判定しています。従って、objの要素数は580個になります。
length(obj) # objの要素数## [1] 580 objの中身を見ると、全体としてはFALSEがほとんどで、TRUEは最初のほうに偏って存在していることがわかります。この理由は、dada2実行結果ファイル（output_JSLAB24.txt）を眺めるとよくわかると思います。つまり、genus中の要素（ASV）の並びは、計12サンプルの出現頻度の総和が高い順になっているということです。どのサンプルでも検出されるようなASV配列は、本来検出したい既知微生物群集HM-280である確率が高いという直感ともよく一致します。
obj # objの中身## [1] TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE
## [25] TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE
## [37] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [85] FALSE FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [97] TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## [133] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [241] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [253] TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [277] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [289] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [301] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [313] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [325] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [337] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [349] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [373] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [385] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [397] FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## [409] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [421] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [433] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [445] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [457] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [469] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [481] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [493] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [505] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [517] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [541] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [553] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [565] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [577] FALSE FALSE FALSE FALSEこの中のTRUEの要素数が「陽性と判定されるべき属のうち、正しく陽性と判定できた属の数（True Positive; TP）」になります。つまり、TP = 51です。
TP # TPの中身## [1] 51注意点というわけではないですが、objの中で見えているFALSEに対応するASV配列の中には、そもそも系統割り当てできなかったものが(580 - 343) = 237個含まれます。これらの系統割り当てされなかったものは、dada2の予測結果は陰性扱いになります。したがって、「陰性と判定されるべき属のうち、間違って陽性と判定した属の数（False Positive; FP）」は、系統割り当てできたASV数（343個）からTP（51個）を引いた数（= 343 - 51 = 292）になります。つまり、FP = 292です。
FP # FPの中身## [1] 292TP / (TP + FP)で定義される適合率（precision）の値は、冒頭でも示していますが0.148688になります。
precision # precisionの中身## [1] 0.148688 ここは少しTips的な話になりますが、以下のようにすれば580個の要素からなるgenusオブジェクトの中から、objがTRUEの要素のみを得ることができます。これは、今回属レベル系統割り当て結果（NA含む）の計580個のASV配列のうち、既知微生物群集HM-280中に存在する計19種類の属名と一致するものを、重複ありでリストアップしていることと同義です。
genus[obj] # genus中のobjがTRUEの要素## [1] "Neisseria" "Enterococcus" "Pseudomonas" "Acinetobacter"
## [5] "Porphyromonas" "Staphylococcus" "Listeria" "Bacteroides"
## [9] "Actinomyces" "Rhodobacter" "Lactobacillus" "Streptococcus"
## [13] "Deinococcus" "Helicobacter" "Pseudomonas" "Bacteroides"
## [17] "Bacteroides" "Bifidobacterium" "Streptococcus" "Streptococcus"
## [21] "Pseudomonas" "Porphyromonas" "Streptococcus" "Pseudomonas"
## [25] "Streptococcus" "Neisseria" "Acinetobacter" "Streptococcus"
## [29] "Pseudomonas" "Neisseria" "Actinomyces" "Acinetobacter"
## [33] "Neisseria" "Neisseria" "Pseudomonas" "Pseudomonas"
## [37] "Enterococcus" "Bacteroides" "Bacteroides" "Neisseria"
## [41] "Actinomyces" "Streptococcus" "Porphyromonas" "Acinetobacter"
## [45] "Deinococcus" "Pseudomonas" "Bacteroides" "Acinetobacter"
## [49] "Helicobacter" "Bacteroides" "Porphyromonas" 従って、これをtable関数の入力として実行すると、属の種類ごとの出現頻度が結果として返されます。結果を眺めることで、例えばPseudomonasが最多の重複ASV数をもつ属だということがわかります。
table(genus[obj]) # genus[obj]中の要素の種類ごとの出現頻度##
## Acinetobacter Actinomyces Bacteroides Bifidobacterium Deinococcus
## 5 3 7 1 2
## Enterococcus Helicobacter Lactobacillus Listeria Neisseria
## 2 2 1 1 6
## Porphyromonas Pseudomonas Rhodobacter Staphylococcus Streptococcus
## 4 8 1 1 7 この結果を入力として、sort関数をdecreasing = Tオプションをつけて実行すれば、出現頻度の多い順（降順）にソートすることができます。
sort(table(genus[obj]), decreasing = T) # table(genus[obj])を降順にソート##
## Pseudomonas Bacteroides Streptococcus Neisseria Acinetobacter
## 8 7 7 6 5
## Porphyromonas Actinomyces Deinococcus Enterococcus Helicobacter
## 4 3 2 2 2
## Bifidobacterium Lactobacillus Listeria Rhodobacter Staphylococcus
## 1 1 1 1 1 decreasing = Tオプションをつけてsort関数を実行した時点で、出力結果の1番目の要素が出現頻度最大のものであることが確定しています。それを利用すれば、以下のように1番右側に[1]を追加するだけで、出現頻度が最大のもののみを抽出可能です。
sort(table(genus[obj]), decreasing = T)[1]# table(genus[obj])中の最多の要素## Pseudomonas
## 8 一連の結果を見ると、W3.2.5で解説したnames属性が最初からつけられていることがわかります。これを利用すれば、以下のように先ほどの結果を入力としてnames関数を実行すれば、出現頻度が最大の属名情報のみを抽出可能です。
names(sort(table(genus[obj]), decreasing = T)[1])# table(genus[obj])中の最多の要素のnames属性## [1] "Pseudomonas"which.max関数 最初から出現頻度が最大のものがターゲットの場合は、which.max関数も利用できます。which.max関数は、数値ベクトルを入力として与えると、最大値をもつ要素の位置（インデックス）を返します。これは数学記号のarg
max（またはargmax）に相当します。「maxならギリギリわかるがarg
max（またはargmax）という記号は見ただけでもうムリ！」というヒトも大勢いらっしゃるとは思います。しかし、which.max関数の入力情報と出力結果を見比べれば、「単にtable(genus[obj])の実行結果中の最大値は12番目の要素ですよ。」だと言っているだけということに気づけます。argmaxやargminはバイオインフォマティクス系の論文でよく出現します（例えばOsabe et al., BMC
Bioinformatics,
2021ではargminが出ます）ので、このような例を突破口として数学記号にも慣れていくとよいです。
table(genus[obj]) # genus[obj]中の要素の種類ごとの出現頻度##
## Acinetobacter Actinomyces Bacteroides Bifidobacterium Deinococcus
## 5 3 7 1 2
## Enterococcus Helicobacter Lactobacillus Listeria Neisseria
## 2 2 1 1 6
## Porphyromonas Pseudomonas Rhodobacter Staphylococcus Streptococcus
## 4 8 1 1 7which.max(table(genus[obj])) # table(genus[obj])中の最大値のインデックス## Pseudomonas
## 12names(which.max(table(genus[obj]))) # which.max(table(genus[obj]))のnames属性## [1] "Pseudomonas" 以下を実行すると、which.maxやwhich.min関数の利用上の注意点もわかると思います。これはwhich.minでの例ですが、実際には最小値の要素は複数存在するのですが、which.minは条件を満たす最初の要素のみを返します。ここで見えている4という数値は、条件を満たす最初の属であるBifidobacteriumの要素のインデックス（位置情報）になります。
table(genus[obj]) # genus[obj]中の要素の種類ごとの出現頻度##
## Acinetobacter Actinomyces Bacteroides Bifidobacterium Deinococcus
## 5 3 7 1 2
## Enterococcus Helicobacter Lactobacillus Listeria Neisseria
## 2 2 1 1 6
## Porphyromonas Pseudomonas Rhodobacter Staphylococcus Streptococcus
## 4 8 1 1 7which.min(table(genus[obj])) # table(genus[obj])中の最小値のインデックス## Bifidobacterium
## 4 この事実を理解した上で?which.minやhelp(which.min)などで関数のマニュアルを見ると、例えばDescription部分で「Determines
the location, i.e., index of the (first) minimum or
maximum of a numeric (or logical)
vector.」と書かれている意味がよくわかると思います。このあたりは通常は読み飛ばしてしまうものですが、いろいろと自分で動作確認しながら進めていく姿勢が大事です。
また、上の実行結果をみると、「which.minの実体は、sort結果の1番目の要素抽出だな」というにも気づけると思います。このようなwhich.minの内部で実際に動いている関数の挙動を想像できるようになれば、プログラム開発者の思考回路に近づいていけます。なお、ここで見えている1という数値は、条件を満たす最初の属であるBifidobacteriumの要素の出現頻度になります。
sort(table(genus[obj]))[1] # table(genus[obj])中の最小の要素## Bifidobacterium
## 1重複ありの場合は、以下のように書けば種レベルの適合率の値までを一気に得ることができます。
obj <- is.element(species, species_truth) # species中の各要素がspecies_truthに含まれるか判定した結果をobjに格納
TP <- sum(obj) # obj中のTRUEの要素数をTPに格納
FP <- length(na.omit(species)) - TP # speciesのNAを除く要素数 - TPをFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
precision # precisionの中身## [1] 0.1694915 W3.2.7との主な違いは、is.element関数の入力として与えている2つのベクトルの順序です。ここでは、580個の要素からなるspecies中の要素1つ1つに対して、それが22個の要素からなるspecies_truth中に存在するかどうかをTRUE
or
FALSEで判定しています。従って、objの要素数は580個になります。
length(obj) # objの要素数## [1] 580 objの中身を見ると、全体としてはFALSEがほとんどで、TRUEは最初のほうに偏って存在していることがわかります。この理由は、dada2実行結果ファイル（output_JSLAB24.txt）を眺めるとよくわかると思います。つまり、species中の要素（ASV）の並びは、計12サンプルの出現頻度の総和が高い順になっているということです。どのサンプルでも検出されるようなASV配列は、本来検出したい既知微生物群集HM-280である確率が高いという直感ともよく一致します。
obj # objの中身## [1] FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE TRUE FALSE
## [13] FALSE FALSE TRUE TRUE FALSE TRUE TRUE FALSE FALSE TRUE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [49] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [97] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [241] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [253] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [277] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [289] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [301] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [313] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [325] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [337] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [349] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [373] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [385] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [397] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [409] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [421] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [433] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [445] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [457] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [469] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [481] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [493] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [505] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [517] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [541] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [553] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [565] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [577] FALSE FALSE FALSE FALSEこの中のTRUEの要素数が「陽性と判定されるべき種のうち、正しく陽性と判定できた種の数（True Positive; TP）」になります。つまり、TP = 10です。
TP # TPの中身## [1] 10注意点というわけではないですが、objの中で見えているFALSEに対応するASV配列の中には、そもそも系統割り当てできなかったものが(580 - 59) = 521個含まれます。これらの系統割り当てされなかったものは、dada2の予測結果は陰性扱いになります。したがって、「陰性と判定されるべき種のうち、間違って陽性と判定した種の数（False Positive; FP）」は、系統割り当てできたASV数（59個）からTP（10個）を引いた数（= 59 - 10 = 49）になります。つまり、FP = 49です。
FP # FPの中身## [1] 49TP / (TP + FP)で定義される適合率（precision）の値は、冒頭でも示していますが0.1694915になります。
precision # precisionの中身## [1] 0.1694915なお、計10個のTPの内訳（どの種が何回割り当てられたか）は、下記の通りです。
sort(table(species[obj]), decreasing = T) # table(species[obj])を降順にソート##
## vulgatus aureus beijerinckii gingivalis mutans
## 3 1 1 1 1
## odontolyticus pylori radiodurans
## 1 1 1 W3.3のスクリプトを以下にまとめます。ただし、属レベルの適合率と種レベルの適合率を区別する必要があるので、最終的なオブジェクト名は変えています。これを実行するために必要な4つのオブジェクト（genus_truth,
genus, species_truth, and
species）は、W2.3を実行すれば得られます。
# 属レベルの適合率
obj <- is.element(genus, genus_truth) # genus中の各要素がgenus_truthに含まれるか判定した結果をobjに格納
TP <- sum(obj) # obj中のTRUEの要素数をTPに格納
FP <- length(na.omit(genus)) - TP # genusのNAを除く要素数 - TPをFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
# オブジェクト名の変更
TP_genus_dup <- TP # TPをTP_genus_dupにコピー
FP_genus_dup <- FP # FPをFP_genus_dupにコピー
precision_genus_dup <- precision # precisionをprecision_genus_dupにコピー
# 種レベルの適合率
obj <- is.element(species, species_truth) # species中の各要素がspecies_truthに含まれるか判定した結果をobjに格納
TP <- sum(obj) # obj中のTRUEの要素数をTPに格納
FP <- length(na.omit(species)) - TP # speciesのNAを除く要素数 - TPをFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
# オブジェクト名の変更
TP_species_dup <- TP # TPをTP_species_dupにコピー
FP_species_dup <- FP # FPをFP_species_dupにコピー
precision_species_dup <- precision # precisionをprecision_species_dupにコピーここでは、重複ありの性能評価として、属レベルと種レベルの適合率を算出しました。dada2実行結果ファイル（output_JSLAB24.txt）では、ASVの配列が違っていても同じ属や種に割り当てられることがあります。それらをそのまま重複ありでカウントして性能評価を行っています。
比較しやすいように、W3.2.9のまとめ情報も以下に再掲します。FNとTNの情報が存在しないため同じフォーマットにはしていませんが、重複ありのTPが確かに重複なしよりも多くなっており妥当だと判断できます。また、重複ありの適合率は、属レベル（0.148688 > 0.0675676）・種レベル（0.1694915 > 0.1481481）ともに重複なしよりも高くなっていることが分かります。
W3.2は同じ属または種単位で、そしてW3.3はASV配列単位で、既知微生物群集HM-280のものと一致するかどうかの評価を行いました。特に適合率は、属レベル、種レベル、重複の有無にかかわらず全て0.2未満であり、随分低い印象を受けます。これはこれで真実ではありますが、実際にはASV配列ごとに出現頻度情報があり、W3.3.1やW3.3.5のobjベクトル中のTRUEの要素が前半部分に偏っていたことからもわかりますが、出現頻度が高いASVはそれなりに検出できていました。ここでは、dada2実行結果ファイル（output_JSLAB24.txt）のSRR列に相当する計12サンプルの出現頻度の総和をASV配列ごとに算出し、出現頻度を含めて属レベルおよび種レベルの適合率を算出します。
ここで用いる入力データは下記の通りです。W3.2やW3.3と比べると、dataオブジェクトが増えただけになります。計5つのオブジェクトを用いますが、いずれもW2.3のコピペ実行で得られます。
data：数値行列
nrow(data)：行数（ASVの総数）は580ncol(data)：列数（サンプル数）は12genus：文字列ベクトル
length(genus)：要素数（ASVの総数）は580length(na.omit(genus))：“NA”を除く要素数（系統割り当てできたASV数）は343length(table(genus))：“NA”を除くユニークな要素数（系統割り当てできたユニークな属の数）は222species：文字列ベクトル
length(species)：要素数（ASVの総数）は580length(na.omit(species))：“NA”を除く要素数（系統割り当てできたASV数）は59length(table(species))：“NA”を除くユニークな要素数（系統割り当てできたユニークな種の数）は54genus_truth：文字列ベクトル
length(genus_truth)：要素数（ユニークな属の数）は19species_truth：文字列ベクトル
length(species_truth)：要素数（ユニークな種の数）は22rowSums関数 rowSums関数は、数値行列を入力として、行ごとの和を返します。
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納 計算手順的に、得られる結果であるfrequencyオブジェクトは、入力行列の行数と同じ数だけの要素からなる数値ベクトルになります。つまり、その要素数は580です。
length(frequency) # frequencyの要素数## [1] 580 frequencyオブジェクトの最初の3つの要素（n = 3）で確認します。dada2実行結果ファイル（output_JSLAB24.txt）のSRR列で手計算されると分かりますが、計算結果は確かに一致しています。
head(frequency, n = 3) # frequency中の最初の3つを表示## ATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCCGGGAACTGCGTTCTGAACTGGGTGACTCGAGTGTGTCAGAGGGAGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCTCCTGGGACAACACTGACGTTCATGCCCGA
## 685510
## TCCGGATTTATTGGGCGTAAAGCGAGCGCAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCCCGGCTCAACCGGGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAGAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGATATATGGAGGAACACCAGTGGCGAAGGCGGCTCTCTGGTCTGTAACTGACGCTGAGGCTCGA
## 493326
## ATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGA
## 484236 なお、rowSums関数は、W3.1.3で紹介したapply関数をMARGIN = 1およびFUN = sumで実行することと同義です。
hoge <- apply(X = data, MARGIN = 1, FUN = sum)# dataの各行に対して和を計算した結果をhogeに格納
head(hoge, n = 3) # hoge中の最初の3つを表示## ATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCCGGGAACTGCGTTCTGAACTGGGTGACTCGAGTGTGTCAGAGGGAGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCTCCTGGGACAACACTGACGTTCATGCCCGA
## 685510
## TCCGGATTTATTGGGCGTAAAGCGAGCGCAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCCCGGCTCAACCGGGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAGAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGATATATGGAGGAACACCAGTGGCGAAGGCGGCTCTCTGGTCTGTAACTGACGCTGAGGCTCGA
## 493326
## ATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGA
## 484236is.na関数 is.na関数は、ベクトルを入力として、要素がNAかどうかをTRUE
or
FALSEの論理値ベクトルとして返します。これを属の予測結果であるgenusオブジェクトを入力として実行すると、予測結果がNAとなった（dada2で系統割り当てができなかった）ASVに相当する要素をTRUE、そうでない要素をFALSEとした結果が返されます。
obj <- is.na(genus) # genusの各要素がNAかどうかの判定結果をobjに格納 従って、obj中のTRUEの要素数はgenus中の系統割り当てできなかったASV数に対応します。
sum(obj) # obj中のTRUEの要素数## [1] 237 obj中のFALSEの要素数は、genus中の系統割り当てできたASV数に対応します。この情報が欲しい場合は、論理値ベクトルのTRUEとFALSEを入れ替えてやればよいです。これを実現するには、オブジェクトの左側に!をつけてやればよいです。つまり、論理値ベクトルのTRUEとFALSEを交換したい場合は、論理演算子である!をつければよいです。!objの時点でTRUEとFALSEが入れ替わっています。
sum(!obj) # obj中のFALSEの要素数## [1] 343 論理演算子&は、要素数が同じ2つ論理値ベクトルを入力として、両方ともTRUEの要素のみをTRUE、それ以外をFALSEとして返します。これを利用すれば、以下の2つの論理値ベクトルを利用して、「系統割り当てはできたが、本来検出したい既知微生物群集HM-280には含まれなかったASVの要素」のみをTRUEとした論理値ベクトルobj_FPを作成することができます。
obj_P：dada2で属レベルの系統割り当てができた（NAではない）ASVをTRUE、それ以外をFALSEとした論理値ベクトルobj_TP：「属レベルの系統割り当て結果（genus）が本来検出したい既知微生物群集HM-280（genus_truth）に含まれるか否かを判定した論理値ベクトル obj_FPは、「属レベルの系統割り当てができた（NAではない）ASVのうち、本来検出したい既知微生物群集HM-280ではないものの位置情報」に相当します。この情報を得るには、obj_P & obj_TPではなくobj_P & !obj_TPとせねばならない点に注意が必要です。「…ではない」なので!obj_TPとせねばならないのです。
obj_P <- !is.na(genus) # genusの各要素がNAでないかどうかの判定結果をobj_Pに格納
obj_TP <- is.element(genus, genus_truth) # genus中の各要素がgenus_truthに含まれるかどうかの判定結果をobj_TPに格納
obj_FP <- obj_P & !obj_TP # obj_Pと!obj_TPの両方でTRUEの要素のみをTRUEとした結果をobj_FPに格納 obj_PのTRUEの要素数は343です。W3.4の冒頭の情報からもある程度想像がつきますが、sum(obj_P)とlength(na.omit(genus))は、予測結果の陽性数（dada2実行結果として得られた属の数）に相当します。
sum(obj_P) # obj_P中のTRUEの要素数## [1] 343 obj_TPのTRUEの要素数は51です。これは、W3.3.7の「属レベル（重複あり）」の中の、「陽性と判定されるべき属のうち、正しく陽性と判定できた属の数（True
Positive; TP）」に相当します。
sum(obj_TP) # obj_TP中のTRUEの要素数## [1] 51 obj_FPのTRUEの要素数は292です。これは、W3.3.7の「属レベル（重複あり）」の中の、「陰性と判定されるべき属のうち、間違って陽性と判定した属の数（False
Positive;
FP）」に相当します。この数自体は、sum(obj_P) - sum(obj_TP)
= 343 - 51からも簡単に算出できます。
sum(obj_FP) # obj_FP中のTRUEの要素数## [1] 292出現頻度を加味した属レベルの適合率は、以下のような感じで算出します。最初の4行は、おさらいです。
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
obj_P <- !is.na(genus) # genusの各要素がNAでないかどうかの判定結果をobj_Pに格納
obj_TP <- is.element(genus, genus_truth) # genus中の各要素がgenus_truthに含まれるかどうかの判定結果をobj_TPに格納
obj_FP <- obj_P & !obj_TP # obj_Pと!obj_TPの両方でTRUEの要素のみをTRUEとした結果をobj_FPに格納
TP <- sum(frequency[obj_TP]) # obj_TPがTRUEのfrequencyの総和を計算した結果をTPに格納
FP <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
precision # precisionの中身## [1] 0.8066153 出現頻度情報を加味した属レベルの適合率は0.8066153となり、W3.3で行った「出現頻度情報を加味しない、重複あり属レベルの適合率（0.148688）」と比べて非常に高いことが分かります。この理由は、「陽性と判定されるべき属のうち、正しく陽性と判定できた属」の出現頻度の総和（TP）が、「陰性と判定されるべき属のうち、間違って陽性と判定した属」の出現頻度の総和（FP）よりもかなり多いためです。
TP # TPの中身## [1] 2057950FP # FPの中身## [1] 493390 W3.3.7やW3.4.4でも示した「陽性と判定されるべき属のうち、正しく陽性と判定できた属」の数（sum(obj_TP)）と、「陰性と判定されるべき属のうち、間違って陽性と判定した属」の数（sum(obj_FP)）の大小関係と比較していただければ、出現頻度情報まで加味することの意味がよくわかると思います。
sum(obj_TP) # obj_TP中のTRUEの要素数## [1] 51sum(obj_FP) # obj_FP中のTRUEの要素数## [1] 292出現頻度を加味した種レベルの適合率は、以下のような感じで算出します。
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
obj_P <- !is.na(species) # speciesの各要素がNAでないかどうかの判定結果をobj_Pに格納
obj_TP <- is.element(species, species_truth) # species中の各要素がspecies_truthに含まれるかどうかの判定結果をobj_TPに格納
obj_FP <- obj_P & !obj_TP # obj_Pと!obj_TPの両方でTRUEの要素のみをTRUEとした結果をobj_FPに格納
TP <- sum(frequency[obj_TP]) # obj_TPがTRUEのfrequencyの総和を計算した結果をTPに格納
FP <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
precision # precisionの中身## [1] 0.2237251 出現頻度情報を加味した種レベルの適合率は0.2237251となり、W3.3で行った「出現頻度情報を加味しない、重複あり種レベルの適合率（0.1694915）」と比べて少し高い程度で、それほど変わらないことが分かります。「陽性と判定されるべき種のうち、正しく陽性と判定できた種」の出現頻度の総和（TP）と「陰性と判定されるべき属のうち、間違って陽性と判定した種」の出現頻度の総和（FP）は、下記の通りです。
TP # TPの中身## [1] 197835FP # FPの中身## [1] 686442 W3.3.7でも示した「陽性と判定されるべき種のうち、正しく陽性と判定できた種」の数（sum(obj_TP)）と、「陰性と判定されるべき種のうち、間違って陽性と判定した種」の数（sum(obj_FP)）は以下の通りでした。種レベルの場合は、出現頻度情報まで加味しても適合率の値はそれほど変わらない理由は、系統割り当てができた種の絶対数が少ないためといえます。
sum(obj_TP) # obj_TP中のTRUEの要素数## [1] 10sum(obj_FP) # obj_FP中のTRUEの要素数## [1] 49 W3.4のスクリプトを以下にまとめます。ただし、W3.2.8およびW3.3.6で得られるオブジェクトと区別する必要があるので、最終的なオブジェクト名は変えています。これを実行するために必要な5つのオブジェクト（data,
genus_truth, genus,
species_truth, and
species）は、W2.3を実行すれば得られます。W3.3.6とノリが若干違う点としては、属レベル・種レベルの適合率計算部分のスクリプトを統一すべく、一旦hogeおよびhoge_truthにしています。作業的に誤差範囲という考え方もあるかとは思いますが、このような積み重ねがオブジェクト名変更ミス発生率を下げます。
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
# 属レベルの適合率
hoge <- genus # genusをhogeにコピー
hoge_truth <- genus_truth # genus_truthをhoge_truthにコピー
obj_P <- !is.na(hoge) # hogeの各要素がNAでないかどうかの判定結果をobj_Pに格納
obj_TP <- is.element(hoge, hoge_truth) # hoge中の各要素がhoge_truthに含まれるかどうかの判定結果をobj_TPに格納
obj_FP <- obj_P & !obj_TP # obj_Pと!obj_TPの両方でTRUEの要素のみをTRUEとした結果をobj_FPに格納
TP <- sum(frequency[obj_TP]) # obj_TPがTRUEのfrequencyの総和を計算した結果をTPに格納
FP <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
# オブジェクト名の変更
TP_genus_freq <- TP # TPをTP_genus_freqにコピー
FP_genus_freq <- FP # FPをFP_genus_freqにコピー
precision_genus_freq <- precision # precisionをprecision_genus_freqにコピー
# 種レベルの適合率
hoge <- species # speciesをhogeにコピー
hoge_truth <- species_truth # species_truthをhoge_truthにコピー
obj_P <- !is.na(hoge) # hogeの各要素がNAでないかどうかの判定結果をobj_Pに格納
obj_TP <- is.element(hoge, hoge_truth) # hoge中の各要素がhoge_truthに含まれるかどうかの判定結果をobj_TPに格納
obj_FP <- obj_P & !obj_TP # obj_Pと!obj_TPの両方でTRUEの要素のみをTRUEとした結果をobj_FPに格納
TP <- sum(frequency[obj_TP]) # obj_TPがTRUEのfrequencyの総和を計算した結果をTPに格納
FP <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をFPに格納
precision <- TP / (TP + FP) # 適合率を計算した結果をprecisionに格納
# オブジェクト名の変更
TP_species_freq <- TP # TPをTP_species_freqにコピー
FP_species_freq <- FP # FPをFP_species_freqにコピー
precision_species_freq <- precision # precisionをprecision_species_freqにコピーここでは、出現頻度を加味した場合の性能評価として、属レベルと種レベルの適合率を算出しました。dada2実行結果ファイル（output_JSLAB24.txt）では、ASVの配列が違っていても同じ属や種に割り当てられることがあります。そして、個々のASV配列は出現頻度情報を持ちます。ここで行ったのは、属名がついたASV配列の総出現頻度（全陽性数 = TP + FP）のうち、計19属（22種）の既知微生物群集HM-280由来の総出現頻度（= TP）がどの程度含まれているのか（つまり適合率）を評価しています。
これまでの解析で、計12個のサンプル全体として、計19属から構成される既知微生物群集HM-280中15属がdada2実行結果として検出されています（W3.2.9）。また、系統割り当てによって343個のASV配列に対してNA以外の属名が付与され、そのうち51個が既知微生物群集由来のASVであることが分かっています（W3.3.7）。各ASV配列がもつ出現頻度情報を加味すると、既知微生物群集由来の51個のASVの総出現頻度（2057950）は、それ以外の292個のASVの総出現頻度（493390）よりも大きな値になることも分かっています（W3.4.8）。ここでは、計19属から構成される既知微生物群集HM-280中の個々の属の総出現頻度の情報を得ます。また、計19属以外の属はothersとして1つにまとめます。
rep関数 rep関数は、xで指定した値を、timesで指定した回数分だけ生成します。xで指定したhoge1は仮想の（dada2で検出された）属名、timesで指定したhoge2は仮想の出現頻度を表します。“B”を2か所に分けて記載しているのは、異なるASVで同じ属に割り当てられるものを念頭においています。
hoge1 <- c("B", "D", "A", "B") # 仮想属名情報をhoge1に格納
hoge2 <- c(2, 4, 3, 1) # 仮想出現頻度をhoge2に格納
freq <- rep(x = hoge1, times = hoge2) # hoge1をhoge2の指定回数分生成した結果をfreqに格納
freq # freqの中身## [1] "B" "B" "D" "D" "D" "D" "A" "A" "A" "B" freqを入力としてtable関数を実行すると、要素の種類ごとの出現頻度が返されます。“B”の出現頻度が意図通り3になっていることがわかります。
table(freq) # freq中の要素の種類ごとの出現頻度## freq
## A B D
## 3 3 4append関数 append関数は、ベクトルに要素を追加したり、ベクトル同士を連結する際に使われます。これは、出現頻度が0の属は0という結果を得たい場合に利用できるテクニックです。例えば、hoge3が出力結果に含めたい、計4種類の属名情報からなる文字列ベクトルだとします。append関数を用いて、実際の属ごとの出現頻度情報に相当するfreqベクトルとhoge3ベクトルを連結した結果をhoge4とします。これにtable関数を実行すると、得られる出現頻度情報は、本当の出現頻度よりも1多い結果となるものの、実際の出現頻度が0の”C”も1足したおかげで結果に含まれるようになります。
hoge3 <- c("A", "B", "C", "D") # 出力結果に含めたい属名情報をhoge3に格納
hoge4 <- append(freq, hoge3) # freqとhoge3を連結した結果をhoge4に格納
table(hoge4) # hoge4中の要素の種類ごとの出現頻度## hoge4
## A B C D
## 4 4 1 5 ここまでできれば、あとは以下のようにtable関数実行結果の各要素を全て-1してやれば、出現頻度が本当は0のものも含めつつ、本来の出現頻度として結果を得ることができます。
table(hoge4) - 1 # hoge4中の要素の種類ごとの出現頻度## hoge4
## A B C D
## 3 3 0 4これは、計12個のサンプルを合わせた結果でも計19属中15属しかdada2実行結果として検出できていないことから、サンプルごとの結果で属ごとの出現頻度を調べた場合に、出現頻度が0の属をどう取り扱うかが課題になることを念頭において、1つの解決策を提示したものになります。
本番に近い解析として、ここではW3.5.1に対応する実データ版を示します。最後のtable(freq)の結果からもわかりますが、計19属からなる本来検出したい既知微生物群集HM-280のうち、dada2で検出できた15属のみの出現頻度情報が得られるやり方になります。実行に必要な3つのオブジェクト（data,
genus, and
genus_truth）は、W2.3を実行すれば得られます。
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
obj_TP <- is.element(genus, genus_truth) # genus中の各要素がgenus_truthに含まれるかどうかの判定結果をobj_TPに格納
hoge1 <- genus[obj_TP] # genus[obj_TP]をhoge1にコピー
hoge2 <- frequency[obj_TP] # frequency[obj_TP]をhoge2にコピー
freq <- rep(x = hoge1, times = hoge2) # hoge1をhoge2の指定回数分生成した結果をfreqに格納
table(freq) # freq中の要素の種類ごとの出現頻度## freq
## Acinetobacter Actinomyces Bacteroides Bifidobacterium Deinococcus
## 178135 22229 44692 2399 11155
## Enterococcus Helicobacter Lactobacillus Listeria Neisseria
## 493338 8434 15355 35818 685642
## Porphyromonas Pseudomonas Rhodobacter Staphylococcus Streptococcus
## 109985 376440 19232 37246 17850これまでの復習にはなりますが、主なオブジェクトについて以下におさらいします：
data：数値行列
nrow(data)：行数（ASVの総数）は580ncol(data)：列数（サンプル数）は12frequency：数値ベクトル
rowSums関数を用いて、ASV配列ごとの（計12サンプルの）出現頻度の総和genus：文字列ベクトル
length(genus)：要素数（ASVの総数）は580length(na.omit(genus))：“NA”を除く要素数（系統割り当てできたASV数）は343length(table(genus))：“NA”を除くユニークな要素数（系統割り当てできたユニークな属の数）は222genus_truth：文字列ベクトル
length(genus_truth)：要素数（ユニークな属の数）は19hoge3に対応obj_TP：論理値ベクトル
genus中の要素1つ1つに対して、それがgenus_truthに含まれるかどうかをTRUE
or FALSEで判定した結果length(obj_TP)：要素数（ASVの総数）は580sum(obj_TP)：TRUEの要素数（本来検出したい既知微生物群集に含まれるASV数）は51genus[obj_TP]：文字列ベクトル
genusのうち、obj_TPがTRUEの要素のみを抽出した結果length(genus[obj_TP])：要素数は51hoge1に対応frequency[obj_TP]：数値ベクトル
frequencyのうち、obj_TPがTRUEの要素のみを抽出した結果length(frequency[obj_TP])：要素数は51hoge2に対応本番の解析として、ここではW3.5.2に対応する実データ版を示します。先ほどのW3.5.3のスクリプトに、最後の2行分を追加しただけです。
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
obj_TP <- is.element(genus, genus_truth) # genus中の各要素がgenus_truthに含まれるかどうかの判定結果をobj_TPに格納
hoge1 <- genus[obj_TP] # genus[obj_TP]をhoge1にコピー
hoge2 <- frequency[obj_TP] # frequency[obj_TP]をhoge2にコピー
freq <- rep(x = hoge1, times = hoge2) # hoge1をhoge2の指定回数分生成した結果をfreqに格納
hoge4 <- append(freq, genus_truth) # freqとgenus_truthを連結した結果をhoge4に格納
table(hoge4) - 1 # hoge4中の要素の種類ごとの出現頻度-1## hoge4
## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135 22229 0 44692
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399 0 11155 493338
## Escherichia Helicobacter Lactobacillus Listeria
## 0 8434 15355 35818
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642 109985 0 376440
## Rhodobacter Staphylococcus Streptococcus
## 19232 37246 17850以下は、既知微生物群集由来の計51配列の出現頻度の総和が、W3.4.8で見えている値（2057950）と同じであることを確認しているだけです。
sum(table(hoge4) - 1) # hoge4中の要素の種類ごとの出現頻度-1## [1] 2057950 W3.5の冒頭で述べたとおり、ここでは「既知微生物群集由来の19属の出現頻度」の数値ベクトルに、「既知微生物群集由来以外の292個のASVの総出現頻度（493390）」を、othersとして追加します。この値は、W3.4.5のFPと同じものになります。従って、ここで示すスクリプトは、先ほどの「W3.5.4で得た計19属分の出現頻度情報を格納した数値ベクトル」に対して、append関数を用いて「W3.4.5のFPに相当する数値」を追加する内容になります。
# W3.5.4
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
obj_TP <- is.element(genus, genus_truth) # genus中の各要素がgenus_truthに含まれるかどうかの判定結果をobj_TPに格納
hoge1 <- genus[obj_TP] # genus[obj_TP]をhoge1にコピー
hoge2 <- frequency[obj_TP] # frequency[obj_TP]をhoge2にコピー
freq <- rep(x = hoge1, times = hoge2) # hoge1をhoge2の指定回数分生成した結果をfreqに格納
hoge4 <- append(freq, genus_truth) # freqとgenus_truthを連結した結果をhoge4に格納
hoge5 <- table(hoge4) - 1 # hoge4中の要素の種類ごとの出現頻度-1をhoge5に格納
# W3.4.5の一部
obj_P <- !is.na(genus) # genusの各要素がNAでないかどうかの判定結果をobj_Pに格納
obj_FP <- obj_P & !obj_TP # obj_Pと!obj_TPの両方でTRUEの要素のみをTRUEとした結果をobj_FPに格納
others <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をothersに格納
names(others) <- "Others" # othersに"Others"をnames属性として付与
# 連結
hoge6 <- append(hoge5, others) # hoge5とothersを連結した結果をhoge6に格納
hoge6 # hoge6の中身## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135 22229 0 44692
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399 0 11155 493338
## Escherichia Helicobacter Lactobacillus Listeria
## 0 8434 15355 35818
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642 109985 0 376440
## Rhodobacter Staphylococcus Streptococcus Others
## 19232 37246 17850 493390以下は、dada2実行結果として得られた計580配列のうち、系統割り当てできた（NAとならなかった）343配列の出現頻度の総和が、W3.4.8で見えている値（2551340）と同じことを確認しているだけです。
sum(hoge6)## [1] 2551340 さきほどまでの結果は、dada2実行結果として得られた計580配列のうち、系統割り当てできた（NAとならなかった）343配列のみの出現頻度をカウントしました。ここでは、系統割り当てできなかった（NAとなった）ものも全て含めて出現頻度情報を得るやり方を示します。W3.5.4に相当する部分は不変で、W3.4.5のobj_FPを得る部分のみが異なっています。
# W3.5.4
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
obj_TP <- is.element(genus, genus_truth) # genus中の各要素がgenus_truthに含まれるかどうかの判定結果をobj_TPに格納
hoge1 <- genus[obj_TP] # genus[obj_TP]をhoge1にコピー
hoge2 <- frequency[obj_TP] # frequency[obj_TP]をhoge2にコピー
freq <- rep(x = hoge1, times = hoge2) # hoge1をhoge2の指定回数分生成した結果をfreqに格納
hoge4 <- append(freq, genus_truth) # freqとgenus_truthを連結した結果をhoge4に格納
hoge5 <- table(hoge4) - 1 # hoge4中の要素の種類ごとの出現頻度-1をhoge5に格納
# W3.4.5の一部
obj_FP <- !obj_TP # obj_TPのTRUEとFALSEを入れ替えた結果をobj_FPに格納
others <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をothersに格納
names(others) <- "Others" # othersに"Others"をnames属性として付与
# 連結
hoge6 <- append(hoge5, others) # hoge5とothersを連結した結果をhoge6に格納
hoge6 # hoge6の中身## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135 22229 0 44692
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399 0 11155 493338
## Escherichia Helicobacter Lactobacillus Listeria
## 0 8434 15355 35818
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642 109985 0 376440
## Rhodobacter Staphylococcus Streptococcus Others
## 19232 37246 17850 786825W3では、サンプル全体の解析を行っています。W4ではサンプルごとの解析を行いますが、計12個あるサンプルごとにW3.5.6のような解析を都度実行するのは面倒です。自作関数についてはW3.1.2でも解説していますが、ここではW3.5.6のスクリプトに酷似した自作関数を作成します。
get_W3.5.6関数 W3.5.6実行時の入力は、3つのオブジェクト（data,
genus, and
genus_truth）でした。したがって、ここで作成する自作関数get_W3.5.6の入力もこの3つを指定するように設計しています。それが1行目の記述内容に相当します。また、W3.5.6の主な出力はhoge6でしたので、自作関数get_W3.5.6もhoge6の中身を返すように設計しています。それが最終行付近のreturn(hoge6)に相当します。
get_W3.5.6 <- function(data, genus, genus_truth){
# W3.5.4
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
obj_TP <- is.element(genus, genus_truth)# genus中の各要素がgenus_truthに含まれるかどうかの判定結果をobj_TPに格納
hoge1 <- genus[obj_TP] # genus[obj_TP]をhoge1にコピー
hoge2 <- frequency[obj_TP] # frequency[obj_TP]をhoge2にコピー
freq <- rep(x = hoge1, times = hoge2) # hoge1をhoge2の指定回数分生成した結果をfreqに格納
hoge4 <- append(freq, genus_truth) # freqとgenus_truthを連結した結果をhoge4に格納
hoge5 <- table(hoge4) - 1 # hoge4中の要素の種類ごとの出現頻度-1をhoge5に格納
# W3.4.5
obj_FP <- !obj_TP # obj_TPのTRUEとFALSEを入れ替えた結果をobj_FPに格納
others <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をothersに格納
names(others) <- "Others" # othersに"Others"をnames属性として付与
# Concatenate
hoge6 <- append(hoge5, others) # hoge5とothersを連結した結果をhoge6に格納
return(hoge6) # hoge6の中身を関数実行結果として返す
} ここで入力として想定している3つのオブジェクト（data,
genus, and
genus_truth）について、復習を兼ねて以下で解説しておきます。
data：数値行列
nrow(data)：行数（ASVの総数）は580ncol(data)：列数（サンプル数）は12genus：文字列ベクトル
length(genus)：要素数（ASVの総数）は580length(na.omit(genus))：“NA”を除く要素数（系統割り当てできたASV数）は343length(table(genus))：“NA”を除くユニークな要素数（系統割り当てできたユニークな属の数）は222genus_truth：文字列ベクトル
length(genus_truth)：要素数（ユニークな属の数）は19get_W3.5.6の実行 get_W3.5.6関数の実行は、以下のような順番で3つのオブジェクトを入力として与えればよいです。ここでは、出力結果をresultオブジェクトに一旦保存してから中身を確認しています。出力結果はW3.5.6と同じですので、確かにこの関数がうまく機能していることが分かると思います。
result <- get_W3.5.6(data, genus, genus_truth) # 属の出現頻度のまとめ情報を取得してresultに格納
result # resultの中身## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135 22229 0 44692
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399 0 11155 493338
## Escherichia Helicobacter Lactobacillus Listeria
## 0 8434 15355 35818
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642 109985 0 376440
## Rhodobacter Staphylococcus Streptococcus Others
## 19232 37246 17850 786825 get_W3.5.6関数は、一見すると属レベルの情報しか受け付けないような印象を受けます。しかし実際には、フォーマットが同じ「種の予測結果に相当するspecies」と「種の正解情報に相当するspecies_truth」を入力として実行することもできます。
result <- get_W3.5.6(data, species, species_truth) # 種の出現頻度のまとめ情報を取得してresultに格納
result # resultの中身## acnes adolescentis aeruginosa agalactiae aureus
## 0 0 0 0 47
## baumannii beijerinckii cereus coli epidermidis
## 0 76 0 0 0
## faecalis gasseri gingivalis meningitides monocytogenes
## 0 0 109189 0 0
## mutans odontolyticus pneumoniae pylori radiodurans
## 2081 22200 0 8430 11149
## sphaeroides vulgatus Others
## 0 44663 2646940 出力結果の出現頻度情報をまとめた数値ベクトルは、計22個の種と1個のothersから構成されます。また、「陽性と判定されるべき種のうち、正しく陽性と判定できた種の数（True
Positive;
TP）」が8であることは既知です（W3.2.6およびW3.2.9）。実際、aureus,
beijerinckii, gingivalis, mutans, odontolyticus, pylori, radiodurans,
vulgatusの8つが0より大きい出現頻度をもつことが確認できます。get_W3.5.6関数が属の情報だけでなく種の情報も入力として受け付け、正しい結果を返すことができる理由は、W3.6.1で定義したこの関数が内部的に利用している、入力として受け付ける引数（ひきすう）の名前が、たまたまりgenusやgenus_truthと属名を連想させるものになっているだけだからです。なお、ここで引数の名前といっているのは、W3.6.1のget_W3.6.1関数定義の先頭行（get_W3.5.6 <- function(data, genus, genus_truth){）で記載しているdata,
genus, genus_truthのことです。
get_W3.5.6_2関数ここでは、「引数の名前」がgenusやgenus_truthのような名前である必要はないということを示すべく、先頭行（1行目）をpredictionおよびtruthに変更しています。これに合わせるため、以下の3～8行目のgenusおよびgenus_truthについても、それぞれpredictionおよびtruthに変更しています。
get_W3.5.6_2 <- function(data, prediction, truth){
# W3.5.4
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
obj_TP <- is.element(prediction, truth) # prediction中の各要素がtruthに含まれるかどうかの判定結果をobj_TPに格納
hoge1 <- prediction[obj_TP] # prediction[obj_TP]をhoge1にコピー
hoge2 <- frequency[obj_TP] # frequency[obj_TP]をhoge2にコピー
freq <- rep(x = hoge1, times = hoge2) # hoge1をhoge2の指定回数分生成した結果をfreqに格納
hoge4 <- append(freq, truth) # freqとtruthを連結した結果をhoge4に格納
hoge5 <- table(hoge4) - 1 # hoge4中の要素の種類ごとの出現頻度-1をhoge5に格納
# W3.4.5
obj_FP <- !obj_TP # obj_TPのTRUEとFALSEを入れ替えた結果をobj_FPに格納
others <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をothersに格納
names(others) <- "Others" # othersに"Others"をnames属性として付与
# Concatenate
hoge6 <- append(hoge5, others) # hoge5とothersを連結した結果をhoge6に格納
return(hoge6) # hoge6の中身を関数実行結果として返す
}get_W3.5.6_2の実行 引数の名前がpredictionおよびtruthに変更されていても、属の場合はgenusおよびgenus_truthに、そして種の場合はspeciesおよびspecies_truthを入力として与えると、適切に処理されることを示します。
get_W3.5.6_2(data, genus, genus_truth) # 属の出現頻度のまとめ情報を取得## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135 22229 0 44692
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399 0 11155 493338
## Escherichia Helicobacter Lactobacillus Listeria
## 0 8434 15355 35818
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642 109985 0 376440
## Rhodobacter Staphylococcus Streptococcus Others
## 19232 37246 17850 786825get_W3.5.6_2(data, species, species_truth)# 種の出現頻度のまとめ情報を取得## acnes adolescentis aeruginosa agalactiae aureus
## 0 0 0 0 47
## baumannii beijerinckii cereus coli epidermidis
## 0 76 0 0 0
## faecalis gasseri gingivalis meningitides monocytogenes
## 0 0 109189 0 0
## mutans odontolyticus pneumoniae pylori radiodurans
## 2081 22200 0 8430 11149
## sphaeroides vulgatus Others
## 0 44663 2646940確かに属の結果はW3.6.2と、そして種の結果はW3.6.3と同じであることがわかります。
get_W3.5.6_2やget_W3.5.6_2関数は、3種類の入力情報（出現頻度情報を含む数値行列、dada2予測結果の文字列ベクトル、既知微生物群集情報からなる文字列ベクトル）を与えて、既知微生物群集中の属または種ごと（とそれ以外をまとめた）の出現頻度をhoge6オブジェクトにまとめ、それを返す役割を持たせたものでした。しかし、hoge6だけではなく、obj_TPのような予測結果のどの要素が既知微生物群集と合致しているかという論理値ベクトル情報も同時に得たいヒトもいるでしょう。ここでは、出力結果としてhoge6とobj_TPの複数の情報をリスト形式で返すget_W3.5.6_3関数を作成します。
get_W3.5.6_3 <- function(data, prediction, truth){
# W3.5.4
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
obj_TP <- is.element(prediction, truth) # prediction中の各要素がtruthに含まれるかどうかの判定結果をobj_TPに格納
hoge1 <- prediction[obj_TP] # prediction[obj_TP]をhoge1にコピー
hoge2 <- frequency[obj_TP] # frequency[obj_TP]をhoge2にコピー
freq <- rep(x = hoge1, times = hoge2) # hoge1をhoge2の指定回数分生成した結果をfreqに格納
hoge4 <- append(freq, truth) # freqとtruthを連結した結果をhoge4に格納
hoge5 <- table(hoge4) - 1 # hoge4中の要素の種類ごとの出現頻度-1をhoge5に格納
# W3.4.5
obj_FP <- !obj_TP # obj_TPのTRUEとFALSEを入れ替えた結果をobj_FPに格納
others <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をothersに格納
names(others) <- "Others" # othersに"Others"をnames属性として付与
# Concatenate
hoge6 <- append(hoge5, others) # hoge5とothersを連結した結果をhoge6に格納
# output
rtn <- list() # 結果をリスト形式で返すためのrtnを作成
rtn$freque <- hoge6 # hoge6をfrequeという名前で追加
rtn$TPposi <- obj_TP # obj_TPをTPposiという名前で追加
return(rtn) # rtnの中身を関数実行結果として返す
} get_W3.5.6_2関数との違いは、最後の5行分のみです。hoge6とobj_TPは同じベクトルではありますが、要素数は異なります。それゆえ、ここでは要素数が異なる場合でも情報を格納できるリスト形式にしています。rtn <-list()は、「listという関数を用いて中身が空のリスト形式のオブジェクトを作成し、それをrtnという名前で保存しているのだ。」と理解すればよいです。以降は後で取り出したい情報を、どのような名前で格納するかを併せて指定して、順次rtnに格納していくだけです。ここでは、hoge6の情報をfrequeという名前で格納しています。また、obj_TPの情報をTPposiという名前で格納しています。rtnというのは、リターン（return）するものというのを念頭においています。W3.6.4で解説したことと関連しますが、関数内部で用いているオブジェクト名は基本的になんでもよいいので、必ずしもrtnである必要はありません。但し、frequeやTPposi`は関数実行後も継承されますのでご注意ください。
get_W3.5.6_3の実行 get_W3.5.6_3関数の実行は、get_W3.5.6やget_W3.5.6_2と同じです。ここでは、属の情報（genusとgenus_truth）を入力として与えて実行し、結果をresというオブジェクトに格納しています。
res <- get_W3.5.6_3(data, genus, genus_truth) # 属の出現頻度のまとめ情報を取得してresに格納
res # resの中身## $freque
## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135 22229 0 44692
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399 0 11155 493338
## Escherichia Helicobacter Lactobacillus Listeria
## 0 8434 15355 35818
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642 109985 0 376440
## Rhodobacter Staphylococcus Streptococcus Others
## 19232 37246 17850 786825
##
## $TPposi
## [1] TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE
## [25] TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE
## [37] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [85] FALSE FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [97] TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## [133] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [241] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [253] TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [277] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [289] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [301] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [313] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [325] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [337] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [349] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [373] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [385] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [397] FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## [409] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [421] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [433] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [445] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [457] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [469] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [481] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [493] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [505] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [517] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [541] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [553] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [565] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [577] FALSE FALSE FALSE FALSE resの出力結果を眺めると、最初にhoge6の中身に相当する「既知微生物群集中の属とそれ以外をまとめた出現頻度」情報が$frequeという部分に、そしてobj_TPの中身に相当する「予測結果のどの要素が既知微生物群集と合致しているかという論理値ベクトル」情報が$TPposiという部分に格納されていることが分かります。これを眺めれば、前述の「但し、frequeやTPposiは関数実行後も継承されますのでご注意ください。」の意味がよくわかると思います。
str関数 str関数は、入力として与えられたオブジェクトの構造（structure）を返します。通常は、ここの表示結果を頼りにして、どのようにして欲しい情報を抽出するかを判断します。
str(res) # resの構造を表示## List of 2
## $ freque: Named num [1:20] 178135 22229 0 44692 2399 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ TPposi: logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ... 例えば、frequeという部分の情報を得たい場合は、res$frequeのようにresとfrequeの間に$を挟んだような形で指定します。出力結果は、get_W3.5.6_2関数の実行結果と同じであることが分かります（W3.6.5）。また、get_W3.5.6_3関数の関数名の由来でもあるW3.5.6の出力結果とも同じです。
res$freque # res$frequeの中身## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135 22229 0 44692
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399 0 11155 493338
## Escherichia Helicobacter Lactobacillus Listeria
## 0 8434 15355 35818
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642 109985 0 376440
## Rhodobacter Staphylococcus Streptococcus Others
## 19232 37246 17850 786825 同様に、TPposiという部分の情報を得たい場合は、res$TPposiのようにresとTPposiの間に$を挟んだような形で指定します。
res$TPposi # res$TPposiの中身## [1] TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE
## [25] TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE
## [37] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [85] FALSE FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [97] TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## [133] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [241] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [253] TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [277] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [289] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [301] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [313] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [325] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [337] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [349] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [373] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [385] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [397] FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## [409] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [421] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [433] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [445] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [457] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [469] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [481] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [493] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [505] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [517] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [541] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [553] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [565] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [577] FALSE FALSE FALSE FALSE 論理値ベクトルであるres$TPposi中のTRUEの要素数は、51です。この数値は、W3.3.1で見えているTPの要素数に相当します。つまり、「陽性（既知微生物群集）と判定されるべき属のうち、正しく陽性と判定できた属の数（True
Positive;
TP）」です。このような情報が得られるもの、「予測結果のどの要素が既知微生物群集と合致しているかという論理値ベクトル」がres$TPposiとして利用可能だからです。
sum(res$TPposi) # res$TPposi中のTRUEの要素数## [1] 51get_W3.6関数 最後に、既知微生物群集中の属または種をどれだけの種類検出できたかという感度（W3.2.2やW3.2.6）、そして予測結果のうち既知微生物群集由来のものがどれだけの割合含まれるかという適合率（W3.4.5やW3.4.6）の情報も出力結果として追加します。関数名は、W3.6の集大成という意味でget_W3.6とします。
get_W3.6 <- function(data, prediction, truth){
# W3.5.4
frequency <- rowSums(data) # dataの行ごとの総和をfrequencyに格納
obj_TP <- is.element(prediction, truth) # prediction中の各要素がtruthに含まれるかどうかの判定結果をobj_TPに格納
hoge1 <- prediction[obj_TP] # prediction[obj_TP]をhoge1にコピー
hoge2 <- frequency[obj_TP] # frequency[obj_TP]をhoge2にコピー
freq <- rep(x = hoge1, times = hoge2) # hoge1をhoge2の指定回数分生成した結果をfreqに格納
hoge4 <- append(freq, truth) # freqとtruthを連結した結果をhoge4に格納
hoge5 <- table(hoge4) - 1 # hoge4中の要素の種類ごとの出現頻度-1をhoge5に格納
# W3.4.5
obj_FP <- !obj_TP # obj_TPのTRUEとFALSEを入れ替えた結果をobj_FPに格納
others <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をothersに格納
names(others) <- "Others" # othersに"Others"をnames属性として付与
# Concatenate
hoge6 <- append(hoge5, others) # hoge5とothersを連結した結果をhoge6に格納
# sensitivity
obj <- is.element(truth, prediction) # truth中の各要素がpredictionに含まれるかどうかの判定結果をobjに格納
sensi <- sum(obj)/length(truth) # 感度の計算結果をsensiに格納
# precision（FPにNA含む）
TP <- sum(frequency[obj_TP]) # obj_TPがTRUEのfrequencyの総和を計算した結果をTPに格納
FP <- others # othersをFPにコピー
preci1 <- TP / (TP + FP) # 適合率を計算した結果をpreci1に格納
# precision（FPにNA含まない）
obj_P <- !is.na(prediction) # predictionの各要素がNAでないかどうかの判定結果をobj_Pに格納
obj_FP <- obj_P & !obj_TP # obj_Pと!obj_TPの両方でTRUEの要素のみをTRUEとした結果をobj_FPに格納
FP <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をFPに格納
preci2 <- TP / (TP + FP) # 適合率を計算した結果をpreci2に格納
# output
rtn <- list() # 結果をリスト形式で返すためのrtnを作成
rtn$freque <- hoge6 # hoge6をfrequeという名前で追加
rtn$TPposi <- obj_TP # obj_TPをTPposiという名前で追加
rtn$sensitivity <- sensi # sensiをsensivitityという名前で追加
rtn$precision_with_NA <- preci1 # preci1をprecision_with_NAという名前で追加
rtn$precision_without_NA <- preci2 # preci2をprecision_without_NAという名前で追加
return(rtn) # rtnの中身を関数実行結果として返す
} TP / (TP +
FP)で定義される適合率は、FPを「TP以外の全てのASV数（NAも含める）」のようにNAを含める場合と、「TP以外の全てのASV数（但し系統割り当てできなかったNAは含めない）」の2つパターンが可能です。それゆえ、ここでは前者をpreci1として計算してprecision_with_NAで出力結果から取り出せるように、そして後者をpreci2として計算してprecision_without_NAで出力結果から取り出せるようにしています。
get_W3.6の実行 get_W3.6関数を実行し、結果をresというオブジェクトに格納しています。
res <- get_W3.6(data, genus, genus_truth) # 属の出現頻度のまとめ情報を取得してresに格納
res # resの中身## $freque
## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135 22229 0 44692
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399 0 11155 493338
## Escherichia Helicobacter Lactobacillus Listeria
## 0 8434 15355 35818
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642 109985 0 376440
## Rhodobacter Staphylococcus Streptococcus Others
## 19232 37246 17850 786825
##
## $TPposi
## [1] TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE
## [25] TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE
## [37] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [85] FALSE FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [97] TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## [133] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [241] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [253] TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [277] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [289] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [301] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [313] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [325] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [337] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [349] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [373] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [385] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [397] FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## [409] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [421] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [433] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [445] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [457] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [469] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [481] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [493] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [505] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [517] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [541] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [553] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [565] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [577] FALSE FALSE FALSE FALSE
##
## $sensitivity
## [1] 0.7894737
##
## $precision_with_NA
## Others
## 0.723414
##
## $precision_without_NA
## [1] 0.8066153 感度の値（=
0.7894737）は、res$sensitivityで抽出できます。W3.2.2で見えているものと同じであることがわかります。
res$sensitivity # res$sensitivityの中身## [1] 0.7894737 NA由来配列の出現頻度を含めて算出した適合率の値（=
0.723414）は、res$precision_with_NAで抽出できます。
res$precision_with_NA # res$precision_with_NAの中身## Others
## 0.723414 NA由来配列の出現頻度を含めずに算出した適合率の値（=
0.8066153）は、res$precision_without_NAで抽出できます。0.723414のres$precision_with_NAよりも値が大きいのは、至極当然です。理由は、NA由来配列の出現頻度を含めなければFPの数が少なくなること、そして適合率はTP
/ (TP + FP)で定義されるからです。
res$precision_without_NA # res$precision_without_NAの中身## [1] 0.8066153source関数 source関数は、入力として指定したファイルの中身を実行します。Linuxが分かるヒト向けの説明としては、シェルスクリプトのようなものになります。ご利益としては、これまでgetNonnaNum、get_W3.5.6、get_W3.5.6_2、get_W3.5.6_3、get_W3.6と5種類の自作関数を示してきました。これらの関数は、一度コンソール（Console）画面上で実行して定義しておくと、ソフトウェアを閉じない限り使えます。しかしRを再起動すると、都度関数を再読み込みせねばなりません。それゆえ、複数の自作関数を定義したファイル（JSLAB_source.R）を予め用意しておき、それを最初に一度にまとめて読み込むことがよく行われます。
source("http://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB_source.R") # 関数の読込このファイルの記述内容を見ていただければわかりますが、関数を定義しているだけなので見た目上何も変化はありません。エラーが出なければうまくいったと思っていただいてかまいません。
sourceのご利益source関数を用いた関数定義ファイルの読込のご利益を実感する1つの方法は、以下の3つの手順を実行し、正しく
「get_W3.6(data, genus, genus_truth) でエラー: 関数 “get_W3.6”
を見つけることができませんでした」
というエラーに遭遇することです。
# W2.1 dada2実行結果ファイルのほう
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t")# in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
obj <- grep(pattern = "SRR", x = colnames(hoge)) # colnames(hoge)の中から"SRR"という文字を含む要素番号を抽出した結果をobjに格納
data <- hoge[, obj] # hoge行列中のobjで指定した列のみ抽出した結果をdataに格納
genus <- hoge$Genus # hoge中のGenus列を抽出した結果をgenusに格納
species <- hoge$Species # hoge中のSpecies列を抽出した結果をspeciesに格納
# W2.2 菌株情報ファイルのほう
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB23_table1.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t") #in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
organism <- hoge$Organism # hoge中のOrganism列を抽出した結果をorganismに格納
hoge <- strsplit(x = organism, split = ",")# organism中の文字列ベクトルをコンマで分割した結果をhogeに格納
hoge2 <- sapply(X = hoge, FUN = `[`, 1) # hoge中のベクトル形式の各要素から、1番目の要素の抽出結果をhoge2に格納
hoge3 <- strsplit(x = hoge2, split = " ") # hoge2中の文字列ベクトルを半角スペースで分割した結果をhoge3に格納
genus_truth <- sapply(X = hoge3, FUN = `[`, 1) # hoge3中のベクトル形式の各要素から、1番目の要素の抽出結果をgenus_truthに格納
species_truth <- sapply(X = hoge3, FUN = `[`, 2)# hoge3中のベクトル形式の各要素から、2番目の要素の抽出結果をspecies_truthに格納
genus_truth <- unique(genus_truth) # genus_truth中のユニークな要素をgenus_truthに格納
species_truth <- unique(species_truth) # species_truth中のユニークな要素をspecies_truthに格納get_W3.6(data, genus, genus_truth) # 属の出現頻度のまとめ情報を取得 無事「get_W3.6(data, genus, genus_truth) でエラー: 関数 “get_W3.6”
を見つけることができませんでした」というエラーに遭遇できたと思います。この状態で、以下を実行すればsource関数のご利益を実感できると思います。
source("http://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB_source.R") # 関数の読込
get_W3.6(data, genus, genus_truth) # 属の出現頻度のまとめ情報を取得## $freque
## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135 22229 0 44692
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399 0 11155 493338
## Escherichia Helicobacter Lactobacillus Listeria
## 0 8434 15355 35818
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642 109985 0 376440
## Rhodobacter Staphylococcus Streptococcus Others
## 19232 37246 17850 786825
##
## $TPposi
## [1] TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE
## [25] TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE
## [37] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [85] FALSE FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [97] TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## [133] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [241] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [253] TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [277] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [289] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [301] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [313] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [325] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [337] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [349] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [373] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [385] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [397] FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## [409] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [421] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [433] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [445] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [457] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [469] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [481] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [493] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [505] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [517] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [541] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [553] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [565] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [577] FALSE FALSE FALSE FALSE
##
## $sensitivity
## [1] 0.7894737
##
## $precision_with_NA
## Others
## 0.723414
##
## $precision_without_NA
## [1] 0.8066153W3.4以降では、dada2実行結果ファイル（output_JSLAB24.txt）のSRR列に相当する計12サンプルの出現頻度の総和をASV配列ごとに算出し、出現頻度を含めた性能評価を行いました。ここでは、サンプルごとに同様の解析を行います。ここで用いる入力データは下記の通りです。計5つのオブジェクトを用いますが、いずれもW2.3のコピペ実行で得られます。
data：数値行列
nrow(data)：行数（ASVの総数）は580ncol(data)：列数（サンプル数）は12genus：文字列ベクトル
length(genus)：要素数（ASVの総数）は580length(na.omit(genus))：“NA”を除く要素数（系統割り当てできたASV数）は343length(table(genus))：“NA”を除くユニークな要素数（系統割り当てできたユニークな属の数）は222species：文字列ベクトル
length(species)：要素数（ASVの総数）は580length(na.omit(species))：“NA”を除く要素数（系統割り当てできたASV数）は59length(table(species))：“NA”を除くユニークな要素数（系統割り当てできたユニークな種の数）は54genus_truth：文字列ベクトル
length(genus_truth)：要素数（ユニークな属の数）は19species_truth：文字列ベクトル
length(species_truth)：要素数（ユニークな種の数）は22第23回の表2と同じものですが、現在取り扱っている計12サンプルの概要を以下に示しておきます。これは、SRP125723 (Hallmaier-Wacker et al., Sci Rep., 2018)のサブセットです。ここでやろうとしていることは、どのDNA抽出キットと保存バッファの組合せが、本来検出したい既知微生物群集HM-280をどれだけうまくカバーしているか？です。
| サンプル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 |
getFreqPerSample関数 getFreqPerSample関数は、W3.6.9で定義したget_W3.6関数をベースに作成しています。get_W3.6関数は数値行列dataを入力として受け付けて、内部的に行ごとの総和をrowSums(data)で計算していました。getFreqPerSample関数は、内部的に行ごとの総和を計算しないようにしています。また、感度の計算部分で、サンプルによって検出できていない（出現頻度が0の）ASV配列も存在しますので、出現頻度が0より大きいものに限定しています。この関数は、W3.6.11で述べたJSLAB_source.Rの中でも定義されています。
getFreqPerSample <- function(data, prediction, truth){
# W3.5.4
frequency <- data # dataをfrequencyにコピー
obj_TP <- is.element(prediction, truth) # prediction中の各要素がtruthに含まれるかどうかの判定結果をobj_TPに格納
hoge1 <- prediction[obj_TP] # prediction[obj_TP]をhoge1にコピー
hoge2 <- frequency[obj_TP] # frequency[obj_TP]をhoge2にコピー
freq <- rep(x = hoge1, times = hoge2) # hoge1をhoge2の指定回数分生成した結果をfreqに格納
hoge4 <- append(freq, truth) # freqとtruthを連結した結果をhoge4に格納
hoge5 <- table(hoge4) - 1 # hoge4中の要素の種類ごとの出現頻度-1をhoge5に格納
# W3.4.5
obj_FP <- !obj_TP # obj_TPのTRUEとFALSEを入れ替えた結果をobj_FPに格納
others <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をothersに格納
names(others) <- "Others" # othersに"Others"をnames属性として付与
# Concatenate
hoge6 <- append(hoge5, others) # hoge5とothersを連結した結果をhoge6に格納
# sensitivity
flag <- frequency > 0 # frequency中の要素が0よりも大きいものをTRUE、それ以外をFALSEとした結果をflagに格納
obj <- is.element(truth, prediction[flag]) # truth中の各要素がprediction (但しflagがTRUEに限定)に含まれるかどうかの判定結果をobjに格納
sensi <- sum(obj)/length(truth) # 感度の計算結果をsensiに格納
# precision（FPにNA含む）
TP <- sum(frequency[obj_TP]) # obj_TPがTRUEのfrequencyの総和を計算した結果をTPに格納
FP <- others # othersをFPにコピー
preci1 <- TP / (TP + FP) # 適合率を計算した結果をpreci1に格納
# precision（FPにNA含まない）
obj_P <- !is.na(prediction) # predictionの各要素がNAでないかどうかの判定結果をobj_Pに格納
obj_FP <- obj_P & !obj_TP # obj_Pと!obj_TPの両方でTRUEの要素のみをTRUEとした結果をobj_FPに格納
FP <- sum(frequency[obj_FP]) # obj_FPがTRUEのfrequencyの総和を計算した結果をFPに格納
preci2 <- TP / (TP + FP) # 適合率を計算した結果をpreci2に格納
# output
rtn <- list() # 結果をリスト形式で返すためのrtnを作成
rtn$freque <- hoge6 # hoge6をfrequeという名前で追加
rtn$TPposi <- obj_TP # obj_TPをTPposiという名前で追加
rtn$sensitivity <- sensi # sensiをsensivitityという名前で追加
rtn$precision_with_NA <- preci1 # preci1をprecision_with_NAという名前で追加
rtn$precision_without_NA <- preci2 # preci2をprecision_without_NAという名前で追加
return(rtn) # rtnの中身を関数実行結果として返す
} getFreqPerSample関数の最初の引数として与える情報はベクトルです。ここでは、計12サンプルから構成される数値行列dataうち、1列目の情報のみを入力として与えています。
j <- 1 # 解析したい列番号情報
result <- getFreqPerSample(data[, j], genus, genus_truth)# j列目の出現頻度のまとめ情報を取得してresultに格納
result # resultの中身## $freque
## Acinetobacter Actinomyces Bacillus Bacteroides
## 10653 2688 0 4577
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 126 0 531 7716
## Escherichia Helicobacter Lactobacillus Listeria
## 0 1501 1173 4145
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 81414 13869 0 51798
## Rhodobacter Staphylococcus Streptococcus Others
## 1589 3105 674 73393
##
## $TPposi
## [1] TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE
## [25] TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE
## [37] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [85] FALSE FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [97] TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## [133] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [241] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [253] TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [277] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [289] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [301] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [313] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [325] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [337] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [349] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [373] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [385] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [397] FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## [409] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [421] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [433] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [445] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [457] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [469] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [481] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [493] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [505] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [517] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [541] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [553] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [565] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [577] FALSE FALSE FALSE FALSE
##
## $sensitivity
## [1] 0.7894737
##
## $precision_with_NA
## Others
## 0.7165768
##
## $precision_without_NA
## [1] 0.8284476 （上記の結果をそのまま眺めてもよいのですが）result$frequeで見られる1番目のサンプルのみの頻度情報は、W3.6.10で見られる12サンプルを合わせた頻度情報よりも明らかに少なく妥当だと判断できます。
result$freque # result$frequeの中身## Acinetobacter Actinomyces Bacillus Bacteroides
## 10653 2688 0 4577
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 126 0 531 7716
## Escherichia Helicobacter Lactobacillus Listeria
## 0 1501 1173 4145
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 81414 13869 0 51798
## Rhodobacter Staphylococcus Streptococcus Others
## 1589 3105 674 73393なお、これは以下と実質的に同じです。
j <- 1 # 解析したい列番号情報
getFreqPerSample(data[, j], genus, genus_truth)$freque# j列目の出現頻度のまとめ情報を取得## Acinetobacter Actinomyces Bacillus Bacteroides
## 10653 2688 0 4577
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 126 0 531 7716
## Escherichia Helicobacter Lactobacillus Listeria
## 0 1501 1173 4145
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 81414 13869 0 51798
## Rhodobacter Staphylococcus Streptococcus Others
## 1589 3105 674 73393かなりローテクな感じですが、以下のようにすれば最初の6サンプル分の属ごとの出現頻度情報をまとめることができます。
freque1 <- getFreqPerSample(data[, 1], genus, genus_truth)$freque# 1列目の出現頻度情報を取得してfreque1に格納
freque2 <- getFreqPerSample(data[, 2], genus, genus_truth)$freque# 2列目の出現頻度情報を取得してfreque2に格納
freque3 <- getFreqPerSample(data[, 3], genus, genus_truth)$freque# 3列目の出現頻度情報を取得してfreque3に格納
freque4 <- getFreqPerSample(data[, 4], genus, genus_truth)$freque# 4列目の出現頻度情報を取得してfreque4に格納
freque5 <- getFreqPerSample(data[, 5], genus, genus_truth)$freque# 5列目の出現頻度情報を取得してfreque5に格納
freque6 <- getFreqPerSample(data[, 6], genus, genus_truth)$freque# 6列目の出現頻度情報を取得してfreque6に格納
result <- cbind(freque1, freque2, freque3,# 6サンプル分の頻度情報を列方向で結合した結果をresultに格納
freque4, freque5, freque6)# 6サンプル分の頻度情報を列方向で結合した結果をresultに格納
result # resultの中身## freque1 freque2 freque3 freque4 freque5 freque6
## Acinetobacter 10653 10257 14988 51921 6146 10441
## Actinomyces 2688 1868 3012 4212 835 1072
## Bacillus 0 0 0 0 0 0
## Bacteroides 4577 3624 6660 9912 1057 630
## Bifidobacterium 126 115 288 247 313 527
## Clostridium 0 0 0 0 0 0
## Deinococcus 531 469 882 1974 1133 2633
## Enterococcus 7716 6148 25697 10822 82833 154386
## Escherichia 0 0 0 0 0 0
## Helicobacter 1501 1379 1344 395 245 164
## Lactobacillus 1173 1058 1488 2685 883 968
## Listeria 4145 3485 5254 5484 1156 865
## Neisseria 81414 81823 107177 33468 27625 11635
## Porphyromonas 13869 12798 14114 4283 4177 2143
## Propionibacterium 0 0 0 0 0 0
## Pseudomonas 51798 48895 62117 16931 16186 5764
## Rhodobacter 1589 1705 2302 5158 883 1216
## Staphylococcus 3105 2677 3893 1488 6586 7184
## Streptococcus 674 620 1109 366 2648 5374
## Others 73393 65793 86215 153638 39158 36743apply関数で実行 apply関数は、行列を入力として、行ごとや列ごとに任意の処理を行うための関数です（W3.1.3）。ここでは、apply関数の中でgetFreqPerSample関数を実行する例を示します。580行
\(\times\)
12列からなる数値行列dataを入力として（X = data）、列ごとに（MARGIN = 2）、getFreqPerSample関数を適用する（FUN = getFreqPerSample）、getFreqPerSample関数で与えないといけない残りの引数を記載（genus, genus_truth）という書き方にしています。getFreqPerSample関数自体は、計3つの引数を入力として与えます（W4.2）。1番目の引数はdata、2番目の引数はgenus、そして3番目の引数はgenus_truthを与えて関数を実行したいので、このような書き方になっています。横長になりすぎるので、3行に分けて記載しています。
result <- apply(X = data, MARGIN = 2, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
FUN = getFreqPerSample, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
genus, genus_truth) # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納 str関数は、入力として与えられたオブジェクトの構造（structure）を返します（W3.6.8）。getFreqPerSample関数は、サンプルごとにfreque,
TPposi, sensitivity,
precision_with_NA, and
precision_without_NAの計5種類の結果をリスト形式で返します（W4.1）。以下で見えているresultオブジェクトは、計12サンプル分の結果がリスト形式でまとめられた形になっていることが（全体的な見た目から）わかります。
str(result) # resultの構造を表示## List of 12
## $ SRR6325813:List of 5
## ..$ freque : Named num [1:20] 10653 2688 0 4577 126 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.717
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.828
## $ SRR6325815:List of 5
## ..$ freque : Named num [1:20] 10257 1868 0 3624 115 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.729
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.834
## $ SRR6325816:List of 5
## ..$ freque : Named num [1:20] 14988 3012 0 6660 288 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.744
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.833
## $ SRR6325817:List of 5
## ..$ freque : Named num [1:20] 51921 4212 0 9912 247 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.493
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.558
## $ SRR6325828:List of 5
## ..$ freque : Named num [1:20] 6146 835 0 1057 313 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.796
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.874
## $ SRR6325829:List of 5
## ..$ freque : Named num [1:20] 10441 1072 0 630 527 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.848
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.926
## $ SRR6325830:List of 5
## ..$ freque : Named num [1:20] 6512 1110 0 751 269 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.825
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.902
## $ SRR6325831:List of 5
## ..$ freque : Named num [1:20] 7193 1156 0 686 106 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.792
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.871
## $ SRR6325859:List of 5
## ..$ freque : Named num [1:20] 11293 2527 0 6465 165 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.764
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.846
## $ SRR6325860:List of 5
## ..$ freque : Named num [1:20] 15193 1740 0 5152 87 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.733
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.824
## $ SRR6325861:List of 5
## ..$ freque : Named num [1:20] 16686 886 0 1230 46 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.557
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.608
## $ SRR6325862:List of 5
## ..$ freque : Named num [1:20] 16852 1123 0 3948 110 ...
## .. ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ TPposi : logi [1:580] TRUE TRUE FALSE TRUE FALSE TRUE ...
## ..$ sensitivity : num 0.789
## ..$ precision_with_NA : Named num 0.678
## .. ..- attr(*, "names")= chr "Others"
## ..$ precision_without_NA: num 0.745 ここでは、W4.3で得たリスト形式のresultオブジェクトを入力として、サンプルごとの出現頻度をまとめて数値行列として得たり、感度や適合率をまとめて数値行列の下部に追加するなどの作業を行います。
ここでは、resultオブジェクトの中から、サンプルごとの結果の1番目の要素であるfreque情報を抽出して、計12サンプル分の属ごとの出現頻度を数値行列として得る作業を行います。入力であるresultがリスト形式ですので、ここではsapply関数を用いて、リスト形式のオブジェクトに対して逐次的な処理を行います。具体的には、W2.2.4と同じような処理を行うべく、オブジェクト中の指定した位置の情報を抽出する[という関数を実行し、1という情報をオプションとして与えています。sapply関数実行結果であるhogeオブジェクトの構造をstr関数で表示させた結果（str(hoge)）より、数値行列形式ではなくリスト形式ではあるものの、「サンプルごとの結果の1番目の要素であるfreque情報」のみの抽出できているっぽいことがわかります。
hoge <- sapply(X = result, FUN = `[`, 1) # result中の各要素から、1番目の要素の抽出結果をhogeに格納
str(hoge) # hogeの構造を表示## List of 12
## $ SRR6325813.freque: Named num [1:20] 10653 2688 0 4577 126 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325815.freque: Named num [1:20] 10257 1868 0 3624 115 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325816.freque: Named num [1:20] 14988 3012 0 6660 288 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325817.freque: Named num [1:20] 51921 4212 0 9912 247 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325828.freque: Named num [1:20] 6146 835 0 1057 313 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325829.freque: Named num [1:20] 10441 1072 0 630 527 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325830.freque: Named num [1:20] 6512 1110 0 751 269 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325831.freque: Named num [1:20] 7193 1156 0 686 106 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325859.freque: Named num [1:20] 11293 2527 0 6465 165 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325860.freque: Named num [1:20] 15193 1740 0 5152 87 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325861.freque: Named num [1:20] 16686 886 0 1230 46 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325862.freque: Named num [1:20] 16852 1123 0 3948 110 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ... リスト形式になっているhoge中の各要素はfreque情報になります。つまり、hoge中の「各要素」は、「サンプルごとの20個の要素からなる数値ベクトル」です。したがって、「リスト形式になっているhoge中の各要素」をcbind関数を用いて列方向で結合してやれば、得られる結果（つまりfreq_mat）は数値行列になります。
freq_mat <- sapply(X = hoge, FUN = cbind) # hoge中の各要素を列方向で結合した結果をfreq_matに格納
str(freq_mat) # freq_matの構造を表示## num [1:20, 1:12] 10653 2688 0 4577 126 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:12] "SRR6325813.freque" "SRR6325815.freque" "SRR6325816.freque" "SRR6325817.freque" ... str(freq_mat)実行結果の最初のほうで見えている「num
[1:20, 1:12]」より、20行 \(\times\)
12列の数値行列であることがなんとなくわかると思います。str関数はリスト形式など比較的複雑なオブジェクトの全体像を把握する際には便利ですが、数値行列のような簡単なオブジェクトのときにはかえってわかりづらいかもしれません。freq_matの場合は、dim(freq_mat)やhead(freq_mat)とかで大まかに把握するほうがわかりやすいかもしれません。最初の6列分については、W4.2の最後のほうで見えている数値行列と一致していることがわかります。W4.2は無駄と思われたかたもいらっしゃるかもしれませんが、このような答え合わせとして使えます。
dim(freq_mat) # freq_matの行数と列数## [1] 20 12head(freq_mat) # freq_mat中の中の最初の6つを表示## SRR6325813.freque SRR6325815.freque SRR6325816.freque SRR6325817.freque
## [1,] 10653 10257 14988 51921
## [2,] 2688 1868 3012 4212
## [3,] 0 0 0 0
## [4,] 4577 3624 6660 9912
## [5,] 126 115 288 247
## [6,] 0 0 0 0
## SRR6325828.freque SRR6325829.freque SRR6325830.freque SRR6325831.freque
## [1,] 6146 10441 6512 7193
## [2,] 835 1072 1110 1156
## [3,] 0 0 0 0
## [4,] 1057 630 751 686
## [5,] 313 527 269 106
## [6,] 0 0 0 0
## SRR6325859.freque SRR6325860.freque SRR6325861.freque SRR6325862.freque
## [1,] 11293 15193 16686 16852
## [2,] 2527 1740 886 1123
## [3,] 0 0 0 0
## [4,] 6465 5152 1230 3948
## [5,] 165 87 46 110
## [6,] 0 0 0 0rownames関数 前述の表示結果からもわかりますが、freq_matには列名情報は含まれていますが、行名情報は含まれていません。このような場合に、rownames関数を用いて行列に行名を付与することができます。行名として付与するのは基本的に属名情報になります。さきほどW4.4.1の最初のほうで得たhogeオブジェクトをよくみると、属名情報が格納されていることがわかります。
hoge <- sapply(X = result, FUN = `[`, 1) # result中の各要素から、1番目の要素の抽出結果をhogeに格納
str(hoge) # hogeの構造を表示## List of 12
## $ SRR6325813.freque: Named num [1:20] 10653 2688 0 4577 126 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325815.freque: Named num [1:20] 10257 1868 0 3624 115 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325816.freque: Named num [1:20] 14988 3012 0 6660 288 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325817.freque: Named num [1:20] 51921 4212 0 9912 247 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325828.freque: Named num [1:20] 6146 835 0 1057 313 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325829.freque: Named num [1:20] 10441 1072 0 630 527 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325830.freque: Named num [1:20] 6512 1110 0 751 269 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325831.freque: Named num [1:20] 7193 1156 0 686 106 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325859.freque: Named num [1:20] 11293 2527 0 6465 165 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325860.freque: Named num [1:20] 15193 1740 0 5152 87 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325861.freque: Named num [1:20] 16686 886 0 1230 46 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## $ SRR6325862.freque: Named num [1:20] 16852 1123 0 3948 110 ...
## ..- attr(*, "names")= chr [1:20] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ... 具体的には、「hogeオブジェクトに含まれる、計12サンプル中のどの要素（サンプル）にも属名情報が含まれている」ことがわかります。例えばリスト形式の1番目の要素を抽出したい場合は、以下のようにします。
hoge[[1]] # hoge中の1番目の要素## Acinetobacter Actinomyces Bacillus Bacteroides
## 10653 2688 0 4577
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 126 0 531 7716
## Escherichia Helicobacter Lactobacillus Listeria
## 0 1501 1173 4145
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 81414 13869 0 51798
## Rhodobacter Staphylococcus Streptococcus Others
## 1589 3105 674 73393 これのnames属性情報が、freq_matに列名として与えたい情報になります。
names(hoge[[1]]) # hoge中の1番目の要素のnames属性## [1] "Acinetobacter" "Actinomyces" "Bacillus"
## [4] "Bacteroides" "Bifidobacterium" "Clostridium"
## [7] "Deinococcus" "Enterococcus" "Escherichia"
## [10] "Helicobacter" "Lactobacillus" "Listeria"
## [13] "Neisseria" "Porphyromonas" "Propionibacterium"
## [16] "Pseudomonas" "Rhodobacter" "Staphylococcus"
## [19] "Streptococcus" "Others" 残りは、freq_matの行名情報としてnames(hoge[[1]])を代入するだけです。
rownames(freq_mat) <- names(hoge[[1]]) # names(hoge[[1]])をfreq_matの行名情報とする念のため確認しているだけです。
freq_mat # freq_matの中身## SRR6325813.freque SRR6325815.freque SRR6325816.freque
## Acinetobacter 10653 10257 14988
## Actinomyces 2688 1868 3012
## Bacillus 0 0 0
## Bacteroides 4577 3624 6660
## Bifidobacterium 126 115 288
## Clostridium 0 0 0
## Deinococcus 531 469 882
## Enterococcus 7716 6148 25697
## Escherichia 0 0 0
## Helicobacter 1501 1379 1344
## Lactobacillus 1173 1058 1488
## Listeria 4145 3485 5254
## Neisseria 81414 81823 107177
## Porphyromonas 13869 12798 14114
## Propionibacterium 0 0 0
## Pseudomonas 51798 48895 62117
## Rhodobacter 1589 1705 2302
## Staphylococcus 3105 2677 3893
## Streptococcus 674 620 1109
## Others 73393 65793 86215
## SRR6325817.freque SRR6325828.freque SRR6325829.freque
## Acinetobacter 51921 6146 10441
## Actinomyces 4212 835 1072
## Bacillus 0 0 0
## Bacteroides 9912 1057 630
## Bifidobacterium 247 313 527
## Clostridium 0 0 0
## Deinococcus 1974 1133 2633
## Enterococcus 10822 82833 154386
## Escherichia 0 0 0
## Helicobacter 395 245 164
## Lactobacillus 2685 883 968
## Listeria 5484 1156 865
## Neisseria 33468 27625 11635
## Porphyromonas 4283 4177 2143
## Propionibacterium 0 0 0
## Pseudomonas 16931 16186 5764
## Rhodobacter 5158 883 1216
## Staphylococcus 1488 6586 7184
## Streptococcus 366 2648 5374
## Others 153638 39158 36743
## SRR6325830.freque SRR6325831.freque SRR6325859.freque
## Acinetobacter 6512 7193 11293
## Actinomyces 1110 1156 2527
## Bacillus 0 0 0
## Bacteroides 751 686 6465
## Bifidobacterium 269 106 165
## Clostridium 0 0 0
## Deinococcus 1110 254 595
## Enterococcus 90159 67629 5505
## Escherichia 0 0 0
## Helicobacter 408 296 1056
## Lactobacillus 1094 905 1460
## Listeria 1454 2027 3202
## Neisseria 45205 52283 88802
## Porphyromonas 6555 6804 18322
## Propionibacterium 0 0 0
## Pseudomonas 21605 22188 49964
## Rhodobacter 673 386 1878
## Staphylococcus 4487 2666 1958
## Streptococcus 2926 1625 528
## Others 38994 43671 59939
## SRR6325860.freque SRR6325861.freque SRR6325862.freque
## Acinetobacter 15193 16686 16852
## Actinomyces 1740 886 1123
## Bacillus 0 0 0
## Bacteroides 5152 1230 3948
## Bifidobacterium 87 46 110
## Clostridium 0 0 0
## Deinococcus 480 558 536
## Enterococcus 5765 19608 17070
## Escherichia 0 0 0
## Helicobacter 973 61 612
## Lactobacillus 1782 744 1115
## Listeria 2619 2176 3951
## Neisseria 74848 11007 70355
## Porphyromonas 14954 1411 10555
## Propionibacterium 0 0 0
## Pseudomonas 39483 7028 34481
## Rhodobacter 1405 935 1102
## Staphylococcus 1431 277 1494
## Streptococcus 464 847 669
## Others 60735 50553 77993 本質的な部分ではありませんでしたが、気になるとずっと気になり続ける部分ですし、。また、最後の最後で自動的に結果をまとめて表として作成したい場合などに、後で苦労することになる場合もあります。実際、W5の作成途中でfreq_mat中に行名情報が存在しないことに気づき、最後の最後ででこの項目を追加する羽目に遭っています。
ここでは、resultオブジェクトの中から、サンプルごとの結果の3番目の要素であるsensitivity情報を抽出して、計12サンプル分の属ごとの感度情報のみを得ます。
hoge <- sapply(X = result, FUN = `[`, 3) # result中の各要素から3番目の要素の抽出結果をhogeに格納
str(hoge) # hogeの構造を表示## List of 12
## $ SRR6325813.sensitivity: num 0.789
## $ SRR6325815.sensitivity: num 0.789
## $ SRR6325816.sensitivity: num 0.789
## $ SRR6325817.sensitivity: num 0.789
## $ SRR6325828.sensitivity: num 0.789
## $ SRR6325829.sensitivity: num 0.789
## $ SRR6325830.sensitivity: num 0.789
## $ SRR6325831.sensitivity: num 0.789
## $ SRR6325859.sensitivity: num 0.789
## $ SRR6325860.sensitivity: num 0.789
## $ SRR6325861.sensitivity: num 0.789
## $ SRR6325862.sensitivity: num 0.789 出力結果であるhogeオブジェクトは、12個の要素からなるリスト形式であることがわかります。1つ1つの要素はスカラー（1つの数字）で、すべて0.789であることがわかります。また、W4.4.1のhead(freq_mat)実行結果で見えている3行目と6行目に0が並んでいる状況を踏まえ、「おそらく全てのサンプルが同じ属を検出できていないのだろう。そのうちの2つは、3行目と6行目の属。15/19
=
0.7894737なので、検出できていないのは4属。」のように解釈します。表示結果の分量もたかが知れているので、以下にhoge実行結果も示します。
hoge # hogeの中身## $SRR6325813.sensitivity
## [1] 0.7894737
##
## $SRR6325815.sensitivity
## [1] 0.7894737
##
## $SRR6325816.sensitivity
## [1] 0.7894737
##
## $SRR6325817.sensitivity
## [1] 0.7894737
##
## $SRR6325828.sensitivity
## [1] 0.7894737
##
## $SRR6325829.sensitivity
## [1] 0.7894737
##
## $SRR6325830.sensitivity
## [1] 0.7894737
##
## $SRR6325831.sensitivity
## [1] 0.7894737
##
## $SRR6325859.sensitivity
## [1] 0.7894737
##
## $SRR6325860.sensitivity
## [1] 0.7894737
##
## $SRR6325861.sensitivity
## [1] 0.7894737
##
## $SRR6325862.sensitivity
## [1] 0.7894737 第24回のW1.6.5でも利用しましたが、上記hogeのような感じで見えているリスト形式のものは、リスト形式を解除するunlist関数を実行してやれば、ベクトル形式に変更することができます。
unlist(hoge) # hogeをベクトル形式に変更## SRR6325813.sensitivity SRR6325815.sensitivity SRR6325816.sensitivity
## 0.7894737 0.7894737 0.7894737
## SRR6325817.sensitivity SRR6325828.sensitivity SRR6325829.sensitivity
## 0.7894737 0.7894737 0.7894737
## SRR6325830.sensitivity SRR6325831.sensitivity SRR6325859.sensitivity
## 0.7894737 0.7894737 0.7894737
## SRR6325860.sensitivity SRR6325861.sensitivity SRR6325862.sensitivity
## 0.7894737 0.7894737 0.7894737ここまでの一連の作業は、以下のような感じでまとめることができます。
sensitivity <- unlist(sapply(X = result, FUN = `[`, 3)) # result中の各要素から3番目の要素のみ抽出した結果をベクトル形式にしてsensitivityに格納
sensitivity # sensitivityの中身## SRR6325813.sensitivity SRR6325815.sensitivity SRR6325816.sensitivity
## 0.7894737 0.7894737 0.7894737
## SRR6325817.sensitivity SRR6325828.sensitivity SRR6325829.sensitivity
## 0.7894737 0.7894737 0.7894737
## SRR6325830.sensitivity SRR6325831.sensitivity SRR6325859.sensitivity
## 0.7894737 0.7894737 0.7894737
## SRR6325860.sensitivity SRR6325861.sensitivity SRR6325862.sensitivity
## 0.7894737 0.7894737 0.7894737 ここでは、resultオブジェクトの中から、サンプルごとの結果の4番目の要素であるprecision_with_NA情報、および5番目の要素であるprecision_without_NA情報を抽出します（W3.6.9）。TP
/ (TP +
FP)で定義される適合率は、FPを「TP以外の全てのASV数（NAも含める）」のようにNAを含める場合が前者（precision_with_NA）、含めない場合が後者（precision_without_NA）の情報になります。
precision_with_NA <- unlist(sapply(X = result, FUN = `[`, 4)) # result中の各要素から4番目の要素のみ抽出した結果をベクトル形式にしてprecision_with_NAに格納
precision_without_NA <- unlist(sapply(X = result, FUN = `[`, 5)) # result中の各要素から5番目の要素のみ抽出した結果をベクトル形式にしてprecision_without_NAに格納 precision_with_NAの中身を念のため示します。
precision_with_NA # precision_with_NAの中身## SRR6325813.precision_with_NA.Others SRR6325815.precision_with_NA.Others
## 0.7165768 0.7289279
## SRR6325816.precision_with_NA.Others SRR6325817.precision_with_NA.Others
## 0.7438195 0.4929171
## SRR6325828.precision_with_NA.Others SRR6325829.precision_with_NA.Others
## 0.7959075 0.8480093
## SRR6325830.precision_with_NA.Others SRR6325831.precision_with_NA.Others
## 0.8253833 0.7919190
## SRR6325859.precision_with_NA.Others SRR6325860.precision_with_NA.Others
## 0.7637025 0.7325757
## SRR6325861.precision_with_NA.Others SRR6325862.precision_with_NA.Others
## 0.5567587 0.6776696 precision_without_NAの中身を念のため示します。
precision_without_NA # precision_without_NAの中身## SRR6325813.precision_without_NA SRR6325815.precision_without_NA
## 0.8284476 0.8337307
## SRR6325816.precision_without_NA SRR6325817.precision_without_NA
## 0.8325379 0.5581588
## SRR6325828.precision_without_NA SRR6325829.precision_without_NA
## 0.8735242 0.9259893
## SRR6325830.precision_without_NA SRR6325831.precision_without_NA
## 0.9019059 0.8711495
## SRR6325859.precision_without_NA SRR6325860.precision_without_NA
## 0.8460054 0.8239576
## SRR6325861.precision_without_NA SRR6325862.precision_without_NA
## 0.6077427 0.7447868 W4.3とW4.4のスクリプトを以下にまとめます。これを実行するために必要な3つのオブジェクト（data,
genus_truth, and
genus）は、W2.3を実行すれば得られます。
source("http://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB_source.R") # 関数の読込
result <- apply(X = data, MARGIN = 2, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
FUN = getFreqPerSample, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
genus, genus_truth) # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
hoge <- sapply(X = result, FUN = `[`, 1) # result中の各要素から、1番目の要素の抽出結果をhogeに格納
freq_mat <- sapply(X = hoge, FUN = cbind) # hoge中の各要素を列方向で結合した結果をfreq_matに格納
rownames(freq_mat) <- names(hoge[[1]]) # names(hoge[[1]])をfreq_matの行名情報とする
sensitivity <- unlist(sapply(X = result, FUN = `[`, 3)) # result中の各要素から3番目の要素のみ抽出した結果をベクトル形式にしてsensitivityに格納
precision_with_NA <- unlist(sapply(X = result, FUN = `[`, 4)) # result中の各要素から4番目の要素のみ抽出した結果をベクトル形式にしてprecision_with_NAに格納
precision_without_NA <- unlist(sapply(X = result, FUN = `[`, 5)) # result中の各要素から5番目の要素のみ抽出した結果をベクトル形式にしてprecision_without_NAに格納以下の4つのオブジェクトが主な結果です。
freq_mat
sensitivity
precision_with_NA
precision_without_NA W4.4は属レベルのものでしたが、種レベルについてもW4.4.5の冒頭部分を変更するだけで簡単に結果が得られます。これを実行するために必要な3つのオブジェクト（data,
species_truth, and
species）は、W2.3を実行すれば得られます。
source("http://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB_source.R") # 関数の読込
result <- apply(X = data, MARGIN = 2, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
FUN = getFreqPerSample, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
species, species_truth) # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
hoge <- sapply(X = result, FUN = `[`, 1) # result中の各要素から、1番目の要素の抽出結果をhogeに格納
freq_mat <- sapply(X = hoge, FUN = cbind) # hoge中の各要素を列方向で結合した結果をfreq_matに格納
rownames(freq_mat) <- names(hoge[[1]]) # names(hoge[[1]])をfreq_matの行名情報とする
sensitivity <- unlist(sapply(X = result, FUN = `[`, 3)) # result中の各要素から3番目の要素のみ抽出した結果をベクトル形式にしてsensitivityに格納
precision_with_NA <- unlist(sapply(X = result, FUN = `[`, 4)) # result中の各要素から4番目の要素のみ抽出した結果をベクトル形式にしてprecision_with_NAに格納
precision_without_NA <- unlist(sapply(X = result, FUN = `[`, 5)) # result中の各要素から5番目の要素のみ抽出した結果をベクトル形式にしてprecision_witho種レベルの出現頻度の数値行列は以下のような感じになります。
freq_mat## SRR6325813.freque SRR6325815.freque SRR6325816.freque
## acnes 0 0 0
## adolescentis 0 0 0
## aeruginosa 0 0 0
## agalactiae 0 0 0
## aureus 0 0 0
## baumannii 0 0 0
## beijerinckii 0 37 0
## cereus 0 0 0
## coli 0 0 0
## epidermidis 0 0 0
## faecalis 0 0 0
## gasseri 0 0 0
## gingivalis 13645 12798 13838
## meningitides 0 0 0
## monocytogenes 0 0 0
## mutans 89 60 164
## odontolyticus 2688 1868 3012
## pneumoniae 0 0 0
## pylori 1501 1379 1344
## radiodurans 531 469 882
## sphaeroides 0 0 0
## vulgatus 4577 3624 6657
## Others 235921 222479 310643
## SRR6325817.freque SRR6325828.freque SRR6325829.freque
## acnes 0 0 0
## adolescentis 0 0 0
## aeruginosa 0 0 0
## agalactiae 0 0 0
## aureus 0 33 14
## baumannii 0 0 0
## beijerinckii 0 0 0
## cereus 0 0 0
## coli 0 0 0
## epidermidis 0 0 0
## faecalis 0 0 0
## gasseri 0 0 0
## gingivalis 4283 4177 2143
## meningitides 0 0 0
## monocytogenes 0 0 0
## mutans 90 249 492
## odontolyticus 4212 835 1072
## pneumoniae 0 0 0
## pylori 395 245 164
## radiodurans 1974 1127 2633
## sphaeroides 0 0 0
## vulgatus 9901 1046 626
## Others 282129 184152 234601
## SRR6325830.freque SRR6325831.freque SRR6325859.freque
## acnes 0 0 0
## adolescentis 0 0 0
## aeruginosa 0 0 0
## agalactiae 0 0 0
## aureus 0 0 0
## baumannii 0 0 0
## beijerinckii 0 0 39
## cereus 0 0 0
## coli 0 0 0
## epidermidis 0 0 0
## faecalis 0 0 0
## gasseri 0 0 0
## gingivalis 6555 6804 18026
## meningitides 0 0 0
## monocytogenes 0 0 0
## mutans 441 245 54
## odontolyticus 1102 1156 2527
## pneumoniae 0 0 0
## pylori 408 296 1056
## radiodurans 1110 254 595
## sphaeroides 0 0 0
## vulgatus 751 686 6465
## Others 212945 200434 224897
## SRR6325860.freque SRR6325861.freque SRR6325862.freque
## acnes 0 0 0
## adolescentis 0 0 0
## aeruginosa 0 0 0
## agalactiae 0 0 0
## aureus 0 0 0
## baumannii 0 0 0
## beijerinckii 0 0 0
## cereus 0 0 0
## coli 0 0 0
## epidermidis 0 0 0
## faecalis 0 0 0
## gasseri 0 0 0
## gingivalis 14954 1411 10555
## meningitides 0 0 0
## monocytogenes 0 0 0
## mutans 75 72 50
## odontolyticus 1719 886 1123
## pneumoniae 0 0 0
## pylori 969 61 612
## radiodurans 480 558 536
## sphaeroides 0 0 0
## vulgatus 5152 1230 3948
## Others 203762 109835 225142種レベルの感度は以下のような感じになります。
sensitivity## SRR6325813.sensitivity SRR6325815.sensitivity SRR6325816.sensitivity
## 0.2727273 0.3181818 0.2727273
## SRR6325817.sensitivity SRR6325828.sensitivity SRR6325829.sensitivity
## 0.2727273 0.3181818 0.3181818
## SRR6325830.sensitivity SRR6325831.sensitivity SRR6325859.sensitivity
## 0.2727273 0.2727273 0.3181818
## SRR6325860.sensitivity SRR6325861.sensitivity SRR6325862.sensitivity
## 0.2727273 0.2727273 0.2727273種レベルの適合率（NAを含めた場合）は以下のような感じになります。
precision_with_NA## SRR6325813.precision_with_NA.Others SRR6325815.precision_with_NA.Others
## 0.08893926 0.08336973
## SRR6325816.precision_with_NA.Others SRR6325817.precision_with_NA.Others
## 0.07695073 0.06883202
## SRR6325828.precision_with_NA.Others SRR6325829.precision_with_NA.Others
## 0.04019514 0.02955180
## SRR6325830.precision_with_NA.Others SRR6325831.precision_with_NA.Others
## 0.04642384 0.04498392
## SRR6325859.precision_with_NA.Others SRR6325860.precision_with_NA.Others
## 0.11338845 0.10280876
## SRR6325861.precision_with_NA.Others SRR6325862.precision_with_NA.Others
## 0.03698281 0.06953043種レベルの適合率（NAを含めない場合）は以下のような感じになります。
precision_without_NA## SRR6325813.precision_without_NA SRR6325815.precision_without_NA
## 0.2204894 0.1981454
## SRR6325816.precision_without_NA SRR6325817.precision_without_NA
## 0.1945461 0.3835332
## SRR6325828.precision_without_NA SRR6325829.precision_without_NA
## 0.2173864 0.3753678
## SRR6325830.precision_without_NA SRR6325831.precision_without_NA
## 0.1863261 0.1528758
## SRR6325859.precision_without_NA SRR6325860.precision_without_NA
## 0.2446185 0.2377287
## SRR6325861.precision_without_NA SRR6325862.precision_without_NA
## 0.2758305 0.1928938 ここでは、取得した出現頻度・感度・適合率の情報をまとめる作業を示します。実行に必要な5つのオブジェクト（data,
genus_truth, genus,
species_truth, and
species）はW2.3を実行すれば得られますが、以下にも再掲します。
# W2.1 dada2実行結果ファイルのほう
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t")# in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
obj <- grep(pattern = "SRR", x = colnames(hoge)) # colnames(hoge)の中から"SRR"という文字を含む要素番号を抽出した結果をobjに格納
data <- hoge[, obj] # hoge行列中のobjで指定した列のみ抽出した結果をdataに格納
genus <- hoge$Genus # hoge中のGenus列を抽出した結果をgenusに格納
species <- hoge$Species # hoge中のSpecies列を抽出した結果をspeciesに格納
# W2.2 菌株情報ファイルのほう
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB23_table1.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t") #in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
organism <- hoge$Organism # hoge中のOrganism列を抽出した結果をorganismに格納
hoge <- strsplit(x = organism, split = ",")# organism中の文字列ベクトルをコンマで分割した結果をhogeに格納
hoge2 <- sapply(X = hoge, FUN = `[`, 1) # hoge中のベクトル形式の各要素から、1番目の要素の抽出結果をhoge2に格納
hoge3 <- strsplit(x = hoge2, split = " ") # hoge2中の文字列ベクトルを半角スペースで分割した結果をhoge3に格納
genus_truth <- sapply(X = hoge3, FUN = `[`, 1) # hoge3中のベクトル形式の各要素から、1番目の要素の抽出結果をgenus_truthに格納
species_truth <- sapply(X = hoge3, FUN = `[`, 2)# hoge3中のベクトル形式の各要素から、2番目の要素の抽出結果をspecies_truthに格納
genus_truth <- unique(genus_truth) # genus_truth中のユニークな要素をgenus_truthに格納
species_truth <- unique(species_truth) # species_truth中のユニークな要素をspecies_truthに格納ここでは、属レベルの抽出結果のまとめを行います。
W4.4.5を再実行し、サンプルごとの属レベルの出現頻度、感度、適合率のオブジェクトを得ておきます。
source("http://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB_source.R") # 関数の読込
result <- apply(X = data, MARGIN = 2, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
FUN = getFreqPerSample, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
genus, genus_truth) # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
hoge <- sapply(X = result, FUN = `[`, 1) # result中の各要素から、1番目の要素の抽出結果をhogeに格納
freq_mat <- sapply(X = hoge, FUN = cbind) # hoge中の各要素を列方向で結合した結果をfreq_matに格納
rownames(freq_mat) <- names(hoge[[1]]) # names(hoge[[1]])をfreq_matの行名情報とする
sensitivity <- unlist(sapply(X = result, FUN = `[`, 3)) # result中の各要素から3番目の要素のみ抽出した結果をベクトル形式にしてsensitivityに格納
precision_with_NA <- unlist(sapply(X = result, FUN = `[`, 4)) # result中の各要素から4番目の要素のみ抽出した結果をベクトル形式にしてprecision_with_NAに格納
precision_without_NA <- unlist(sapply(X = result, FUN = `[`, 5)) # result中の各要素から5番目の要素のみ抽出した結果をベクトル形式にしてprecision_without_NAに格納 ここでは、出現頻度、感度、適合率の結果を、行方向で結合させます。rbind関数
は、同じ要素数の複数のベクトルなどを行（row）方向で結合（bind）させたい場合に利用します。
result_pers <- rbind(freq_mat, # 行方向で結合した結果をresult_persに格納
sensitivity, # 行方向で結合した結果をresult_persに格納
precision_with_NA, # 行方向で結合した結果をresult_persに格納
precision_without_NA)# 行方向で結合した結果をresult_persに格納得られるオブジェクトは以下のような感じになります。
str(result_pers) # result_persの構造を表示## num [1:23, 1:12] 10653 2688 0 4577 126 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:23] "Acinetobacter" "Actinomyces" "Bacillus" "Bacteroides" ...
## ..$ : chr [1:12] "SRR6325813.freque" "SRR6325815.freque" "SRR6325816.freque" "SRR6325817.freque" ...dim(result_pers) # result_persの行数と列数## [1] 23 12 result_persは、23行 \(\times\) 12列からなることがわかります。20行
\(\times\)
12列からなる行列freq_matの下部に、列数と同じ要素からなる3つのベクトル情報をそれぞれ結合して得られたものなので妥当ですね。
まず、result_persの中身を示します。指数表記で若干見づらいかもしれませんが、数値自体は妥当です。
result_pers # result_persの中身を表示## SRR6325813.freque SRR6325815.freque SRR6325816.freque
## Acinetobacter 1.065300e+04 1.025700e+04 1.498800e+04
## Actinomyces 2.688000e+03 1.868000e+03 3.012000e+03
## Bacillus 0.000000e+00 0.000000e+00 0.000000e+00
## Bacteroides 4.577000e+03 3.624000e+03 6.660000e+03
## Bifidobacterium 1.260000e+02 1.150000e+02 2.880000e+02
## Clostridium 0.000000e+00 0.000000e+00 0.000000e+00
## Deinococcus 5.310000e+02 4.690000e+02 8.820000e+02
## Enterococcus 7.716000e+03 6.148000e+03 2.569700e+04
## Escherichia 0.000000e+00 0.000000e+00 0.000000e+00
## Helicobacter 1.501000e+03 1.379000e+03 1.344000e+03
## Lactobacillus 1.173000e+03 1.058000e+03 1.488000e+03
## Listeria 4.145000e+03 3.485000e+03 5.254000e+03
## Neisseria 8.141400e+04 8.182300e+04 1.071770e+05
## Porphyromonas 1.386900e+04 1.279800e+04 1.411400e+04
## Propionibacterium 0.000000e+00 0.000000e+00 0.000000e+00
## Pseudomonas 5.179800e+04 4.889500e+04 6.211700e+04
## Rhodobacter 1.589000e+03 1.705000e+03 2.302000e+03
## Staphylococcus 3.105000e+03 2.677000e+03 3.893000e+03
## Streptococcus 6.740000e+02 6.200000e+02 1.109000e+03
## Others 7.339300e+04 6.579300e+04 8.621500e+04
## sensitivity 7.894737e-01 7.894737e-01 7.894737e-01
## precision_with_NA 7.165768e-01 7.289279e-01 7.438195e-01
## precision_without_NA 8.284476e-01 8.337307e-01 8.325379e-01
## SRR6325817.freque SRR6325828.freque SRR6325829.freque
## Acinetobacter 5.192100e+04 6.146000e+03 1.044100e+04
## Actinomyces 4.212000e+03 8.350000e+02 1.072000e+03
## Bacillus 0.000000e+00 0.000000e+00 0.000000e+00
## Bacteroides 9.912000e+03 1.057000e+03 6.300000e+02
## Bifidobacterium 2.470000e+02 3.130000e+02 5.270000e+02
## Clostridium 0.000000e+00 0.000000e+00 0.000000e+00
## Deinococcus 1.974000e+03 1.133000e+03 2.633000e+03
## Enterococcus 1.082200e+04 8.283300e+04 1.543860e+05
## Escherichia 0.000000e+00 0.000000e+00 0.000000e+00
## Helicobacter 3.950000e+02 2.450000e+02 1.640000e+02
## Lactobacillus 2.685000e+03 8.830000e+02 9.680000e+02
## Listeria 5.484000e+03 1.156000e+03 8.650000e+02
## Neisseria 3.346800e+04 2.762500e+04 1.163500e+04
## Porphyromonas 4.283000e+03 4.177000e+03 2.143000e+03
## Propionibacterium 0.000000e+00 0.000000e+00 0.000000e+00
## Pseudomonas 1.693100e+04 1.618600e+04 5.764000e+03
## Rhodobacter 5.158000e+03 8.830000e+02 1.216000e+03
## Staphylococcus 1.488000e+03 6.586000e+03 7.184000e+03
## Streptococcus 3.660000e+02 2.648000e+03 5.374000e+03
## Others 1.536380e+05 3.915800e+04 3.674300e+04
## sensitivity 7.894737e-01 7.894737e-01 7.894737e-01
## precision_with_NA 4.929171e-01 7.959075e-01 8.480093e-01
## precision_without_NA 5.581588e-01 8.735242e-01 9.259893e-01
## SRR6325830.freque SRR6325831.freque SRR6325859.freque
## Acinetobacter 6.512000e+03 7.193000e+03 1.129300e+04
## Actinomyces 1.110000e+03 1.156000e+03 2.527000e+03
## Bacillus 0.000000e+00 0.000000e+00 0.000000e+00
## Bacteroides 7.510000e+02 6.860000e+02 6.465000e+03
## Bifidobacterium 2.690000e+02 1.060000e+02 1.650000e+02
## Clostridium 0.000000e+00 0.000000e+00 0.000000e+00
## Deinococcus 1.110000e+03 2.540000e+02 5.950000e+02
## Enterococcus 9.015900e+04 6.762900e+04 5.505000e+03
## Escherichia 0.000000e+00 0.000000e+00 0.000000e+00
## Helicobacter 4.080000e+02 2.960000e+02 1.056000e+03
## Lactobacillus 1.094000e+03 9.050000e+02 1.460000e+03
## Listeria 1.454000e+03 2.027000e+03 3.202000e+03
## Neisseria 4.520500e+04 5.228300e+04 8.880200e+04
## Porphyromonas 6.555000e+03 6.804000e+03 1.832200e+04
## Propionibacterium 0.000000e+00 0.000000e+00 0.000000e+00
## Pseudomonas 2.160500e+04 2.218800e+04 4.996400e+04
## Rhodobacter 6.730000e+02 3.860000e+02 1.878000e+03
## Staphylococcus 4.487000e+03 2.666000e+03 1.958000e+03
## Streptococcus 2.926000e+03 1.625000e+03 5.280000e+02
## Others 3.899400e+04 4.367100e+04 5.993900e+04
## sensitivity 7.894737e-01 7.894737e-01 7.894737e-01
## precision_with_NA 8.253833e-01 7.919190e-01 7.637025e-01
## precision_without_NA 9.019059e-01 8.711495e-01 8.460054e-01
## SRR6325860.freque SRR6325861.freque SRR6325862.freque
## Acinetobacter 1.519300e+04 1.668600e+04 1.685200e+04
## Actinomyces 1.740000e+03 8.860000e+02 1.123000e+03
## Bacillus 0.000000e+00 0.000000e+00 0.000000e+00
## Bacteroides 5.152000e+03 1.230000e+03 3.948000e+03
## Bifidobacterium 8.700000e+01 4.600000e+01 1.100000e+02
## Clostridium 0.000000e+00 0.000000e+00 0.000000e+00
## Deinococcus 4.800000e+02 5.580000e+02 5.360000e+02
## Enterococcus 5.765000e+03 1.960800e+04 1.707000e+04
## Escherichia 0.000000e+00 0.000000e+00 0.000000e+00
## Helicobacter 9.730000e+02 6.100000e+01 6.120000e+02
## Lactobacillus 1.782000e+03 7.440000e+02 1.115000e+03
## Listeria 2.619000e+03 2.176000e+03 3.951000e+03
## Neisseria 7.484800e+04 1.100700e+04 7.035500e+04
## Porphyromonas 1.495400e+04 1.411000e+03 1.055500e+04
## Propionibacterium 0.000000e+00 0.000000e+00 0.000000e+00
## Pseudomonas 3.948300e+04 7.028000e+03 3.448100e+04
## Rhodobacter 1.405000e+03 9.350000e+02 1.102000e+03
## Staphylococcus 1.431000e+03 2.770000e+02 1.494000e+03
## Streptococcus 4.640000e+02 8.470000e+02 6.690000e+02
## Others 6.073500e+04 5.055300e+04 7.799300e+04
## sensitivity 7.894737e-01 7.894737e-01 7.894737e-01
## precision_with_NA 7.325757e-01 5.567587e-01 6.776696e-01
## precision_without_NA 8.239576e-01 6.077427e-01 7.447868e-01 指数表記になる原因は、最初の20行分に相当する出現頻度（freq_mat）のほうは～6桁の整数である一方、残りの3行分の感度や適合率のほうは1未満の実数だからです。つまり、全体として表示桁数が大きくなってしまうためです。指数表記を回避したい場合は、以下に示すようにoptions(scipen=100)というコマンドを最初に実行しておけばよいです。
options(scipen=100) # 指数表記を回避
result_pers # result_persの中身を表示## SRR6325813.freque SRR6325815.freque SRR6325816.freque
## Acinetobacter 10653.0000000 10257.0000000 14988.0000000
## Actinomyces 2688.0000000 1868.0000000 3012.0000000
## Bacillus 0.0000000 0.0000000 0.0000000
## Bacteroides 4577.0000000 3624.0000000 6660.0000000
## Bifidobacterium 126.0000000 115.0000000 288.0000000
## Clostridium 0.0000000 0.0000000 0.0000000
## Deinococcus 531.0000000 469.0000000 882.0000000
## Enterococcus 7716.0000000 6148.0000000 25697.0000000
## Escherichia 0.0000000 0.0000000 0.0000000
## Helicobacter 1501.0000000 1379.0000000 1344.0000000
## Lactobacillus 1173.0000000 1058.0000000 1488.0000000
## Listeria 4145.0000000 3485.0000000 5254.0000000
## Neisseria 81414.0000000 81823.0000000 107177.0000000
## Porphyromonas 13869.0000000 12798.0000000 14114.0000000
## Propionibacterium 0.0000000 0.0000000 0.0000000
## Pseudomonas 51798.0000000 48895.0000000 62117.0000000
## Rhodobacter 1589.0000000 1705.0000000 2302.0000000
## Staphylococcus 3105.0000000 2677.0000000 3893.0000000
## Streptococcus 674.0000000 620.0000000 1109.0000000
## Others 73393.0000000 65793.0000000 86215.0000000
## sensitivity 0.7894737 0.7894737 0.7894737
## precision_with_NA 0.7165768 0.7289279 0.7438195
## precision_without_NA 0.8284476 0.8337307 0.8325379
## SRR6325817.freque SRR6325828.freque SRR6325829.freque
## Acinetobacter 51921.0000000 6146.0000000 10441.0000000
## Actinomyces 4212.0000000 835.0000000 1072.0000000
## Bacillus 0.0000000 0.0000000 0.0000000
## Bacteroides 9912.0000000 1057.0000000 630.0000000
## Bifidobacterium 247.0000000 313.0000000 527.0000000
## Clostridium 0.0000000 0.0000000 0.0000000
## Deinococcus 1974.0000000 1133.0000000 2633.0000000
## Enterococcus 10822.0000000 82833.0000000 154386.0000000
## Escherichia 0.0000000 0.0000000 0.0000000
## Helicobacter 395.0000000 245.0000000 164.0000000
## Lactobacillus 2685.0000000 883.0000000 968.0000000
## Listeria 5484.0000000 1156.0000000 865.0000000
## Neisseria 33468.0000000 27625.0000000 11635.0000000
## Porphyromonas 4283.0000000 4177.0000000 2143.0000000
## Propionibacterium 0.0000000 0.0000000 0.0000000
## Pseudomonas 16931.0000000 16186.0000000 5764.0000000
## Rhodobacter 5158.0000000 883.0000000 1216.0000000
## Staphylococcus 1488.0000000 6586.0000000 7184.0000000
## Streptococcus 366.0000000 2648.0000000 5374.0000000
## Others 153638.0000000 39158.0000000 36743.0000000
## sensitivity 0.7894737 0.7894737 0.7894737
## precision_with_NA 0.4929171 0.7959075 0.8480093
## precision_without_NA 0.5581588 0.8735242 0.9259893
## SRR6325830.freque SRR6325831.freque SRR6325859.freque
## Acinetobacter 6512.0000000 7193.0000000 11293.0000000
## Actinomyces 1110.0000000 1156.0000000 2527.0000000
## Bacillus 0.0000000 0.0000000 0.0000000
## Bacteroides 751.0000000 686.0000000 6465.0000000
## Bifidobacterium 269.0000000 106.0000000 165.0000000
## Clostridium 0.0000000 0.0000000 0.0000000
## Deinococcus 1110.0000000 254.0000000 595.0000000
## Enterococcus 90159.0000000 67629.0000000 5505.0000000
## Escherichia 0.0000000 0.0000000 0.0000000
## Helicobacter 408.0000000 296.0000000 1056.0000000
## Lactobacillus 1094.0000000 905.0000000 1460.0000000
## Listeria 1454.0000000 2027.0000000 3202.0000000
## Neisseria 45205.0000000 52283.0000000 88802.0000000
## Porphyromonas 6555.0000000 6804.0000000 18322.0000000
## Propionibacterium 0.0000000 0.0000000 0.0000000
## Pseudomonas 21605.0000000 22188.0000000 49964.0000000
## Rhodobacter 673.0000000 386.0000000 1878.0000000
## Staphylococcus 4487.0000000 2666.0000000 1958.0000000
## Streptococcus 2926.0000000 1625.0000000 528.0000000
## Others 38994.0000000 43671.0000000 59939.0000000
## sensitivity 0.7894737 0.7894737 0.7894737
## precision_with_NA 0.8253833 0.7919190 0.7637025
## precision_without_NA 0.9019059 0.8711495 0.8460054
## SRR6325860.freque SRR6325861.freque SRR6325862.freque
## Acinetobacter 15193.0000000 16686.0000000 16852.0000000
## Actinomyces 1740.0000000 886.0000000 1123.0000000
## Bacillus 0.0000000 0.0000000 0.0000000
## Bacteroides 5152.0000000 1230.0000000 3948.0000000
## Bifidobacterium 87.0000000 46.0000000 110.0000000
## Clostridium 0.0000000 0.0000000 0.0000000
## Deinococcus 480.0000000 558.0000000 536.0000000
## Enterococcus 5765.0000000 19608.0000000 17070.0000000
## Escherichia 0.0000000 0.0000000 0.0000000
## Helicobacter 973.0000000 61.0000000 612.0000000
## Lactobacillus 1782.0000000 744.0000000 1115.0000000
## Listeria 2619.0000000 2176.0000000 3951.0000000
## Neisseria 74848.0000000 11007.0000000 70355.0000000
## Porphyromonas 14954.0000000 1411.0000000 10555.0000000
## Propionibacterium 0.0000000 0.0000000 0.0000000
## Pseudomonas 39483.0000000 7028.0000000 34481.0000000
## Rhodobacter 1405.0000000 935.0000000 1102.0000000
## Staphylococcus 1431.0000000 277.0000000 1494.0000000
## Streptococcus 464.0000000 847.0000000 669.0000000
## Others 60735.0000000 50553.0000000 77993.0000000
## sensitivity 0.7894737 0.7894737 0.7894737
## precision_with_NA 0.7325757 0.5567587 0.6776696
## precision_without_NA 0.8239576 0.6077427 0.7447868ここまででサンプルごとの結果のまとめが完了になります。
次に、W3.6.10を再実行します。
res <- get_W3.6(data, genus, genus_truth) # 属の出現頻度のまとめ情報を取得してresに格納
res # resの中身## $freque
## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135 22229 0 44692
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399 0 11155 493338
## Escherichia Helicobacter Lactobacillus Listeria
## 0 8434 15355 35818
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642 109985 0 376440
## Rhodobacter Staphylococcus Streptococcus Others
## 19232 37246 17850 786825
##
## $TPposi
## [1] TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE
## [25] TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE
## [37] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [85] FALSE FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [97] TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## [133] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [241] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [253] TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [277] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [289] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [301] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [313] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [325] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [337] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [349] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [373] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [385] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [397] FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## [409] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [421] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [433] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [445] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [457] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [469] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [481] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [493] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [505] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [517] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [541] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [553] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [565] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [577] FALSE FALSE FALSE FALSE
##
## $sensitivity
## [1] 0.7894737
##
## $precision_with_NA
## Others
## 0.723414
##
## $precision_without_NA
## [1] 0.8066153 resオブジェクトには、サンプル全体の属レベルの出現頻度、感度、適合率の情報が含まれていることがわかります。これをc関数を用いてW5.1.2と同様のノリで以下のように連結します。
overall <- c(res$freque, # ベクトルとして連結した結果をoverallに格納
res$sensitivity, # ベクトルとして連結した結果をoverallに格納
res$precision_with_NA, # ベクトルとして連結した結果をoverallに格納
res$precision_without_NA) # ベクトルとして連結した結果をoverallに格納
overall # overallの中身## Acinetobacter Actinomyces Bacillus Bacteroides
## 178135.0000000 22229.0000000 0.0000000 44692.0000000
## Bifidobacterium Clostridium Deinococcus Enterococcus
## 2399.0000000 0.0000000 11155.0000000 493338.0000000
## Escherichia Helicobacter Lactobacillus Listeria
## 0.0000000 8434.0000000 15355.0000000 35818.0000000
## Neisseria Porphyromonas Propionibacterium Pseudomonas
## 685642.0000000 109985.0000000 0.0000000 376440.0000000
## Rhodobacter Staphylococcus Streptococcus Others
## 19232.0000000 37246.0000000 17850.0000000 786825.0000000
## Others
## 0.7894737 0.7234140 0.8066153 overallの要素数は、サンプルごとのまとめ情報であるresult_persの行数と同じになります。
length(overall) # overallの要素数## [1] 23dim(result_pers) # result_persの行数と列数## [1] 23 12 サンプル全体の結果であるresから、出現頻度、感度、適合率の情報のみを抽出して1つのベクトルとしてまとめます。そしてそのベクトルをサンプルごとの結果のまとめであるresult_persの右側に列方向で連結します。
result_all <- cbind(result_pers, overall) # result_persとoverallを列方向で連結してresult_allに格納ファイル出力します。第24回のW4.4.3と同様のテクニックを用いています。念のため、ここでもタブ区切りテキストファイル版（output_JSLAB25_genus.txt）とxlsx版（output_JSLAB25_genus.xlsx）をダウンロードできるようにしておきます。
tmp <- cbind(rownames(result_all), result_all)# result_allの行名と、result_allを列方向で連結してtmpに格納
write.table(x = tmp, file = "output_JSLAB25_genus.txt", sep = "\t", row.names = F)# tmpの中身を指定したファイル名でタブ区切りテキストファイルとして、行名をつけずに保存W5.1のスクリプトを以下にまとめます。
source("http://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB_source.R") # 関数の読込
# サンプルごと
result <- apply(X = data, MARGIN = 2, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
FUN = getFreqPerSample, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
genus, genus_truth) # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
hoge <- sapply(X = result, FUN = `[`, 1) # result中の各要素から、1番目の要素の抽出結果をhogeに格納
freq_mat <- sapply(X = hoge, FUN = cbind) # hoge中の各要素を列方向で結合した結果をfreq_matに格納
rownames(freq_mat) <- names(hoge[[1]]) # names(hoge[[1]])をfreq_matの行名情報とする
sensitivity <- unlist(sapply(X = result, FUN = `[`, 3)) # result中の各要素から3番目の要素のみ抽出した結果をベクトル形式にしてsensitivityに格納
precision_with_NA <- unlist(sapply(X = result, FUN = `[`, 4)) # result中の各要素から4番目の要素のみ抽出した結果をベクトル形式にしてprecision_with_NAに格納
precision_without_NA <- unlist(sapply(X = result, FUN = `[`, 5)) # result中の各要素から5番目の要素のみ抽出した結果をベクトル形式にしてprecision_without_NAに格納
result_pers <- rbind(freq_mat, # 行方向で結合した結果をresult_persに格納
sensitivity, # 行方向で結合した結果をresult_persに格納
precision_with_NA, # 行方向で結合した結果をresult_persに格納
precision_without_NA)# 行方向で結合した結果をresult_persに格納
# サンプル全体
res <- get_W3.6(data, genus, genus_truth) # 属の出現頻度のまとめ情報を取得してresに格納
overall <- c(res$freque, # ベクトルとして連結した結果をoverallに格納
res$sensitivity, # ベクトルとして連結した結果をoverallに格納
res$precision_with_NA, # ベクトルとして連結した結果をoverallに格納
res$precision_without_NA) # ベクトルとして連結した結果をoverallに格納
# ファイル出力
result_all <- cbind(result_pers, overall) # result_persとoverallを列方向で連結してresult_allに格納
tmp <- cbind(rownames(result_all), result_all)# result_allの行名と、result_allを列方向で連結してtmpに格納
write.table(x = tmp, file = "output_JSLAB25_genus.txt", sep = "\t", row.names = F)# tmpの中身を指定したファイル名でタブ区切りテキストファイルとして、行名をつけずに保存上記のスクリプトは、W5.1で使ってきたものをそのままペタペタ貼り付けていっただけですが、この後に行う種レベルの解析を意識して、変更可能性のあるオブジェクト名を最初のほうに集約させたものも以下に示しておきます。
source("http://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB_source.R") # 関数の読込
# パラメータ
prediction <- genus # dada2実行結果のオブジェクトをpredictionにコピー
truth <- genus_truth # 正解情報をtruthにコピー
out_f <- "output_JSLAB25_genus.txt" # 出力ファイル名をout_fにコピー
# サンプルごと
result <- apply(X = data, MARGIN = 2, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
FUN = getFreqPerSample, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
prediction, truth) # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
hoge <- sapply(X = result, FUN = `[`, 1) # result中の各要素から、1番目の要素の抽出結果をhogeに格納
freq_mat <- sapply(X = hoge, FUN = cbind) # hoge中の各要素を列方向で結合した結果をfreq_matに格納
rownames(freq_mat) <- names(hoge[[1]]) # names(hoge[[1]])をfreq_matの行名情報とする
sensitivity <- unlist(sapply(X = result, FUN = `[`, 3)) # result中の各要素から3番目の要素のみ抽出した結果をベクトル形式にしてsensitivityに格納
precision_with_NA <- unlist(sapply(X = result, FUN = `[`, 4)) # result中の各要素から4番目の要素のみ抽出した結果をベクトル形式にしてprecision_with_NAに格納
precision_without_NA <- unlist(sapply(X = result, FUN = `[`, 5)) # result中の各要素から5番目の要素のみ抽出した結果をベクトル形式にしてprecision_without_NAに格納
result_pers <- rbind(freq_mat, # 行方向で結合した結果をresult_persに格納
sensitivity, # 行方向で結合した結果をresult_persに格納
precision_with_NA, # 行方向で結合した結果をresult_persに格納
precision_without_NA)# 行方向で結合した結果をresult_persに格納
# サンプル全体
res <- get_W3.6(data, prediction, truth) # 出現頻度のまとめ情報を取得してresに格納
overall <- c(res$freque, # ベクトルとして連結した結果をoverallに格納
res$sensitivity, # ベクトルとして連結した結果をoverallに格納
res$precision_with_NA, # ベクトルとして連結した結果をoverallに格納
res$precision_without_NA) # ベクトルとして連結した結果をoverallに格納
# ファイル出力
result_all <- cbind(result_pers, overall) # result_persとoverallを列方向で連結してresult_allに格納
tmp <- cbind(rownames(result_all), result_all)# result_allの行名と、result_allを列方向で連結してtmpに格納
write.table(x = tmp, file = out_f, sep = "\t", row.names = F)# tmpの中身を指定したファイル名でタブ区切りテキストファイルとして、行名をつけずに保存W5.1.7のスクリプトの最初のほうを少し変更するだけで、種レベルの結果も簡単に得られます。念のため、ここでもタブ区切りテキストファイル版（output_JSLAB25_species.txt）とxlsx版（output_JSLAB25_species.xlsx）をダウンロードできるようにしておきます。
source("http://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB_source.R") # 関数の読込
# パラメータ
prediction <- species # dada2実行結果のオブジェクトをpredictionにコピー
truth <- species_truth # 正解情報をtruthにコピー
out_f <- "output_JSLAB25_species.txt" # 出力ファイル名をout_fにコピー
# サンプルごと
result <- apply(X = data, MARGIN = 2, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
FUN = getFreqPerSample, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
prediction, truth) # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
hoge <- sapply(X = result, FUN = `[`, 1) # result中の各要素から、1番目の要素の抽出結果をhogeに格納
freq_mat <- sapply(X = hoge, FUN = cbind) # hoge中の各要素を列方向で結合した結果をfreq_matに格納
rownames(freq_mat) <- names(hoge[[1]]) # names(hoge[[1]])をfreq_matの行名情報とする
sensitivity <- unlist(sapply(X = result, FUN = `[`, 3)) # result中の各要素から3番目の要素のみ抽出した結果をベクトル形式にしてsensitivityに格納
precision_with_NA <- unlist(sapply(X = result, FUN = `[`, 4)) # result中の各要素から4番目の要素のみ抽出した結果をベクトル形式にしてprecision_with_NAに格納
precision_without_NA <- unlist(sapply(X = result, FUN = `[`, 5)) # result中の各要素から5番目の要素のみ抽出した結果をベクトル形式にしてprecision_without_NAに格納
result_pers <- rbind(freq_mat, # 行方向で結合した結果をresult_persに格納
sensitivity, # 行方向で結合した結果をresult_persに格納
precision_with_NA, # 行方向で結合した結果をresult_persに格納
precision_without_NA)# 行方向で結合した結果をresult_persに格納
# サンプル全体
res <- get_W3.6(data, prediction, truth) # 出現頻度のまとめ情報を取得してresに格納
overall <- c(res$freque, # ベクトルとして連結した結果をoverallに格納
res$sensitivity, # ベクトルとして連結した結果をoverallに格納
res$precision_with_NA, # ベクトルとして連結した結果をoverallに格納
res$precision_without_NA) # ベクトルとして連結した結果をoverallに格納
# ファイル出力
result_all <- cbind(result_pers, overall) # result_persとoverallを列方向で連結してresult_allに格納
tmp <- cbind(rownames(result_all), result_all)# result_allの行名と、result_allを列方向で連結してtmpに格納
write.table(x = tmp, file = out_f, sep = "\t", row.names = F)# tmpの中身を指定したファイル名でタブ区切りテキストファイルとして、行名をつけずに保存 菌株情報ファイル中のスペルミスを修正したうえでW5をやり直します。具体的には、meningitidesをmeningitidisに修正したJSLAB23_table1_corrected.txtを入力ファイルとして実行します。W2.3で得られるオブジェクトとの違いは、species_truthのみになります。
# W2.1 dada2実行結果ファイルのほう
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/output_JSLAB24.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t")# in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
obj <- grep(pattern = "SRR", x = colnames(hoge)) # colnames(hoge)の中から"SRR"という文字を含む要素番号を抽出した結果をobjに格納
data <- hoge[, obj] # hoge行列中のobjで指定した列のみ抽出した結果をdataに格納
genus <- hoge$Genus # hoge中のGenus列を抽出した結果をgenusに格納
species <- hoge$Species # hoge中のSpecies列を抽出した結果をspeciesに格納
# W2.2 菌株情報ファイルのほう
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB23_table1_corrected.txt" # 読み込みたいファイル情報をin_fに格納
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t") #in_fを、1行目を列名、1列目を行名として読み込んだ結果をhogeに格納
organism <- hoge$Organism # hoge中のOrganism列を抽出した結果をorganismに格納
hoge <- strsplit(x = organism, split = ",")# organism中の文字列ベクトルをコンマで分割した結果をhogeに格納
hoge2 <- sapply(X = hoge, FUN = `[`, 1) # hoge中のベクトル形式の各要素から、1番目の要素の抽出結果をhoge2に格納
hoge3 <- strsplit(x = hoge2, split = " ") # hoge2中の文字列ベクトルを半角スペースで分割した結果をhoge3に格納
genus_truth <- sapply(X = hoge3, FUN = `[`, 1) # hoge3中のベクトル形式の各要素から、1番目の要素の抽出結果をgenus_truthに格納
species_truth <- sapply(X = hoge3, FUN = `[`, 2)# hoge3中のベクトル形式の各要素から、2番目の要素の抽出結果をspecies_truthに格納
genus_truth <- unique(genus_truth) # genus_truth中のユニークな要素をgenus_truthに格納
species_truth <- unique(species_truth) # species_truth中のユニークな要素をspecies_truthに格納次に、W5.2に相当する内容を再度実行します。W5.2との違いは、出力ファイル名のみです。念のため、ここでもタブ区切りテキストファイル版（output_JSLAB25_species_corrected.txt）とxlsx版（output_JSLAB25_species_corrected.xlsx）をダウンロードできるようにしておきます
source("http://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB_source.R") # 関数の読込
# パラメータ
prediction <- species # dada2実行結果のオブジェクトをpredictionにコピー
truth <- species_truth # 正解情報をtruthにコピー
out_f <- "output_JSLAB25_species_corrected.txt" # 出力ファイル名をout_fにコピー
# サンプルごと
result <- apply(X = data, MARGIN = 2, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
FUN = getFreqPerSample, # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
prediction, truth) # dataの列ごとに出現頻度のまとめ情報を取得した結果をresultに格納
hoge <- sapply(X = result, FUN = `[`, 1) # result中の各要素から、1番目の要素の抽出結果をhogeに格納
freq_mat <- sapply(X = hoge, FUN = cbind) # hoge中の各要素を列方向で結合した結果をfreq_matに格納
rownames(freq_mat) <- names(hoge[[1]]) # names(hoge[[1]])をfreq_matの行名情報とする
sensitivity <- unlist(sapply(X = result, FUN = `[`, 3)) # result中の各要素から3番目の要素のみ抽出した結果をベクトル形式にしてsensitivityに格納
precision_with_NA <- unlist(sapply(X = result, FUN = `[`, 4)) # result中の各要素から4番目の要素のみ抽出した結果をベクトル形式にしてprecision_with_NAに格納
precision_without_NA <- unlist(sapply(X = result, FUN = `[`, 5)) # result中の各要素から5番目の要素のみ抽出した結果をベクトル形式にしてprecision_without_NAに格納
result_pers <- rbind(freq_mat, # 行方向で結合した結果をresult_persに格納
sensitivity, # 行方向で結合した結果をresult_persに格納
precision_with_NA, # 行方向で結合した結果をresult_persに格納
precision_without_NA)# 行方向で結合した結果をresult_persに格納
# サンプル全体
res <- get_W3.6(data, prediction, truth) # 出現頻度のまとめ情報を取得してresに格納
overall <- c(res$freque, # ベクトルとして連結した結果をoverallに格納
res$sensitivity, # ベクトルとして連結した結果をoverallに格納
res$precision_with_NA, # ベクトルとして連結した結果をoverallに格納
res$precision_without_NA) # ベクトルとして連結した結果をoverallに格納
# ファイル出力
result_all <- cbind(result_pers, overall) # result_persとoverallを列方向で連結してresult_allに格納
tmp <- cbind(rownames(result_all), result_all)# result_allの行名と、result_allを列方向で連結してtmpに格納
write.table(x = tmp, file = out_f, sep = "\t", row.names = F)# tmpの中身を指定したファイル名でタブ区切りテキストファイルとして、行名をつけずに保存