Commits

Marinka Zitnik committed 04516b0

saving graphs to subdirectories, produce graphs for time subgroups, invalid patches stats, plate csv

Comments (0)

Files changed (3)

     Combine these object measurements with those of the previous object?:No
     File name:\\g<Movie>-TrackedPatch.csv
     Use the object name for the file name?:No
+    Data to export:Cell
+    Combine these object measurements with those of the previous object?:No
+    File name:\\g<Movie>-TrackedCells.csv
+    Use the object name for the file name?:No

yp/cell_filter.py

 
 def main(argv):
     if len(sys.argv) < 4:
-        sys.exit("Usage: %s <DefaultOUT_Patch.csv path> <DefaultOUT_FilteredCell.csv path> <res dir>" % sys.argv[0])
+        sys.exit("Usage: python %s <DefaultOUT_Patch.csv path> <DefaultOUT_FilteredCell.csv path> <res dir>" % sys.argv[0])
     if not os.path.exists(sys.argv[1]):
         sys.exit("ERROR: Patch data file %s was not found!" % sys.argv[1])
     if not os.path.exists(sys.argv[2]):

yp/patch_tracking.py

     
     :param filepath_patch: Full path to the CP exported data (module ExportToSpreadsheet) about patches. 
     :type filepath_patch: `str`
+    :param filepath_cell: Full path to the CP exported data about cells. 
+    :type filepath_cell: `str`
     :param res_dir: Full path to directory with results.
     :type res_dir: `str`
     """
     
     namesc = fc.readline().strip().split(",")
     print "No. columns %d: %s ..." % (len(namesc), namesc[:min(3, len(namesc))])
-    obj = "TrackObjects_Label_20"
-    assert obj in namesc, "ERROR: Feature '%s' not found." % obj
+    obj_c = "TrackObjects_Label_20"
+    assert obj_c in namesc, "ERROR: Feature '%s' not found." % obj_c
     
     io2t = {str(i): {} for i in xrange(1, 51)} #image x objectNumber -> trackedObjects_Label
     for line in fc:
         dline = dict(zip(namesc, line.strip().split(",")))
-        io2t[dline["ImageNumber"]][dline["ObjectNumber"]] = dline[obj]
+        io2t[dline["ImageNumber"]][dline["ObjectNumber"]] = dline[obj_c]
     fc.close()
     ###
 
     patch2dist = defaultdict(list) #distance traveled by patch from first match (>=1) to end of its lifetime (<50)
     patch2startstop = {} #patch start stop locations
     gl_distance = {} #patch minimum distance to perimeter point
+    invalid_times = []
     
     done = set()
     for k, tracked_patch in enumerate(dbp[obj]):
-        if tracked_patch in done:
-            continue
-        done.add(tracked_patch)        
-        if tracked_patch in it2d["50"]:
-            #za dolociti premik od membrane so v redu krpice, ki izginejo znotraj opazovanega casa
-            continue 
+        if tracked_patch in done: continue
+        done.add(tracked_patch)
+        img = dbp["ImageNumber"][k]
+        #save invalid patches
+        if tracked_patch in it2i["1"] or tracked_patch in it2i["50"]:
+            invalid_times.append(np.sum([1 for img1 in xrange(int(img), 51) if tracked_patch in it2d[str(img1)]]))
+        #patch move relative to cell membrane: valid patches disappear during observation time 
+        if tracked_patch in it2d["50"]: continue
         time_dist = []
-        img = dbp["ImageNumber"][k]
         time_dist.append((int(img), it2d[img][tracked_patch]))
         patch2startstop[tracked_patch] = [it2loc[img][tracked_patch]]*2
     
         d = zip(*time_dist)[1]
         if not np.isnan(d).any(): 
             gl_distance[tracked_patch] = d
-            
-        if tracked_patch in it2i["1"]:
-            #za dolocitev zivljenjskega casa pridejo v postev kripice, ki se zacnejo in koncajo znotraj opazovanega casa
-            continue
+
+        #patch lifetime and intensity: valid patches appear and disappear during observation time
+        if tracked_patch in it2i["1"]: continue
         time_inten = []
         yerr = []
         img = dbp["ImageNumber"][k]
         iimg = int(img)
         time_inten.append((iimg, it2i[img][tracked_patch][0]))
         yerr.append(it2i[img][tracked_patch][1])
-        global_inten_start[0].append(it2i[img][tracked_patch][0])
-        tmp_val = [it2i[img][tracked_patch][0]]
         for img1 in xrange(iimg+1, 50):
             if not tracked_patch in it2i[str(img1)]:
                 break
             time_inten.append((img1, it2i[str(img1)][tracked_patch][0]))
             yerr.append(it2i[str(img1)][tracked_patch][1])
-            global_inten_start[img1-iimg].append(it2i[str(img1)][tracked_patch][0])
-            tmp_val.append(it2i[str(img1)][tracked_patch][0])
           
         if len(time_inten)>=5:
+            tmp_val = zip(*time_inten)[1]
             gis.append(tmp_val)  
+            for img_idx, inten in enumerate(tmp_val): 
+                global_inten_start[img_idx].append(inten)
+            _align_max(global_inten_max, tmp_val)
              
-            global_inten_max[0].append(it2i[img][tracked_patch][0])
-            max_idx = tmp_val.index(max(tmp_val))
-            global_inten_max[0].append(tmp_val[max_idx])
-            for i in xrange(1, len(tmp_val)-max_idx):
-                global_inten_max[i].append(tmp_val[max_idx+i])
-            for i in xrange(1, max_idx):
-                global_inten_max[-i].append(tmp_val[max_idx-i])
-            
         gl_lifetime.append((tracked_patch, len(time_inten)))
         
         #do not plot and save membrane distance if object lifetime in less than 5 frames
     gl_patches, gl_times = zip(*gl_lifetime)  
     patch2startstop = {p: math.sqrt( (v[0][0]-v[1][0])**2 + (v[0][1]-v[1][1])**2) for p, v in patch2startstop.iteritems()}       
     
+    feat_vals = []
+    #saving basic stats to file
     fname = res_dir + "/patches.txt"
     fres = open(fname, "w")
     fres.write("Cell file: %s\n" % filepath_cell)
     fres.write("Patch file: %s\n" % filepath_patch)
     fres.write("=====================================\n")
-    fres.write("N(patches all tracked from 1st including 50th frame): %d\n" % len(done))
-    fres.write("Avg(N(cells per frame)): %5.3f" % np.mean([len(io2t[img]) for img in io2t]))    
-    fres.write("N(cells in 1st frame): %d" % len(io2t["1"]))
-    fres.write("N(cells in 50th frame): %d" % len(io2t["50"]))
+    feat_vals.append(np.mean([len(io2t[img]) for img in io2t]))
+    fres.write("Avg(N(cell per frame)): %5.3f\n" % feat_vals[-1])
+    feat_vals.append(len(io2t["1"]))    
+    fres.write("N(cell in 1st frame): %d\n" % feat_vals[-1])
+    feat_vals.append(len(io2t["50"]))
+    fres.write("N(cell in 50th frame): %d\n" % feat_vals[-1])
+    fres.write("=====================================\n")
+    feat_vals.append(len(done))
+    fres.write("N(tracked patch): %d\n" % feat_vals[-1])
+    feat_vals.append(len(done)-len(invalid_times))
+    fres.write("N(valid patch for lifetime determination): %d\n" % feat_vals[-1])
+    feat_vals.append(len(invalid_times))
+    fres.write("N(invalid patch=exist in 1st or 50th frame): %d\n" % feat_vals[-1])
+    for step in xrange(1, int(np.max(invalid_times)), 2):
+        feat_vals.append(np.sum([1 for e in invalid_times if e>=step]))
+        fres.write("\tN(invalid patch with >=%d measur.): %d\n" % (step, feat_vals[-1]))
+    fres.write("=====================================\n")
     for k in xrange(0, 46, 5):
         ptch = [p for p, e in gl_lifetime if e >= k and e < k+5]
-        fres.write("[%d,%d)\tN: %4d\n" % (k, k+5, len(ptch)))
-        fres.write("\tAvg(Avg-distance-to-clos.-point-on-perim.-par.-cell for patch in Patches): %5.3f\n" % np.mean([np.mean(v) for p, v in gl_distance.iteritems() if p in ptch]))
-        fres.write("\tAvg(Max-distance-to-clos.-point-on-perim.-par.-cell for patch in Patches): %5.3f\n" % np.mean([np.max(v) for p, v in gl_distance.iteritems() if p in ptch]))
-        fres.write("\tAvg(Trajectory-length for patch in Patches): %5.3f\n" % np.mean([np.sum(d) for p, d in patch2dist.iteritems() if p in ptch]))
-        fres.write("\tAvg(Avg-move-per-frame  for patch in Patches): %5.3f\n" % np.mean([np.mean(d) for p, d in patch2dist.iteritems() if p in ptch]))
-        fres.write("\tAvg(Max-move-per-frame for patch in Patches): %5.3f\n" % np.mean([np.max(d) for p, d in patch2dist.iteritems() if p in ptch]))
-        fres.write("\tAvg(Start-stop-distance for patch in Patches): %5.3f\n" % np.mean([v for p, v in patch2startstop.iteritems() if p in ptch]))
+        feat_vals.append(len(ptch))
+        fres.write("[%d,%d)\tN: %4d\n" % (k, k+5, feat_vals[-1]))
+        feat_vals.append(np.mean([np.mean(v) for p, v in gl_distance.iteritems() if p in ptch]))
+        fres.write("\tAvg(Avg-distance-to-clos.-point-on-perim.-par.-cell for patch in Patches): %5.3f\n" % feat_vals[-1])
+        feat_vals.append(np.mean([np.max(v) for p, v in gl_distance.iteritems() if p in ptch]))
+        fres.write("\tAvg(Max-distance-to-clos.-point-on-perim.-par.-cell for patch in Patches): %5.3f\n" % feat_vals[-1])
+        feat_vals.append(np.mean([np.sum(d) for p, d in patch2dist.iteritems() if p in ptch]))
+        fres.write("\tAvg(Trajectory-length for patch in Patches): %5.3f\n" % feat_vals[-1])
+        feat_vals.append(np.mean([np.mean(d) for p, d in patch2dist.iteritems() if p in ptch]))
+        fres.write("\tAvg(Avg-move-per-frame  for patch in Patches): %5.3f\n" % feat_vals[-1])
+        feat_vals.append(np.mean([np.max(d) for p, d in patch2dist.iteritems() if p in ptch]))
+        fres.write("\tAvg(Max-move-per-frame for patch in Patches): %5.3f\n" % feat_vals[-1])
+        feat_vals.append(np.mean([v for p, v in patch2startstop.iteritems() if p in ptch]))
+        fres.write("\tAvg(Start-stop-distance for patch in Patches): %5.3f\n" % feat_vals[-1])
         
     fres.write("=====================================\n")
-    fres.write("Mean life. (all): %5.3f\n" % np.mean(gl_times))
+    feat_vals.append(np.mean(gl_times))
+    fres.write("Mean life. (all): %5.3f\n" % feat_vals[-1])
     mx = np.max(gl_times)
     vals = np.array(gl_times)
     for step in xrange(1, mx+1, 2):
-        fres.write("Mean life. (>=%d): %5.3f\n" % (step, np.mean(vals[vals >= step])))
+        feat_vals.append(np.mean(vals[vals >= step]))
+        fres.write("Mean life. (>=%d): %5.3f\n" % (step, feat_vals[-1]))
     fres.write("===============Total=================\n")
-    fres.write("Avg(Avg-distance-to-clos.-point-on-perim.-par.-cell for patch in Patches): %5.3f\n" % np.mean([np.mean(v) for v in gl_distance.values()]))
-    fres.write("Avg(Max-distance-to-clos.-point-on-perim.-par.-cell for patch in Patches): %5.3f\n" % np.mean([np.max(v) for v in gl_distance.values()]))
-    fres.write("Avg(Trajectory-length for patch in Patches): %5.3f\n" % np.mean([np.sum(d) for d in patch2dist.values()]))
-    fres.write("Avg(Avg-move-per-frame  for patch in Patches): %5.3f\n" % np.mean([np.mean(d) for d in patch2dist.values()]))
-    fres.write("Avg(Max-move-per-frame for patch in Patches): %5.3f\n" % np.mean([np.max(d) for d in patch2dist.values()]))
-    fres.write("Avg(Start-stop-distance for patch in Patches): %5.3f\n" % np.mean(patch2startstop.values()))
+    feat_vals.append(np.mean([np.mean(v) for v in gl_distance.values()]))
+    fres.write("Avg(Avg-distance-to-clos.-point-on-perim.-par.-cell for patch in Patches): %5.3f\n" % feat_vals[-1])
+    feat_vals.append(np.mean([np.max(v) for v in gl_distance.values()]))
+    fres.write("Avg(Max-distance-to-clos.-point-on-perim.-par.-cell for patch in Patches): %5.3f\n" % feat_vals[-1])
+    feat_vals.append(np.mean([np.sum(d) for d in patch2dist.values()]))
+    fres.write("Avg(Trajectory-length for patch in Patches): %5.3f\n" % feat_vals[-1])
+    feat_vals.append(np.mean([np.mean(d) for d in patch2dist.values()]))
+    fres.write("Avg(Avg-move-per-frame  for patch in Patches): %5.3f\n" % feat_vals[-1])
+    feat_vals.append(np.mean([np.max(d) for d in patch2dist.values()]))
+    fres.write("Avg(Max-move-per-frame for patch in Patches): %5.3f\n" % feat_vals[-1])
+    feat_vals.append(np.mean(patch2startstop.values()))
+    fres.write("Avg(Start-stop-distance for patch in Patches): %5.3f\n" % feat_vals[-1])
     fres.write("=====================================\n")
-    
     fres.close()
     print "Processing info. saved to file: %s" % fname
+
+    #adding new observation to plate specific file
+    tmp = os.path.dirname(os.path.dirname(res_dir))   
+    plate_name = tmp.split(os.path.sep)[-1][:-11]
+    plate_fname = tmp + "%s_patch_tracking.csv" % plate_name
+    if not os.path.exists(plate_fname):
+        fplate = open(plate_fname, "w")
+        fplate.write("%s\n" % ",".join(["feat_%d" % i for i in xrange(len(feat_vals))]))
+        fplate.close()
     
-    _plot_gl_int_profile_start(global_inten_start, gis, res_dir)
-    _plot_gl_int_profile_max(global_inten_max, gis, res_dir)
-    _plot_gl_int_profile_end(gis, res_dir)
+    fplate = open(plate_fname, "a")
+    fplate.write("%s\n" % ",".join(["%5.4f" % v for v in feat_vals]))
+    fplate.close()
+    print "New observation added to file: %s" % plate_fname
+    
+    #global intensity profile per time-series group
+    for k in xrange(0, 46, 5):
+        type_gis = [e for e in gis if len(e) >= k and len(e) < k+5]
+        if len(type_gis) == 0: continue
+        typ = "%s-%s" % (k, k+5)
+        mx = np.max([len(e) for e in type_gis])
+        type_inten_start = {i: [e[i] for e in type_gis if len(e)>i] for i in xrange(mx+1)}
+        type_inten_max = defaultdict(list)
+        for e in type_gis:
+            _align_max(type_inten_max, e) 
+        _plot_gl_int_profile_start(type_inten_start, type_gis, res_dir + "/%s" % typ, typ)
+        _plot_gl_int_profile_max(type_inten_max, type_gis, res_dir + "/%s" % typ, typ)
+        _plot_gl_int_profile_end(type_gis, res_dir + "/%s" % typ, typ)
+    
+    #global intensity profile for all time-series groups
+    _plot_gl_int_profile_start(global_inten_start, gis, res_dir, "global")
+    _plot_gl_int_profile_max(global_inten_max, gis, res_dir, "global")
+    _plot_gl_int_profile_end(gis, res_dir, "global")
+    
+def _align_max(inten_max, tmp_val):
+    max_idx = tmp_val.index(max(tmp_val))
+    inten_max[0].append(tmp_val[max_idx])
+    for i in xrange(1, len(tmp_val)-max_idx):
+        inten_max[i].append(tmp_val[max_idx+i])
+    for i in xrange(1, max_idx+1):
+        inten_max[-i].append(tmp_val[max_idx-i])
     
 def _plot_profile(time_dist, time_inten, yerr, tracked_patch, res_dir, lim_dist, lim_inten):
     """Compute and save distance-intensity profile of a patch and its spline approximation."""
     plt.savefig(path + "/intensity_dist_patch%s.pdf" % tracked_patch)
     plt.close()
             
-def _plot_gl_int_profile_start(global_inten_start, gis, res_dir):
+def _plot_gl_int_profile_start(global_inten_start, gis, res_dir, type):
     """Compute global plot of patch intensity profiles aligned by start of lifetime."""
     plt.clf()
     plt.figure(1, figsize = (8,8))
     plt.setp(ax_histx.get_xticklabels(), visible = False)
     ax_histx.set_xlim(ax_scatter.get_xlim())
     
-    plt.savefig(res_dir + "/global_intensity_profile_start.pdf")
+    if not os.path.exists(res_dir):
+        os.makedirs(res_dir)
+    plt.savefig(res_dir + "/%s_intensity_profile_start.pdf" % type)
     plt.close()
     
-def _plot_gl_int_profile_max(global_inten_max, gis, res_dir):
+def _plot_gl_int_profile_max(global_inten_max, gis, res_dir, type):
     """
     Compute global plot of patch intensity profiles aligned by max intensity and 
     its spline approximation.
         ax.plot(x_sp[j:j+2], y_sp[j:j+2], 'o-', color = "red", linewidth = 2, alpha = 0.2 + (weights[j] + weights[j+1])/2.)
     
     ax.set_ylabel("Mean intensity")
-    ax.set_xlabel("Intensity alignment centered at maximum")
+    ax.set_xlabel("Lifetime centered at maximum intensity")
     ax.set_xlim((x_sp[0]-1, x_sp[-1]+1))
     
-    plt.savefig(res_dir + "/global_intensity_profile_max.pdf")
+    if not os.path.exists(res_dir):
+        os.makedirs(res_dir)
+    plt.savefig(res_dir + "/%s_intensity_profile_max.pdf" % type)
     plt.close()
     
-def _plot_gl_int_profile_end(gis, res_dir):
+def _plot_gl_int_profile_end(gis, res_dir, type):
     """Compute global plot of patch intensity profiles aligned by end of lifetime."""
     plt.clf()
     fig = plt.figure(1, figsize = (8,8))
     for v in gis:
         ax.plot((-1)*np.arange(len(v))[::-1], v, 'o-', alpha = 0.5)
     ax.set_ylabel("Mean intensity")
-    ax.set_xlabel("Intensity alignment aligned by end")
+    ax.set_xlabel("Alignment by end of lifetime")
     ax.set_xlim((-np.max([len(v) for v in gis])-1, 1))
     
-    plt.savefig(res_dir + "/global_intensity_profile_end.pdf")
+    if not os.path.exists(res_dir):
+        os.makedirs(res_dir)
+    plt.savefig(res_dir + "/%s_intensity_profile_end.pdf" % type)
     plt.close()
 
 def main(argv):
     if len(sys.argv) < 3:
-        sys.exit("Usage: %s <DefaultOUT_TrackedPatch.csv path> <DefaultOUT_TrackedCell.csv path> <res dir>" % sys.argv[0])
+        sys.exit("Usage: python %s <DefaultOUT_TrackedPatch.csv path> <DefaultOUT_TrackedCell.csv path> <res dir>" % sys.argv[0])
     if not os.path.exists(sys.argv[1]):
         sys.exit("ERROR: Patch data file %s was not found!" % sys.argv[1])
     if not os.path.exists(sys.argv[2]):