Reproducible workflow for … In this workflow, ….
Hit the Hide Code button to hide the R code (shown by default).
In order to run the workflow, you either need to first run the DADA2 Workflow for 2018 High Temp samples and then the Data Preparation workflow or begin with the output files from the Data Preparation and PIME workflows. See the Data Availability page for complete details.
In this workflow we cluster ASVs into operational taxonomic units (OTUs) using the otu function from the kmer package. Briefly, the workflow consists of XXX steps.
fasta files for each object of interest. While we’re at we also save ASV, taxonomy, and sample metadata tables for each object.otu. There are several parameters we can set to control how OTUs are picked including the threshold (value between 0 and 1 giving the OTU identity cutoff) and the method (central, centroid, farthest). The output will be a table that maps ASVs to OTUs.Here we perform OTU clustering for the 16S rRNA FULL data set and the merged set but the code will work for any number of phyloseq objects.
# LOAD ONLY REMOVE WHEN WORKFLOW FINISHED
remove(list = ls())
load("page_build/trepo/otu_wf.rdata")ps_list <- c("ssu_ps_work_merge", "ssu_ps_work")
for (i in ps_list) {
## Get ASV data
tmp_get <- get(i)
tmp_asv <- data.frame(t(otu_table(tmp_get)))
tmp_asv$total <- rowSums(tmp_asv)
tmp_asv <- tmp_asv[order(tmp_asv$total, decreasing = TRUE),]
tmp_asv$total <- NULL
tmp_asv <- tmp_asv %>% tibble::rownames_to_column("ASV_ID")
tmp_asv <- tmp_asv %>% tidyr::separate(ASV_ID, c(NA, "tmp_value"),
sep = "ASV", remove = FALSE)
tmp_asv$tmp_value <- as.numeric(as.character(tmp_asv$tmp_value))
tmp_asv <- tmp_asv[order(tmp_asv$tmp_value, decreasing = FALSE), ]
tmp_asv$tmp_value <- NULL
## Get TAX data
tmp_tax <- data.frame(tax_table(tmp_get))
tmp_tax$ASV_ID <- NULL
tmp_tax <- tmp_tax %>% tibble::rownames_to_column("ASV_ID")
tmp_tax <- tmp_tax[base::match(tmp_asv$ASV_ID, tmp_tax$ASV_ID), ]
## Fasta file
tmp_fasta_df <- tmp_tax[, c(1,8)] %>% tibble::remove_rownames()
tmp_fasta <- dataframe2fas(tmp_fasta_df)
## NAMING
tmp_asv_name <- purrr::map_chr(i, ~ paste0(., "_asv"))
assign(tmp_asv_name, tmp_asv)
tmp_tax_name <- purrr::map_chr(i, ~ paste0(., "_tax"))
assign(tmp_tax_name, tmp_tax)
tmp_fasta_name <- purrr::map_chr(i, ~ paste0(., ".fasta"))
assign(tmp_fasta_name, tmp_fasta)
tmp_path <- file.path("files/trepo/otu/tables/")
write.fasta(tmp_fasta, paste(tmp_path, tmp_fasta_name, sep = ""))
write.table(tmp_asv, paste(tmp_path, tmp_asv_name, ".txt", sep = ""),
quote = FALSE, sep = "\t", row.names = FALSE)
write.table(tmp_tax, paste(tmp_path, tmp_tax_name, ".txt", sep = ""),
quote = FALSE, sep = "\t", row.names = FALSE)
rm(list = ls(pattern = "tmp_"))
rm(list = ls(pattern = ".fasta"))
}Code for clustering OTUs is based on the kmer-vignette.Important paramets are threshold (value between 0 and 1 giving the OTU identity cutoff) and method (central, centroid, farthest)
# THIS GIVES WEIRD RESULTS WHEN RUN IN LOOP
for (i in ps_list) {
tmp_fasta_name <- purrr::map_chr(i, ~paste0(., ".fasta"))
tmp_path <- file.path("files/trepo/otu/tables/")
tmp_asvs <- read.dna(paste(tmp_path, tmp_fasta_name, sep = ""),
format = "fasta")
# CANT use if different lengths
#tmp_dna <- tmp_dna[, apply(tmp_dna, 2, function(v) !any(v == 0xf0))]
tmp_otus <- kmer::otu(tmp_asvs, k = 5, threshold = 0.97,
method = "centroid", nstart = 20)
tmp_df <- data.frame(tmp_otus)
tmp_otu_name <- purrr::map_chr(i, ~paste0(., "_cluster_results"))
assign(tmp_otu_name, tmp_df)
rm(list = ls(pattern = "tmp_"))
}
ssu_ps_work_cluster_results
ssu_ps_work_merge_cluster_resultsfor (i in ps_list) {
tmp_df <- get(purrr::map_chr(i, ~paste0(., "_cluster_results")))
## Make lookup table from results
tmp_df <- tmp_df %>% tibble::rownames_to_column("ASV_ID")
tmp_df[[2]] <- paste0("OTU", tmp_df[[2]])
tmp_df <- tmp_df %>% dplyr::rename("OTU_ID" = 2)
tmp_df$temp_rep <- tmp_df$ASV_ID
tmp_df$ASV_ID <- tmp_df$ASV_ID %>% stringr::str_replace("\\*", "")
tmp_df$temp_rep <- tmp_df$temp_rep %>% stringr::str_replace("\\*", "_rep")
tmp_df <- tmp_df %>% tidyr::separate(temp_rep, c(NA, "is_asv_rep"), sep = "_", fill = "left")
tmp_df[[3]][tmp_df[[3]] != "rep"] <- "FALSE"
tmp_df[[3]][tmp_df[[3]] == "rep"] <- "TRUE"
tmp_path <- file.path("files/trepo/otu/tables/")
tmp_lookup_name <- purrr::map_chr(i, ~ paste0(., "_asv_to_otu_map"))
assign(tmp_lookup_name, tmp_df)
write.table(tmp_df, paste(tmp_path, tmp_lookup_name, ".txt", sep = ""),
quote = FALSE, sep = "\t", row.names = FALSE)
tmp_l_asv <- length(base::unique(tmp_df[[1]]))
tmp_l_otu <- length(base::unique(tmp_df[[2]]))
print(c(tmp_l_asv, tmp_l_otu))
print(tmp_l_otu - tmp_l_asv)
rm(list = ls(pattern = "tmp_"))
}
objects()#https://riptutorial.com/data-table/example/13084/using--sd-and--sdcols
for (i in ps_list) {
tmp_asv <- get(purrr::map_chr(i, ~ paste0(., "_asv")))
tmp_asv_to_otu <- get(purrr::map_chr(i, ~ paste0(., "_asv_to_otu_map")))
tmp_asv_otu <- dplyr::left_join(tmp_asv_to_otu, tmp_asv)
tmp_asv_otu <- arrange(tmp_asv_otu, OTU_ID, desc(is_asv_rep))
tmp_asv_otu <- data.table(tmp_asv_otu)
tmp_cols <- which(sapply(tmp_asv_otu, is.numeric))
tmp_otu <- tmp_asv_otu[order(OTU_ID), lapply(.SD, sum),
by = .(OTU_ID), .SDcols = tmp_cols]
tmp_otu <- tmp_otu %>% tidyr::separate(OTU_ID, c(NA, "tmp_value"),
sep = "OTU", remove = FALSE)
tmp_otu$tmp_value <- as.numeric(as.character(tmp_otu$tmp_value))
tmp_otu <- tmp_otu[order(tmp_otu$tmp_value, decreasing = FALSE), ]
tmp_otu$tmp_value <- NULL
tmp_otu <- tmp_otu %>% tibble::column_to_rownames("OTU_ID")
tmp_otu_t <- transpose(tmp_otu)
colnames(tmp_otu_t) <- rownames(tmp_otu)
rownames(tmp_otu_t) <- colnames(tmp_otu)
print(identical(colnames(tmp_otu), rownames(tmp_otu_t)))
print(identical(colnames(tmp_otu_t), rownames(tmp_otu)))
tmp_otu_t <- as.matrix(tmp_otu_t)
tmp_otu_name <- purrr::map_chr(i, ~ paste0(., "_otu_tab"))
assign(tmp_otu_name, tmp_otu_t)
rm(list = ls(pattern = "tmp_"))
}rm(list = ls(pattern = "tmp_"))
rm(list = ls(pattern = "_otu_tax"))
for (i in ps_list) {
tmp_lookup <- get(purrr::map_chr(i, ~ paste0(., "_asv_to_otu_map")))
tmp_lookup <- tmp_lookup %>% dplyr::filter(is_asv_rep == "TRUE")
tmp_tax <- get(purrr::map_chr(i, ~ paste0(., "_tax")))
tmp_tax_merge <- dplyr::left_join(tmp_lookup, tmp_tax)
tmp_tax_merge[, c(1,3)] <- list(NULL)
tmp_tax_merge <- tmp_tax_merge %>% tidyr::separate(OTU_ID, c(NA, "tmp_value"),
sep = "OTU", remove = FALSE)
tmp_tax_merge$tmp_value <- as.numeric(as.character(tmp_tax_merge$tmp_value))
tmp_tax_merge <- tmp_tax_merge[order(tmp_tax_merge$tmp_value, decreasing = FALSE), ]
tmp_tax_merge$tmp_value <- NULL
tmp_tax_merge <- tmp_tax_merge %>% tibble::remove_rownames()
tmp_tax_merge <- tmp_tax_merge %>% tibble::column_to_rownames("OTU_ID")
tmp_tax_merge <- as.matrix(tmp_tax_merge)
tmp_tax_name <- purrr::map_chr(i, ~ paste0(., "_otu_tax"))
assign(tmp_tax_name, tmp_tax_merge)
rm(list = ls(pattern = "tmp_"))
}
objects1 <- objects()Now we can create phyloseq objects for the new OTU data. This is a little tricky because in the original ps objects, the ASVs are named in order of total abundance. However, the new OTU designations are not named in this fashion.
for (i in ps_list) {
tmp_otu <- get(purrr::map_chr(i, ~ paste0(., "_otu_tab")))
tmp_tax <- get(purrr::map_chr(i, ~ paste0(., "_otu_tax")))
tmp_samp <- data.frame(sample_data(get(i)))
tmp_samp <- tmp_samp %>% dplyr::select(1:6)
# Reorder by OTU abundance
tmp_otu <- data.frame(tmp_otu)
tmp_otu <- tmp_otu %>% select(order(-colSums(tmp_otu)))
tmp_tax <- data.frame(tmp_tax)
tmp_tax <- tmp_tax %>% tibble::rownames_to_column("TEMP")
tmp_tax$NAME <- tmp_tax$TEMP
tmp_tax <- tmp_tax %>% tibble::column_to_rownames("TEMP")
tmp_tax <- tmp_tax[match(colnames(tmp_otu), tmp_tax$NAME),]
print(identical(colnames(tmp_otu), row.names(tmp_tax)))
# Save original name
tmp_org <- data.frame(colnames(tmp_otu))
tmp_org <- tmp_org %>% dplyr::rename("OTU_ID" = 1)
tmp_org_names <- purrr::map_chr(i, ~ paste0(., "_org_names"))
assign(tmp_org_names, tmp_org)
# Rename to sequential
colnames(tmp_otu) <- c(paste0("OTU", 1:ncol(tmp_otu)))
row.names(tmp_tax) <- c(paste0("OTU", 1:nrow(tmp_tax)))
tmp_tax$NAME <- NULL
print(identical(colnames(tmp_otu), row.names(tmp_tax)))
# MAke PS object & save
tmp_tax <- as.matrix(tmp_tax)
tmp_otu <- as.matrix(tmp_otu)
tmp_ps <- phyloseq(otu_table(tmp_otu, taxa_are_rows = FALSE),
sample_data(tmp_samp), tax_table(tmp_tax))
# Adding tree reorganizes OTU order
#tmp_tree <- rtree(ntaxa(tmp_ps), rooted = TRUE,
# tip.label = taxa_names(tmp_ps))
#tmp_ps <- merge_phyloseq(tmp_ps, sample_data, tmp_tree)
tmp_ps_name <- purrr::map_chr(i, ~ paste0(., "_otu"))
assign(tmp_ps_name, tmp_ps)
print(identical(readcount(tmp_ps), readcount(get(i))))
# SAVE stuff
tmp_path <- file.path("files/trepo/otu/rdata/")
saveRDS(tmp_ps, paste(tmp_path, tmp_ps_name, ".rds", sep = ""))
tmp_path2 <- file.path("files/trepo/otu/tables/")
write.table(tmp_org, paste(tmp_path2, tmp_org_names, ".txt", sep = ""),
quote = FALSE, sep = "\t", row.names = FALSE)
rm(list = ls(pattern = "tmp_"))
}
objects2 <- objects()OK, now we need to make sure that all of this fiddling with OTU names didn’t mess anything up. So we need to compare the following to make sure they are all identical:
Specifically, we will compare the lineage and the DNA sequence from all three data frames. If they are equal, we are good to go. First the full data set.
tmp_rn <- data.frame(row.names(t(otu_table(ssu_ps_work_otu))))
tmp_bind <- cbind(ssu_ps_work_org_names, tmp_rn)
tmp_bind <- tmp_bind %>% dplyr::rename("NEW_ID" = 2)
temp_map_true <- ssu_ps_work_asv_to_otu_map[ssu_ps_work_asv_to_otu_map$is_asv_rep == 'TRUE',]
tmp_join <- dplyr::left_join(temp_map_true, tmp_bind, keep = FALSE)
ssu_ps_work_asv_to_otu_to_new_otu_map <- tmp_join
rm(list = ls(pattern = "tmp_"))
temp_map <- ssu_ps_work_asv_to_otu_to_new_otu_map
temp_map[,3] <- NULL
## ORIGINAL ASV
temp_org_tax <- data.frame(tax_table(ssu_ps_work))
temp_org_tax <- temp_org_tax[temp_org_tax$ASV_ID %in% temp_map$ASV_ID,]
## ORIGINAL OTU
temp_org_otu_tax <- data.frame(ssu_ps_work_otu_tax)
temp_org_otu_tax <- temp_org_otu_tax %>% tibble::add_column(row.names(temp_org_otu_tax)) %>%
dplyr::rename(., "OTU_ID" = 8)
#RENAMED OTU
temp_rn_otu_tax <- data.frame(tax_table(ssu_ps_work_otu))
temp_rn_otu_tax <- temp_rn_otu_tax %>% tibble::add_column(row.names(temp_rn_otu_tax)) %>%
dplyr::rename(., "OTU_ID" = 8)
ssu_results <- data.frame()
for (i in row.names(temp_map)) {
tmp_x <- temp_org_tax[temp_org_tax$ASV_ID %in% temp_map[i,1],]
tmp_y <- temp_org_otu_tax[temp_org_otu_tax$OTU_ID %in% temp_map[i,2],]
tmp_z <- temp_rn_otu_tax[temp_rn_otu_tax$OTU_ID %in% temp_map[i,3],]
check_true <- all_equal(tmp_x[, 1:7], tmp_z[, 1:7], tmp_y[, 1:7])
ssu_results <- rbind(ssu_results, check_true)
rm(list = ls(pattern = "tmp_"))
}
rm(list = ls(pattern = "temp_")) [1] TRUE.
<0 rows> (or 0-length row.names)And then the merged data.
tmp_rn <- data.frame(row.names(t(otu_table(ssu_ps_work_merge_otu))))
tmp_bind <- cbind(ssu_ps_work_merge_org_names, tmp_rn)
tmp_bind <- tmp_bind %>% dplyr::rename("NEW_ID" = 2)
temp_map_true <- ssu_ps_work_merge_asv_to_otu_map[ssu_ps_work_merge_asv_to_otu_map$is_asv_rep == 'TRUE',]
tmp_join <- dplyr::left_join(temp_map_true, tmp_bind, keep = FALSE)
ssu_ps_work_merge_asv_to_otu_to_new_otu_map <- tmp_join
rm(list = ls(pattern = "tmp_"))
temp_map <- ssu_ps_work_merge_asv_to_otu_to_new_otu_map
temp_map[,3] <- NULL
## ORIGINAL ASV
temp_org_tax <- data.frame(tax_table(ssu_ps_work_merge))
temp_org_tax <- temp_org_tax[temp_org_tax$ASV_ID %in% temp_map$ASV_ID,]
## ORIGINAL OTU
temp_org_otu_tax <- data.frame(ssu_ps_work_merge_otu_tax)
temp_org_otu_tax <- temp_org_otu_tax %>%
tibble::add_column(row.names(temp_org_otu_tax)) %>%
dplyr::rename(., "OTU_ID" = 8)
#RENAMED OTU
temp_rn_otu_tax <- data.frame(tax_table(ssu_ps_work_merge_otu))
temp_rn_otu_tax <- temp_rn_otu_tax %>%
tibble::add_column(row.names(temp_rn_otu_tax)) %>%
dplyr::rename(., "OTU_ID" = 8)
its_results <- data.frame()
for (i in row.names(temp_map)) {
tmp_x <- temp_org_tax[temp_org_tax$ASV_ID %in% temp_map[i,1],]
tmp_y <- temp_org_otu_tax[temp_org_otu_tax$OTU_ID %in% temp_map[i,2],]
tmp_z <- temp_rn_otu_tax[temp_rn_otu_tax$OTU_ID %in% temp_map[i,3],]
check_true <- all_equal(tmp_x[, 1:7], tmp_z[, 1:7], tmp_y[, 1:7])
its_results <- rbind(its_results, check_true)
rm(list = ls(pattern = "tmp_"))
}
rm(list = ls(pattern = "temp_")) [1] TRUE.
<0 rows> (or 0-length row.names)The last step is to add an OTU_ID column to the new ps objects so they are the same as their ASV counterparts.
tax_table(ssu_ps_work_otu) <- cbind(tax_table(ssu_ps_work_otu),
rownames(tax_table(ssu_ps_work_otu)))
colnames(tax_table(ssu_ps_work_otu)) <- c("Kingdom", "Phylum", "Class",
"Order", "Family", "Genus",
"OTU_SEQ", "OTU_ID")
tax_table(ssu_ps_work_merge_otu) <- cbind(tax_table(ssu_ps_work_merge_otu),
rownames(tax_table(ssu_ps_work_merge_otu)))
colnames(tax_table(ssu_ps_work_merge_otu)) <- c("Kingdom", "Phylum", "Class",
"Order", "Family", "Genus",
"OTU_SEQ", "OTU_ID") for (i in ps_list){
tmp_ps <- get(purrr::map_chr(i, ~ paste0(., "_otu")))
tmp_otu <- data.frame(t(otu_table(tmp_ps)))
tmp_otu[] <- lapply(tmp_otu, as.numeric)
tmp_otu <- as.matrix(tmp_otu)
tmp_tax <- as.matrix(data.frame(tax_table(tmp_ps)))
tmp_samples <- data.frame(sample_data(tmp_ps))
tmp_amp <- merge_phyloseq(otu_table(tmp_otu, taxa_are_rows = TRUE),
tax_table(tmp_tax),
sample_data(tmp_samples))
tmp_amp_name <- purrr::map_chr(i, ~ paste0(., "_ps_otu_amp"))
assign(tmp_amp_name, tmp_amp)
rm(list = ls(pattern = "tmp_"))
}for (i in ps_list){
tmp_ps <- get(purrr::map_chr(i, ~ paste0(., "_ps_otu_amp")))
tmp_samp <- data.frame(sample_data(tmp_ps))
tmp_asv <- data.frame(otu_table(tmp_ps))
tmp_asv <- tmp_asv %>% tibble::rownames_to_column("OTU")
tmp_tax <- data.frame(tax_table(tmp_ps))
tmp_tax <- tmp_tax %>% tibble::rownames_to_column("OTU")
tmp_tax$OTU_SEQ <- NULL
#tmp_tax$OTU_ID <- tmp_tax$OTU
colnames(tmp_tax)[colnames(tmp_tax) == "OTU_ID"] <- "Species"
tmp_asv_tax <- left_join(tmp_asv, tmp_tax, by = "OTU")
tmp_amp <- amp_load(tmp_asv_tax, metadata = tmp_samp, tree = tmp_ps@phy_tree)
tmp_amp_name <- purrr::map_chr(i, ~ paste0(., "_otu_amp"))
assign(tmp_amp_name, tmp_amp)
tmp_path <- file.path("files/trepo/otu/rdata/")
saveRDS(tmp_amp, paste(tmp_path, tmp_amp_name, ".rds", sep = ""))
rm(list = ls(pattern = "tmp_"))
}
rm(list = ls(pattern = "_ps_otu_amp"))The ssu_ps_work_merge ASV phyloseq object has 6457054 reads and 47725 ASVs.
The ssu_ps_work_merge_otu OTU phyloseq object has 6457054 reads and 21910 OTUs.
The ssu_ps_work ASV phyloseq object has 6457054 reads and 47725 ASVs.
The ssu_ps_work_otu OTU phyloseq object has 6457054 reads and 21906 OTUs. save.image("page_build/trepo/otu_wf.rdata")The source code for this page can be accessed on GitHub by clicking this link.
If you see mistakes or want to suggest changes, please create an issue on the source repository.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. Source code is available at https://github.com/tropical-repo/web/, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".