---

title: "3. Filtering"

description: |

Reproducible workflow for Arbitrary, PERfect, & PIME filtering of community data.

format:

html:

toc: true

toc-depth: 3

---

```{r}

#| message: false

#| results: hide

#| code-fold: true

#| code-summary: "Click here for setup information"

remove(list = ls())

knitr::opts_chunk$set(echo = TRUE, eval = FALSE)

set.seed(119)

#library(conflicted)

#pacman::p_depends(PERFect, local = TRUE)

#pacman::p_depends_reverse(PERFect, local = TRUE)

library(phyloseq); packageVersion("phyloseq")

library(Biostrings); packageVersion("Biostrings")

pacman::p_load(tidyverse, patchwork,

agricolae, labdsv, naniar, PERFect, pime, ape, gdata,

pairwiseAdonis, microbiome, seqRFLP, microbiomeMarker,

reactable, downloadthis, captioner, jamba,

install = FALSE, update = FALSE)

library(ggstatsplot)

library(ggside)

options(scipen=999)

knitr::opts_current$get(c(

"cache",

"cache.path",

"cache.rebuild",

"dependson",

"autodep"

))

source(file.path("assets", "functions.R"))

```

```{r}

#| include: false

#| eval: true

## Load to build page

load("page_build/filtering_wf.rdata")

ssu18_ps_work <- readRDS("files/data-prep/rdata/ssu18_ps_work.rds")

its18_ps_work <- readRDS("files/data-prep/rdata/its18_ps_work.rds")

```

```{r}

#| include: false

#| eval: true

#| results: hide

# Create the caption(s) with captioner

caption_tab_ssu <- captioner(prefix = "(16S rRNA) Table", suffix = " |", style = "b")

caption_fig_ssu <- captioner(prefix = "(16S rRNA) Figure", suffix = " |", style = "b")

caption_tab_its <- captioner(prefix = "(ITS) Table", suffix = " |", style = "b")

caption_fig_its <- captioner(prefix = "(ITS) Figure", suffix = " |", style = "b")

# Create a function for referring to the tables in text

ref <- function(x) str_extract(x, "[^|]*") %>%

trimws(which = "right", whitespace = "[ ]")

```

</details>

```{r}

#| echo: false

#| eval: true

source("assets/captions/captions_filtering.R")

```

In [Part A](#a.-arbitrary-filtering), we apply arbitrary filtering to the 16S rRNA and ITS data sets. In [Part B](#b.-perfect-filtering) we use PERFect [(PERmutation Filtering test for microbiome data)](https://github.com/katiasmirn/PERFect) [@smirnova2019perfect] to filter the data sets. And in [Part C](#c.-pime-filtering) of this workflow, we use PIME [(Prevalence Interval for Microbiome Evaluation)](https://github.com/microEcology/pime) [@roesch2020pime] to filter the FULL 16S rRNA and ITS data sets.

## Workflow Input

Files needed to run this workflow can be downloaded from figshare.

<iframe src="https://widgets.figshare.com/articles/14690739/embed?show_title=1" width="100%" height="301" allowfullscreen frameborder="0"></iframe>

# Filtering Results

We begin by summarizing the results of each filtering method on the 16S and ITS data sets. Below you can find the complete workflow for each filtering method.

## 16S rRNA

<small>`r caption_tab_ssu("ssu_all_filt_summary")`</small>

```{r}

#| echo: false

ssu18_filtering_results <- dplyr::left_join(ssu18_sum_arbitrary, ssu18_sum_perfect,

by = c("Metric", "Full data")) %>%

dplyr::left_join(., ssu18_sum_pime,

by = c("Metric", "Full data")) %>%

dplyr::rename("Arbitrary" = 3, "PERfect" = 4, "PIME" = 5)

```

```{r}

#| echo: false

#| eval: true

seq_table <- ssu18_filtering_results

write.table(seq_table, "files/filtering/tables/ssu_all_filt_summary.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = FALSE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 100

),

columns = list(

Metric = colDef(name = "Metric",

sticky = "left",

style = list(borderRight = "1px solid #eee"),

headerStyle = list(borderRight = "1px solid #eee"),

align = "left",

minWidth = 200)

),

searchable = FALSE, defaultPageSize = 10, showPagination = FALSE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = TRUE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = FALSE, showSortable = FALSE, fullWidth = TRUE,

theme = reactableTheme(style = list(fontSize = "0.8em")))

rm(seq_table)

```

{{< downloadthis files/filtering/tables/ssu_all_filt_summary.txt dname="ssu_all_filt_summary" label="Download table as text file" icon="table" type="primary" >}}

```{r}

#| echo: false

tmp_read_count <- data.frame(sample_sums(ssu18_ps_work)) %>%

tibble::rownames_to_column("SamName")

tmp_asv_count <- estimate_richness(ssu18_ps_work, split = TRUE, measures = "Observed") %>%

tibble::rownames_to_column("SamName")

tmp_samp_data <- data.frame(sample_data(ssu18_ps_work))

tmp_samp_data <- tmp_samp_data[order(tmp_samp_data$PLOT, tmp_samp_data$TEMP),]

tmp_join <- dplyr::left_join(tmp_samp_data, tmp_read_count, by = "SamName") %>%

dplyr::left_join(., tmp_asv_count, by = "SamName") %>%

dplyr::rename("Sample_ID" = 1) %>%

dplyr::rename("total_reads" = 7) %>%

dplyr::rename("total_asvs" = 8)

tmp_join[, 3] <- NULL

ssu_work_samp_summary <- tmp_join

rm(list = ls(pattern = "tmp_"))

```

```{r}

#| echo: false

ssu18_filtering_results <- dplyr::left_join(ssu18_sum_arbitrary, ssu18_sum_perfect,

by = c("Metric", "Full data")) %>%

dplyr::left_join(., ssu18_sum_pime,

by = c("Metric", "Full data")) %>%

dplyr::rename("Arbitrary" = 3, "PERfect" = 4, "PIME" = 5)

```

```{r}

#| echo: false

ssu_all_filt_samp_summary <- dplyr::left_join(ssu_work_samp_summary, ssu_filt_samp_summary,

by = c("Sample_ID", "PLOT", "TREAT", "TEMP", "PAIR")) %>%

dplyr::left_join(., ssu_perfect_samp_summary,

by = c("Sample_ID", "PLOT", "TREAT", "TEMP", "PAIR")) %>%

dplyr::left_join(., ssu_pime_samp_summary,

by = c("Sample_ID", "PLOT", "TREAT", "TEMP", "PAIR")) %>%

dplyr::relocate(c("total_reads.x", "total_reads.y",

"total_reads.x.x", "total_reads.y.y"), .after = "PAIR")

ssu_all_filt_samp_summary[,2:5] <- NULL

```

<small>`r caption_tab_ssu("ssu_all_filt_samp_summary")`</small>

```{r}

#| echo: false

#| eval: true

seq_table <- ssu_all_filt_samp_summary

write.table(seq_table, "files/filtering/tables/ssu_all_filt_samp_summary.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

columns = list(

Sample_ID = colDef(name = "Sample_ID",

sticky = "left",

style = list(borderRight = "5px solid #eee"),

headerStyle = list(borderRight = "5px solid #eee"),

align = "left",

minWidth = 150),

total_reads.x = colDef(name = "FULL"),

total_reads.y = colDef(name = "Arbitrary"),

total_reads.x.x = colDef(name = "PERfect"),

total_reads.y.y = colDef(name = "PIME",

style = list(borderRight = "5px solid #eee"),

headerStyle = list(borderRight = "5px solid #eee")),

total_asvs.x = colDef(name = "FULL"),

total_asvs.y = colDef(name = "Arbitrary"),

total_asvs.x.x = colDef(name = "PERfect"),

total_asvs.y.y = colDef(name = "PIME")

),

columnGroups = list(

colGroup(name = "Total Reads", columns = c("total_reads.x", "total_reads.y", "total_reads.x.x", "total_reads.y.y"),

headerStyle = list(borderRight = "5px solid #eee", fontSize = "1.5em")),

colGroup(name = "Total ASVs", columns = c("total_asvs.x", "total_asvs.y", "total_asvs.x.x", "total_asvs.y.y"),

headerStyle = list(fontSize = "1.5em"))

),

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = TRUE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 100

),

searchable = TRUE, defaultPageSize = 5, showPagination = TRUE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = TRUE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = FALSE, showSortable = TRUE, fullWidth = TRUE,

theme = reactableTheme(style = list(fontSize = "0.8em")))

rm(seq_table)

```

{{< downloadthis files/filtering/tables/ssu_all_filt_samp_summary.txt dname="ssu_all_filt_samp_summary" label="Download table as text file" icon="table" type="primary" >}}

## ITS

<small>`r caption_tab_its("its_all_filt_summary")`</small>

```{r}

#| echo: false

its18_filtering_results <- dplyr::left_join(its18_sum_arbitrary, its18_sum_perfect,

by = c("Metric", "Full data")) %>%

dplyr::left_join(., its18_sum_pime,

by = c("Metric", "Full data")) %>%

dplyr::rename("Arbitrary" = 3, "PERfect" = 4, "PIME" = 5)

```

```{r}

#| echo: false

#| eval: true

seq_table <- its18_filtering_results

write.table(seq_table, "files/filtering/tables/its_all_filt_summary.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = FALSE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 100

),

columns = list(

Metric = colDef(name = "Metric",

sticky = "left",

style = list(borderRight = "1px solid #eee"),

headerStyle = list(borderRight = "1px solid #eee"),

align = "left",

minWidth = 200)

),

searchable = FALSE, defaultPageSize = 10, showPagination = FALSE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = TRUE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = FALSE, showSortable = FALSE, fullWidth = TRUE,

theme = reactableTheme(style = list(fontSize = "0.8em")))

rm(seq_table)

```

{{< downloadthis files/filtering/tables/its_all_filt_summary.txt dname="its_all_filt_summary" label="Download table as text file" icon="table" type="primary" >}}

```{r}

#| echo: false

tmp_read_count <- data.frame(sample_sums(its18_ps_work)) %>%

tibble::rownames_to_column("SamName")

tmp_asv_count <- estimate_richness(its18_ps_work, split = TRUE, measures = "Observed") %>%

tibble::rownames_to_column("SamName")

tmp_samp_data <- data.frame(sample_data(its18_ps_work))

tmp_samp_data <- tmp_samp_data[order(tmp_samp_data$PLOT, tmp_samp_data$TEMP),]

tmp_join <- dplyr::left_join(tmp_samp_data, tmp_read_count, by = "SamName") %>%

dplyr::left_join(., tmp_asv_count, by = "SamName") %>%

dplyr::rename("Sample_ID" = 1) %>%

dplyr::rename("total_reads" = 7) %>%

dplyr::rename("total_asvs" = 8)

tmp_join[, 3] <- NULL

its_work_samp_summary <- tmp_join

rm(list = ls(pattern = "tmp_"))

```

```{r}

#| echo: false

its18_filtering_results <- dplyr::left_join(its18_sum_arbitrary, its18_sum_perfect,

by = c("Metric", "Full data")) %>%

dplyr::left_join(., its18_sum_pime,

by = c("Metric", "Full data")) %>%

dplyr::rename("Arbitrary" = 3, "PERfect" = 4, "PIME" = 5)

```

```{r}

#| echo: false

its_all_filt_samp_summary <- dplyr::left_join(its_work_samp_summary, its_filt_samp_summary,

by = c("Sample_ID", "PLOT", "TREAT", "TEMP", "PAIR")) %>%

dplyr::left_join(., its_perfect_samp_summary,

by = c("Sample_ID", "PLOT", "TREAT", "TEMP", "PAIR")) %>%

dplyr::left_join(., its_pime_samp_summary,

by = c("Sample_ID", "PLOT", "TREAT", "TEMP", "PAIR")) %>%

dplyr::relocate(c("total_reads.x", "total_reads.y",

"total_reads.x.x", "total_reads.y.y"), .after = "PAIR")

its_all_filt_samp_summary[,2:5] <- NULL

```

<small>`r caption_tab_its("its_all_filt_samp_summary")`</small>

```{r}

#| echo: false

#| eval: true

seq_table <- its_all_filt_samp_summary

write.table(seq_table, "files/filtering/tables/its_all_filt_samp_summary.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

columns = list(

Sample_ID = colDef(name = "Sample_ID",

sticky = "left",

style = list(borderRight = "5px solid #eee"),

headerStyle = list(borderRight = "5px solid #eee"),

align = "left",

minWidth = 150),

total_reads.x = colDef(name = "FULL"),

total_reads.y = colDef(name = "Arbitrary"),

total_reads.x.x = colDef(name = "PERfect"),

total_reads.y.y = colDef(name = "PIME",

style = list(borderRight = "5px solid #eee"),

headerStyle = list(borderRight = "5px solid #eee")),

total_asvs.x = colDef(name = "FULL"),

total_asvs.y = colDef(name = "Arbitrary"),

total_asvs.x.x = colDef(name = "PERfect"),

total_asvs.y.y = colDef(name = "PIME")

),

columnGroups = list(

colGroup(name = "Total Reads", columns = c("total_reads.x", "total_reads.y", "total_reads.x.x", "total_reads.y.y"),

headerStyle = list(borderRight = "5px solid #eee", fontSize = "1.5em")),

colGroup(name = "Total ASVs", columns = c("total_asvs.x", "total_asvs.y", "total_asvs.x.x", "total_asvs.y.y"),

headerStyle = list(fontSize = "1.5em"))

),

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = TRUE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 100

),

searchable = TRUE, defaultPageSize = 5, showPagination = TRUE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = TRUE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = FALSE, showSortable = TRUE, fullWidth = TRUE,

theme = reactableTheme(style = list(fontSize = "0.8em")))

rm(seq_table)

```

{{< downloadthis files/filtering/tables/its_all_filt_samp_summary.txt dname="its_all_filt_samp_summary" label="Download table as text file" icon="table" type="primary" >}}

# A. Arbitrary Filtering

For low-count arbitrary filtering, we set the minimum read count to `5` and the prevalence to 20%. Then we apply another filter to remove ASVs with a low variance (`0.2`).

```{r}

#| include: false

## Initial Load for ARBITRARY ANALYSIS #1

set.seed(119)

ssu18_ps_work <- readRDS("files/data-prep/rdata/ssu18_ps_work.rds")

its18_ps_work <- readRDS("files/data-prep/rdata/its18_ps_work.rds")

```

## 16S rRNA

```{r}

#| code-fold: true

tmp_low_count <- phyloseq::genefilter_sample(

ssu18_ps_work,

filterfun_sample(function(x) x >= 5),

A = 0.2 * nsamples(ssu18_ps_work))

tmp_low_count <- phyloseq::prune_taxa(tmp_low_count, ssu18_ps_work)

tmp_low_var <- phyloseq::filter_taxa(tmp_low_count,

function(x) var(x) > 0.2, prune = TRUE)

ssu18_ps_filt <- tmp_low_var

ssu18_ps_filt@phy_tree <- NULL

tmp_tree <- rtree(ntaxa(ssu18_ps_filt), rooted = TRUE,

tip.label = taxa_names(ssu18_ps_filt))

ssu18_ps_filt <- merge_phyloseq(ssu18_ps_filt,

sample_data,

tmp_tree)

rm(list = ls(pattern = "tmp_"))

ssu18_ps_filt

```

Here is the filtered 16S rRNA phyloseq object.

```{r}

#| echo: false

#| eval: true

ssu18_ps_filt

```

```{r}

#| echo: false

## This code is gross. I spent hours trying to find a better

## solution but failed :(

rm(list = ls(pattern = "tmp_"))

tmp_samples <- c("ssu18_ps_work", "ssu18_ps_filt")

tmp_objects <- c(

"Total number of reads",

"Min. number of reads",

"Max. number of reads",

"Average number of reads",

"Median number of reads",

"Total ASVs",

"Min. number of ASVs",

"Max. number of ASVs",

"Average number of ASVs",

"Median number of ASVs")

tmp_total_reads <- c()

for (i in tmp_samples) {

tmp_get <- sum(readcount(get(i)))

tmp_total_reads <- c(append(tmp_total_reads, tmp_get))

}

tmp_min_read <- c()

for (i in tmp_samples) {

tmp_get <- min(readcount(get(i)))

tmp_min_read <- c(append(tmp_min_read, tmp_get))

}

tmp_max_read <- c()

for (i in tmp_samples) {

tmp_get <- max(readcount(get(i)))

tmp_max_read <- c(append(tmp_max_read, tmp_get))

}

tmp_mean_reads <- c()

for (i in tmp_samples) {

tmp_get <- round(mean(readcount(get(i))), digits = 0)

tmp_mean_reads <- c(append(tmp_mean_reads, tmp_get))

}

tmp_median_reads <- c()

for (i in tmp_samples) {

tmp_get <- median(readcount(get(i)))

tmp_median_reads <- c(append(tmp_median_reads, tmp_get))

}

tmp_total_asvs <- c()

for (i in tmp_samples) {

tmp_get <- ntaxa(get(i))

tmp_total_asvs <- c(append(tmp_total_asvs, tmp_get))

}

tmp_min_asvs <- c()

for (i in tmp_samples) {

tmp_get <- min(richness(get(i), index = "observed"))

tmp_min_asvs <- c(append(tmp_min_asvs, tmp_get))

}

tmp_max_asvs <- c()

for (i in tmp_samples) {

tmp_get <- max(richness(get(i), index = "observed"))

tmp_max_asvs <- c(append(tmp_max_asvs, tmp_get))

}

tmp_mean_asvs <- c()

for (i in tmp_samples) {

tmp_get <- round(mean(richness(get(i), index = "observed")$observed), digits = 0)

tmp_mean_asvs <- c(append(tmp_mean_asvs, tmp_get))

}

tmp_med_asvs <- c()

for (i in tmp_samples) {

tmp_get <- median(richness(get(i), index = "observed")$observed)

tmp_med_asvs <- c(append(tmp_med_asvs, tmp_get))

}

tmp_bind <- data.frame(do.call("rbind", list(

tmp_total_reads, tmp_min_read, tmp_max_read, tmp_mean_reads, tmp_median_reads,

tmp_total_asvs, tmp_min_asvs, tmp_max_asvs, tmp_mean_asvs, tmp_med_asvs)))

ssu18_sum_arbitrary <- dplyr::bind_cols(tmp_objects, tmp_bind) %>%

dplyr::rename("Metric" = 1, "Full data" = 2, "Arbitrary filter" = 3)

rm(list = ls(pattern = "tmp_"))

```

<small>`r caption_tab_ssu("ssu_arbitrary_filt_summary")`</small>

```{r}

#| echo: false

#| eval: true

seq_table <- ssu18_sum_arbitrary

write.table(seq_table, "files/filtering/tables/ssu_arbitrary_filt_summary.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = FALSE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 150

),

columns = list(

Metric = colDef(name = "Metric",

sticky = "left",

style = list(borderRight = "1px solid #eee"),

headerStyle = list(borderRight = "1px solid #eee"),

align = "left",

minWidth = 200)

),

searchable = FALSE, defaultPageSize = 10, showPagination = FALSE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = TRUE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = FALSE, showSortable = FALSE, fullWidth = TRUE,

theme = reactableTheme(style = list(fontSize = "0.8em")))

rm(seq_table)

```

{{< downloadthis files/filtering/tables/ssu_arbitrary_filt_summary.txt dname="ssu_arbitrary_filt_summary" label="Download table as text file" icon="table" type="primary" >}}

We can look at how filtering affected total reads and ASVs for each sample.

```{r}

#| echo: false

tmp_read_count <- data.frame(sample_sums(ssu18_ps_filt)) %>%

tibble::rownames_to_column("SamName")

tmp_asv_count <- estimate_richness(ssu18_ps_filt, split = TRUE, measures = "Observed") %>%

tibble::rownames_to_column("SamName")

tmp_samp_data <- data.frame(sample_data(ssu18_ps_filt))

tmp_samp_data <- tmp_samp_data[order(tmp_samp_data$PLOT, tmp_samp_data$TEMP),]

tmp_join <- dplyr::left_join(tmp_samp_data, tmp_read_count, by = "SamName") %>%

dplyr::left_join(., tmp_asv_count, by = "SamName") %>%

dplyr::rename("Sample_ID" = 1) %>%

dplyr::rename("total_reads" = 7) %>%

dplyr::rename("total_asvs" = 8)

tmp_join[, 3] <- NULL

ssu_filt_samp_summary <- tmp_join

rm(list = ls(pattern = "tmp_"))

```

<small>`r caption_tab_ssu("ssu_filt_samp_summary")`</small>

```{r}

#| echo: false

#| eval: true

seq_table <- ssu_filt_samp_summary

write.table(seq_table, "files/filtering/tables/ssu_filt_samp_summary.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = TRUE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 100

),

columns = list(

Sample_ID = colDef(name = "Sample_ID",

sticky = "left",

style = list(borderRight = "1px solid #eee"),

headerStyle = list(borderRight = "1px solid #eee"),

align = "left",

minWidth = 150)

),

searchable = TRUE, defaultPageSize = 5, showPagination = TRUE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = TRUE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = FALSE, showSortable = TRUE, fullWidth = FALSE,

theme = reactableTheme(style = list(fontSize = "0.8em")))

rm(seq_table)

```

{{< downloadthis files/filtering/tables/ssu_filt_samp_summary.txt dname="ssu_filt_samp_summary" label="Download table as text file" icon="table" type="primary" >}}

## ITS

```{r}

#| code-fold: true

tmp_low_count <- phyloseq::genefilter_sample(

its18_ps_work,

filterfun_sample(function(x) x >= 5),

A = 0.2 * nsamples(its18_ps_work))

tmp_low_count <- phyloseq::prune_taxa(tmp_low_count, its18_ps_work)

tmp_low_var <- phyloseq::filter_taxa(tmp_low_count,

function(x) var(x) > 0.2, prune = TRUE)

its18_ps_filt <- tmp_low_var

rm(list = ls(pattern = "tmp_"))

its18_ps_filt

```

And the filtered ITS phyloseq object.

```{r}

#| echo: false

#| eval: true

its18_ps_filt

```

```{r}

#| echo: false

## This code is gross. I spent hours trying to find a better

## solution but failed :(

rm(list = ls(pattern = "tmp_"))

tmp_samples <- c("its18_ps_work", "its18_ps_filt")

tmp_objects <- c(

"Total number of reads",

"Min. number of reads",

"Max. number of reads",

"Average number of reads",

"Median number of reads",

"Total ASVs",

"Min. number of ASVs",

"Max. number of ASVs",

"Average number of ASVs",

"Median number of ASVs")

tmp_total_reads <- c()

for (i in tmp_samples) {

tmp_get <- sum(readcount(get(i)))

tmp_total_reads <- c(append(tmp_total_reads, tmp_get))

}

tmp_min_read <- c()

for (i in tmp_samples) {

tmp_get <- min(readcount(get(i)))

tmp_min_read <- c(append(tmp_min_read, tmp_get))

}

tmp_max_read <- c()

for (i in tmp_samples) {

tmp_get <- max(readcount(get(i)))

tmp_max_read <- c(append(tmp_max_read, tmp_get))

}

tmp_mean_reads <- c()

for (i in tmp_samples) {

tmp_get <- round(mean(readcount(get(i))), digits = 0)

tmp_mean_reads <- c(append(tmp_mean_reads, tmp_get))

}

tmp_median_reads <- c()

for (i in tmp_samples) {

tmp_get <- median(readcount(get(i)))

tmp_median_reads <- c(append(tmp_median_reads, tmp_get))

}

tmp_total_asvs <- c()

for (i in tmp_samples) {

tmp_get <- ntaxa(get(i))

tmp_total_asvs <- c(append(tmp_total_asvs, tmp_get))

}

tmp_min_asvs <- c()

for (i in tmp_samples) {

tmp_get <- min(richness(get(i), index = "observed"))

tmp_min_asvs <- c(append(tmp_min_asvs, tmp_get))

}

tmp_max_asvs <- c()

for (i in tmp_samples) {

tmp_get <- max(richness(get(i), index = "observed"))

tmp_max_asvs <- c(append(tmp_max_asvs, tmp_get))

}

tmp_mean_asvs <- c()

for (i in tmp_samples) {

tmp_get <- round(mean(richness(get(i), index = "observed")$observed), digits = 0)

tmp_mean_asvs <- c(append(tmp_mean_asvs, tmp_get))

}

tmp_med_asvs <- c()

for (i in tmp_samples) {

tmp_get <- median(richness(get(i), index = "observed")$observed)

tmp_med_asvs <- c(append(tmp_med_asvs, tmp_get))

}

tmp_bind <- data.frame(do.call("rbind", list(

tmp_total_reads, tmp_min_read, tmp_max_read, tmp_mean_reads, tmp_median_reads,

tmp_total_asvs, tmp_min_asvs, tmp_max_asvs, tmp_mean_asvs, tmp_med_asvs)))

its18_sum_arbitrary <- dplyr::bind_cols(tmp_objects, tmp_bind) %>%

dplyr::rename("Metric" = 1, "Full data" = 2, "Arbitrary filter" = 3)

rm(list = ls(pattern = "tmp_"))

```

<small>`r caption_tab_its("its_arbitrary_filt_summary")`</small>

```{r}

#| echo: false

#| eval: true

seq_table <- its18_sum_arbitrary

write.table(seq_table, "files/filtering/tables/its_arbitrary_filt_summary.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = FALSE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 150

),

columns = list(

Metric = colDef(name = "Metric",

sticky = "left",

style = list(borderRight = "1px solid #eee"),

headerStyle = list(borderRight = "1px solid #eee"),

align = "left",

minWidth = 200)

),

searchable = FALSE, defaultPageSize = 10, showPagination = FALSE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = FALSE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = FALSE, showSortable = FALSE, fullWidth = TRUE,

theme = reactableTheme(style = list(fontSize = "0.8em")))

rm(seq_table)

```

{{< downloadthis files/filtering/tables/its_arbitrary_filt_summary.txt dname="its_arbitrary_filt_summary" label="Download table as text file" icon="table" type="primary" >}}

And again, let's look how filtering affected total reads and ASVs for each sample.

```{r}

#| echo: false

tmp_read_count <- data.frame(sample_sums(its18_ps_filt)) %>%

tibble::rownames_to_column("SamName")

tmp_asv_count <- estimate_richness(its18_ps_filt, split = TRUE, measures = "Observed") %>%

tibble::rownames_to_column("SamName")

tmp_samp_data <- data.frame(sample_data(its18_ps_filt))

tmp_samp_data <- tmp_samp_data[order(tmp_samp_data$PLOT, tmp_samp_data$TEMP),]

tmp_join <- dplyr::left_join(tmp_samp_data, tmp_read_count, by = "SamName") %>%

dplyr::left_join(., tmp_asv_count, by = "SamName") %>%

dplyr::rename("Sample_ID" = 1) %>%

dplyr::rename("total_reads" = 7) %>%

dplyr::rename("total_asvs" = 8)

tmp_join[, 3] <- NULL

its_filt_samp_summary <- tmp_join

rm(list = ls(pattern = "tmp_"))

```

<small>`r caption_tab_its("its_filt_samp_summary")`</small>

```{r}

#| echo: false

#| eval: true

seq_table <- its_filt_samp_summary

write.table(seq_table, "files/filtering/tables/its_filt_samp_summary.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = TRUE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 100

),

columns = list(

Sample_ID = colDef(name = "Sample_ID",

sticky = "left",

style = list(borderRight = "1px solid #eee"),

headerStyle = list(borderRight = "1px solid #eee"),

align = "left",

minWidth = 150)

),

searchable = TRUE, defaultPageSize = 5, showPagination = TRUE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = TRUE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = FALSE, showSortable = TRUE, fullWidth = FALSE,

theme = reactableTheme(style = list(fontSize = "0.8em")))

rm(seq_table)

```

{{< downloadthis files/filtering/tables/its_filt_samp_summary.txt dname="its_filt_samp_summary" label="Download table as text file" icon="table" type="primary" >}}

```{r}

#| include: false

saveRDS(ssu18_ps_filt, "files/filtering/arbitrary/rdata/ssu18_ps_filt.rds")

saveRDS(its18_ps_filt, "files/filtering/arbitrary/rdata/its18_ps_filt.rds")

```

And that's it for Arbitrary filtering. Moving on.

# B. PERfect Filtering

To run PERFect, we need the ASV tables in data frame format with samples as rows and ASVs as columns. PERFect is sensitive to the order of ASVs, so here we test **a**) the default order in the phyloseq object and **b**) ASVs ordered by decreasing abundance.

## Setup

```{r}

#| code-fold: true

samp_ps_filt <- c("ssu18_ps_work", "its18_ps_work")

for (i in samp_ps_filt) {

tmp_get <- get(i)

tmp_get_tab <- data.frame(t(otu_table(tmp_get)))

tmp_get_tab <- tmp_get_tab %>% tibble::rownames_to_column("ID")

tmp_get_tab <- jamba::mixedSortDF(tmp_get_tab, decreasing = FALSE,

useRownames = FALSE, byCols = 1)

tmp_get_tab <- tmp_get_tab %>% tibble::remove_rownames()

tmp_get_tab <- tmp_get_tab %>% tibble::column_to_rownames("ID")

tmp_get_tab <- data.frame(t(tmp_get_tab))

tmp_get_tab_ord <- data.frame(t(otu_table(tmp_get)))

tmp_get_tab_ord <- tmp_get_tab_ord %>% tibble::rownames_to_column("ID")

tmp_get_tab_ord <- jamba::mixedSortDF(tmp_get_tab_ord, decreasing = TRUE,

useRownames = FALSE, byCols = 1)

tmp_get_tab_ord <- tmp_get_tab_ord %>% tibble::remove_rownames()

tmp_get_tab_ord <- tmp_get_tab_ord %>% tibble::column_to_rownames("ID")

tmp_get_tab_ord <- data.frame(t(tmp_get_tab_ord))

tmp_tab_name <- purrr::map_chr(i, ~ paste0(., "_perfect"))

assign(tmp_tab_name, tmp_get_tab)

tmp_tab_ord_name <- purrr::map_chr(i, ~ paste0(., "_ord_perfect"))

assign(tmp_tab_ord_name, tmp_get_tab_ord)

rm(list = ls(pattern = "tmp_"))

}

objects()

```

## Filter

Next we run the filtering analysis using [`PERFect_sim`](https://rdrr.io/github/katiasmirn/PERFect/man/PERFect_sim.html). The other option is [`PERFect_perm`](https://rdrr.io/github/katiasmirn/PERFect/man/PERFect_perm.html) however I could not get `PERFect_perm` to work as of this writing. The process never finished :/

We need to set an initial p-value cutoff. For 16S rRNA, we use `0.05` and for ITS we use `0.1`.

```{r}

# Set a pvalue cutoff

ssu_per_pval <- 0.05

its_per_pval <- 0.10

```

```{r}

#| code-fold: true

for (i in samp_ps_filt) {

tmp_get <- get(purrr::map_chr(i, ~ paste0(., "_perfect")))

tmp_pval <- gsub("18.*", "_per_pval", i)

tmp_get_ord <- get(purrr::map_chr(i, ~ paste0(., "_ord_perfect")))

tmp_sim <- PERFect_sim(X = tmp_get, alpha = get(tmp_pval), Order = "NP", center = FALSE)

dim(tmp_sim$filtX)

tmp_sim_ord <- PERFect_sim(X = tmp_get_ord, alpha = get(tmp_pval), Order = "NP", center = FALSE)

dim(tmp_sim_ord$filtX)

tmp_sim_name <- purrr::map_chr(i, ~ paste0(., "_perfect_sim"))

assign(tmp_sim_name, tmp_sim)

tmp_sim_ord_name <- purrr::map_chr(i, ~ paste0(., "_ord_perfect_sim"))

assign(tmp_sim_ord_name, tmp_sim_ord)

tmp_path <- file.path("files/filtering/perfect/rdata/")

saveRDS(tmp_sim, paste(tmp_path, tmp_sim_name, ".rds", sep = ""))

saveRDS(tmp_sim_ord, paste(tmp_path, tmp_sim_ord_name, ".rds", sep = ""))

rm(list = ls(pattern = "tmp_"))

}

objects(pattern = "_sim")

```

```{r}

#| echo: false

## FOR TESTING

dim(ssu18_ps_work_perfect_sim$filtX)

dim(ssu18_ps_work_ord_perfect_sim$filtX)

dim(its18_ps_work_perfect_sim$filtX)

dim(its18_ps_work_ord_perfect_sim$filtX)

###########

tmp_obj <- ssu18_ps_work_ord_perfect_sim

tmp_ps <- ssu18_ps_work

###########

tmp_filt <- data.frame(t(tmp_obj$filtX))

tmp_filt <- tmp_filt %>% tibble::rownames_to_column("ID")

tmp_pval <- data.frame(tmp_obj$pvals)

tmp_pval <- tmp_pval %>% tibble::rownames_to_column("ID")

tmp_org <- data.frame(t(otu_table(tmp_ps)))

tmp_org <- tmp_org %>% tibble::rownames_to_column("ID")

tmp_com <- dplyr::left_join(tmp_pval, tmp_org, by = "ID") %>%

dplyr::left_join(., tmp_filt, by = "ID")

write.table(tmp_com, "tmp_com.txt", quote = FALSE, sep = "\t", row.names = FALSE)

```

```{r}

#| echo: false

ps_sim <- c("ssu18_ps_work_perfect_sim", "ssu18_ps_work_ord_perfect_sim",

"its18_ps_work_perfect_sim", "its18_ps_work_ord_perfect_sim", )

for (i in ps_sim) {

tmp_get <- get(i)

tmp_name <- i

tmp_get$info <- NULL

tmp_get$DFL <- NULL

tmp_get$est <- NULL

tmp_get$fit <- NULL

tmp_get$hist <- NULL

tmp_get$pDFL <- NULL

assign(tmp_name, tmp_get)

rm(list = ls(pattern = "tmp_"))

}

objects()

```

```{r}

#| echo: false

#| message: false

#| warning: false

#Save some deets so we do not need to save large objects

ssu_default_num_asvs <- dim(ssu18_ps_work_perfect_sim$filtX)[2]

ssu_ord_num_asvs <- dim(ssu18_ps_work_ord_perfect_sim$filtX)[2]

its_default_num_asvs <- dim(its18_ps_work_perfect_sim$filtX)[2]

its_ord_num_asvs <- dim(its18_ps_work_ord_perfect_sim$filtX)[2]

ssu_default_pvals <- ssu18_ps_work_perfect_sim$pvals

ssu_ord_pvals <- ssu18_ps_work_ord_perfect_sim$pvals

its_default_pvals <- its18_ps_work_perfect_sim$pvals

its_ord_pvals <- its18_ps_work_ord_perfect_sim$pvals

```

How many ASVs were retained after filtering?

First the 16S rRNA data set. Default ordering resulted in

**`r ssu_default_num_asvs`** ASVs

and reordering the data resulted in

**`r ssu_ord_num_asvs`** ASVs.

And then the ITS data set. Default ordering resulted in

**`r its_default_num_asvs`** ASVs

and reordering the data resulted in

**`r its_ord_num_asvs`** ASVs.

For some reason, the package does not remove based on the p value cutoff that we set earlier (`0.05`). So we need to filter out the ASVs that have a higher p-value than the cutoff.

```{r echo=FALSE, results='hold', eval=TRUE, render=pander::pander}

cat("Total 16S rRNA ASVs with p-value less than", ssu_per_pval[1], "\n")

tmp_df <- ssu_default_pvals

tmp_df <- data.frame(tmp_df)

pval_asv <- tmp_df %>% dplyr::summarise(count = sum(tmp_df <= ssu_per_pval))

print(paste("default order: ASVs before checking p value was", ssu_default_num_asvs, "and after was", pval_asv$count[1]), quote = FALSE)

tmp_df <- ssu_ord_pvals

tmp_df <- data.frame(tmp_df)

pval_asv_ord <- tmp_df %>% dplyr::summarise(count = sum(tmp_df <= ssu_per_pval))

print(paste("decreasing order: ASVs before checking p value was", ssu_ord_num_asvs, "and after was", pval_asv_ord$count[1]), quote = FALSE)

print("--------------------------------------", quote = FALSE)

cat("Total ITS ASVs with p-value less than", its_per_pval[1], "\n")

tmp_df <- its_default_pvals

tmp_df <- data.frame(tmp_df)

pval_asv <- tmp_df %>% dplyr::summarise(count = sum(tmp_df <= its_per_pval))

print(paste("default order: ASVs before checking p-value was", its_default_num_asvs, "and after was", pval_asv$count[1]), quote = FALSE)

tmp_df <- its_ord_pvals

tmp_df <- data.frame(tmp_df)

pval_asv_ord <- tmp_df %>% dplyr::summarise(count = sum(tmp_df <= its_per_pval))

print(paste("decreasing order: ASVs before checking p-value was", its_ord_num_asvs, "and after was", pval_asv_ord$count[1]), quote = FALSE)

```

Now we can make phyloseq objects. Manual inspection of the results from `PERFect_sim` for the 16S rRNA data indicated that using the **decreasing** order and filtering p-values less than `r ssu_per_pval` yielded in the best results. Manual inspection of the results from `PERFect_sim` for the ITS data indicated that using the **default** order and filtering p-values less than `r its_per_pval` yielded in the best results. These approaches limited the number of ASVs found in only 1 or 2 samples. So first we filter out ASVs with p-values lower than the defined cutoff and then make the objects.

```{r}

# pvalue cutoffs set earlier