from Bio import SeqIO
from sys import stderr
import csv
import re
from crc64iso.crc64iso import crc64

info_file = str(snakemake.input['info'])
a2m_file = str(snakemake.input['a2m'])
pos_file = str(snakemake.input['pos'])
usearch_file = str(snakemake.input['usearch'])
blasso_file  = str(snakemake.input['blasso'])
out_file = str(snakemake.output)

infos = dict()
data_type = 'JGI'
header = []
id_cols = { 'OM-RGC_ID', 'genome_id', 'gene_oid', 'ID' }
with open(info_file) as fh:
    tsv = csv.reader(fh, delimiter = '\t')
    header = next(tsv)
    header = [v + str(header[:i].count(v) + 1) if header.count(v) > 1 else v for i, v in enumerate(header)]
    id_cols = set(header) & id_cols
    assert len(id_cols) > 0, "Id column not found"
    id_col = id_cols.pop()
    has_sequence = 'sequence' in header
    for line in tsv:
        row = dict(zip(header, line))
        record_id = row.pop(id_col).replace(' ', '_')
        if has_sequence:
            del row['sequence']
        infos[record_id] = row
    header.remove(id_col)
    if has_sequence:
        header.remove('sequence')

positions = []
groups = []
with open(pos_file) as fh:
    for line in fh:
        pos, group = line.rstrip().split()
        pos0 = int(pos) - 1
        positions.append((pos0, group))
        if len(group) > 1 and group not in groups:
            groups.append(group)

blassos = dict()
with open(blasso_file) as fh:
    for line in fh:
        record_id, wl_mean, wl_sd = line.split()
        blassos[record_id] = { 'wl_mean': wl_mean, 'wl_sd': wl_sd }

matches = dict()
with open(usearch_file) as fh:
    for line in fh:
        query, target, ident = line.split('\t')
        if query not in matches:
            target_id, *rest = target.split()
            search = re.search('{(.+?)}', target)
            assert search, "Unexpected reference header format: %s" % target
            family = search.group(1)
            matches[query] = { 'target': target, 'family': family, 'ident': float(ident) }

aminoacids = [ 'A','C','D','E','F','G','H','I','K','L','M','N','P','Q','R','S','T','V','W','Y' ]
with open(out_file, 'w') as out_fh:
    tsv = csv.writer(out_fh, delimiter = "\t")
    tsv.writerow([ 'record_id', 'checksum', 'target', 'family', 'ident', 'wl_mean', 'wl_sd' ] + groups + header)
    with open(a2m_file) as a2m_fh:
        a2m = SeqIO.parse(a2m_fh, 'fasta')
        ref = next(a2m)
        ref_pos = 0
        aln_pos = []
        for pos in range(len(ref.seq)):
            if ref.seq[pos] not in [ '-', '.' ]:
                for pos0, group in positions:
                    if ref_pos == pos0:
                        aln_pos.append((pos, group))
                ref_pos += 1
        for rec in a2m:
            if rec.id in matches:
                seq = re.sub('[^A-Z]', '', str(rec.seq))
                cks = crc64(seq)
                match  = matches[rec.id]
                blasso = blassos[rec.id]
                info_id, *rest = rec.id.split('/')
                info   = infos[info_id]
                res = { group: [] for group in groups }
                for pos, group in aln_pos:
                    if group in groups:
                        res[group].append(rec.seq[pos])
                    elif rec.seq[pos] != group:
                        break
                else:
                    row = [ rec.id, cks, match['target'], match['family'], match['ident'], blasso['wl_mean'], blasso['wl_sd'] ]
                    row += [ ''.join(res[x]) for x in groups ]
                    row += info.values()
                    tsv.writerow(row)

snakemake@source("plot_functions.R")

with(snakemake@input, {
    metadata_file <<- metadata
    tsv_file  <<- tsv
    crt_files <<- crt
    colors_file <<- colors
})
with(snakemake@output, {
    tax_file <<- tax
    sum_file <<- sum
    car_file <<- car
})

aliases <- list(
    BCD = "BLH",
    Lycopene_cycl = "CrtY",
    crtI_fam = "CrtI"
)

metadata <- read.table(metadata_file, header = T, sep = "\t", comment.char = "") %>%
    filter(ecosystem_category == "Aquatic", habitat != "Groundwater") %>%
    mutate(Ecosystem.Type = ecosystem_type, Specific.Ecosystem = habitat, Ecosystem.Subtype = habitat, Isolation = habitat, x = longitude, y = latitude) %>%
    my_ecosystem

taxa <- read.table(colors_file, sep = "\t", comment.char = "") %>%
    arrange(V1) %>%
    with(setNames(V2,V1))

genes <- lapply(crt_files, read.table) %>%
    bind_rows %>%
    select(1, 3, 5) %>%
    setNames(c("record_id", "gene", "E_value")) %>%
    filter(E_value < 1e-15) %>%
    separate(record_id, into = c("genome_id", "Gene.ID"), sep = "/") %>%
    mutate(gene = recode(gene, !!!aliases)) %>%
    distinct(genome_id, gene)

all_data <- read.table(tsv_file, header = T, sep = "\t", quote = "", comment.char = "") %>%
    rename(Taxonomy = "ecosystem") %>% # fix bug
    filter(ecosystem_category == "Aquatic", !grepl("-", motif)) %>%
    mutate(Ecosystem.Type = ecosystem_type, Specific.Ecosystem = habitat, Ecosystem.Subtype = habitat, Isolation = habitat, x = longitude, y = latitude, Abundance = 1) %>%
    my_ecosystem %>% 
    my_rhodopsin_class %>%
    separate(record_id, into = c("genome_id", "Gene.ID"), sep = "/") %>%
    mutate(window_type = case_when(
        window %in% c("W", "F") ~ "FW",
        window %in% c("G") ~ "G",
        T ~ "other"
    )) %>%
    distinct(genome_id, class, window, .keep_all = T)

data <- complete(all_data, class, window, nesting(x, y, Ecosystem), fill = list(Abundance = 0)) %>%
    group_by(x, y, Ecosystem, class, window) %>%
    summarize(Abundance = sum(Abundance), .groups = "drop")

by_ecosystem <- group_by(data, class, window, Ecosystem) %>%
    summarize(Abundance = sum(Abundance), .groups = "drop") %>%
    group_by(Ecosystem) %>%
    group_split %>%
    lapply(my_logo_matrix)

aliases <- list(Cyanobacteriota = "Cyanobacteria")
by_taxa <- distinct(all_data, otu_id, .keep_all = T) %>%
    separate_rows(Taxonomy, sep = ";") %>%
    filter(class %in% c("x", "p"), window_type %in% c("FW", "G")) %>%
    separate(Taxonomy, into = c("Rank", "Name"), sep = "__") %>%
    filter(Rank == "p") %>%
    group_by(class, window_type, Ecosystem, Name) %>%
    summarize(Abundance = sum(Abundance), .groups = "drop") %>%
    group_by(Name) %>%
    mutate(Name = ifelse(sum(Abundance) > 10, Name, NA)) %>%
    mutate(Name = recode(Name, !!!aliases))

tmp <- left_join(metadata, all_data, by = "genome_id", suffix = c("", ".all_data")) %>%
    distinct(genome_id, .keep_all = T) %>%
    group_by(Ecosystem) %>%
    summarize(num_px = sum(class %in% c("x", "p") & window_type %in% c("FW", "G")), num_all = n())

all_taxa <- unique(na.omit(by_taxa$Name))
missing_taxa <- all_taxa[! all_taxa %in% names(taxa)]
if (length(missing_taxa) > 0) {
    write("Taxa not in the color file:", stderr())
    write(missing_taxa, stderr())
    q()
}

by_gene <- all_data %>%
    filter(class %in% c("x", "p"), window_type %in% c("WF", "G")) %>%
    left_join(genes, by = "genome_id") %>%
    mutate(present = T) %>%
    complete(nesting(genome_id, class, window_type), gene, fill = list(present = F)) %>%
    filter(!is.na(gene)) %>%
    group_by(class, window_type, gene, present) %>%
    summarize(n = n())
p <- ggplot(by_gene, aes(x = gene, fill = present, y = n)) +
    geom_bar(position = "stack", stat = "identity") +
    facet_grid(class ~ window_type) +
    theme_bw()
ggsave(car_file, p)

p <- ggplot(by_taxa, aes(fill = Name, x = Ecosystem, y = Abundance)) +
    geom_bar(position = "stack", stat = "identity") +
    scale_fill_manual(values = taxa) +
    facet_grid(class ~ window_type) +
    theme_bw()
ggsave(tax_file, p, width = 7, height = 5)

p <- my_ggplot_ecosystems()
p <- Reduce(my_add_grob_to_ecosystem, by_ecosystem, p)
ggsave(sum_file, p, height = 8, width = 16)

snakemake@source("plot_functions.R")
library(readxl)
library(car)
library(lme4)

with(snakemake@input, {
    tsv_file       <<- tsv
    metadata_file  <<- metadata
    abundance_file <<- abundance
    crtZ_file      <<- crtZ
    crtZ_abundance_file <<- crtZ_abundance
    single_copy_files <<- single_copy
})
with(snakemake@params, {
    sheet_name <<- sheet
})
with(snakemake@output, {
    map_file <<- map
    map_percell_file <<- map_percell
    sum_file <<- sum
    dep1_file <<- dep1
    lat1_file <<- lat1
    test1_file <<- test1
    dep2_file <<- dep2
    lat2_file <<- lat2
    test2_file <<- test2
})

metadata <- read_excel(metadata_file, sheet_name, .name_repair = "universal") %>%
    mutate(x = round(Longitude), y = round(Latitude)) %>%
    mutate(Ecosystem = factor("Marine"))
abundance <- read.table(abundance_file, header = T) %>%
    gather(Sample, Abundance, -OMRGC_ID) %>%
    left_join(metadata, by = c(Sample = "PANGAEA.sample.id"))
all.data <- read.table(tsv_file, header = T, sep = "\t", na.strings = c("NA", "")) %>%
    filter(!is.na(wl_mean)) %>%
    my_rhodopsin_class %>%
    left_join(abundance, by = c(record_id = "OMRGC_ID")) %>%
    mutate(Depth = suppressWarnings(as.numeric(Depth..nominal)), fenestrated = case_when(window == "G" ~ T, window %in% c("W","F") ~ F, T ~ NA))
single_copy <- lapply(single_copy_files, read.table, header = T) %>%
    lapply(select, -OMRGC_ID) %>%
    lapply(colSums) %>%
    bind_rows %>%
    colMeans %>%
    data.frame(Abundance = ., Sample = names(.)) %>%
    filter(Sample %in% all.data$Sample) %>%
    left_join(metadata, by = c(Sample = "PANGAEA.sample.id"))

motifs_counts <- all.data %>%
    distinct(checksum, .keep_all = T) %>%
    group_by(family, motif) %>%
    summarize(Abundance = n())
motifs_weighted <- all.data %>%
    group_by(family, motif) %>%
    summarize(Abundance = sum(Abundance))

# NB: these plots are not saved
#p <- ggplot(motifs_counts, aes(fill = motif, x = family, y = Abundance)) +
#    geom_bar(position = "stack", stat = "identity")
#p <- ggplot(motifs_weighted, aes(fill = motif, x = family, y = Abundance)) +
#    geom_bar(position = "stack", stat = "identity")

ratio.data <- filter(all.data, class %in% c("x", "p"), !is.na(fenestrated), !is.na(Depth)) %>%
    # filter(Polar == "Non polar") %>%
    group_by(Sample, x, y, Ecosystem, Depth, Station, Polar, fenestrated) %>%
    summarize(Abundance = sum(Abundance), .groups = "drop") %>%
    spread(fenestrated, Abundance, fill = 0) %>%
    mutate(Abundance = `TRUE` + `FALSE`, ratio = `TRUE` / Abundance)
model <- lmer(ratio ~ log10(Depth + 1) + abs(y) + (1 | Station), data = ratio.data)
capture.output(Anova(model, type = "II"), file = test1_file)

p <- ggplot(ratio.data, aes(x = Depth, y = ratio, size = Abundance)) +
    geom_point() +
    geom_smooth(method = 'lm', formula = y ~ x) +
    scale_x_continuous(trans = 'log10')
ggsave(dep1_file, p)
p <- ggplot(ratio.data, aes(x = abs(y), y = ratio, size = Abundance)) +
    geom_point() +
    geom_smooth(method = 'lm', formula = y ~ x)
ggsave(lat1_file, p)

weighted_mean <- function(x, w) {
    sum(w * x) / sum(w)
}
weighted_sd <- function(x, w) {
    x_mean <- weighted_mean(x, w)
    M <- sum(w > 0)
    sqrt( sum(w * (x - x_mean)^2) * M / (M-1) / sum(w) )
}

wl.data <- filter(all.data, class %in% c("x", "p"), !is.na(fenestrated), !is.na(Depth)) %>%
    group_by(x, y, Ecosystem, Depth, Station, Polar) %>%
    summarize(wl_avg = weighted_mean(wl_mean, Abundance), wl_sd = weighted_sd(wl_mean, Abundance), Abundance = sum(Abundance), .groups = "drop")
model <- lmer(wl_avg ~ log10(Depth + 1) + abs(y) + (1 | Station), data = wl.data)
capture.output(Anova(model, type = "II"), file = test2_file)

p <- ggplot(wl.data, aes(x = Depth, y = wl_avg, size = Abundance)) +
    geom_pointrange(aes(ymin = wl_avg - wl_sd, ymax = wl_avg + wl_sd)) +
    geom_smooth(method = 'lm', formula = y ~ x) +
    scale_size_continuous(range = c(0.1,0.5)) +
    scale_x_continuous(trans = 'log10')
ggsave(dep2_file, p)
p <- ggplot(wl.data, aes(x = abs(y), y = wl_avg, size = Abundance)) +
    geom_pointrange(aes(ymin = wl_avg - wl_sd, ymax = wl_avg + wl_sd)) +
    scale_size_continuous(range = c(0.1,0.5)) +
    geom_smooth(method = 'lm', formula = y ~ x)
ggsave(lat2_file, p)

crtZ.abundance <- read.table(crtZ_abundance_file, header = T) %>%
    gather(Sample, Abundance, -OMRGC_ID) %>%
    left_join(metadata, by = c(Sample = "PANGAEA.sample.id"))
crtZ.data <- read.table(crtZ_file) %>%
    select(OMRGC_ID = "V1", evalue = "V5") %>%
    filter(evalue < 1e-15) %>%
    left_join(crtZ.abundance, by = "OMRGC_ID") %>%
    group_by(Sample, x, y, Station) %>%
    summarize(CrtZ.abundance = sum(Abundance), .groups = "drop") %>%
    left_join(select(ratio.data, Sample, Depth, Abundance, `TRUE`, `FALSE`), by = "Sample") %>%
    mutate(ratio = CrtZ.abundance / `TRUE`)

model <- lmer(ratio ~ log10(Depth + 1) + abs(y) + (1 | Station), data = crtZ.data)
capture.output(Anova(model, type = "II"), file = test2_file)

p <- ggplot(crtZ.data, aes(x = CrtZ.abundance, y = `TRUE`)) +
    geom_point() +
    geom_smooth(method = 'lm', formula = y ~ x) +
    scale_x_continuous(trans = 'log10') +
    scale_y_continuous(trans = 'log10')
ggsave("TRUE.pdf", p)

p <- ggplot(crtZ.data, aes(x = CrtZ.abundance, y = `FALSE`)) +
    geom_point() +
    geom_smooth(method = 'lm', formula = y ~ x) +
    scale_x_continuous(trans = 'log10') +
    scale_y_continuous(trans = 'log10')
ggsave("FALSE.pdf", p)

p <- ggplot(crtZ.data, aes(x = Depth, y = ratio, size = Abundance)) +
    geom_point() +
    geom_smooth(method = 'lm', formula = y ~ x) +
    scale_x_continuous(trans = 'log10') +
    scale_y_continuous(trans = 'log10')
ggsave("tmp.pdf", p)

# prots <- distinct(all.data, record_id, .keep_all=T)
# p <- ggplot(prots, aes(x = window, y = wl_mean)) + geom_boxplot()
# ggsave(win_wl_file, p)

data <- group_by(all.data, x, y, Ecosystem, class, window) %>%
    summarize(Abundance = sum(Abundance), .groups = "drop") %>%
    complete(class, window, nesting(x, y, Ecosystem), fill = list(Abundance = 0))
single_copy_data <- group_by(single_copy, x, y, Ecosystem) %>%
    summarize(Abundance = sum(Abundance), .groups = "drop")

by_location <- left_join(data, single_copy_data, by = c("x", "y", "Ecosystem"), suffix = c("", ".SC")) %>%
    group_by(x, y, Ecosystem) %>%
    filter(class %in% c("p", "x")) %>%
    summarize(
        Total.abundance = sum(Abundance),
        Total.percell = Total.abundance / first(Abundance.SC),
        G  = sum(Abundance[window %in% "G"],         na.rm = T),
        WF = sum(Abundance[window %in% c("W", "F")], na.rm = T),
        .groups = "drop"
    )

by_ecosystem <- group_by(data, class, window, Ecosystem) %>%
    summarize(Abundance = sum(Abundance), .groups = "drop") %>%
    my_logo_matrix

width <- 16
height <- 8

p <- mutate(by_location, Total = log10(Total.abundance)) %>%
    my_pies(width, height, labeller = function(r) sprintf('1e%d', r))
ggsave(map_file, p, height = height, width = width)

scale <- 10
p <- mutate(by_location, Total = Total.percell * scale) %>%
    my_pies(width, height, labeller = function(r) sprintf('%d%%', r / scale * 100))
ggsave(map_percell_file, p, height = height, width = width)

p <- my_ggplot_ecosystems()
p <- my_add_grob_to_ecosystem(p, by_ecosystem)
ggsave(sum_file, p, height = 3, width = 4)

#### Calculation of % of fenestrated rhods per cell in a sample ####
#### not yet integrated ####
data <- group_by(all.data, Sample, Layer, class, window) %>%
    summarize(Abundance = sum(Abundance), .groups = "drop") %>%
    complete(class, window, Sample, fill = list(Abundance = 0))
single_copy_data <- group_by(single_copy, Sample) %>%
    summarize(Abundance = sum(Abundance), .groups = "drop")

by_location <- left_join(data, single_copy_data, by = "Sample", suffix = c("", ".SC")) %>%
    filter(class %in% c("p", "x"), window %in% c("G","W","F")) %>%
    group_by(Sample, Layer) %>%
    summarize(
        Total.abundance = sum(Abundance),
        Total.percell = Total.abundance / first(Abundance.SC),
        G  = sum(Abundance[window %in% "G"],         na.rm = T),
        WF = sum(Abundance[window %in% c("W", "F")], na.rm = T),
        .groups = "drop"
    ) %>%
    mutate(G.percell = Total.percell * G / Total.abundance) %>%
    group_by(Layer) %>%
    summarize(G.percell.mean = mean(G.percell), G.percell.sd = sd(G.percell))

library(treeio)
library(ggtree)
library(ape)
library(phytools)
library(dplyr)
library(tidyr)
library(ggplot2)
library(phangorn)
library(ggnewscale)
library(castor)
library(seqinr)

if (interactive()) {
    Snakemake <- setClass("Snakemake", slots = list(input = "list", output = "list"))
    snakemake <- Snakemake(
        input = list(outgroup = "input/outgroups.fasta", tree = "analysis/phylogeny/rhodopsins.treefile", metadata = "analysis/parse/phylogeny.tsv"),
        output = list("tmp.pdf")
    )
}

with(snakemake@input, {
    outgroup_file <<- outgroup
    tree_file     <<- tree
    tsv_file      <<- tsv
    colors_file   <<- colors
    a2m_file      <<- a2m
})
with(snakemake@output, {
    output_file_small <<- small
    output_file_big   <<- big
    output_jtree      <<- jtree
})

taxa <- read.table("metadata/taxa.txt", sep = "\t", comment.char = "") %>%
    arrange(1) %>%
    with(setNames(V2, V1))

outgroups <- names(read.fasta(outgroup_file))
tsv <- read.table(tsv_file, header = T, sep = "\t", na.strings = "", fill = T) %>%
    select(-target)
metadata <- read.fasta(a2m_file, seqtype = "AA", as.string = T) %>%
    {data.frame(label = names(.), sequence = as.character(.))} %>%
    left_join(tsv, by = c(label = "record_id")) %>%
    mutate(is_outgroup = label %in% outgroups) %>%
    mutate(Alias = gsub(",.+", "", Alias)) %>%
    mutate(Activity = ifelse(grepl("\\]$", Activity), NA, gsub("[][]", "", Activity))) %>%
    mutate(Highlight = !is.na(Highlight)) %>%
    mutate(D85 = substr(motif, 1, 1), T89 = substr(motif, 2, 2), D96 = substr(motif, 3, 3), G156 = window)

tree <- read.tree(tree_file)

tree.unrooted <- as_tibble(tree) %>%
    left_join(metadata, by = "label") %>%
    `class<-`(c("tbl_tree", "data.frame")) %>%
    as.treedata
write.jtree(tree.unrooted, file = output_jtree)

tree.tib <- ape::root(tree, outgroups, edgelabel = T, resolve.root = T) %>%
    drop.tip(outgroups) %>%
    as_tibble %>%
    mutate(support = suppressWarnings(as.numeric(label))) %>%
    left_join(metadata, by = "label") %>%
    mutate(is_outgroup = ifelse(label %in% outgroups, T, NA)) %>%
    `class<-`(c("tbl_tree", "data.frame"))
tree.phylo <- as.treedata(tree.tib)@phylo

clustalx <- c(
    A = "BLUE",
    I = "BLUE",
    L = "BLUE",
    M = "BLUE",
    F = "BLUE",
    W = "BLUE",
    V = "BLUE",
    C = "BLUE",
    K = "RED",
    R = "RED",
    E = "MAGENTA",
    D = "MAGENTA",
    N = "GREEN",
    Q = "GREEN",
    S = "GREEN",
    T = "GREEN",
    G = "ORANGE",
    P = "YELLOW",
    H = "CYAN",
    Y = "CYAN"
)

clades <- filter(tree.tib, ! node %in% parent) %>%
    pull(Clade) %>%
    as.factor
hsp <- hsp_max_parsimony(tree.phylo, as.numeric(clades)) %>%
    `$`("likelihoods") %>%
    data.frame(check.names = F) %>%
    mutate(node = row_number()) %>%
    gather(hsp, prob, -node) %>%
    filter(prob > 0.9) %>%
    mutate(hsp = levels(clades)[as.numeric(hsp)])

tree <- left_join(tree.tib, hsp, by = "node") %>%
    mutate(hsp.parent = .[match(parent, node),"hsp"]) %>%
    mutate(hsp = ifelse(!is.na(hsp.parent) & hsp == hsp.parent, NA, hsp)) %>%
    as.treedata

p_small <- ggtree(tree, layout = "circular") +
    geom_highlight(mapping = aes(subset = !is.na(hsp), fill = hsp), alpha = 0.1) +
    geom_treescale() +
    geom_tiplab(mapping = aes(subset = !is.na(Alias), label = Alias, size = Highlight), offset = 0.1) +
    # geom_tippoint(mapping = aes(subset = !is.na(Activity), color = Activity), size = 1) +
    scale_size_manual(values = c(2.5, 5)) + new_scale("size") + # labels
    geom_point2(aes(subset = !is.na(support) & support >= 90, size = support >= 95), color = "darkgray") +
    scale_size_manual(values = c(0.5, 1)) + # support values
    new_scale("color") +
    geom_text2(aes(label = G156, color = G156, angle = angle - 90, x = 4.3), size = 3) +
    scale_color_manual(values = clustalx) # residues

p_big <- ggtree(tree, aes(color = Taxon), layout = "rectangular") +
    geom_highlight(mapping = aes(subset = !is.na(hsp), fill = hsp), alpha = 0.1) +
    scale_color_manual(values = taxa) + new_scale("color") +
    geom_treescale() +
    geom_tiplab(mapping = aes(label = sprintf("%s [%s]", ifelse(is.na(Alias), label, Alias), Organism)), size = 2, offset = 0.1) +
    geom_tippoint(mapping = aes(subset = !is.na(Activity), color = Activity), size = 1) +
    geom_point2(aes(subset = !is.na(support) & support >= 90, size = support >= 95, x = branch), shape = 15, color = "darkgray") +
    scale_size_manual(values = c(0.5, 1)) + # support values
    new_scale("color") +
    geom_text2(aes(label = D85, color = D85, x = 4.8), size = 1.5) +
    geom_text2(aes(label = T89, color = T89, x = 5.0), size = 1.5) +
    geom_text2(aes(label = D96, color = D96, x = 5.2), size = 1.5) +
    geom_text2(aes(label = G156, color = G156, x = 5.4), size = 2) +
    scale_color_manual(values = clustalx) # residues

ggsave(output_file_small, p_small, height = 7, width = 8)
ggsave(output_file_big,   p_big,   height = 6.5, width = 5)