Cosmetic changes

growthAnalysis/growth_control.py

 import os
 import math
-from itertools import permutations, product, repeat, takewhile, groupby
+from itertools import permutations, product, repeat, takewhile
 from functools import partial, reduce
 from operator import add, itemgetter
 from collections import defaultdict
 from subprocess import call
-from typing import Dict, List, Tuple
+from typing import Dict, List, Tuple, Set
 
 import ast
 import astor
@@ -18,6 +18,7 @@ from scipy.signal import argrelmax
 from scipy.stats import gaussian_kde
 from scipy.stats import ttest_ind
 from more_itertools import interleave_longest
+from mpl_toolkits.mplot3d import Axes3D
 
 import common.statesN as f
 import common.seg as seg
@@ -59,14 +60,16 @@ sts_ts = f.states_per_t()
 fst = itemgetter(0)
 snd = itemgetter(1)
 
+
 def set_plot_params(fontsize=11):
-    sns.set_style("whitegrid", {'grid.linestyle':'--'})
+    sns.set_style("whitegrid", {'grid.linestyle': '--'})
     plt.rc('xtick', labelsize=fontsize)
     plt.rc('ytick', labelsize=fontsize)
     plt.rc('axes', labelsize=fontsize)
     params = {'legend.fontsize': fontsize}
     plt.rcParams.update(params)
 
+
 class BExpr2():
 
     def __init__(self, e_str):
@@ -134,19 +137,40 @@ class GPairRegions():
 
         for c in self.ts:
             if seg.evalB(getTopExpr(he), c.exprs):
-                d[c.cid] = 100.0
+                d[c.cid] = 1.0
             else:
                 d[c.cid] = 0.5
 
         for c in self.ts:
             if seg.evalB(getTopExpr(le), c.exprs):
-                d[c.cid] = -10.0
+                d[c.cid] = 0.0
+
+        f = plot_ts_q_(self.ts, d,
+                        lb="{g1}-{g2}".format(g1=self.g1, g2=self.g2),
+                        bounds=(0.1, 0.9), txt=txt)
+
+        return f, d
+
+    def get_regions(self) -> Tuple[Set[int], Set[int]]:
+        grs_state = get_grs_state()
+        sts = sts_ts[self.t]
+        sts_e1 = get_sts_expr(self.ts, self.e1, sts)
+        sts_e2 = get_sts_expr(self.ts, self.e2, sts)
+        gr_e1 = meanor0(reduce(add, [grs_state[st] for st in sts_e1], []))
+        gr_e2 = meanor0(reduce(add, [grs_state[st] for st in sts_e2], []))
 
-        f = plot_ts_q(self.ts, d,
-                      lb="{g1}-{g2}".format(g1=self.g1, g2=self.g2),
-                      bounds=(0, 1), txt=txt)
+        he, le = [e for e, gr in sorted([(self.e1, gr_e1), (self.e2, gr_e2)],
+                                        key=snd, reverse=True)]
+
+        r1 = list()
+        r2 = list()
+        for c in self.ts:
+            if seg.evalB(getTopExpr(he), c.exprs):
+                r1.append(c.cid)
+            elif seg.evalB(getTopExpr(le), c.exprs):
+                r2.append(c.cid)
 
-        return f
+        return set(r1), set(r2)
 
     def plot_distr(self):
         set_plot_params(fontsize=12)
@@ -178,7 +202,9 @@ class GPairRegions():
         sns.distplot(l_grs, bins=10, color=cls[0], kde=True, ax=ax,
                      label=astor.to_source(le).replace("'", ""))
         ax.set_xlabel(r'$\mu$m/h')
-        plt.legend(fontsize=10, loc='upper right', frameon=False)
+        #plt.legend(fontsize=10, loc='upper right', frameon=False)
+        gns = "{g1}-{g2}".format(g1=self.g1, g2=self.g2)
+        print(gns, gratio(h_grs, l_grs))
         plt.savefig("{g1}-{g2}-regions.png".format(g1=self.g1, g2=self.g2), dpi=300)
         plt.show()
 
@@ -200,15 +226,18 @@ def groupByRegion(t, ts, es: List[str]) -> Dict[int, List[str]]:
 
     return res
 
+
 def mean_gratio(ps):
     return meanor0([abs(v) for p, v in ps])
 
+
 def mean_gratio_(ps):
     return meanor0([v for p, v in ps])
 
+
 def groupByRegionGPairs(t: int,
                         ts: STissue,
-                        ress: Dict[Tuple[str, str], float],
+                        ress: List[Tuple[Tuple[str, str], float]],
                         n=6,
                         prc=0.5):
     d = defaultdict(list)
@@ -229,6 +258,7 @@ def groupByRegionGPairs(t: int,
 
     return [(i+1, v) for i, (k, v) in zip(range(len(d_)), d_)]
 
+
 def plotRegions(t, ts, d):
     fns = list()
     for rid, ps in d:
@@ -244,6 +274,7 @@ def plotRegions(t, ts, d):
 
     return "t{t}_regions.png".format(t=t)
 
+
 def cluster_states(grs_state: Dict[int, List[float]]):
     mgrs = {st:np.mean(grs) for st, grs in grs_state.items()}
     g = gaussian_kde(np.array(list(mgrs.values())))
@@ -261,6 +292,7 @@ def cluster_states(grs_state: Dict[int, List[float]]):
 
     return g, xsm, ysm, cls
 
+
 def plot_clustering(g, xsm, ysm):
     xs = np.arange(-0.05, 0.2, 0.001)
     ys = g.evaluate(xs)
@@ -273,6 +305,7 @@ def plot_clustering(g, xsm, ysm):
     plt.savefig("clusters_grates_states.png", dpi=300)
     plt.show()
 
+
 def get_grs_state() -> Dict[int, List[float]]:
     tss, linss = lin.mkSeries1(d="../data/FM1/tv/",
       	                       dExprs="../data/geneExpression/",
@@ -287,11 +320,14 @@ def get_grs_state() -> Dict[int, List[float]]:
 
     return grates_state
 
+
 def ge(s: str): return ast.parse(s)
 
+
 def construct_expr(nms, fns):
     return ge(" ".join(interleave_longest(nms, fns)))
 
+
 def mk_bexprs(k, geneNms):
     not_genes = ["not {gn}".format(gn=gn) for gn in geneNms]
     not_gene_nms = geneNms + not_genes
@@ -305,17 +341,20 @@ def mk_bexprs(k, geneNms):
 
     return
 
+
 def search_bexprs(ts, k, obj):
     ress = dict()
     for i, bexpr in enumerate(mk_bexprs(k, geneNms)):
-        #cids = ts.filterGBExpr(bexpr)
+        # cids = ts.filterGBExpr(bexpr)
         ress[bexpr] = obj(cids2=bexpr)
 
     return ress
 
+
 def hash_expr(ts, e):
     return np.sum(ts.filterGBExpr(e))
 
+
 def get_cids_highg(ts):
     sts = f.statesI()
     grs_state = get_grs_state()
@@ -326,10 +365,12 @@ def get_cids_highg(ts):
 
     return cids
 
+
 def get_cids_sts(ts, sts):
     st_genes = f.statesI()
     return reduce(add, [ts.filterGs(st_genes[st]) for st in sts])
 
+
 def bacc(ts, cids1, cids2):
     scids1 = set(cids1)
     scids2 = set(cids2)
@@ -354,18 +395,21 @@ def bacc(ts, cids1, cids2):
 
     return 0.5*(tpr + tnr)
 
+
 def meanor0(xs):
     if xs:
         return np.mean(xs)
     else:
         return 0.0
 
+
 def avg_grate(grates, cids2):
     grs_cids = [grates.get(cid, None) for cid in cids2]
     grs_cids_ = [gr for gr in grs_cids if gr]
 
     return meanor0(grs_cids_)
 
+
 def gratio(xs, xs1):
     if not xs or not xs1:
         return 0.0
@@ -375,6 +419,7 @@ def gratio(xs, xs1):
 
     return (m1 - m2) / (m1 + m2)
 
+
 def states_d(ts, cids2, sts, grs_state_m):
     sts_envs = f.statesI_env(seg.geneNms)
     e = getTopExpr(cids2)
@@ -385,6 +430,7 @@ def states_d(ts, cids2, sts, grs_state_m):
 
     return gratio(mgrs_sts_t, mgrs_sts_f)
 
+
 def go_search(t, k):
     tss, linss = lin.mkSeries()
     lin.filterL1_st(tss)
@@ -560,8 +606,9 @@ def vis_stripplot_hue2(d, i):
 
     plt.show()
 
+
 def cellsToDat(fn, cells):
-    #so we can use newman
+    # so we can use newman
     nvars = len(cells[0].exprs.keys()) + 5
     ncells = len(cells)
     ntpoints = 1
@@ -572,9 +619,11 @@ def cellsToDat(fn, cells):
     with open(fn, 'w+') as fout:
         fout.write("\n".join([header] + [cell.toOrganism() for cell in cells]))
 
+
 def plot_ts_binary(ts_, cids, lb="region", txt=""):
     from copy import deepcopy
     ts = deepcopy(ts_)
+
     def toNewmanFn(fn):
         (fn, ext) = os.path.splitext(fn)
         return fn + "000" + ".tif"
@@ -631,7 +680,54 @@ def plot_ts_binary(ts_, cids, lb="region", txt=""):
     call(["open",
           "{lb}.png".format(lb=lb)])
 
-def plot_ts_q(ts_, d, lb="vals", bounds=(-0.03, 0.05), txt="", interactive=False, size=100, cm=3):
+
+def plot_ts_q_(ts, d, lb="vals", bounds=(0, 1), txt=""):
+    vmin, vmax = bounds
+
+    xs = [c.pos.x for c in ts]
+    ys = [c.pos.y for c in ts]
+    zs = [c.pos.z for c in ts]
+    cs = [d[c.cid] for c in ts]
+
+    fig = plt.figure(figsize=(1.5, 1.5))
+    ax = fig.add_subplot(111, projection='3d')
+    ax.view_init(elev=104, azim=-89)
+    ax.scatter(xs, ys, zs, c=cs, cmap="coolwarm", alpha=0.9, s=7, vmin=vmin,
+               vmax=vmax, edgecolors='black', linewidths=0.5)
+    ax.set_title(txt)
+    plt.axis('off')
+
+    fout = "{lb}.png".format(lb=lb)
+    plt.savefig(fout, dpi=300, bbox_inches='tight')
+
+    plt.show()
+
+
+def plot_distrs(grs, r1, r2, lb):
+    grs1 = [grs.get(cid, np.nan) for cid in r1]
+    grs2 = [grs.get(cid, np.nan) for cid in r2]
+
+    grs1_ = [gr for gr in grs1 if not math.isnan(gr)]
+    grs2_ = [gr for gr in grs2 if not math.isnan(gr)]
+
+    cls = sns.color_palette("coolwarm")
+
+    fig = plt.figure(figsize=(3.5, 3.7))
+    ax = fig.add_subplot('111')
+
+    ax.set_xlim(-0.05, 0.18)
+
+    sns.distplot(grs1_, bins=10, color=cls[-1], kde=True, ax=ax,
+                 label="region1")
+    sns.distplot(grs2_, bins=10, color=cls[0], kde=True, ax=ax,
+                 label="region2")
+    ax.set_xlabel(r'$\mu$m/h'.format())
+    plt.savefig("zone{i}.png".format(i=lb), dpi=300, bbox_inches='tight')
+    plt.show()
+
+
+def plot_ts_q(ts_, d, lb="vals", bounds=(0, 1), txt="",
+              interactive=False, size=100, cm=1):
     from copy import deepcopy
     ts = deepcopy(ts_)
 
@@ -1154,3 +1250,66 @@ def go():
 
     g1.plot()
     g1.plot_distr()
+
+
+def go_():
+    tss, linss = lin.mkSeriesIm0(dataDir="../data/",
+                                 ft=lambda t: t in {10, 40, 96, 120, 132})
+    lin.filterL1_st(tss)
+
+    t = 132
+    ts = tss[t]
+    g = GPairRegions(t, ts, 'ETTIN', 'LFY')  # zone 1
+    g1 = GPairRegions(t, ts, 'SEP1', 'LFY')  # zone 2
+    g2 = GPairRegions(t, ts, 'AP1', 'LFY')   # zone 3
+
+    g.plot()
+    g.plot_distr()
+
+    g1.plot()
+    g1.plot_distr()
+
+    g2.plot()
+    g2.plot_distr()
+
+def go_t120():
+    tss, linss = lin.mkSeriesIm0(dataDir="../data/",
+                                 ft=lambda t: t in {10, 40, 96, 120, 132})
+    lin.filterL1_st(tss)
+
+    t = 120
+    ts = tss[t]
+    g = GPairRegions(t, ts, 'ETTIN', 'LFY')  # zone 1
+    g1 = GPairRegions(t, ts, 'AP1', 'LFY')   # zone 2
+
+    g.plot()
+    g.plot_distr()
+
+    g1.plot()
+    g1.plot_distr()
+
+
+def go__():
+    tss, linss = lin.mkSeriesIm0(dataDir="../data/",
+                                 ft=lambda t: t in {10, 40, 96, 120, 132})
+    lin.filterL1_st(tss)
+
+    t = 132
+    ts = tss[t]
+
+    r1 = [c.cid for c in ts if c.exprs['LFY']]
+    r2 = [c.cid for c in ts if c.exprs['CUC1_2_3']]
+
+    d = dict()
+    for c in ts:
+        if c.cid in r1:
+            d[c.cid] = 1.0
+        elif c.cid in r2:
+            d[c.cid] = 0.0
+        else:
+            d[c.cid] = 0.5
+
+    grates = grs.grates_avg_cons(tss, linss)
+
+    plot_ts_q_(ts, d)
+    plot_distrs(grates[t], r1, r2, lb="boundary")