--- a
+++ b/singlecellmultiomics/utils/bdbplot.py
@@ -0,0 +1,3199 @@
+from lxml import etree
+import math
+import collections
+from collections import Counter
+from collections import OrderedDict
+import numpy as np
+from singlecellmultiomics.utils import bdbbio
+import os
+import matplotlib.cm
+from Bio import SeqIO
+from Bio.Seq import Seq
+import matplotlib.pyplot as plt
+from colorama import Fore #,Back, Style
+from colorama import Back
+from colorama import Style
+from colorama import init
+import scipy
+import scipy.cluster
+import time
+import itertools
+
+init(autoreset=True)
+
+
+
+#Convert a nested dictionary to a matrix
+# ({'A':{'1':2}, 'B':{'1':3, '2':4}}) will become
+#(array([[  2.,  nan],
+#        [  3.,   4.]]), ['A', 'B'], ['1', '2'])
+
+
+def interpolateBezier( points, steps=10, t=None):
+    if len(points)==3:
+        mapper = lambda t,p: (1-t)**2 * p[0] + 2*(1-t)*t*p[1] + t**2*p[2]
+    elif len(points)==4:
+        mapper = lambda t,p: (np.power( (1-t),3)*p[0] +\
+         3* np.power((1-t),2) *t *p[1] +\
+         3*(1-t)*np.power(t,2)*p[2] +\
+         np.power(t,3)*p[3])
+
+    if t is not None:
+        return   mapper(t, [q[0] for q in points]), mapper(t, [q[1] for q in points])
+    xGen = ( mapper(t, [q[0] for q in points]) for t in np.linspace(0, 1, steps) )
+    yGen = ( mapper(t, [q[1] for q in points]) for t in np.linspace(0, 1, steps) )
+
+    return zip(xGen, yGen)
+
+def interpolateBezierAngle(points, t, ds=0.001):
+    x0, y0 = interpolateBezier(points, t=t-ds)
+    x1, y1 = interpolateBezier(points, t=t+ds)
+    return np.arctan2( y1-y0, x0-x1)
+
+
+def initMatrix(rowNames,columnNames, mtype="obj"):
+    if mtype=="obj":
+        matrix = np.empty( (len(rowNames), len(columnNames)), dtype=object)
+    elif mtype=="npzeros":
+        matrix = np.zeros( (len(rowNames), len(columnNames)))
+    return(matrix)
+
+def nestedDictionaryToNumpyMatrix( nestedDictionary, setNan=True, mtype="obj", transpose=False, indicateProgress=False):
+
+    rowNames = sorted( list(nestedDictionary.keys()), key=int )
+    columnNames = set()
+    for key in nestedDictionary:
+        columnNames.update( set(nestedDictionary[key].keys() ))
+    columnNames = sorted( list(columnNames) )
+
+    keys = list(columnNames)
+
+    if ':' in keys[0]:
+        sargs = np.argsort( [ int(k.split(':')[1]) for k in keys] )
+        print(sargs)
+
+        columnNames = [ keys[index] for index in sargs ]
+        print(columnNames)
+
+    matrix = initMatrix(rowNames,columnNames, mtype)
+
+    if setNan:
+        matrix[:] = np.nan
+
+    prevTime = time.time()
+    for rowIndex,rowName in enumerate(rowNames):
+        if indicateProgress and (time.time()-prevTime)>1:
+            prevTime = time.time()
+            print("Matrix creation progress: %.2f%%" % (100.0*rowIndex/len(rowNames)))
+
+        for colIndex,colName in enumerate(columnNames):
+            try:
+                matrix[rowIndex,colIndex] = nestedDictionary[rowName][colName]
+            except:
+                pass
+
+    if transpose:
+        matrix = matrix.transpose()
+        columnNames, rowNames =  rowNames, columnNames
+
+    if indicateProgress:
+        print("Matrix finished")
+    return( (matrix, rowNames, columnNames) )
+
+def pruneNonUniqueColumnsFromMatrix(matrix,rows,columns, minInstances=1,minOccurence=1):
+    colsToKeep = []
+    for columnIndex in range(matrix.shape[1]):
+        if len( set(np.unique( matrix[:,columnIndex].astype(str) ))-set( ["nan"]) )>minInstances:
+
+            cnts = Counter( list(matrix[:,columnIndex].astype(str)) )
+
+            counts = Counter({k: cnts[k] for k in  cnts if cnts[k] >= minOccurence})
+            del counts['nan']
+
+            if len(counts.values())>minInstances:
+                colsToKeep.append(columnIndex)
+
+    matrix = matrix[:,colsToKeep]
+    columns = np.array(columns)[colsToKeep]
+    return(matrix, rows, columns)
+
+
+
+# Convert dictionary of tuples to a numpy matrix
+def tupleAnnotationsToNumpyMatrix( originRowNames, originColNames, tuples, setNan=True, mtype="obj" ):
+
+    m = initMatrix(originRowNames,originColNames, mtype) #np.zeros( (len(originRowNames), len(originColNames)))
+    if setNan:
+        m[:] = np.nan
+    for tup in tuples:
+        value = tuples[tup]
+        if tup[0] in originRowNames and tup[1] in originColNames:
+         m[originRowNames.index(tup[0]), originColNames.index(tup[1])  ] = value
+    return(m)
+
+def dictAnnotationsToNumpyMatrix( originRowNames, originColNames, dictionary, mtype="obj" ):
+    tuples = {}
+    for rowKey in dictionary:
+        for columnKey in dictionary[rowKey]:
+            tuples[ (rowKey, columnKey) ] = dictionary[rowKey][columnKey]
+    return(tupleAnnotationsToNumpyMatrix(originRowNames, originColNames, tuples, mtype=mtype))
+
+
+def getSomeColors(n):
+    return( plt.cm.Set1(np.linspace(0, 1, n)) )
+
+
+
+def _ipol(a, b, first, last,  interpolateValue):
+    #Due to floating point rounding errors the interpolate value can be very close to last,
+    # it is ok to return last in those cases
+    if last>first and interpolateValue>=last:
+        return(b)
+    if last<first and interpolateValue>=first:
+        return(a)
+
+    y_interp = scipy.interpolate.interp1d([first, last], [a,b])
+    return( y_interp(interpolateValue) )
+
+def interpolate(interpolateValue,  colorScaleKeys, nodeColorMapping):
+        #Seek positions around value to interpolate
+        first = colorScaleKeys[0]
+        index = 0
+        last = first
+        for value in colorScaleKeys:
+            if value>=interpolateValue:
+                last = value
+                break
+            else:
+                first = value
+            index+=1
+        if value==interpolateValue:
+            return(nodeColorMapping[value])
+
+        #Do interpolation
+        colorA = nodeColorMapping[first]
+        colorB = nodeColorMapping[last]
+        dx = last-first
+
+        # Check out of bounds condition
+        if interpolateValue< first:
+            return(colorA)
+        if interpolateValue>last:
+            return(colorB)
+
+
+
+        return( _ipol(colorA[0], colorB[0], first, last, interpolateValue), _ipol(colorA[1], colorB[1], first, last, interpolateValue), _ipol(colorA[2], colorB[2], first, last, interpolateValue))
+
+
+
+def plotFeatureSpace(features, classLabels, featureNames, path=None, bins = 50, title=None):
+
+    print(Fore.GREEN + "Feature space plotter:")
+    #1d mode:
+    print(features.shape)
+    classAbundance = Counter(classLabels)
+    print(classAbundance)
+
+    classes = set(classLabels)
+    classColors = ['#FF6A00','#0066FF','#FF33FF','#666666']
+    classColors = [ tuple(float(int(hexColor.replace('#','')[i:i+2], 16))/255.0 for i in (0, 2 ,4)) for hexColor in classColors]
+
+    if features.shape[1]==1:
+        print("Performing 1-D density plot of %s samples" % ( features.shape[0]))
+
+        plt.close('all')
+
+        histStart = features.min()
+        histEnd = features.max()
+        if histStart == histEnd:
+            histEnd += 1
+            histStart -= 1
+        precision = (histEnd-histStart)/bins
+
+
+        print("Histogram will be plotted from %s to %s " % (histStart, histEnd))
+        fig, ax = plt.subplots() #figsize=(120, 10))
+
+        print(classColors)
+        for classIndex,className in enumerate(list(classes)):
+            boolList = np.array(classLabels)==np.array(className)
+            classSize = classAbundance[className]
+            ax.hist(
+                features[boolList,0],
+                np.arange(histStart,histEnd+precision,precision),
+                normed=False, fc= (classColors[classIndex]+ (0.5,)),
+                ec=classColors[classIndex],
+                lw=1.5, histtype='stepfilled',
+                label='%s[%s]' % (className,classSize)
+                )
+
+        plt.ylabel("Density")
+        plt.xlabel(featureNames[0])
+        if title is not None:
+            plt.title(title)
+        ax.legend(loc='upper right')
+        #plt.yscale('log', nonposy='clip')
+        if path is None:
+            plt.show()
+        else:
+            plt.savefig(path)
+        return(True)
+
+
+    if features.shape[1]==2:
+        #2d
+        print("Performing 2-D density plot of %s samples" % ( features.shape[0]))
+        fig, ax = plt.subplots()
+
+        for classIndex,className in enumerate(list(classes)):
+            #print(np.where( classLabels==className, features ))
+            boolList = np.array(classLabels)==np.array(className)
+            plt.plot( features[boolList,0], features[boolList,1], ".",label='%s[%s]' % (className,sum(boolList)), c= (classColors[classIndex] + (0.5,)))
+
+        plt.xlabel(featureNames[0])
+        plt.ylabel(featureNames[1])
+        plt.legend(loc="lower right")
+        plt.tight_layout()
+        try:
+            plt.savefig(path)
+
+        except Exception as e:
+            print(e)
+        return(True)
+    print(Fore.RED + "Invalid amount of dimensions for feature space plotting (%s)" % features.shape[1])
+
+def matplotHeatmap( D, YC, figsize=(10,10), clust=True, xLab=None, yLab=None, show=True, colormap=plt.cm.YlGnBu_r, colorbarLabel=None ):
+    plt.rcParams["axes.grid"] = False
+    import scipy
+    import scipy.cluster.hierarchy as sch
+    # Compute and plot first dendrogram.
+    fig = plt.figure(figsize=figsize)
+    if not clust:
+        idx1 = range(0, len(YC))
+        idx1 = range(0, len(YC))
+
+    if clust:
+        ax1 = fig.add_axes([0.09,0.1,0.2,0.6])
+
+        L = sch.linkage(D, method='centroid')
+        Z1 = sch.dendrogram(L, orientation='right')
+        ax1.set_xticks([])
+        ax1.set_yticks([])
+
+        # Compute and plot second dendrogram.
+        ax2 = fig.add_axes([0.3,0.71,0.6,0.2])
+        Z2 = sch.dendrogram(L)
+        ax2.set_xticks([])
+        ax2.set_yticks([])
+        idx1 = Z1['leaves']
+        idx2 = Z2['leaves']
+        D = D[idx1,:]
+        D = D[:,idx1]
+
+    # Plot distance matrix.
+    axmatrix = fig.add_axes([0.3,0.1,0.6,0.6])
+
+    im = axmatrix.matshow(D, aspect='auto', origin='lower', cmap=colormap)
+    axmatrix.set_xticks([])
+    axmatrix.set_yticks([])
+    ###
+
+    if xLab is None:
+        axmatrix.set_xticks(range(len(YC)))
+        axmatrix.set_xticklabels(YC[idx1])
+    else:
+        axmatrix.set_xticks(range(len(xLab)))
+        axmatrix.set_xticklabels(xLab)
+    axmatrix.xaxis.set_label_position('bottom')
+    axmatrix.xaxis.tick_bottom()
+
+    plt.xticks(rotation=-90)
+
+    if yLab is None:
+        axmatrix.set_yticks(range(len(YC)))
+        axmatrix.set_yticklabels(YC[idx1], minor=False)
+    else:
+        axmatrix.set_yticks(range(len(yLab)))
+        axmatrix.set_yticklabels(yLab, minor=False)
+    axmatrix.yaxis.set_label_position('right')
+    axmatrix.yaxis.tick_right()
+
+
+    # Plot colorbar.
+    axcolor = fig.add_axes([1.05,0.1,0.02,0.6])
+    cbar = fig.colorbar(im, cax=axcolor)
+    if colorbarLabel is not None:
+        cbar.set_label(colorbarLabel, rotation=270,  labelpad=15)
+    if show:
+        fig.savefig('dendrogram.png')
+
+def tsnePlot(data, labels=None, components=2, perplexity=30.0, iterations=1000):
+    from sklearn.manifold import TSNE
+    #from MulticoreTSNE import MulticoreTSNE as TSNE
+    model = TSNE(n_components=components, perplexity=perplexity, n_iter=iterations ) #random_state=0, n_jobs=8,
+    transformedPoints = model.fit_transform(data.astype(np.float64))
+
+
+    classes = list(set(labels))
+    classColors = getSomeColors(len(classes))
+    color = np.array([ classColors[classes.index(label)] for label in labels ])
+    nplabels = np.array(labels)
+    print("TSNE input:")
+    print(data.shape)
+    print("TSNE plotting for %s classes " %  len(classes))
+
+    #print("Color mapping is %s" % ",".join(color))
+    #Plot the data:
+    if components==2:
+        fig = plt.figure()
+        #plt.style.use('ggplot')
+        ax = fig.add_subplot(111)
+
+        print(color)
+        #plt.scatter(transformedPoints[:, 0], transformedPoints[:, 1], c=color, cmap=plt.cm.Spectral, s=1, alpha=0.5) #labels=labels,
+        for classIndex, className in enumerate(classes):
+            classColor = classColors[classIndex]
+            print(className)
+            print(classColor)
+            plt.scatter(transformedPoints[className==nplabels, 0], transformedPoints[className==nplabels, 1],   s=3, alpha=0.9, label=className)  #c=classColor,, cmap=plt.cm.Spectral,
+
+        #plt.axis('tight')
+        box = ax.get_position()
+        ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
+
+        # Put a legend to the right of the current axis
+        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
+
+        plt.show()
+
+
+    elif components==3:
+        from mpl_toolkits.mplot3d import Axes3D
+        fig = plt.figure()
+        ax = fig.add_subplot(111, projection='3d')
+
+        for classIndex,className in enumerate(classes):
+            samplesForClass = (className==nplabels)
+            ax.scatter(transformedPoints[samplesForClass,0], transformedPoints[samplesForClass,1], transformedPoints[samplesForClass,2], c=classColors[classIndex], label=classIndex)
+        ax.legend()
+        plt.show()
+
+
+class BDBcolor():
+
+    def __init__(self, r=0, g=0, b=0, a=1.0 ):
+
+        if str(r)[0]=='#':
+            #parse hex colour:
+            #parts = r.replace('#','').replace('(','').replace(')','').split(',')
+            cleaned = r.replace('#','')
+
+            #parts = [int(i) for i in parts]
+            r =  int(cleaned[0:2], 16)
+            g =  int(cleaned[2:4], 16)
+            b =  int(cleaned[4:6], 16)
+
+        self.r = max(0,min(255,r))
+        self.g = max(0,min(255,g))
+        self.b = max(0,min(255,b))
+        self.a = max(0,min(1.0,a))
+
+    def getRGBStr(self):
+        return('rgb(%s,%s,%s)' % (int(self.r), int(self.g), int(self.b)))
+
+    def getRGBAStr(self):
+        return('rgba(%s,%s,%s,%s)' % (self.r, self.g, self.b, self.a))
+
+
+    def getReadableInverted(self):
+
+        hsv = self.getHSV()
+        hsv['v'] = (  255.0-( hsv['v'] ) )
+        rgb = self.HSVtoRGB(0,0,hsv['v'])
+        return( BDBcolor( rgb['r'],rgb['g'],rgb['b'] ))
+
+    def getHSV(self):
+        h=0
+        s=0
+        v=0
+        minV = min( self.r, self.g, self.b )
+        maxV = max( self.r, self.g, self.b )
+        v = maxV
+        delta = maxV - minV
+        if maxV != 0:
+            s = delta / float(maxV)
+        else:
+            s = 0
+            h = -1
+            return({'h':h, 's':s, 'v':v})
+
+        if delta==0:
+                h = 255
+        else:
+                if self.r == maxV:
+                    h = ( self.g - self.b ) / float(delta)
+
+                else:
+                        if self.g == maxV:
+
+                                h = 2.0 + float( self.b - self.r ) / float(delta)
+                        else:
+                                h = 4.0 + float( self.r - self.g ) / float(delta)
+
+
+        h *= 60
+        if h < 0:
+            h += 360
+
+        return({'h':h, 's':s, 'v':v})
+
+
+    def HSVtoRGB( self, h,s,v ):
+
+
+        i=0
+        f=0
+        p=0
+        q=0
+        t = 0
+
+        if s == 0:
+            #Grey
+            r = g = b = v
+            return({'r':round(r), 'g':round(g), 'b':round(b)})
+
+        h /= 60         # sector 0 to 5
+        i = math.floor( h )
+        f = h - i           # factorial part of h
+        p = v * ( 1 - s )
+        q = v * ( 1 - s * f )
+        t = v * ( 1 - s * ( 1 - f ) )
+
+        if i==0:
+            r = v
+            g = t
+            b = p
+        elif i==1:
+            r = q
+            g = v
+            b = p
+        elif i==2:
+            r = p
+            g = v
+            b = t
+        elif i==3:
+            r = p
+            g = q
+            b = v
+        elif i==4:
+            r = t
+            g = p
+            b = v
+        else:
+            r = v
+            g = p
+            b = q
+
+        return({'r':round(r), 'g':round(g), 'b':round(b)})
+
+
+class BDBPlot():
+
+
+    def __init__(self):
+
+        # We need to declare the xlink namespace, to create references to things in our own file
+        self.xlink =  'http://www.w3.org/1999/xlink'
+
+        NSMAP = {'xlink':self.xlink }
+
+        self.svgTree = etree.Element('svg',nsmap = NSMAP)
+        self.svgTree.set('xmlns','http://www.w3.org/2000/svg')
+        self.svgTree.set('version','1.2')
+
+
+
+        self.root = self.svgTree.getroottree()
+
+        self.nextFilterId = 0
+        self.nextObjId = 0
+        self.nextTspanId = 0
+        #Create definition element
+        self.svgTree.append( self.getDefinitionBlock() )
+        self.debug = 2 # 2 all
+
+        self.xMin = 0
+        self.xMax = 10
+        self.yMin = 0
+        self.yMax = 10
+        self.plotStartX = 30
+        self.plotStartY = 30
+        self.plotHeight = 400
+        self.plotWidth = 600
+
+        self.setWidth(800)
+        self.setHeight(1000)
+        self.script = ""
+
+    def clear(self):
+
+        toRm = []
+        for child in self.svgTree:
+            toRm.append(child)
+        for child in toRm:
+            self.svgTree.remove(child)
+
+        self.svgTree.append( self.getDefinitionBlock() )
+        self.nextFilterId = 0
+        self.nextObjId = 0
+
+    def getGroup(self, identifier, zIndex=0 ):
+        g = etree.Element('g')
+        g.set('id', str(identifier))
+        self.svgTree.append(g)
+        return(g)
+
+
+    def getTspan(self):
+        tspan = etree.Element('tspan')
+        tspan.set('id', str(self.nextTspanId))
+        self.nextTspanId+=1
+        return(tspan)
+
+    def addLegend(self,colorMapping):
+
+        y = 0
+        c = self.getYLabelCoord(y)
+        yp = c[1]+10 + 80
+        for color in colorMapping:
+
+
+            text = self.getText(str(colorMapping[color]), c[0]+self.plotStartX, yp,BDBcolor(80,80,80,1))
+            text.set('text-anchor','begin')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+            text.set('fill','%s' % color)
+            self.svgTree.append( text )
+            yp+=20
+
+
+
+    def getGroupColors(self, n):
+
+        if n==1:
+            return(['#3770C4'])
+        if n==2:
+            return(['#3770C4','#66A43E'])
+        if n==3:
+            return(['#3770C4','#66A43E','#F6853A'])
+        if n==4:
+            return(['#3770C4','#A43E3E','#66A43E','#F6853A'])
+        if n==5:
+            return(['#3770C4','#A43E3E','#66A43E','#F6853A','#A33DA2'])
+        if n==6:
+            return(['#3770C4','#A43E3E','#66A43E','#F6853A','#A33DA2','#AAD400'])
+        if n==7:
+            return(['#3770C4','#A43E3E','#66A43E','#F6853A','#A33DA2','#AAD400','#9DAC93'])
+        if n==8:
+            return(['#3770C4','#A43E3E','#66A43E','#F6853A','#A33DA2','#AAD400','#9DAC93','#7FCADF'])
+        if n==9:
+            return(['#3770C4','#A43E3E','#66A43E','#F6853A','#A33DA2','#AAD400','#9DAC93','#7FCADF','#D1AC17'])
+
+
+
+        return(['#3770C4','#A43E3E','#66A43E','#F6853A','#A33DA2','#AAD400','#9DAC93','#7FCADF','#D1AC17','#000080','#FF0066','#6C5D53'] + ['#333333']*n)
+
+
+    #Macro to set a title quickly
+    def setTitle(self, string, x=None, y=10, size=25, fill='#333333' ):
+
+        centerX = x is None
+        if x is None:
+            x,_= self.getPlottingCoord(self.xMin + 0.5*(self.xMax - self.xMin), 0)
+
+
+        text = self.getText(str(string), x,y, fill=fill)
+        if centerX:
+            text.set('text-anchor','middle')
+        else:
+            text.set('text-anchor','begin')
+        text.set('dominant-baseline','central')
+        text.set('font-family','Gill Sans MT')
+        text.set('font-size', str(size))
+        self.svgTree.append(text)
+        return(self)
+
+    def setSubtitle(self, string, x=None, y=40, size=15, fill='#222222'):
+        self.setTitle(string,x,y,size,fill)
+
+    def setWidth(self, width):
+        self.width = width
+        self.svgTree.set('width','%s' % width)
+
+    def setHeight(self, height):
+        self.height = height
+        self.svgTree.set('height','%s' % height)
+
+
+    def getDx(self):
+        return( float(self.plotWidth )/float((self.xMax - self.xMin)))
+
+
+    def getDy(self):
+        return( float(self.plotHeight )/float((self.yMax - self.yMin)))
+
+
+    def getPlottingCoord(self, x,y,z=0):
+
+        return( (self.plotStartX + (float(x)/(self.xMax - self.xMin))*self.plotWidth, self.plotHeight+self.plotStartY -(( float(y)/(self.yMax - self.yMin)))*self.plotHeight))
+
+
+    def getXLabelCoord(self, x):
+        return( (self.plotStartX + (float(x)/(self.xMax - self.xMin))*self.plotWidth, self.plotHeight+self.plotStartY+2 ))
+
+
+    def getYLabelCoord(self, y):
+        return( (self.plotStartX, self.plotHeight+self.plotStartY -(( float(y)/(self.yMax - self.yMin)))*self.plotHeight ))
+
+    def getNextObjId(self):
+        self.nextObjId+=1
+        return(str(self.nextObjId))
+
+
+    def getDefinitionBlock(self):
+
+        self.defs = etree.Element('defs')
+        self.defs.set('id','defs0')
+        return(self.defs)
+
+
+    def filter(self):
+        filterDef = etree.Element('filter')
+        filterDef.set('id','filter_%s' % (self.nextFilterId))
+        self.nextFilterId+=1
+        return(filterDef)
+
+    def getAxis(self,hv=0):
+
+        if hv==1:
+            p = self.getPath(self.getPathDefinition([self.getPlottingCoord(self.xMin, self.yMin), self.getPlottingCoord(self.xMax, self.yMin)]))
+        elif hv==2:
+            p = self.getPath(self.getPathDefinition([self.getPlottingCoord(self.xMin, self.yMax),self.getPlottingCoord(self.xMin, self.yMin)]))
+        else:
+            p = self.getPath(self.getPathDefinition([self.getPlottingCoord(self.xMin, self.yMax),self.getPlottingCoord(self.xMin, self.yMin), self.getPlottingCoord(self.xMax, self.yMin)]))
+        return(p)
+
+    def getPathDefinition(self, coordinates, preventAliasing=False ):
+
+        definition = []
+        for idx,coordinateTuple in enumerate(coordinates):
+
+            if preventAliasing:
+                coordinateTuple = ( round(coordinateTuple[0])+0.5, round(coordinateTuple[1])+0.5 )
+
+            if idx==0:
+                definition.append('M%s,%s' % (coordinateTuple[0],coordinateTuple[1]))
+            else:
+                definition.append('L%s,%s' % (coordinateTuple[0],coordinateTuple[1]))
+        return(' '.join(definition))
+
+
+    def getLinearGradientDefinition(self, tuplesWithStops): #Format: %x, color
+
+        definitionElement =  etree.Element('linearGradient')
+        definitionElement.set('id', self.getNextObjId()) #not really needed
+        definitionElement.set('x1', tuplesWithStops[0][0])
+        definitionElement.set('y1', tuplesWithStops[0][0])
+        definitionElement.set('x2', tuplesWithStops[-1][0])
+        definitionElement.set('y2', tuplesWithStops[-1][0])
+
+        for i, tup in enumerate(tuplesWithStops):
+
+            stop = etree.SubElement(definitionElement, 'stop')
+            stop.set('offset',tup[0])
+            stop.set('stop-color',tup[1])
+            if i==0:
+                stop.set('class','start')
+        stop.set('class','stop')
+        return(definitionElement)
+
+    def shadow(self, dy=2, dx=2, gaussStd=2,color='rgb(0,0,0)', floodOpacity=0.9,
+     width=None, # Width: set a pixel region around the filter to prevent clipping (https://stackoverflow.com/questions/17883655/svg-shadow-cut-off)
+     height=None):
+        f = self.filter()
+
+        if width is not None:
+            f.set('width', str(width))
+            f.set('x', '-%s' % (width*0.5))
+        if height is not None:
+            f.set('height', str(height))
+            f.set('y', '-%s' % (height*0.5))
+
+        f.set('color-interpolation-filters','sRGB')
+        self.nextFilterId+=1
+        #Flood
+        flood = etree.SubElement(f, 'feFlood')
+        flood.set('result','flood')
+        flood.set('flood-color',color)
+        flood.set('flood-opacity','%s' % floodOpacity)
+
+        #Composite filter
+        composite1 = etree.SubElement(f, 'feComposite')
+        composite1.set('in2','SourceGraphic')
+        composite1.set('operator','in')
+        composite1.set('in','flood')
+        composite1.set('result','composite1')
+
+        #Gaussian blur
+        gauss = etree.SubElement(f, 'feGaussianBlur')
+        gauss.set('result','blur')
+        gauss.set('stdDeviation','%s' % gaussStd)
+
+        #Shadow offset
+        offset = etree.SubElement(f, 'feOffset')
+        offset.set('result','offset')
+        offset.set('dy','%s' % dy)
+        offset.set('dx','%s' % dx)
+
+        #Final composite filter
+        composite2 = etree.SubElement(f, 'feComposite')
+        composite2.set('in2','offset')
+        composite2.set('operator','over')
+        composite2.set('in','SourceGraphic')
+        composite2.set('result','composite2')
+        return(f)
+
+
+    def makeInnerShadow(self, shadow ):
+
+        i = etree.SubElement(shadow, 'feComposite')
+        i.set('operator','in')
+        i.set('in2','SourceGraphic')
+
+
+
+    def addDef(self, filterDef, defId=None):
+        if defId is not None:
+            filterDef.set('id',defId)
+        self.defs.append(filterDef)
+        return( filterDef.get('id') )
+
+    def hasDef(self, defId):
+        return( len(self.defs.findall(".//*[@id='%s']" % defId))>0 )
+
+    def getDef(self,defId):
+        if not self.hasDef(defId):
+            print(('Definition %s was not found' % defId))
+            exit()
+        return( self.defs.findall(".//*[@id='%s']" % defId)[0] )
+
+    def warn(self, msg):
+        print(('[WARN] %s' % msg))
+
+    def getRectangle(self, x,y, width, height):
+        rectangle = etree.Element('rect')
+        rectangle.set('id', self.getNextObjId())
+        rectangle.set('x',str(x))
+        rectangle.set('y',str(y))
+        rectangle.set('width',str(width))
+        rectangle.set('height',str(height))
+        rectangle.set('style',"fill:rgba(100,100,100,1);stroke:#1b1b1b")
+        return(rectangle)
+
+    def getImage( self, path, x=None, y=None, width=None, height=None, preserveAspectRatio=None):
+        image = etree.Element('image')
+        image.set('id', self.getNextObjId())
+        if x is not None:
+            image.set('x',str(x))
+        if y is not None:
+            image.set('y',str(y))
+        if width is not None:
+            image.set('width',str(width))
+
+        if preserveAspectRatio is not None:
+            image.set('preserveAspectRatio',str(preserveAspectRatio))
+
+        if height is not None:
+            image.set('height',str(height))
+
+        image.set('{%s}href'% self.xlink ,str(path))
+
+        return(image)
+
+
+
+    #Modify attribute in style
+    def modifyStyleString(self, style, setAttr={},remove=[]):
+
+        attributes = {}
+
+        if style is not None and style.strip()!='':
+            parts = style.split(';')
+
+            for part in parts:
+                kvPair = part.split(':')
+                if len(kvPair)==2:
+                    key = kvPair[0]
+                    value = kvPair[1]
+
+                    if key not in remove:
+                        attributes[kvPair[0]] = kvPair[1]
+                else:
+                    self.warn('Style parsing %s failed (ignoring)' % part)
+
+
+
+            if self.debug>=3:
+                print('Style decomposition')
+                for key in attributes:
+                    print(('%s\t:\t%s' % (key, attributes[key]) ))
+
+        #Roll changes:
+        for attribute in setAttr:
+            attributes[attribute] = setAttr[attribute]
+
+        #Create new style string
+        newStyle = []
+        for attr in attributes:
+            newStyle.append('%s:%s' % (attr, attributes[attr]))
+
+        return(';'.join(newStyle))
+
+    def modifyStyle(self, element,setAttr={},remove=[]):
+
+        if not 'style' in element.attrib:
+            element.set('style','')
+
+        element.set('style', self.modifyStyleString(element.get('style'), setAttr, remove ))
+
+
+    def setTextRotation(self, element, angle ):
+        element.set('transform','rotate(%s, %s, %s)'%(angle,element.get('x'), element.get('y')))
+
+
+
+    def humanReadable(self, value, targetDigits=2,fp=0):
+
+        #Float:
+        if value<1 and value>0:
+            return('%.2f' % value )
+
+        if value == 0.0:
+            return('0')
+
+        baseId = int(math.floor( math.log10(float(value))/3.0 ))
+        suffix = ""
+        if baseId==0:
+            sVal =  str(round(value,targetDigits))
+            if len(sVal)>targetDigits and sVal.find('.'):
+                sVal = sVal.split('.')[0]
+
+        elif baseId>0:
+
+            sStrD = max(0,targetDigits-len(str( '{:.0f}'.format((value/(math.pow(10,baseId*3)))) )))
+
+
+            sVal = ('{:.%sf}' % min(fp, sStrD)).format((value/(math.pow(10,baseId*3))))
+            suffix = 'kMGTYZ'[baseId-1]
+        else:
+
+            sStrD = max(0,targetDigits-len(str( '{:.0f}'.format((value*(math.pow(10,-baseId*3)))) )))
+            sVal = ('{:.%sf}' %  min(fp, sStrD)).format((value*(math.pow(10,-baseId*3))))
+            suffix = 'mnpf'[-baseId-1]
+
+            if len(sVal)+1>targetDigits:
+                # :(
+                sVal = str(round(value,fp))[1:]
+                suffix = ''
+
+
+        return('%s%s' % (sVal,suffix))
+
+
+
+    def getText(self, text, x=0, y=0, fill=BDBcolor(), pathId=None):
+
+
+        textElement =  etree.Element('text')
+        textElement.set('id', self.getNextObjId())
+
+        if pathId != None:
+            tp =  etree.Element('textPath')
+            tp.text = text
+            tp.set('{%s}href'%self.xlink,  '#%s' % (pathId))
+            tp.set("startOffset", "50%")
+
+            textElement.append(tp)
+        else:
+            textElement.text = str(text)
+
+        textElement.set('x',str(x))
+        textElement.set('y',str(y))
+
+        if type(fill) is str:
+            textElement.set('fill',fill)
+        else:
+            textElement.set('fill',fill.getRGBStr())
+        textElement.set('shape-rendering','crispEdges')
+
+
+        return(textElement)
+
+    def getCenteredText(self,text, x,y, bold=False, fontSize=14, fill='rgba(50,50,50,1)', **kwargs):
+        text = self.getText(text,x,y,**kwargs)
+        text.set('text-anchor','middle')
+        text.set('dominant-baseline','middle')
+        text.set('font-family','Helvetica')
+        #text.set('font-family','Cambria Math')
+        if bold:
+            text.set('font-weight', 'bold')
+        text.set('font-size', str(fontSize))
+        text.set('fill', fill)
+        return text
+
+    def addTspan(self, textObject, text=None):
+        tspan = self.getTspan()
+        textObject.append(tspan)
+
+        #Fill with text if supplied:
+        if text is not None:
+            tspan.text = text
+
+        return(tspan)
+    #superscript
+    def addSuper(self, text, superText, offset=-10):
+
+        superElement =  etree.Element('tspan')
+        superElement.text = superText
+        superElement.set('dy', '%s' % offset)
+        text.append(superElement)
+
+
+
+
+    def polarToCartesian(self, centerX, centerY, radius, angleInDegrees = None, angleInRadians = None):
+
+        angleInRadians = (angleInDegrees-90) * math.pi / 180.0 if angleInDegrees is not None else angleInRadians
+        return({
+          'x': centerX + (radius * math.cos(angleInRadians)),
+          'y': centerY + (radius * math.sin(angleInRadians))
+        })
+
+    def describeArc(self, x, y, radius, startAngle, endAngle, sweep=0, largeArcFlag=None):
+
+        start = self.polarToCartesian(x, y, radius, endAngle)
+        end = self.polarToCartesian(x, y, radius, startAngle)
+
+        if largeArcFlag==None:
+            if endAngle - startAngle <= 180:
+                largeArcFlag  = "0"
+            else:
+                largeArcFlag  = "1"
+
+        d = " ".join([str(x) for x in [
+            "M", start['x'], start['y'],
+            "A", radius, radius, 0, largeArcFlag, sweep, end['x'], end['y']
+        ]])
+
+        return(d)
+
+    def describeArcRad(self, x, y, radius, startAngle, endAngle, sweep=0, largeArcFlag=None):
+
+        start = self.polarToCartesian(x, y, radius, angleInRadians=startAngle )
+        end = self.polarToCartesian(x, y, radius, angleInRadians=endAngle  )
+
+        if largeArcFlag==None:
+            if endAngle - startAngle <= math.pi:
+                largeArcFlag  = "0"
+            else:
+                largeArcFlag  = "1"
+
+        d = " ".join([str(x) for x in [
+            "M", start['x'], start['y'],
+            "A", radius, radius, 0, largeArcFlag, sweep, end['x'], end['y']
+        ]])
+
+        return(d)
+
+
+    def getCircle(self, centerX, centerY, radius):
+        circle = etree.Element('circle')
+        circle.set('id', self.getNextObjId())
+        circle.set('cx',str(centerX))
+        circle.set('cy',str(centerY))
+        circle.set('r',str(radius))
+        circle.set('style',"fill:none;stroke:#1b1b1b;stroke-width:1.29999995;stroke-linecap:round;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:5.2, 5.2;stroke-dashoffset:0")
+        return(circle)
+
+
+    def getPath(self, pathDef):
+        path = etree.Element('path')
+        path.set('id', self.getNextObjId())
+        path.set('d', pathDef)
+        path.set('style',"fill:none;stroke:#1b1b1b;stroke-width:1;stroke-linecap:round")
+        return(path)
+
+
+    def dump(self):
+        print(( etree.toString(self.root, pretty_print=True)))
+
+
+
+
+    def write(self, path, pretty=False, htmlCallback=None, bodyCallback=None):
+
+        try:
+            os.makedirs(os.path.dirname(path),exist_ok=True)
+        except:
+            pass
+
+        if len(self.script)>0:
+            html = etree.Element('html')
+
+
+            head = etree.Element('head')
+            body = etree.Element('body')
+
+            s = etree.Element('script')
+            s.set('type','text/javascript')
+            s.text = self.script
+
+            jquery = etree.Element('script')
+            jquery.set('src','https://ajax.googleapis.com/ajax/libs/jquery/2.2.3/jquery.min.js')
+            jquery.text = ' '
+            head.append(jquery)
+
+
+            jqueryUi = etree.Element('link')
+            jqueryUi.set('rel', 'stylesheet')
+            jqueryUi.set('href','https://ajax.googleapis.com/ajax/libs/jqueryui/1.12.1/themes/smoothness/jquery-ui.css')
+            jqueryUi.text = ' '
+            head.append(jqueryUi)
+
+            jqueryUi = etree.Element('script')
+            jqueryUi.set('src','https://ajax.googleapis.com/ajax/libs/jqueryui/1.12.1/jquery-ui.min.js')
+            jqueryUi.text = ' '
+            head.append(jqueryUi)
+
+            jqueryColor = etree.Element('script')
+            jqueryColor.set('src','http://code.jquery.com/color/jquery.color-2.1.2.js')
+            jqueryColor.text = ' '
+            head.append(jqueryColor)
+
+            body.append(self.svgTree)
+            html.append(head)
+
+            body.append(s)
+
+            if htmlCallback is not None:
+                htmlCallback(html)
+
+            if bodyCallback is not None:
+                bodyCallback(body)
+            html.append(body)
+
+            import html as pyhtml
+            if pretty:
+                import xml.dom.minidom as minidom
+                with open(path, 'w') as f:
+                    f.write( pyhtml.unescape(minidom.parseString(etree.tostring(html.getroottree())).toprettyxml(indent=" ").decode('utf-8') ))
+            else:
+
+                #html.getroottree().write(path)
+                with open(path, 'w') as f:
+                    f.write( pyhtml.unescape( etree.tostring(html.getroottree()).decode('utf-8') ) )
+        else:
+            try:
+                self.svgTree.getroottree().write(path)
+            except:
+                print("failed saving %s" % path)
+        return(path)
+
+    def SVGtoPNG(self,svgPath, pngPath, width=None, inkscapePath="C:\Program Files (x86)\Inkscape\inkscape.exe"):
+
+        if width is not None:
+            pass
+        else:
+            width = self.width
+
+
+        #cmd = '"%(INKSCAPE_PATH)s" -z --verb=org.ekips.filter.embedimage --verb=FileSave --verb=FileClose -f %(source_svg)s -w %(width)s -j -e %(dest_png)s' %  {'INKSCAPE_PATH':inkscapePath, 'source_svg':svgPath, 'dest_png':pngPath, 'width':width}
+        cmd = '"%(INKSCAPE_PATH)s" -z  --verb=FileSave --verb=FileClose -f %(source_svg)s -w %(width)s -e %(dest_png)s' %  {'INKSCAPE_PATH':inkscapePath, 'source_svg':svgPath, 'dest_png':pngPath, 'width':width}
+        os.system('%s' % cmd)
+        os.system('%s' % cmd)
+
+
+
+
+
+#circle = bdbplot.getCircle(200,200,50)
+#bdbplot.modifyStyle(circle, {'filter':'url(#%s)'%shadow.get('id')})
+#bdbplot.svgTree.append( circle )
+
+#circle = bdbplot.getCircle(250,250,50)
+#bdbplot.modifyStyle(circle, {'filter':'url(#%s)'%shadow.get('id')})
+#bdbplot.svgTree.append( circle )
+
+#path = bdbplot.getPath('M100 100 L300 100 L300 300 L300 100')
+#bdbplot.svgTree.append( path )
+
+
+##
+# Spaghettogram
+##
+
+
+
+##
+# Histogram
+##
+
+#
+# dictionary of read abundace->freq
+def readCountHistogram(abundanceFreqDict, logAbundance=True):
+
+    #We expect a distribution which is very steep.
+    lookAhead = 3
+
+
+    if logAbundance:
+        f = abundanceFreqDict
+        abundanceFreqDict = Counter({})
+        for a in f:
+            #print(("%s %s" % (a, f[a])))
+            try:
+                logA = int(round( math.log10(int(a)*100), 0))
+            except Exception as e:
+                logA = 0
+                print(e)
+            abundanceFreqDict[logA] += f[a]
+           # print(("%s %s -> %s %s" % (a, f[a], logA, abundanceFreqDict[logA] )))
+
+    #Find the highest abundant read:
+    hfreq = max([n for n in abundanceFreqDict])
+
+    #Find closed distribution:
+    #Mapping from abundance value to plotting X coordinate
+    xxMapping = {}
+
+    closedEnd = 1
+    perBin = 100
+    for c in range(1,perBin):
+
+        if abundanceFreqDict[c]==0 and 0==sum(abundanceFreqDict[q] for q in range(c,c+lookAhead+1)):
+            closedEnd = c-1
+            break
+
+        else:
+            xxMapping[c] = c-0.5
+
+    #check how many extra blocks we need
+
+    extraBlocks = 0
+    prevX = closedEnd-0.5
+    maxX = 1
+    prevAbundance = closedEnd
+    extraBlockCoords = []
+    extraBlockContinuous = {}
+    for abundance in sorted(abundanceFreqDict.keys()):
+        if abundance > closedEnd:
+
+
+            if (abundance-prevAbundance) > 1:
+                xxMapping[abundance] = prevX+2
+                extraBlockContinuous[extraBlocks] = True
+            else:
+                xxMapping[abundance] = prevX+1
+                extraBlockContinuous[extraBlocks] = False
+
+            prevAbundance=abundance
+            extraBlockCoords.append(abundance)
+            prevX = xxMapping[abundance]
+            maxX= xxMapping[abundance]+1
+            extraBlocks+=1
+    #print(xxMapping)
+    bdbplot = BDBPlot()
+    bdbplot.plotStartX = 100
+    bdbplot.plotStartY = 150
+
+    bdbplot.plotHeight =400
+    bdbplot.plotWidth = max(600, (maxX) * 25)
+
+    bdbplot.setWidth(bdbplot.plotWidth+bdbplot.plotStartX+10)
+    bdbplot.setHeight(800)
+
+
+    bdbplot.xMax = max(1,maxX) # prevent 0 (breaks everything, 0 divisions and such)
+    bdbplot.yMax = max(1,int(math.log10(hfreq)))
+
+    axis = bdbplot.getAxis(2)
+    bdbplot.svgTree.append( axis )
+
+    #Draw specialised x-axis
+
+    p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(bdbplot.xMin, bdbplot.yMin),bdbplot.getPlottingCoord(xxMapping[closedEnd]+1, bdbplot.yMin)]))
+    bdbplot.svgTree.append( p )
+
+    for extraBlock in range(0,extraBlocks):
+        x = xxMapping[extraBlockCoords[extraBlock]]
+
+        #p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(x, bdbplot.yMin),bdbplot.getPlottingCoord(x+1, bdbplot.yMin)]))
+
+        if  extraBlockContinuous[extraBlock]:
+            d = 0.15
+            p = bdbplot.getPath(bdbplot.getPathDefinition([
+                bdbplot.getPlottingCoord(x-1, bdbplot.yMin),
+                bdbplot.getPlottingCoord(x-0.75, bdbplot.yMin+d),
+                bdbplot.getPlottingCoord(x-0.25, bdbplot.yMin-d),
+                bdbplot.getPlottingCoord(x, bdbplot.yMin)
+                ]))
+
+            bdbplot.modifyStyle(p, {'stroke-width':'1', 'stroke-linecap':'round', 'stroke-dasharray':'2 2','stroke-dashoffset':'0'} )
+            bdbplot.svgTree.append( p )
+
+
+        p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(x, bdbplot.yMin),bdbplot.getPlottingCoord(x+1, bdbplot.yMin)]))
+        bdbplot.svgTree.append( p )
+
+
+
+        #p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(x+2, bdbplot.yMin),bdbplot.getPlottingCoord(x+3, bdbplot.yMin)]))
+        #bdbplot.modifyStyle(p, {'stroke-width':'1', 'stroke-linecap':'round', 'stroke-dasharray':'2 2','stroke-dashoffset':'0'} )
+        #bdbplot.svgTree.append( p )
+
+
+    #Draw fine grid
+    for y in range(1,bdbplot.yMax+1):
+        p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(bdbplot.xMin, y),bdbplot.getPlottingCoord(bdbplot.xMax, y)]))
+
+        bdbplot.modifyStyle(p, {'stroke-width':'0.5', 'stroke-linecap':'round', 'stroke-dasharray':'2 2','stroke-dashoffset':'0'} )
+        bdbplot.svgTree.append( p )
+
+
+    ### Block plotting ###
+    rectangles = []
+    barShadow = bdbplot.shadow(1,1)
+    bdbplot.addDef(barShadow)
+
+    for abundance in range(1,int(hfreq)+1):
+
+        if abundanceFreqDict[abundance]>0:
+            plotX = xxMapping[abundance]
+            frequency = abundanceFreqDict[abundance]
+
+            if frequency==1:
+                value=0.20
+            else:
+                value = math.log10(frequency)
+
+            c = bdbplot.getPlottingCoord(plotX,value)
+            origin = bdbplot.getPlottingCoord(plotX,0)
+
+            barWidth = float(bdbplot.plotWidth)/(bdbplot.xMax+1)
+
+            rectangleParams = (c[0], c[1], barWidth,  (float(value)/bdbplot.yMax) * bdbplot.plotHeight-3)
+            rectangles.append(rectangleParams)
+            bar = bdbplot.getRectangle( *rectangleParams )
+            bdbplot.modifyStyle(bar, {'filter':'url(#%s)'%barShadow.get('id'),'fill':'rgba(255,255,255,1)'})
+            bdbplot.svgTree.append( bar )
+
+            text = bdbplot.getText(str( bdbplot.humanReadable(frequency,1 ) ), c[0]+0.5*barWidth, c[1]-10,BDBcolor(0,0,0,1))
+            text.set('text-anchor','middle')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Gill Sans MT')
+            #text.set('font-family','Cambria Math')
+            text.set('font-size', '14')
+            text.set('font-weight', 'bold')
+            text.set('fill', 'rgba(50,50,50,1)')
+            bdbplot.svgTree.append( text )
+
+            #AXIS LABEL
+            c = bdbplot.getXLabelCoord(plotX+0.5)
+            text = bdbplot.getText(str(bdbplot.humanReadable(abundance)), c[0], c[1]+15,BDBcolor(80,80,80,1))
+            text.set('text-anchor','middle')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+            bdbplot.svgTree.append( text )
+
+
+
+
+    for rect in rectangles:
+        bar = bdbplot.getRectangle( *rect )
+        bdbplot.modifyStyle(bar, {'fill':bdbplot.getGroupColors(1)[0], 'stroke':'#FFFFFF','stroke-width':'1.5'})
+        bdbplot.svgTree.append( bar )
+
+
+
+    #Y axis label
+    for y in range(1,bdbplot.yMax+1):
+        c = bdbplot.getYLabelCoord(y)
+
+        value = math.pow(10,y)
+
+        text = bdbplot.getText(str(10), c[0]-10, c[1],BDBcolor(80,80,80,1))
+        text.set('text-anchor','end')
+        text.set('dominant-baseline','middle')
+        text.set('font-family','Cambria Math')
+        bdbplot.addSuper(text,str(y))
+        bdbplot.svgTree.append( text )
+
+
+
+    c = bdbplot.getYLabelCoord( (bdbplot.yMax/2))
+    text = bdbplot.getText('Frequency', c[0]-60, c[1]-30,BDBcolor(0,0,0,1))
+    text.set('text-anchor','middle')
+    text.set('dominant-baseline','middle')
+    text.set('font-family','Gill Sans MT')
+    text.set('font-size', '25')
+    #bdbplot.modifyStyle(text, {'font-size': '20'})
+    bdbplot.setTextRotation(text,270)
+    bdbplot.svgTree.append( text )
+
+    c = bdbplot.getXLabelCoord(bdbplot.xMax/2)
+    text = bdbplot.getText('Read abundance', c[0], c[1]+50,BDBcolor(0,0,0,1))
+    text.set('text-anchor','middle')
+    text.set('dominant-baseline','middle')
+    text.set('font-family','Gill Sans MT')
+    text.set('font-size', '25')
+    bdbplot.svgTree.append( text )
+
+
+    return(bdbplot)
+
+
+
+
+class subdividedHistClass():
+
+    def __init__(self, name, dataPoints, logTransform=False, offset=0):
+        self.logTransform = logTransform
+        self.totalValue = 0
+        self.barSpacerWidth = 5
+        self.barWidth = 40
+        self.maxValue = 0
+        self.bars = []
+        self.name = name
+        self.startX = offset
+        x = offset
+        self.width = 0
+        for dName,count in dataPoints.most_common():
+            self.totalValue+=count
+            x+=self.barSpacerWidth
+            self.bars.append({'x':x, 'y':count,'name':dName})
+            x+=self.barWidth
+
+            self.maxValue = max(self.maxValue, count)
+
+        self.width = x + self.barSpacerWidth - offset
+
+    def plot(self, bdbplot, scarAliases,subClassColors ):
+        #Draw full class rectangle:
+        value = 1
+        if self.logTransform and self.totalValue>0:
+            value = math.log(self.totalValue)
+
+        if not self.logTransform:
+            value = self.totalValue
+
+        c = bdbplot.getPlottingCoord(self.startX,value)
+        origin = bdbplot.getPlottingCoord(self.startX,0)
+
+        rectangleParams = (c[0], c[1], self.width,  (float(value)/bdbplot.yMax) * bdbplot.plotHeight)
+        bar = bdbplot.getRectangle( *rectangleParams )
+        bdbplot.modifyStyle(bar, {'fill':'rgba(150,150,150,0.8)','stroke-width':'0'})
+        bdbplot.svgTree.append( bar )
+
+
+        text = bdbplot.getText(self.name, c[0]+0.5*self.width, c[1]- 10,BDBcolor(0,0,0,1))
+        text.set('text-anchor','middle')
+        text.set('dominant-baseline','middle')
+        text.set('font-family','Gill Sans MT')
+        #text.set('font-family','Cambria Math')
+        text.set('font-size', '14')
+        text.set('font-weight', 'bold')
+        text.set('fill', 'rgba(50,50,50,1)')
+        bdbplot.svgTree.append( text )
+
+        barShadow = bdbplot.shadow(1,1)
+        bdbplot.addDef(barShadow)
+        rectangles= []
+        for bar in self.bars:
+            #Add bar:
+            plotX = bar['x']
+            frequency= bar['y']
+            className = scarAliases[ bar['name'] ]
+            barColor = subClassColors[bar['name']]
+
+            #Add class label to X axis:
+            c = bdbplot.getXLabelCoord(plotX+self.barWidth*0.5)
+            text = bdbplot.getText(className, c[0], c[1]+15,BDBcolor(80,80,80,1))
+            text.set('text-anchor','middle')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+            bdbplot.svgTree.append( text )
+
+            if self.logTransform:
+                if frequency==1:
+                    value=0.20
+                else:
+                    value = math.log10(frequency)
+            else:
+                value = frequency
+
+            c = bdbplot.getPlottingCoord(plotX,value)
+            origin = bdbplot.getPlottingCoord(plotX,0)
+
+            #barWidth = float(bdbplot.plotWidth)/(bdbplot.xMax+1)
+            barWidth = self.barWidth
+            rectangleParams = (c[0], c[1], barWidth,  (float(value)/bdbplot.yMax) * bdbplot.plotHeight)
+            rectangles.append(rectangleParams)
+            bar = bdbplot.getRectangle( *rectangleParams )
+            bdbplot.modifyStyle(bar, {'stroke-width':'0','filter':'url(#%s)'%barShadow.get('id'),'fill':barColor})
+            bdbplot.svgTree.append( bar )
+
+            text = bdbplot.getText(str( bdbplot.humanReadable(frequency,3,2 ) ), c[0]+0.5*barWidth, c[1]-10,BDBcolor(0,0,0,1))
+            text.set('text-anchor','middle')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Gill Sans MT')
+            #text.set('font-family','Cambria Math')
+            text.set('font-size', '14')
+            text.set('font-weight', 'bold')
+            text.set('fill', 'rgba(50,50,50,1)')
+
+            bdbplot.svgTree.append( text )
+
+            #Percentile:
+            text = bdbplot.getText(str( bdbplot.humanReadable(100*(float(frequency)/self.totalValue),3,2 )+'%' ), c[0]+0.5*barWidth, c[1]+10,BDBcolor(255,255,255,1))
+            text.set('text-anchor','middle')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Gill Sans MT')
+            #text.set('font-family','Cambria Math')
+            text.set('font-size', '14')
+            text.set('font-weight', 'bold')
+            text.set('fill', 'rgba(255,255,255,1)')
+
+            bdbplot.svgTree.append( text )
+
+
+
+def subdividedClassHistogram( classes, logTransform= False, scarAliases={}):
+
+    classSpacerWidth = 20
+
+    currentX = classSpacerWidth
+    classIndex = 0
+    maxValue = 0
+
+    subHistClassList = []
+
+    for className in classes:
+        shc = subdividedHistClass(className,classes[className], logTransform, currentX)
+        currentX += shc.width + classSpacerWidth
+
+        if logTransform:
+            if shc.totalValue>0:
+                maxValue = max(maxValue,math.log(shc.totalValue))
+        else:
+            maxValue = max(maxValue,shc.totalValue)
+
+        subHistClassList.append( shc )
+
+
+    bdbplot = BDBPlot()
+
+    ## color list:
+    subClassColors = {}
+    idx = 0
+    #print(len(scarAliases))
+    gc = bdbplot.getGroupColors( len(scarAliases) )
+    for scar in scarAliases:
+        subClassColors[scar] = gc[idx]
+        print(('%s -> %s' %(scar, gc[idx])))
+        idx+=1
+
+    ## Plot area preparation
+    bdbplot.plotStartX = 100
+    bdbplot.plotStartY = 100
+
+    bdbplot.plotHeight =400
+    bdbplot.plotWidth = currentX
+
+    bdbplot.setWidth(bdbplot.plotWidth+bdbplot.plotStartX+10)
+    bdbplot.setHeight(800)
+    bdbplot.xMax = max(1,currentX) # prevent 0 (breaks everything, 0 divisions and such)
+    if logTransform:
+        bdbplot.yMax = max(1,int(maxValue)+1)
+    else:
+        bdbplot.yMax = max(1,int(maxValue+1))
+
+    axis = bdbplot.getAxis()
+    bdbplot.svgTree.append( axis )
+    if logTransform:
+        for y in range(1,bdbplot.yMax+1):
+            p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(bdbplot.xMin, y),bdbplot.getPlottingCoord(bdbplot.xMax, y)], True))
+
+            bdbplot.modifyStyle(p, {'stroke-width':'0.5', 'stroke-linecap':'round', 'stroke-dasharray':'2 2','stroke-dashoffset':'0'} )
+            bdbplot.svgTree.append( p )
+
+    #Draw fine grid
+    if logTransform:
+        for y in range(1,bdbplot.yMax+1):
+            p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(bdbplot.xMin, y),bdbplot.getPlottingCoord(bdbplot.xMax, y)], True))
+
+            bdbplot.modifyStyle(p, {'stroke-width':'0.5', 'stroke-linecap':'round', 'stroke-dasharray':'2 2','stroke-dashoffset':'0'} )
+            bdbplot.svgTree.append( p )
+
+    else:
+        stepSize = 50000
+        if bdbplot.yMax<101:
+            stepSize=10
+        for y in range(0,bdbplot.yMax+1,stepSize):
+            p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(bdbplot.xMin, y),bdbplot.getPlottingCoord(bdbplot.xMax, y)], True))
+
+            bdbplot.modifyStyle(p, {'stroke-width':'0.5', 'stroke-linecap':'round', 'stroke-dasharray':'2 2','stroke-dashoffset':'0'} )
+            bdbplot.svgTree.append( p )
+
+            c = bdbplot.getYLabelCoord(y)
+
+            text = bdbplot.getText(bdbplot.humanReadable(y,2,3 ), c[0]-10, c[1],BDBcolor(80,80,80,1))
+            text.set('text-anchor','end')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+            bdbplot.svgTree.append( text )
+
+
+    if logTransform:
+        for y in range(1,bdbplot.yMax+1):
+            c = bdbplot.getYLabelCoord(y)
+
+
+            value = math.pow(10,y)
+
+            text = bdbplot.getText(str(10), c[0]-10, c[1],BDBcolor(80,80,80,1))
+            text.set('text-anchor','end')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+            bdbplot.addSuper(text,str(y))
+            bdbplot.svgTree.append( text )
+
+
+    for shc in subHistClassList:
+        shc.plot(bdbplot,scarAliases,subClassColors)
+
+
+    c = bdbplot.getYLabelCoord( (bdbplot.yMax/2))
+    text = bdbplot.getText('Reads', c[0]-60, c[1]-30,BDBcolor(0,0,0,1))
+    text.set('text-anchor','middle')
+    text.set('dominant-baseline','middle')
+    text.set('font-family','Gill Sans MT')
+    text.set('font-size', '25')
+    #bdbplot.modifyStyle(text, {'font-size': '20'})
+    bdbplot.setTextRotation(text,270)
+    bdbplot.svgTree.append( text )
+
+    c = bdbplot.getXLabelCoord(bdbplot.xMax/2)
+    text = bdbplot.getText('Sample', c[0], c[1]+50,BDBcolor(0,0,0,1))
+    text.set('text-anchor','middle')
+    text.set('dominant-baseline','middle')
+    text.set('font-family','Gill Sans MT')
+    text.set('font-size', '25')
+    bdbplot.svgTree.append( text )
+
+
+    return({'plot':bdbplot,'colorMapping':subClassColors})
+
+
+
+def classHistogram( classCountMapping, logTransform = False, classColors=None,  placeLeft=None,placeRight=None, reverseOrder=True, height=400, xLabel='Sample', yLabel='Reads', yStepper = 50000, classWidth=50, freqFontSize=14, xLabelFontSize=10,  rotateClassLabels=0,zebraFillMode=False, defaultFillColor='#404040', barSpacing=5, xLabelOffset=50, showZeros=False, axisLabelFontSize=25, drawFreqLabels=True, title=None, freqMethod='humanReadable' ): # freqmethod 'humanReadable' ,'float'
+
+    amountOfClasses = len(classCountMapping)
+    #classWidth = 50
+    maxValue = 0
+    for className in classCountMapping:
+        maxValue = max(classCountMapping[className],maxValue)
+
+
+    bdbplot = BDBPlot()
+    bdbplot.plotStartX = 100
+    bdbplot.plotStartY = 100
+
+    bdbplot.plotHeight = height
+    bdbplot.plotWidth = max(600, (amountOfClasses) * classWidth)
+
+    bdbplot.setWidth(bdbplot.plotWidth+bdbplot.plotStartX+10)
+    bdbplot.setHeight(700)
+
+    if title is not None:
+        text = bdbplot.getText(title, 10,20, fill='#666666')
+        text.set('text-anchor','begin')
+        text.set('dominant-baseline','central')
+        text.set('font-family','Gill Sans MT')
+        text.set('font-size', '25')
+        bdbplot.svgTree.append(text)
+
+
+    bdbplot.xMax = max(1,amountOfClasses) # prevent 0 (breaks everything, 0 divisions and such)
+
+    if logTransform:
+        bdbplot.yMax = max(1,int(math.log10(maxValue)+1))
+    else:
+        bdbplot.yMax = max(1,int(maxValue+1))
+
+    axis = bdbplot.getAxis(2)
+    bdbplot.svgTree.append( axis )
+
+    classIndex = 0
+
+    rectangles = []
+    barShadow = bdbplot.shadow(1,1)
+    bdbplot.addDef(barShadow)
+
+    whiteShadow =  bdbplot.shadow(1,1,3,'rgb(255,255,255)')
+    bdbplot.addDef(whiteShadow)
+
+    #Draw fine grid
+    if logTransform:
+        for y in range(1,bdbplot.yMax+1):
+            p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(bdbplot.xMin, y),bdbplot.getPlottingCoord(bdbplot.xMax, y)], True))
+
+            bdbplot.modifyStyle(p, {'stroke-width':'0.5', 'stroke-linecap':'round', 'stroke-dasharray':'2 2','stroke-dashoffset':'0'} )
+            bdbplot.svgTree.append( p )
+
+    else:
+        for y in np.arange(0,bdbplot.yMax+yStepper,yStepper):
+            p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(bdbplot.xMin, y),bdbplot.getPlottingCoord(bdbplot.xMax, y)], True))
+
+            bdbplot.modifyStyle(p, {'stroke-width':'0.5', 'stroke-linecap':'round', 'stroke-dasharray':'2 2','stroke-dashoffset':'0'} )
+            bdbplot.svgTree.append( p )
+
+            c = bdbplot.getYLabelCoord(y)
+            if(freqMethod=='humanReadable'):
+                text = bdbplot.getText(bdbplot.humanReadable(y,2,3 ), c[0]-10, c[1],BDBcolor(80,80,80,1))
+            else:
+                text = bdbplot.getText(float(y), c[0]-10, c[1],BDBcolor(80,80,80,1))
+            text.set('text-anchor','end')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+
+            bdbplot.svgTree.append( text )
+
+
+
+    if isinstance(classCountMapping, collections.OrderedDict):
+
+        if reverseOrder:
+            classOrderKeys = list(reversed(list(classCountMapping.keys())))
+        else:
+            classOrderKeys = list(classCountMapping.keys())
+
+        classOrder = [
+            (key, classCountMapping[key]) for key in classOrderKeys
+        ]
+
+    else:
+        if reverseOrder:
+            classOrder = list(reversed(classCountMapping.most_common()))
+        else:
+            classOrder = list(classCountMapping.most_common())
+    classOrderKeys = {}
+    for className,freq in classOrder:
+        classOrderKeys[className] = (className, freq)
+    #Prepend left desired classes to the left
+    if placeLeft is not None:
+        for className in placeLeft:
+            #print(className)
+            tup = classOrderKeys[className]
+            classOrder.remove(tup)
+            classOrder.insert(0, tup)
+    if placeRight is not None:
+        for className in placeRight:
+            #print(className)
+            tup = classOrderKeys[className]
+            classOrder.remove(tup)
+            classOrder.append( tup)
+
+    barWidth = (float(bdbplot.plotWidth)/(bdbplot.xMax+1))
+    for className, frequency in classOrder:
+
+        #Add class label to X axis:
+        c = bdbplot.getXLabelCoord(classIndex)
+        text = bdbplot.getText(str(className), c[0] + ( 0.5*(barWidth)), c[1]+15,BDBcolor(80,80,80,1))
+        text.set('text-anchor','middle')
+        text.set('dominant-baseline','middle')
+        text.set('font-family','Cambria Math')
+        text.set('font-size',str(xLabelFontSize))
+
+        if rotateClassLabels!=0:
+            bdbplot.setTextRotation(text, rotateClassLabels)
+            text.set('text-anchor','start')
+        if rotateClassLabels>90:
+            bdbplot.setTextRotation(text, rotateClassLabels)
+            text.set('text-anchor','end')
+
+
+        bdbplot.svgTree.append( text )
+
+        #Add bar:
+        plotX = classIndex
+
+        if logTransform:
+
+
+            if frequency==1:
+                value=0.20
+            elif frequency==0:
+                value= 0
+            else:
+                value = math.log10(frequency)
+        else:
+            value = frequency
+
+        c = bdbplot.getPlottingCoord(plotX,value)
+        origin = bdbplot.getPlottingCoord(plotX,0)
+
+
+
+
+        #barWidth = float(bdbplot.plotWidth)/(bdbplot.xMax+1)
+
+
+        rectangleParams = (c[0]+0.5*barSpacing, c[1], barWidth-barSpacing,  (float(value)/bdbplot.yMax) * bdbplot.plotHeight-3)
+        rectangles.append(rectangleParams)
+        bar = bdbplot.getRectangle( *rectangleParams )
+        bdbplot.modifyStyle(bar, {'filter':'url(#%s)'%barShadow.get('id'),'fill':'rgba(255,255,255,1)'})
+        bdbplot.svgTree.append( bar )
+
+        if (showZeros or frequency>0) and drawFreqLabels:
+            if(freqMethod=='humanReadable'):
+                text = bdbplot.getText(str( bdbplot.humanReadable(frequency,2,1 ) ), c[0]+0.5*barWidth, c[1]-10,BDBcolor(0,0,0,1))
+            else:
+                text = bdbplot.getText('%.2f' % frequency  , c[0]+0.5*barWidth, c[1]-10,BDBcolor(0,0,0,1))
+            text.set('text-anchor','middle')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Gill Sans MT')
+            #text.set('font-family','Cambria Math')
+            text.set('font-size', str(freqFontSize))
+            text.set('font-weight', 'bold')
+            text.set('fill', 'rgba(50,50,50,1)')
+            bdbplot.modifyStyle(text, {'filter':'url(#%s)'%whiteShadow.get('id')})
+
+            bdbplot.svgTree.append( text )
+
+        classIndex += 1
+
+
+    if logTransform:
+        for y in range(1,bdbplot.yMax+1):
+            c = bdbplot.getYLabelCoord(y)
+
+
+            value = math.pow(10,y)
+
+            text = bdbplot.getText(str(10), c[0]-10, c[1],BDBcolor(80,80,80,1))
+            text.set('text-anchor','end')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+            bdbplot.addSuper(text,str(y))
+            bdbplot.svgTree.append( text )
+
+
+    for idx,rect in enumerate(rectangles):
+        bar = bdbplot.getRectangle( *rect )
+
+        if zebraFillMode:
+            fillColor = bdbplot.getGroupColors(3)[idx%2]
+        else:
+            fillColor = defaultFillColor
+        if classColors is not None:
+            if classOrder[idx][0] in classColors:
+                fillColor = classColors[classOrder[idx][0]]
+            else:
+                #print(('Setting %s to default color' % classOrder[idx][0]))
+                pass
+
+        bdbplot.modifyStyle(bar, {'fill':fillColor, 'stroke':'#FFFFFF','stroke-width':'1.75'})
+        bdbplot.svgTree.append( bar )
+
+
+
+
+    c = bdbplot.getYLabelCoord( (bdbplot.yMax/2))
+    text = bdbplot.getText(yLabel, c[0]-60, c[1]-30,BDBcolor(0,0,0,1))
+    text.set('text-anchor','middle')
+    text.set('dominant-baseline','middle')
+    text.set('font-family','Gill Sans MT')
+    text.set('font-size', str(axisLabelFontSize))
+    #bdbplot.modifyStyle(text, {'font-size': '20'})
+    bdbplot.setTextRotation(text,270)
+    bdbplot.svgTree.append( text )
+
+    c = bdbplot.getXLabelCoord(bdbplot.xMax/2)
+    text = bdbplot.getText(xLabel, c[0], c[1]+xLabelOffset,BDBcolor(0,0,0,1))
+    text.set('text-anchor','middle')
+    text.set('dominant-baseline','middle')
+    text.set('font-family','Gill Sans MT')
+    text.set('font-size', str(axisLabelFontSize))
+    bdbplot.svgTree.append( text )
+
+
+    return(bdbplot)
+
+
+
+class Heatmap(object):
+
+    def __init__(self, npMatrix, colorMatrix=None, rowNames=None,  rowColors=None, columnColors=None, cellFormat=None, cellIdentifiers=None, cellRotations=None, cellStrings=None, columnNames=None, cellSize=25, title=None, subtitle=None, nominalColoring=False, rotateColumnLabels=90, cellAnnot=None, cellAnnotFormat=None, metaDataMatrix=None, cluster=False, groupSize=10 ):
+
+        self.nominalColoring = nominalColoring
+        self.nominalColoringMapping = None
+        self.zeroColor = None
+        self.colormap = matplotlib.cm.get_cmap('plasma')
+        self.NanColor = (0.9,0.9,0.9,1)
+        self.rotateColumnLabels = rotateColumnLabels
+        self.footerHeight = 400
+        self.groupSpacerSize = 10
+        self.cellFont = 'Cambria'
+        self.labelFont = 'Cambria'
+
+        print("Plotting %s by %s matrix" % npMatrix.shape)
+        if rowNames is not None:
+            print("Supplied %s rownames" % len(rowNames))
+            rowNames = [ "%s: %s"%t for t in enumerate(rowNames) ]
+
+        if columnNames is not None:
+            print("Supplied %s column names" % len(columnNames))
+
+        if cluster:
+            print("Clustering")
+            self.matrix = npMatrix
+            clusterMatrix = np.zeros( npMatrix.shape )
+            for (column,row), value in np.ndenumerate(npMatrix):
+                l = self.nominalIndex(value)
+                clusterMatrix[column,row] = l
+
+            distances = scipy.spatial.distance.pdist( np.nan_to_num(clusterMatrix.transpose()), 'cityblock' )
+            mdistMatrix = scipy.spatial.distance.squareform(distances)
+            clustering = scipy.cluster.hierarchy.linkage( mdistMatrix, 'ward' )
+            leavesList = list( scipy.cluster.hierarchy.leaves_list(clustering) )
+            #npMatrix  = clusterMatrix
+            print(leavesList)
+        else:
+            leavesList = list(range(npMatrix.shape[1]))
+        print("Ranging %s" % len(leavesList))
+
+        self.matrix = npMatrix[:,leavesList]
+        if colorMatrix is not None:
+            self.colorMatrix = colorMatrix[:,leavesList] #Values between zero and one
+        else:
+            self.colorMatrix = None
+        self.rowColors = np.array(rowColors)[leavesList] if rowColors is not None else []
+        self.columnColors = columnColors if columnColors is not None else []
+        self.cellIdentifiers = cellIdentifiers if cellIdentifiers is not None else None
+        self.cellAnnot = cellAnnot
+        self.rowNames =np.array(rowNames)[leavesList] if rowNames is not None else []
+        self.columnNames = columnNames if columnNames is not None else []
+        self.cellSize = cellSize
+        self.cellFormat = cellFormat if cellFormat is not None else lambda x: x
+        self.cellAnnotFormat  = cellAnnotFormat if cellAnnotFormat is not None else lambda x: x
+        self.cellStrings = cellStrings[:,leavesList] if cellStrings is not None else []
+        self.cellRotations = cellRotations[:,leavesList] if cellRotations is not None else None
+
+        self.metaDataMatrix = metaDataMatrix[:,leavesList] if metaDataMatrix is not None else None
+        self.title = title
+        self.subtitle = subtitle
+        self.leftMargin = 80
+        self.labelWidth = 20
+        self.topMargin = 150
+        self.cellSpacing = 1
+        self.cellFontSize = 10
+        self.labelFontSize = 15
+        self.groupSize = groupSize
+        #self.colormap = matplotlib.cm.get_cmap('inferno')
+        print(self.rowNames)
+
+
+    def getRowName(self, index):
+        try:
+            return(str(self.rowNames[index]))
+        except:
+            return('')
+    def getColName(self, index):
+        try:
+            return(str(self.columnNames[index]))
+        except:
+            return('')
+
+    def getRowColor(self,index):
+        try:
+            c = self.rowColors[index]
+
+            return( c )
+        except:
+            return( BDBcolor( 50,50,50, 1) )
+
+    def getColumnColor(self,index):
+        try:
+            c = self.columnColors[index]
+            return( c )
+        except:
+            return( BDBcolor( 50,50,50, 1) )
+
+    def getCellString(self, row,column):
+        try:
+            return( self.cellFormat(self.cellStrings[column,row]))
+        except:
+            return('')
+
+    def getCellAnnotString(self, row,column):
+        try:
+            return( str(self.cellAnnotFormat(self.cellAnnot[column, row])))
+        except:
+            return('')
+
+    def getCellId(self, row,column):
+        try:
+            return( str(self.cellIdentifiers[column, row]) )
+        except:
+            return(None)
+    def getColor(self, value):
+
+        theValueIsNan = self.isnan(value)
+
+        if theValueIsNan:
+            r,g,b,a  = self.NanColor
+        elif self.zeroColor is not None and value==0:
+            r,g,b,a  = self.zeroColor
+        else:
+
+            if self.nominalColoring:
+                r,g,b,a = self.nominalColor(value)
+            else:
+                try:
+                    r,g,b,a = self.colormap(value)
+                except:
+
+                    print("Reverted to nominal coloring mode")
+                    self.nominalColoring = True
+                    r,g,b,a = self.nominalColor(value)
+        return(r,g,b,a)
+
+
+    def nominalIndex(self,value):
+        #Force build of nominal matrix
+        r,g,b,a = self.nominalColor(value)
+        #get index:
+        try:
+            idx = list(self.nominalColoringMapping.keys()).index(value)
+        except:
+            idx=0.01
+        return(idx / len(self.nominalColoringMapping.keys()))
+
+
+    def getCellMetaData(self, row,column):
+        try:
+            return( str(self.metaDataMatrix[column, row]) )
+        except:
+            return(None)
+
+    def addTitle(self, plot):
+        if self.title is not None:
+            text = plot.getText(self.title, 10,20, fill='#666666')
+            text.set('text-anchor','begin')
+            text.set('dominant-baseline','central')
+            text.set('font-family','Gill Sans MT')
+            text.set('font-size', '25')
+            plot.svgTree.append(text)
+
+    def getCellRotation(self, row, column):
+        try:
+            if self.cellRotations[column, row] == np.nan:
+                return(None)
+
+            return( str(self.cellRotations[column, row]) )
+        except:
+            return(None)
+
+
+    def addSubtitle(self, plot):
+        if self.subtitle is not None:
+            text = plot.getText(self.subtitle, 10,40, fill='#222222')
+            text.set('text-anchor','begin')
+            text.set('dominant-baseline','central')
+            text.set('font-family','Cambria')
+            text.set('font-size', '15')
+            plot.svgTree.append(text)
+
+    def addGroupTitle(self, plot, title, indexStart, indexEnd):
+
+        matrixGroup = plot.svgTree.findall(".//g[@id='matrix']")[0]
+        #Calculate x-starting coordinate
+        xStart = self.columnIndexToXCoordinate(indexStart)
+        #Calculate x-ending coordinate
+        xEnd = self.columnIndexToXCoordinate(indexEnd)+self.cellSize
+        yStart = self.rowIndexToYCoordinate(0)-self.cellSize
+        yEnd = yStart - self.cellSize*0.5
+
+        p = plot.getPath(plot.getPathDefinition([(xStart, yStart),(xStart,yEnd),(xEnd, yEnd), (xEnd, yStart)]))
+        plot.modifyStyle(p, {'stroke-width':'1', 'stroke-linecap':'round','stroke-dashoffset':'0', 'stroke':'#333333'} )
+        matrixGroup.append( p )
+
+        text = plot.getText(title, xStart+0.5*( (xEnd-xStart)), yEnd - 8, fill='#000000')
+        text.set('text-anchor','middle')
+        text.set('dominant-baseline','central')
+        text.set('font-family','Gill Sans MT')
+        #stext.set('font-family','Cambria')
+        matrixGroup.append(text)
+
+
+    def columnIndexToXCoordinate(self, columnIndex, objSize=None):
+        objSize = self.cellSize if objSize is None else objSize
+        return(  (self.cellSize + (self.cellSize-objSize)*0.5)  + (columnIndex-1) * self.cellSize+ self.cellSpacing*columnIndex + int(columnIndex/self.groupSize)*self.groupSpacerSize)
+        #return(columnIndex * (self.cellSize+self.cellSpacing) + int(columnIndex/4)*self.groupSpacerSize)
+
+
+    def rowIndexToYCoordinate(self, rowIndex, objSize=None):
+        objSize = self.cellSize if objSize is None else objSize
+        return( (self.cellSize + (self.cellSize-objSize)*0.5) + (rowIndex-1) * (self.cellSize) + self.cellSpacing*rowIndex + int(rowIndex/self.groupSize)*self.groupSpacerSize)
+
+
+    def nominalColor(self, value):
+        if self.nominalColoringMapping == None:
+            #Find all unique values in the matrix:
+
+            uniqueValues = list(set( v for _,v in np.ndenumerate( self.matrix.astype(str ) ) ))
+            try:
+                uniqueValues = sorted(uniqueValues)
+            except:
+                pass
+
+            if len(uniqueValues)!=0:
+                self.nominalColoringMapping = {val:self.colormap(float(index)/float(len(uniqueValues))) for index,val in enumerate(uniqueValues)}
+            else:
+                return(self.NanColor)
+            print(self.nominalColoringMapping)
+
+
+        return( self.nominalColoringMapping.get(value, self.NanColor) )
+
+    def isnan(self,value):
+
+        theValueIsNan = False
+        if value is None:
+            theValueIsNan = True
+        else:
+            try:
+                theValueIsNan = np.isnan(value)
+            except:
+                theValueIsNan = False
+        return(theValueIsNan)
+
+    def getPlot(self):
+
+        plot = BDBPlot()
+
+        matrixGroup = plot.getGroup('matrix')
+        matrixGroup.set('transform', 'translate(%s, %s)' % (self.leftMargin,self.topMargin ))
+        columnCount, rowCount = self.matrix.shape
+
+        plotWidth = self.leftMargin + (columnCount*(self.cellSize+self.cellSpacing )) + self.groupSpacerSize*(columnCount/4)
+        plotHeight = self.topMargin + self.rowIndexToYCoordinate(rowCount) +self.cellSize+self.cellSpacing  + self.footerHeight
+        plot.setWidth(plotWidth)
+        plot.setHeight(plotHeight)
+        ySlack = self.cellSize/2
+
+        cellShadow = plot.shadow(0.5,0.5,1, 'rgb(0,0,0)', 0.98)
+
+        plot.addDef(cellShadow)
+
+
+        foreGroundTilesGroup = plot.getGroup('foreGroundTiles')
+        foreGroundTilesGroup.set('style', 'filter:url(#%s);'%cellShadow.get('id') )
+        matrixGroup.append(foreGroundTilesGroup)
+        for (column,row), value in np.ndenumerate(self.matrix):
+
+            zIndex = 0
+            cellSize = self.cellSize*0.75 if self.getCellRotation(row, column)!=None and self.getCellRotation(row, column)!=0 and self.getCellRotation(row, column)!='nan' else self.cellSize
+            x = self.columnIndexToXCoordinate(column,cellSize)
+            y = self.rowIndexToYCoordinate(row,cellSize)
+            rect = plot.getRectangle(x,y,cellSize, cellSize)
+
+            try:
+                cVal = self.colorMatrix[column,row]
+            except:
+                cVal = value
+
+
+            theValueIsNan = self.isnan(cVal)
+            r,g,b,a = self.getColor(value)
+
+            r = int(r*255.0)
+            g = int(g*255.0)
+            b = int(b*255.0)
+            #print('rgb(%s,%s,%s)' %  (r,g,b))
+
+
+            rect.set( 'fill','rgb(%s,%s,%s)' %  (r,g,b))
+            plot.modifyStyle(rect, {'fill': 'rgba(%s,%s,%s,1)' % (r,g,b), 'stroke-width':'0', 'stroke':'rgba(%s,%s,%s,1)' % (0,0,0)})
+
+            if self.getCellId(row,column) is not None:
+                rect.set('cell_id',self.getCellId(row,column))
+            if self.getCellMetaData(row,column) is not None:
+                rect.set('meta',self.getCellMetaData(row,column))
+
+            if self.getCellRotation(row, column)!=None and self.getCellRotation(row, column)!=0 and self.getCellRotation(row, column)!='nan':
+                rect.set('transform','rotate(%s, %s, %s)'%(self.getCellRotation(row, column),x+cellSize*0.5, y+cellSize*0.5))
+                zIndex = len(matrixGroup)-1
+
+            if theValueIsNan:
+                matrixGroup.insert( 0, rect)
+            else:
+                foreGroundTilesGroup.insert( zIndex, rect)
+
+            brightness = 255 - ((r+g+b)/3.0)
+            if ((r+g+b)/3.0)<100:
+                c = BDBcolor( brightness, brightness, brightness, 1)
+            else:
+                c = BDBcolor( 0, 0, 0, 1)
+
+            if theValueIsNan==False:
+                if self.cellAnnot is None:
+                    text = plot.getText(str(self.getCellString(row,column)), x+0.5*cellSize, y+0.5*cellSize, fill=c.getRGBStr())
+                    text.set('text-anchor','middle')
+                    text.set('dominant-baseline','central')
+                    #text.set('font-family','Gill Sans MT')
+                    text.set('font-family',self.cellFont)
+                    text.set('font-size', str(self.cellFontSize))
+                    matrixGroup.append(text)
+                else:
+                    text = plot.getText(str(self.getCellString(row,column)), x+0.5*cellSize, y+0.3*cellSize, fill=c.getRGBStr())
+                    text.set('text-anchor','middle')
+                    text.set('dominant-baseline','central')
+                    #text.set('font-family','Gill Sans MT')
+                    text.set('font-family',self.cellFont)
+                    text.set('font-size', str(self.cellFontSize))
+                    #text.set('font-weight','bold')
+                    matrixGroup.append(text)
+
+                    text = plot.getText(str(self.getCellAnnotString(row,column)), x+0.5*cellSize, y+0.75*cellSize, fill=c.getRGBStr())
+                    text.set('text-anchor','middle')
+                    text.set('dominant-baseline','central')
+                    #text.set('font-family','Gill Sans MT')
+                    text.set('font-family',self.cellFont)
+                    text.set('font-size', str(self.cellFontSize*0.90))
+                    matrixGroup.append(text)
+
+
+            if column==0:
+                #x -= self.labelWidth
+
+                self.leftMargin = max( len( self.getRowName(row) )*8, self.leftMargin )
+
+                text = plot.getText(self.getRowName(row), self.columnIndexToXCoordinate(column)-self.labelWidth,  self.rowIndexToYCoordinate(row)+ySlack, fill=self.getRowColor(row))
+                text.set('text-anchor','end')
+                text.set('dominant-baseline','central')
+                #text.set('font-family','Gill Sans MT')
+                text.set('font-family',self.labelFont)
+                text.set('font-size', str(self.labelFontSize))
+
+                matrixGroup.append(text)
+            if row==(rowCount-1) or row==0:
+                offset = self.labelWidth+self.cellSize if row==(rowCount-1) else -self.cellSize*0.5
+                text = plot.getText(self.getColName(column), self.columnIndexToXCoordinate(column)+self.cellSize*0.5,  self.rowIndexToYCoordinate(row)+offset, fill=self.getColumnColor(column))
+                text.set('text-anchor','middle')
+                text.set('dominant-baseline','central')
+                #text.set('font-family','Gill Sans MT')
+                text.set('font-family',self.labelFont)
+                text.set('font-size', str(self.labelFontSize))
+                if self.rotateColumnLabels!=0:
+                    plot.setTextRotation(text, self.rotateColumnLabels)
+                    if self.rotateColumnLabels>90:
+                        text.set('text-anchor','end' if row==(rowCount-1) else 'start')
+                    else:
+                        text.set('text-anchor','start' if row==(rowCount-1) else 'end')
+                matrixGroup.append(text)
+
+        plot.svgTree.append(matrixGroup)
+        self.addTitle(plot)
+        self.addSubtitle(plot)
+
+        matrixGroup.set('transform', 'translate(%s, %s)' % (self.leftMargin,self.topMargin ))
+        plotWidth = self.leftMargin + (columnCount*(self.cellSize+self.cellSpacing )) + self.groupSpacerSize*(columnCount/4)
+        plotHeight = self.topMargin + self.rowIndexToYCoordinate(rowCount) +self.cellSize+self.cellSpacing  + self.footerHeight
+        plot.setWidth(plotWidth)
+        plot.setHeight(plotHeight)
+        return(plot)
+
+
+
+def testHeatmap():
+    xmin = 0.0
+    xmax = 10.0
+    dx = 1.0
+    ymin=0.0
+    ymax=5.0
+    dy = 1.0
+    x,y = np.meshgrid(np.arange(xmin,xmax,dx),np.arange(ymin,ymax,dy))
+    npMat = (x*y)
+    m = npMat.max()
+    print(m)
+    npMat /= m
+    print(npMat)
+    xlabels = ["X"+str(x) for x in range(npMat.shape[1])]
+    ylabels = ["Y"+str(y) for y in range(npMat.shape[0])]
+    h = Heatmap(npMat, npMat, rowNames=xlabels, columnNames=ylabels)
+    p = h.getPlot()
+    p.write('test.svg')
+
+
+def histogram(values = [1,7,3,2,1,0,0,0,1], rebin=False, binCount=9, reScale=False, logScale=False, logScaleData=False):
+
+    if rebin:
+        newBars = {}
+        bars = [0]*(binCount+1)
+        frequencies = dict()
+        for v in values:
+            if not v in frequencies:
+                frequencies[v]=1
+            else:
+                frequencies[v]+=1
+
+        #Take log of frequencies
+        if logScale:
+            for q in frequencies:
+
+                if frequencies[q]>0:
+                    frequencies[q] = math.log10(frequencies[q])
+                else:
+                    frequencies[q] = -1
+
+
+        minValue = min([float(x) for x in list(frequencies.keys())])
+        maxValue = max([float(x) for x in list(frequencies.keys())])
+        binSize = (maxValue - minValue)/binCount
+
+        sampleTotal = sum(frequencies.values())
+
+        for binIndex in range(0,binCount+1):
+            binStart = minValue+ binIndex*binSize
+            binEnd = binStart+binSize
+            binTotal = 0
+            for d in frequencies.irange(binStart, binEnd, (True,False)):
+                binTotal+=frequencies[d]
+                if reScale:
+                    newBars[binStart+0.5*binSize] =  float(binTotal)/float(sampleTotal)
+                    bars[binIndex] = float(binTotal)/float(sampleTotal)
+                else:
+                    newBars[binStart+0.5*binSize] =  binTotal
+                    bars[binIndex] = binTotal
+
+
+    else:
+        bars=values
+
+
+    bdbplot = BDBPlot()
+
+
+    bdbplot.plotStartX = 100
+    bdbplot.plotStartY = 100
+    bdbplot.xMax = max(1,int(math.ceil(len(bars)))) # prevent 0 (breaks everything, 0 divisions and such)
+    bdbplot.yMax = max(1,int(math.ceil(max(bars))))
+
+
+    #plotWall = bdbplot.getRectangle(bdbplot.plotStartX,bdbplot.plotStartY,bdbplot.plotWidth,bdbplot.plotHeight)
+    #bdbplot.svgTree.append( plotWall )
+    #bdbplot.modifyStyle(plotWall, {'filter':'url(#%s)'%shadow.get('id'),'fill':'rgba(255,255,255,1)'})
+    #bdbplot.modifyStyle(plotWall, {'fill':'rgba(255,255,255,1)'})
+
+    axis = bdbplot.getAxis()
+    bdbplot.svgTree.append( axis )
+
+    #Draw fine grid
+    for y in range(1,bdbplot.yMax+1):
+        p = bdbplot.getPath(bdbplot.getPathDefinition([bdbplot.getPlottingCoord(bdbplot.xMin, y),bdbplot.getPlottingCoord(bdbplot.xMax, y)]))
+
+        bdbplot.modifyStyle(p, {'stroke-width':'0.5', 'stroke-linecap':'round', 'stroke-dasharray':'2 2','stroke-dashoffset':'0'} )
+        bdbplot.svgTree.append( p )
+
+
+    # Add axis labels:
+    if logScaleData:
+        for x in range(0,bdbplot.xMax+1):
+            c = bdbplot.getXLabelCoord(x)
+
+            text = bdbplot.getText(str(10), c[0], c[1]+15,BDBcolor(80,80,80,1))
+            text.set('text-anchor','middle')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+            bdbplot.addSuper(text,str(x))
+            bdbplot.svgTree.append( text )
+
+    else:
+        for x in range(0,bdbplot.xMax+1):
+            c = bdbplot.getXLabelCoord(x)
+            text = bdbplot.getText(str(x), c[0], c[1]+15,BDBcolor(80,80,80,1))
+            text.set('text-anchor','middle')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+            bdbplot.svgTree.append( text )
+    if logScale:
+        for y in range(1,bdbplot.yMax+1):
+            c = bdbplot.getYLabelCoord(y)
+
+            value = math.pow(10,y)
+
+            text = bdbplot.getText(str(10), c[0]-10, c[1],BDBcolor(80,80,80,1))
+            text.set('text-anchor','end')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+            bdbplot.addSuper(text,str(y))
+            bdbplot.svgTree.append( text )
+
+    else:
+        for y in range(0,bdbplot.yMax+1):
+            c = bdbplot.getYLabelCoord(y)
+            text = bdbplot.getText(str(y), c[0]-10, c[1],BDBcolor(80,80,80,1))
+            text.set('text-anchor','end')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Cambria Math')
+            bdbplot.svgTree.append( text )
+
+
+
+    c = bdbplot.getYLabelCoord( (bdbplot.yMax/2))
+    text = bdbplot.getText('Frequency', c[0]-60, c[1]-30,BDBcolor(0,0,0,1))
+    text.set('text-anchor','middle')
+    text.set('dominant-baseline','middle')
+    text.set('font-family','Gill Sans MT')
+    text.set('font-size', '25')
+    #bdbplot.modifyStyle(text, {'font-size': '20'})
+    bdbplot.setTextRotation(text,270)
+    bdbplot.svgTree.append( text )
+
+    c = bdbplot.getXLabelCoord(bdbplot.xMax/2)
+    text = bdbplot.getText('Read abundance', c[0], c[1]+50,BDBcolor(0,0,0,1))
+    text.set('text-anchor','middle')
+    text.set('dominant-baseline','middle')
+    text.set('font-family','Gill Sans MT')
+    text.set('font-size', '25')
+    bdbplot.svgTree.append( text )
+
+
+    barShadow = bdbplot.shadow(1,1)
+    bdbplot.addDef(barShadow)
+
+    for barIndex,barValue in enumerate(bars):
+        c = bdbplot.getPlottingCoord(barIndex,barValue)
+        origin = bdbplot.getPlottingCoord(barIndex,0)
+        bar = bdbplot.getRectangle(  c[0], c[1], float(bdbplot.plotWidth)/(len(bars)+1),  (float(barValue)/bdbplot.yMax) * bdbplot.plotHeight )
+
+        bdbplot.modifyStyle(bar, {'filter':'url(#%s)'%barShadow.get('id'),'fill':'#FFFFFF'})
+        bdbplot.svgTree.append( bar )
+
+    for barIndex,barValue in enumerate(bars):
+        c = bdbplot.getPlottingCoord(barIndex,barValue)
+        origin = bdbplot.getPlottingCoord(barIndex,0)
+
+        barWidth = float(bdbplot.plotWidth)/(len(bars)+1)
+
+        bar = bdbplot.getRectangle(  c[0], c[1], barWidth,  (float(barValue)/bdbplot.yMax) * bdbplot.plotHeight )
+        bdbplot.modifyStyle(bar, {'fill':bdbplot.getGroupColors(1)[0], 'stroke':'#FFFFFF','stroke-width':'1.5'})
+        bdbplot.svgTree.append( bar )
+
+        if barValue!=-1:
+            text = bdbplot.getText(str( bdbplot.humanReadable( int(math.pow(10,barValue)) ) ), c[0]+0.5*barWidth, c[1]-10,BDBcolor(0,0,0,1))
+            text.set('text-anchor','middle')
+            text.set('dominant-baseline','middle')
+            text.set('font-family','Gill Sans MT')
+            #text.set('font-family','Cambria Math')
+            text.set('font-size', '14')
+            text.set('font-weight', 'bold')
+            text.set('fill', 'rgba(50,50,50,1)')
+            bdbplot.svgTree.append( text )
+
+    return(bdbplot)
+
+def densityXY(scatterData, plotPath, xlabel='x', ylabel='y', logX=False, forceShow=False, logY=False):
+
+
+    from scipy.stats import gaussian_kde
+
+    if len(scatterData['x'])==0:
+        #self.warn('No datapoints left for comparison')
+        return(False)
+
+    #@todo: cast this earlier.
+    scatterData['x'] = np.array(scatterData['x'])
+    scatterData['y'] = np.array(scatterData['y'])
+
+    xy = np.vstack([scatterData['x'],scatterData['y']])
+    z = gaussian_kde(xy)(xy)
+
+    # Sort the points by density, so that the densest points are plotted last
+    idx = z.argsort()
+    scatterData['x'] = scatterData['x'][idx]
+    scatterData['y'] = scatterData['y'][idx]
+    z = z[idx]
+
+    plt.close('all')
+    fig, ax = plt.subplots()
+    ax.scatter(scatterData['x'], scatterData['y'], c=z, s=50, edgecolor='')
+    if logX:
+        ax.set_xscale('log')
+    if logY:
+        ax.set_yscale('log')
+
+    plt.ylabel(ylabel)
+    plt.xlabel(xlabel)
+    if plotPath and forceShow:
+        plt.show()
+    if plotPath is None:
+        plt.show()
+    else:
+        plt.savefig(plotPath, bbox_inches='tight')
+        plt.close('all')
+
+
+##
+# SIMPLE X Y PLOT
+##
+def simpleXY():
+
+    bdbplot = BDBPlot()
+
+    shadow = bdbplot.shadow()
+    bdbplot.addDef(shadow)
+    bdbplot.plotStartX = 100
+    bdbplot.plotStartY = 100
+    bdbplot.xMax = 10
+    bdbplot.yMax = 10
+
+    #plotWall = bdbplot.getRectangle(bdbplot.plotStartX,bdbplot.plotStartY,bdbplot.plotWidth,bdbplot.plotHeight)
+    #bdbplot.svgTree.append( plotWall )
+    #bdbplot.modifyStyle(plotWall, {'filter':'url(#%s)'%shadow.get('id'),'fill':'rgba(255,255,255,1)'})
+    #bdbplot.modifyStyle(plotWall, {'fill':'rgba(255,255,255,1)'})
+
+    axis = bdbplot.getAxis()
+    bdbplot.svgTree.append( axis )
+
+    # Add axis labels:
+    for x in range(0,11):
+        c = bdbplot.getXLabelCoord(x)
+        text = bdbplot.getText(str(x), c[0], c[1]+15,BDBcolor(80,80,80,1))
+        text.set('text-anchor','middle')
+        text.set('dominant-baseline','middle')
+        text.set('font-family','Cambria Math')
+        bdbplot.addSuper(text,'y')
+        bdbplot.svgTree.append( text )
+
+    for y in range(1,11):
+        c = bdbplot.getYLabelCoord(y)
+        text = bdbplot.getText(str(y), c[0]-15, c[1],BDBcolor(80,80,80,1))
+        text.set('text-anchor','middle')
+        text.set('dominant-baseline','middle')
+        text.set('font-family','Cambria Math')
+        bdbplot.svgTree.append( text )
+
+
+
+    c = bdbplot.getYLabelCoord(5)
+    text = bdbplot.getText('Y Axis', c[0]-40, c[1],BDBcolor(0,0,0,1))
+    text.set('text-anchor','middle')
+    text.set('dominant-baseline','middle')
+    text.set('font-family','Gill Sans MT')
+    bdbplot.setTextRotation(text,270)
+    bdbplot.svgTree.append( text )
+
+    c = bdbplot.getXLabelCoord(5)
+    text = bdbplot.getText('X Axis', c[0], c[1]+40,BDBcolor(0,0,0,1))
+    text.set('text-anchor','middle')
+    text.set('dominant-baseline','middle')
+    text.set('font-family','Gill Sans MT')
+    bdbplot.svgTree.append( text )
+
+
+
+
+    for x in range(0,10):
+
+        c = bdbplot.getPlottingCoord(x,x)
+        circle = bdbplot.getCircle(c[0],c[1],2)
+        #bdbplot.modifyStyle(circle, {'filter':'url(#%s)'%shadow.get('id')})
+        bdbplot.svgTree.append( circle )
+
+    #
+    #for i in range(0,10):
+    #
+    #   a = ((math.pi*2)/10.0) * i
+    #   text = bdbplot.getText('%s' % i,250 + 50*math.cos(a),250 + 50*math.sin(a),BDBcolor(80,80,80,1))
+    #   text.set('text-anchor','middle')
+    #   text.set('dominant-baseline','middle')
+    #   text.set('font-family','Cambria Math')
+    #   bdbplot.svgTree.append( text )
+
+
+    bdbplot.dump()
+    bdbplot.write('test.svg')
+
+
+#vals = [0,11]
+#
+#import random
+#for i in range(0,16):
+#   vals += [i]* int(math.ceil(math.exp(i/2+1)))
+#
+#
+##print('c(%s)' % ','.join(str(i) for i in vals))
+#
+#plot = histogram(vals, True,15, False, True)
+#
+#
+#plot.write('test.svg')
+
+#
+#d = Counter({1:100,2:50,3:10,5:5,6:10,  1000:1, 500:2, 10000:3})
+#plot = readCountHistogram(d)
+#text = plot.getText('Embryo 1, component 1',10,40)
+#text.set('font-family','Gill Sans MT')
+#text.set('font-size', '42')
+#plot.svgTree.append( text )
+#
+#totalCount = sum( [v*d[v]for v in d] )
+#text = plot.getText( '%s reads total' % plot.humanReadable(totalCount),10,75, BDBcolor(77,77,77,1))
+#text.set('font-family','Gill Sans MT')
+#text.set('font-size', '23')
+#plot.svgTree.append( text )
+#
+#
+#plot.write('test.svg')
+
+import networkx as nx
+import scipy.interpolate
+from Bio import pairwise2
+
+
+class GraphRenderer():
+
+    def interpolate(self, interpolateValue,  colorScaleKeys, nodeColorMapping):
+
+        #Seek positions around value to interpolate
+        first = colorScaleKeys[0]
+        index = 0
+        last = first
+        for value in colorScaleKeys:
+
+            if value>=interpolateValue:
+                last = value
+                break
+            else:
+                first = value
+            index+=1
+        if value==interpolateValue:
+            return(nodeColorMapping[value])
+
+        #Do interpolation
+        colorA = nodeColorMapping[first]
+        colorB = nodeColorMapping[last]
+        dx = last-first
+
+        return( self._ipol(colorA[0], colorB[0], first, last, interpolateValue), self._ipol(colorA[1], colorB[1], first, last, interpolateValue), self._ipol(colorA[2], colorB[2], first, last, interpolateValue))
+
+
+    def _ipol(self,a, b, first, last,  interpolateValue):
+        #Due to floating point rounding errors the interpolate value can be very close to last,
+        # it is ok to return last in those cases
+        if last>first and interpolateValue>=last:
+            return(b)
+        if last<first and interpolateValue>=first:
+            return(a)
+
+        y_interp = scipy.interpolate.interp1d([first, last], [a,b])
+        return( y_interp(interpolateValue) )
+
+
+    def sortByIndexAndBase(self, value):
+
+        parts = value.split('_')
+        pos = ['A','T','C','G','N'].index(parts[1])
+        if pos==None:
+            pos=0
+        else:
+            pos+=1
+        return( int(parts[0]) + pos*0.1 )
+
+
+    def __init__(self, nxGraph, coloringMode = 'nodeRGB', coloringAttribute='confidence', performDistanceMeasure=True, performFrequencyMeasure=True, alias='none'):
+        self.g = nxGraph
+        self.undirectedG = self.g.to_undirected()
+        self.plot = BDBPlot()
+
+        self.nodeShadow = self.plot.shadow(0.5,0.5,1, 'rgb(0,0,0)', 0.98)
+        self.plot.addDef(self.nodeShadow)
+
+        minX = 0
+        maxX = 0
+        minY = 0
+        maxY = 0
+        for nodeName in self.g:
+            node = self.g.node[nodeName]
+            if 'x' in node and 'y' in node and 'size' in node :
+                if (node['x']-node['size'])<minX:
+                    minX = node['x']-node['size']
+                if (node['x']+node['size'])>maxX:
+                    maxX = node['x']+node['size']
+
+                if (node['y']-node['size'])<minY:
+                    minY = node['y']-node['size']
+                if (node['y']+node['size'])>maxY:
+                    maxY = node['y']+node['size']
+
+
+
+        #Estimate color scale: #############
+        colorScaleY = 0
+        colorScaleHeight = 0
+        colorScaleGraphSpacing = 0
+        createdColorScale = False
+        if coloringAttribute is not None and coloringMode == 'nodeRGB':
+            createdColorScale= True
+            colorScaleParts = 10
+            colorScaleWidth = 400
+            colorScaleHeight = 35
+            colorScaleSpacing = 5
+            colorScaleShadow = self.plot.shadow(1,1,1)
+            colorScaleDeltaX = float(colorScaleWidth-colorScaleSpacing)/colorScaleParts
+            self.plot.addDef(colorScaleShadow)
+            labelHeight = 15
+            colorScaleX = 5
+            colorScaleGraphSpacing = 10
+            colorScaleY = (maxY-minY)+colorScaleGraphSpacing
+
+            nodeColorMapping = {}
+            abundanceMapping = Counter({})
+            lowestValue = 100000
+            highestValue = -lowestValue
+            for nodeName in self.g:
+                node = self.g.node[nodeName]
+                if coloringAttribute in node and 'r' in node:
+                    if node[coloringAttribute]==1:
+                        value = 0
+                    else:
+                        value = -math.log( 1.0 - node[coloringAttribute],10 )
+
+                    abundanceMapping[value]+= node['abundance']
+                    lowestValue = min(lowestValue, value)
+                    highestValue = max(highestValue, value)
+                    nodeColorMapping[value] = (node['r'], node['g'], node['b'])
+
+            #Create color scale; ##########
+            lowestValue = 2.5
+            highestValue = 3.4
+
+            #lowestValue = 3.0
+            #highestValue = 3.6
+            print((lowestValue, highestValue))
+            colorScale = self.plot.getGroup('colorScale')
+
+            self.plot.svgTree.append(colorScale)
+            r = self.plot.getRectangle(colorScaleX,colorScaleY,colorScaleWidth,colorScaleHeight+labelHeight+colorScaleSpacing )
+            colorScale.append(r)
+
+            r.set('style', 'fill:#FFFFFF' )
+            #r.set('style', 'filter:url(#%s);fill:#FFFFFF'%colorScaleShadow.get('id') )
+            colorScaleKeys = sorted(OrderedDict(sorted(nodeColorMapping.items())) )
+            #print(nodeColorMapping)
+            deltaValue = (highestValue-lowestValue) / (colorScaleParts-1)
+            currentValue = lowestValue
+            x = colorScaleSpacing+colorScaleX
+            y = colorScaleY + colorScaleSpacing
+
+            idx = 0
+            while idx<colorScaleParts:
+                print(('Interpolating for %s' % currentValue ))
+                c = self.interpolate(currentValue, colorScaleKeys, nodeColorMapping)
+                r = self.plot.getRectangle(x,y,colorScaleDeltaX-colorScaleSpacing,colorScaleHeight-2*colorScaleSpacing )
+                r.set('fill',  'rgb(%s, %s, %s)' % c)
+                r.set('style', 'filter:url(#%s);'%  self.nodeShadow.get('id') )
+                r.set('stroke',  'None')
+                colorScale.append( r )
+
+                text = self.plot.getText('%.1f' %  (10.0*currentValue) ,x+0.5*(colorScaleDeltaX-colorScaleSpacing),y+colorScaleHeight+colorScaleSpacing, BDBcolor(0,0,0,1))
+                #self.plot.setTextRotation(text, 90)
+                text.set('text-anchor','middle')
+                text.set('dominant-baseline','central')
+                text.set('font-family','Gill Sans MT')
+                #text.set('font-family','Cambria Math')
+                text.set('font-size', '10')
+                #text.set('font-weight', 'bold')
+                text.set('fill', 'rgba(0,0,0,1)')
+                self.plot.svgTree.append( text )
+
+
+                currentValue+=deltaValue
+                x += colorScaleDeltaX
+                idx+=1
+
+
+        #######################
+
+
+        #Colorize nodes by distance to center
+
+
+        distancesFound = Counter({})
+        distancesReads = Counter({})
+        for node in self.g:
+            if 'idStr' in self.g.node[node] and  'Wt' == self.g.node[node]['idStr']:
+                centerNode = node
+                break
+
+        ldistThreshold = 4
+        maxHamming = 8
+        classColors = {'H0':'#0000DD','H1': '#66A43E', 'H2': '#0D40DB', 'H3': '#3970DD', 'H4': '#769AE0', 'H5': '#A6B9DD', 'H6': '#D7DAE0', 'H7': '#FFFFFF', 'N1':'#FFCC00','N2':'#FF6600', 'N3':'#C83737','N4':'#800000'}
+        if performDistanceMeasure:
+
+            print('Estimating all distances...')
+            weirdSequences = []
+            for nodeIndex,node in enumerate(self.undirectedG):
+
+                if nodeIndex%100==0:
+                    completion = 100.0*(float(nodeIndex)/len(self.undirectedG))
+                print('\rcompletion %s    ' % completion, end=' ')
+
+                abundance = self.g.node[node]['abundance']
+                hammingDistance = bdbbio.getHammingDistance(node, centerNode)
+                unformattedAlignments = pairwise2.align.localxx(node,centerNode) #bdbbio.getLevenshteinDistance(node,centerNode)
+                ldist = len(node)
+
+                self.g.node[node]['exactHammingDistance'] =  hammingDistance
+                self.g.node[node]['exactNWDistance'] =  ldist
+                if len(unformattedAlignments)>0:
+                    ldist = int(round(len(node)-float(unformattedAlignments[0][2])))
+                    #pairwise2.format_alignment(*unformattedAlignments[0])
+
+                if hammingDistance<maxHamming and ldist<=hammingDistance:
+                    distancesReads['H%s'%hammingDistance]+=abundance
+
+                    if hammingDistance==0:
+                        self.g.node[node]['color'] = '#66A43E'
+                        distancesFound['H0']+=1
+                    elif hammingDistance==1:
+                        self.g.node[node]['color'] = '#0D40DB'
+                        distancesFound['H1']+=1
+                    elif hammingDistance==2:
+                        self.g.node[node]['color'] = '#2362E0'
+                        distancesFound['H2']+=1
+                    elif hammingDistance==3:
+                        self.g.node[node]['color'] = '#769AE0'
+                        distancesFound['H3']+=1
+                    elif hammingDistance==4:
+                        self.g.node[node]['color'] = '#A6B9DD'
+                        distancesFound['H4']+=1
+                    elif hammingDistance==5:
+                        self.g.node[node]['color'] = '#D7DAE0'
+                        distancesFound['H5']+=1
+                    elif hammingDistance==6:
+                        self.g.node[node]['color'] = '#FFFFFF'
+                        distancesFound['H6']+=1
+                else:
+
+
+                    if ldist<=ldistThreshold:
+
+                        distancesFound['N%s'%ldist]+=1
+                        distancesReads['N%s'%ldist]+=abundance
+                        if 'N%s'%ldist not in classColors:
+
+                            brightness = 255-int(round((float(ldist)/ldistThreshold)*100))
+                            self.g.node[node]['color'] = 'rgb(%s,%s,%s)' % (brightness,0,0)
+                            classColors['N%s'%ldist] = self.g.node[node]['color']
+                        self.g.node[node]['color'] = classColors['N%s'%ldist]
+                    else:
+                        self.g.node[node]['color'] = '#404040'
+                        #distancesFound['l%s'%ldist]+=1
+                        #distancesReads['l%s'%ldist]+=abundance
+                        distancesFound['N>%s'%ldistThreshold]+=1
+                        distancesReads['N>%s'%ldistThreshold]+=abundance
+                        weirdSequences.append(SeqIO.SeqRecord(Seq(node), 'NW%s-a%s-%s' % (str(ldist),abundance,str(nodeIndex))))
+
+            #classHistogram( classCountMapping, logTransform = False, classColors=None, placeLeft=None ):
+
+            classHistogram(distancesFound, True, classColors, None,['N>%s'%ldistThreshold], False).write('%s_distancesByNodes.svg' % alias)
+            classHistogram(distancesReads, True, classColors, None,['N>%s'%ldistThreshold], False).write('%s_distancesByCountB.svg' % alias)
+            nx.write_graphml( self.g, './%s-distanceAnnotated.graphml' % alias)
+            fastaPath = './%s-weirdSequences.fa' % alias
+            SeqIO.write(weirdSequences, fastaPath, "fasta")
+
+
+
+
+        if performFrequencyMeasure:
+
+            sequenceColors = {'rest':'#404040'}
+            sequenceFrequencies = Counter({})
+            frequencyTable = []
+            frequencyCounter = Counter({})
+            for nodeIndex,node in enumerate(self.undirectedG):
+                abundance = self.g.node[node]['abundance']
+                nodeName = nodeIndex
+                if 'idStr' in self.g.node[node]:
+                    nodeName = self.g.node[node]['idStr']
+                #if 'r' in self.g.node[node]:
+                #   sequenceColors[str(nodeName)] = 'rgb(%s, %s, %s)' % (self.g.node[node]['r'],self.g.node[node]['g'],self.g.node[node]['b'])
+                if 'color' in self.g.node[node]:
+                    sequenceColors[str(nodeName)] = self.g.node[node]['color']
+
+                if abundance>3:
+                    sequenceFrequencies[str(nodeName)] += abundance
+
+                else:
+                    #sequenceFrequencies['rest'] += abundance
+                    pass
+                for index,base in enumerate(str(node)):
+
+                    if index>=len(frequencyTable):
+                        frequencyTable.append(Counter({}))
+                        for b in ['A','T','C','G']:
+                            frequencyCounter['%s_%s' % (index, b)] = 0
+                    frequencyCounter['%s_%s' % (index, base)]+=abundance
+                    frequencyTable[index][base]+= abundance
+
+
+            freqKeys = sorted(list(frequencyCounter.keys()), key=lambda x: self.sortByIndexAndBase(x))
+            colors = {}
+            for key in freqKeys:
+                parts = key.split('_')
+                i = ['N','A','T','C','G'].index(parts[1])
+                if i==None:
+                    i = 0
+                colors[key] = ['#404040', '#336bbd','#ff6600','#aa0000','#5aa02c'][i]
+
+
+            classHistogram(sequenceFrequencies, True, sequenceColors, None,[], False, classWidth=40, rotateClassLabels=90).write('%s_sequenceFrequenciesLog.svg' % alias)
+            classHistogram(sequenceFrequencies, False, sequenceColors, None,[], False, classWidth=40, rotateClassLabels=90).write('%s_sequenceFrequenciesLin.svg' % alias)
+            classHistogram(frequencyCounter, True, colors, None,freqKeys, False, classWidth=30, rotateClassLabels=90).write('%s_baseFrequencies.svg' % alias)
+
+
+
+        print(('Offsetting %s %s' % (minX, minY)))
+        self.plot.setWidth( maxX-minX )
+        if createdColorScale:
+            self.plot.setHeight( ( (colorScaleY+colorScaleHeight+colorScaleGraphSpacing)-minY ))
+        else:
+            self.plot.setHeight( maxY - minY + 10)
+        h = maxY-minY
+        ## plotting
+        edgeGroup = self.plot.getGroup('edges')
+        nodeGroup = self.plot.getGroup('nodes')
+        labelGroup = self.plot.getGroup('labels')
+
+        print('Adding edges')
+
+        for fromNode,toNode,data in self.g.edges(data=True):
+
+            if 'hdist' in data and data['hdist']==1:
+                fn = self.g.node[fromNode]
+                tn = self.g.node[toNode]
+                p = self.plot.getPath( self.plot.getPathDefinition([ (int(round(fn['x']-minX)), h-int(round(fn['y']-minY))), (int(round(tn['x']-minX)), h-int(round(tn['y']-minY))) ]) )
+
+                etree.strip_attributes(p,'style')
+                if coloringMode == 'nodeRGB':
+                    if 'r' in tn:
+                        p.set('stroke', 'rgb(%s, %s, %s)' % (tn['r'],tn['g'],tn['b']))
+                else:
+                    if 'color' in tn:
+                        p.set('stroke', '%s' % (tn['color']))
+                p.set('stroke-width', '0.75')
+                p.set('stroke-opacity', '0.8')
+                #self.plot.modifyStyle(p, { 'stroke-width':'0.5', 'stroke-opacity':"0.5"  } ) #,'stroke-width':'0.5' 'stroke':'rgba(80,80,80,0.8)'})
+                #self.plot.modifyStyle(p, { 'stroke':'rgba(80,80,80,0.8)' } ) #,'stroke-width':'0.5' 'stroke':'rgba(80,80,80,0.8)'})
+
+                edgeGroup.append(p)
+
+        #self.plot.modifyStyle(edgeGroup, {'stroke-width':'0.5', 'stroke':'rgba(80,80,80,0.8)'})
+        self.plot.svgTree.append(edgeGroup)
+        self.plot.svgTree.append(nodeGroup)
+        self.plot.svgTree.append(labelGroup)
+        print('Adding nodes')
+
+
+
+
+
+        for componentIndex,connectedComponent in enumerate(nx.connected_component_subgraphs(self.undirectedG)):
+            componentGroup = self.plot.getGroup('component_%s' % componentIndex)
+            smallNodes = self.plot.getGroup('component_%s_smallNodes' % componentIndex)
+            bigNodes = self.plot.getGroup('component_%s_bigNodes' % componentIndex)
+            componentGroup.append(smallNodes)
+            componentGroup.append(bigNodes)
+            nodeGroup.append(componentGroup)
+            for nodeName in connectedComponent:
+                node = self.g.node[nodeName]
+                if 'x' in node and 'y' in node and 'size' in node :
+
+
+                    if node['size']>0:
+                        circle = self.plot.getCircle(int(round(node['x']-minX)), h-int(round(node['y']-minY)), int(round(node['size'])))
+                        circle.set('style','fill:none')
+                        if coloringMode == 'nodeRGB' and 'r' in node:
+                            self.plot.modifyStyle(circle, {'fill':'rgb(%s,%s,%s)' % (node['r'], node['g'], node['b']), 'stroke':'rgba(0,0,247,0.8)'})
+                        else:
+                            if 'color' in node:
+                                self.plot.modifyStyle(circle, {'fill':'%s' % (node['color']), 'stroke':'rgba(0,0,247,0.8)'})
+
+                        if node['size']>0.0:
+                            bigNodes.append(circle)
+                            bigNodes.append(circle)
+                            if node['abundance']>500:
+                                if 'idStr' in node:
+                                    label = node['idStr']
+                                else:
+                                    label= ''
+
+                                text = self.plot.getText('%s %s' % (label,node['abundance']) ,int(round(node['x']-minX)), h-int(round(node['y']-minY))+4,BDBcolor(80,80,80,1))
+                                text.set('text-anchor','middle')
+                                text.set('dominant-baseline','central')
+                                text.set('font-family','Gill Sans MT')
+                                #text.set('font-family','Cambria Math')
+                                text.set('font-size', '14')
+                                #text.set('font-weight', 'bold')
+                                text.set('fill', 'rgba(50,50,50,1)')
+                                labelGroup.append( text )
+
+
+                        else:
+                            smallNodes.append(circle)
+                            smallNodes.append(circle)
+
+                else:
+                    print('Skipped a node; could not find coordinates')
+
+            bigNodes.set('style', 'filter:url(#%s);'%self.nodeShadow.get('id') )
+            #self.plot.modifyStyle(bigNodes, {'filter':'url(#%s)'%self.nodeShadow.get('id')})
+
+
+
+
+
+
+
+
+def testGraphRenderer():
+    p = GraphRenderer( nx.read_graphml("C:\\Users\BuysDB\Desktop\Control1-g3.graphml"),'distances',alias='controlDist')
+    p.plot.write('control_spaghettogramRenderDistancesB.svg')
+    p.plot.SVGtoPNG('control_spaghettogramRenderDistancesB.svg', 'control_spaghettogramRenderDistancesB.png',2048)
+
+    p = GraphRenderer( nx.read_graphml("C:\\Users\BuysDB\Desktop\Control1-g3.graphml"),'nodeRGB',alias='controlConf')
+    p.plot.write('control_spaghettogramRenderConfidenceB.svg')
+    p.plot.SVGtoPNG('control_spaghettogramRenderConfidenceB.svg', 'control_spaghettogramRenderConfidenceB.png',2048)
+
+
+def artGraphRenderer():
+    p = GraphRenderer( nx.read_graphml("C:\\Users\BuysDB\Desktop\ArtSimulatedGraphHC26k.graphml"),'nodeRGB','confidence',False, True, 'simulated')
+    p.plot.write('ArtSimulatedGraphHC26k.svg')
+    p.plot.SVGtoPNG('ArtSimulatedGraphHC26k.svg', 'ArtSimulatedGraphHC26k.png',2048)
+
+def embryoGraphRenderer():
+    p = GraphRenderer( nx.read_graphml("C:\\Users\BuysDB\Desktop\embryo7remapped2.graphml"),'nodeRGB','confidence',False, True, 'embryo7')
+    p.plot.write('embryo7remapped3.svg')
+    p.plot.SVGtoPNG('embryo7remapped3.svg', 'embryo7remapped3.png',2048)
+
+def midbrainGraphRenderer():
+    p = GraphRenderer( nx.read_graphml("C:\\Users\BuysDB\Desktop\midbrain.graphml"),'nodeRGB','confidence',False, True, 'brain')
+    p.plot.write('midbrain.svg')
+    p.plot.SVGtoPNG('midbrain.svg', 'midbrain.png',2048)
+
+class SequenceBin():
+
+    def __init__(self, sequence, abundance):
+        self.sequence = sequence
+        self.abundance = abundance
+        self.confidences = []
+        self.diffIndices = []
+        self.x = 0
+        self.y = 0
+
+
+
+
+class HammingBin():
+
+    def __init__(self, index, hammingDistance=0):
+        self.hammingDistance = hammingDistance
+        self.index = index
+        self.sequences = []
+
+    def addSequence(self, sequence, abundance):
+        self.sequences.append(sequence)
+
+
+
+#SVG table renderer
+class SVGTable(BDBPlot):
+    #@param datamatrix list of lists containing values
+    #@param header list containing column names
+    def __init__(self, dataMatrix, header):
+        BDBPlot.__init__(self)
+        self.data = dataMatrix
+        self.header = header
+        self.cellPointer = 0
+
+    #Render a cell in the matrix
+    def cellRenderFunction(self,x,y):
+        self.svgTree.getGroup()
+
+
+
+class SpaghettoPlot():
+
+    def __init__(self, networkxGraph):
+        self.g = networkxGraph
+        self.minRadius = 1
+        self.nucleotideColours = {'A':'#FF2222','T':'#22FF22', 'G':'#2222FF','C':'#FFFF22'}
+
+
+    def getSurfaceBasedNodeRadius(self, abundance, maxRadius):
+
+        #return(math.log(abundance+1)+1)
+        maxO = math.pi*math.pow(float(maxRadius),2)
+        return( max(self.minRadius,math.sqrt( float(abundance*maxO)/math.pi) ) )
+
+
+
+    def layout(self):
+
+
+        maxRadius = 100
+        minDistance = 15 #Distance between the nodes
+
+        components = []
+        self.readLen = 0
+
+        toRemove = []
+        for nodeA,nodeB,d in self.g.edges_iter(data='hdist'):
+            if d!=1:
+                toRemove.append( (nodeA, nodeB) )
+        self.g.remove_edges_from(toRemove)
+
+        for componentIndex,connectedComponent in enumerate(nx.connected_component_subgraphs(self.g)):
+            if len(connectedComponent)>3:
+                #Find center(s)
+                plot = BDBPlot()
+
+                centerNode = None
+                centerAbundance = 0
+                for node in connectedComponent:
+                    a = connectedComponent.node[node]['abundance']
+                    if a > centerAbundance:
+                        centerNode = node
+                        centerAbundance = a
+                    #Todo: compat for more centers
+                    self.readLen = max(self.readLen, len(node))
+                #Estimate radial size of connected component:
+                #longest path from center to member
+                radialSize = 1
+
+                distanceMap = {0:[centerNode]}
+                for targetNode in connectedComponent:
+                    if targetNode!=centerNode:
+
+                        pathLen = nx.shortest_path_length(self.g,source=centerNode,target=targetNode)
+                        radialSize = max( radialSize, pathLen)
+                        if not pathLen in distanceMap:
+                            distanceMap[pathLen] = []
+                        distanceMap[pathLen].append(targetNode)
+
+
+                #print(distanceMap)
+
+                #Read index to angle mapping
+                phis = []
+                for index in range(0,self.readLen):
+                    phis.append( -math.pi*0.5 + float(math.pi*2) * (float(index)/self.readLen) )
+
+
+
+                currentRadius = self.getSurfaceBasedNodeRadius(float(self.g.node[centerNode]['abundance'])/centerAbundance,maxRadius)
+                centerRadius = currentRadius
+                # Construct coordinates
+                coordinates = {}
+                coordinates[centerNode] = {'x':0,'y':0,'r':currentRadius}
+
+                hammingRadials = []
+                for distance in range(1,radialSize+1):
+
+                    print(('Distance %s radius: %s' % (distance, currentRadius)))
+                    if distance in distanceMap:
+                        #Find the radius of this circle
+                        maxNodeRadius = 0
+                        for node in distanceMap[distance]:
+                            maxNodeRadius = max(maxNodeRadius, self.getSurfaceBasedNodeRadius(float(self.g.node[node]['abundance'])/centerAbundance,maxRadius))
+
+                        currentRadius += (maxNodeRadius+ minDistance) #added *.5!
+
+                        #Calculate x and y coordinates:
+
+                        for node in distanceMap[distance]:
+                            r = self.getSurfaceBasedNodeRadius(float(self.g.node[node]['abundance'])/centerAbundance,maxRadius)
+
+                            #Retrieve the hamming index of the node
+                            indices = bdbbio.getHammingIndices(node,centerNode) #[ int(x) for x in self.g.node[node]['indices'].split(',') ]
+
+                            #Take the first index as base
+                            for index in indices:
+                                angle = phis[ index ]
+
+                                coordinates[node] = {'x':math.cos(angle)*currentRadius, 'y':math.sin(angle)*currentRadius, 'r':r, 'a':angle }
+                            #print(coordinates[node])
+                    hammingRadials.append(currentRadius)
+
+                translateX = -currentRadius
+                translateY = -currentRadius
+                #for nodeName in coordinates:
+                #   translateX = min(translateX, coordinates[nodeName]['x']-coordinates[nodeName]['r'] )
+                #   translateY = min(translateY, coordinates[nodeName]['y']-coordinates[nodeName]['r'] )
+                #Draw component
+
+                for i,r in enumerate(hammingRadials):
+
+                    circle = plot.getCircle( -translateX,  -translateY, r)
+                    plot.modifyStyle(circle, {'stroke-width':'0.3', 'stroke-miterlimit':'4', 'stroke-dasharray':'2.08,2.08'})
+                    plot.svgTree.append( circle )
+
+                    text = plot.getText('h'+str(i+1),-translateX, -r-translateY+5, fill=BDBcolor(30,30,30,0))
+                    text.set('text-anchor','middle')
+                    text.set('dominant-baseline','central')
+                    text.set('font-family','Cambria')
+                    text.set('font-size','6')
+                    plot.svgTree.append(text)
+
+
+                    if i==0:
+                        r=centerRadius*2 + 10
+                        for index, angle in enumerate(phis):
+                            text = plot.getText(centerNode[index], math.cos(angle)*r*0.5 - translateX, math.sin(angle)*r*0.5-translateY, fill=BDBcolor(30,30,30,0))
+                            text.set('text-anchor','middle')
+                            text.set('dominant-baseline','central')
+                            text.set('font-family','Gill Sans MT')
+                            text.set('font-size','10')
+                            plot.svgTree.append(text)
+
+                            p = r - 35
+                            text = plot.getText(str(index), math.cos(angle)*p*0.5 - translateX, math.sin(angle)*p*0.5-translateY, fill=BDBcolor(30,30,30,0))
+                            text.set('text-anchor','middle')
+                            text.set('dominant-baseline','central')
+                            text.set('font-family','Cambria')
+                            text.set('font-size','8')
+                            plot.svgTree.append(text)
+
+
+
+                for nodeName in coordinates:
+
+                    circle = plot.getCircle( coordinates[nodeName]['x']-translateX,  coordinates[nodeName]['y']-translateY, coordinates[nodeName]['r'])
+                    plot.modifyStyle(circle, {'stroke-width':'0', 'fill':'#FF6655', 'fill-opacity':'0.30'})
+                    plot.svgTree.append( circle )
+
+                plot.write('./components/component%s.svg' % componentIndex )