shell: """
    (bowtie2 --minins {params.min_fraglength} --maxins {params.max_fraglength} --fr --no-mixed --no-discordant --al-conc-gz fastq/aligned/{wildcards.sample}_{FACTOR}-chipseq-aligned.fastq.gz --un-conc-gz fastq/unaligned/{wildcards.sample}_{FACTOR}-chipseq-unaligned.fastq.gz -p {threads} -x {params.idx_path}/{basename} -1 {input.r1} -2 {input.r2}  | \
     samtools view -buh -q {params.minmapq} - | \
     samtools sort -T .{wildcards.sample} -@ {threads} -o {output.bam} -) &> {output.log}
    mv fastq/aligned/{wildcards.sample}_{FACTOR}-chipseq-aligned.fastq.1.gz {output.aligned_r1}
    mv fastq/aligned/{wildcards.sample}_{FACTOR}-chipseq-aligned.fastq.2.gz {output.aligned_r2}
    mv fastq/unaligned/{wildcards.sample}_{FACTOR}-chipseq-unaligned.fastq.1.gz {output.unaligned_r1}
    mv fastq/unaligned/{wildcards.sample}_{FACTOR}-chipseq-unaligned.fastq.2.gz {output.unaligned_r2}
    """

shell: """
    bedtools makewindows -g <(faidx {input.fasta} -i chromsizes) -w {wildcards.windowsize} | \
    awk 'BEGIN{{FS=OFS="\t"}}{{print $1, $2, $3, ".", 0, "."}}' | \
    LC_COLLATE=C sort -k1,1 -k2,2n > {output}
    """

shell: """
    (mkdir -p qual_ctrl/fastqc/{wildcards.fqtype}
    fastqc --adapters <(echo -e "adapter\t{params.adapter}") --nogroup --noextract -t {threads} -o qual_ctrl/fastqc/{wildcards.fqtype} {input.fastq}
    unzip -p qual_ctrl/fastqc/{wildcards.fqtype}/{params.fname}_fastqc.zip {params.fname}_fastqc/fastqc_data.txt > {output}) &> {log}
    """

import argparse
import numpy as np
import pandas as pd
import pyBigWig as pybw

#given paths to bigwig files representing replicates,
#return a dictionary where keys are chromosome names and
#values are the average coverage across replicates
def average_bigwigs(coverage_paths):
    coverage = {}
    for index, path in enumerate(coverage_paths):
        bw = pybw.open(path)
        chroms = bw.chroms()
        for chrom in chroms:
            if index==0:
                coverage[chrom] = bw.values(chrom, 0, chroms[chrom], numpy=True)
            else:
                coverage[chrom] = np.add(bw.values(chrom, 0, chroms[chrom], numpy=True), coverage[chrom])
            if index==len(coverage_paths):
                coverage[chrom] = np.divide(coverage[chrom], index)
        bw.close()
    return coverage

#given "unstranded" chromosome coordinates,
#return 0-based offset of summit position from start.
#If multiple positions have the same max signal, return the mean position
def get_summit(row, coverage):
    local_coverage = coverage[row['chrom']][row['start']:row['end']]
    if not np.any(np.isfinite(local_coverage)):
        return int(len(local_coverage) / 2)
    return int(np.mean(np.argwhere(local_coverage==np.amax(local_coverage[np.isfinite(local_coverage)]))))

def main(condition_paths,
        control_paths,
        diffexp_path,
        narrowpeak_out,
        bed_out):

    #condition and control coverage are imported separately and
    #averaged across replicates in case the number of samples
    #in each group is different
    condition_coverage = average_bigwigs(condition_paths)
    coverage = average_bigwigs(control_paths)
    for chrom in coverage:
        coverage[chrom] = np.add(coverage[chrom], condition_coverage[chrom])

    #we only need to perform operations using start and end as integers,
    #so everything else can be treated as an object to avoid reformatting
    diffexp_df = pd.read_csv(diffexp_path, sep="\t",
                             dtype={'chrom':str,
                                    'start':np.uint32,
                                    'end':np.uint32,
                                    'name':str,
                                    'score':str,
                                    'strand':str,
                                    'log2FC_enrichment':str,
                                    'lfc_SE':str,
                                    'stat':str,
                                    'log10_pval':str,
                                    'log10_padj':str,
                                    'mean_counts':str,
                                    'condition_enrichment':str,
                                    'condition_enrichment_SE':str,
                                    'control_enrichment':str,
                                    'control_enrichment_SE':str})

    if diffexp_df.shape[0] > 0:
        diffexp_df['summit'] = diffexp_df.apply(get_summit, coverage=coverage, axis=1)
        diffexp_df = diffexp_df.assign(summit_start = diffexp_df['start'] + diffexp_df['summit'])
        diffexp_df = diffexp_df.assign(summit_end = diffexp_df['summit_start'] + 1)

    #NOTE: we convert NAs (found in pvalue and score columns) to zero for narrowpeak compatibility
    diffexp_df.to_csv(narrowpeak_out,
                      sep="\t",
                      columns=(['chrom', 'start', 'end', 'name', 'score', 'strand',
                               'log2FC_enrichment', 'log10_pval', 'log10_padj', 'summit'] if
                               diffexp_df.shape[0] > 0 else []),
                      header=False,
                      index=False,
                      float_format="%.3f",
                      encoding='utf-8',
                      na_rep="0")
    diffexp_df.to_csv(bed_out,
                      sep="\t",
                      columns=(['chrom', 'summit_start', 'summit_end', 'name', 'score', 'strand'] if
                          diffexp_df.shape[0] > 0 else []),
                      header=False,
                      index=False,
                      float_format="%.3f",
                      encoding='utf-8',
                      na_rep="0")

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Add back summit information to ChIP-seq differential binding results.')
    parser.add_argument('-i', dest = 'condition_paths', type=str, nargs='+', help='BigWigs for all condition samples')
    parser.add_argument('-j', dest = 'control_paths', type=str, nargs='+', help='BigWigs for all control samples')
    parser.add_argument('-d', dest = 'diffexp_path', type=str, help='differential binding results file')
    parser.add_argument('-n', dest = 'narrowpeak_out', type=str, help='output path for narrowPeak file')
    parser.add_argument('-b', dest = 'bed_out', type=str, help='output path for BED file of summit positions')
    args = parser.parse_args()

    main(args.condition_paths,
         args.control_paths,
         args.diffexp_path,
         args.narrowpeak_out,
         args.bed_out)

library(tidyverse)
library(magrittr)
library(DESeq2)
library(gridExtra)

get_countdata = function(path, samples){
    df = read_tsv(path) %>%
        select(samples) %>%
        rowid_to_column(var="index") %>%
        column_to_rownames(var="index") %>%
        as.data.frame()
    df = df[rowSums(df)>1,]
    return(df)
}

initialize_dds = function(data_path,
                          samples,
                          conditions,
                          sample_type,
                          condition_id,
                          control_id){
    dds = DESeqDataSetFromMatrix(countData = get_countdata(data_path, samples),
                                 colData = data.frame(condition = factor(conditions,
                                                                         levels = c(control_id,
                                                                                    condition_id)),
                                                      sample_type = factor(sample_type,
                                                                          levels = c("input",
                                                                                     "ChIP")),
                                                      row.names = samples),
                                 design = ~ sample_type + condition + sample_type:condition)
    return(dds)
}

extract_normalized_counts = function(dds){
    dds %>%
        counts(normalized=TRUE) %>%
        as.data.frame() %>%
        rownames_to_column(var="index") %>%
        as_tibble() %>%
        return()
}

extract_rlog_counts = function(dds){
    dds %>%
        rlog(blind=FALSE) %>%
        assay() %>%
        as.data.frame() %>%
        rownames_to_column(var="index") %>%
        as_tibble() %>%
        return()
}

build_mean_sd_df_pre = function(dds){
     dds %>%
        normTransform() %>%
        assay() %>%
        as_tibble() %>%
        rowid_to_column(var="index") %>%
        gather(sample, signal, -index) %>%
        group_by(index) %>%
        summarise(mean = mean(signal),
                  sd = sd(signal)) %>%
        mutate(rank = min_rank(dplyr::desc(mean))) %>%
        return()
}

build_mean_sd_df_post = function(counts){
    counts %>%
        gather(sample, signal, -index) %>%
        group_by(index) %>%
        summarise(mean = mean(signal),
                  sd = sd(signal)) %>%
        mutate(rank = min_rank(dplyr::desc(mean))) %>%
        return()
}

reverselog_trans <- function(base = exp(1)) {
    trans <- function(x) -log(x, base)
    inv <- function(x) base^(-x)
    scales::trans_new(paste0("reverselog-", format(base)), trans, inv,
              scales::log_breaks(base = base),
              domain = c(1e-100, Inf))
}

mean_sd_plot = function(df, ymax, title){
    ggplot(data = df, aes(x=rank, y=sd)) +
        geom_hex(aes(fill=..count.., color=..count..), bins=100, size=0) +
        geom_smooth(color="#4292c6") +
        scale_fill_viridis_c(option="inferno", name=expression(log[10](count)), guide=FALSE) +
        scale_color_viridis_c(option="inferno", guide=FALSE) +
        scale_x_continuous(trans = reverselog_trans(10),
                           name="rank(mean enrichment)",
                           expand = c(0,0)) +
        scale_y_continuous(limits = c(NA, ymax),
                           name = "SD") +
        theme_light() +
        ggtitle(title) +
        theme(text = element_text(size=8))
}

extract_deseq_results = function(dds,
                                 annotations,
                                 alpha,
                                 lfc){
    control_enrichment = results(dds,
            contrast=c(0,1,0,0),
            tidy=TRUE) %>%
        as_tibble() %>%
        select(row,
               control_enrichment = log2FoldChange,
               control_enrichment_SE = lfcSE)
    condition_enrichment = results(dds,
            contrast=c(0,1,0,1),
            tidy=TRUE) %>%
        as_tibble() %>%
        select(row,
               condition_enrichment = log2FoldChange,
               condition_enrichment_SE = lfcSE)

    results(dds,
            alpha=alpha,
            lfcThreshold=lfc,
            altHypothesis="greaterAbs",
            tidy=TRUE) %>%
        as_tibble() %>%
        left_join(control_enrichment,
                  by="row") %>%
        left_join(condition_enrichment,
                  by="row") %>%
        left_join(annotations, ., by=c("index"="row")) %>%
        arrange(padj) %>%
        mutate(name = if_else(name==".",
                              paste0("peak_", row_number()),
                              name),
               score = as.integer(pmin(-125*log2(padj), 1000))) %>%
        mutate_at(vars(pvalue, padj), ~(-log10(.))) %>%
        mutate_if(is.double, round, 3) %>%
        select(index, chrom, start, end, name, score, strand,
               log2FC_enrichment=log2FoldChange, lfc_SE=lfcSE,
               stat, log10_pval=pvalue, log10_padj=padj, mean_counts=baseMean,
               condition_enrichment, condition_enrichment_SE,
               control_enrichment, control_enrichment_SE) %>%
        return()
}

write_counts_table = function(results_df,
                              annotations,
                              counts_df,
                              output_path){
    results_df %>%
        select(1:7) %>%
        right_join(annotations %>% select(-c(name, score)),
                   by = c("index", "chrom", "start", "end", "strand")) %>%
        left_join(counts_df, by="index") %>%
        select(-index) %>%
        write_tsv(output_path) %>%
        return()
}

plot_ma = function(df_sig = results_df_filtered_significant,
                   df_nonsig = results_df_filtered_nonsignificant,
                   xvar = mean_expr,
                   yvar = log2_enrichment,
                   lfc,
                   condition,
                   control){
    xvar = enquo(xvar)
    yvar = enquo(yvar)
    ggplot() +
        geom_hline(yintercept = 0, color="black", linetype="dashed") +
        geom_hline(yintercept = c(-lfc, lfc), color="grey70", linetype="dashed") +
        stat_bin_hex(data = df_nonsig,
                     geom="point",
                     aes(x=!!xvar, y=!!yvar, alpha=..count..),
                     binwidth = c(.01, 0.01),
                     color="black", stroke=0, size=0.7) +
        stat_bin_hex(data = df_sig,
                     geom="point",
                     aes(x=!!xvar, y=!!yvar, alpha=..count..),
                     binwidth = c(.01, 0.01),
                     color="red", stroke=0, size=0.7) +
        scale_x_log10(name="mean of normalized counts") +
        scale_alpha_continuous(range = c(0.5, 1)) +
        ylab(bquote(log[2]~frac("enrichment in" ~ .(condition),
                                "enrichment in" ~ .(control)))) +
        theme_light() +
        theme(text = element_text(size=8, color="black"),
              axis.text = element_text(color = "black"),
              axis.title.y = element_text(angle=0, hjust=1, vjust=0.5),
              legend.position = "none")
}

plot_volcano = function(df = results_df_filtered,
                        xvar = log2_enrichment,
                        yvar = log10_padj,
                        lfc,
                        alpha,
                        condition,
                        control){
    xvar = enquo(xvar)
    yvar = enquo(yvar)
    ggplot() +
        geom_vline(xintercept = 0, color="black", linetype="dashed") +
        geom_vline(xintercept = c(-lfc, lfc), color="grey70", linetype="dashed") +
        stat_bin_hex(data = df,
                     geom = "point",
                     aes(x = !!xvar, y = !!yvar, color=log10(..count..)),
                     binwidth = c(0.01, 0.1),
                     alpha=0.8, stroke=0, size=0.7) +
        geom_hline(yintercept = -log10(alpha), color="red", linetype="dashed") +
        xlab(bquote(log[2] ~ frac("enrichment in" ~ .(condition),
                                  "enrichment in" ~ .(control)))) +
        ylab(expression(-log[10] ~ FDR)) +
        scale_color_viridis_c(option="inferno") +
        theme_light() +
        theme(text = element_text(size=8),
              axis.title.y = element_text(angle=0, hjust=1, vjust=0.5),
              legend.position = "none")
}

main = function(exp_table="depleted-v-non-depleted_allsamples-experimental-Rpb1-chipseq-counts-verified-coding-genes.tsv.gz",
                spike_table="depleted-v-non-depleted_allsamples-spikein-Rpb1-chipseq-counts-peaks.tsv.gz",
                samples=read_tsv(exp_table) %>% select(-c(1:6)) %>% names(),
                conditions=rep(c(rep("non-depleted",4), rep("depleted",4)), 2),
                sample_type=c(rep("input",8), rep("ChIP", 8)),
                # batches = rep(c(rep(1,2), rep(2,2)), 4),
                norm="spikenorm",
                condition="depleted",
                control="non-depleted",
                alpha=0.1,
                lfc=0,
                counts_norm_out="counts_norm.tsv",
                counts_rlog_out="counts_rlog.tsv",
                results_all_out="results_all.tsv",
                results_up_out="results_up.tsv",
                results_down_out="results_down.tsv",
                results_unchanged_out="results_unch.tsv",
                # bed_all_out="all.bed",
                # bed_up_out="up.bed",
                # bed_down_out="down.bed",
                # bed_unchanged_out="nonsignificant.bed",
                qc_plots_out="qcplots.png"){

    annotations = read_tsv(exp_table) %>%
        select(1:6) %>%
        rownames_to_column(var="index") %>%
        mutate(chrom = str_replace(chrom, "-minus$|-plus$", ""))

    dds = initialize_dds(data_path = exp_table,
                         samples = samples,
                         conditions = conditions,
                         sample_type = sample_type,
                         condition_id = condition,
                         control_id = control)

    if (norm=="spikenorm"){
        dds_spike = initialize_dds(data_path = spike_table,
                                   samples = samples,
                                   conditions = conditions,
                                   sample_type = sample_type,
                                   condition_id = condition,
                                   control_id = control) %>%
            estimateSizeFactors()
        sizeFactors(dds) = sizeFactors(dds_spike)
    } else {
        dds %<>% estimateSizeFactors()
    }
    dds %<>% estimateDispersions() %>% nbinomWaldTest()

    #extract normalized counts and write to file
    counts_norm = extract_normalized_counts(dds = dds)
    counts_rlog = extract_rlog_counts(dds = dds)

    mean_sd_df_pre = build_mean_sd_df_pre(dds)
    mean_sd_df_post = build_mean_sd_df_post(counts_rlog)

    sd_max = max(c(mean_sd_df_pre[["sd"]],
                   mean_sd_df_post[["sd"]]),
                 na.rm=TRUE)*1.01

    mean_sd_plot_pre = mean_sd_plot(df = mean_sd_df_pre,
                                    ymax = sd_max,
                                    title = expression(log[2] ~ "counts," ~ "pre-shrinkage"))
    mean_sd_plot_post = mean_sd_plot(df = mean_sd_df_post,
                                     ymax = sd_max,
                                     title = expression(regularized ~ log[2] ~ "counts"))

    results_df = extract_deseq_results(dds = dds,
                                       annotations = annotations,
                                       alpha = alpha,
                                       lfc = lfc) %>%
        mutate(chrom = str_replace(chrom, "-minus$|-plus$", ""))

    write_counts_table(results_df = results_df,
                       annotations = annotations,
                       counts_df = counts_norm,
                       output_path = counts_norm_out)
    write_counts_table(results_df = results_df,
                       annotations = annotations,
                       counts_df = counts_rlog,
                       output_path = counts_rlog_out)

    results_df %<>%
        select(-index) %>%
        write_tsv(results_all_out)
    # results_df %>%
    #     select(1:6) %>%
    #     write_tsv(bed_all_out, col_names=FALSE)

    results_df_significant = results_df %>%
        filter(log10_padj > -log10(alpha))
    results_df_nonsignificant = results_df %>%
        filter(log10_padj <= -log10(alpha)) %>%
        write_tsv(results_unchanged_out)
    # results_df_nonsignificant %>%
    #     select(1:6) %>%
    #     write_tsv(bed_unchanged_out, col_names=FALSE)

    results_df_significant %>%
        filter(log2FC_enrichment >= 0) %>%
        write_tsv(results_up_out)
        # write_tsv(results_up_out) %>%
        # select(1:6) %>%
        # write_tsv(bed_up_out, col_names=FALSE)

    results_df_significant %>%
        filter(log2FC_enrichment < 0) %>%
        write_tsv(results_down_out)
        # write_tsv(results_down_out) %>%
        # select(1:6) %>%
        # write_tsv(bed_down_out, col_names=FALSE)

    maplot = plot_ma(df_sig = results_df_significant,
                     df_nonsig = results_df_nonsignificant,
                     xvar = mean_counts,
                     yvar = log2FC_enrichment,
                     lfc = lfc,
                     condition = condition,
                     control = control)

    volcano = plot_volcano(df = results_df,
                           xvar = log2FC_enrichment,
                           yvar = log10_padj,
                           lfc = lfc,
                           alpha = alpha,
                           condition = condition,
                           control = control)

    qc_plots = arrangeGrob(mean_sd_plot_pre,
                           mean_sd_plot_post,
                           maplot,
                           volcano,
                           ncol=2)

    ggsave(qc_plots_out,
           plot = qc_plots,
           width = 16*1.5,
           height = 9*1.5,
           units="cm")
}

main(exp_table = snakemake@input[["exp_counts"]],
     spike_table = snakemake@input[["spike_counts"]],
     samples = snakemake@params[["samples"]],
     conditions = snakemake@params[["conditions"]],
     sample_type = snakemake@params[["sampletypes"]],
     norm = snakemake@wildcards[["norm"]],
     condition = snakemake@wildcards[["condition"]],
     control = snakemake@wildcards[["control"]],
     alpha = snakemake@params[["alpha"]],
     lfc = snakemake@params[["lfc"]],
     counts_norm_out = snakemake@output[["counts_norm"]],
     counts_rlog_out = snakemake@output[["counts_rlog"]],
     results_all_out = snakemake@output[["results_all"]],
     results_up_out = snakemake@output[["results_up"]],
     results_down_out = snakemake@output[["results_down"]],
     results_unchanged_out = snakemake@output[["results_nonsig"]],
     # bed_all_out = snakemake@output[["bed_all"]],
     # bed_up_out = snakemake@output[["bed_up"]],
     # bed_down_out = snakemake@output[["bed_down"]],
     # bed_unchanged_out = snakemake@output[["bed_nonsig"]],
     qc_plots_out = snakemake@output[["qc_plots"]])

library(tidyverse)
library(forcats)
library(viridis)
library(ggthemes)
library(ggrepel)

import = function(path){
    read_tsv(path) %>%
        mutate_at(vars(sample, status), funs(fct_inorder(., ordered=TRUE))) %>%
        return()
}

main = function(seq_len_dist_in, per_tile_in, per_base_qual_in,
                per_base_seq_in, per_base_n_in, per_seq_gc_in,
                per_seq_qual_in, adapter_content_in, seq_dup_in,
                seq_len_dist_out, per_tile_out, per_base_qual_out,
                per_base_seq_out, per_seq_gc_out, per_seq_qual_out,
                adapter_content_out, seq_dup_out){
    #damnit fastqc...why bin some of the data and then output in this shite format...
    length_distribution = import(seq_len_dist_in) %>%
        separate(length, into=c('a','b'), sep="-", fill="right", convert=TRUE) %>%
        mutate_at(vars(count), funs(if_else(is.na(b), ., ./2))) %>%
        gather(key, length, c(a,b)) %>%
        filter(!is.na(length)) %>%
        select(-key)

    nsamples = n_distinct(length_distribution$sample)

    per_tile_quality = import(per_tile_in) %>%
        mutate_at(vars(tile), funs(fct_inorder(as.character(.), ordered=TRUE)))

    per_base_qual = import(per_base_qual_in) %>%
        left_join(length_distribution, by=c("base"="length", "sample", "status")) %>%
        group_by(sample, status) %>%
        mutate_at(vars(count), funs(if_else(is.na(.), 0, .))) %>%
        mutate(n = lag(sum(count)-cumsum(count), default=sum(count))) %>%
        mutate_at(vars(n), funs(./max(n)))

    adapter_content = import(adapter_content_in)

    per_base_seq = import(per_base_seq_in) %>%
        left_join(import(per_base_n_in), by=c("base", "sample","status")) %>%
        rename(position=base, n=n_count) %>%
        gather(base, pct, -c(position, sample, status)) %>%
        mutate_at(vars(base), funs(toupper(.))) %>%
        mutate_at(vars(base), funs(fct_inorder(., ordered=TRUE)))

    per_seq_gc = import(per_seq_gc_in) %>%
        filter(count > 0) %>%
        group_by(sample, status) %>%
        mutate(norm_count = count/sum(count))

    per_seq_qual = import(per_seq_qual_in) %>%
        filter(count > 0) %>%
        group_by(sample, status) %>%
        mutate(norm_count = count/sum(count))

    duplication_levels = import(seq_dup_in) %>%
        mutate_at(vars(duplication_level), funs(fct_inorder(., ordered=TRUE)))

    #kmer_content = import(kmer_in)

    theme_standard = theme_light() +
        theme(text = element_text(size=12, color="black", face="bold"),
              axis.text = element_text(size=12, color="black"),
              axis.title = element_text(size=12, color="black", face="bold"),
              strip.placement = "outside",
              strip.background = element_blank(),
              strip.text = element_text(size=12, color="black", face="bold"),
              strip.text.y = element_text(angle=-180, hjust=1))

    length_dist_plot = ggplot(data = length_distribution %>%
                                     group_by(sample, status) %>%
                                     mutate(normcount = count/max(count)),
                                 aes(x=length, y=normcount)) +
        geom_col(fill="#114477") +
        scale_x_continuous(breaks=scales::pretty_breaks(n=6), name="read length (nt)") +
        scale_y_continuous(breaks=scales::pretty_breaks(n=2), name="normalized counts") +
        facet_grid(sample~status, scales="free_y", switch="y") +
        ggtitle("read length distributions") +
        theme_standard +
        theme(axis.text.y = element_text(size=10, face="plain"))

    ggsave(seq_len_dist_out, plot=length_dist_plot, width=26, height=2+2*nsamples, units="cm", limitsize=FALSE)

    tile_quality_plot = ggplot(data = per_tile_quality %>% filter(status=="raw"),
                               aes(x=base, y=tile, fill=mean)) +
        geom_raster() +
        scale_fill_viridis(direction=-1, guide=guide_colorbar(title="mean\nquality\nscore", barheight=10)) +
        scale_x_continuous(expand=c(0,0), name="cycle number", breaks=scales::pretty_breaks(n=6)) +
        ylab("flow cell tile") +
        ggtitle("per tile sequencing quality") +
        facet_grid(sample~status, scales="free_y", switch="y") +
        theme_standard +
        theme(axis.text.y = element_blank(),
              axis.ticks.y = element_blank(),
              panel.grid.major.y = element_blank(),
              strip.text.x = element_blank())

    ggsave(per_tile_out, plot=tile_quality_plot, width=24, height=2+2*nsamples, units="cm", limitsize=FALSE)

    per_base_qual_plot = ggplot(data=per_base_qual, aes(x=base, y=fct_rev(sample),
                                               height=n, fill=mean, color=mean)) +
        geom_tile() +
        scale_color_viridis(guide=FALSE, direction=-1) +
        scale_fill_viridis(guide=guide_colorbar(title="mean quality", barwidth=12,
                                                barheight=1, title.position = "top",
                                                title.hjust=0.5),
                           direction=-1) +
        scale_x_continuous(expand=c(0,0), name="read length (nt)", breaks=scales::pretty_breaks(n=6)) +
        facet_grid(sample~status, scales="free_y", switch="y") +
        ggtitle("per base sequencing quality",
                subtitle = expression("bar height " %prop% " fraction of reads")) +
        theme_standard +
        theme(axis.text.y = element_blank(),
              axis.title.y = element_blank(),
              plot.subtitle = element_text(size=12),
              legend.position="top",
              legend.margin = margin(0,0,0,0))

    ggsave(per_base_qual_out, plot=per_base_qual_plot, width=26, height=2.5+1.25*nsamples, units="cm", limitsize=FALSE)

    adapter_plot = ggplot(data = adapter_content, aes(x=position, y=0, fill=pct, color=pct)) +
        geom_raster() +
        scale_color_viridis(guide=FALSE) +
        scale_fill_viridis(guide=guide_colorbar(title="% reads with adapter", barwidth=12,
                                                barheight=1, title.position = "top",
                                                title.hjust=0.5)) +
        scale_x_continuous(expand=c(0,0), name="read length (nt)") +
        scale_y_continuous(expand=c(0,0), name=NULL, breaks=0, labels=NULL) +
        facet_grid(sample~status, switch="y") +
        ggtitle("adapter content") +
        theme_standard +
        theme(legend.position="top",
              legend.margin = margin(0,0,0,0))

    ggsave(adapter_content_out, plot=adapter_plot, width=32, height=2+1.25*nsamples, units="cm", limitsize=FALSE)

    per_base_seq_plot = ggplot(data = per_base_seq, aes(x=position, y=pct, color=base)) +
        geom_line() +
        scale_color_ptol(guide=guide_legend(label.position="top", label.hjust=0.5,
                                            keyheight=0.2)) +
        scale_x_continuous(expand=c(0,1), name="position in read", breaks=scales::pretty_breaks(n=6)) +
        scale_y_continuous(name="% of reads", breaks=scales::pretty_breaks(n=2)) +
        facet_grid(sample~status, switch="y") +
        ggtitle("per base sequence content") +
        theme_standard +
        theme(legend.position="top",
              legend.title = element_blank(),
              legend.margin = margin(0,0,0,0),
              legend.key.size = unit(1, "cm"),
              legend.text = element_text(size=12, face="bold"),
              axis.text.y = element_text(size=10, face="plain"))

    ggsave(per_base_seq_out, plot=per_base_seq_plot, width=32, height=2+2.25*nsamples, units="cm", limitsize=FALSE)

    per_seq_gc_plot = ggplot(data = per_seq_gc, aes(x=gc_content, y=norm_count)) +
        geom_line(color="#114477") +
        scale_x_continuous(expand=c(0,0), name="GC%") +
        #xlab("GC%") +
        scale_y_continuous(breaks=scales::pretty_breaks(n=2), name="normalized counts") +
        facet_grid(sample~status, switch="y") +
        ggtitle("per sequence GC content") +
        theme_standard +
        theme(axis.text.y = element_text(size=10, face="plain"),
              panel.spacing.x = unit(1, "cm"),
              plot.margin = margin(5.5, 12, 5.5, 5.5, unit="pt"))

    ggsave(per_seq_gc_out, plot=per_seq_gc_plot, width=26, height=2+2*nsamples, units="cm", limitsize=FALSE)

    per_seq_qual_plot = ggplot(data = per_seq_qual, aes(x=quality, y=norm_count)) +
        geom_col(fill="#114477") +
        scale_x_continuous(breaks=scales::pretty_breaks(n=5), name="quality score") +
        scale_y_continuous(breaks=scales::pretty_breaks(n=2), name="normalized counts") +
        facet_grid(sample~status, switch="y") +
        ggtitle("per sequence quality scores") +
        theme_standard +
        theme(axis.text.y = element_text(size=10, face="plain"))

    ggsave(per_seq_qual_out, plot=per_seq_qual_plot, width=26, height=2+1.5*nsamples, units="cm", limitsize=FALSE)

    dup_level_plot = ggplot(data = duplication_levels, aes(x=duplication_level, y=pct_of_total)) +
        geom_col(fill="#114477") +
        xlab("duplication level") +
        ylab("% of total reads") +
        facet_grid(sample~status, switch="y") +
        ggtitle("sequence duplication levels") +
        theme_standard +
        theme(axis.text.x = element_text(size=10, face="plain", angle=60, hjust=1),
              axis.text.y = element_text(size=10, face="plain"))

    ggsave(seq_dup_out, plot=dup_level_plot, width=26, height=2+1.5*nsamples, units="cm", limitsize=FALSE)

    ##ermmm...no obvious way to make this one look nice, but then it doesn't really need to
    #kmer_content_plot = ggplot(data = kmer_content, aes(x=max_position, y=log2(obs_over_exp_max), label=sequence)) +
    #    geom_point(shape=16, stroke=0, size=1, alpha=0.5) +
    #    geom_label_repel(size=2, label.size=unit(0.05, "pt"), label.padding=unit(0.1, "pt"), label.r=unit(0,"pt"), segment.size=0.1,
    #                     box.padding=unit(0.05,"pt"), segment.alpha=0.4) +
    #    xlab("position in read") +
    #    ylab(expression(bold(log[2]~ frac("observed", "expected")))) +
    #    ggtitle("k-mer content",
    #            subtitle = "top 20 overrepresented k-mers") +
    #    facet_grid(sample~status, switch="y", scales="free_y") +
    #    theme_standard + theme(plot.subtitle = element_text(size=12, face="plain"))
    #
    #ggsave(kmer_out, plot=kmer_content_plot, width=35, height=2+5*nsamples, units="cm", limitsize=FALSE)
}

main(seq_len_dist_in = snakemake@input[["seq_length_dist"]],
     per_tile_in = snakemake@input[["per_tile_qual"]],
     per_base_qual_in = snakemake@input[["per_base_qual"]],
     per_base_seq_in = snakemake@input[["per_base_seq_content"]],
     per_base_n_in = snakemake@input[["per_base_n"]],
     per_seq_gc_in = snakemake@input[["per_seq_gc"]],
     per_seq_qual_in = snakemake@input[["per_seq_qual"]],
     adapter_content_in = snakemake@input[["adapter_content"]],
     seq_dup_in = snakemake@input[["seq_duplication"]],
     #kmer_in = snakemake@input[["kmer"]],
     seq_len_dist_out = snakemake@output[["seq_length_dist"]],
     per_tile_out = snakemake@output[["per_tile_qual"]],
     per_base_qual_out = snakemake@output[["per_base_qual"]],
     per_base_seq_out = snakemake@output[["per_base_seq_content"]],
     per_seq_gc_out = snakemake@output[["per_seq_gc"]],
     per_seq_qual_out = snakemake@output[["per_seq_qual"]],
     adapter_content_out = snakemake@output[["adapter_content"]],
     seq_dup_out = snakemake@output[["seq_duplication"]])
     #kmer_out = snakemake@output[["kmer"]])