script:
    "../scripts/test_carnival.R"

script:
    "../scripts/test_carnival.R"

script:
    "../scripts/RMD_scripts/carnival_downstream.Rmd"

script:
    "../scripts/RMD_scripts/carnival_join.Rmd"

script:
    "../scripts/transcriptutorial/sample_resolution_dorothea.R"

script:
    "../scripts/transcriptutorial/sample_resolution_progeny.R"

script:
    "../scripts/transcriptutorial/sample_resolution_carnival.R"

script:
    "../scripts/transcriptutorial/06_analysis_CARNIVAL_results.Rmd"

script:
    "../scripts/plot-pca.R"

script:
    "../scripts/deseq2.R"

script:
    "../scripts/rlog_transform.R"

script:
    "../scripts/DESeq2_analysis.Rmd"

script:
    "../scripts/RMD_scripts/gsea_report.Rmd"

script:
    "../scripts/DEseq2_cohort.Rmd"

script:
    "../scripts/run_mitch.R"

script:
    "../scripts/export_diffexp_xlsx.R"

script:
    "../scripts/dorothea.Rmd"

script:
    "../scripts/RMD_scripts/dorothea_diffexp.Rmd"

script:
    "../scripts/RMD_scripts/progeny_analysis.Rmd"

script:
    "../scripts/count-matrix.py"

script:
    "../scripts/gtf2bed.py"

shell:
    "junction_annotation.py {params.extra} -i {input.bam} -r {input.bed} -o {params.prefix} "
    "> {log[0]} 2>&1"

shell:
    "junction_saturation.py {params.extra} -i {input.bam} -r {input.bed} -o {params.prefix} "
    "> {log} 2>&1"

shell:
    "bam_stat.py -i {input} > {output} 2> {log}"

shell:
    "infer_experiment.py -r {input.bed} -i {input.bam} > {output} 2> {log}"

shell:
    "inner_distance.py -r {input.bed} -i {input.bam} -o {params.prefix} > {log} 2>&1"

shell:
    "read_distribution.py -r {input.bed} -i {input.bam} > {output} 2> {log}"

shell:
    "read_duplication.py -i {input} -o {params.prefix} > {log} 2>&1"

shell:
    "read_GC.py -i {input} -o {params.prefix} > {log} 2>&1"

script:
    "../scripts/export_fpkm.py"

body .main-container {
  max-width: 1800px !important;
  width: 1800px !important;
}
body {
  max-width: 1800px !important;
}

if (!is.null(conf$diffexp$custom_model[[cond_id]])) {
  model_string <- conf$diffexp$custom_model[[cond_id]]
  # Rerun deseq since it was intialized
  for (x in names(conf[["diffexp"]][["contrasts"]])) {
    colData(dds_obj)[, x] <- as.factor(colData(dds_obj)[, x])
  }
  colData(dds_obj)[, cond_id] <- as.factor(colData(dds_obj)[, cond_id])
  design(dds_obj) <- as.formula(model_string)
  dds_obj <- DESeq(dds_obj)
  dds_obj@colData[, cond_id] <- relevel(dds_obj@colData[, cond_id],
    ref = contrast_groups[2]
  )
}
if (lfc_shrink) {
  dds_obj@colData[, cond_id] <- relevel(dds_obj@colData[, cond_id],
    ref = contrast_groups[2]
  )
  design(dds_obj) <- as.formula(model_string)
  dds_obj <- DESeq(dds_obj)
  write(coef_string, file = stderr())
  res <- lfcShrink(dds_obj,
    coef = stringr::str_replace_all(coef_string, "-", "."),
    type = "apeglm"
  )
} else {
  res <- results(dds_obj, contrast = c(cond_id, contrast_groups))
}
DESeq2::plotMA(res)

sample_overview <- readr::read_tsv(samp_map)

deseq_tb <- readr::read_tsv(diffexp_tb_path,
  col_names = c(
    "gene_id",
    "baseMean",
    "logFoldChange",
    "lfcSE", "stat",
    "pvalue", "padj"
  ),
  skip = 1
)

sample_overview <- sample_overview %>%
  dplyr::filter(sample %in% rownames(colData(dds_obj)))

filer <- fpkm %>% dplyr::filter(gene %in% deseq_tb$gene_id)

joined_df <- join_tables(deseq_tb, filer)

joined_df <- joined_df %>%
  dplyr::mutate(overexpressed_in = ifelse(logFoldChange > 0,
    contrast_groups[1], contrast_groups[2]
  ))

mean_tb <- joined_df %>%
  dplyr::filter(padj < pvalue_threshold &
    abs(logFoldChange) > lfc_threshold) %>%
  RNAscripts::mean_tibble_from_mat(., "logFoldChange",
    contrast_groups = contrast_groups,
    s_map = sample_overview, cond_id = cond_id
  )
DT::datatable(mean_tb)

diff_exp_trans <- SummarizedExperiment::assay(rld)[joined_df$gene, ]
rownames(diff_exp_trans) <- joined_df$gname
write(cond_id, file = stderr())
df_col <- data.frame(SummarizedExperiment::colData(rld)[, c(cond_id)])

if (!is.null(group_colors)) {
  col_annotation <- HeatmapAnnotation(
    condition = dds_obj@colData[, cond_id],
    col = list(condition = group_colors)
  )
} else {
  col_annotation <- HeatmapAnnotation(condition = dds_obj@colData[, cond_id])
}
rownames(df_col) <- colnames(SummarizedExperiment::assay(rld))
colnames(df_col) <- cond_id
index_vec <- which(joined_df$padj < pvalue_threshold &
  abs(joined_df$logFoldChange) > lfc_threshold)
diff_exp_genes <- diff_exp_trans[index_vec, ]
small_joined_df <- joined_df[joined_df$padj < pvalue_threshold &
  abs(joined_df$logFoldChange) > lfc_threshold, ]
small_joined_df <- na.omit(small_joined_df)

if (nrow(small_joined_df) > 1) {
  scaled_diffexp <- diff_exp_genes
  diffexp_heatmap <- Heatmap(
    head(scaled_diffexp[order(
      abs(small_joined_df$logFoldChange * -log10(small_joined_df$padj)),
      decreasing = TRUE
    ), ], 50),
    top_annotation = col_annotation
  )
  save_cheatmap_svg(
    x = diffexp_heatmap,
    filename = file.path(
      plot_path,
      "diffexp_heatmap.svg"
    )
  )
  diffexp_heatmap
}

#' Uniquify a vector
#'
#' @param x Input factor to check
#' @return corrected facotr value
#' @examples
#' NULL
uniquify <- function(x) {
  if (length(x) == 1) {
    x
  } else {
    sprintf("%s_%02d", x, seq_along(x))
  }
}

#' Disambiguate a vector
#'
#' Add numbered suffixes to redundant entries in a vector
#'
#' @param in_vector Vector to be disambiguated
#' @importFrom stats ave
#' @return The disambiguated vector.
#' @export
#'
#' @examples
#' NULL
#'
disambiguate_vector <- function(in_vector) {
  ave(in_vector, in_vector, FUN = uniquify)
}
diff_exp_heatmaps <- list()
diff_norm_heatmaps <- list()


if (!is.null(group_colors)) {
  col_annotation <-
    HeatmapAnnotation(
      condition = dds_obj@colData[dds_obj@colData[, cond_id] %in%
        contrast_groups, cond_id],
      col = list(condition = group_colors)
    )
} else {
  col_annotation <- HeatmapAnnotation(
    condition = dds_obj@colData[dds_obj@colData[, cond_id] %in%
      contrast_groups, cond_id]
  )
}

filtered_join_df <-
  joined_df[] %>% dplyr::filter(padj < pvalue_threshold &
    abs(logFoldChange) > lfc_threshold)
# Since it can happen that duplicate gene names appear -> uplicates are marked

filtered_join_df$gname <- disambiguate_vector(filtered_join_df$gname)

s_table <-
  filter_split_table(rld, contrast_groups, filtered_join_df,
    reorder_genes = FALSE,
    cond_id = cond_id
  )
r_means <- rowMeans(s_table)
r_sds <- apply(s_table, 1, sd)
s_table_centered <- s_table - r_means
s_table_norm <- t(scale(t(s_table)))



diff_exp_heatmaps <- ComplexHeatmap::Heatmap(
  s_table_centered,
  cluster_rows = TRUE,
  top_annotation = col_annotation,
  name = "Centered exp.",
  row_split = filtered_join_df$overexpressed_in,
  column_names_rot = 45,
  show_row_names = nrow(s_table_norm) < 100,
)
diff_norm_heatmaps <- ComplexHeatmap::Heatmap(
  s_table_norm,
  cluster_rows = TRUE,
  top_annotation = col_annotation,
  name = "z-score exp.",
  column_names_rot = 45,
  show_row_names = nrow(s_table_norm) < 100
)

if (nrow(filtered_join_df) > 0) {
  save_cheatmap_svg(diff_exp_heatmaps,
    file.path(plot_path, "split_heatmap.svg"),
    width = 14, height = 12
  )

  save_cheatmap_svg(diff_norm_heatmaps,
    file.path("standard_redblue_heatmap.svg"),
    width = 14, height = 12
  )

  plot(diff_exp_heatmaps)
  plot(diff_norm_heatmaps)
}

genereate_gsea_plots <- function(gsea_obj, gsea_name) {
  cat("\n")
  cat("### ", gsea_name , " \n")
  cat("\n")


}
msig_enrichment <- enrich_list$msig_enrichment

msig_gsea <- enrich_list$msig_gsea

kegg <- enrich_list$kegg

reactome_stuff <- enrich_list$reactome_stuff

ensembl_gene_list <- joined_df %>% dplyr::select(c(gene, stat))
gene_list <- joined_df %>% dplyr::select(c(gname, stat))

better_dotplot <- function(gset, c_groups = contrast_groups) {
  pos_gsea <- gset %>% dplyr::filter(NES > 0.5) %>% dplyr::arrange(desc(NES))

  if (nrow(pos_gsea) > 1) {
    dp_pos_nes <-
      dotplot(
        pos_gsea,
        size = "NES",
        color = "p.adjust",
        showCategory = 20,
        title = glue::glue("Pathways enriched in {contrast_groups[1]}")
      ) +
      scale_size(range = c(1, 7), limits = c(1, max(gset@result$NES)))
  } else {
    dp_pos_nes <- NULL
  }

  neg_gsea <- gset %>% dplyr::filter(NES < -0.5) %>% dplyr::arrange(NES)

  if (nrow(neg_gsea) > 1) {
    dp_neg_nes <-
      dotplot(
        neg_gsea,
        size = "NES",
        color = "p.adjust",
        showCategory = 20,
        title = glue::glue("Pathways enriched in {contrast_groups[2]}")
      ) +
      scale_size(range = c(7, 1), limits = c(min(gset@result$NES), -1))
  } else {
    dp_neg_nes <- NULL
  }

  return(list(dp_pos_nes, dp_neg_nes))
}
#' Title
#'
#' @param d_plot
#' @param p_group
#' @param p_path
#'
#' @return
#' @export
#'
#' @examples
save_dotplots <- function(d_plot, p_group, p_path = plot_path, gsea_type) {
  ggplot2::ggsave(
    filename = file.path(
      p_path,
      glue::glue("{p_group}_{gsea_type}_dplot.svg")
    ),
    d_plot, width = 10, height = 10
  )
}

org_db <- get_org_db(organism)

do_gseaplots <- function(enrich_res, gsea_name) {
  cat("\n")
  cat("### ", gsea_name , " \n")
  cat("\n")

  id_class <- conf[["gsea"]][[gsea_name]][["id_class"]]
  all_plots <- list()
  if (!is.null(enrich_res)) {
    if (nrow(enrich_res) > 1) {
      enrich_res <- clusterProfiler::setReadable(enrich_res, org_db, id_class)
      z <- enrichplot::upsetplot(enrich_res)
      all_plots$upset <- z
      if (conf[["gsea"]][[gsea_name]][["use_gsea"]]) {
        enrich_plots <- RNAscripts::plot_enrichment(enrich_res, 
          X = "ID", Y = "NES",
          pval = "p.adjust"
        )
        ggsave(filename = file.path(plot_path, paste0(
          contrast_groups[1], "_vs_",
          contrast_groups[2],
          "_", gsea_name, "_enrichplot.svg"
        )), plot = enrich_plots)
        print(enrich_plots)
        all_plots$enrichplots <- enrich_plots

        dplot <- better_dotplot(enrich_res, c_groups = contrast_groups)
        x<- purrr::map(dplot, print)
        purrr::map2(dplot, contrast_groups, save_dotplots, gsea_type = gsea_name )
        all_plots$dotplots <- dplot      
        cnetplot <- cnetplot(enrich_res,
                                foldChange = set_names(joined_df$logFoldChange, gene_list$gname),
                                node_label = "all",
                                layout = "nicely",
                                cex_label_category = 0.8,
                                cex_label_gene = 0.6
                                )
        print(cnetplot)

        ggsave(
          filename = file.path(plot_path, paste0(
            contrast_groups[1],
            "_vs_", contrast_groups[2], "_", gsea_name,
            "_cnet.svg"
          )),
          cnetplot, width = 16, height = 16
        )
        all_plots <- cnetplot
      } else {
        enrich_plots <- NULL
        all_plots$dotplots <- dotplot(msig_enrichment) 
        print(dotplot)
      }
    } else {
      all_plots <- NULL
    }
  }
  cat("\n \n")
  return(all_plots)
}


all_plots <- purrr::map2(enrich_list, names(enrich_list), do_gseaplots)

## READ all DIFFexp tables
diff_exp_tables <- purrr::map(diffexp_tables_paths,
  readr::read_tsv,
  col_names = c(
    "gene_id",
    "baseMean",
    "logFoldChange",
    "lfcSE",
    "stat",
    "pvalue",
    "padj"
  ), skip = 1
)
names(diff_exp_tables) <- contrast_list
all_ensembl <- diff_exp_tables[[1]][, 1]
### Set Rownames for data frame
diff_exp_frames <- purrr::map(diff_exp_tables, function(x) {
  gene_ids <- x$gene_id
  diff_exp_fr <- as.data.frame(x[, -1])
  rownames(diff_exp_fr) <- gene_ids
  diff_exp_fr
})
### Read FPKM Table
fpkm <- readr::read_tsv(fpkm_path)
all_genes <- fpkm %>%
  dplyr::filter(gene %in% all_ensembl$gene_id) %>%
  dplyr::select(gname)

Big_tables <- purrr::map(diff_exp_tables, join_tables, fpkm = fpkm)

#' Create dynamic table from significant gene names
#'
#' @param sig_gene_names import list of genes to put in matrix
#'
#' @return DT::datatable
#' @export
#'
#' @examples
create_overlap_matrix <- function(sig_gene_names) {
  gene_name_index <- unlist(sig_gene_names) %>%
    unique() %>%
    sort()

  #' Takes two character vecotrs (info and index) and checks which index fields are present in info
  #'
  #' @return returns a vecor of length(index) with 1 or 0 depending on True or False
  check_against_index <- function(info, index) {
    return(as.numeric(index %in% info))
  }

  gene_matrix <- sapply(sig_gene_names, check_against_index, index = gene_name_index)
  rownames(gene_matrix) <- gene_name_index


  gene_matrix <- cbind(gene_matrix, total = rowSums(gene_matrix))

  ordered_matrix <- gene_matrix[order(gene_matrix[, "total"], decreasing = T), ]

  DT::datatable(ordered_matrix)
}

ifelse(all(purrr::map(
  sig_gene_names,
  is.null
)),
create_overlap_matrix(sig_gene_names),
print("No enrichments present")
)

plot_gene_sets <- function(gset_res, plot_title) {
  print(plot_title)
  if (!is.null(gset_res)) {
    if (any(gset_res@result$p.adjust < 0.05)) {
      cnet_plot <- enrichplot::cnetplot(gset_res) + ggplot2::ggtitle(plot_title)
      return(cnet_plot)
    }
  } else {
    return(NULL)
  }
}
if (exists("H_res")) {
  cnet_plots <- purrr::pmap(list(H_res, names(H_res)), plot_gene_sets)
}

H_barplots <- purrr::map(H_res, barplot)
H_heatmaps <- purrr::map(H_res, clusterProfiler::heatplot, showCategory = 10)


for (set_n in names(H_res)) {
  if (any(H_res[[set_n]]@result$p.adjust < 0.05)) {
    cat("### ", set_n, "H set \n")
    y <- H_barplots[[set_n]]
    ifelse(nrow(y$data) > 0, plot(y), "")
    ifelse(nrow(H_res[[set_n]]) > 0,
      plot(dotplot(H_res[[set_n]], orderBy = "x")),
      ""
    )
    cat("\n")
    cat("\n")
    ifelse(nrow(cnet_plots[[set_n]]$data) > 5,
      plot(cnet_plots[[set_n]]), ""
    )
    cat("\n \n")
  }
}

dds_obj <- readRDS(dds_path)

rld <- readRDS(rld_path)

sig_transcripts <- purrr::map(sig_gene_map, function(x) {
  dplyr::pull(x, gene)
}) %>%
  unlist() %>%
  unique()

# select <- order(var(counts(dds_obj,normalized=TRUE)),
#                decreasing=TRUE)[1:400]
df <- data.frame(SummarizedExperiment::colData(dds_obj)[, c(cond_id)])
rownames(df) <- colnames(SummarizedExperiment::assay(rld))
colnames(df) <- cond_id

annotation_cols <- list(cond_id = group_colors)

ha <- ComplexHeatmap::HeatmapAnnotation(condition = rld@colData[, cond_id], col = annotation_cols)
normalized_expression <- SummarizedExperiment::assay(rld[sig_transcripts, ]) - rowMeans(SummarizedExperiment::assay(rld[sig_transcripts, ]))
ComplexHeatmap::Heatmap(normalized_expression,
  cluster_rows = TRUE, show_row_names = FALSE,
  # col = viridis::viridis(100),
  top_annotation = ha,
  name = "Normalized gene expression in differentially expressed genes",
  #  km = 4
)

assay_data <- SummarizedExperiment::assay(rld)
metadata <- SummarizedExperiment::colData(rld)

pca_data <- PCAtools::pca(assay_data, metadata = metadata, removeVar = 0.1)
scree <- screeplot(pca_data)
scree

PCAtools::biplot(pca_data, colby = cond_id, legendPosition = "right")
plot <- PCAtools::pairsplot(pca_data,
  components = PCAtools::getComponents(pca_data, c(1:5)),
  triangle = F,
  hline = 0, vline = 0,
  pointSize = 0.8,
  gridlines.major = FALSE, gridlines.minor = FALSE,
  colby = cond_id,
  title = "Pairs plot", plotaxes = T,
  margingaps = ggplot2::unit(c(-0.01, -0.01, -0.01, -0.01), "cm"),
  legendPosition = "none"
)
plot