Datasets
Models
Downloads: 1
Diff of
/singlecellmultiomics/utils/bdbplot.py
[000000]
..
[2c420a]
Switch to side-by-side view
--- 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 )