<commit_before>__author__ = 'boris'
from scipy.stats.mstats import winsorize
import scipy.stats as stats
import matplotlib.pyplot as plt
from collections import defaultdict
import numpy
import math
import mod_lib
import mod_utils
import operator
import os
import uniform_colormaps
from statsmodels.nonparametric.smoothers_lowess import lowess
from matplotlib import rc
plt.rcParams['pdf.fonttype'] = 42 #leaves most text as actual text in PDFs, not outlines
#plt.rcParams["font.family"] = "sans-serif"
rc('font',**{'family':'sans-serif','sans-serif':['Bitstream Vera Sans']})
#rc('text', usetex=True)

def plot_mutated_nts_pie(libraries, out_prefix, subtract_background=False, subtract_control=False, exclude_constitutive=False):
    #Makes an array of pie charts, 1 per library
    if subtract_background:
        #if subtracting background, need to only look at those which have a defined control
        libraries = [library for library in libraries if (library.lib_settings.sample_name in
                     library.experiment_settings.get_property('experimentals')) or (library.lib_settings.sample_name in
                     library.experiment_settings.get_property('with_mod_controls'))]
    elif subtract_control:
        #if subtracting background, need to only look at those which have a defined control
        libraries = [library for library in libraries if library.lib_settings.sample_name in
                     library.experiment_settings.get_property('experimentals')]
    num_subplots = len(libraries)
    num_plots_wide = math.ceil(math.sqrt(num_subplots))
    num_plots_high = num_plots_wide
    fig = plt.figure(figsize=(4*num_plots_wide, 4*num_plots_high))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    plot_index =1
    for library in libraries:
        plot = fig.add_subplot(num_plots_high, num_plots_wide, plot_index)
        mutated_nts_count = library.count_mutation_rates_by_nucleotide(subtract_background=subtract_background, subtract_control=subtract_control,
                                                                       exclude_constitutive=exclude_constitutive)
        labels = sorted(mutated_nts_count.keys())
        sizes = numpy.array([mutated_nts_count[nt] for nt in labels])
        total = float(sum(sizes))
        sizes = sizes/total
        merged_labels = ['%s %.3f' % (labels[i], sizes[i]) for i in range(len(sizes))]
        plot.pie(sizes, labels = merged_labels, colors = mod_utils.rainbow)
        plot.set_title(library.lib_settings.sample_name)
        plot_index += 1
    if subtract_background:
        plt.suptitle('background-subtracted mutation rate fractions')
    if subtract_control:
        plt.suptitle('control-subtracted mutation rate fractions')
    else:
        plt.suptitle('mutation rate fractions')
    plt.savefig(out_prefix + '.pdf', transparent='True', format='pdf')
    plt.clf()

def plot_mutation_breakdown_pie(libraries, out_prefix, exclude_constitutive=False):
    #Makes an array of pie charts, 4 pers library, with types of mutations for each nt
    num_subplots = len(libraries)*4
    num_plots_wide = 4
    num_plots_high = len(libraries)
    fig = plt.figure(figsize=(16, 4*num_plots_high))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    plot_index =1
    for library in libraries:
        mutated_nts_count = library.count_mutation_types_by_nucleotide(exclude_constitutive=exclude_constitutive)
        for nt in 'ATCG':
            plot = fig.add_subplot(num_plots_high, num_plots_wide, plot_index)
            sorted_muts = sorted(mutated_nts_count[nt].items(), key=operator.itemgetter(1), reverse=True)
            labels = [pair[0] for pair in sorted_muts[:4]]
            others = [pair[0] for pair in sorted_muts[4:]]
            sizes = [mutated_nts_count[nt][mut_type] for mut_type in labels]
            labels.append('other')
            other_sum = sum([mutated_nts_count[nt][mut_type] for mut_type in others])
            sizes.append(other_sum)
            sizes = numpy.array(sizes)
            total = float(sum(sizes))
            sizes = sizes/total
            merged_labels = ['%s %.3f' % (labels[i], sizes[i]) for i in range(len(sizes))]
            plot.pie(sizes, labels = merged_labels, colors = mod_utils.rainbow)
            plot.set_title(library.lib_settings.sample_name)
            plot_index += 1
    plt.suptitle('mutation rate type fractions')
    plt.savefig(out_prefix + '.pdf', transparent='True', format='pdf')
    plt.clf()

def plot_mutation_rate_cdfs(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False):
    #Makes 2 CDF plots. One of all libraries, showing the coverage-normalized mutation rates
    # and one showing background-subtracted mutation rates

    fig = plt.figure(figsize=(24, 16))
    plots = []
    plot = fig.add_subplot(231)
    plots.append(plot)
    colormap = plt.get_cmap('spectral')
    colorindex = 0
    for library in libraries:
        all_mutation_rates = [val for val in
                              library.list_mutation_rates(subtract_background=False, subtract_control=False,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive)]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("mutation rate")
    plot.set_title('raw mutation rates')
    lg=plt.legend(loc=4,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-0.001, 0.02)

    plot = fig.add_subplot(232)
    plots.append(plot)
    colorindex = 0
    libraries_to_plot = [library for library in libraries if library.lib_settings.sample_name in
             library.experiment_settings.get_property('experimentals')]
    for library in libraries_to_plot:
        all_mutation_rates = [val for val in
                              library.list_mutation_rates(subtract_background=True, subtract_control=False,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive)]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("background-subtracted mutation rate")
    plot.set_title('normalized mutation rates')
    lg=plt.legend(loc=4,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-0.001, 0.02)

    plot = fig.add_subplot(233)
    plots.append(plot)
    colorindex = 0
    libraries_to_plot = [library for library in libraries if library.lib_settings.sample_name in
             library.experiment_settings.get_property('experimentals')]
    for library in libraries_to_plot:
        all_mutation_rates = [val for val in
                              library.list_mutation_rates(subtract_background=False, subtract_control=True,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive)]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("control-subtracted mutation rate")
    plot.set_title('control normalized mutation rates')
    lg=plt.legend(loc=4,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-0.001, 0.02)

    plot = fig.add_subplot(234)
    plots.append(plot)
    colormap = plt.get_cmap('spectral')
    colorindex = 0
    for library in libraries:
        all_mutation_rates = [math.log(val, 10) for val in
                              library.list_mutation_rates(subtract_background=False, subtract_control=False,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive) if val>0]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("log10 mutation rate")
    plot.set_title('raw mutation rates')
    lg=plt.legend(loc=2,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-5, -1)

    plot = fig.add_subplot(235)
    plots.append(plot)
    colorindex = 0
    libraries_to_plot = [library for library in libraries if library.lib_settings.sample_name in
             library.experiment_settings.get_property('experimentals')]
    for library in libraries_to_plot:
        all_mutation_rates = [math.log(val, 10) for val in
                              library.list_mutation_rates(subtract_background=True, subtract_control=False,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive) if val>0]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("background-subtracted log10 mutation rate")
    plot.set_title('normalized mutation rates')
    lg=plt.legend(loc=2,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-5, -1)

    plot = fig.add_subplot(236)
    plots.append(plot)
    colorindex = 0
    libraries_to_plot = [library for library in libraries if library.lib_settings.sample_name in
             library.experiment_settings.get_property('experimentals')]
    for library in libraries_to_plot:
        all_mutation_rates = [math.log(val, 10) for val in
                              library.list_mutation_rates(subtract_background=False, subtract_control=True,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive) if val>0]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("control-subtracted log10 mutation rate")
    plot.set_title('control normalized mutation rates')
    lg=plt.legend(loc=2,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-5, -1)

    for plot in plots:
        plot.set_ylabel("cumulative fraction of %s nucleotides" % (nucleotides_to_count))
        plot.set_ylim(0, 1)
    plt.tight_layout()
    plt.savefig(out_prefix + '.pdf', transparent='True', format='pdf')
    plt.clf()

def plot_mutation_rate_violins(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False):
    #Makes violin plots of raw mutation rates
    data = []
    labels = []
    for library in libraries:
        labels.append(library.lib_settings.sample_name)
        data.append([math.log10(val) for val in library.list_mutation_rates(subtract_background=False, subtract_control=False,
                                                        nucleotides_to_count=nucleotides_to_count,
                                                        exclude_constitutive=exclude_constitutive) if val>0])

    colormap = uniform_colormaps.viridis
    fig = plt.figure(figsize=(5,8))
    ax1 = fig.add_subplot(111)

    # Hide the grid behind plot objects
    ax1.yaxis.grid(True, linestyle='-', which='major', color='lightgrey', alpha=0.5)
    ax1.set_axisbelow(True)

    #ax1.set_xlabel(ylabel)
    plt.subplots_adjust(left=0.1, right=0.95, top=0.9, bottom=0.25)

    pos = range(1,len(libraries)+1)  # starts at 1 to play nice with boxplot
    dist = max(pos)-min(pos)
    w = min(0.15*max(dist,1.0),0.5)
    for library,p in zip(libraries,pos):
        d = [math.log10(val) for val in library.list_mutation_rates(subtract_background=False, subtract_control=False,
                                                        nucleotides_to_count=nucleotides_to_count,
                                                        exclude_constitutive=exclude_constitutive) if val>0]
        k = stats.gaussian_kde(d) #calculates the kernel density
        m = k.dataset.min() #lower bound of violin
        M = k.dataset.max() #upper bound of violin
        x = numpy.arange(m,M,(M-m)/100.) # support for violin
        v = k.evaluate(x) #violin profile (density curve)
        v = v/v.max()*w #scaling the violin to the available space
        plt.fill_betweenx(x,p,v+p,facecolor=colormap((p-1)/float(len(libraries))),alpha=0.3)
        plt.fill_betweenx(x,p,-v+p,facecolor=colormap((p-1)/float(len(libraries))),alpha=0.3)
    if True:
        bplot = plt.boxplot(data,notch=1)
        plt.setp(bplot['boxes'], color='black')
        plt.setp(bplot['whiskers'], color='black')
        plt.setp(bplot['fliers'], color='red', marker='.')

    per50s = []
    i = 1
    for datum in data:
        #per50s.append(stats.scoreatpercentile(datum, 50))
        t = stats.scoreatpercentile(datum, 50)

        per50s.append(t)
        #ax1.annotate(str(round(t,3)), xy=(i+0.1, t), xycoords='data', arrowprops=None, fontsize='small', color='black')
        i+= 1
    #ax1.set_xticks([0.0, 0.5, 1.0, 1.5])
    #ax1.set_yscale('log')
    ax1.set_ylabel('log10 mutation rate')
    ax1.set_ylim(-5, 0)
    xtickNames = plt.setp(ax1, xticklabels=labels)
    plt.setp(xtickNames, rotation=90, fontsize=6)
    plt.savefig(out_prefix+'_logviolin.pdf', transparent='True', format='pdf')
    plt.clf()

def ma_plots_interactive(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False,
                         max_fold_reduction=0.001, max_fold_increase=100):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :return: for each library use bokeh to plot an interactive plot of average magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    from bokeh.plotting import figure, output_file, show, save, ColumnDataSource, gridplot
    from bokeh.models import Range1d
    from bokeh.models import HoverTool
    from collections import OrderedDict

    # output to static HTML file
    output_file("%s.html" % (out_prefix))
    plot_figs=[]

    for library in libraries:
        mag, fold_change, annotation = [], [], []
        prot_mag, prot_fold_change, prot_annotation = [], [], []
        deprot_mag, deprot_fold_change, deprot_annotation = [], [], []
        for rRNA_name in library.rRNA_mutation_data:
            for position in library.rRNA_mutation_data[rRNA_name].nucleotides:
                nucleotide = library.rRNA_mutation_data[rRNA_name].nucleotides[position]
                if (exclude_constitutive and nucleotide.exclude_constitutive)or nucleotide.identity not in nucleotides_to_count:
                    continue
                else:
                    protection_call = nucleotide.determine_protection_status(confidence_interval=library.experiment_settings.get_property('confidence_interval_cutoff'),
                                                                   fold_change_cutoff=library.experiment_settings.get_property('fold_change_cutoff'))
                    control_fold_change = nucleotide.get_control_fold_change_in_mutation_rate()
                    avg_mutation_rate = (nucleotide.mutation_rate+nucleotide.get_control_nucleotide().mutation_rate)/2.0
                    if control_fold_change == 0:
                        control_fold_change = max_fold_reduction
                    elif control_fold_change == float('inf'):
                        control_fold_change = max_fold_increase
                    if protection_call == 'no_change':
                        mag.append(avg_mutation_rate)
                        fold_change.append(control_fold_change)
                        annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'deprotected':
                        deprot_mag.append(avg_mutation_rate)
                        deprot_fold_change.append(control_fold_change)
                        deprot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'protected':
                        prot_mag.append(avg_mutation_rate)
                        prot_fold_change.append(control_fold_change)
                        prot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
        source = ColumnDataSource(data=dict(x = mag, y = fold_change, label = annotation))
        prot_source = ColumnDataSource(data=dict(x = prot_mag, y = prot_fold_change, label = prot_annotation))
        deprot_source = ColumnDataSource(data=dict(x = deprot_mag, y = deprot_fold_change,
                                                   label = deprot_annotation))
        TOOLS = "pan,wheel_zoom,reset,save,hover"
        PlotFig = figure(x_axis_label = "avg signal ([%s] + [%s])/2" % (library.lib_settings.sample_name, library.get_normalizing_lib_with_mod().lib_settings.sample_name),
                         y_axis_label = "fold change [%s]/[%s]" % (library.lib_settings.sample_name, library.get_normalizing_lib_with_mod().lib_settings.sample_name),
                         y_axis_type="log", x_axis_type="log", tools=TOOLS, toolbar_location="right")
        PlotFig.circle("x", "y", size = 5, source=source, color=mod_utils.bokeh_black)
        PlotFig.circle("x", "y", size = 5, source=prot_source, color=mod_utils.bokeh_vermillion)
        PlotFig.circle("x", "y", size = 5, source=deprot_source, color=mod_utils.bokeh_bluishGreen)
        PlotFig.x_range = Range1d(start=0.00001, end=1)
        PlotFig.y_range = Range1d(start=.001, end=100)

        #adjust what information you get when you hover over it
        Hover = PlotFig.select(dict(type=HoverTool))
        Hover.tooltips = OrderedDict([("nuc", "@label")])
        plot_figs.append([PlotFig])
    p = gridplot(plot_figs)
    save(p)


def ma_plots_interactive_by_count(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False,
                         max_fold_reduction=0.001, max_fold_increase=100, lowess_correct=False):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :param max_fold_reduction:
    :param max_fold_increase:
    :param lowess_correct: wether to use lowess regression to correct fold_change data
    :return: for each library use bokeh to plot an interactive plot of average magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    from bokeh.plotting import figure, output_file, show, save, ColumnDataSource, gridplot
    from bokeh.models import Range1d
    from bokeh.models import HoverTool
    from collections import OrderedDict

    # output to static HTML file
    output_file("%s.html" % (out_prefix))
    plot_figs=[]

    for library in libraries:
        mag, fold_change, annotation = [], [], []
        prot_mag, prot_fold_change, prot_annotation = [], [], []
        deprot_mag, deprot_fold_change, deprot_annotation = [], [], []
        if lowess_correct:
            library.lowess_correct_fold_changes(nucleotides_to_count = nucleotides_to_count,
                                                exclude_constitutive=exclude_constitutive,
                                                max_fold_reduction=max_fold_reduction, max_fold_increase=max_fold_increase)
        for rRNA_name in library.rRNA_mutation_data:
            for position in library.rRNA_mutation_data[rRNA_name].nucleotides:
                nucleotide = library.rRNA_mutation_data[rRNA_name].nucleotides[position]
                if (exclude_constitutive and nucleotide.exclude_constitutive)or nucleotide.identity not in nucleotides_to_count:
                    continue
                else:
                    protection_call = nucleotide.determine_protection_status(confidence_interval=library.experiment_settings.get_property('confidence_interval_cutoff'),
                                                                   fold_change_cutoff=library.experiment_settings.get_property('fold_change_cutoff'), lowess_correct=lowess_correct)
                    if lowess_correct:
                        control_fold_change = nucleotide.lowess_fc
                    else:
                        control_fold_change = nucleotide.get_control_fold_change_in_mutation_rate()
                    avg_mutation_rate = (nucleotide.total_mutation_counts+nucleotide.get_control_nucleotide().total_mutation_counts)/2.0
                    if control_fold_change == 0:
                        control_fold_change = max_fold_reduction
                    elif control_fold_change == float('inf'):
                        control_fold_change = max_fold_increase
                    if protection_call == 'no_change':
                        mag.append(avg_mutation_rate)
                        fold_change.append(control_fold_change)
                        annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'deprotected':
                        deprot_mag.append(avg_mutation_rate)
                        deprot_fold_change.append(control_fold_change)
                        deprot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'protected':
                        prot_mag.append(avg_mutation_rate)
                        prot_fold_change.append(control_fold_change)
                        prot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
        source = ColumnDataSource(data=dict(x = mag, y = fold_change, label = annotation))
        prot_source = ColumnDataSource(data=dict(x = prot_mag, y = prot_fold_change, label = prot_annotation))
        deprot_source = ColumnDataSource(data=dict(x = deprot_mag, y = deprot_fold_change, label = deprot_annotation))
        fc_log = [math.log(fc, 10) for fc in fold_change]
        mag_log = [math.log(m, 10) if m>0 else -1. for m in mag]
        lowess_fc_log = lowess(fc_log, mag_log, return_sorted=False)
        lowess_fc = 10**lowess_fc_log
        lowess_fc_sort_by_mag = [x for (y,x) in sorted(zip(mag,lowess_fc), key=lambda pair: pair[0])]
        lowess_line = ColumnDataSource(data=dict(x = sorted(mag), y = lowess_fc_sort_by_mag))
        TOOLS = "pan,wheel_zoom,reset,save,hover"
        PlotFig = figure(x_axis_label = "average mutation counts",
                         y_axis_label = "fold change [%s]/[%s]" % (library.lib_settings.sample_name, library.get_normalizing_lib_with_mod().lib_settings.sample_name),
                         y_axis_type="log", x_axis_type="log", tools=TOOLS, toolbar_location="right")
        PlotFig.circle("x", "y", size = 5, source=source, color=mod_utils.bokeh_black)
        PlotFig.circle("x", "y", size = 5, source=prot_source, color=mod_utils.bokeh_vermillion)
        PlotFig.circle("x", "y", size = 5, source=deprot_source, color=mod_utils.bokeh_bluishGreen)
        PlotFig.line("x", "y", source=lowess_line, color='red', line_dash='dashed')
        PlotFig.x_range = Range1d(start=.1, end=100000)
        PlotFig.y_range = Range1d(start=.001, end=100)

        #adjust what information you get when you hover over it
        Hover = PlotFig.select(dict(type=HoverTool))
        Hover.tooltips = OrderedDict([("nuc", "@label")])
        plot_figs.append([PlotFig])
    p = gridplot(plot_figs)
    save(p)

def ma_plots(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False,
             max_fold_reduction=0.001, max_fold_increase=100):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :return: for each library use bokeh to plot an interactive plot of magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    output_file = "%s.pdf" % (out_prefix)
    plot_figs=[]

    num_subplots = len(libraries)
    num_plots_wide = math.ceil(math.sqrt(num_subplots))
    num_plots_high = num_plots_wide
    fig = plt.figure(figsize=(4*num_plots_wide, 4*num_plots_high))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    plot_index =1
    for library in libraries:
        plot = fig.add_subplot(num_plots_high, num_plots_wide, plot_index)
        mag, fold_change, annotation = [], [], []
        prot_mag, prot_fold_change, prot_annotation = [], [], []
        deprot_mag, deprot_fold_change, deprot_annotation = [], [], []
        for rRNA_name in library.rRNA_mutation_data:
            for position in library.rRNA_mutation_data[rRNA_name].nucleotides:
                nucleotide = library.rRNA_mutation_data[rRNA_name].nucleotides[position]
                if (exclude_constitutive and nucleotide.exclude_constitutive)or nucleotide.identity not in nucleotides_to_count:
                    continue
                else:
                    protection_call = nucleotide.determine_protection_status(confidence_interval=library.experiment_settings.get_property('confidence_interval_cutoff'),
                                                                   fold_change_cutoff=library.experiment_settings.get_property('fold_change_cutoff'))
                    control_fold_change = nucleotide.get_control_fold_change_in_mutation_rate()
                    avg_mutation_rate = (nucleotide.mutation_rate+nucleotide.get_control_nucleotide().mutation_rate)/2.0
                    if control_fold_change == 0:
                        control_fold_change = max_fold_reduction
                    elif control_fold_change == float('inf'):
                        control_fold_change = max_fold_increase
                    if protection_call == 'no_change':
                        mag.append(avg_mutation_rate)
                        fold_change.append(control_fold_change)
                        annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'deprotected':
                        deprot_mag.append(avg_mutation_rate)
                        deprot_fold_change.append(control_fold_change)
                        deprot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'protected':
                        prot_mag.append(avg_mutation_rate)
                        prot_fold_change.append(control_fold_change)
                        prot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
        plot.set_xlabel("([%s] + [%s])/2" % (library.lib_settings.sample_name, library.get_normalizing_lib_with_mod().lib_settings.sample_name), fontsize = 8)
        plot.set_ylabel("[%s]/[%s]" % (library.lib_settings.sample_name, library.get_normalizing_lib_with_mod().lib_settings.sample_name), fontsize = 8)
        plot.set_yscale('log')
        plot.set_xscale('log')
        plot.scatter(mag, fold_change, color=mod_utils.black, s=3)
        plot.scatter(prot_mag, prot_fold_change, color=mod_utils.vermillion, s=5)
        plot.scatter(deprot_mag, deprot_fold_change, color=mod_utils.bluishGreen, s=5)
        plot.set_xlim(0.00001,1)
        plot.set_ylim(.001,100)
        plot_figs.append(plot)
        plot_index+=1
    plt.savefig(output_file, transparent='True', format='pdf')

def mutation_rate_scatter(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :return: for each library use bokeh to plot an interactive plot of magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    output_file = "%s.pdf" % (out_prefix)
    plot_figs=[]

    num_subplots = len(libraries)
    num_plots_wide = math.ceil(math.sqrt(num_subplots))
    num_plots_high = num_plots_wide
    fig = plt.figure(figsize=(4*num_plots_wide, 4*num_plots_high))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    plot_index =1
    for library in libraries:
        plot = fig.add_subplot(num_plots_high, num_plots_wide, plot_index)
        x, y, annotation = [], [], []
        for rRNA_name in library.rRNA_mutation_data:
            for position in library.rRNA_mutation_data[rRNA_name].nucleotides:
                nucleotide = library.rRNA_mutation_data[rRNA_name].nucleotides[position]
                if (exclude_constitutive and nucleotide.exclude_constitutive)or nucleotide.identity not in nucleotides_to_count:
                    continue
                else:
                    x.append(nucleotide.mutation_rate)
                    y.append(nucleotide.get_control_nucleotide().mutation_rate)
        plot.set_xlabel("%s mismatch_fraction" % (library.lib_settings.sample_name), fontsize = 8)
        plot.set_ylabel("%s mismatch_fraction" % (library.get_normalizing_lib_with_mod().lib_settings.sample_name), fontsize = 8)
        plot.set_yscale('log')
        plot.set_xscale('log')
        x=numpy.array(x)
        y=numpy.array(y)
        plot.scatter(x, y, color=mod_utils.black, s=3)
        logx=numpy.log10(x)
        logy=numpy.log10(y)
        logx, logy = mod_utils.filter_x_y_pairs(logx, logy)
        r, p = stats.pearsonr(logx, logy)
        plot.annotate("r^2 = {:.2f}".format(r**2),
                    xy=(.1, .9), xycoords=plot.transAxes)
        plot.set_xlim(0.00001,1)
        plot.set_ylim(0.00001,1)
        plot_figs.append(plot)
        plot_index+=1
    plt.savefig(output_file, transparent='True', format='pdf')

def ma_plots_by_count(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False,
             max_fold_reduction=0.001, max_fold_increase=100, lowess_correct=False):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :return: for each library use bokeh to plot an interactive plot of magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    output_file = "%s.pdf" % (out_prefix)
    plot_figs=[]

    num_subplots = len(libraries)
    num_plots_wide = math.ceil(math.sqrt(num_subplots))
    num_plots_high = num_plots_wide
    fig = plt.figure(figsize=(4*num_plots_wide, 4*num_plots_high))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    plot_index =1
    for library in libraries:
        if lowess_correct:
            library.lowess_correct_fold_changes(nucleotides_to_count = nucleotides_to_count,
                                                exclude_constitutive=exclude_constitutive,
                                                max_fold_reduction=max_fold_reduction, max_fold_increase=max_fold_increase)
        plot = fig.add_subplot(num_plots_high, num_plots_wide, plot_index)<commit_msg>changed some labels. added interactive scatterplot<commit_after>__author__ = 'boris'
from scipy.stats.mstats import winsorize
import scipy.stats as stats
import matplotlib.pyplot as plt
from collections import defaultdict
import numpy
import math
import mod_lib
import mod_utils
import operator
import os
import uniform_colormaps
from statsmodels.nonparametric.smoothers_lowess import lowess
from matplotlib import rc
plt.rcParams['pdf.fonttype'] = 42 #leaves most text as actual text in PDFs, not outlines
#plt.rcParams["font.family"] = "sans-serif"
rc('font',**{'family':'sans-serif','sans-serif':['Bitstream Vera Sans']})
#rc('text', usetex=True)

def plot_mutated_nts_pie(libraries, out_prefix, subtract_background=False, subtract_control=False, exclude_constitutive=False):
    #Makes an array of pie charts, 1 per library
    if subtract_background:
        #if subtracting background, need to only look at those which have a defined control
        libraries = [library for library in libraries if (library.lib_settings.sample_name in
                     library.experiment_settings.get_property('experimentals')) or (library.lib_settings.sample_name in
                     library.experiment_settings.get_property('with_mod_controls'))]
    elif subtract_control:
        #if subtracting background, need to only look at those which have a defined control
        libraries = [library for library in libraries if library.lib_settings.sample_name in
                     library.experiment_settings.get_property('experimentals')]
    num_subplots = len(libraries)
    num_plots_wide = math.ceil(math.sqrt(num_subplots))
    num_plots_high = num_plots_wide
    fig = plt.figure(figsize=(4*num_plots_wide, 4*num_plots_high))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    plot_index =1
    for library in libraries:
        plot = fig.add_subplot(num_plots_high, num_plots_wide, plot_index)
        mutated_nts_count = library.count_mutation_rates_by_nucleotide(subtract_background=subtract_background, subtract_control=subtract_control,
                                                                       exclude_constitutive=exclude_constitutive)
        labels = sorted(mutated_nts_count.keys())
        sizes = numpy.array([mutated_nts_count[nt] for nt in labels])
        total = float(sum(sizes))
        sizes = sizes/total
        merged_labels = ['%s %.3f' % (labels[i], sizes[i]) for i in range(len(sizes))]
        plot.pie(sizes, labels = merged_labels, colors = mod_utils.rainbow)
        plot.set_title(library.lib_settings.sample_name)
        plot_index += 1
    if subtract_background:
        plt.suptitle('background-subtracted mutation rate fractions')
    if subtract_control:
        plt.suptitle('control-subtracted mutation rate fractions')
    else:
        plt.suptitle('mutation rate fractions')
    plt.savefig(out_prefix + '.pdf', transparent='True', format='pdf')
    plt.clf()

def plot_mutation_breakdown_pie(libraries, out_prefix, exclude_constitutive=False):
    #Makes an array of pie charts, 4 pers library, with types of mutations for each nt
    num_subplots = len(libraries)*4
    num_plots_wide = 4
    num_plots_high = len(libraries)
    fig = plt.figure(figsize=(16, 4*num_plots_high))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    plot_index =1
    for library in libraries:
        mutated_nts_count = library.count_mutation_types_by_nucleotide(exclude_constitutive=exclude_constitutive)
        for nt in 'ATCG':
            plot = fig.add_subplot(num_plots_high, num_plots_wide, plot_index)
            sorted_muts = sorted(mutated_nts_count[nt].items(), key=operator.itemgetter(1), reverse=True)
            labels = [pair[0] for pair in sorted_muts[:4]]
            others = [pair[0] for pair in sorted_muts[4:]]
            sizes = [mutated_nts_count[nt][mut_type] for mut_type in labels]
            labels.append('other')
            other_sum = sum([mutated_nts_count[nt][mut_type] for mut_type in others])
            sizes.append(other_sum)
            sizes = numpy.array(sizes)
            total = float(sum(sizes))
            sizes = sizes/total
            merged_labels = ['%s %.3f' % (labels[i], sizes[i]) for i in range(len(sizes))]
            plot.pie(sizes, labels = merged_labels, colors = mod_utils.rainbow)
            plot.set_title(library.lib_settings.sample_name)
            plot_index += 1
    plt.suptitle('mutation rate type fractions')
    plt.savefig(out_prefix + '.pdf', transparent='True', format='pdf')
    plt.clf()

def plot_mutation_rate_cdfs(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False):
    #Makes 2 CDF plots. One of all libraries, showing the coverage-normalized mutation rates
    # and one showing background-subtracted mutation rates

    fig = plt.figure(figsize=(24, 16))
    plots = []
    plot = fig.add_subplot(231)
    plots.append(plot)
    colormap = plt.get_cmap('spectral')
    colorindex = 0
    for library in libraries:
        all_mutation_rates = [val for val in
                              library.list_mutation_rates(subtract_background=False, subtract_control=False,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive)]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("mutation rate")
    plot.set_title('raw mutation rates')
    lg=plt.legend(loc=4,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-0.001, 0.02)

    plot = fig.add_subplot(232)
    plots.append(plot)
    colorindex = 0
    libraries_to_plot = [library for library in libraries if library.lib_settings.sample_name in
             library.experiment_settings.get_property('experimentals')]
    for library in libraries_to_plot:
        all_mutation_rates = [val for val in
                              library.list_mutation_rates(subtract_background=True, subtract_control=False,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive)]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("background-subtracted mutation rate")
    plot.set_title('normalized mutation rates')
    lg=plt.legend(loc=4,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-0.001, 0.02)

    plot = fig.add_subplot(233)
    plots.append(plot)
    colorindex = 0
    libraries_to_plot = [library for library in libraries if library.lib_settings.sample_name in
             library.experiment_settings.get_property('experimentals')]
    for library in libraries_to_plot:
        all_mutation_rates = [val for val in
                              library.list_mutation_rates(subtract_background=False, subtract_control=True,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive)]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("control-subtracted mutation rate")
    plot.set_title('control normalized mutation rates')
    lg=plt.legend(loc=4,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-0.001, 0.02)

    plot = fig.add_subplot(234)
    plots.append(plot)
    colormap = plt.get_cmap('spectral')
    colorindex = 0
    for library in libraries:
        all_mutation_rates = [math.log(val, 10) for val in
                              library.list_mutation_rates(subtract_background=False, subtract_control=False,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive) if val>0]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("log10 mutation rate")
    plot.set_title('raw mutation rates')
    lg=plt.legend(loc=2,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-5, -1)

    plot = fig.add_subplot(235)
    plots.append(plot)
    colorindex = 0
    libraries_to_plot = [library for library in libraries if library.lib_settings.sample_name in
             library.experiment_settings.get_property('experimentals')]
    for library in libraries_to_plot:
        all_mutation_rates = [math.log(val, 10) for val in
                              library.list_mutation_rates(subtract_background=True, subtract_control=False,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive) if val>0]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("background-subtracted log10 mutation rate")
    plot.set_title('normalized mutation rates')
    lg=plt.legend(loc=2,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-5, -1)

    plot = fig.add_subplot(236)
    plots.append(plot)
    colorindex = 0
    libraries_to_plot = [library for library in libraries if library.lib_settings.sample_name in
             library.experiment_settings.get_property('experimentals')]
    for library in libraries_to_plot:
        all_mutation_rates = [math.log(val, 10) for val in
                              library.list_mutation_rates(subtract_background=False, subtract_control=True,
                                                          nucleotides_to_count=nucleotides_to_count, exclude_constitutive=exclude_constitutive) if val>0]
        plot.hist(all_mutation_rates, 10000, normed=1, cumulative=True, histtype='step', color=colormap(colorindex/float(len(libraries))),
                  label=library.lib_settings.sample_name, lw=2)
        colorindex += 1
    plot.set_xlabel("control-subtracted log10 mutation rate")
    plot.set_title('control normalized mutation rates')
    lg=plt.legend(loc=2,prop={'size':6}, labelspacing=0.2)
    lg.draw_frame(False)
    plot.set_xlim(-5, -1)

    for plot in plots:
        plot.set_ylabel("cumulative fraction of %s nucleotides" % (nucleotides_to_count))
        plot.set_ylim(0, 1)
    plt.tight_layout()
    plt.savefig(out_prefix + '.pdf', transparent='True', format='pdf')
    plt.clf()

def plot_mutation_rate_violins(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False):
    #Makes violin plots of raw mutation rates
    data = []
    labels = []
    for library in libraries:
        labels.append(library.lib_settings.sample_name)
        data.append([math.log10(val) for val in library.list_mutation_rates(subtract_background=False, subtract_control=False,
                                                        nucleotides_to_count=nucleotides_to_count,
                                                        exclude_constitutive=exclude_constitutive) if val>0])

    colormap = uniform_colormaps.viridis
    fig = plt.figure(figsize=(5,8))
    ax1 = fig.add_subplot(111)

    # Hide the grid behind plot objects
    ax1.yaxis.grid(True, linestyle='-', which='major', color='lightgrey', alpha=0.5)
    ax1.set_axisbelow(True)

    #ax1.set_xlabel(ylabel)
    plt.subplots_adjust(left=0.1, right=0.95, top=0.9, bottom=0.25)

    pos = range(1,len(libraries)+1)  # starts at 1 to play nice with boxplot
    dist = max(pos)-min(pos)
    w = min(0.15*max(dist,1.0),0.5)
    for library,p in zip(libraries,pos):
        d = [math.log10(val) for val in library.list_mutation_rates(subtract_background=False, subtract_control=False,
                                                        nucleotides_to_count=nucleotides_to_count,
                                                        exclude_constitutive=exclude_constitutive) if val>0]
        k = stats.gaussian_kde(d) #calculates the kernel density
        m = k.dataset.min() #lower bound of violin
        M = k.dataset.max() #upper bound of violin
        x = numpy.arange(m,M,(M-m)/100.) # support for violin
        v = k.evaluate(x) #violin profile (density curve)
        v = v/v.max()*w #scaling the violin to the available space
        plt.fill_betweenx(x,p,v+p,facecolor=colormap((p-1)/float(len(libraries))),alpha=0.3)
        plt.fill_betweenx(x,p,-v+p,facecolor=colormap((p-1)/float(len(libraries))),alpha=0.3)
    if True:
        bplot = plt.boxplot(data,notch=1)
        plt.setp(bplot['boxes'], color='black')
        plt.setp(bplot['whiskers'], color='black')
        plt.setp(bplot['fliers'], color='red', marker='.')

    per50s = []
    i = 1
    for datum in data:
        #per50s.append(stats.scoreatpercentile(datum, 50))
        t = stats.scoreatpercentile(datum, 50)

        per50s.append(t)
        #ax1.annotate(str(round(t,3)), xy=(i+0.1, t), xycoords='data', arrowprops=None, fontsize='small', color='black')
        i+= 1
    #ax1.set_xticks([0.0, 0.5, 1.0, 1.5])
    #ax1.set_yscale('log')
    ax1.set_ylabel('log10 mutation rate')
    ax1.set_ylim(-5, 0)
    xtickNames = plt.setp(ax1, xticklabels=labels)
    plt.setp(xtickNames, rotation=90, fontsize=6)
    plt.savefig(out_prefix+'_logviolin.pdf', transparent='True', format='pdf')
    plt.clf()

def ma_plots_interactive(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False,
                         max_fold_reduction=0.001, max_fold_increase=100):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :return: for each library use bokeh to plot an interactive plot of average magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    from bokeh.plotting import figure, output_file, show, save, ColumnDataSource, gridplot
    from bokeh.models import Range1d
    from bokeh.models import HoverTool
    from collections import OrderedDict

    # output to static HTML file
    output_file("%s.html" % (out_prefix))
    plot_figs=[]

    for library in libraries:
        mag, fold_change, annotation = [], [], []
        prot_mag, prot_fold_change, prot_annotation = [], [], []
        deprot_mag, deprot_fold_change, deprot_annotation = [], [], []
        for rRNA_name in library.rRNA_mutation_data:
            for position in library.rRNA_mutation_data[rRNA_name].nucleotides:
                nucleotide = library.rRNA_mutation_data[rRNA_name].nucleotides[position]
                if (exclude_constitutive and nucleotide.exclude_constitutive)or nucleotide.identity not in nucleotides_to_count:
                    continue
                else:
                    protection_call = nucleotide.determine_protection_status(confidence_interval=library.experiment_settings.get_property('confidence_interval_cutoff'),
                                                                   fold_change_cutoff=library.experiment_settings.get_property('fold_change_cutoff'))
                    control_fold_change = nucleotide.get_control_fold_change_in_mutation_rate()
                    avg_mutation_rate = (nucleotide.mutation_rate+nucleotide.get_control_nucleotide().mutation_rate)/2.0
                    if control_fold_change == 0:
                        control_fold_change = max_fold_reduction
                    elif control_fold_change == float('inf'):
                        control_fold_change = max_fold_increase
                    if protection_call == 'no_change':
                        mag.append(avg_mutation_rate)
                        fold_change.append(control_fold_change)
                        annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'deprotected':
                        deprot_mag.append(avg_mutation_rate)
                        deprot_fold_change.append(control_fold_change)
                        deprot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'protected':
                        prot_mag.append(avg_mutation_rate)
                        prot_fold_change.append(control_fold_change)
                        prot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
        source = ColumnDataSource(data=dict(x = mag, y = fold_change, label = annotation))
        prot_source = ColumnDataSource(data=dict(x = prot_mag, y = prot_fold_change, label = prot_annotation))
        deprot_source = ColumnDataSource(data=dict(x = deprot_mag, y = deprot_fold_change,
                                                   label = deprot_annotation))
        TOOLS = "pan,wheel_zoom,reset,save,hover"
        PlotFig = figure(x_axis_label = "avg signal ([%s] + [%s])/2" % (library.lib_settings.sample_name, library.get_normalizing_lib_with_mod().lib_settings.sample_name),
                         y_axis_label = "fold change [%s]/[%s]" % (library.lib_settings.sample_name, library.get_normalizing_lib_with_mod().lib_settings.sample_name),
                         y_axis_type="log", x_axis_type="log", tools=TOOLS, toolbar_location="right")
        PlotFig.circle("x", "y", size = 5, source=source, color=mod_utils.bokeh_black)
        PlotFig.circle("x", "y", size = 5, source=prot_source, color=mod_utils.bokeh_vermillion)
        PlotFig.circle("x", "y", size = 5, source=deprot_source, color=mod_utils.bokeh_bluishGreen)
        PlotFig.x_range = Range1d(start=0.00001, end=1)
        PlotFig.y_range = Range1d(start=.001, end=100)

        #adjust what information you get when you hover over it
        Hover = PlotFig.select(dict(type=HoverTool))
        Hover.tooltips = OrderedDict([("nuc", "@label")])
        plot_figs.append([PlotFig])
    p = gridplot(plot_figs)
    save(p)


def ma_plots_interactive_by_count(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False,
                         max_fold_reduction=0.001, max_fold_increase=100, lowess_correct=False):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :param max_fold_reduction:
    :param max_fold_increase:
    :param lowess_correct: wether to use lowess regression to correct fold_change data
    :return: for each library use bokeh to plot an interactive plot of average magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    from bokeh.plotting import figure, output_file, show, save, ColumnDataSource, gridplot
    from bokeh.models import Range1d
    from bokeh.models import HoverTool
    from collections import OrderedDict

    # output to static HTML file
    output_file("%s.html" % (out_prefix))
    plot_figs=[]

    for library in libraries:
        mag, fold_change, annotation = [], [], []
        prot_mag, prot_fold_change, prot_annotation = [], [], []
        deprot_mag, deprot_fold_change, deprot_annotation = [], [], []
        if lowess_correct:
            library.lowess_correct_fold_changes(nucleotides_to_count = nucleotides_to_count,
                                                exclude_constitutive=exclude_constitutive,
                                                max_fold_reduction=max_fold_reduction, max_fold_increase=max_fold_increase)
        for rRNA_name in library.rRNA_mutation_data:
            for position in library.rRNA_mutation_data[rRNA_name].nucleotides:
                nucleotide = library.rRNA_mutation_data[rRNA_name].nucleotides[position]
                if (exclude_constitutive and nucleotide.exclude_constitutive)or nucleotide.identity not in nucleotides_to_count:
                    continue
                else:
                    protection_call = nucleotide.determine_protection_status(confidence_interval=library.experiment_settings.get_property('confidence_interval_cutoff'),
                                                                   fold_change_cutoff=library.experiment_settings.get_property('fold_change_cutoff'), lowess_correct=lowess_correct)
                    if lowess_correct:
                        control_fold_change = nucleotide.lowess_fc
                    else:
                        control_fold_change = nucleotide.get_control_fold_change_in_mutation_rate()
                    avg_mutation_rate = (nucleotide.total_mutation_counts+nucleotide.get_control_nucleotide().total_mutation_counts)/2.0
                    if control_fold_change == 0:
                        control_fold_change = max_fold_reduction
                    elif control_fold_change == float('inf'):
                        control_fold_change = max_fold_increase
                    if protection_call == 'no_change':
                        mag.append(avg_mutation_rate)
                        fold_change.append(control_fold_change)
                        annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'deprotected':
                        deprot_mag.append(avg_mutation_rate)
                        deprot_fold_change.append(control_fold_change)
                        deprot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'protected':
                        prot_mag.append(avg_mutation_rate)
                        prot_fold_change.append(control_fold_change)
                        prot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
        source = ColumnDataSource(data=dict(x = mag, y = fold_change, label = annotation))
        prot_source = ColumnDataSource(data=dict(x = prot_mag, y = prot_fold_change, label = prot_annotation))
        deprot_source = ColumnDataSource(data=dict(x = deprot_mag, y = deprot_fold_change, label = deprot_annotation))
        fc_log = [math.log(fc, 10) for fc in fold_change]
        mag_log = [math.log(m, 10) if m>0 else -1. for m in mag]
        lowess_fc_log = lowess(fc_log, mag_log, return_sorted=False)
        lowess_fc = 10**lowess_fc_log
        lowess_fc_sort_by_mag = [x for (y,x) in sorted(zip(mag,lowess_fc), key=lambda pair: pair[0])]
        lowess_line = ColumnDataSource(data=dict(x = sorted(mag), y = lowess_fc_sort_by_mag))
        TOOLS = "pan,wheel_zoom,reset,save,hover"
        PlotFig = figure(x_axis_label = "average mutation counts",
                         y_axis_label = "fold change [%s]/[%s]" % (library.lib_settings.sample_name, library.get_normalizing_lib_with_mod().lib_settings.sample_name),
                         y_axis_type="log", x_axis_type="log", tools=TOOLS, toolbar_location="right")
        PlotFig.circle("x", "y", size = 5, source=source, color=mod_utils.bokeh_black)
        PlotFig.circle("x", "y", size = 5, source=prot_source, color=mod_utils.bokeh_vermillion)
        PlotFig.circle("x", "y", size = 5, source=deprot_source, color=mod_utils.bokeh_bluishGreen)
        PlotFig.line("x", "y", source=lowess_line, color='red', line_dash='dashed')
        PlotFig.x_range = Range1d(start=.1, end=100000)
        PlotFig.y_range = Range1d(start=.001, end=100)

        #adjust what information you get when you hover over it
        Hover = PlotFig.select(dict(type=HoverTool))
        Hover.tooltips = OrderedDict([("nuc", "@label")])
        plot_figs.append([PlotFig])
    p = gridplot(plot_figs)
    save(p)

def ma_plots(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False,
             max_fold_reduction=0.001, max_fold_increase=100):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :return: for each library use bokeh to plot an interactive plot of magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    output_file = "%s.pdf" % (out_prefix)
    plot_figs=[]

    num_subplots = len(libraries)
    num_plots_wide = math.ceil(math.sqrt(num_subplots))
    num_plots_high = num_plots_wide
    fig = plt.figure(figsize=(4*num_plots_wide, 4*num_plots_high))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    plot_index =1
    for library in libraries:
        plot = fig.add_subplot(num_plots_high, num_plots_wide, plot_index)
        mag, fold_change, annotation = [], [], []
        prot_mag, prot_fold_change, prot_annotation = [], [], []
        deprot_mag, deprot_fold_change, deprot_annotation = [], [], []
        for rRNA_name in library.rRNA_mutation_data:
            for position in library.rRNA_mutation_data[rRNA_name].nucleotides:
                nucleotide = library.rRNA_mutation_data[rRNA_name].nucleotides[position]
                if (exclude_constitutive and nucleotide.exclude_constitutive)or nucleotide.identity not in nucleotides_to_count:
                    continue
                else:
                    protection_call = nucleotide.determine_protection_status(confidence_interval=library.experiment_settings.get_property('confidence_interval_cutoff'),
                                                                   fold_change_cutoff=library.experiment_settings.get_property('fold_change_cutoff'))
                    control_fold_change = nucleotide.get_control_fold_change_in_mutation_rate()
                    avg_mutation_rate = (nucleotide.mutation_rate+nucleotide.get_control_nucleotide().mutation_rate)/2.0
                    if control_fold_change == 0:
                        control_fold_change = max_fold_reduction
                    elif control_fold_change == float('inf'):
                        control_fold_change = max_fold_increase
                    if protection_call == 'no_change':
                        mag.append(avg_mutation_rate)
                        fold_change.append(control_fold_change)
                        annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'deprotected':
                        deprot_mag.append(avg_mutation_rate)
                        deprot_fold_change.append(control_fold_change)
                        deprot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'protected':
                        prot_mag.append(avg_mutation_rate)
                        prot_fold_change.append(control_fold_change)
                        prot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
        plot.set_xlabel("([%s] + [%s])/2" % (library.lib_settings.sample_name, library.get_normalizing_lib_with_mod().lib_settings.sample_name), fontsize = 8)
        plot.set_ylabel("[%s]/[%s]" % (library.lib_settings.sample_name, library.get_normalizing_lib_with_mod().lib_settings.sample_name), fontsize = 8)
        plot.set_yscale('log')
        plot.set_xscale('log')
        plot.scatter(mag, fold_change, color=mod_utils.black, s=3)
        plot.scatter(prot_mag, prot_fold_change, color=mod_utils.vermillion, s=5)
        plot.scatter(deprot_mag, deprot_fold_change, color=mod_utils.bluishGreen, s=5)
        plot.set_xlim(0.00001,1)
        plot.set_ylim(.001,100)
        plot_figs.append(plot)
        plot_index+=1
    plt.savefig(output_file, transparent='True', format='pdf')

def mutation_rate_scatter(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :return: for each library use bokeh to plot an interactive plot of magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    output_file = "%s.pdf" % (out_prefix)
    plot_figs=[]

    num_subplots = len(libraries)
    num_plots_wide = math.ceil(math.sqrt(num_subplots))
    num_plots_high = num_plots_wide
    fig = plt.figure(figsize=(4*num_plots_wide, 4*num_plots_high))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    plot_index =1
    for library in libraries:
        plot = fig.add_subplot(num_plots_high, num_plots_wide, plot_index)
        x, y, annotation = [], [], []
        for rRNA_name in library.rRNA_mutation_data:
            for position in library.rRNA_mutation_data[rRNA_name].nucleotides:
                nucleotide = library.rRNA_mutation_data[rRNA_name].nucleotides[position]
                if (exclude_constitutive and nucleotide.exclude_constitutive)or nucleotide.identity not in nucleotides_to_count:
                    continue
                else:
                    x.append(nucleotide.mutation_rate)
                    y.append(nucleotide.get_control_nucleotide().mutation_rate)
        plot.set_xlabel("%s misincorporation fraction" % (library.lib_settings.sample_name), fontsize = 8)
        plot.set_ylabel("%s misincorporation fraction" % (library.get_normalizing_lib_with_mod().lib_settings.sample_name), fontsize = 8)
        plot.set_yscale('log')
        plot.set_xscale('log')
        x=numpy.array(x)
        y=numpy.array(y)
        plot.scatter(x, y, color=mod_utils.black, s=3)
        logx=numpy.log10(x)
        logy=numpy.log10(y)
        logx, logy = mod_utils.filter_x_y_pairs(logx, logy)
        r, p = stats.pearsonr(logx, logy)
        plot.annotate("r^2 = {:.2f}".format(r**2),
                    xy=(.1, .9), xycoords=plot.transAxes)
        plot.set_xlim(0.00001,1)
        plot.set_ylim(0.00001,1)
        plot_figs.append(plot)
        plot_index+=1
    plt.savefig(output_file, transparent='True', format='pdf')

def scatter_interactive(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False,
                         max_fold_reduction=0.001, max_fold_increase=100):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :return: for each library use bokeh to plot an interactive plot of average magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    from bokeh.plotting import figure, output_file, show, save, ColumnDataSource, gridplot
    from bokeh.models import Range1d
    from bokeh.models import HoverTool
    from collections import OrderedDict

    # output to static HTML file
    output_file("%s.html" % (out_prefix))
    plot_figs=[]

    for library in libraries:
        x, y, annotation = [], [], []
        prot_x, prot_y, prot_annotation = [], [], []
        deprot_x, deprot_y, deprot_annotation = [], [], []
        for rRNA_name in library.rRNA_mutation_data:
            for position in library.rRNA_mutation_data[rRNA_name].nucleotides:
                nucleotide = library.rRNA_mutation_data[rRNA_name].nucleotides[position]
                if (exclude_constitutive and nucleotide.exclude_constitutive)or nucleotide.identity not in nucleotides_to_count:
                    continue
                else:
                    protection_call = nucleotide.determine_protection_status(confidence_interval=library.experiment_settings.get_property('confidence_interval_cutoff'),
                                                                   fold_change_cutoff=library.experiment_settings.get_property('fold_change_cutoff'))
                    control_fold_change = nucleotide.get_control_fold_change_in_mutation_rate()
                    avg_mutation_rate = (nucleotide.mutation_rate+nucleotide.get_control_nucleotide().mutation_rate)/2.0
                    if control_fold_change == 0:
                        control_fold_change = max_fold_reduction
                    elif control_fold_change == float('inf'):
                        control_fold_change = max_fold_increase
                    if protection_call == 'no_change':
                        x.append(nucleotide.get_control_nucleotide().mutation_rate)
                        y.append(nucleotide.mutation_rate)
                        annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'deprotected':
                        deprot_x.append(nucleotide.get_control_nucleotide().mutation_rate)
                        deprot_y.append(nucleotide.mutation_rate)
                        deprot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
                    elif protection_call == 'protected':
                        prot_x.append(nucleotide.get_control_nucleotide().mutation_rate)
                        prot_y.append(nucleotide.mutation_rate)
                        prot_annotation.append('%s_%s%d' %(rRNA_name,nucleotide.identity,position))
        source = ColumnDataSource(data=dict(x = x, y = y, label = annotation))
        prot_source = ColumnDataSource(data=dict(x = prot_x, y = prot_y, label = prot_annotation))
        deprot_source = ColumnDataSource(data=dict(x = deprot_x, y = deprot_y,
                                                   label = deprot_annotation))
        TOOLS = "pan,wheel_zoom,reset,save,hover"
        PlotFig = figure(x_axis_label = "%s misincorporation fraction" % (library.get_normalizing_lib_with_mod().lib_settings.sample_name),
                         y_axis_label = "%s misincorporation fraction" % (library.lib_settings.sample_name),
                         y_axis_type="log", x_axis_type="log", tools=TOOLS, toolbar_location="right")
        PlotFig.circle("x", "y", size = 5, source=source, color=mod_utils.bokeh_black)
        PlotFig.circle("x", "y", size = 5, source=prot_source, color=mod_utils.bokeh_vermillion)
        PlotFig.circle("x", "y", size = 5, source=deprot_source, color=mod_utils.bokeh_bluishGreen)
        PlotFig.x_range = Range1d(start=0.00001, end=1)
        PlotFig.y_range = Range1d(start=0.00001, end=1)

        #adjust what information you get when you hover over it
        Hover = PlotFig.select(dict(type=HoverTool))
        Hover.tooltips = OrderedDict([("nuc", "@label")])
        plot_figs.append([PlotFig])
    p = gridplot(plot_figs)
    save(p)


def ma_plots_by_count(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False,
             max_fold_reduction=0.001, max_fold_increase=100, lowess_correct=False):
    """

    :param libraries:
    :param out_prefix:
    :param nucleotides_to_count:
    :param exclude_constitutive:
    :return: for each library use bokeh to plot an interactive plot of magnitude of signal (experimental+control)/2
            vs log10 fold change (experimental/control).
            Protected and de-protected calls will be colored, based on a fold change cutoff and confidence interval.
            All nucleotides will be labelled on mouseover.
    """
    output_file = "%s.pdf" % (out_prefix)
    plot_figs=[]

    num_subplots = len(libraries)
    num_plots_wide = math.ceil(math.sqrt(num_subplots))
    num_plots_high = num_plots_wide
    fig = plt.figure(figsize=(4*num_plots_wide, 4*num_plots_high))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    plot_index =1
    for library in libraries:
        if lowess_correct:
            library.lowess_correct_fold_changes(nucleotides_to_count = nucleotides_to_count,
                                                exclude_constitutive=exclude_constitutive,
                                                max_fold_reduction=max_fold_reduction, max_fold_increase=max_fold_increase)
        plot = fig.add_subplot(num_plots_high, num_plots_wide, plot_index)
{'fraction_non_alphanumeric': 0.09734487692774654, 'fraction_numerical': 0.012791420353534623, 'mean_word_length': 2.8170583691934237, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 82, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '755597', 'n_tokens_mistral': 22231, 'n_tokens_neox': 20364, 'n_words': 3853}