In sponslerdb/vsearchR:

Load vsearch output

Philly 2018 blanks

trnL_PH_blanks_2018 <- load_blast6("~/vsearchr/inst/extdata/trnL_vsearch_philly_2018.output/blanks/")
ITS1_PH_blanks_2018 <- load_blast6("~/vsearchr/inst/extdata/ITS1_vsearch_philly_2018.output/blanks/")
ITS2_PH_blanks_2018 <- load_blast6("~/vsearchr/inst/extdata/ITS2_vsearch_philly_2018.output/blanks/")

Load taxonomy tables

Philly 2018

trnL_PH_tax <- load_MTXA("~/vsearchr/inst/extdata/trnL_PH_Amplicons2.tax")
ITS1_PH_tax <- load_MTXA("~/vsearchr/inst/extdata/ITS1_PH_Amplicons.tax")
ITS2_PH_tax <- load_MTXA("~/vsearchr/inst/extdata/ITS2_PH_Amplicons.tax")

Join vsearch output to taxonomy tables

Philly 2018 blanks

trnL_PH_blanks_2018_join <- tax_join(trnL_PH_blanks_2018, trnL_PH_tax, min_id = 97.0, min_len = 150)
ITS1_PH_blanks_2018_join <- tax_join(ITS1_PH_blanks_2018, ITS1_PH_tax, min_id = 95.0, min_len = 300)
ITS2_PH_blanks_2018_join <- tax_join(ITS2_PH_blanks_2018, ITS2_PH_tax, min_id = 95.0, min_len = 300)

Tally genera by sample

Philly 2018 blanks

trnL_PH_blanks_2018_tally <- tally_gen(trnL_PH_blanks_2018_join) %>%
  select(sample, genus, read_count) 

ITS1_PH_blanks_2018_tally <- tally_gen(ITS1_PH_blanks_2018_join) %>%
 select(sample, genus, read_count) 

ITS2_PH_blanks_2018_tally <- tally_gen(ITS2_PH_blanks_2018_join) %>%
 select(sample, genus, read_count) 


trnL_PH_blanks_2018_read_summary <- trnL_PH_blanks_2018_tally %>%
  group_by(sample) %>%
  summarize(total_reads = sum(read_count)) %>%
  mutate(class = case_when(
    grepl("library", sample) ~ "lib",
    grepl("template", sample) ~ "temp"
  ))

ITS1_PH_blanks_2018_read_summary <- ITS1_PH_blanks_2018_tally %>%
  group_by(sample) %>%
  summarize(total_reads = sum(read_count)) %>%
  mutate(class = case_when(
    grepl("library", sample) ~ "lib",
    grepl("template", sample) ~ "temp"
  ))

ITS2_PH_blanks_2018_read_summary <- ITS2_PH_blanks_2018_tally %>%
  group_by(sample) %>%
  summarize(total_reads = sum(read_count)) %>%
  mutate(class = case_when(
    grepl("library", sample) ~ "lib",
    grepl("template", sample) ~ "temp"
  ))

write_csv(trnL_PH_blanks_2018_tally, "~/vsearchr/inst/extdata/output_2018/trnL_control_tally_2018.csv")
write_csv(ITS1_PH_blanks_2018_tally, "~/vsearchr/inst/extdata/output_2018/ITS1_control_tally_2018.csv")
write_csv(ITS2_PH_blanks_2018_tally, "~/vsearchr/inst/extdata/output_2018/ITS2_control_tally_2018.csv")

write_csv(trnL_PH_blanks_2018_read_summary, "~/vsearchr/inst/extdata/output_2018/trnL_control_read_summary_2018.csv")
write_csv(ITS1_PH_blanks_2018_read_summary, "~/vsearchr/inst/extdata/output_2018/ITS1_control_read_summary_2018.csv")
write_csv(ITS2_PH_blanks_2018_read_summary, "~/vsearchr/inst/extdata/output_2018/ITS2_control_read_summary_2018.csv")

Plot

```r ggplot(trnL_PH_blanks_2018_read_summary, aes(x = sample, y = sqrt(total_reads))) + geom_bar(stat = "identity") + theme(axis.text.x = element_text(angle = 90)) + ggtitle("trnL")

ggplot(ITS1_PH_blanks_2018_read_summary, aes(x = sample, y = sqrt(total_reads))) + geom_bar(stat = "identity") + theme(axis.text.x = element_text(angle = 90)) + ggtitle("ITS1")

ggplot(ITS2_PH_blanks_2018_read_summary, aes(x = sample, y = sqrt(total_reads))) + geom_bar(stat = "identity") + theme(axis.text.x = element_text(angle = 90)) + ggtitle("ITS2")

ggplot(ITS1_PH_blanks_2018_read_summary, aes(x = class, y = total_reads)) + geom_boxplot()

ggplot(ITS2_PH_blanks_2018_read_summary, aes(x = class, y = total_reads)) + geom_boxplot()

ggplot(trnL_PH_blanks_2018_read_summary, aes(x = class, y = total_reads)) + geom_boxplot()

summary(aov(total_reads ~ class, data = ITS1_PH_blanks_2018_read_summary)) summary(aov(total_reads ~ class, data = ITS2_PH_blanks_2018_read_summary)) summary(aov(total_reads ~ class, data = trnL_PH_blanks_2018_read_summary))

`````

sponslerdb/vsearchR documentation built on June 19, 2020, 10:01 a.m.