diff --git a/code/00-functions.r b/code/00-functions.r
index 8662799..2898c0f 100644
--- a/code/00-functions.r
+++ b/code/00-functions.r
@@ -51,9 +51,11 @@ denoise <- function(marker){
     err.fwd <- learnErrors(marker.fwd, multithread = threads, randomize = T)
     err.fwd.plot <- plotErrors(err.fwd, nominalQ = T) + ggtitle(paste(insert, 'forward read error model'))
     file.path(logs, paste0(marker, '-error-fwd.png')) |> ggsave(err.fwd.plot, width = 12, height = 9)
+    file.path(logs, paste0(marker, '-error-fwd.rds')) |> saveRDS(err.fwd.plot, file = _)
     err.rev <- learnErrors(marker.rev, multithread = threads, randomize = T)
     err.rev.plot <- plotErrors(err.rev, nominalQ = T) + ggtitle(paste(insert, 'reverse read error model'))
     file.path(logs, paste0(marker, '-error-rev.png')) |> ggsave(err.rev.plot, width = 12, height = 9)
+    file.path(logs, paste0(marker, '-error-rev.rds')) |> saveRDS(err.rev.plot, file = _)
     
     # Denoise reads in both directions ####
     dada.fwd <- dada(derep.fwd, err = err.fwd, multithread = threads, pool = 'pseudo')
@@ -112,11 +114,10 @@ denoise <- function(marker){
     in.start <- which(shift$pair == start) + 1
     in.sum <- shift[in.start:nrow(shift), 'depth'] |> sum()
     
-    paste0(correct, ' / ', total, '\n', insert, '-trimmed reads\ncorrectly demultiplexed', '\n\n',
-           del.sum, ' reads with deletions in ', stop, '\n',
-           in.sum, ' reads with insertions in ', start, '\n\n') |>
-        cat(file = file.path(logs, paste0(marker, '-indel.txt')
-        ))
+    indel <- data.frame(primers = insert, total = total, correct = correct,
+                        deletions = del.sum, insertions = in.sum)
+    
+    file.path(logs, paste0(marker, '-indel.txt')) |> write.table(indel, file = _, row.names = F, quote = F)
     
     # Plot frameshift errors ####
     frameshift <- ggplot(shift, aes(x = pair, y = depth)) +
@@ -132,7 +133,8 @@ denoise <- function(marker){
               axis.title.x = element_text(face = 'bold'),
               axis.title.y = element_text(face = 'bold'))
     
-    file.path(logs, paste0(marker, '-frameshift.png')) |> ggsave(frameshift, width = 12, height = 9)
+    file.path(out, paste0(marker, '-frameshift.png')) |> ggsave(frameshift, width = 12, height = 9)
+    file.path(out, paste0(marker, '-frameshift.rds')) |> saveRDS(frameshift, file = _)
 }
 
 # 04-compile.r ####
diff --git a/code/06-analyze.r b/code/06-analyze.r
index cb55a54..4759973 100644
--- a/code/06-analyze.r
+++ b/code/06-analyze.r
@@ -1,7 +1,6 @@
 # Determine the viability of hamPCR ####
 
 # Load packages ####
-library(patchwork)
 library(ggplot2)
 library(caret)
 library(nlme)
@@ -128,6 +127,7 @@ load <- ggplot(standard, aes(x = log10(dilution), y = log10(mean_noga_load)
           legend.title = element_text(face = 'bold'))
 
 file.path(out, 'load.png') |> ggsave(load, width = 12, height = 9)
+file.path(out, 'load.rds') |> saveRDS(load, file = _)
 
 # Plot the relationship between frameshift pairs and model residuals ####
 set.seed(666)
@@ -145,6 +145,7 @@ shifts <- ggplot(standard, aes(x = fun_n, y = gi_n, color = res)) +
           legend.title = element_text(face = 'bold'))
 
 file.path(out, 'residuals.png') |> ggsave(shifts, width = 9, height = 6)
+file.path(out, 'residuals.rds') |> saveRDS(shifts, file = _)
 
 # Extract OTU tables from the entire dataset ####
 fun.tab <- subset(meta, template == 'dfssmt',
@@ -195,6 +196,7 @@ taxa <- ggplot(long, aes(x = id, y = value, fill = OTU)) +
           axis.ticks.x = element_blank())
 
 file.path(out, 'taxa.png') |> ggsave(taxa, width = 6, height = 6)
+file.path(out, 'taxa.rds') |> saveRDS(taxa, file = _)
 
 # Compare relative abundance and load hierarchical clustering assignments ####
 rownames(fun.ra) <- 1:nrow(fun.ra)