Merge pull request #2 from mmolari/feat/improve-clips

Feat/improve clips

Author: mmolari

nextflow.config

@@ -43,7 +43,7 @@ profiles {
             withLabel: 'q6h_1core' {
                 time = '06:00:00'
                 cpus = 1
-                memory = '16G'
+                memory = '32G'
                 clusterOptions = '--qos=6hours'
             }
 

scripts/plot_clips.py

@@ -128,17 +128,17 @@ def plot_joint_distr(df):
     bins = np.arange(max([f.max(), r.max()]) + 2) - 0.5
     norm = mpl.colors.LogNorm()
     m = ax.hist2d(f, r, bins=bins, norm=norm)
-    ax.set_xlabel("n. insertions forward")
-    ax.set_ylabel("n. insertions reverse")
+    ax.set_xlabel("n. clips forward")
+    ax.set_ylabel("n. clips reverse")
     plt.colorbar(m[3], ax=ax, label="n. sites")
 
     ax = axs[1]
     f, r = df["Ff"], df["Fr"]
     bins = np.linspace(0, 1, 25)
     norm = mpl.colors.LogNorm()
     m = ax.hist2d(f, r, bins=bins, norm=norm)
-    ax.set_xlabel("freq insertions forward")
-    ax.set_ylabel("freq insertions reverse")
+    ax.set_xlabel("freq clips forward")
+    ax.set_ylabel("freq clips reverse")
     plt.colorbar(m[3], ax=ax, label="n. sites")
 
     ax = axs[2]
@@ -148,8 +148,8 @@ def plot_joint_distr(df):
     bins = np.logspace(np.log10(m) - 0.1, np.log10(M) + 0.1, 25)
     norm = mpl.colors.LogNorm()
     m = ax.hist2d(f, r, bins=bins, norm=norm)
-    ax.set_xlabel("avg. insertion length forward")
-    ax.set_ylabel("avg. insertion length reverse")
+    ax.set_xlabel("avg. clip length forward")
+    ax.set_ylabel("avg. clip length reverse")
     ax.set_xscale("log")
     ax.set_yscale("log")
     plt.colorbar(m[3], ax=ax, label="n. sites")
@@ -179,7 +179,7 @@ def plot_n_clips_genome(dfs, step=5000):
         ax.legend()
         ax.set_xlim(bins[0] - step * 5, bins[-1] + step * 5)
         ax.set_xlabel("genome position (bp)")
-        ax.set_ylabel(f"n. insertions per {step//1000} kbp")
+        ax.set_ylabel(f"n. clips per {step//1000} kbp")
         ax.grid(alpha=0.3)
         ax.set_title(f"t = {t}")
         ytl = [int(y) for y in ax.get_yticks()]
@@ -260,7 +260,7 @@ def plot_trajectories(sel_df, Ts, threshold=5, Nx=3):
 
     # get vial number
     vial = re.search("vial_(\d+)/?$", str(data_path)).groups()[0]
-    print(f"preparing insertion plots for vial {vial}")
+    print(f"preparing clip plots for vial {vial}")
 
     # load clips and consensus frequency
     clips, clip_seqs = load_clips(data_path)

scripts/plot_insertions.py

@@ -39,7 +39,7 @@ def L_tot(x):
     dfs = {}
     for t in Ts:
 
-        ins = insertions[str(t)]
+        ins = insertions[t]
         pos = np.array(sorted(ins.keys()))
         df = {}
 

scripts/plot_utils.py

@@ -49,6 +49,7 @@ def load_time_dict(data_path, filename):
     dicts = {}
     for tpf in timepoint_fld:
         time_id = m.search(str(tpf)).groups()[0]
+        time_id = int(time_id)
         tp_file = tpf / "pileup" / filename
         with gzip.open(tp_file, "r") as f:
             dicts[time_id] = pkl.load(f)
@@ -138,7 +139,9 @@ def color_dict(timepoints):
 def dict_histograms(items_dict, ax, colors, plotmeans=False, **kwargs):
     """plot the histogram of a dictionary of items"""
     means = {tp: arr.mean() for tp, arr in items_dict.items()}
-    for tp, arr in items_dict.items():
+    ordered_Ts = sorted(items_dict.keys())
+    for tp in ordered_Ts:
+        arr = items_dict[tp]
         ax.hist(arr, label=tp, histtype="step", color=colors[tp], **kwargs)
         if plotmeans:
             ax.axvline(means[tp], ls=":", color=colors[tp])