Merge branch 'Joint'

Author: Peter Combs

BindUtils.py

@@ -0,0 +1,62 @@
+import pandas as pd
+from bisect import bisect
+from glob import glob
+from collections import defaultdict
+
+bind_dist = 5000
+
+try:
+    tss_dict = locals()['tss_dict']
+except KeyError:
+    tss_dict = defaultdict(list)
+    tss_per_gene = defaultdict(list)
+    for i, row in pd.read_table('Reference/tss').iterrows():
+        tss_dict[row['chr'].replace('dmel_', '')].append((row['TSS_start'],
+                                                          row['gene_name']))
+        tss_per_gene[row['gene_name']].append((row['chr'].replace('dmel_', ''),
+                                              row['TSS_start']))
+
+def find_near(chrom, coord, dist, nearest_only=False):
+    out = set()
+    center = bisect([i[0] for i in chrom], coord)
+    for (coord2, gene) in reversed(chrom[:center]):
+        if coord - coord2 > dist: break
+        out.add(gene)
+        if nearest_only:
+            break
+
+    for (coord2, gene) in chrom[center:]:
+        if coord2 - coord > dist: break
+        out.add(gene)
+        if nearest_only:
+            break
+
+    return out
+
+def find_nearest_dist(chrom, coord):
+    center = bisect([i[0] for i in chrom, coord])
+    dist = abs(chrom[center][0] - coord)
+    distp1 = abs(coord - chrom[center+1][0])
+    distm1 = abs(coord - chrom[center-1][0])
+    # I should only need one of these, but I'm too lazy to read the
+    # documentation to see whether it's the +1 or -1
+    return min(dist, distp1, distm1)
+
+def get_binding_matrix(genes):
+    binding_matrix = pd.DataFrame(index=genes,
+                                  columns=tfs, data=0.0, dtype=float)
+    for tf in tfs:
+        for gene in binding_matrix.index:
+            binding_matrix.ix[gene, tf] = float(gene in has_tfs[tf])
+    binding_matrix.sort_index(axis=1, inplace=True)
+    return binding_matrix
+
+tf_files = glob('Reference/*_peaks')
+tfs = ([i.split('/')[1].split('_')[0] for i in tf_files])
+has_tfs = {}
+for tf, tf_file in zip(tfs, tf_files):
+    has_tf = set()
+    for coord in pd.read_table(tf_file).NewPeak:
+        chr, coord = coord.split(':')
+        has_tf.update(find_near(tss_dict[chr], int(coord), bind_dist))
+    has_tfs[tf] = has_tf

DistributionDifference.py

@@ -0,0 +1,265 @@
+import numpy as np
+from scipy import interpolate
+
+import collections
+import functools
+import emd
+
+class memoized(object):
+   '''Decorator. Caches a function's return value each time it is called.
+   If called later with the same arguments, the cached value is returned
+   (not reevaluated).
+   '''
+   def __init__(self, func):
+      self.func = func
+      self.cache = {}
+   def __call__(self, *args):
+      if not isinstance(args, collections.Hashable):
+         # uncacheable. a list, for instance.
+         # better to not cache than blow up.
+         return self.func(*args)
+      if tuple(*args) in self.cache:
+         return self.cache[tuple(*args)]
+      else:
+         value = self.func(*args)
+         self.cache[tuple(*args)] = value
+         return value
+   def __repr__(self):
+      '''Return the function's docstring.'''
+      return self.func.__doc__
+   def __get__(self, obj, objtype):
+      '''Support instance methods.'''
+      return functools.partial(self.__call__, obj)
+
+@memoized
+def convert_to_distribution(points):
+    """Converts the given data points to a smoothed distribution from 0-100%
+
+    """
+
+    x = np.linspace(0,1, len(points), endpoint=True)
+    f = interpolate.interp1d(x, points, kind='cubic')
+    retval =  np.cumsum(f(np.linspace(0, 1, 30, endpoint=True)).clip(0,1e5))
+    return retval / (sum(retval)+1e-10)
+
+def diff_stat(points1, points2):
+
+    dist1 = convert_to_distribution(points1)
+    dist2 = convert_to_distribution(points2)
+
+
+    normfac = np.log(max(max(points1), max(points2)) + 1)
+
+    return np.max(np.abs(dist1 - dist2)) * normfac
+
+divmat = np.zeros([0,0])
+
+
+def tang_stat(points1, points2):
+    assert len(points1) == len(points2)
+    points1 = np.array(points1 / np.mean(points1))
+    points2 = np.array(points2 / np.mean(points2))
+
+    va = np.reshape(np.repeat(points1, len(points2)), (len(points2), -1),
+                    order='C')
+    vb = np.reshape(np.repeat(points2, len(points1)), (-1, len(points2)),
+                    order='F')
+
+    global divmat
+    if np.shape(divmat) != (len(points1), len(points2)):
+        x, y = np.mgrid[0:len(points1), 0:len(points2)]
+        divmat = 1/(np.abs(x - y) + 1)
+    return np.sqrt(np.sum(np.triu((va - vb)**2 * divmat)))
+
+#    stat = 0
+#    for i in range(len(points1)):
+#        for j in range(len(points2)):
+#            stat += (points1[i] - points2[j])**2 / (np.abs(i - j)+1)
+#            if i == j:
+#                break
+#
+#    return np.sqrt(stat)
+
+def earth_mover(points1, points2):
+    xs1 = np.linspace(0,1,len(points1),
+                      endpoint=True)[np.array(np.isfinite(points1))]
+    xs2 = np.linspace(0,1,len(points2),
+                      endpoint=True)[np.array(np.isfinite(points2))]
+    points1 = points1[np.isfinite(points1)]
+    points2 = points2[np.isfinite(points2)]
+    return emd.emd(xs1, xs2,
+                   points1/np.sum(points1),
+                   points2/np.sum(points2))
+
+startswith = lambda x: lambda y: y.startswith(x)
+
+def earth_mover_multi(points1, points2):
+    dist = 0.0
+    embs = {col.split('sl')[0] for col in points1.index}
+    sums = [[],[]]
+    for emb in embs:
+        dist += earth_mover(points1.select(startswith(emb))+1e-5,
+                            points2.select(startswith(emb))+1e-5)**2
+        sums[0].append(points1.select(startswith(emb)).mean())
+        sums[1].append(points2.select(startswith(emb)).mean())
+    dist += earth_mover(np.array(sums[0]), np.array(sums[1]))
+    return dist**.5
+
+def mp_earth_mover(args):
+    i, j = args
+    return earth_mover(i, j)
+
+def mp_earth_mover_multi(args):
+    i, j = args
+    return earth_mover_multi(i, j)
+
+import progressbar as pb
+def pdist(X, metric, p=2, w=None, V=None, VI=None):
+    X = np.asarray(X, order='c')
+
+
+    s = X.shape
+    if len(s) != 2:
+        raise ValueError('A 2-dimensional array must be passed.')
+
+    m, n = s
+    dm = np.zeros((m * (m - 1) / 2,), dtype=np.double)
+
+
+
+    k = 0
+    prog = pb.ProgressBar(widgets=['calculating distances', pb.Bar(),
+                                   pb.Percentage(), pb.ETA()])
+    for i in prog(range(0, m - 1)):
+        for j in range(i + 1, m):
+            dm[k] = metric(X[i], X[j])
+            k = k + 1
+    return dm
+
+def mp_mapped(args):
+    manager, X, i, j = args
+    metric = manager.get_metric()
+    return metric(X[i], X[j])
+
+def mp_pdist(X, metric, p=2, w=None, V=None, VI=None):
+    import multiprocessing
+    from multiprocessing.managers import BaseManager
+    X = np.asarray(X, order='c')
+
+
+    s = X.shape
+    if len(s) != 2:
+        raise ValueError('A 2-dimensional array must be passed.')
+
+    m, n = s
+    dm = np.zeros((m * (m - 1) / 2,), dtype=np.double)
+
+    pool = multiprocessing.Pool(10)
+    func = globals()["mp_"+metric.__name__]
+
+    k = 0
+    prog = pb.ProgressBar(widgets=['calculating distances', pb.Bar(),
+                                   pb.Percentage(), pb.ETA()])
+    for i in prog(range(0, m - 1)):
+        ks = np.arange(k, k + m - i - 1)
+        inputs = [(X[i], X[j]) for j in range(i+1, m)]
+        dm[ks] = pool.map(func, inputs)
+        k  = ks[-1] + 1
+    return dm
+
+def mp_pandas_pdist(X, metric, p=2, w=None, V=None, VI=None):
+    import multiprocessing
+
+    s = X.shape
+    if len(s) != 2:
+        raise ValueError('A 2-dimensional array must be passed.')
+
+    m, n = s
+    dm = np.zeros((m * (m - 1) / 2,), dtype=np.double)
+
+    pool = multiprocessing.Pool()
+    if metric.__name__.endswith('multi'):
+        func = globals()["mp_"+metric.__name__]
+    else:
+        func = globals()["mp_"+metric.__name__+"_multi"]
+
+    k = 0
+    prog = pb.ProgressBar(widgets=['calculating distances', pb.Bar(),
+                                   pb.Percentage(), pb.ETA()])
+    for i in prog(range(0, m - 1)):
+        ks = np.arange(k, k + m - i - 1)
+        inputs = [(X.ix[i], X.ix[j]) for j in range(i+1, m)]
+        dm[ks] = pool.map(func, inputs)
+        k  = ks[-1] + 1
+    return dm
+
+def pandas_pdist(X, metric, p=2, w=None, V=None, VI=None):
+    s = X.shape
+    if len(s) != 2:
+        raise ValueError('A 2-dimensional array must be passed.')
+
+    m, n = s
+    dm = np.zeros((m * (m - 1) / 2,), dtype=np.double)
+
+    if metric.__name__.endswith('multi'):
+        func = globals()["mp_"+metric.__name__]
+    else:
+        func = globals()["mp_"+metric.__name__+"_multi"]
+
+    k = 0
+    prog = pb.ProgressBar(widgets=['calculating distances', pb.Bar(),
+                                   pb.Percentage(), pb.ETA()])
+    for i in prog(range(0, m - 1)):
+        ks = np.arange(k, k + m - i - 1)
+        inputs = [(X.ix[i], X.ix[j]) for j in range(i+1, m)]
+        print(len(inputs))
+        print(ks)
+        dm[ks] = map(func, inputs)
+        k  = ks[-1] + 1
+    return dm
+
+
+if __name__ == "__main__":
+    import pandas as pd
+    import matplotlib.pyplot as mpl
+
+    zld_exp = pd.read_table('analysis/summary.tsv', index_col=0).sort_index()
+    wt_exp = pd.read_table('prereqs/WT5.53_summary.tsv', index_col=0).sort_index()
+
+    zld_bind = pd.read_table('journal.pgen.1002266.s005.xls', skiprows=1)
+    zld_bind.TSS_gene = zld_bind.TSS_gene.apply(str.strip)
+    by_gene = zld_bind.groupby('TSS_gene')
+
+
+    types = {'Intergenic':'N', 'Intronic':'I', 'Promoter':'P', 'UTR5':'5',
+             'CDS':'C', 'UTR3':'3'}
+
+    zld_comp = zld_exp.select(lambda x: 'cyc14A' in x, axis=1)
+    wt_comp = wt_exp.select(lambda x: 'cyc14A' in x, axis=1)
+
+    diff_col = pd.Series(index=zld_comp.index)
+    for gene in wt_exp.index:
+        assert gene in zld_exp.index
+        diff_col[gene] = diff_stat(zld_comp.ix[gene], wt_comp.ix[gene])
+
+    zld_exp['diff_col'] = diff_col
+    wt_exp['diff_col'] = diff_col
+
+    zld_exp.sort(column='diff_col', ascending=False, inplace=True)
+    wt_exp.sort(column='diff_col', ascending=False, inplace=True)
+
+    zld_fig_genes = zld_exp.select(lambda x: '14A' in x or '11' in x, axis=1)
+    wt_fig_genes = wt_exp.select(lambda x: '14A' in x or '11' in x, axis=1)
+
+    zld_fig_genes = zld_fig_genes[wt_fig_genes.max(axis=1) > 10][:120]
+    wt_fig_genes = wt_fig_genes[wt_fig_genes.max(axis=1) > 10][:120]
+
+    assert (zld_fig_genes.index == wt_fig_genes.index).all()
+    import PlotUtils
+
+    PlotUtils.svg_heatmap((wt_fig_genes, zld_fig_genes),
+                            'analysis/results/cyc13diff.svg',
+                            norm_rows_by=wt_fig_genes.max(axis=1),
+                            draw_row_labels=True,
+                            cmap = (mpl.cm.Blues, mpl.cm.Reds),
+                            box_size=15, total_width=150)