Create sequence logo showing amino-acid or nucleotide preferences.

    The heights of each letter is equal to the preference of
    that site for that amino acid or nucleotide.

    Note that stop codons may or may not be included in the logo
    depending on whether they are present in *pi_d*.  

    CALLING VARIABLES:

    * *sites* is a list of all of the sites that are being included
      in the logo, as strings. They must be in natural sort or an error
      will be raised **unless** *allowunsorted* is *True*. The sites
      in the plot are ordered in the same arrangement
      listed in *sites*. These should be **strings**, not integers.

    * *datatype* should be one of the following strings:
    
        * 'prefs' for preferences
        
        * 'diffprefs' for differential preferences
        
        * 'diffsel' for differential selection

    * *data* is a dictionary that has a key for every entry in
      *sites*. For every site *r* in *sites*, *sites[r][x]*
      is the value for character *x*. 
      Preferences must sum to one; differential preferences to zero.
      All sites must have the same set of characters. The characters
      must be the set of nucleotides or amino acids with or without
      stop codons.

    * *plotfile* is a string giving the name of the created PDF file 
      of the logo plot.
      It must end in the extension ``.pdf``.

    * *nperline* is the number of sites per line. Often 40 to 80 are good values.

    * *numberevery* is specifies how frequently we put labels for the sites on
      x-axis.

    * *allowunsorted* : if *True* then we allow the entries in *sites* to 
      **not** be sorted. This means that the logo plot will **not** have
      sites in sorted order.

    * *ydatamax* : meaningful only if *datatype* is 'diffprefs'. In this case, it gives
      the maximum that the logo stacks extend in the positive and negative directions.
      Cannot be smaller than the maximum extent of the differential preferences.

    * *ylimits*: is **mandatory** if *datatype* is 'diffsel', and meaningless 
      otherwise. It is *(ymin, ymax)* where *ymax > 0 > ymin*, and gives extent 
      of the data in the positive and negative directions. Must encompass the 
      actual maximum and minimum of the data.

    * *overlay* : make overlay bars that indicate other properties for
      the sites. If you set to something other than `None`, it should be
      a list giving one to three properties. Each property is a tuple:
      *(prop_d, shortname, longname)* where:

        - *prop_d* is a dictionary keyed by site numbers that are in *sites*.
          For each *r* in *sites*, *prop_d[r]* gives the value of the property,
          or if there is no entry in *prop_d* for *r*, then the property
          is undefined and is colored white. Properties can either be:

            * continuous: in this case, all of the values should be numbers.

            * discrete : in this case, all of the values should be strings.
              While in practice, if you have more than a few discrete
              categories (different strings), the plot will be a mess.

        - *shortname* : short name for the property; will not format well
          if more than 4 or 5 characters.

        - *longname* : longer name for property used on axes label. Can be the
          same as *shortname* if you don't need a different long name.

        - In the special case where both *shortname* and *longname* are 
          the string `wildtype`, then rather than an overlay bar we
          right the one-character wildtype identity in `prop_d` for each
          site.

    * *fix_limits* is only meaningful if *overlay* is being used. In this case, for any
      *shortname* in *overlay* that also keys an entry in *fix_limits*, we use
      *fix_limits[shortname]* to set the limits for *shortname*. Specifically,
      *fix_limits[shortname]* should be the 2-tuple *(ticks, ticknames)*. *ticks*
      should be a list of tick locations (numbers) and *ticknames* should be a list of
      the corresponding tick label for that tick.

    * If *fixlongname* is *True*, then we use the *longname* in *overlay* exactly as written;
      otherwise we add a parenthesis indicating the *shortname* for which this *longname*
      stands.

    * *overlay_cmap* can be the name of a valid *matplotlib.colors.Colormap*, such as the
      string *jet* or *bwr*. Otherwise, it can be *None* and a (hopefully) good choice will 
      be made for you.

    * *custom_cmap* can be the name of a valid *matplotlib.colors.Colormap* which will be
      used to color amino-acid one-letter codes in the logoplot by the *map_metric* when
      either 'kd' or 'mw' is used as *map_metric*.

    * *relativestackheight* indicates how high the letter stack is relative to
      the default. The default is multiplied by this number, so make it > 1
      for a higher letter stack.

    * *map_metric* specifies the amino-acid property metric used to map colors to amino-acid
      letters. Valid options are 'kd' (Kyte-Doolittle hydrophobicity scale, default), 'mw' 
      (molecular weight), 'functionalgroup' (functional groups: small, nucleophilic, hydrophobic,
      aromatic, basic, acidic, and amide), and 'charge' (charge at neutral pH). If 'charge' is used, then the
      argument for *custom_cmap* will no longer be meaningful, since 'charge' uses its own 
      blue/black/red colormapping. Similarly, 'functionalgroup' uses its own colormapping.

    * *noseparator* is only meaningful if *datatype* is 'diffsel' or 'diffprefs'.
      If it set to *True*, then we do **not** print a black horizontal line to
      separate positive and negative values.

    * *underlay* if `True` then make an underlay rather than an overlay.

    * *scalebar*: show a scale bar. If `False`, no scale bar shown. Otherwise
      should be a 2-tuple of `(scalebarlen, scalebarlabel)`. Currently only
      works when data is `diffsel`.

def LogoPlot(sites, datatype, data, plotfile, nperline,
        numberevery=10, allowunsorted=False, ydatamax=1.01,
        overlay=None, fix_limits={}, fixlongname=False,
        overlay_cmap=None, ylimits=None, relativestackheight=1,
        custom_cmap='jet', map_metric='kd', noseparator=False,
        underlay=False, scalebar=False):
    """Create sequence logo showing amino-acid or nucleotide preferences.

    The heights of each letter is equal to the preference of
    that site for that amino acid or nucleotide.

    Note that stop codons may or may not be included in the logo
    depending on whether they are present in *pi_d*.  

    CALLING VARIABLES:

    * *sites* is a list of all of the sites that are being included
      in the logo, as strings. They must be in natural sort or an error
      will be raised **unless** *allowunsorted* is *True*. The sites
      in the plot are ordered in the same arrangement
      listed in *sites*. These should be **strings**, not integers.

    * *datatype* should be one of the following strings:
    
        * 'prefs' for preferences
        
        * 'diffprefs' for differential preferences
        
        * 'diffsel' for differential selection

    * *data* is a dictionary that has a key for every entry in
      *sites*. For every site *r* in *sites*, *sites[r][x]*
      is the value for character *x*. 
      Preferences must sum to one; differential preferences to zero.
      All sites must have the same set of characters. The characters
      must be the set of nucleotides or amino acids with or without
      stop codons.

    * *plotfile* is a string giving the name of the created PDF file 
      of the logo plot.
      It must end in the extension ``.pdf``.

    * *nperline* is the number of sites per line. Often 40 to 80 are good values.

    * *numberevery* is specifies how frequently we put labels for the sites on
      x-axis.

    * *allowunsorted* : if *True* then we allow the entries in *sites* to 
      **not** be sorted. This means that the logo plot will **not** have
      sites in sorted order.

    * *ydatamax* : meaningful only if *datatype* is 'diffprefs'. In this case, it gives
      the maximum that the logo stacks extend in the positive and negative directions.
      Cannot be smaller than the maximum extent of the differential preferences.

    * *ylimits*: is **mandatory** if *datatype* is 'diffsel', and meaningless 
      otherwise. It is *(ymin, ymax)* where *ymax > 0 > ymin*, and gives extent 
      of the data in the positive and negative directions. Must encompass the 
      actual maximum and minimum of the data.

    * *overlay* : make overlay bars that indicate other properties for
      the sites. If you set to something other than `None`, it should be
      a list giving one to three properties. Each property is a tuple:
      *(prop_d, shortname, longname)* where:

        - *prop_d* is a dictionary keyed by site numbers that are in *sites*.
          For each *r* in *sites*, *prop_d[r]* gives the value of the property,
          or if there is no entry in *prop_d* for *r*, then the property
          is undefined and is colored white. Properties can either be:

            * continuous: in this case, all of the values should be numbers.

            * discrete : in this case, all of the values should be strings.
              While in practice, if you have more than a few discrete
              categories (different strings), the plot will be a mess.

        - *shortname* : short name for the property; will not format well
          if more than 4 or 5 characters.

        - *longname* : longer name for property used on axes label. Can be the
          same as *shortname* if you don't need a different long name.

        - In the special case where both *shortname* and *longname* are 
          the string `wildtype`, then rather than an overlay bar we
          right the one-character wildtype identity in `prop_d` for each
          site.

    * *fix_limits* is only meaningful if *overlay* is being used. In this case, for any
      *shortname* in *overlay* that also keys an entry in *fix_limits*, we use
      *fix_limits[shortname]* to set the limits for *shortname*. Specifically,
      *fix_limits[shortname]* should be the 2-tuple *(ticks, ticknames)*. *ticks*
      should be a list of tick locations (numbers) and *ticknames* should be a list of
      the corresponding tick label for that tick.

    * If *fixlongname* is *True*, then we use the *longname* in *overlay* exactly as written;
      otherwise we add a parenthesis indicating the *shortname* for which this *longname*
      stands.

    * *overlay_cmap* can be the name of a valid *matplotlib.colors.Colormap*, such as the
      string *jet* or *bwr*. Otherwise, it can be *None* and a (hopefully) good choice will 
      be made for you.

    * *custom_cmap* can be the name of a valid *matplotlib.colors.Colormap* which will be
      used to color amino-acid one-letter codes in the logoplot by the *map_metric* when
      either 'kd' or 'mw' is used as *map_metric*.

    * *relativestackheight* indicates how high the letter stack is relative to
      the default. The default is multiplied by this number, so make it > 1
      for a higher letter stack.

    * *map_metric* specifies the amino-acid property metric used to map colors to amino-acid
      letters. Valid options are 'kd' (Kyte-Doolittle hydrophobicity scale, default), 'mw' 
      (molecular weight), 'functionalgroup' (functional groups: small, nucleophilic, hydrophobic,
      aromatic, basic, acidic, and amide), and 'charge' (charge at neutral pH). If 'charge' is used, then the
      argument for *custom_cmap* will no longer be meaningful, since 'charge' uses its own 
      blue/black/red colormapping. Similarly, 'functionalgroup' uses its own colormapping.

    * *noseparator* is only meaningful if *datatype* is 'diffsel' or 'diffprefs'.
      If it set to *True*, then we do **not** print a black horizontal line to
      separate positive and negative values.

    * *underlay* if `True` then make an underlay rather than an overlay.

    * *scalebar*: show a scale bar. If `False`, no scale bar shown. Otherwise
      should be a 2-tuple of `(scalebarlen, scalebarlabel)`. Currently only
      works when data is `diffsel`.
    """
    assert datatype in ['prefs', 'diffprefs', 'diffsel'], "Invalid datatype {0}".format(datatype)

    # check data, and get characters
    assert sites, "No sites specified"
    assert set(sites) == set(data.keys()), "Not a match between sites and the keys of data"
    characters = list(data[sites[0]].keys())
    aas = sorted(AA_TO_INDEX.keys())
    if set(characters) == set(NT_TO_INDEX.keys()):
        alphabet_type = 'nt'
    elif set(characters) == set(aas) or set(characters) == set(aas + ['*']):
        alphabet_type = 'aa'
    else:
        raise ValueError("Invalid set of characters in data. Does not specify either nucleotides or amino acids:\n%s" % str(characters))
    for r in sites:
        if set(data[r].keys()) != set(characters):
            raise ValueError("Not all sites in data have the same set of characters")

    firstblankchar = 'B' # character for first blank space for diffprefs / diffsel
    assert firstblankchar not in characters, "firstblankchar in characters"
    lastblankchar = 'b' # character for last blank space for diffprefs / diffsel
    assert lastblankchar not in characters, "lastblankchar in characters"
    separatorchar = '-' # separates positive and negative for diffprefs / diffsel
    assert separatorchar not in characters, "lastblankchar in characters"
    if noseparator:
        separatorheight = 0
    else:
        separatorheight = 0.02 # height of separator as frac of total for diffprefs / diffsel

    if os.path.splitext(plotfile)[1].lower() != '.pdf':
        raise ValueError("plotfile must end in .pdf: %s" % plotfile)
    if os.path.isfile(plotfile):
        os.remove(plotfile) # remove existing plot

    if not allowunsorted:
        sorted_sites = natsort.natsorted([r for r in sites])
        if sorted_sites != sites:
            raise ValueError("sites is not properly sorted")

    # Following are specifications of weblogo sizing taken from its documentation
    stackwidth = 9.5 # stack width in points, not default size of 10.8, but set to this in weblogo call below
    barheight = 5.5 # height of bars in points if using overlay
    barspacing = 2.0 # spacing between bars in points if using overlay
    stackaspectratio = 4.4 # ratio of stack height:width, doesn't count part going over maximum value of 1
    assert relativestackheight > 0, "relativestackheight must be > 0"
    stackaspectratio *= relativestackheight
    if overlay:
        if len(overlay) > 3:
            raise ValueError("overlay cannot have more than 3 entries")
        ymax = (stackaspectratio * stackwidth + len(overlay) * (barspacing + barheight)) / float(stackaspectratio * stackwidth)
        aspectratio = ymax * stackaspectratio # effective aspect ratio for full range
    else:
        ymax = 1.0 
        aspectratio = stackaspectratio
    rmargin = 11.5 # right margin in points, fixed by weblogo
    stackheightmargin = 16 # margin between stacks in points, fixed by weblogo

    showscalebar = False
    try:
        # write data into transfacfile (a temporary file)
        (fd, transfacfile) = tempfile.mkstemp()
        f = os.fdopen(fd, 'w')
        ordered_alphabets = {} # keyed by site index (0, 1, ...) with values ordered lists for characters from bottom to top
        if datatype == 'prefs':
            chars_for_string = characters
            f.write('ID ID\nBF BF\nP0 %s\n' % ' '.join(chars_for_string))
            for (isite, r) in enumerate(sites):
                f.write('%d %s\n' % (isite, ' '.join([str(data[r][x]) for x in characters])))
                pi_r = [(data[r][x], x) for x in characters]
                pi_r.sort()
                ordered_alphabets[isite] = [tup[1] for tup in pi_r] # order from smallest to biggest
        elif datatype == 'diffprefs':
            chars_for_string = characters + [firstblankchar, lastblankchar, separatorchar]
            ydatamax *= 2.0 # maximum possible range of data, multiply by two for range
            f.write('ID ID\nBF BF\nP0 %s\n' % ' '.join(chars_for_string))
            for (isite, r) in enumerate(sites):
                positivesum = sum([data[r][x] for x in characters if data[r][x] > 0]) + separatorheight / 2.0
                negativesum = sum([data[r][x] for x in characters if data[r][x] < 0]) - separatorheight / 2.0
                if abs(positivesum + negativesum) > 1.0e-3:
                    raise ValueError("Differential preferences sum of %s is not close to zero for site %s" % (positivesum + negativesum, r))
                if 2.0 * positivesum > ydatamax:
                    raise ValueError("You need to increase ydatamax: the total differential preferences sum to more than the y-axis limits. Right now, ydatamax is %.3f while the total differential preferences are %.3f" % (ydatamax, 2.0 * positivesum))
                f.write('%d' % isite)
                deltapi_r = []
                for x in characters:
                    deltapi_r.append((data[r][x], x))
                    f.write(' %s' % (abs(data[r][x]) / float(ydatamax)))
                deltapi_r.sort()
                firstpositiveindex = 0
                while deltapi_r[firstpositiveindex][0] < 0:
                    firstpositiveindex += 1
                ordered_alphabets[isite] = [firstblankchar] + [tup[1] for tup in deltapi_r[ : firstpositiveindex]] + [separatorchar] + [tup[1] for tup in deltapi_r[firstpositiveindex : ]] + [lastblankchar] # order from most negative to most positive with blank characters and separators
                f.write(' %g %g %g\n' % (0.5 * (ydatamax + 2.0 * negativesum) / ydatamax, 0.5 * (ydatamax + 2.0 * negativesum) / ydatamax, separatorheight)) # heights for blank charactors and separators
        elif datatype == 'diffsel':
            assert ylimits, "You must specify ylimits if using diffsel"
            (dataymin, dataymax) = ylimits
            assert dataymax > 0 > dataymin, "Invalid ylimits of {0}".format(ylimits)
            yextent = float(dataymax - dataymin)
            separatorheight *= yextent
            chars_for_string = characters + [firstblankchar, lastblankchar, separatorchar]
            f.write('ID ID\nBF BF\nP0 {0}\n'.format(' '.join(chars_for_string)))
            for (isite, r) in enumerate(sites):
                positivesum = sum([data[r][x] for x in characters if data[r][x] > 0]) + separatorheight / 2.0
                negativesum = sum([data[r][x] for x in characters if data[r][x] < 0]) - separatorheight / 2.0
                assert positivesum <= dataymax, "Data exceeds ylimits in positive direction"
                assert negativesum >= dataymin, "Data exceeds ylimits in negative direction"
                f.write('{0}'.format(isite))
                diffsel_r = []
                for x in characters:
                    diffsel_r.append((data[r][x], x))
                    f.write(' {0}'.format(abs(data[r][x]) / yextent))
                diffsel_r.sort()
                firstpositiveindex = 0
                while diffsel_r[firstpositiveindex][0] < 0:
                    firstpositiveindex += 1
                ordered_alphabets[isite] = [firstblankchar] + [tup[1] for tup in diffsel_r[ : firstpositiveindex]] + [separatorchar] + [tup[1] for tup in diffsel_r[firstpositiveindex : ]] + [lastblankchar] # order from most negative to most positive with blank characters and separators
                f.write(' %g %g %g\n' % ((negativesum - dataymin) / yextent, (dataymax - positivesum) / yextent, separatorheight / yextent)) # heights for blank charactors and separators
            # height of one unit on y-axis in points
            heightofone = stackwidth * stackaspectratio / yextent
            assert heightofone > 0
            if scalebar:
                showscalebar = (heightofone * scalebar[0], scalebar[1])
        else:
            raise ValueError("Invalid datatype of %s" % datatype)
        f.close()

        # create web logo
        charstring = ''.join(chars_for_string)
        assert len(charstring) == len(chars_for_string), "Length of charstring doesn't match length of chars_for_string. Do you have unallowable multi-letter characters?\n%s" % (str(chars_for_string))
        logoprior = weblogolib.parse_prior('equiprobable', charstring, 0)
        motif = _my_Motif.read_transfac(open(transfacfile), charstring)
        logodata = weblogolib.LogoData.from_counts(motif.alphabet, motif, logoprior)
        logo_options = weblogolib.LogoOptions()
        logo_options.fineprint = None
        logo_options.stacks_per_line = nperline
        logo_options.stack_aspect_ratio = aspectratio
        logo_options.stack_width = stackwidth
        logo_options.unit_name = 'probability'
        logo_options.show_yaxis = False
        logo_options.yaxis_scale = ymax 
        if alphabet_type == 'aa':
            map_functions = {'kd':KyteDoolittleColorMapping,
                             'mw': MWColorMapping,
                             'charge' : ChargeColorMapping,
                             'functionalgroup':FunctionalGroupColorMapping}
            map_fcn = map_functions[map_metric]
            (cmap, colormapping, mapper) = map_fcn(maptype=custom_cmap)
        elif alphabet_type == 'nt':
            colormapping = {}
            colormapping['A'] = '#008000'
            colormapping['T'] = '#FF0000'
            colormapping['C'] = '#0000FF'
            colormapping['G'] = '#FFA500'
        else:
            raise ValueError("Invalid alphabet_type %s" % alphabet_type)
        colormapping[firstblankchar] = colormapping[lastblankchar] = '#000000' # black, but color doesn't matter as modified weblogo code replaces with empty space
        colormapping[separatorchar] = '#000000' # black
        color_scheme = weblogolib.colorscheme.ColorScheme()
        for x in chars_for_string:
            if hasattr(color_scheme, 'rules'):
                color_scheme.rules.append(weblogolib.colorscheme.SymbolColor(x, colormapping[x], "'%s'" % x))
            else:
                # this part is needed for weblogo 3.4
                color_scheme.groups.append(weblogolib.colorscheme.ColorGroup(x, colormapping[x], "'%s'" % x))
        logo_options.color_scheme = color_scheme
        logo_options.annotate = [{True:r, False:''}[0 == isite % numberevery] for (isite, r) in enumerate(sites)]
        logoformat = weblogolib.LogoFormat(logodata, logo_options)
        # _my_pdf_formatter is modified from weblogo version 3.4 source code
        # to allow custom ordering of the symbols.
        pdf = _my_pdf_formatter(logodata, logoformat, ordered_alphabets) 
        with open(plotfile, 'wb') as f:
            f.write(pdf)
        assert os.path.isfile(plotfile), "Failed to find expected plotfile %s" % plotfile
    finally:
        # close if still open
        try:
            f.close()
        except:
            pass
        # remove temporary file
        if os.path.isfile(transfacfile):
            os.remove(transfacfile)

    # now build the overlay
    if overlay or showscalebar:
        try:
            (fdoverlay, overlayfile) = tempfile.mkstemp(suffix='.pdf')
            (fdmerged, mergedfile) = tempfile.mkstemp(suffix='.pdf')
            foverlay = os.fdopen(fdoverlay, 'wb')
            foverlay.close() # close, but we still have the path overlayfile...
            fmerged = os.fdopen(fdmerged, 'wb')
            logoheight = stackwidth * stackaspectratio + stackheightmargin
            LogoOverlay(sites, overlayfile, overlay, nperline, sitewidth=stackwidth, rmargin=rmargin, logoheight=logoheight, barheight=barheight, barspacing=barspacing, fix_limits=fix_limits, fixlongname=fixlongname, overlay_cmap=overlay_cmap, underlay=underlay, scalebar=showscalebar)
            plotfile_f = open(plotfile, 'rb')
            plot = PyPDF2.PdfFileReader(plotfile_f).getPage(0)
            overlayfile_f = open(overlayfile, 'rb')
            overlaypdf = PyPDF2.PdfFileReader(overlayfile_f).getPage(0)
            xshift = overlaypdf.artBox[2] - plot.artBox[2]
            yshift = (barheight + barspacing) * len(overlay) - 0.5 * barspacing
            overlaypdf.mergeTranslatedPage(plot, xshift, 
                    yshift * int(underlay), expand=True)
            overlaypdf.compressContentStreams() 
            output = PyPDF2.PdfFileWriter()
            output.addPage(overlaypdf)
            output.write(fmerged)
            fmerged.close()
            shutil.move(mergedfile, plotfile)
        finally:
            try:
                plotfile_f.close()
            except:
                pass
            try:
                overlayfile_f.close()
            except:
                pass
            try:
                foverlay.close()
            except:
                pass
            try:
                fmerged.close()
            except:
                pass
            for fname in [overlayfile, mergedfile]:
                if os.path.isfile(fname):
                    os.remove(fname)