gene_vuln_comparison_plots_parallel.R

# Vulnerability modeling - parallel

# Imports ----
source("vulnerability_modeling_tools.R")
source("vulnerability_plotting_tools.R")

require("optparse")

# Options ----
Sys.setenv(LOCAL_CPPFLAGS = '-march=corei7 -mtune=corei7')
options(mc.cores = parallel::detectCores())



# Main ----

args = commandArgs(trailingOnly=TRUE)

# test if there is at least one argument: if not, return an error

option_list = list(
    # Data stuff:
    make_option("--strain1", type="character", default="",
              help="Strain for condition 1. Used as reference."),
    make_option("--label1", type="character", default="",
              help="Label for condition 1."),
    make_option("--data1", type="character", default="",
              help="Data for condition 1."),

    make_option("--strain2", type="character", default="",
              help="Strain for condition 2. Used as reference."),
    make_option("--label2", type="character", default="",
              help="Label for condition 2."),
    make_option("--data2", type="character", default="",
              help="Data for condition 2."),

    make_option("--dir", type="character", default="./",
              help="Base data dir"),
    make_option("--genes", type="character", default="",
              help="Gene path"),


    #System stuff:
    make_option("--cores", type="numeric", default=15,
              help="Number of cores to use"),


    # plotting stuff:
    make_option("--ymax", type="numeric", default=0.05,
              help="Max value for the y-axis"),
    make_option("--ymin", type="numeric", default=-1.0,
              help="minimum for the y-axis"),

    make_option("--border1", type="character", default="#262626",
              help="Border color for condition 1"),
    make_option("--fill1", type="character", default="#262626",
              help="Fill color for condition 1"),

    make_option("--border2", type="character", default="#609393",
              help="Border color for condition 2"),
    make_option("--fill2", type="character", default="#008B8B",
              help="Fill color for condition 2")

);

opt_parser = OptionParser(option_list=option_list);
opt = parse_args(opt_parser);
  
desired_strain1 = opt$strain1
label1 = opt$label1

desired_strain2 = opt$strain2
label2 = opt$label2

desired_gene_path = opt$genes

if (file.exists(opt$data1)){
    DATA_PATH1 = opt$data1
} else {
   DATA_PATH1 = file.path(opt$dir, paste0("data_", label1, ".txt"))
}

if (file.exists(opt$data2)){
    DATA_PATH2 = opt$data2
} else {
   DATA_PATH2 = file.path(opt$dir, paste0("data_", label2, ".txt"))
}



read_data = function(path){
    XFULL = read.table(path, sep=",", stringsAsFactors = FALSE, header=T)
    ii_bad = is.na(XFULL$ESS)
    XFULL$ESS[ii_bad] = "N/A"
    return(XFULL)
}

message(sprintf("Reading %s data...", label1))
data1 = read_data(DATA_PATH1)
message(sprintf("Reading %s data...", label2))
data2 = read_data(DATA_PATH2)


if (file.exists(desired_gene_path)){
        unique_genes = scan(desired_gene_path, what="character")
} else {
    unique_genes = unique(c(unique(data1$ORF), unique(data2$ORF)))
}
 

N = length(unique_genes)
message(sprintf("Plotting comparison plots on all %s genes...", N))
# run_model("MSMEG4621:proB", data, desired_strain, folder=label)



output_prefix = file.path(opt$dir, "plots/comparisons", paste(label2, "vs", label1, sep="_"))


gene_comparisons = list()
for(i in 1:N){

  gene = unique_genes[i]

  gene_comparisons[[i]] = list(
        list(strain=desired_strain1, gene=gene, label = label1,
            folder=label1,
            fill_palette = opt$fill1, border_palette=opt$border1),

        list(strain=desired_strain2, gene=gene, label = label2,
            folder=label2,
            fill_palette = opt$fill2, border_palette=opt$border2)
    )
  
}         

if(opt$cores > 1){
    mclapply(gene_comparisons, plot_comparisons, mc.cores=opt$cores,
        prefix=opt$dir,
        output_prefix=output_prefix,
        width=13, height=13, band_alpha=0.1, style="lines",
        str_resolution=0.01, summary_size=4, samples=1000,
        )
} else {
    for(gene_comp in gene_comparisons){
        message("Comparing:")
        message(gene_comp)
        message("\n\n")
        plot_comparisons(gene_comp, prefix=opt$dir,
            output_prefix=output_prefix,
            width=13, height=13, band_alpha=0.1, style="lines",
            str_resolution=0.01, summary_size=4, samples=1000)

    }
}