suppressPackageStartupMessages(library(tidyverse))

dir <- ifelse(!is.na(snakemake@params[['out_dir']]), snakemake@params[['out_dir']], '.')

filename_assess_assembly_all_scores_tsv <- paste0(dir, "/pomoxis/assess_assembly_all_scores.tsv")
filename_assess_assembly_all_meanQ_pdf <- "assess_assembly_all_meanQ.pdf"

# plot assess_assembly mean Q scores for each assembly and reference

df <- read_tsv(filename_assess_assembly_all_scores_tsv,
  col_names = c("assembly", "reference", "Qscore", "percErr"),
  col_types = "ffdc"
)

for (i_ref in unique(df$reference)) {
  df.plot <- df %>%
    filter(reference == i_ref)

  p <- df.plot %>%
    ggplot(aes(y = reorder(assembly, Qscore, FUN = max), x = Qscore)) +
    geom_point(shape = 21, size = 2, fill = "deepskyblue3") +
    theme_bw() +
    theme(legend.position = "bottom") +
    labs(title = "pomoxis / assess-assembly", subtitle = paste("Reference:", i_ref)) +
    ylab("") +
    xlab("Q-score ") +
    theme(axis.text = element_text(size = rel(0.75)))

  filename_out <- paste0(dir, "/pomoxis/", i_ref, "_", filename_assess_assembly_all_meanQ_pdf)
  writeLines(paste("Saving plot to", filename_out))
  ggsave(p, filename = filename_out)
}

suppressPackageStartupMessages(library(tidyverse))

dir <- ifelse(!is.na(snakemake@params[['out_dir']]), snakemake@params[['out_dir']], '.')

filename_busco_stats_tsv <- paste0(dir, "/busco/all_stats.tsv")
filename_busco_pdf <- paste0(dir, "/busco/busco_stats.pdf")

# plot busco

df <- read_tsv(filename_busco_stats_tsv,
  col_names = c("assembly", "Complete", "Fragmented", "Missing", "n"),
  col_types = "fdddd"
) %>%
  gather(key = "measure", value = "percent", -assembly, -n)

p <- df %>%
  ggplot(aes(y = reorder(assembly, percent, max), x = percent)) +
  geom_text(data = filter(df, measure == "Complete"), aes(color = measure, label = percent), size = 2.5, hjust=1.4) +
  geom_point(size = 2, aes(color = measure)) +
  theme_bw() +
  theme(legend.position = "bottom") +
  labs(title = "BUSCO") +
  xlab("Percent") +
  ylab("") +
  xlim(0, 100) +
  theme(axis.text.y = element_text(size = rel(0.75))) +
  scale_color_discrete(guide = guide_legend(title = ""))

filename_out <- filename_busco_pdf
writeLines(paste("Saving plot to", filename_out))
ggsave(p, filename = filename_out)

suppressPackageStartupMessages(library(tidyverse))

dir <- ifelse(!is.na(snakemake@params[['out_dir']]), snakemake@params[['out_dir']], '.')

# diamond for uniprot -------------------------------------------------------------------------------------
dir_ideel_diamond <- paste0(dir, "/ideel/diamond/")
filename_ideel_hist_pdf <- paste0(dir, "/ideel/ideel_uniprot_histograms.pdf")
filename_ideel_box_pdf <- paste0(dir, "/ideel/ideel_uniprot_boxplots.pdf")

filelist <- list.files(path = dir_ideel_diamond, pattern = ".*\\.tsv", recursive = FALSE, full.names = FALSE)
df_list <- vector("list", length(filelist))

for (i in seq_along(filelist)) {
  filename <- paste0(dir_ideel_diamond, filelist[[i]])
	df <- read_tsv(filename,
									 col_names = c("qseqid", "sseqid", "qlen", "slen"),
									 col_types = "ccii"
    ) %>%
    mutate(assembly = filelist[[i]]) %>%
    mutate(assembly = str_remove(assembly, "\\.tsv")) %>%
    mutate(SCov = qlen / slen)

  df_list[[i]] <- df
}
df.all <- dplyr::bind_rows(df_list) %>% mutate(assembly = factor(assembly))

p <- ggplot(df.all, aes(x = SCov)) +
  geom_histogram(aes(color=assembly, fill=assembly), alpha = 0.8, position = "identity", binwidth = 0.01, size=0.2) +
  #geom_density(aes(fill = assembly), alpha = 0.4, color = "grey40") +
  facet_wrap(~assembly) +
  theme_bw() +
  ggtitle(paste("ideel with Uniprot Sprot")) +
  ylab("") + xlab("qlen / slen") +
  xlim(0.5, 1.5) +
  scale_color_discrete(guide = FALSE) +
  scale_fill_discrete(guide = FALSE) +
  theme(strip.text.x = element_text(size = 6)) +
  theme(axis.text = element_text(size = rel(0.75)))

filename_out <- filename_ideel_hist_pdf
writeLines(paste("Saving plot to", filename_out))
ggsave(p, filename = filename_out)


p.box <- df.all %>%
  mutate(assembly = fct_reorder(assembly, SCov, mean)) %>%
  ggplot(aes(y = assembly, x = SCov)) +
  geom_boxplot(outlier.shape = NA) +
  #geom_histogram(aes(color=assembly, fill=assembly), alpha = 0.8, position = "identity", binwidth = 0.01, size=0.2) +
  #geom_density(aes(fill = assembly), alpha = 0.4, color = "grey40") +
  #facet_wrap(~assembly) +
  theme_bw() +
  ggtitle(paste("ideel with Uniprot Sprot")) +
  ylab("") + xlab("qlen / slen") +
  xlim(0.5, 1.5) +
  theme(axis.text = element_text(size = rel(0.75)))

filename_out <- filename_ideel_box_pdf
writeLines(paste("Saving plot to", filename_out))
ggsave(p.box, filename = filename_out)


# diamond for references --------------------------------------------------------------------------------
dir_ideel_diamond_ref <- paste0(dir, "/ideel/diamond-ref/")

ref_list <- list.dirs(path = dir_ideel_diamond_ref, recursive = FALSE, full.names = FALSE)

for (r in seq_along(ref_list)) {
  filelist <- list.files(path = paste0(dir_ideel_diamond_ref,"/", ref_list[[r]]), pattern = ".*\\.tsv", recursive = FALSE, full.names = FALSE)
  df_list <- vector("list", length(filelist))

  for (i in seq_along(filelist)) {
    filename <- paste0(dir_ideel_diamond_ref, "/", ref_list[[r]], "/", filelist[[i]])
    df <- read_tsv(filename, col_names = c("qseqid", "sseqid", "qlen", "slen"), col_types = "ccii") %>%
      mutate(assembly = filelist[[i]]) %>%
      mutate(assembly = str_remove(assembly, "\\.tsv")) %>%
      mutate(assembly = str_remove(assembly, paste0("_",ref_list[[r]]))) %>%
      mutate(SCov = qlen / slen)

    df_list[[i]] <- df
  }
  df.all <- dplyr::bind_rows(df_list) %>% mutate(assembly = factor(assembly))

	p <- ggplot(df.all, aes(x = SCov)) +
		geom_histogram(aes(color=assembly, fill=assembly), alpha = 0.8, position = "identity", binwidth = 0.01, size=0.2) +
		#geom_density(aes(fill = assembly), alpha = 0.4, color = "grey40") +
		facet_wrap(~assembly, scales = "free_x") +
		theme_bw() +
		labs(title = "ideel", subtitle = paste("Reference:", ref_list[[r]])) +
		ylab("") + xlab("qlen / slen") +
		xlim(0.5, 1.5) +
		scale_color_discrete(guide = FALSE) +
		scale_fill_discrete(guide = FALSE) +
		theme(strip.text.x = element_text(size = 6)) +
		theme(axis.text = element_text(size = rel(0.75)))

	filename_pdf <- paste0(dir, "/ideel/", ref_list[[r]], "_ideel_histograms.pdf")
	writeLines(paste("Saving plot to", filename_pdf))
	ggsave(p, filename = filename_pdf)

	p.box <- df.all %>%
	  mutate(assembly = fct_reorder(assembly, SCov, mean)) %>%
	  ggplot(aes(y = assembly, x = SCov)) +
	  geom_boxplot(outlier.shape = NA) +
	  theme_bw() +
	  labs(title = "ideel", subtitle = paste("Reference:", ref_list[[r]])) +
	  ylab("") + xlab("qlen / slen") +
	  xlim(0.5, 1.5) +
	  theme(axis.text.y = element_text(size = rel(0.75)))

	filename_pdf <- paste0(dir, "/ideel/", ref_list[[r]], "_ideel_boxplots.pdf")
	writeLines(paste("Saving plot to", filename_pdf))
	ggsave(p.box, filename = filename_pdf)

}

suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(tidytext))

dir <- ifelse(!is.na(snakemake@params[['out_dir']]), snakemake@params[['out_dir']], '.')

filename_nucdiff_stats_tsv <- paste0(dir, "/nucdiff/all_stats.tsv")
filename_suffix_nucdiff_stats_tsv <- "_nucdiff_stats.pdf"

# plot nucdiff stats for each reference

df <- read_tsv(filename_nucdiff_stats_tsv,
  col_names = c("assembly", "reference", "measure", "value"),
  col_types = "fffd"
) %>%
  mutate(measure = fct_relevel(measure, "Substitutions"))

for (i_ref in unique(df$reference)) {
  df.plot <- df %>%
    filter(reference == i_ref)

  p <- df.plot %>%
    filter(reference == i_ref) %>%
    group_by(measure) %>% 
      mutate(assembly = fct_reorder(assembly, -value)) %>%
    ungroup() %>%
    ggplot(aes(y = assembly, x = value)) +
    #geom_line(aes(group = reference), color = "grey70", size = 1) +
    geom_point(shape = 21, size = 2, fill = "deepskyblue3") +
    facet_wrap(~measure, scales = "free_x") +
    theme_bw() +
    theme(legend.position = "bottom") +
    labs(title = "nucdiff", subtitle = paste("Reference:", i_ref)) +
    ylab("") +
    xlab("") +
    theme(
      panel.grid.major.x = element_blank(),
      axis.text = element_text(size = rel(0.75))
    )
  filename_out <- paste0(dir, "/nucdiff/", i_ref, filename_suffix_nucdiff_stats_tsv)
  writeLines(paste("Saving plot to", filename_out))
  ggsave(p, filename = filename_out)
}

shell:
    """
    minimap2 -x asm5 -t {threads} --MD -a {input.ref} {input.assembly} 2>{log} | samtools sort -o {output} - >>{log} 2>&1
    """

shell:
    """
    quast -t {threads} --glimmer -o {out_dir}/quast/{wildcards.ref} -r {input.reference} {input.fa} >{log} 2>&1
    """

shell:
    """
    prodigal -a {output} -i {input} >{log} 2>&1
    # remove stop codon (*) from end of sequences
    sed -i 's/*$//' {output}
    """

shell:
    """
    seqkit stats -T -t protein {input} >{output} 2>{log}
    """

shell:
    """
    wget -N -P {out_dir}/bakta https://zenodo.org/record/7669534/files/db-light.tar.gz >{log} 2>&1
    tar --directory {out_dir}/bakta -xf {out_dir}/bakta/db-light.tar.gz >{log} 2>&1
    amrfinder_update --force_update --database {out_dir}/bakta/db-light/amrfinderplus-db/ >{log} 2>&1
    """