From c0a4b5cd0b8c519529f6780f88a01f42d8c2d70b Mon Sep 17 00:00:00 2001 From: AlexanderButyaev Date: Wed, 17 Sep 2025 16:51:22 -0400 Subject: [PATCH 1/3] test MPLBACKEND Agg before everything --- DashML/GUI/DT.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/DashML/GUI/DT.py b/DashML/GUI/DT.py index bacaf755..c65ed274 100644 --- a/DashML/GUI/DT.py +++ b/DashML/GUI/DT.py @@ -1,4 +1,10 @@ -import os, re +import os + +# Must be set before importing matplotlib or anything that might import pyplot +os.environ.setdefault("MPLBACKEND", "Agg") +os.environ.setdefault("QT_API", "pyqt6") + +import re import sys import numpy as np import pandas as pd From 28588dfacc7518152f2d59ba04f7cad4b74e4329 Mon Sep 17 00:00:00 2001 From: AlexanderButyaev Date: Fri, 19 Sep 2025 11:10:33 -0400 Subject: [PATCH 2/3] removed plt/global change of mpl backend + use Figure / FigureCanvasAgg in thread +local RC context (vs global) --- DashML/Basecall/Basecall_Plot.py | 425 ++++++++++++++++++++----------- 1 file changed, 281 insertions(+), 144 deletions(-) diff --git a/DashML/Basecall/Basecall_Plot.py b/DashML/Basecall/Basecall_Plot.py index 09793cec..2ef6f947 100644 --- a/DashML/Basecall/Basecall_Plot.py +++ b/DashML/Basecall/Basecall_Plot.py @@ -1,118 +1,201 @@ -import os.path +from pathlib import Path import math -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt + import numpy as np import pandas as pd -def plot_modification(mods, dir_name="Default", save_path="./", mod_type="Modifications", seq_name="Sequence"): - ### deprecated #### +import matplotlib +from matplotlib.figure import Figure +from matplotlib.backends.backend_agg import FigureCanvasAgg + +from contextlib import contextmanager + + +@contextmanager +def _mpl_context(params: dict | None = None): + with matplotlib.rc_context(rc=params or {}): + yield + + +def _ensure_dir(p: str | Path) -> Path: + p = Path(p) + p.mkdir(parents=True, exist_ok=True) + return p + + +def plot_modification( + mods, + dir_name: str = "Default", + save_path: str = "./", + mod_type: str = "Modifications", + seq_name: str = "Sequence", +): + #### deprecated #### return ##### get positions and mods ############ - modifications = np.array(mods, dtype=int) + # ----- data prep (no pyplot) ----- + modifications = np.asarray(mods, dtype=int) positions = modifications[:, 0].tolist() - mods = modifications[:, 1].tolist() - moz = [i for i, n in enumerate(mods) if n > 1] - mods = [mods[i] for i in moz] - positions = [positions[i] for i in moz] - positions_size = 50 - num_plots = math.floor(len(positions)/positions_size) + # ALBU: I assume the second column is counts + counts = modifications[:, 1].tolist() + + # keep only counts > 1 + idx = [i for i, n in enumerate(counts) if n > 1] + positions = [positions[i] for i in idx] + counts = [counts[i] for i in idx] + positions_size = 50 + num_plots = math.floor(len(positions) / positions_size) - #### set color based on mod type ####### + # ----- color selection ----- bar_color = "b" if mod_type == "Deletions": bar_color = "r" elif mod_type == "Insertions": bar_color = "g" - #### set figure length ##### - figlen = len(mods) * .3 - + # ----- figure length ----- + figlen = len(counts) * 0.3 if figlen < 6.4: figlen = 6.4 elif figlen > 100: figlen = 100 - #### plot modifications ##### - #fig, ax = plt.figure(figsize=(figlen, 4.8)) #6.4, 4.8 default size - #fig, ax = plt.subplots(num_plots) - #fig.set_size_inches(10, 50) - #fig.set_figheight(figlen) - #fig.suptitle(seq_name +' '+ dir_name +' ' + mod_type, fontsize="xx-large") + # fig.set_figwidth(10) - #fig.set_figwidth(10) - - dir_name = save_path + dir_name + '_Modification_Plots' - if not os.path.exists(dir_name): - os.makedirs(dir_name) + # ----- output dir ----- + outdir = Path(save_path) / f"{dir_name}_Modification_Plots" + outdir.mkdir(parents=True, exist_ok=True) + outputs = [] curr_pos = 0 - for i in range(num_plots): - mod_pos = np.array(positions[curr_pos:curr_pos+positions_size]) - mod_num = np.array(mods[curr_pos:curr_pos+positions_size]) - df = pd.DataFrame({'Reference Nucleotide Position': mod_pos, - 'Number of Modifications': mod_num}) - ax = df.plot.bar(x='Reference Nucleotide Position', y='Number of Modifications', rot=90, figsize=(figlen, 10), - color=bar_color, xlabel='Reference Nucleotide Position', ylabel='Number of ' + mod_type, - title=seq_name +' ' + mod_type, legend=False) - # pps = ax.bar(X, np.array(mods[curr_pos:curr_pos+positions_size]), - # color=bar_color, align="center", edgecolor="black", width=.4) - #ax.bar_label(ax, label_type='edge') - curr_pos = curr_pos + positions_size - fig = plt.gcf() - #plt.showblock=False) - figname = os.path.join(dir_name + '/' + seq_name + '_' + mod_type + '_' + str(i) + '.png') - fig.savefig(figname) - - -def plot_modification_summary(mods, dir_name="Default", save_path="./", \ - mod_type="Modification Summary", seq_name="Sequence"): + with _mpl_context({"font.size": 20}): + for i in range(num_plots): + mod_pos = np.array(positions[curr_pos : curr_pos + positions_size]) + mod_num = np.array(counts[curr_pos : curr_pos + positions_size]) + + df = pd.DataFrame( + { + "Reference Nucleotide Position": mod_pos, + "Number of Modifications": mod_num, + } + ) + + # isolated Agg figure/canvas + fig = Figure(figsize=(figlen, 10), dpi=100) + FigureCanvasAgg(fig) # attach Agg renderer + ax = fig.add_subplot(111) + + df.plot.bar( + x="Reference Nucleotide Position", + y="Number of Modifications", + rot=90, + ax=ax, + color=bar_color, + xlabel="Reference Nucleotide Position", + ylabel=f"Number of {mod_type}", + title=f"{seq_name} {mod_type}", + legend=False, + ) + + #### OLD plot modifications (just for idea - don't use) ##### + # fig, ax = plt.figure(figsize=(figlen, 4.8)) #6.4, 4.8 default size + # fig, ax = plt.subplots(num_plots) + # fig.set_size_inches(10, 50) + # fig.set_figheight(figlen) + # fig.suptitle(seq_name +' '+ dir_name +' ' + mod_type, fontsize="xx-large") + + fig.tight_layout() + + figname = outdir / f"{seq_name}_{mod_type}_{i}.png" + fig.savefig(figname, bbox_inches="tight") + outputs.append(str(figname.resolve())) + + curr_pos += positions_size + + return outputs + + +def plot_modification_summary( + mods, + dir_name: str = "Default", + save_path: str | Path = "./", + mod_type: str = "Modification Summary", + seq_name: str = "Sequence", +): ##### deprecated for now ######### return - ##### get positions and mods ############ - modifications = np.array(mods, dtype=int) + # ----- data prep ----- + modifications = np.asarray(mods, dtype=int) positions = modifications[:, 0].tolist() - dels = (modifications[:, 1]/np.linalg.norm(modifications[:, 1])).tolist() - ins = (modifications[:, 2]/np.linalg.norm(modifications[:,2])).tolist() - mismatch = (modifications[:, 3]/np.linalg.norm(modifications[:,3])).tolist() + + # Safe norm (avoid divide-by-zero) + def _safe_norm(v: np.ndarray) -> np.ndarray: + n = np.linalg.norm(v) + return (v / n) if n > 0 else np.zeros_like(v, dtype=float) + + dels = _safe_norm(modifications[:, 1]).tolist() + ins = _safe_norm(modifications[:, 2]).tolist() + mismatch = _safe_norm(modifications[:, 3]).tolist() + positions_size = 50 num_plots = math.floor(len(positions) / positions_size) - #### set figure length ##### - figlen = len(mods) * .3 + # ----- figure length ----- + figlen = len(mods) * 0.3 if figlen < 6.4: figlen = 6.4 elif figlen > 6.4: figlen = 30 - dir_name = save_path + dir_name + '_Modification_Plots' - if not os.path.exists(dir_name): - os.makedirs(dir_name) + # ----- output dir ----- + outdir = Path(save_path) / f"{dir_name}_Modification_Plots" + outdir.mkdir(parents=True, exist_ok=True) + outputs: list[str] = [] curr_pos = 0 - for i in range(num_plots): - mod_pos = np.array(positions[curr_pos:curr_pos + positions_size]) - mod_del = np.array(dels[curr_pos:curr_pos + positions_size]) - mod_ins = np.array(ins[curr_pos:curr_pos + positions_size]) - mod_mis = np.array(mismatch[curr_pos:curr_pos + positions_size]) - - df = pd.DataFrame({'Reference Nucleotide Position': mod_pos, - 'Deletions': mod_del, 'Insertions': mod_ins, - 'Mismatches': mod_mis}) - ax = df.plot.bar(x='Reference Nucleotide Position', rot=90, figsize=(figlen, 10), - xlabel='Reference Nucleotide Position', ylabel='Number of ' + mod_type, - title=seq_name + ' ' + mod_type, legend=True, ylim=(0,.01)) - # pps = ax.bar(X, np.array(mods[curr_pos:curr_pos+positions_size]), - # color=bar_color, align="center", edgecolor="black", width=.4) - # ax.bar_label(ax, label_type='edge') - curr_pos = curr_pos + positions_size - fig = plt.gcf() - #plt.showblock=False) - figname = os.path.join(dir_name + '/' + seq_name + '_' + mod_type + '_' + str(i) + '.png') - fig.savefig(figname) + with _mpl_context({"font.size": 20}): + for i in range(num_plots): + mod_pos = np.array(positions[curr_pos : curr_pos + positions_size]) + mod_del = np.array(dels[curr_pos : curr_pos + positions_size]) + mod_ins = np.array(ins[curr_pos : curr_pos + positions_size]) + mod_mis = np.array(mismatch[curr_pos : curr_pos + positions_size]) + + df = pd.DataFrame( + { + "Reference Nucleotide Position": mod_pos, + "Deletions": mod_del, + "Insertions": mod_ins, + "Mismatches": mod_mis, + } + ) + + # isolated Agg figure/canvas + fig = Figure(figsize=(figlen, 10), dpi=100) + FigureCanvasAgg(fig) # attach Agg renderer + ax = fig.add_subplot(111) + + df.plot.bar( + x="Reference Nucleotide Position", + rot=90, + ax=ax, + xlabel="Reference Nucleotide Position", + ylabel=f"Number of {mod_type}", + title=f"{seq_name} {mod_type}", + legend=True, + ylim=(0, 0.01), + ) + + fig.tight_layout() + + figname = outdir / f"{seq_name}_{mod_type}_{i}.png" + fig.savefig(figname, bbox_inches="tight") + outputs.append(str(figname.resolve())) + + curr_pos += positions_size + + return outputs # #### plot modifications ##### # fig = plt.figure(figsize=(figlen, 4.8)) #6.4, 4.8 default size @@ -141,86 +224,140 @@ def plot_modification_summary(mods, dir_name="Default", save_path="./", \ # plt.tight_layout() # #plt.show) -def plot_mismatch(mismatches, seq_name="Sequence", dir_name="Default", save_path="./", mod_type="Mismatched Bases"): + +def plot_mismatch( + mismatches, + seq_name: str = "Sequence", + dir_name: str = "Default", + save_path: str = "./", + mod_type: str = "Mismatched Bases", +): #### deprecated for now #### return ##### get positions and mods ############ - data = np.array(mismatches) - positions = data[:,0].tolist() - num_mismatch = data[:,1].tolist() - nuc_mismatch = data[:,2].tolist() - ref_seq_label = data[:,3].tolist() + # ----- data prep (no pyplot) ----- + data = np.asarray(mismatches) + positions = data[:, 0].tolist() + num_mismatch = data[:, 1].astype(float).tolist() + nuc_mismatch = data[:, 2].tolist() + ref_seq_label = data[:, 3].tolist() - ref_label = [] - for i in range(len(ref_seq_label)): - ref_label.append(str(positions[i]) +" "+ str(ref_seq_label[i])) + ref_label = [ + f"{positions[i]} {ref_seq_label[i]}" for i in range(len(ref_seq_label)) + ] print(seq_name) #### set figure length ##### - figlen = len(positions) * .3 + figlen = len(positions) * 0.3 if figlen < 6.4: figlen = 6.4 elif figlen > 100: figlen = 100 - #### plot modifications ##### - fig = plt.figure(figsize=(figlen, 6)) # 6.4, 4.8 default size - ax = fig.add_axes([.1, .2, .8, .7]) - ax.margins(.005, .2) - ax.set_xticks(range(0, len(positions))) - ax.set_xticklabels(ref_label) - ax.tick_params(direction='out', length=6, rotation=90) - ax.set_xlabel('Reference Nucleotide Position', fontsize="xx-large") - ax.set_ylabel('Number of ' + mod_type, fontsize="xx-large") - ax.set_title(seq_name + ' ' + mod_type, fontsize="xx-large") - ax.grid(axis='y') - # proper ticks - X = np.arange(len(num_mismatch)) - pps = ax.bar(X, np.array(num_mismatch).astype(float), color="b", align="center", edgecolor="black", width=.2) - ax.bar_label(pps, nuc_mismatch, label_type='edge') - #### save plot ##### - fig = plt.gcf() - dir_name = save_path + dir_name + '_Modification_Plots' - if not os.path.exists(dir_name): - os.makedirs(dir_name) - figname = os.path.join(dir_name + '/' + seq_name + '_Mismatched_Bases' + '.png') - fig.savefig(figname) - #plt.showblock=False) + # ----- output dir ----- + outdir = Path(save_path) / f"{dir_name}_Modification_Plots" + outdir.mkdir(parents=True, exist_ok=True) + + # ----- render without pyplot/global backend ----- + with matplotlib.rc_context({"font.size": 20}): + fig = Figure(figsize=(figlen, 6), dpi=100) + FigureCanvasAgg(fig) # attach Agg renderer + ax = fig.add_axes([0.1, 0.2, 0.8, 0.7]) + + ax.margins(0.005, 0.2) + ax.set_xticks(range(len(positions))) + ax.set_xticklabels(ref_label) + ax.tick_params(direction="out", length=6, rotation=90) + ax.set_xlabel("Reference Nucleotide Position", fontsize="xx-large") + ax.set_ylabel(f"Number of {mod_type}", fontsize="xx-large") + ax.set_title(f"{seq_name} {mod_type}", fontsize="xx-large") + ax.grid(axis="y") + + X = np.arange(len(num_mismatch)) + bars = ax.bar( + X, + np.asarray(num_mismatch, dtype=float), + color="b", + align="center", + edgecolor="black", + width=0.2, + ) + + # Ensure labels length matches the bars length + if len(nuc_mismatch) >= len(bars): + labels = nuc_mismatch[: len(bars)] + else: + labels = nuc_mismatch + [""] * (len(bars) - len(nuc_mismatch)) + + ax.bar_label(bars, labels, label_type="edge") + + figname = outdir / f"{seq_name}_Mismatched_Bases.png" + fig.savefig(figname, bbox_inches="tight") + + return str(figname.resolve()) + def plot_average_mod_rate(df, dir_name="Default", save_path="./"): - plt.rcParams.update({'font.size': 20}) - df = df.drop("Condition", axis=1) - #x = df.values # returns a numpy array - #min_max_scaler = preprocessing.MinMaxScaler() - #x_scaled = min_max_scaler.fit_transform(x) - #df = pd.DataFrame(x_scaled) - df = df.iloc[:,:] * 100 - fig = df.plot(kind="bar", rot=30, title="Overall Modification Rates", figsize=(12, 12)) - #### save plot ##### - fig = plt.gcf() - dir_name = save_path + dir_name + '_Modification_Plots' - if not os.path.exists(dir_name): - os.makedirs(dir_name) - figname = os.path.join(dir_name + '/Average_Modification_Rate' + '.png') - fig.savefig(figname) - #plt.showblock=False) + print("plot_average_mod_rate") + + # --- Data prep --- + df2 = df.drop(columns=["Condition"], errors="ignore") + df2 = df2 * 100.0 # percent + + # --- Build figure without pyplot/global backend --- + with _mpl_context({"font.size": 20}): + fig = Figure(figsize=(12, 12), dpi=100) # independent of GUI backend + FigureCanvasAgg(fig) # attach Agg renderer to this Figure + ax = fig.add_subplot(111) + + # Use pandas plotting but direct it to our Axes (no pyplot) + df2.plot(kind="bar", rot=30, title="Overall Modification Rates", ax=ax) + + fig.tight_layout() + + outdir = _ensure_dir(Path(save_path) / f"{dir_name}_Modification_Plots") + figname = outdir / "Average_Modification_Rate.png" + + # Save via the Figure object (thread-safe) + fig.savefig(figname, bbox_inches="tight") + + # Return absolute path for the GUI thread to load + return str(figname.resolve()) + def plot_average_mod_by_pos_rate(df, dir_name="Default", save_path="./"): - plt.rcParams.update({'font.size': 20}) - #x = df.values # returns a numpy array - #min_max_scaler = preprocessing.MinMaxScaler() - #x_scaled = min_max_scaler.fit_transform(x) - #df = pd.DataFrame(x_scaled) - df = df.iloc[:, :] * 100 - fig = df.plot(kind="line", rot=45, title="Modification Rates by Position", subplots=True, figsize=(12,8), - ylim=(0,100), color='purple') - #fig = df.plot(kind="bar", rot=90, title="Modification Rates by Position", figsize=(12,8), stacked=True) - #### save plot ##### - fig = plt.gcf() - dir_name = save_path + dir_name + '_Modification_Plots' - if not os.path.exists(dir_name): - os.makedirs(dir_name) - figname = os.path.join(dir_name + '/Position_Modification_Rate' + '.png') - fig.savefig(figname) - #plt.showblock=False) + print("plot_average_mod_by_pos_rate") + + # --- Data prep (no plotting side-effects) --- + df2 = df * 100.0 # percent + + # x = df.values # returns a numpy array + # min_max_scaler = preprocessing.MinMaxScaler() + # x_scaled = min_max_scaler.fit_transform(x) + # df = pd.DataFrame(x_scaled) + # df = df.iloc[:, :] * 100 + # --- Build figure without pyplot/global backend --- + with _mpl_context({"font.size": 20}): + fig = Figure(figsize=(12, 8), dpi=100) # independent of GUI backend + FigureCanvasAgg(fig) # attach Agg renderer to this Figure + ax = fig.add_subplot(111) + df2.plot( + kind="line", + rot=45, + title="Modification Rates by Position", + ylim=(0, 100), + color="purple", + ax=ax, + ) + + fig.tight_layout() + + outdir = _ensure_dir(Path(save_path) / f"{dir_name}_Modification_Plots") + figname = outdir / "Position_Modification_Rate.png" + + # Save via the Figure object (thread-safe) + fig.savefig(figname, bbox_inches="tight") + + # Return absolute path for the GUI thread to load + return str(figname.resolve()) From 7d3aaf4b683f15941ea2a0400e82d11610b31cd0 Mon Sep 17 00:00:00 2001 From: AlexanderButyaev Date: Mon, 22 Sep 2025 13:59:05 -0400 Subject: [PATCH 3/3] reworked almost all plot related stuff --- DashML/Basecall/Basecall_Bias.py | 276 ++-- DashML/Basecall/Basecall_Plot.py | 11 +- DashML/GUI/DT.py | 1127 +++++++++++++---- DashML/Landscape/Cluster/Centroid_Analysis.py | 87 +- DashML/Landscape/Cluster/Centroid_Fold.py | 342 +++-- DashML/Landscape/Cluster/Centroid_MFE.py | 228 ++-- DashML/Landscape/Cluster/Cluster_Num.py | 297 +++-- DashML/Landscape/Cluster/Cluster_means.py | 282 +++-- DashML/Landscape/Cluster/Cluster_mode.py | 262 ++-- DashML/Landscape/Cluster/Cluster_native.py | 218 ++-- DashML/Predict/Gmm_Analysis.py | 109 +- DashML/Predict/Gmm_Analysis2.py | 129 +- DashML/Predict/Lof.py | 69 +- DashML/Predict/Peak_Analysis_BC.py | 107 +- DashML/Predict/Predict_BPP.py | 213 ++-- DashML/Predict/Predicts.py | 175 ++- DashML/Predict/run_predict.py | 268 ++-- 17 files changed, 2768 insertions(+), 1432 deletions(-) diff --git a/DashML/Basecall/Basecall_Bias.py b/DashML/Basecall/Basecall_Bias.py index 9f50d684..78a6bee8 100644 --- a/DashML/Basecall/Basecall_Bias.py +++ b/DashML/Basecall/Basecall_Bias.py @@ -1,30 +1,39 @@ import os.path import re +from matplotlib.figure import Figure +from matplotlib.backends.backend_agg import FigureCanvasAgg import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt import numpy as np import pandas as pd + class Modification_Bias: - def __init__(self,reference_sequences, f_path, structure="", loop_structure="", dir_name="Default"): + def __init__( + self, + reference_sequences, + f_path, + structure="", + loop_structure="", + dir_name="Default", + ): self.dir_name = dir_name self.reference_sequences = reference_sequences self.f_path = f_path - #get base pairing data + # get base pairing data if structure: self.structure = structure - #get loop structure data - if 'single_loop_structure' in loop_structure: - self.single_loop_structure = loop_structure.get('single_loop_structure') - if 'paired_loop_structure' in loop_structure: - self.paired_loop_structure = loop_structure.get('paired_loop_structure') + # get loop structure data + if "single_loop_structure" in loop_structure: + self.single_loop_structure = loop_structure.get("single_loop_structure") + if "paired_loop_structure" in loop_structure: + self.paired_loop_structure = loop_structure.get("paired_loop_structure") ###### Create Output Folder ifne ########## - self.save_path = self.f_path + 'Out_Structure/' + self.dir_name + '_Modification_Rates' + self.save_path = ( + self.f_path + "Out_Structure/" + self.dir_name + "_Modification_Rates" + ) print("Modification Output To: " + self.save_path) - ########### Plot Modifications ################# def parse_modifications(self, position, indel): df = pd.DataFrame() @@ -46,26 +55,44 @@ def parse_modifications(self, position, indel): if indel[key][i][2] > 0: mismatch.append((position[key][i], indel[key][i][2])) if len(indel[key][i][3]) > 0: - #pos, number of mismatches,mismatch nucleotide,refseq nucleotide - mismatches.append((position[key][i], indel[key][i][3][0][1], indel[key][i][3][0][0],sequence[i])) - if (indel[key][i][0] > 0 or indel[key][i][1] > 0 or indel[key][i][2] > 0): - summary.append((position[key][i], indel[key][i][0], indel[key][i][1], indel[key][i][2])) + # pos, number of mismatches,mismatch nucleotide,refseq nucleotide + mismatches.append( + ( + position[key][i], + indel[key][i][3][0][1], + indel[key][i][3][0][0], + sequence[i], + ) + ) + if ( + indel[key][i][0] > 0 + or indel[key][i][1] > 0 + or indel[key][i][2] > 0 + ): + summary.append( + ( + position[key][i], + indel[key][i][0], + indel[key][i][1], + indel[key][i][2], + ) + ) ### if modifications then plot ######## - #if len(deletions) > 0: - #pltmod.plot_modification(deletions, seq_name=key, mod_type="Deletions", dir_name=dir_name) - #if len(insertions) > 0: - #pltmod.plot_modification(insertions, seq_name=key, mod_type="Insertions", dir_name=dir_name) - #if len(mismatch) > 0: - #pltmod.plot_modification(mismatch, seq_name=key, mod_type="Mismatch", dir_name=dir_name) + # if len(deletions) > 0: + # pltmod.plot_modification(deletions, seq_name=key, mod_type="Deletions", dir_name=dir_name) + # if len(insertions) > 0: + # pltmod.plot_modification(insertions, seq_name=key, mod_type="Insertions", dir_name=dir_name) + # if len(mismatch) > 0: + # pltmod.plot_modification(mismatch, seq_name=key, mod_type="Mismatch", dir_name=dir_name) if len(mismatches) > 0: - #### calculate mismatch bias #### - pos_data, df_tmp = self.mismatch_bias(mismatches, key) - df = df.append(df_tmp) - pymol[key] = self.pymol_parse(data=pos_data) - #pltmod.plot_mismatch(mismatches, seq_name=key, dir_name=dir_name) + #### calculate mismatch bias #### + pos_data, df_tmp = self.mismatch_bias(mismatches, key) + df = df.append(df_tmp) + pymol[key] = self.pymol_parse(data=pos_data) + # pltmod.plot_mismatch(mismatches, seq_name=key, dir_name=dir_name) ##### Plot Summary of Modifications ############ - #if len(summary) > 0: - #pltmod.plot_modification_summary(summary, dir_name=dir_name, seq_name=key) + # if len(summary) > 0: + # pltmod.plot_modification_summary(summary, dir_name=dir_name, seq_name=key) deletions = [] insertions = [] mismatch = [] @@ -79,14 +106,31 @@ def parse_modifications(self, position, indel): ####### Modifications by position ############ df["Position"] = pd.to_numeric(df["Position"]) df = df.reset_index(drop=True) - #print(df.dtypes) #print(df.columns) - df = df.loc[:, ["Sequence", "Position", "Mismatches","Number of Mismatches","Total Mismatches in Aligned Sequences", "Modified Bp Freq"]] - df.sort_values(by=['Sequence', 'Position'], ascending=[True, True], inplace=True) + # print(df.dtypes) #print(df.columns) + df = df.loc[ + :, + [ + "Sequence", + "Position", + "Mismatches", + "Number of Mismatches", + "Total Mismatches in Aligned Sequences", + "Modified Bp Freq", + ], + ] + df.sort_values( + by=["Sequence", "Position"], ascending=[True, True], inplace=True + ) self.print_df(df, f_suffix="_combined") ####### Modifications in sequence grouped by mismatch ########## df = df.groupby("Mismatches", as_index=False)["Number of Mismatches"].sum() - df.insert(2, "Modified Bp Freq", (df.iloc[:,1]/df.iloc[:,1].sum()) * 100, allow_duplicates=True) - #print(df) + df.insert( + 2, + "Modified Bp Freq", + (df.iloc[:, 1] / df.iloc[:, 1].sum()) * 100, + allow_duplicates=True, + ) + # print(df) ######## parse pymol commmands for modified positions ########## def pymol_parse(self, data): @@ -106,24 +150,39 @@ def print_pymol_modifications(self, pymol): pymol_df.columns = ["Sequence", "Mismatch Positions"] self.print_df(pymol_df, f_suffix="_pymol") - #### calculate mismatch bias per sequence #### def mismatch_bias(self, mismatches, curr_seq): - base = re.compile('[A-Z]', re.IGNORECASE) + base = re.compile("[A-Z]", re.IGNORECASE) bias = {} pymol = [] data = np.array(mismatches) df = pd.DataFrame(mismatches) - df.columns = ["Position", "Number of Mismatches", "Called Base", "Reference Base"] + df.columns = [ + "Position", + "Number of Mismatches", + "Called Base", + "Reference Base", + ] df["Position"] = pd.to_numeric(df["Position"]) - df['Called Base'] = df['Called Base'].str.replace('[^A-Z]', '', regex=True) - df['Reference Base'] = df['Reference Base'].str.replace('[^A-Z]', '', regex=True) - df['Mismatches'] = df['Reference Base'] + df['Called Base'] - df = df.drop(columns=['Called Base', 'Reference Base']) - df.insert(3, "Unmodified Structure", df['Position'].apply(lambda x: 'S' if x in self.structure['single_strand_positions'] else 'BP'), allow_duplicates=True) + df["Called Base"] = df["Called Base"].str.replace("[^A-Z]", "", regex=True) + df["Reference Base"] = df["Reference Base"].str.replace( + "[^A-Z]", "", regex=True + ) + df["Mismatches"] = df["Reference Base"] + df["Called Base"] + df = df.drop(columns=["Called Base", "Reference Base"]) + df.insert( + 3, + "Unmodified Structure", + df["Position"].apply( + lambda x: ( + "S" if x in self.structure["single_strand_positions"] else "BP" + ) + ), + allow_duplicates=True, + ) - #print(df) + # print(df) # if len(self.single_loop_structure) > 1: # df_single_loop_structure = pd.DataFrame(self.single_loop_structure, columns=["Position", "Reference Base", "Loop Structure"]) # df_single_loop_structure = df_single_loop_structure.drop(["Reference Base"], axis =1) @@ -131,7 +190,9 @@ def mismatch_bias(self, mismatches, curr_seq): # df = df.fillna(0) # #print(df_single_loop_structure) # #print(df) - df.sort_values(by=['Position', 'Mismatches', 'Unmodified Structure'], inplace=True) + df.sort_values( + by=["Position", "Mismatches", "Unmodified Structure"], inplace=True + ) ######### Print Tables TODO: why??? ############ self.print_df(df) ########## Plot Modification Data ############# @@ -142,24 +203,25 @@ def mismatch_bias(self, mismatches, curr_seq): #### total number of mismatches in aligned sequences * number of aligned reads # TODO: may want to include number of aligned reads tot = np.sum(data[:, 1].astype(np.int)) - #print("Tot: " + str(tot)) - #print(data[:,2]) + # print("Tot: " + str(tot)) + # print(data[:,2]) ######## compile string for pymol command ##### self.pymol_parse(data[:, 0]) - #TODO: add position for strand to loop data[:, 0] + # TODO: add position for strand to loop data[:, 0] for i in range(len(data[:, 3])): - #remove artifacts from pysam indel notation + # remove artifacts from pysam indel notation ref_base = re.search(base, np.array2string(data[i, 3])).group() called_base = re.search(base, np.array2string(data[i, 2])).group() if ref_base != called_base: key = ref_base + called_base - value = data[i,1].astype(np.int) + value = data[i, 1].astype(np.int) if bias.get(key): curr_val = bias.get(key)[0] - else: curr_val = 0 - pos = data[i,0] + else: + curr_val = 0 + pos = data[i, 0] if curr_val == None: bias[key] = (value, pos) else: @@ -168,52 +230,108 @@ def mismatch_bias(self, mismatches, curr_seq): ######## Prepare Data Frame ############## df = pd.DataFrame.from_dict(bias, orient="index") df = df.reset_index() - df.columns = ["Mismatches", "Number of Mismatches","Position"] - df.insert(0, "Sequence", [curr_seq]*len(df), allow_duplicates=True) - #print(df) - df.insert(3, "Total Mismatches in Aligned Sequences", df[ "Number of Mismatches"].sum(), allow_duplicates=True) + df.columns = ["Mismatches", "Number of Mismatches", "Position"] + df.insert(0, "Sequence", [curr_seq] * len(df), allow_duplicates=True) + # print(df) + df.insert( + 3, + "Total Mismatches in Aligned Sequences", + df["Number of Mismatches"].sum(), + allow_duplicates=True, + ) ### TODO: total mismatches at position percentage would be better - df.insert(4, "Modified Bp Freq", df["Number of Mismatches"]/df["Total Mismatches in Aligned Sequences"] * 100, allow_duplicates=True) + df.insert( + 4, + "Modified Bp Freq", + df["Number of Mismatches"] + / df["Total Mismatches in Aligned Sequences"] + * 100, + allow_duplicates=True, + ) df = df.reset_index(drop=True) self.print_df(df, curr_seq, f_suffix="_modification_bias") return data[:, 0], df # TODO: add whether a modification occurs on single strand or on known base pairs, then later in/out helix - #TODO: add function calculates bias vs guassian, vs unmodified tetrahymena - #TODO: finally bias calculation for all modifications + # TODO: add function calculates bias vs guassian, vs unmodified tetrahymena + # TODO: finally bias calculation for all modifications ######## Print Data Frame to CSV ############## - def print_df(self,df, curr_seq="", f_suffix=""): + def print_df(self, df, curr_seq="", f_suffix=""): if not os.path.exists(self.save_path): os.makedirs(self.save_path) ###### output modification bias for file ########### if curr_seq: - df.to_csv(self.save_path + "/" + curr_seq + "_" + self.dir_name + f_suffix + ".csv", index=False) + df.to_csv( + self.save_path + + "/" + + curr_seq + + "_" + + self.dir_name + + f_suffix + + ".csv", + index=False, + ) else: - df.to_csv(self.save_path + "/" + self.dir_name + f_suffix + ".csv", index=False) + df.to_csv( + self.save_path + "/" + self.dir_name + f_suffix + ".csv", index=False + ) - def plot_df(self, df, curr_seq): - df['Number of Mismatches'] = (df['Number of Mismatches']/df['Number of Mismatches'].sum() ) * 100 - df = df.loc[df['Number of Mismatches'] > df['Number of Mismatches'].mean()] - df['Mismatches'] = df['Mismatches'] + df['Unmodified Structure'] - #print(df) + def plot_df(self, df: pd.DataFrame, curr_seq: str): + df = df.copy() + total = float(df["Number of Mismatches"].sum()) + if total > 0: + df["Number of Mismatches"] = (df["Number of Mismatches"] / total) * 100 + else: + df["Number of Mismatches"] = 0.0 # Avoid division by zero - fig = df.plot(x='Position', y = 'Number of Mismatches', kind="bar", rot=45, title= curr_seq + " Mismatches by Structure", figsize=(12,8)) - # fig = df.plot(kind="bar", rot=90, title="Modification Rates by Position", figsize=(12,8), stacked=True) - #### save plot ##### - i = 0 - score = df['Mismatches'].to_numpy() - for p in fig.patches: - fig.annotate(score[i], xy=(p.get_x(), p.get_height())) - i += 1 + df = df.loc[df["Number of Mismatches"] > df["Number of Mismatches"].mean()] + # ensure numeric then combine labels (or keep as strings) + df["Mismatches"] = df["Mismatches"].astype(str) + df[ + "Unmodified Structure" + ].astype(str) - fig = plt.gcf() + outdir = self.save_path + os.makedirs(outdir, exist_ok=True) - if not os.path.exists(self.save_path): - os.makedirs(self.save_path) + with matplotlib.rc_context({"font.size": 12}): + fig = Figure(figsize=(12, 8), dpi=100) + FigureCanvasAgg(fig) # attach Agg renderer + ax = fig.add_subplot(111) + + df.plot( + x="Position", + y="Number of Mismatches", + kind="bar", + rot=45, + title=f"{curr_seq} Mismatches by Structure", + ax=ax, + legend=False, + ) + + labels = df["Mismatches"].to_numpy(dtype=str) + bars = [p for p in ax.patches] + for i, p in enumerate(bars[: len(labels)]): + # place label slightly above the bar + ax.annotate( + labels[i], + xy=(p.get_x() + p.get_width() / 2.0, p.get_height()), + xytext=(0, 3), + textcoords="offset points", + ha="center", + va="bottom", + fontsize=8, + ) + + fig.tight_layout() + + figname = os.path.join( + outdir, f"{curr_seq}_Position_Structure_Modification_Rate.png" + ) + fig.savefig(figname, bbox_inches="tight") + fig.clear() + del fig - figname = os.path.join(self.save_path + '/'+ curr_seq + '_' + 'Position_Structure_Modification_Rate' + '.png') - fig.savefig(figname) - #plt.show) + return figname diff --git a/DashML/Basecall/Basecall_Plot.py b/DashML/Basecall/Basecall_Plot.py index 2ef6f947..bb84ff86 100644 --- a/DashML/Basecall/Basecall_Plot.py +++ b/DashML/Basecall/Basecall_Plot.py @@ -33,7 +33,7 @@ def plot_modification( #### deprecated #### return ##### get positions and mods ############ - # ----- data prep (no pyplot) ----- + # ----- data prep ----- modifications = np.asarray(mods, dtype=int) positions = modifications[:, 0].tolist() # ALBU: I assume the second column is counts @@ -235,7 +235,7 @@ def plot_mismatch( #### deprecated for now #### return ##### get positions and mods ############ - # ----- data prep (no pyplot) ----- + # ----- data prep ----- # data = np.asarray(mismatches) positions = data[:, 0].tolist() num_mismatch = data[:, 1].astype(float).tolist() @@ -259,7 +259,7 @@ def plot_mismatch( outdir = Path(save_path) / f"{dir_name}_Modification_Plots" outdir.mkdir(parents=True, exist_ok=True) - # ----- render without pyplot/global backend ----- + # ----- render without global backend ----- with matplotlib.rc_context({"font.size": 20}): fig = Figure(figsize=(figlen, 6), dpi=100) FigureCanvasAgg(fig) # attach Agg renderer @@ -305,13 +305,12 @@ def plot_average_mod_rate(df, dir_name="Default", save_path="./"): df2 = df.drop(columns=["Condition"], errors="ignore") df2 = df2 * 100.0 # percent - # --- Build figure without pyplot/global backend --- + # --- Build figure without global backend --- with _mpl_context({"font.size": 20}): fig = Figure(figsize=(12, 12), dpi=100) # independent of GUI backend FigureCanvasAgg(fig) # attach Agg renderer to this Figure ax = fig.add_subplot(111) - # Use pandas plotting but direct it to our Axes (no pyplot) df2.plot(kind="bar", rot=30, title="Overall Modification Rates", ax=ax) fig.tight_layout() @@ -337,7 +336,7 @@ def plot_average_mod_by_pos_rate(df, dir_name="Default", save_path="./"): # x_scaled = min_max_scaler.fit_transform(x) # df = pd.DataFrame(x_scaled) # df = df.iloc[:, :] * 100 - # --- Build figure without pyplot/global backend --- + # --- Build figure without global backend --- with _mpl_context({"font.size": 20}): fig = Figure(figsize=(12, 8), dpi=100) # independent of GUI backend FigureCanvasAgg(fig) # attach Agg renderer to this Figure diff --git a/DashML/GUI/DT.py b/DashML/GUI/DT.py index c65ed274..546beca7 100644 --- a/DashML/GUI/DT.py +++ b/DashML/GUI/DT.py @@ -1,8 +1,6 @@ import os -# Must be set before importing matplotlib or anything that might import pyplot -os.environ.setdefault("MPLBACKEND", "Agg") -os.environ.setdefault("QT_API", "pyqt6") +os.environ.setdefault("QT_API", "PyQt6") import re import sys @@ -17,18 +15,28 @@ import DashML.Database_fx.Insert_DB as dbins import DashML.Database_fx.Select_DB as dbsel from PyQt6.QtWidgets import ( - QApplication, QWidget, QFrame, QVBoxLayout, QHBoxLayout, - QLabel, QPushButton, QCheckBox, QLineEdit, QFileDialog, QSizePolicy, - QScrollArea, QMessageBox, QComboBox + QApplication, + QWidget, + QFrame, + QVBoxLayout, + QHBoxLayout, + QLabel, + QPushButton, + QCheckBox, + QLineEdit, + QFileDialog, + QSizePolicy, + QScrollArea, + QMessageBox, + QComboBox, ) -from PyQt6.QtCore import Qt,QTimer, QObject, pyqtSignal, QThread -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt +from PyQt6.QtCore import Qt, QTimer, QObject, pyqtSignal, QThread + import matplotlib.image as mpimg from matplotlib.backends.backend_qtagg import FigureCanvasQTAgg as FigureCanvas from matplotlib.figure import Figure + class SelectSampleSection(QFrame): sample_changed = pyqtSignal() @@ -57,11 +65,19 @@ def init_ui(self): self.id_input.setReadOnly(True) self.id_input.setStyleSheet("background-color: #D3D3D3;") self.contig_combo = QComboBox() - self.df['combo'] = self.df["ID"].astype(str) + " " + self.df["contig"] + " " + self.df["type1"] + " " + self.df["type2"] - self.contig_combo.addItems(self.df['combo']) + self.df["combo"] = ( + self.df["ID"].astype(str) + + " " + + self.df["contig"] + + " " + + self.df["type1"] + + " " + + self.df["type2"] + ) + self.contig_combo.addItems(self.df["combo"]) # Select first item properly if not self.df.empty: - self.contig_combo.setCurrentText(self.df['combo'].iloc[0]) + self.contig_combo.setCurrentText(self.df["combo"].iloc[0]) self.lid = self.df["ID"].iloc[0] self.contig = self.df["contig"].iloc[0] @@ -73,43 +89,73 @@ def init_ui(self): id_contig_layout.addWidget(QLabel("ID:")) id_contig_layout.addWidget(self.id_input) id_contig_layout.addSpacing(20) - id_contig_layout.addWidget(QLabel("Contig (Select Existing Contig to Load Data):")) + id_contig_layout.addWidget( + QLabel("Contig (Select Existing Contig to Load Data):") + ) id_contig_layout.addWidget(self.contig_combo) id_contig_layout.addWidget(self.contig_input) layout.addLayout(id_contig_layout) # Sequence name comes first layout.addWidget(QLabel("Sequence Name (Reference Name in fasta file)")) - self.sequence_name_input = QLineEdit(self.df["sequence_name"].iloc[0] if not self.df.empty else "") + self.sequence_name_input = QLineEdit( + self.df["sequence_name"].iloc[0] if not self.df.empty else "" + ) layout.addWidget(self.sequence_name_input) # Then nucleotide sequence layout.addWidget(QLabel("Sequence (Nucleotide Sequence)")) - self.sequence_input = QLineEdit(self.df["sequence"].iloc[0] if not self.df.empty else "") + self.sequence_input = QLineEdit( + self.df["sequence"].iloc[0] if not self.df.empty else "" + ) self.sequence_input.setMinimumWidth(400) layout.addWidget(self.sequence_input) # structure and secondary - layout.addWidget(QLabel("Secondary Structure (Optional Confirmed Secondary Structure w/o Psuedoknots)")) - self.secondary_input = QLineEdit(self.df["secondary"].iloc[0] if not self.df.empty else "") + layout.addWidget( + QLabel( + "Secondary Structure (Optional Confirmed Secondary Structure w/o Psuedoknots)" + ) + ) + self.secondary_input = QLineEdit( + self.df["secondary"].iloc[0] if not self.df.empty else "" + ) self.secondary_input.setMinimumWidth(400) layout.addWidget(self.secondary_input) - layout.addWidget(QLabel("Experiment (Optional Experimental Conf. eg Xray Crystallography)")) - self.experiment_input = QLineEdit(self.df["experiment"].iloc[0] if not self.df.empty else "") + layout.addWidget( + QLabel("Experiment (Optional Experimental Conf. eg Xray Crystallography)") + ) + self.experiment_input = QLineEdit( + self.df["experiment"].iloc[0] if not self.df.empty else "" + ) self.experiment_input.setMinimumWidth(400) layout.addWidget(self.experiment_input) - self.temperature_input = QLineEdit(str(self.df["temp"].iloc[0]) if not self.df.empty else "") - self.condition_1_input = QLineEdit(self.df["type1"].iloc[0] if not self.df.empty else "") - self.condition_2_input = QLineEdit(self.df["type2"].iloc[0] if not self.df.empty else "") - self.experimental_run_input = QLineEdit(str(self.df["run"].iloc[0]) if not self.df.empty else "") + self.temperature_input = QLineEdit( + str(self.df["temp"].iloc[0]) if not self.df.empty else "" + ) + self.condition_1_input = QLineEdit( + self.df["type1"].iloc[0] if not self.df.empty else "" + ) + self.condition_2_input = QLineEdit( + self.df["type2"].iloc[0] if not self.df.empty else "" + ) + self.experimental_run_input = QLineEdit( + str(self.df["run"].iloc[0]) if not self.df.empty else "" + ) self.fields = [ ("Temperature (Default temperature is 37)", self.temperature_input), ("Condition 1 (Unmodified eg DMSO)", self.condition_1_input), - ("Condition 2 (Probe Reagent eg 1M7, repeat Condition 1 if unmodified)", self.condition_2_input), - ("Experimental Run (Multiple experiments on same sequence should have separate run numbers.)", self.experimental_run_input) + ( + "Condition 2 (Probe Reagent eg 1M7, repeat Condition 1 if unmodified)", + self.condition_2_input, + ), + ( + "Experimental Run (Multiple experiments on same sequence should have separate run numbers.)", + self.experimental_run_input, + ), ] for label, widget in self.fields: layout.addWidget(QLabel(label)) @@ -160,19 +206,27 @@ def validate_inputs_and_collect_errors(self): secondary = self.secondary_input.text().strip() experiment = self.experiment_input.text().strip() - if secondary != '': + if secondary != "": if not re.fullmatch(r"^[().]+$", secondary): - errors.append("Secondary structure must be in valid dot-bracket notation (only '.', '(', ')').") + errors.append( + "Secondary structure must be in valid dot-bracket notation (only '.', '(', ')')." + ) if not experiment: - errors.append("If specifying a control structure, an experiment must also be provided.") + errors.append( + "If specifying a control structure, an experiment must also be provided." + ) if len(self.sequence_input.text().strip()) != len(secondary): - errors.append("Sequence and secondary structure must be the same length.") + errors.append( + "Sequence and secondary structure must be the same length." + ) # Check if experiment is given without secondary (if needed) if experiment and not secondary: - errors.append("If specifying an experiment, a secondary control structure must also be provided.") + errors.append( + "If specifying an experiment, a secondary control structure must also be provided." + ) return errors @@ -211,9 +265,14 @@ def cancel_edit(self): def set_fields_enabled(self, enabled): for widget in [ - self.sequence_input, self.sequence_name_input, self.temperature_input, - self.condition_1_input, self.condition_2_input, - self.experimental_run_input, self.secondary_input, self.experiment_input + self.sequence_input, + self.sequence_name_input, + self.temperature_input, + self.condition_1_input, + self.condition_2_input, + self.experimental_run_input, + self.secondary_input, + self.experiment_input, ]: widget.setReadOnly(not enabled) @@ -226,40 +285,66 @@ def clear_fields(self): widget.clear() def validate_inputs(self): - widgets = [self.sequence_input, self.sequence_name_input, self.temperature_input, - self.condition_1_input, self.condition_2_input, - self.experimental_run_input] + widgets = [ + self.sequence_input, + self.sequence_name_input, + self.temperature_input, + self.condition_1_input, + self.condition_2_input, + self.experimental_run_input, + ] return all(w.text().strip() != "" for w in widgets) def add_seq(self): - contig_value = self.contig_input.text().strip() if self.contig_input.isVisible() else self.contig_combo.currentText() + contig_value = ( + self.contig_input.text().strip() + if self.contig_input.isVisible() + else self.contig_combo.currentText() + ) try: - dtr = pd.DataFrame.from_dict({ - "contig": [self.sequence_name_input.text()], - "sequence": [self.sequence_input.text()], - "secondary": [self.secondary_input.text()], - "experiment": [self.experiment_input.text()], - "sequence_name": [self.sequence_name_input.text()], - "sequence_len": [len(self.sequence_input.text().strip())], - "temp": [int(self.temperature_input.text())], - "is_modified": [0 if self.condition_1_input.text() == self.condition_2_input.text() else 1], - "type1": [self.condition_1_input.text()], - "type2": [self.condition_2_input.text()], - "complex": [0], # TODO: Add complex input if needed - "run": [int(self.experimental_run_input.text())] - }, orient='columns') + dtr = pd.DataFrame.from_dict( + { + "contig": [self.sequence_name_input.text()], + "sequence": [self.sequence_input.text()], + "secondary": [self.secondary_input.text()], + "experiment": [self.experiment_input.text()], + "sequence_name": [self.sequence_name_input.text()], + "sequence_len": [len(self.sequence_input.text().strip())], + "temp": [int(self.temperature_input.text())], + "is_modified": [ + ( + 0 + if self.condition_1_input.text() + == self.condition_2_input.text() + else 1 + ) + ], + "type1": [self.condition_1_input.text()], + "type2": [self.condition_2_input.text()], + "complex": [0], # TODO: Add complex input if needed + "run": [int(self.experimental_run_input.text())], + }, + orient="columns", + ) lid = dbins.insert_library(dtr) - if lid==None: + if lid == None: raise Exception("Error processing sequence. Please try again.") dtr["ID"] = lid - #print(dtr.head()) + # print(dtr.head()) self.df = dbsel.select_library_full() # Compose new combo string - contig_value = (str(dtr["ID"][0]) + " " + str(dtr["contig"][0]) + " " + - str(dtr["type1"][0]) + " " + str(dtr["type2"][0])) + contig_value = ( + str(dtr["ID"][0]) + + " " + + str(dtr["contig"][0]) + + " " + + str(dtr["type1"][0]) + + " " + + str(dtr["type2"][0]) + ) self.contig_combo.addItem(contig_value) self.contig_combo.setCurrentText(contig_value) @@ -282,7 +367,7 @@ def update_fields_from_contig(self, contig): except Exception: return if lid in self.df["ID"].values: - row = self.df[self.df['ID'] == lid].iloc[0] + row = self.df[self.df["ID"] == lid].iloc[0] self.lid = row["ID"] self.contig = row["contig"] self.id_input.setText(str(row["ID"])) @@ -324,19 +409,34 @@ def __init__(self, source, lid, contig, signal_path, modification, modification2 def run(self): try: # Select only current contig for upload from tx - cols = ['contig', 'position', 'reference_kmer', 'read_index', - 'event_level_mean', 'event_length', 'event_stdv'] - df = pd.read_csv(self.signal_path, sep='\t', usecols=cols) - df = df[df['contig'] == self.contig] - df['LID'] = self.lid - df['type1'] = self.modification - df['type2'] = self.modification2 + cols = [ + "contig", + "position", + "reference_kmer", + "read_index", + "event_level_mean", + "event_length", + "event_stdv", + ] + df = pd.read_csv(self.signal_path, sep="\t", usecols=cols) + df = df[df["contig"] == self.contig] + df["LID"] = self.lid + df["type1"] = self.modification + df["type2"] = self.modification2 dbins.insert_signal(df) # groups for display plotting - df1 = df.groupby('position').agg(mean_val=('event_level_mean', 'mean')).reset_index() - df2 = df.groupby('position').agg(mean_val=('event_length', 'mean')).reset_index() + df1 = ( + df.groupby("position") + .agg(mean_val=("event_level_mean", "mean")) + .reset_index() + ) + df2 = ( + df.groupby("position") + .agg(mean_val=("event_length", "mean")) + .reset_index() + ) self.finished.emit("signal", df1, df2) except Exception as e: @@ -357,8 +457,9 @@ def __init__(self, source, lid, contig, basecall_path, modification, modificatio def run(self): try: - gpath1, gpath2 = run_basecall.get_modification(self.lid, self.contig, - self.basecall_path, self.modification, plot=True) + gpath1, gpath2 = run_basecall.get_modification( + self.lid, self.contig, self.basecall_path, self.modification, plot=True + ) self.finished.emit("basecall", gpath1, gpath2) except Exception as e: self.error.emit(str(e)) @@ -382,8 +483,14 @@ def __init__(self, sample_section): self.basecall_thread = None self.basecall_worker = None + # Keep strong refs to canvases so we can close/cleanup explicitly + self._canvases = [] + self.init_ui() + # Ensure figures/canvases are cleaned if this widget goes away + self.destroyed.connect(self._cleanup_all) + def init_ui(self): self.setFrameShape(QFrame.Shape.Box) self.setMinimumHeight(400) @@ -434,19 +541,23 @@ def init_ui(self): self.basecall_graphs_layout = QHBoxLayout() default_basecall1 = "default1.png" default_basecall2 = "default2.png" - self.basecall_graph1 = self.create_sample_graph("Basecall: Avg Modification Rates", default_basecall1) - self.basecall_graph2 = self.create_sample_graph("Basecall: Modification by Position", default_basecall2) + self.basecall_graph1 = self.create_sample_graph( + "Basecall: Avg Modification Rates", default_basecall1 + ) + self.basecall_graph2 = self.create_sample_graph( + "Basecall: Modification by Position", default_basecall2 + ) self.basecall_graphs_layout.addWidget(self.basecall_graph1) self.basecall_graphs_layout.addWidget(self.basecall_graph2) layout.addLayout(self.basecall_graphs_layout) - self.setLayout(layout) + self.setLayout(layout) # WHY? def select_basecall_file(self): directory = QFileDialog.getExistingDirectory( self, "Select Basecall Alignment Directory", - options=QFileDialog.Option.ShowDirsOnly + options=QFileDialog.Option.ShowDirsOnly, ) if directory: self.basecall_data_input.setText(directory) @@ -454,22 +565,19 @@ def select_basecall_file(self): def reset_ui(self): # Clear input field self.basecall_data_input.clear() - - # Replace graphs with placeholders - for old_graph in [self.basecall_graph1, self.basecall_graph2]: - self.basecall_graphs_layout.removeWidget(old_graph) - old_graph.setParent(None) - - # Add default graphs - default1 = "default1.png" - default2 = "default2.png" - self.basecall_graph1 = self.create_sample_graph("Basecall: Avg Modification Rates", default1) - self.basecall_graph2 = self.create_sample_graph("Basecall: Modification by Position", default2) - self.basecall_graphs_layout.addWidget(self.basecall_graph1) - self.basecall_graphs_layout.addWidget(self.basecall_graph2) - - # Optionally disable load button until input is provided + # Replace graphs with placeholders (safely) + graph1 = self.create_sample_graph( + "Basecall: Avg Modification Rates", "default1.png", parent=self + ) + graph2 = self.create_sample_graph( + "Basecall: Modification by Position", "default2.png", parent=self + ) + self._replace_graphs( + (self.basecall_graphs_layout, "basecall_graph1", graph1), + (self.basecall_graphs_layout, "basecall_graph2", graph2), + ) self.load_basecall_button.setEnabled(True) + def update_sample_info(self): self.lid = self.sample_section.get_lid() self.contig = self.sample_section.get_contig() @@ -479,6 +587,7 @@ def update_sample_info(self): self.id_display.setText(f"{self.lid} {self.contig}") self.id_display.setReadOnly(True) self.reset_ui() + def load_basecall(self): self.load_basecall_button.setEnabled(False) self.lid = self.sample_section.get_lid() @@ -496,7 +605,7 @@ def load_basecall(self): self.sample_section.get_contig(), self.path, self.modification, - self.modification2 + self.modification2, ) self.basecall_worker.moveToThread(self.basecall_thread) self.basecall_thread.started.connect(self.basecall_worker.run) @@ -509,8 +618,12 @@ def load_basecall(self): def on_worker_finished(self, source, gpath1, gpath2): # Replace old graphs with new ones - new_graph1 = self.create_sample_graph(f"{source.capitalize()}: Avg Modification Rates", gpath1) - new_graph2 = self.create_sample_graph(f"{source.capitalize()}: Modification by Position", gpath2) + new_graph1 = self.create_sample_graph( + f"{source.capitalize()}: Avg Modification Rates", gpath1, parent=self + ) + new_graph2 = self.create_sample_graph( + f"{source.capitalize()}: Modification by Position", gpath2, parent=self + ) # Remove old widgets for old_graph in [self.basecall_graph1, self.basecall_graph2]: @@ -532,12 +645,88 @@ def on_worker_error(self, message): QMessageBox.critical(self, "Error", message) self.load_basecall_button.setEnabled(True) - def create_sample_graph(self, label, gpath=None): - fig = Figure(figsize=(4, 3)) + # ----------------- helpers ----------------- + + def _replace_graphs(self, *triples): + """ + Replace graphs given as (layout, attr_name, new_canvas) triples. + Ensures old canvas is closed/deleted before adding the new one. + """ + for layout, attr_name, new_canvas in triples: + old_canvas = getattr(self, attr_name, None) + if old_canvas is not None: + self._remove_canvas_widget(old_canvas) + setattr(self, attr_name, new_canvas) + layout.addWidget(new_canvas) + + def _remove_canvas_widget(self, canvas: FigureCanvas): + """Close its figure, remove from layout, and delete the widget.""" + try: + self.basecall_graphs_layout.removeWidget(canvas) + except Exception: + pass + + fig = canvas.figure + if fig is not None: + try: + fig.clear() + except Exception: + pass + + try: + import matplotlib.pyplot as plt + + plt.close(fig) + except Exception: + pass + + # schedule widget destruction + try: + canvas.deleteLater() + except Exception: + pass + + # drop our strong reference (if we kept it) + try: + self._canvases.remove(canvas) + except ValueError: + pass + + def cleanup(self): + """Public: call this when explicitly dispose the widget.""" + self._cleanup_all() + + def _cleanup_all(self): + """Call on widget destroy.""" + # Close & delete any known canvases + for canvas in list(self._canvases): + self._remove_canvas_widget(canvas) + self._canvases.clear() + + # Also close currently attached graphs if not already cleared + for graph in [ + getattr(self, "basecall_graph1", None), + getattr(self, "basecall_graph2", None), + ]: + if graph is not None: + self._remove_canvas_widget(graph) + + # Null out just in case + try: + self.basecall_graph1 = None + self.basecall_graph2 = None + except Exception: + pass + + def create_sample_graph(self, label, gpath=None, parent=None): + parent = parent or self # ensure Qt owns it + fig = Figure( + figsize=(4, 3), constrained_layout=True + ) ## TODO: verify constrained_layout canvas = FigureCanvas(fig) ax = fig.add_subplot(111) ax.set_title(label, fontsize=10) - ax.axis('off') + ax.axis("off") image = None @@ -547,10 +736,14 @@ def create_sample_graph(self, label, gpath=None): # Check if it's a default packaged image if basename.startswith("default"): try: - with importlib.resources.files("DashML.GUI").joinpath(basename).open("rb") as f: - image = mpimg.imread(f, format='png') + with importlib.resources.files("DashML.GUI").joinpath( + basename + ).open("rb") as f: + image = mpimg.imread(f, format="png") except FileNotFoundError: - print(f"Default image '{basename}' not found in DashML.GUI package.") + print( + f"Default image '{basename}' not found in DashML.GUI package." + ) except Exception as e: print(f"Error loading default image '{basename}': {e}") elif os.path.exists(gpath): @@ -560,13 +753,15 @@ def create_sample_graph(self, label, gpath=None): print(f"Error loading external image '{gpath}': {e}") if image is not None: - ax.imshow(image, aspect='auto') + ax.imshow(image, aspect="auto") else: - ax.text(0.5, 0.5, "No Data", ha='center', va='center') + ax.text(0.5, 0.5, "No Data", ha="center", va="center") - fig.tight_layout() + # fig.tight_layout() + self._canvases.append(canvas) # track for cleanup return canvas + class LoadSignal(QFrame): request_start_worker = pyqtSignal(dict) # emit parameters as dict @@ -584,8 +779,12 @@ def __init__(self, sample_section): self.signal_thread = None self.signal_worker = None + self._canvases = [] + self.init_ui() + self.destroyed.connect(self._cleanup_all) + def init_ui(self): self.setFrameShape(QFrame.Shape.Box) self.setMinimumHeight(400) @@ -637,8 +836,12 @@ def init_ui(self): self.signal_graphs_layout = QHBoxLayout() default_signal1 = "default3.png" default_signal2 = "default4.png" - self.signal_graph1 = self.create_sample_graph("Signal: Avg Signal", default_signal1) - self.signal_graph2 = self.create_sample_graph("Signal: Avg Dwell Time", default_signal2) + self.signal_graph1 = self.create_sample_graph( + "Signal: Avg Signal", default_signal1 + ) + self.signal_graph2 = self.create_sample_graph( + "Signal: Avg Dwell Time", default_signal2 + ) self.signal_graphs_layout.addWidget(self.signal_graph1) self.signal_graphs_layout.addWidget(self.signal_graph2) @@ -667,10 +870,18 @@ def reset_ui(self): # Add default graphs default1 = "default3.png" default2 = "default4.png" - self.signal_graph1 = self.create_sample_graph("Signal: Avg Signal", gpath=default1) - self.signal_graph2 = self.create_sample_graph("Signal: Avg Dwell Time", gpath=default2) - self.signal_graphs_layout.addWidget(self.signal_graph1) - self.signal_graphs_layout.addWidget(self.signal_graph2) + graph1 = self.create_sample_graph( + "Signal: Avg Signal", gpath=default1, parent=self + ) + graph2 = self.create_sample_graph( + "Signal: Avg Dwell Time", gpath=default2, parent=self + ) + + # safely close/remove old, add new + self._replace_graphs( + (self.signal_graphs_layout, "signal_graph1", graph1), + (self.signal_graphs_layout, "signal_graph2", graph2), + ) # Optionally disable load button until input is provided self.load_signal_button.setEnabled(True) @@ -684,6 +895,7 @@ def update_sample_info(self): self.id_display.setText(f"{self.lid} {self.contig}") self.id_display.setReadOnly(True) self.reset_ui() + def load_signal(self): self.load_signal_button.setEnabled(False) self.lid = self.sample_section.get_lid() @@ -694,7 +906,14 @@ def load_signal(self): self.path = self.signal_data_input.text() self.signal_thread = QThread() - self.signal_worker = SignalWorker(self.source, self.lid, self.contig, self.path, self.modification, self.modification2) + self.signal_worker = SignalWorker( + self.source, + self.lid, + self.contig, + self.path, + self.modification, + self.modification2, + ) self.signal_worker.moveToThread(self.signal_thread) self.signal_thread.started.connect(self.signal_worker.run) self.signal_worker.finished.connect(self.on_worker_finished) @@ -706,21 +925,18 @@ def load_signal(self): def on_worker_finished(self, source, gpath1, gpath2): # Replace old graphs with new ones - new_graph1 = self.create_sample_graph(f"{source.capitalize()}: Avg Signal", gpath1) - new_graph2 = self.create_sample_graph(f"{source.capitalize()}: Avg Dwell Time", gpath2) - - # Remove old widgets - for old_graph in [self.signal_graph1, self.signal_graph2]: - self.signal_graphs_layout.removeWidget(old_graph) - old_graph.setParent(None) - - # Add new widgets - self.signal_graphs_layout.addWidget(new_graph1) - self.signal_graphs_layout.addWidget(new_graph2) + new_graph1 = self.create_sample_graph( + f"{source.capitalize()}: Avg Signal", gpath1, parent=self + ) + new_graph2 = self.create_sample_graph( + f"{source.capitalize()}: Avg Dwell Time", gpath2, parent=self + ) - # Save references - self.signal_graph1 = new_graph1 - self.signal_graph2 = new_graph2 + # safely swap graphs + self._replace_graphs( + (self.signal_graphs_layout, "signal_graph1", new_graph1), + (self.signal_graphs_layout, "signal_graph2", new_graph2), + ) QMessageBox.information(self, "Success", "Sample added successfully.") self.load_signal_button.setEnabled(True) @@ -729,9 +945,12 @@ def on_worker_error(self, message): QMessageBox.critical(self, "Error", message) self.load_signal_button.setEnabled(True) - def create_sample_graph(self, label, gpath=None): - fig = Figure(figsize=(4, 3)) + def create_sample_graph(self, label, gpath=None, parent=None): + parent = parent or self # ensure Qt owns it + + fig = Figure(figsize=(4, 3), constrained_layout=True) canvas = FigureCanvas(fig) + canvas.setParent(parent) ax = fig.add_subplot(111) ax.set_title(label, fontsize=10) @@ -740,18 +959,20 @@ def create_sample_graph(self, label, gpath=None): # Handle DataFrame case first if isinstance(gpath, pd.DataFrame): if not gpath.empty: - if 'position' in gpath.columns and 'mean_val' in gpath.columns: - ax.scatter(gpath['position'], gpath['mean_val'], s=10, alpha=0.7) + if "position" in gpath.columns and "mean_val" in gpath.columns: + ax.scatter(gpath["position"], gpath["mean_val"], s=10, alpha=0.7) ax.set_xlabel("Position", fontsize=10) ax.set_ylabel("Mean Value", fontsize=10) - ax.tick_params(axis='both', which='major', labelsize=10) - ax.tick_params(axis='both', which='minor', labelsize=5) + ax.tick_params(axis="both", which="major", labelsize=10) + ax.tick_params(axis="both", which="minor", labelsize=5) else: - ax.text(0.5, 0.5, "Invalid DataFrame format", ha='center', va='center') + ax.text( + 0.5, 0.5, "Invalid DataFrame format", ha="center", va="center" + ) ax.set_xticks([]) ax.set_yticks([]) else: - ax.text(0.5, 0.5, "Empty DataFrame", ha='center', va='center') + ax.text(0.5, 0.5, "Empty DataFrame", ha="center", va="center") ax.set_xticks([]) ax.set_yticks([]) elif isinstance(gpath, str): @@ -760,10 +981,14 @@ def create_sample_graph(self, label, gpath=None): # Check if default image if basename.startswith("default"): try: - with importlib.resources.files("DashML.GUI").joinpath(basename).open("rb") as f: - image = mpimg.imread(f, format='png') + with importlib.resources.files("DashML.GUI").joinpath( + basename + ).open("rb") as f: + image = mpimg.imread(f, format="png") except FileNotFoundError: - print(f"Default image '{basename}' not found in DashML.GUI package.") + print( + f"Default image '{basename}' not found in DashML.GUI package." + ) except Exception as e: print(f"Error loading default image '{basename}': {e}") elif os.path.exists(gpath): @@ -773,20 +998,101 @@ def create_sample_graph(self, label, gpath=None): print(f"Error loading external image '{gpath}': {e}") if image is not None: - ax.imshow(image, aspect='auto') - ax.axis('off') + ax.imshow(image, aspect="auto") + ax.axis("off") else: - ax.text(0.5, 0.5, "Image not found", ha='center', va='center') + ax.text(0.5, 0.5, "Image not found", ha="center", va="center") ax.set_xticks([]) ax.set_yticks([]) else: - ax.text(0.5, 0.5, "No Data", ha='center', va='center') + ax.text(0.5, 0.5, "No Data", ha="center", va="center") ax.set_xticks([]) ax.set_yticks([]) - fig.tight_layout() + # fig.tight_layout() + + self._canvases.append(canvas) # track for cleanup return canvas + def _remove_canvas_widget(self, canvas: FigureCanvas): + """Close its figure, remove from layout, and delete the widget.""" + try: + self.signal_graphs_layout.removeWidget(canvas) + except Exception: + pass + + fig = canvas.figure + if fig is not None: + try: + fig.clear() + except Exception: + pass + try: + import matplotlib.pyplot as plt + + plt.close(fig) + except Exception: + pass + + # schedule widget destruction + try: + canvas.deleteLater() + except Exception: + pass + + try: + self._canvases.remove(canvas) + except Exception: + pass + + def _replace_graphs(self, *triples): + """ + Replace graphs given as (layout, attr_name, new_canvas) triples. + Ensures old canvas is closed/deleted before adding the new one. + """ + for layout, attr_name, new_canvas in triples: + old_canvas = getattr(self, attr_name, None) + if old_canvas is not None: + self._remove_canvas_widget(old_canvas) + setattr(self, attr_name, new_canvas) + layout.addWidget(new_canvas) + + def cleanup(self): + """Public: call this when explicitly dispose the widget.""" + self._cleanup_all() + + def _cleanup_all(self): + """Close figures, delete canvases, and stop threads on widget destruction.""" + # Stop/cleanup thread safely + try: + if self.signal_thread and self.signal_thread.isRunning(): + # If worker supports interruption, request it here: + # try: self.signal_worker.requestInterruption() except: pass + self.signal_thread.quit() + self.signal_thread.wait(2000) + except Exception: + pass + + # Remove and delete known canvases + for c in list(self._canvases): + self._remove_canvas_widget(c) + self._canvases.clear() + + # Also remove currently referenced graphs if still present + for graph in [ + getattr(self, "signal_graph1", None), + getattr(self, "signal_graph2", None), + ]: + if graph is not None: + self._remove_canvas_widget(graph) + + # Null out just in case + try: + self.signal_graph1 = None + self.signal_graph2 = None + except Exception: + pass + class PredictSection(QFrame): def __init__(self): @@ -799,9 +1105,12 @@ def __init__(self): self.contig_unmod = None self.lid_mod = None self.contig_mod = None - self.df= None + self.df = None + + self._canvases = [] self.init_ui() + self.destroyed.connect(self._cleanup_all) def init_ui(self): layout = QVBoxLayout() @@ -814,31 +1123,41 @@ def init_ui(self): # Dropdowns dropdown_layout = QHBoxLayout() self.unmod_combo = QComboBox() - self.mod_combo= QComboBox() - self.unmod_combo.setSizePolicy(QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Fixed) - self.mod_combo.setSizePolicy(QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Fixed) + self.mod_combo = QComboBox() + self.unmod_combo.setSizePolicy( + QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Fixed + ) + self.mod_combo.setSizePolicy( + QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Fixed + ) self.df = dbsel.select_library_full() - self.df['combo'] = (self.df["ID"].astype(str) + " " + self.df["contig"] + " " - + self.df["type1"] + " " + self.df["type2"]) - self.df_unmod = self.df.loc[self.df['is_modified'] == 0] - self.df_mod = self.df.loc[self.df['is_modified'] == 1] - self.unmod_combo.addItems(self.df_unmod['combo']) + self.df["combo"] = ( + self.df["ID"].astype(str) + + " " + + self.df["contig"] + + " " + + self.df["type1"] + + " " + + self.df["type2"] + ) + self.df_unmod = self.df.loc[self.df["is_modified"] == 0] + self.df_mod = self.df.loc[self.df["is_modified"] == 1] + self.unmod_combo.addItems(self.df_unmod["combo"]) # Select first item properly if not self.df_unmod.empty: - self.unmod_combo.setCurrentText(self.df_unmod['combo'].iloc[0]) + self.unmod_combo.setCurrentText(self.df_unmod["combo"].iloc[0]) self.lid_unmod = self.df_unmod["ID"].iloc[0] self.contig_unmod = self.df_unmod["contig"].iloc[0] - self.mod_combo.addItems(self.df_mod['combo']) + self.mod_combo.addItems(self.df_mod["combo"]) # Select first item properly if not self.df_mod.empty: - self.mod_combo.setCurrentText(self.df_mod['combo'].iloc[0]) + self.mod_combo.setCurrentText(self.df_mod["combo"].iloc[0]) self.lid_mod = self.df_mod["ID"].iloc[0] self.contig_mod = self.df_mod["contig"].iloc[0] - dropdown_layout.addWidget(QLabel("Select Unmodified:")) dropdown_layout.addWidget(self.unmod_combo) dropdown_layout.addWidget(QLabel("Select Modified:")) @@ -863,8 +1182,8 @@ def init_ui(self): self.graph1_layout = QHBoxLayout() self.graph2_layout = QHBoxLayout() - self.graph1 = self.create_sample_graph("default5.png") - self.graph2 = self.create_sample_graph("default6.png") + self.graph1 = self.create_sample_graph("default5.png", parent=self) + self.graph2 = self.create_sample_graph("default6.png", parent=self) self.graph1_layout.addWidget(self.graph1) self.graph2_layout.addWidget(self.graph2) @@ -893,14 +1212,19 @@ def populate_dropdowns(self): self.mod_combo.update() # force repaint def refresh_data(self): - #print("PredictSection: refresh_data() called") + # print("PredictSection: refresh_data() called") self.df = dbsel.select_library_full() - self.df['combo'] = ( - self.df["ID"].astype(str) + " " + self.df["contig"] + " " - + self.df["type1"] + " " + self.df["type2"] + self.df["combo"] = ( + self.df["ID"].astype(str) + + " " + + self.df["contig"] + + " " + + self.df["type1"] + + " " + + self.df["type2"] ) - self.df_unmod = self.df.loc[self.df['is_modified'] == 0] - self.df_mod = self.df.loc[self.df['is_modified'] == 1] + self.df_unmod = self.df.loc[self.df["is_modified"] == 0] + self.df_mod = self.df.loc[self.df["is_modified"] == 1] # Clear and repopulate unmodified combo self.unmod_combo.clear() @@ -918,7 +1242,8 @@ def refresh_data(self): self.lid_mod = self.df_mod["ID"].iloc[0] self.contig_mod = self.df_mod["contig"].iloc[0] - def create_sample_graph(self, image_path, label=None): + def create_sample_graph(self, image_path, label=None, parent=None): + parent = parent or self image = None aspect_ratio = 1.5 # Default fallback @@ -930,12 +1255,16 @@ def create_sample_graph(self, image_path, label=None): if basename.startswith("default"): try: # Use importlib.resources to get the image inside the package - with importlib.resources.files(DashML.GUI).joinpath(basename).open("rb") as f: - image = mpimg.imread(f, format='png') + with importlib.resources.files(DashML.GUI).joinpath(basename).open( + "rb" + ) as f: + image = mpimg.imread(f, format="png") height, width = image.shape[:2] aspect_ratio = width / height if height != 0 else 1.5 except FileNotFoundError: - print(f"Default image '{basename}' not found in DashML.GUI package.") + print( + f"Default image '{basename}' not found in DashML.GUI package." + ) except Exception as e: print(f"Error loading default image '{basename}': {e}") else: @@ -953,22 +1282,25 @@ def create_sample_graph(self, image_path, label=None): # Dynamically set figure size based on aspect ratio fig = Figure(figsize=(6 * aspect_ratio, 6)) canvas = FigureCanvas(fig) + canvas.setParent(parent) canvas.setSizePolicy(QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Expanding) ax = fig.add_subplot(111) - ax.axis('off') + ax.axis("off") if label is not None: ax.set_title(label, fontsize=12, pad=10) if image is not None: - ax.imshow(image, aspect='auto') + ax.imshow(image, aspect="auto") else: - ax.text(0.5, 0.5, "No image", ha='center', va='center') + ax.text(0.5, 0.5, "No image", ha="center", va="center") fig.tight_layout() + self._canvases.append(canvas) return canvas + def run_mod(self): self.run_prediction("Predict Modifications") @@ -984,8 +1316,12 @@ def run_prediction(self, label): self.lid_mod = self.mod_combo.currentText().strip().split()[0] self.contig_mod = self.mod_combo.currentText().strip().split()[1] - seq_unmod = self.df_unmod.loc[self.df_unmod['ID'] == int(self.lid_unmod), 'sequence'].unique()[0] - seq_mod = self.df_mod.loc[self.df_mod['ID'] == int(self.lid_mod), 'sequence'].unique()[0] + seq_unmod = self.df_unmod.loc[ + self.df_unmod["ID"] == int(self.lid_unmod), "sequence" + ].unique()[0] + seq_mod = self.df_mod.loc[ + self.df_mod["ID"] == int(self.lid_mod), "sequence" + ].unique()[0] if label == "Predict Modifications": if str(seq_mod) == str(seq_unmod): @@ -993,27 +1329,31 @@ def run_prediction(self, label): unmod_lids=self.lid_unmod, lids=self.lid_mod, continue_reads=False, - vienna=use_base_pairing + vienna=use_base_pairing, ) if image1 is not None: # Use singleShot if you want delayed UI update - QTimer.singleShot(3000, lambda: self.update_graphs(image1, image2, label)) + QTimer.singleShot( + 3000, lambda: self.update_graphs(image1, image2, label) + ) else: QMessageBox.critical( self, "Prediction Error", - "Prediction failed. Please check the input or try again." + "Prediction failed. Please check the input or try again.", ) else: QMessageBox.critical( self, "Prediction Error", - "Please select the unmod and modified versions of the same sequence and try again." + "Please select the unmod and modified versions of the same sequence and try again.", ) except Exception as e: print(f"Prediction failed: {e}") - QMessageBox.critical(self, "Prediction Error", f"An unexpected error occurred: {e}") + QMessageBox.critical( + self, "Prediction Error", f"An unexpected error occurred: {e}" + ) finally: # Always re-enable UI controls self.btn_mod.setEnabled(True) @@ -1023,33 +1363,88 @@ def run_prediction(self, label): def update_graphs(self, image_path, image_path2, label): if label == "Predict Modifications": # Replace old graphs with new ones - new_graph1 = self.create_sample_graph(image_path, label="Predicted Modifications") - - # Remove old graph1 from its layout - self.graph1_layout.removeWidget(self.graph1) - self.graph1.setParent(None) - - # Add the new graph1 to the same layout - self.graph1_layout.addWidget(new_graph1) - self.graph1 = new_graph1 - + new_graph1 = self.create_sample_graph( + image_path, label="Predicted Modifications", parent=self + ) # Replace old graphs with new ones - new_graph2 = self.create_sample_graph(image_path2, label="Predicted Modifications") - - # Remove old graph2 from its layout - self.graph2_layout.removeWidget(self.graph2) - self.graph2.setParent(None) + new_graph2 = self.create_sample_graph( + image_path2, label="Predicted Modifications", parent=self + ) - # Add the new graph2 to the same layout - self.graph2_layout.addWidget(new_graph2) - self.graph2 = new_graph2 + # Swap safely + self._replace_graphs( + (self.graph1_layout, "graph1", new_graph1), + (self.graph2_layout, "graph2", new_graph2), + ) # Show success message - QMessageBox.information(self, "Prediction Complete", f"{label} prediction finished.") + QMessageBox.information( + self, "Prediction Complete", f"{label} prediction finished." + ) # Re-enable both buttons self.btn_mod.setEnabled(True) + def _remove_canvas_widget(self, canvas: FigureCanvas): + """Close its figure, remove from layout, and delete the widget.""" + # remove from any known layouts + try: + self.graph1_layout.removeWidget(canvas) + except Exception: + pass + try: + self.graph2_layout.removeWidget(canvas) + except Exception: + pass + + # close its figure + fig = getattr(canvas, "figure", None) + if fig is not None: + try: + fig.clear() + except Exception: + pass + try: + import matplotlib.pyplot as plt + + plt.close(fig) + except Exception: + pass + + # schedule widget deletion & drop references + try: + canvas.deleteLater() + except Exception: + pass + try: + self._canvases.remove(canvas) + except ValueError: + pass + + def _replace_graphs(self, *triples): + """ + Replace graphs given as (layout, attr_name, new_canvas) triples. + Ensures old canvas is closed/deleted before adding the new one. + """ + for layout, attr_name, new_canvas in triples: + old_canvas = getattr(self, attr_name, None) + if old_canvas is not None: + self._remove_canvas_widget(old_canvas) + setattr(self, attr_name, new_canvas) + layout.addWidget(new_canvas) + + def _cleanup_all(self): + """Close figures and delete canvases when this widget is destroyed.""" + for c in list(self._canvases): + self._remove_canvas_widget(c) + self._canvases.clear() + + # Also clear currently referenced graph widgets if still present + for name in ("graph1", "graph2"): + c = getattr(self, name, None) + if c is not None: + self._remove_canvas_widget(c) + setattr(self, name, None) class CreateLandscapeSection(QFrame): @@ -1064,10 +1459,20 @@ def __init__(self, parent=None): self.lid_mod = None self.contig_mod = None self.df = dbsel.select_library_full() - self.df['combo'] = (self.df["ID"].astype(str) + " " + self.df["contig"] + " " - + self.df["type1"] + " " + self.df["type2"]) - self.df_unmod = self.df.loc[self.df['is_modified'] == 0] - self.df_mod = self.df.loc[self.df['is_modified'] == 1] + self.df["combo"] = ( + self.df["ID"].astype(str) + + " " + + self.df["contig"] + + " " + + self.df["type1"] + + " " + + self.df["type2"] + ) + self.df_unmod = self.df.loc[self.df["is_modified"] == 0] + self.df_mod = self.df.loc[self.df["is_modified"] == 1] + + self._canvases = [] + self.destroyed.connect(self._cleanup_all) self.init_ui() @@ -1085,19 +1490,23 @@ def init_ui(self): dropdown_layout = QHBoxLayout() self.unmod_combo = QComboBox() self.mod_combo = QComboBox() - self.unmod_combo.setSizePolicy(QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Fixed) - self.mod_combo.setSizePolicy(QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Fixed) - self.unmod_combo.addItems(self.df_unmod['combo']) + self.unmod_combo.setSizePolicy( + QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Fixed + ) + self.mod_combo.setSizePolicy( + QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Fixed + ) + self.unmod_combo.addItems(self.df_unmod["combo"]) # Select first item properly if not self.df_unmod.empty: - self.unmod_combo.setCurrentText(self.df_unmod['combo'].iloc[0]) + self.unmod_combo.setCurrentText(self.df_unmod["combo"].iloc[0]) self.lid_unmod = self.df_unmod["ID"].iloc[0] self.contig_unmod = self.df_unmod["contig"].iloc[0] - self.mod_combo.addItems(self.df_mod['combo']) + self.mod_combo.addItems(self.df_mod["combo"]) # Select first item properly if not self.df_mod.empty: - self.mod_combo.setCurrentText(self.df_mod['combo'].iloc[0]) + self.mod_combo.setCurrentText(self.df_mod["combo"].iloc[0]) self.lid_mod = self.df_mod["ID"].iloc[0] self.contig_mod = self.df_mod["contig"].iloc[0] @@ -1126,15 +1535,17 @@ def init_ui(self): self.graphs_layout = QHBoxLayout() self.graphs_layout.setSpacing(6) - self.graph1 = self.create_sample_graph("default7.png") - self.graph2 = self.create_sample_graph("default7.png") - self.graph3 = self.create_sample_graph("default7.png") - self.graph4 = self.create_sample_graph("default8.png") + graph1 = self.create_sample_graph("default7.png", parent=self) + graph2 = self.create_sample_graph("default7.png", parent=self) + graph3 = self.create_sample_graph("default7.png", parent=self) + graph4 = self.create_sample_graph("default8.png", parent=self) - self.graphs_layout.addWidget(self.graph1) - self.graphs_layout.addWidget(self.graph2) - self.graphs_layout.addWidget(self.graph3) - self.graphs_layout.addWidget(self.graph4) + self._replace_graphs( + (self.graphs_layout, "graph1", graph1), + (self.graphs_layout, "graph2", graph2), + (self.graphs_layout, "graph3", graph3), + (self.graphs_layout, "graph4", graph4), + ) layout.addLayout(self.graphs_layout) self.setLayout(layout) @@ -1154,12 +1565,17 @@ def populate_dropdowns(self): def refresh_data(self): # print("PredictSection: refresh_data() called") self.df = dbsel.select_library_full() - self.df['combo'] = ( - self.df["ID"].astype(str) + " " + self.df["contig"] + " " - + self.df["type1"] + " " + self.df["type2"] + self.df["combo"] = ( + self.df["ID"].astype(str) + + " " + + self.df["contig"] + + " " + + self.df["type1"] + + " " + + self.df["type2"] ) - self.df_unmod = self.df.loc[self.df['is_modified'] == 0] - self.df_mod = self.df.loc[self.df['is_modified'] == 1] + self.df_unmod = self.df.loc[self.df["is_modified"] == 0] + self.df_mod = self.df.loc[self.df["is_modified"] == 1] # Clear and repopulate unmodified combo self.unmod_combo.clear() @@ -1185,14 +1601,18 @@ def create_landscape(self): self.contig_unmod = self.unmod_combo.currentText().strip().split()[1] self.lid_mod = self.mod_combo.currentText().strip().split()[0] self.contig_mod = self.mod_combo.currentText().strip().split()[1] - seq_unmod = self.df_unmod.loc[self.df_unmod['ID'] == int(self.lid_unmod), 'sequence'].unique()[0] - seq_mod = self.df_mod.loc[self.df_mod['ID'] == int(self.lid_mod), 'sequence'].unique()[0] + seq_unmod = self.df_unmod.loc[ + self.df_unmod["ID"] == int(self.lid_unmod), "sequence" + ].unique()[0] + seq_mod = self.df_mod.loc[ + self.df_mod["ID"] == int(self.lid_mod), "sequence" + ].unique()[0] if str(seq_mod) == str(seq_unmod): images = landscape.run_landscape( unmod_lid=self.lid_unmod, lid=self.lid_mod, - optimize_clusters=optimize_clusters + optimize_clusters=optimize_clusters, ) if images is not None: @@ -1201,49 +1621,53 @@ def create_landscape(self): QMessageBox.critical( self, "Landscape Error", - "Landscape failed. Please check the input or try again." + "Landscape failed. Please check the input or try again.", ) else: QMessageBox.critical( self, "Landscape Error", - "Please select the unmod and modified versions of the same sequence and try again." + "Please select the unmod and modified versions of the same sequence and try again.", ) except Exception as e: print(f"Landscape creation failed: {e}") - QMessageBox.critical(self, "Landscape Error", f"An unexpected error occurred: {e}") + QMessageBox.critical( + self, "Landscape Error", f"An unexpected error occurred: {e}" + ) finally: self.create_button.setEnabled(True) def update_graphs(self, images): - # Remove old graphs from layout and UI - for graph in [self.graph1, self.graph2, self.graph3, self.graph4]: - self.graphs_layout.removeWidget(graph) - graph.setParent(None) - - # Create new graphs using updated images - self.graph1 = self.create_sample_graph(images[0], label="HeatMap") - self.graph2 = self.create_sample_graph(images[1], label="Dendrogram") - self.graph3 = self.create_sample_graph(images[2], label="Read Corr") - self.graph4 = self.create_sample_graph(images[3], label="Cluster Opt.") - - # Add them back to the layout - self.graphs_layout.addWidget(self.graph1) - self.graphs_layout.addWidget(self.graph2) - self.graphs_layout.addWidget(self.graph3) - self.graphs_layout.addWidget(self.graph4) + # Build new canvases + graph1 = self.create_sample_graph(images[0], label="HeatMap", parent=self) + graph2 = self.create_sample_graph(images[1], label="Dendrogram", parent=self) + graph3 = self.create_sample_graph(images[2], label="Read Corr", parent=self) + graph4 = self.create_sample_graph(images[3], label="Cluster Opt.", parent=self) + + # Swap safely + self._replace_graphs( + (self.graphs_layout, "graph1", graph1), + (self.graphs_layout, "graph2", graph2), + (self.graphs_layout, "graph3", graph3), + (self.graphs_layout, "graph4", graph4), + ) # Done - QMessageBox.information(self, "Landscape Complete", "Landscape analysis finished.") + QMessageBox.information( + self, "Landscape Complete", "Landscape analysis finished." + ) self.create_button.setEnabled(True) - def create_sample_graph(self, image_path, label=None): + def create_sample_graph(self, image_path, label=None, parent=None): + parent = parent or self # Check if this is a default image if os.path.basename(image_path).startswith("default"): try: # Load from DashML.GUI package - with importlib.resources.files("DashML.GUI").joinpath(os.path.basename(image_path)).open("rb") as f: - image = mpimg.imread(f, format='png') + with importlib.resources.files("DashML.GUI").joinpath( + os.path.basename(image_path) + ).open("rb") as f: + image = mpimg.imread(f, format="png") except FileNotFoundError: print(f"Default image {image_path} not found in DashML.GUI package.") image = None @@ -1258,25 +1682,85 @@ def create_sample_graph(self, image_path, label=None): # Fixed figure size to ensure all graphs have the same size fig = Figure(figsize=(6, 6)) # 6x6 inches fixed size canvas = FigureCanvas(fig) + canvas.setParent(parent) canvas.setSizePolicy(QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Expanding) ax = fig.add_subplot(111) - ax.axis('off') # Hide axes lines and ticks + ax.axis("off") # Hide axes lines and ticks if label is not None: - ax.set_title(label, fontsize=12, fontweight='bold', pad=10) + ax.set_title(label, fontsize=12, fontweight="bold", pad=10) if image is not None: - ax.imshow(image, aspect='auto') # Scale naturally + ax.imshow(image, aspect="auto") # Scale naturally else: - ax.text(0.5, 0.5, "No image", ha='center', va='center') + ax.text(0.5, 0.5, "No image", ha="center", va="center") # Adjust layout to make room for title fig.tight_layout() fig.subplots_adjust(top=0.85) # Adjust top margin so title doesn't get cut off + self._canvases.append(canvas) return canvas + def _remove_canvas_widget(self, canvas: FigureCanvas): + """Close its figure, remove from layout, and delete the widget.""" + try: + self.graphs_layout.removeWidget(canvas) + except Exception: + pass + + fig = getattr(canvas, "figure", None) + if fig is not None: + try: + fig.clear() + except Exception: + pass + try: + import matplotlib.pyplot as plt + + plt.close(fig) + except Exception: + pass + + try: + canvas.deleteLater() + except Exception: + pass + try: + self._canvases.remove(canvas) + except ValueError: + pass + + def _replace_graphs(self, *triples): + """ + Replace graphs given as (layout, attr_name, new_canvas) triples. + Ensures old canvas is closed/deleted before adding the new one. + """ + for layout, attr_name, new_canvas in triples: + old_canvas = getattr(self, attr_name, None) + if old_canvas is not None: + self._remove_canvas_widget(old_canvas) + setattr(self, attr_name, new_canvas) + layout.addWidget(new_canvas) + + def cleanup(self): + """Public: call this when explicitly dispose the widget.""" + self._cleanup_all() + + def _cleanup_all(self): + """Close figures and delete canvases when this widget is destroyed.""" + for c in list(self._canvases): + self._remove_canvas_widget(c) + self._canvases.clear() + + # Also clear currently referenced graph widgets if still present + for name in ("graph1", "graph2", "graph3", "graph4"): + c = getattr(self, name, None) + if c is not None: + self._remove_canvas_widget(c) + setattr(self, name, None) + class MainApp(QWidget): def __init__(self): @@ -1289,17 +1773,19 @@ def __init__(self): # Pass shared sample section to LoadBasecall self.load_basecall_section = LoadBasecall(self.sample_section) - self.load_basecall_section.request_start_worker.connect(self.handle_start_worker) + self.load_basecall_section.request_start_worker.connect( + self.handle_start_worker + ) # Pass shared sample section to LoadSignal self.load_signal_section = LoadSignal(self.sample_section) self.load_signal_section.request_start_worker.connect(self.handle_start_worker) - #predict + # predict self.predict_section = PredictSection() self.sample_section.sample_changed.connect(self.predict_section.refresh_data) - #landscape + # landscape self.landscape_section = CreateLandscapeSection() self.sample_section.sample_changed.connect(self.landscape_section.refresh_data) @@ -1311,6 +1797,9 @@ def __init__(self): # self.signal_thread = None # self.signal_worker = None + self._active_threads = set() # track active threads + self._active_workers = set() # track active workers + def init_ui(self): scroll = QScrollArea() container = QWidget() @@ -1326,7 +1815,9 @@ def init_ui(self): container.setLayout(layout) # Apply responsive sizing - container.setSizePolicy(QSizePolicy.Policy.Preferred, QSizePolicy.Policy.Maximum) + container.setSizePolicy( + QSizePolicy.Policy.Preferred, QSizePolicy.Policy.Maximum + ) container.adjustSize() scroll.setWidget(container) @@ -1339,31 +1830,54 @@ def init_ui(self): def handle_start_worker(self, params): thread = QThread() - source = params.get('source') - if source == 'basecall': + source = params.get("source") + if source == "basecall": worker = BasecallWorker( - params['lid'], params['contig'], params['basecall_path'], - params['modification'], params['modification2'] + source, # preserve signature + params["lid"], + params["contig"], + params["basecall_path"], + params["modification"], + params["modification2"], ) - elif source == 'signal': + elif source == "signal": worker = SignalWorker( - params['lid'], params['contig'], params['signal_path'], - params['modification'], params['modification2'] + source, # preserve signature + params["lid"], + params["contig"], + params["signal_path"], + params["modification"], + params["modification2"], ) else: print(f"Unknown source: {source}") return - self._start_worker(thread, worker) + self._start_worker(thread, worker, source) - def _start_worker(self, thread, worker): + def _start_worker(self, thread, worker, source): worker.moveToThread(thread) thread.started.connect(worker.run) worker.finished.connect(self.on_worker_finished) worker.error.connect(self.on_worker_error) worker.finished.connect(thread.quit) thread.finished.connect(thread.deleteLater) + worker.finished.connect(worker.deleteLater) + + # Track active threads/workers for Cleanup + self._active_threads.add(thread) + self._active_workers.add(worker) + + # Remove from tracking once finished + def _cleanup_refs(): + self._active_threads.discard(thread) + self._active_workers.discard(worker) + + thread.finished.connect(_cleanup_refs) + + worker._source = source # for debug + thread.start() def on_worker_finished(self, source, gpath1, gpath2): @@ -1372,28 +1886,45 @@ def on_worker_finished(self, source, gpath1, gpath2): if source == "basecall": layout = self.load_basecall_section.basecall_graphs_layout - old1, old2 = self.load_basecall_section.basecall_graph1, self.load_basecall_section.basecall_graph2 - new_graph1 = self.load_basecall_section.create_sample_graph(f"{source.capitalize()}: Avg Modification Rates", - gpath1) - new_graph2 = self.load_basecall_section.create_sample_graph(f"{source.capitalize()}: Modification by Position", - gpath2) - self.load_basecall_section.basecall_graph1 = new_graph1 - self.load_basecall_section.basecall_graph2 = new_graph2 - else: + + new_graph1 = self.load_basecall_section.create_sample_graph( + f"{source.capitalize()}: Avg Modification Rates", + gpath1, + parent=self.load_basecall_section, + ) + new_graph2 = self.load_basecall_section.create_sample_graph( + f"{source.capitalize()}: Modification by Position", + gpath2, + parent=self.load_basecall_section, + ) + self.load_basecall_section._replace_graphs( + (layout, "basecall_graph1", new_graph1), + (layout, "basecall_graph2", new_graph2), + ) + elif source == "signal": layout = self.load_signal_section.signal_graphs_layout - old1, old2 = self.load_signal_section.signal_graph1, self.load_basecall_section.signal_graph2 + new_graph1 = self.load_signal_section.create_sample_graph( - f"{source.capitalize()}: Average Signal Rates by Position", gpath1) + f"{source.capitalize()}: Average Signal Rates by Position", + gpath1, + parent=self.load_signal_section, + ) new_graph2 = self.load_signal_section.create_sample_graph( - f"{source.capitalize()}: Average Dwell Time by Position", gpath2) - self.load_signal_section.signal_graph1 = new_graph1 - self.load_signal_section.signal_graph2 = new_graph2 - - layout.replaceWidget(old1, new_graph1) - layout.replaceWidget(old2, new_graph2) + f"{source.capitalize()}: Average Dwell Time by Position", + gpath2, + parent=self.load_signal_section, + ) + # Swap safely + self.load_signal_section._replace_graphs( + (layout, "signal_graph1", new_graph1), + (layout, "signal_graph2", new_graph2), + ) + else: + print(f"[WARNING] Unknown source on finish: {source}") + return - old1.setParent(None) - old2.setParent(None) + # Do NOT call layout.replaceWidget + setParent(None); the section + # helper already closes figures, deletes widgets, and updates refs. # Clear references to allow garbage collection if source == "basecall": @@ -1407,12 +1938,50 @@ def on_worker_error(self, message): self.load_basecall_section.load_basecall_button.setEnabled(True) self.load_signal_section.load_signal_button.setEnabled(True) + # Add more explicit error handling if needed + try: + QMessageBox.critical(self, "Worker Error", message) + except Exception: + print(f"Worker error: {message}") + + def closeEvent(self, event): + """ + Ask sections to cleanup canvases/threads if they expose cleanup() + Gracefully stop any running threads + """ + for section in [ + self.load_basecall_section, + self.load_signal_section, + self.predict_section, + self.landscape_section, + ]: + try: + if hasattr(section, "cleanup"): + section.cleanup() + except Exception: + pass + + # Stop any still-running threads + for t in list(self._active_threads): + try: + if t.isRunning(): + t.quit() + t.wait(2000) + except Exception: + pass + self._active_threads.clear() + self._active_workers.clear() + + super().closeEvent(event) + + def main(): - print("Launching GUI...") - app = QApplication(sys.argv) - window = MainApp() - window.show() - sys.exit(app.exec()) + print("Launching GUI...") + app = QApplication(sys.argv) + window = MainApp() + window.show() + sys.exit(app.exec()) + if __name__ == "__main__": main() diff --git a/DashML/Landscape/Cluster/Centroid_Analysis.py b/DashML/Landscape/Cluster/Centroid_Analysis.py index bfaf7aef..d0a0e3b6 100644 --- a/DashML/Landscape/Cluster/Centroid_Analysis.py +++ b/DashML/Landscape/Cluster/Centroid_Analysis.py @@ -2,9 +2,7 @@ import platform import numpy as np import pandas as pd -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt + import seaborn as sns from scipy.cluster.hierarchy import dendrogram, linkage from sklearn.cluster import AgglomerativeClustering, KMeans @@ -13,57 +11,80 @@ from kmodes.kmodes import KModes -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) -if platform.system() == 'Linux': +if platform.system() == "Linux": ##### server ##### data_path = "/home/jwbear/projects/def-jeromew/jwbear/dendrogram/Out/" - save_path = "/home/jwbear/projects/def-jeromew/jwbear/dendrogram/Dendrogram/Dendrogram_Out/" + save_path = ( + "/home/jwbear/projects/def-jeromew/jwbear/dendrogram/Dendrogram/Dendrogram_Out/" + ) else: data_path = sys.path[1] + "/DashML/Deconvolution/Dendrogram/Putative_Structures/" out_path = sys.path[1] + "/DashML/Deconvolution/Out/" save_path = sys.path[1] + "/DashML/Deconvolution/Dendrogram/Putative_Structures/" -def get_mfe_distances(seq='HCV'): +def get_mfe_distances(seq="HCV"): # native conformation mfes from RNAeval - df_native = pd.read_csv(data_path + "native_mfes.txt", names=['sequence', 'mfe'], - dtype={'sequence':str, 'mfe': float}, - skipinitialspace=True) + df_native = pd.read_csv( + data_path + "native_mfes.txt", + names=["sequence", "mfe"], + dtype={"sequence": str, "mfe": float}, + skipinitialspace=True, + ) # putative conformations mfes from RNAeval, multiple can take average? - df_putative = pd.read_csv(data_path + "putative_mfes.txt", names=['sequence', 'cluster_type', 'cluster_number', 'mfe'], - dtype={'sequence':str, 'cluster_type': str, 'cluster_number':int, 'mfe': float}, - skipinitialspace=True) - #cluster sizes - df_counts = pd.read_csv("/Users/timshel/structure_landscapes/DashML/Deconvolution/Dendrogram/Clusters/" - "cluster_counts.csv", names=['sequence', 'cluster_type', 'cluster_number', 'cluster_size'], - dtype={'sequence': str, 'cluster_type': str, 'cluster_number': int, 'cluster_size': int}, - skipinitialspace=True) + df_putative = pd.read_csv( + data_path + "putative_mfes.txt", + names=["sequence", "cluster_type", "cluster_number", "mfe"], + dtype={ + "sequence": str, + "cluster_type": str, + "cluster_number": int, + "mfe": float, + }, + skipinitialspace=True, + ) + # cluster sizes + df_counts = pd.read_csv( + "/Users/timshel/structure_landscapes/DashML/Deconvolution/Dendrogram/Clusters/" + "cluster_counts.csv", + names=["sequence", "cluster_type", "cluster_number", "cluster_size"], + dtype={ + "sequence": str, + "cluster_type": str, + "cluster_number": int, + "cluster_size": int, + }, + skipinitialspace=True, + ) - #merge data - df_counts['mfe'] = 0 - df_putative['cluster_size'] = 0 - df = df_putative.merge(df_counts, on=['sequence', 'cluster_type', 'cluster_number'], how='left') - df.drop(columns=['mfe_y', 'cluster_size_x'], inplace=True) - df.rename(columns={'mfe_x':'mfe', 'cluster_size_y': 'cluster_size'}, inplace=True) - df = df.sort_values(by=['sequence', 'cluster_type', 'cluster_number']) + # merge data + df_counts["mfe"] = 0 + df_putative["cluster_size"] = 0 + df = df_putative.merge( + df_counts, on=["sequence", "cluster_type", "cluster_number"], how="left" + ) + df.drop(columns=["mfe_y", "cluster_size_x"], inplace=True) + df.rename(columns={"mfe_x": "mfe", "cluster_size_y": "cluster_size"}, inplace=True) + df = df.sort_values(by=["sequence", "cluster_type", "cluster_number"]) # get delta(native mfe, cluster mfe) # todo add seq variable loop - native_mfe = df_native.loc[df_native['sequence']==seq, 'mfe'][0] - df['native_mfe'] = native_mfe + native_mfe = df_native.loc[df_native["sequence"] == seq, "mfe"][0] + df["native_mfe"] = native_mfe # negative means mfe of cluster is more stable than native conformation # a more negative mfe is more stable # a more positive mfe is less stable # native should be most negative, more negative the better in clusters - df['delta_native'] = np.subtract(np.abs(df['native_mfe']),np.abs(df['mfe'])) + df["delta_native"] = np.subtract(np.abs(df["native_mfe"]), np.abs(df["mfe"])) print(df) -#get_mfe_distances() -#sys.exit(0) +# get_mfe_distances() +# sys.exit(0) -#get_cluster_hamming_distance() -#distance_matrix() +# get_cluster_hamming_distance() +# distance_matrix() diff --git a/DashML/Landscape/Cluster/Centroid_Fold.py b/DashML/Landscape/Cluster/Centroid_Fold.py index e527cfe1..6f7180b8 100644 --- a/DashML/Landscape/Cluster/Centroid_Fold.py +++ b/DashML/Landscape/Cluster/Centroid_Fold.py @@ -12,162 +12,216 @@ # 900 pairings in 30 x 30 x seqlen # used more for interactions to capture hidden pairings, pomdp -#RNA basepairing with reactivity is probably better here +# RNA basepairing with reactivity is probably better here # RNAfold -p -d2 --noLP --MEA --shape=HCV_rnafold2.dat < hcv.fa > hcv_bp.out # RNAcofold -a -d2 --noLP < sequences.fa > cofold.out # bp percentages where predict is true but over 95% can be unmodified, may be due to cluster effects # ignore non-predicted or missing values -#RNAfold -p -d2 --noLP < test_sequenc.fa > test_sequenc.out +# RNAfold -p -d2 --noLP < test_sequenc.fa > test_sequenc.out # RNAcofold -a -d2 --noLP < sequences.fa > cofold.out # $ RNAfold --shape=reactivities.dat < sequence.fa # where the file reactivities.dat is a two column text file with sequence positions (1-based) # normalized reactivity values (usually between 0 and 2. Missing values may be left out, or assigned a negative score: -MAX_THREADS=100 +MAX_THREADS = 100 -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) gen = SnowflakeGenerator(42) + def get_tmp_dir(): - package_dir = os.path.dirname(os.path.abspath(__file__)) # Folder containing this file + package_dir = os.path.dirname( + os.path.abspath(__file__) + ) # Folder containing this file tmp_dir = os.path.join(package_dir, "TMP") os.makedirs(tmp_dir, exist_ok=True) return tmp_dir + +# reactivity format for RNAcofold def scale_reactivities(reactivities): try: - min = reactivities.min() - max = reactivities.max() - smin = 0 - smax = 2 + rmin = float(np.nanmin(reactivities)) + rmax = float(np.nanmax(reactivities)) + smin, smax = 0.0, 2.0 - reactivities = ((reactivities - min) / (max - min)) * (smax - smin) + smin + reactivities = ((reactivities - rmin) / (rmax - rmin)) * (smax - smin) + smin except ValueError: # raised if `y` is empty. pass - return reactivities -#extract putative mfes from output + +# extract putative mfes from output def extract_mfes(output, dt, dt2): - contig1 = re.sub('\[|\]|\\s|\\n+|\|\"|\'', '', str(dt['contig'].unique())) - lid = re.sub('\[|\]|\\s|\\n+|\|\"|\'', '', str(dt['LID'].unique())) - ssid = re.sub('(UUID)|\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(dt['SSID'].unique())) - cluster = int(re.sub('\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(dt['cluster'].unique()))) - method = re.sub('\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(dt['method'].unique())) - - contig2 = re.sub('\[|\]|\\s|\\n+|\|\"|\'', '', str(dt2['contig'].unique())) - lid2 = re.sub('\[|\]|\\s|\\n+|\|\"|\'', '', str(dt2['LID'].unique())) - ssid2 = re.sub('(UUID)|\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(dt2['SSID'].unique())) - cluster2 = int(re.sub('\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(dt2['cluster'].unique()))) - method2 = re.sub('\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(dt2['method'].unique())) - - df = pd.DataFrame(columns=['SSID', 'LID1', 'contig1', 'LID2', 'contig2', 'cluster1', 'cluster2', 'secondary', 'mfe', - 'frequency', 'deltag', 'type', 'mea', 'method']) + contig1 = re.sub("\[|\]|\\s|\\n+|\|\"|'", "", str(dt["contig"].unique())) + lid = re.sub("\[|\]|\\s|\\n+|\|\"|'", "", str(dt["LID"].unique())) + ssid = re.sub("(UUID)|\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(dt["SSID"].unique())) + cluster = int( + re.sub("\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(dt["cluster"].unique())) + ) + method = re.sub("\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(dt["method"].unique())) + + contig2 = re.sub("\[|\]|\\s|\\n+|\|\"|'", "", str(dt2["contig"].unique())) + lid2 = re.sub("\[|\]|\\s|\\n+|\|\"|'", "", str(dt2["LID"].unique())) + ssid2 = re.sub("(UUID)|\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(dt2["SSID"].unique())) + cluster2 = int( + re.sub("\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(dt2["cluster"].unique())) + ) + method2 = re.sub("\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(dt2["method"].unique())) + + df = pd.DataFrame( + columns=[ + "SSID", + "LID1", + "contig1", + "LID2", + "contig2", + "cluster1", + "cluster2", + "secondary", + "mfe", + "frequency", + "deltag", + "type", + "mea", + "method", + ] + ) # Define the conversion dictionary - convert_dict = {'SSID': int, 'LID1': int, 'contig1': str, 'LID2': int, 'contig2': str, 'cluster1': int, 'cluster2': int, - 'secondary': str, 'mfe' : float,'frequency' : float, 'deltag' : float, - 'type': str, 'mea': float, 'method': str} + convert_dict = { + "SSID": int, + "LID1": int, + "contig1": str, + "LID2": int, + "contig2": str, + "cluster1": int, + "cluster2": int, + "secondary": str, + "mfe": float, + "frequency": float, + "deltag": float, + "type": str, + "mea": float, + "method": str, + } # Convert columns using the dictionary df = df.astype(convert_dict) mfe, freq, deltag, mea = 0, 0, 0, 0 - ab, aa, bb, a, b = "", "","","","" + ab, aa, bb, a, b = "", "", "", "", "" for i, line in enumerate(output.split("\\n")): - #print(line) - secondary, type = "", "" - if re.search('(frequency)', line): - l = re.split(";", line,1) - #print(l) - f = re.search('([0-9]+.[0-9]+e*-*[0-9]*)', l[0]) + # print(line) + secondary, type = "", "" + if re.search("(frequency)", line): + l = re.split(";", line, 1) + # print(l) + f = re.search("([0-9]+.[0-9]+e*-*[0-9]*)", l[0]) freq = f.group() - #print("frequency", freq) - d = re.search('-*([0-9]+\.[0-9]+e*-*[0-9]*)', l[1]) + # print("frequency", freq) + d = re.search("-*([0-9]+\.[0-9]+e*-*[0-9]*)", l[1]) deltag = d.group() - #print("deltag", deltag) - elif re.search(r'(-*[0-9]+\.[0-9]+\\t)', line): - l = re.split(r'\\t', line) + # print("deltag", deltag) + elif re.search(r"(-*[0-9]+\.[0-9]+\\t)", line): + l = re.split(r"\\t", line) ab, aa, bb, a, b = l[0], l[1], l[2], l[3], l[4] - #print("energies", ab, aa, bb, a, b) - elif re.search('^((\.+)|(\(+|\)+))\W', line): - l = re.split("\s", line,1) + # print("energies", ab, aa, bb, a, b) + elif re.search("^((\.+)|(\(+|\)+))\W", line): + l = re.split("\s", line, 1) secondary = l[0] - #print("structure", secondary) + # print("structure", secondary) energy = l[1] - if re.search('^(\(|-[0-9])', energy): + if re.search("^(\(|-[0-9])", energy): type = "MFE" - mfe = re.sub('\(|\)', '', energy) - #print("mfe", mfe) - elif re.search('^(\[|-[0-9])', energy): + mfe = re.sub("\(|\)", "", energy) + # print("mfe", mfe) + elif re.search("^(\[|-[0-9])", energy): type = "PROB" - mfe = re.sub('\[|\]','', energy) - #print("prob mfe", mfe) - elif re.search('^(\{|-[0-9])',energy): + mfe = re.sub("\[|\]", "", energy) + # print("prob mfe", mfe) + elif re.search("^(\{|-[0-9])", energy): i = 0 e = re.split("\s", energy) if len(e[0]) <= 1: i = 1 - if re.search('(d\=)', energy): + if re.search("(d\=)", energy): type = "CENTROID" - mfe = re.sub('\{|\}|d|\=', '', e[i]) - #print("centroid mfe",mfe) - distance = re.sub('\{|\}|d|\=', '', e[i+1]) - #print("centroid distance", distance) - elif re.search('(MEA\=)', energy): + mfe = re.sub("\{|\}|d|\=", "", e[i]) + # print("centroid mfe",mfe) + distance = re.sub("\{|\}|d|\=", "", e[i + 1]) + # print("centroid distance", distance) + elif re.search("(MEA\=)", energy): type = "MEA" - mfe = re.sub('\{|\}|(MEA)|\=', '', e[i]) - #print("mea mfe", mfe) - mea = re.sub('\{|\}|(MEA)|\=', '', e[i + 1]) - #print("MEA", mea) + mfe = re.sub("\{|\}|(MEA)|\=", "", e[i]) + # print("mea mfe", mfe) + mea = re.sub("\{|\}|(MEA)|\=", "", e[i + 1]) + # print("MEA", mea) #'LID', 'contig', 'secondary', 'mfe', 'cluster', 'frequency', 'diversity', 'type', 'distance' if secondary != "": if len(mfe) == 0: mfe = 0 - df.loc[len(df)] = [ssid, lid, contig1, lid2, contig2, cluster, cluster2, secondary, mfe, freq, deltag, type, mea] - - df.loc[df['type'] == 'MFE', 'frequency'] = float(freq) - df.loc[df['type'] == 'MFE', 'deltag'] = float(deltag) - #print(df) - - #insert into db + df.loc[len(df)] = [ + ssid, + lid, + contig1, + lid2, + contig2, + cluster, + cluster2, + secondary, + mfe, + freq, + deltag, + type, + mea, + ] + + df.loc[df["type"] == "MFE", "frequency"] = float(freq) + df.loc[df["type"] == "MFE", "deltag"] = float(deltag) + # print(df) + + # insert into db dbins.insert_centroid_secondary_intrx(df) - return {'AB': ab, 'AA': aa, 'BB' : bb, 'A': a, 'B': b} + return {"AB": ab, "AA": aa, "BB": bb, "A": a, "B": b} -def get_bpfiles(contig, contig2, cluster, cluster2,method): + +def get_bpfiles(contig, contig2, cluster, cluster2, method): file_list = { - 'AA' : "AA{0},{1}{2}{3}{4}_dp5.ps".format(contig, contig2, cluster, cluster2), - 'AB' :"AB{0},{1}{2}{3}{4}_dp5.ps".format(contig, contig2, cluster, cluster2), - 'BB' : "BB{0},{1}{2}{3}{4}_dp5.ps".format(contig, contig2, cluster, cluster2), - 'A' : "A{0},{1}{2}{3}{4}_dp5.ps".format(contig, contig2, cluster, cluster2), - 'B' : "B{0},{1}{2}{3}{4}_dp5.ps".format(contig, contig2, cluster, cluster2) + "AA": "AA{0},{1}{2}{3}{4}_dp5.ps".format(contig, contig2, cluster, cluster2), + "AB": "AB{0},{1}{2}{3}{4}_dp5.ps".format(contig, contig2, cluster, cluster2), + "BB": "BB{0},{1}{2}{3}{4}_dp5.ps".format(contig, contig2, cluster, cluster2), + "A": "A{0},{1}{2}{3}{4}_dp5.ps".format(contig, contig2, cluster, cluster2), + "B": "B{0},{1}{2}{3}{4}_dp5.ps".format(contig, contig2, cluster, cluster2), } return file_list -#extract base pairing probabilities from output save to Structure_BPP + +# extract base pairing probabilities from output save to Structure_BPP # ubox The upper right triangle displays all predicted base pairs with not more than two inconsistent sequences. # lbox The lower left triangle contains only the secondary structure formed by the most believable base pairs. # It still contains pairs that are not in the final prediction, lbox all equal def extract_bpps(df, df2, mfes): - contig = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df['contig'].unique())) - ssid = re.sub('(UUID)|\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(df['SSID'].unique())) - seqlen = len(df2['sequence'].unique()) - cluster = re.sub('\[|\]|\\s|\\n+|\|\"|\'', '', str(df['cluster'].unique())) - method = re.sub('\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(df['method'].unique())) + contig = re.sub("\[|\]|\\s|\\n+|'|\"", "", str(df["contig"].unique())) + ssid = re.sub("(UUID)|\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(df["SSID"].unique())) + seqlen = len(df2["sequence"].unique()) + cluster = re.sub("\[|\]|\\s|\\n+|\|\"|'", "", str(df["cluster"].unique())) + method = re.sub("\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(df["method"].unique())) - contig2 = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df2['contig'].unique())) - seqlen2 = len(df2['sequence'].unique()) - cluster2 = re.sub('\[|\]|\\s|\\n+|\|\"|\'', '', str(df2['cluster'].unique())) - method2 = re.sub('\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(df2['method'].unique())) + contig2 = re.sub("\[|\]|\\s|\\n+|'|\"", "", str(df2["contig"].unique())) + seqlen2 = len(df2["sequence"].unique()) + cluster2 = re.sub("\[|\]|\\s|\\n+|\|\"|'", "", str(df2["cluster"].unique())) + method2 = re.sub("\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(df2["method"].unique())) # get file list by creating file names bp_files = get_bpfiles(contig, contig2, cluster, cluster2, method) # for each file loop through and extract bpp and probabilities for rx, bpf in bp_files.items(): - #print(rx) + # print(rx) mfe = mfes[rx] base1, base2, prob = [], [], [] + def get_bases(): try: tmp_dir = get_tmp_dir() @@ -176,7 +230,9 @@ def get_bases(): with open(tmp_file_path, "r+") as f: lines = f.readlines() for i, line in enumerate(lines): - if re.search("([0-9]+ [0-9]+ [0-9]+\.*[0-9]* (ubox))", line.strip()): + if re.search( + "([0-9]+ [0-9]+ [0-9]+\.*[0-9]* (ubox))", line.strip() + ): for l in lines[i:]: if re.search("(showpage)|(end)|(%%EOF)", l): return @@ -185,7 +241,7 @@ def get_bases(): base1.append(int(bases[0])) base2.append(int(bases[1])) prob.append(float(bases[2])) - #print(bases) + # print(bases) except Exception as e: raise Exception(str(e)) print(e) @@ -195,15 +251,31 @@ def get_bases(): get_bases() bppid = next(gen) - df_bpp = pd.DataFrame({'BPPID': bppid, 'SSID': ssid, 'type_interaction': rx, 'mfe':float(mfe)}, index=[0]) - df_prob = pd.DataFrame({'BPPID': bppid, 'SSID': ssid, 'base1':base1, 'base2':base2, 'probability':prob}) + df_bpp = pd.DataFrame( + {"BPPID": bppid, "SSID": ssid, "type_interaction": rx, "mfe": float(mfe)}, + index=[0], + ) + df_prob = pd.DataFrame( + { + "BPPID": bppid, + "SSID": ssid, + "base1": base1, + "base2": base2, + "probability": prob, + } + ) # fix bases RNAcofold extends base numbers for different molecules, AA - df_prob['base1'] = np.where(df_prob['base1'] > seqlen, df_prob['base1'] - seqlen, df_prob['base1']) - df_prob['base2'] = np.where(df_prob['base2'] > seqlen2, df_prob['base2'] - seqlen2, df_prob['base2']) - #insert into db + df_prob["base1"] = np.where( + df_prob["base1"] > seqlen, df_prob["base1"] - seqlen, df_prob["base1"] + ) + df_prob["base2"] = np.where( + df_prob["base2"] > seqlen2, df_prob["base2"] - seqlen2, df_prob["base2"] + ) + # insert into db dbins.insert_centroid_bpp(df_bpp, df_prob) return + def get_reactive_dimers(df, df2): try: print(df.head()) @@ -212,34 +284,40 @@ def get_reactive_dimers(df, df2): ## vienna rna does not always follow naming conventions if names exceed length tmp = get_tmp_dir() - contig = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df['contig'].unique())) - cluster = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df['cluster'].unique())) - sequence = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df['sequence'].unique())) - method = re.sub('\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(df['method'].unique())) - reactivity = scale_reactivities(df['centroid'].to_numpy()) + contig = re.sub("\[|\]|\\s|\\n+|'|\"", "", str(df["contig"].unique())) + cluster = re.sub("\[|\]|\\s|\\n+|'|\"", "", str(df["cluster"].unique())) + sequence = re.sub("\[|\]|\\s|\\n+|'|\"", "", str(df["sequence"].unique())) + method = re.sub("\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(df["method"].unique())) + reactivity = scale_reactivities(df["centroid"].to_numpy()) - contig2 = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df2['contig'].unique())) - cluster2 = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df2['cluster'].unique())) - sequence2 = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df2['sequence'].unique())) - method2 = re.sub('\(|\)|\[|\]|\\s|\\n+|\|\"|\'', '', str(df2['method'].unique())) - reactivity2 = scale_reactivities(df2['centroid'].to_numpy()) + contig2 = re.sub("\[|\]|\\s|\\n+|'|\"", "", str(df2["contig"].unique())) + cluster2 = re.sub("\[|\]|\\s|\\n+|'|\"", "", str(df2["cluster"].unique())) + sequence2 = re.sub("\[|\]|\\s|\\n+|'|\"", "", str(df2["sequence"].unique())) + method2 = re.sub("\(|\)|\[|\]|\\s|\\n+|\|\"|'", "", str(df2["method"].unique())) + reactivity2 = scale_reactivities(df2["centroid"].to_numpy()) if len(sequence + sequence2) >= 5000: - sys.error('Length Exceeds Vienna CoFold') - #TODO pesky switch to RNAFold only to calculate putative structure, no base pairing info + sys.error("Length Exceeds Vienna CoFold") + # TODO pesky switch to RNAFold only to calculate putative structure, no base pairing info # limit is 5k nt - #RNAfold -p -d2 --noLP < sequence1.fa > sequence1.out + # RNAfold -p -d2 --noLP < sequence1.fa > sequence1.out return None # create input file - sf = "{}{}{}{}{}{}_sequence.fa".format(tmp, contig, contig2, cluster, cluster2, method) + sf = "{}{}{}{}{}{}_sequence.fa".format( + tmp, contig, contig2, cluster, cluster2, method + ) f = open(sf, "w") - lines = "> {},{}{}{}{}\n{}&{}\n".format(contig, contig2, cluster, cluster2, sequence, sequence2, method) + lines = "> {},{}{}{}{}\n{}&{}\n".format( + contig, contig2, cluster, cluster2, sequence, sequence2, method + ) f.write(lines) f.close() # create dat file of reactivities - datf = "{}{}{}{}{}{}.dat".format(tmp, contig, contig2, cluster, cluster2, method) + datf = "{}{}{}{}{}{}.dat".format( + tmp, contig, contig2, cluster, cluster2, method + ) f = open(datf, "w") i = 0 for i in range(0, len(reactivity)): @@ -254,14 +332,18 @@ def get_reactive_dimers(df, df2): in_path = sf # print(in_path) # out_path = save_path_dmso + seqname + "_" + seqname + "_intra.out" - p1 = subprocess.Popen(["RNAcofold", "-p", "-a", "-d2", "--noLP", "--MEA", dat_file, in_path], - stdin=subprocess.PIPE, stdout=subprocess.PIPE, cwd=tmp) + p1 = subprocess.Popen( + ["RNAcofold", "-p", "-a", "-d2", "--noLP", "--MEA", dat_file, in_path], + stdin=subprocess.PIPE, + stdout=subprocess.PIPE, + cwd=tmp, + ) output = str(p1.communicate(timeout=None)) - #print(output) + # print(output) # generate unique run number ssid = next(gen) - df['SSID'] = ssid - df2['SSID'] = ssid + df["SSID"] = ssid + df2["SSID"] = ssid # extract mfes and base pairing data mfes = extract_mfes(output, df, df2) # mfes = {'AB': .5, 'AA': .7, 'BB': .8, 'A': .9, 'B': .11} @@ -276,9 +358,9 @@ def get_reactive_dimers(df, df2): raise Exception(str(e)) return None finally: - if 'f' in locals(): + if "f" in locals(): f.close() - if 'p1' in locals(): + if "p1" in locals(): p1.stdout.close() p1.stdin.close() @@ -290,21 +372,21 @@ def get_probabilities(lids): def get_probpool(): with ThreadPoolExecutor(max_workers=MAX_THREADS) as executor: # intra & inter for interactions - for li in df['LID'].unique(): - d1 = df.loc[(df['LID'] == li) & (df['method'] == 'kmeans')] - for li2 in df['LID'].unique(): - d2 = df.loc[(df['LID'] == li2) & (df['method'] == 'kmeans')] - for cluster in d1['cluster'].unique(): - dc1 = d1.loc[(d1['cluster'] == cluster)] - for cluster2 in d2['cluster'].unique(): - dc2 = d2.loc[(d2['cluster'] == cluster2)] - if len(dc1)>0 and len(dc2)>0: + for li in df["LID"].unique(): + d1 = df.loc[(df["LID"] == li) & (df["method"] == "kmeans")] + for li2 in df["LID"].unique(): + d2 = df.loc[(df["LID"] == li2) & (df["method"] == "kmeans")] + for cluster in d1["cluster"].unique(): + dc1 = d1.loc[(d1["cluster"] == cluster)] + for cluster2 in d2["cluster"].unique(): + dc2 = d2.loc[(d2["cluster"] == cluster2)] + if len(dc1) > 0 and len(dc2) > 0: executor.submit(get_reactive_dimers, dc1, dc2) - get_probpool() return -#get_probabilities('36') -#sys.exit(0) + +# get_probabilities('36') +# sys.exit(0) diff --git a/DashML/Landscape/Cluster/Centroid_MFE.py b/DashML/Landscape/Cluster/Centroid_MFE.py index 2f1c8b29..a8bbd919 100644 --- a/DashML/Landscape/Cluster/Centroid_MFE.py +++ b/DashML/Landscape/Cluster/Centroid_MFE.py @@ -13,101 +13,145 @@ # Calculates putative structure and MFE of cluster centroids using reactivities # Used to measure distance in MFE from native structure -#RNA basepairing with reactivity is probably better here +# RNA basepairing with reactivity is probably better here # RNAfold -p -d2 --noLP --MEA --shape=HCV_rnafold2.dat < hcv.fa > hcv_bp.out # RNAcofold -a -d2 --noLP < sequences.fa > cofold.out # todo bp percentages where predict is true but over 95% can be unmodified # ignore non-predicted or missing values -#RNAfold -p -d2 --noLP < test_sequenc.fa > test_sequenc.out +# RNAfold -p -d2 --noLP < test_sequenc.fa > test_sequenc.out # RNAcofold -a -d2 --noLP < sequences.fa > cofold.out # $ RNAfold --shape=reactivities.dat < sequence.fa # where the file reactivities.dat is a two column text file with sequence positions (1-based) # normalized reactivity values (usually between 0 and 2. Missing values may be left out, or assigned a negative score: -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) + def get_tmp_dir(): - package_dir = os.path.dirname(os.path.abspath(__file__)) # Folder containing this file + package_dir = os.path.dirname( + os.path.abspath(__file__) + ) # Folder containing this file tmp_dir = os.path.join(package_dir, "TMP") os.makedirs(tmp_dir, exist_ok=True) return tmp_dir + +# reactivity format for RNAcofold def scale_reactivities(reactivities): - min = reactivities.min() - max = reactivities.max() - smin = 0 - smax = 2 + try: + rmin = float(np.nanmin(reactivities)) + rmax = float(np.nanmax(reactivities)) + smin, smax = 0.0, 2.0 - reactivities = ((reactivities - min) / (max - min)) * (smax - smin) + smin + reactivities = ((reactivities - rmin) / (rmax - rmin)) * (smax - smin) + smin + except ValueError: # raised if `y` is empty. + pass return reactivities -#extract putative mfes from output -def extract_mfes(lid, seq, cluster, clust_type,output): - df = pd.DataFrame(columns=['LID', 'contig', 'secondary', 'mfe', 'cluster', - 'frequency', 'diversity', 'type', 'distance', 'mea', 'method']) + +# extract putative mfes from output +def extract_mfes(lid, seq, cluster, clust_type, output): + df = pd.DataFrame( + columns=[ + "LID", + "contig", + "secondary", + "mfe", + "cluster", + "frequency", + "diversity", + "type", + "distance", + "mea", + "method", + ] + ) # Define the conversion dictionary - convert_dict = {'LID': int, 'contig': str, 'secondary': str, 'mfe' : float, 'cluster' : int, - 'frequency' : float, 'diversity' : float, 'type' : str, 'distance': float, - 'mea': float, 'method': str} + convert_dict = { + "LID": int, + "contig": str, + "secondary": str, + "mfe": float, + "cluster": int, + "frequency": float, + "diversity": float, + "type": str, + "distance": float, + "mea": float, + "method": str, + } # Convert columns using the dictionary df = df.astype(convert_dict) for i, line in enumerate(output.split("\\n")): - secondary, type = "", "" - mfe, distance, freq, diversity, distance, mea = 0, 0 ,0, 0, 0, 0 - if re.search('(frequency)', line): - l = re.split(";", line,1) - #print(l) - f = re.search('([0-9]+.[0-9]+e*-*[0-9]*)', l[0]) + secondary, type = "", "" + mfe, distance, freq, diversity, distance, mea = 0, 0, 0, 0, 0, 0 + if re.search("(frequency)", line): + l = re.split(";", line, 1) + # print(l) + f = re.search("([0-9]+.[0-9]+e*-*[0-9]*)", l[0]) freq = f.group() - #print("frequency", freq) - d = re.search('([0-9]+.[0-9]+e*-*[0-9]*)', l[1]) + # print("frequency", freq) + d = re.search("([0-9]+.[0-9]+e*-*[0-9]*)", l[1]) diversity = d.group() - #print("diversity", diversity) - df.loc[df['type']=='MFE', 'frequency'] = float(freq) - df.loc[df['type'] == 'MFE', 'diversity'] = float(diversity) - elif re.search('(\.+)|(\(+|\)+)\W', line): - l = re.split("\s", line,1) + # print("diversity", diversity) + df.loc[df["type"] == "MFE", "frequency"] = float(freq) + df.loc[df["type"] == "MFE", "diversity"] = float(diversity) + elif re.search("(\.+)|(\(+|\)+)\W", line): + l = re.split("\s", line, 1) secondary = l[0] - #print("structure", secondary) + # print("structure", secondary) energy = l[1] - if re.search('^(\(|-[0-9])', energy): + if re.search("^(\(|-[0-9])", energy): type = "MFE" - mfe = re.sub('\(|\)', '', energy) - #print("mfe", mfe) - elif re.search('^(\[|-[0-9])', energy): + mfe = re.sub("\(|\)", "", energy) + # print("mfe", mfe) + elif re.search("^(\[|-[0-9])", energy): type = "PROB" - mfe = re.sub('\[|\]','', energy) - #print("prob mfe", mfe) - elif re.search('^(\{|-[0-9])',energy): + mfe = re.sub("\[|\]", "", energy) + # print("prob mfe", mfe) + elif re.search("^(\{|-[0-9])", energy): i = 0 e = re.split("\s", energy) if len(e[0]) <= 1: i = 1 - if re.search('(d\=)', energy): + if re.search("(d\=)", energy): type = "CENTROID" - mfe = re.sub('\{|\}|d|\=', '', e[i]) - #print("centroid mfe",mfe) - distance = re.sub('\{|\}|d|\=', '', e[i+1]) - #print("centroid distance", distance) - elif re.search('(MEA\=)', energy): + mfe = re.sub("\{|\}|d|\=", "", e[i]) + # print("centroid mfe",mfe) + distance = re.sub("\{|\}|d|\=", "", e[i + 1]) + # print("centroid distance", distance) + elif re.search("(MEA\=)", energy): type = "MEA" - mfe = re.sub('\{|\}|(MEA)|\=', '', e[i]) - #print("mea mfe", mfe) - mea = re.sub('\{|\}|(MEA)|\=', '', e[i + 1]) - #print("MEA", mea) + mfe = re.sub("\{|\}|(MEA)|\=", "", e[i]) + # print("mea mfe", mfe) + mea = re.sub("\{|\}|(MEA)|\=", "", e[i + 1]) + # print("MEA", mea) #'LID', 'contig', 'secondary', 'mfe', 'cluster', 'frequency', 'diversity', 'type', 'distance' if secondary != "": if len(mfe) == 0: mfe = 0 - df.loc[len(df)] = [lid, seq, secondary, mfe, cluster, freq, diversity, type, distance, mea, clust_type] + df.loc[len(df)] = [ + lid, + seq, + secondary, + mfe, + cluster, + freq, + diversity, + type, + distance, + mea, + clust_type, + ] - #insert into db + # insert into db dbins.insert_centroid_secondary(df) return + ##### get putative structures for cluster centroids ##### # RNAfold -p -a -d2 --noLP --MEA < sequence.fa > sequence_cofold.out def getRNAfold(lid, seqname, sequence, clust_num, clust_rx, clust_type): @@ -122,7 +166,7 @@ def getRNAfold(lid, seqname, sequence, clust_num, clust_rx, clust_type): # Compose full file path tmp_file_path = os.path.join(tmp_dir, "sequence.fa") - #create input file + # create input file f = open(tmp_file_path, "w") f.write(">" + seqname + "-" + str(clust_num) + "\n" + sequence + "\n") f.close() @@ -135,18 +179,31 @@ def getRNAfold(lid, seqname, sequence, clust_num, clust_rx, clust_type): f.write(str(i + 1) + "\t" + str(clust_rx[i]) + "\n") f.close() dat_file = "--shape=" + tmp_file_path - #print("dat file: ", dat_file) - - #send to rnacofold - #in_path = save_path_dir + "/sequence.fa" - #print(in_path) - #out_path = save_path_dir + '/' + seqname + "_" + str(clust_num) + ".out" - p1 = subprocess.Popen(["RNAfold", "-p","-S 1.2", "-d2", "--noLP", "--MEA", dat_file, "sequence.fa"], - stdin=subprocess.PIPE, stdout=subprocess.PIPE, cwd=get_tmp_dir()) - #reval = sequence + "\n" + ss + "\n@" - #p1.communicate(input=reval.encode()) + # print("dat file: ", dat_file) + + # send to rnacofold + # in_path = save_path_dir + "/sequence.fa" + # print(in_path) + # out_path = save_path_dir + '/' + seqname + "_" + str(clust_num) + ".out" + p1 = subprocess.Popen( + [ + "RNAfold", + "-p", + "-S 1.2", + "-d2", + "--noLP", + "--MEA", + dat_file, + "sequence.fa", + ], + stdin=subprocess.PIPE, + stdout=subprocess.PIPE, + cwd=get_tmp_dir(), + ) + # reval = sequence + "\n" + ss + "\n@" + # p1.communicate(input=reval.encode()) output = str(p1.communicate(timeout=300)) - extract_mfes(lid, seqname, clust_num, clust_type,output) + extract_mfes(lid, seqname, clust_num, clust_type, output) except subprocess.CalledProcessError as e: print(e) return None @@ -155,7 +212,7 @@ def getRNAfold(lid, seqname, sequence, clust_num, clust_rx, clust_type): raise Exception(str(e)) return None finally: - #f.close() + # f.close() p1.stdout.close() p1.stdin.close() @@ -165,29 +222,38 @@ def getRNAfold(lid, seqname, sequence, clust_num, clust_rx, clust_type): def get_putative_structure(lids): # dataframes of centroids df = dbsel.select_centroidz(lids) - clusters = df['cluster'].unique() + clusters = df["cluster"].unique() with ThreadPoolExecutor(max_workers=MAX_THREADS) as executor: for cluster in clusters: - dt = df.loc[(df['cluster']==cluster) & (df['method']=='kmeans')] - for li in dt['LID'].unique(): - seq_name = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df['contig'].unique())) - sequence = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df['sequence'].unique())) - reactivities = scale_reactivities(dt['centroid'].to_numpy()) - executor.submit(getRNAfold, li, seq_name, sequence, cluster, reactivities, 'kmeans') + dt = df.loc[(df["cluster"] == cluster) & (df["method"] == "kmeans")] + for li in dt["LID"].unique(): + seq_name = re.sub("\[|\]|\\s|\\n+|'|\"", "", str(df["contig"].unique())) + sequence = re.sub( + "\[|\]|\\s|\\n+|'|\"", "", str(df["sequence"].unique()) + ) + reactivities = scale_reactivities(dt["centroid"].to_numpy()) + executor.submit( + getRNAfold, li, seq_name, sequence, cluster, reactivities, "kmeans" + ) for cluster in clusters: - dt = df.loc[(df['cluster']==cluster) & (df['method']=='hamming')] - for li in dt['LID'].unique(): - seq_name = re.sub('\[|\]|\\s|\\n+|\|\"|\'', '', str(df['contig'].unique())) - sequence = re.sub('\[|\]|\\s|\\n+|\'|\"', '', str(df['sequence'].unique())) - dt.loc[dt['centroid'] == 1, 'centorid'] = 0 - dt.loc[dt['centroid'] == -1, 'centroid'] = 2 - reactivities = dt['centroid'].to_numpy() - executor.submit(getRNAfold,li, seq_name, sequence, cluster, reactivities, 'hamming') - + dt = df.loc[(df["cluster"] == cluster) & (df["method"] == "hamming")] + for li in dt["LID"].unique(): + seq_name = re.sub( + "\[|\]|\\s|\\n+|\|\"|'", "", str(df["contig"].unique()) + ) + sequence = re.sub( + "\[|\]|\\s|\\n+|'|\"", "", str(df["sequence"].unique()) + ) + dt.loc[dt["centroid"] == 1, "centorid"] = 0 + dt.loc[dt["centroid"] == -1, "centroid"] = 2 + reactivities = dt["centroid"].to_numpy() + executor.submit( + getRNAfold, li, seq_name, sequence, cluster, reactivities, "hamming" + ) #### TODO check reactivities (kmeans) from centroid files -#get_putative_structure(37) -#sys.exit(0) +# get_putative_structure(37) +# sys.exit(0) diff --git a/DashML/Landscape/Cluster/Cluster_Num.py b/DashML/Landscape/Cluster/Cluster_Num.py index 9b4a73cf..7c9fdef0 100644 --- a/DashML/Landscape/Cluster/Cluster_Num.py +++ b/DashML/Landscape/Cluster/Cluster_Num.py @@ -4,10 +4,12 @@ import math import numpy as np import pandas as pd + +from pathlib import Path import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt -import seaborn as sns +from matplotlib.figure import Figure +from matplotlib.backends.backend_agg import FigureCanvasAgg + from scipy.cluster.hierarchy import dendrogram, linkage, cophenet, fcluster from sklearn.cluster import AgglomerativeClustering, KMeans from sklearn.metrics import silhouette_score @@ -18,163 +20,252 @@ #### Get Average Number of Optimal Clusters for Predictions #### -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) -save_path=None +save_path: Path | None = None -#reactivity format for RNAcofold + +# reactivity format for RNAcofold def scale_reactivities(reactivities): - min = reactivities.min() - max = reactivities.max() - smin = 0 - smax = 2 + try: + rmin = float(np.nanmin(reactivities)) + rmax = float(np.nanmax(reactivities)) + smin, smax = 0.0, 2.0 - reactivities = ((reactivities - min) / (max - min)) * (smax - smin) + smin + reactivities = ((reactivities - rmin) / (rmax - rmin)) * (smax - smin) + smin + except ValueError: # raised if `y` is empty. + pass return reactivities + # Create linkage matrix and then plot the dendrogram -def distancematrix_modes(df, seq='hcv'): +def distancematrix_modes(df, seq="hcv"): ##### Compute Distance Between Reads ##### # reactivity kmode hamming - dx = df.sort_values(by=['read_id', 'position']) - read_list = dx['read_id'].unique() - dx = (dx.pivot(index=["read_id"], columns="position", values="Predict") - .rename_axis(columns=None) - .reset_index()) + dx = df.sort_values(by=["read_id", "position"]) + read_list = dx["read_id"].unique() + dx = ( + dx.pivot(index=["read_id"], columns="position", values="Predict") + .rename_axis(columns=None) + .reset_index() + ) dx.fillna(value=0, inplace=True) - dx = dx.drop(columns=['read_id']) + dx = dx.drop(columns=["read_id"]) mm = dx.to_numpy() - #clustermap - ham = distance.cdist(mm, mm, 'hamming') - #ham = np.divide(ham,np.linalg.norm(ham)) - #print(ham) + # clustermap + ham = distance.cdist(mm, mm, "hamming") + # ham = np.divide(ham,np.linalg.norm(ham)) + # print(ham) return ham, mm, read_list -def distancematrix_means(df, seq='hcv'): + +def distancematrix_means(df, seq="hcv"): ##### Compute Distance Between Reads ##### - #reactivity kmeans distance - dx = df.sort_values(by=['read_id', 'position']) - read_list = dx['read_id'].unique() - dx = (dx.pivot(index=["read_id"], columns="position", values="Reactivity") - .rename_axis(columns=None) - .reset_index()) + # reactivity kmeans distance + dx = df.sort_values(by=["read_id", "position"]) + read_list = dx["read_id"].unique() + dx = ( + dx.pivot(index=["read_id"], columns="position", values="Reactivity") + .rename_axis(columns=None) + .reset_index() + ) dx.fillna(value=0, inplace=True) - dx = dx.drop(columns=['read_id']) + dx = dx.drop(columns=["read_id"]) ## create matrix for euclidean distance reads = dx.to_numpy() - #mm = np.divide(mm, np.linalg.norm(mm)) + # mm = np.divide(mm, np.linalg.norm(mm)) ## calculate distances - distk = distance.cdist(reads, reads, 'euclidean') + distk = distance.cdist(reads, reads, "euclidean") return distk, reads, read_list -def optimal_dendrogram(X, seq, tit, linkage_method='complete', threshold='auto', percent=.2): + +def optimal_dendrogram( + X, seq, tit, linkage_method="complete", threshold="auto", percent=0.2 +): distance_matrix = pdist(X) Z = linkage(distance_matrix, method=linkage_method) - plt.figure(figsize=(8, 5)) - dendrogram(Z) + with matplotlib.rc_context({"font.size": 9}): + fig = Figure(figsize=(8, 5), dpi=100) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + dendrogram(Z, ax=ax) - if 'kmeans' in tit: - percent = .45 + # special percent tweak for kmeans titles + if "kmeans" in (tit or "").lower(): + percent = 0.45 + + if threshold == "auto": + max_height = float(np.max(Z[:, 2])) + threshold_val = percent * max_height + num_clusters = int(np.sum(Z[:, 2] > threshold_val) + 1) + ax.axhline( + y=threshold_val, + color="r", + linestyle="--", + label=f"{int(percent * 100)}% Height = {threshold_val:.2f}", + ) + ax.legend() + print( + f"Estimated number of clusters at {percent * 100}% height threshold: {num_clusters}" + ) + elif isinstance(threshold, (int, float)): + threshold_val = float(threshold) + num_clusters = int(np.sum(Z[:, 2] > threshold_val) + 1) + ax.axhline( + y=threshold_val, + color="r", + linestyle="--", + label=f"Threshold = {threshold_val:g}", + ) + ax.legend() + print( + f"Estimated number of clusters at threshold {threshold_val}: {num_clusters}" + ) + else: + threshold_val = None + num_clusters = None - if threshold == 'auto': - max_height = np.max(Z[:, 2]) - threshold_val = percent * max_height - num_clusters = np.sum(Z[:, 2] > threshold_val) + 1 - plt.axhline(y=threshold_val, color='r', linestyle='--', - label=f'{int(percent * 100)}% Height = {threshold_val:.2f}') - plt.legend() - print(f"Estimated number of clusters at {percent * 100}% height threshold: {num_clusters}") + ax.set_title("Dendrogram") + ax.set_xlabel("Read Index") + ax.set_ylabel("Distance") - elif isinstance(threshold, (int, float)): - num_clusters = np.sum(Z[:, 2] > threshold) + 1 - plt.axhline(y=threshold, color='r', linestyle='--', label=f'Threshold = {threshold}') - plt.legend() - print(f"Estimated number of clusters at threshold {threshold}: {num_clusters}") + assert save_path is not None, "save_path is not set" + img_path = save_path + f"{seq}{tit}_dendrogram.png" + fig.tight_layout() + fig.savefig(img_path, dpi=300, bbox_inches="tight") + fig.clear() + del fig - plt.title('Dendrogram') - plt.xlabel('Read Index') - plt.ylabel('Distance') - plt.savefig(save_path + seq + tit + '_dendrogram.png', dpi=300) - #plt.showblock=False) + return str(img_path), num_clusters, threshold_val def get_optimal_clusters(X, seq, read_list, tit, plot=False): - # Heirarchical Clustering - # perform clustering + """ + Returns (img_sil, img_wss). + """ s_scores = [] c_scores = [] wss_values = [] + n_samples = len(X) - max = min(50, n_samples) - for k in np.arange(2,max, 1): - model = AgglomerativeClustering(n_clusters=k, linkage='complete', metric='precomputed', compute_distances=True, - compute_full_tree=True, distance_threshold=None) - model = model.fit_predict(X) - c_scores.append(1 + np.amax(model)) - #print(np.unique(model)) - #print(f"Number of clusters = {1 + np.amax(model)}") - - ### cluster evaluation - # Silhouette Score maximal score - if (1 + np.amax(model)) >= 2 and (k < len(np.unique(model)) < len(X)): - silhouette = silhouette_score(X, model, metric='hamming') - prev = silhouette - s_scores.append(silhouette) + k_max = min(50, n_samples) + + # Range starts at 2 clusters + for k in range(2, max(3, k_max)): # ensure at least try k=2 + model = AgglomerativeClustering( + n_clusters=k, + linkage="complete", + metric="precomputed", + compute_distances=True, + compute_full_tree=True, + distance_threshold=None, + ) + labels = model.fit_predict(X) + k_found = int(labels.max()) + 1 + c_scores.append(k_found) + + # Silhouette score: + if k_found >= 2 and (k < len(np.unique(labels)) < len(X)): + try: + sil = silhouette_score(X, labels, metric="hamming") + except Exception: + sil = float("nan") else: - s_scores.append(0) + sil = float("nan") + s_scores.append(sil) - # Total Within-Cluster Sum of Squares (WSS) - wss = np.sum([np.sum((X[model == c] - X[model == c].mean(axis=0)) ** 2) for c in np.unique(model)]) + # Total Within-Cluster Sum of Squares (approximate in distance space): + # Move back to a rough feature embedding via classical MDS is ideal, + # but to keep parity we approximate using mean of cluster distances. + wss = 0.0 + for c in np.unique(labels): + idx = labels == c + if idx.sum() > 1: + D = X[np.ix_(idx, idx)] + # sum of squared distances to "centroid" approximated by mean + mu = D.mean(axis=0) + wss += np.sum((D - mu) ** 2) wss_values.append(wss) - #print('Silhouette Score:', silhouette, "WSS:", wss) - print('Silhouette Score:', len(s_scores), "WSS:", len(wss_values)) + print("Silhouette Score points:", len(s_scores), "WSS points:", len(wss_values)) + + img_sil = None + img_wss = None + best_k = None + if plot: - # Plot Elbow Curve - plt.figure(figsize=(6, 4)) - plt.plot( c_scores, s_scores, marker='o', label='Silhouette Score') - plt.title('Elbow Method using Silhouette Score') - plt.xlabel('Number of Clusters') - plt.ylabel('Score') - img_sil = save_path + seq + tit + '_Silhouette.png' - plt.savefig(img_sil, dpi=300) - #plt.showblock=False) - - # Plot Elbow Curve - plt.figure(figsize=(6, 4)) - plt.plot( c_scores, wss_values, marker='*', label='WSS') - plt.title('Elbow Method using WSS') - plt.xlabel('Number of Clusters') - plt.ylabel('Score') - img_wss = save_path + seq + tit + '_WSS.png' - plt.savefig(img_wss, dpi=300) - #plt.showblock=False) + assert save_path is not None, "save_path is not set" + # Plot Silhouette elbow + with matplotlib.rc_context({"font.size": 9}): + fig1 = Figure(figsize=(6, 4), dpi=100) + FigureCanvasAgg(fig1) + ax1 = fig1.add_subplot(111) + ax1.plot(c_scores, s_scores, marker="o", label="Silhouette Score") + ax1.set_title("Elbow Method using Silhouette Score") + ax1.set_xlabel("Number of Clusters") + ax1.set_ylabel("Score") + if np.isfinite(np.nanmax(s_scores)): + best_k = c_scores[int(np.nanargmax(s_scores))] + ax1.axvline(best_k, linestyle="--", color="r", label=f"Best k={best_k}") + ax1.legend() + img_sil = save_path + seq + tit + "_Silhouette.png" + fig1.tight_layout() + fig1.savefig(img_sil, dpi=300, bbox_inches="tight") + fig1.clear() + del fig1 + + # Plot WSS elbow + with matplotlib.rc_context({"font.size": 9}): + fig2 = Figure(figsize=(6, 4), dpi=100) + FigureCanvasAgg(fig2) + ax2 = fig2.add_subplot(111) + ax2.plot(c_scores, wss_values, marker="*", label="WSS") + ax2.set_title("Elbow Method using WSS") + ax2.set_xlabel("Number of Clusters") + ax2.set_ylabel("Score") + if len(wss_values) > 1 and np.isfinite(np.nanmin(wss_values)): + elbow_k = c_scores[int(np.nanargmin(wss_values))] + ax2.axvline( + elbow_k, linestyle="--", color="r", label=f"Elbow k={elbow_k}" + ) + ax2.legend() + img_wss = save_path + seq + tit + "_WSS.png" + fig2.tight_layout() + fig2.savefig(img_wss, dpi=300, bbox_inches="tight") + fig2.clear() + del fig2 + return img_sil, img_wss + def den_array(df, plot=True): - seq = str(df['contig'].unique()) - seq = re.sub('\[|\]|\\s|\\n+|\'|\"', '', seq) + seq = str(df["contig"].unique()) + seq = re.sub("\[|\]|\\s|\\n+|'|\"", "", seq) ####### evaluate clusters ### modes - tit = '_hamming_' + str(df['LID'].unique()) + tit = "_hamming_" + str(df["LID"].unique()) distk, reads, read_list = distancematrix_modes(df, seq=seq) - img_sil, img_wss = get_optimal_clusters(X=distk, read_list=read_list, seq=seq, tit=tit, plot=plot) + img_sil, img_wss = get_optimal_clusters( + X=distk, read_list=read_list, seq=seq, tit=tit, plot=plot + ) optimal_dendrogram(X=distk, seq=seq, tit=tit) #### means - tit = '_kmeans_' + str(df['LID'].unique()) + tit = "_kmeans_" + str(df["LID"].unique()) distk, reads, read_list = distancematrix_means(df, seq=seq) # evaluate clusters - img_sil_k, img_wss_k = get_optimal_clusters(X=distk, read_list=read_list, seq=seq, tit=tit, plot=plot) + img_sil_k, img_wss_k = get_optimal_clusters( + X=distk, read_list=read_list, seq=seq, tit=tit, plot=plot + ) optimal_dendrogram(X=distk, seq=seq, tit=tit) del df diff --git a/DashML/Landscape/Cluster/Cluster_means.py b/DashML/Landscape/Cluster/Cluster_means.py index 62b9618d..b27a5f92 100644 --- a/DashML/Landscape/Cluster/Cluster_means.py +++ b/DashML/Landscape/Cluster/Cluster_means.py @@ -5,8 +5,8 @@ import numpy as np import pandas as pd import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt +from matplotlib.figure import Figure +from matplotlib.backends.backend_agg import FigureCanvasAgg import seaborn as sns from scipy.cluster.hierarchy import dendrogram, linkage from sklearn.cluster import AgglomerativeClustering, KMeans @@ -15,45 +15,61 @@ from scipy.spatial import distance import DashML.Database_fx.Insert_DB as dbins +from pathlib import Path -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) -save_path=None +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) + +save_path: Path | None = None + + +def _ensure_dir(p: Path) -> Path: + p.mkdir(parents=True, exist_ok=True) + return p -#reactivity format for RNAcofold -def scale_reactivities(reactivities): - min = reactivities.min() - max = reactivities.max() - smin = 0 - smax = 2 - reactivities = ((reactivities - min) / (max - min)) * (smax - smin) + smin - return reactivities +def set_save_path(path: str | os.PathLike) -> None: + """Call this once from your app to set a base output path.""" + global save_path + save_path = Path(path) + + +# reactivity format for RNAcofold +def scale_reactivities(reactivities): + rmin = float(np.nanmin(reactivities)) + rmax = float(np.nanmax(reactivities)) + smin, smax = 0.0, 2.0 + if rmax <= rmin: + return np.full_like(reactivities, smin, dtype=float) + return ((reactivities - rmin) / (rmax - rmin)) * (smax - smin) + smin # Create linkage matrix and then plot the dendrogram -def distancematrix(df, seq='hcv'): +def distancematrix(df, seq="hcv"): ##### Compute Distance Between Reads ##### - #reactivity kmeans distance - dx = df.sort_values(by=['read_id', 'position']) - read_list = dx['read_id'].unique() - dx = (dx.pivot(index=["read_id"], columns="position", values="Reactivity") - .rename_axis(columns=None) - .reset_index()) + # reactivity kmeans distance + dx = df.sort_values(by=["read_id", "position"]) + read_list = dx["read_id"].unique() + dx = ( + dx.pivot(index=["read_id"], columns="position", values="Reactivity") + .rename_axis(columns=None) + .reset_index() + ) dx.fillna(value=0, inplace=True) - dx = dx.drop(columns=['read_id']) + dx = dx.drop(columns=["read_id"]) ## create matrix for euclidean distance reads = dx.to_numpy() - #mm = np.divide(mm, np.linalg.norm(mm)) + # mm = np.divide(mm, np.linalg.norm(mm)) ## calculate distances - distk = distance.cdist(reads, reads, 'euclidean') + distk = distance.cdist(reads, reads, "euclidean") return distk, reads, read_list + def plot_dendrogram(model, **kwargs): # Create linkage matrix and then plot the dendrogram @@ -76,124 +92,186 @@ def plot_dendrogram(model, **kwargs): # Plot the corresponding dendrogram dendrogram(linkage_matrix, **kwargs) + #### method ultimately decided by which method most closely approximates the desired function def get_clusters(X, seq, read_list, tit, clust_num): # Heirarchical Clustering - #perform clustering - model = AgglomerativeClustering(n_clusters=clust_num, linkage='complete', metric='precomputed', - compute_distances=True, compute_full_tree=True) - model = model.fit(X) - df = pd.DataFrame(columns=['read_id', 'cluster']) - df['read_id'] = read_list - df['cluster'] = model.labels_ - # plot the top three levels of the dendrogram - plt.figure(figsize=(10, 7)) - plt.title("KMeans Hierarchical Clustering Dendrogram") - plt.xlabel("Number of points in node (or index of point).") - ax = plt.gca() - plot_dendrogram(model, truncate_mode="level", p=5, ax=ax) - img_path = save_path + seq + tit + "_dendrogram.png" - plt.savefig(img_path, dpi=600) - #plt.showblock=False) - return df, img_path + # perform clustering + model = AgglomerativeClustering( + n_clusters=clust_num, + linkage="complete", + metric="precomputed", + compute_distances=True, + compute_full_tree=True, + ).fit(X) + + df = pd.DataFrame({"read_id": read_list, "cluster": model.labels_}) + + with matplotlib.rc_context({"font.size": 9}): + fig = Figure(figsize=(10, 7), dpi=100) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + ax.set_title("KMeans Hierarchical Clustering Dendrogram") + ax.set_xlabel("Number of points in node (or index of point).") + plot_dendrogram(model, truncate_mode="level", p=5, ax=ax) + assert save_path is not None, "save_path is not set — call set_save_path()" + img_path = save_path + seq + tit + "_dendrogram.png" + # out = _ensure_dir(save_path) / f"{seq}{tit}_dendrogram.png" + fig.tight_layout() + fig.savefig(img_path, dpi=600, bbox_inches="tight") + + return df, str(img_path) # hybrid correlation matrix using distance instance of correlation -def clustermap(X, seq='HCV', tit='Kmeans'): +def clustermap(X, seq="HCV", tit="Kmeans"): # cannot contain nan X = np.nan_to_num(X) - z_col = linkage(np.transpose(squareform(X)), method='complete') + z_col = linkage(np.transpose(squareform(X)), method="complete") print(z_col.shape) - z_row = linkage(squareform(X), method='complete') + z_row = linkage(squareform(X), method="complete") print(z_row.shape) - g = sns.clustermap(data=X, row_linkage=z_row, col_linkage=z_col, figsize = (8,8)) - #g.ax_heatmap.set_xlabel("Cluster Labels", labelpad=1, fontweight='extra bold') - #g.ax_heatmap.set_ylabel("Clusters Labels", labelpad=1, fontweight='extra bold') - #plt.title(tit + " Inter-Cluster Distances ") - plt.ylabel("Distance (Min-Max)") - g.savefig(save_path + seq + tit + '_clustermap.png', dpi=600) - #plt.showblock=False) - #print(dir(g)) - #print(g.data2d) + g = sns.clustermap(data=X, row_linkage=z_row, col_linkage=z_col, figsize=(8, 8)) + # g.ax_heatmap.set_xlabel("Cluster Labels", labelpad=1, fontweight='extra bold') + # g.ax_heatmap.set_ylabel("Clusters Labels", labelpad=1, fontweight='extra bold') + # plt.title(tit + " Inter-Cluster Distances ") + # plt.ylabel("Distance (Min-Max)") + g.ax_heatmap.set_ylabel("Distance (Min-Max)") + g.savefig(save_path + seq + tit + "_clustermap.png", dpi=600) + # plt.showblock=False) + # print(dir(g)) + # print(g.data2d) + # g.figure.clear() return -def corrmatrix(seq='HCV', data=None, tit="Kmeans"): - g = sns.heatmap(np.corrcoef(data)) - g.set_xlabel("Reads") - g.set_ylabel("Reads") - g.set_title( "Kmeans Read Correlations") - g.get_figure().savefig(save_path + seq + tit + '_read_corr_matrix.png', dpi=600) - img_path = save_path + seq + tit + '_read_corr_matrix.png' - #plt.showblock=False) - #print(np.corrcoef(data)) +def corrmatrix(seq="HCV", data=None, tit="Kmeans"): + assert data is not None, "data must be provided" + assert save_path is not None, "save_path is not set" + img_path = save_path + seq + tit + "_read_corr_matrix.png" + + with matplotlib.rc_context({"font.size": 9}): # CHANGES HERE + fig = Figure(figsize=(8, 6), dpi=600) # CHANGES HERE + FigureCanvasAgg(fig) # CHANGES HERE + ax = fig.add_subplot(111) # CHANGES HERE + + sns.heatmap(np.corrcoef(data), ax=ax) # CHANGES HERE + ax.set_xlabel("Reads") + ax.set_ylabel("Reads") + ax.set_title("Kmeans Read Correlations") + + fig.tight_layout() # CHANGES HERE + fig.savefig(img_path, dpi=600, bbox_inches="tight") # CHANGES HERE + fig.clear() # CHANGES HERE + del fig + return img_path + #### kmeans to get centroids #### +def corrmatrix(seq="HCV", data=None, tit="Kmeans"): + assert save_path is not None, "save_path is not set — call set_save_path()" + + # DIY Figure/Axes and pass to seaborn + with matplotlib.rc_context({"font.size": 9}): + fig = Figure(figsize=(8, 6), dpi=600) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + sns.heatmap(np.corrcoef(data), ax=ax) + ax.set_xlabel("Reads") + ax.set_ylabel("Reads") + ax.set_title("Kmeans Read Correlations") + fig.tight_layout() + img_path = save_path + seq + tit + "_read_corr_matrix.png" + fig.savefig(img_path, dpi=600, bbox_inches="tight") + return str(img_path) + + def get_centroids(X, clusters, seq): # set cluster labels to reads - df = X.merge(clusters, on=['read_id'], how='left') - df['method'] = 'kmeans' - seq_ids = df['LID'].unique().flatten() + df = X.merge(clusters, on=["read_id"], how="left") + df["method"] = "kmeans" + seq_ids = df["LID"].unique().flatten() # save all clusters dbins.insert_clusters(df, seq_ids) - df_centroid = pd.DataFrame(columns=['LID', 'cluster', 'position', 'centroid']) - seqlen = int(df['position'].max()) - for i in clusters['cluster'].unique(): - dt = df.loc[df['cluster']==i] - lids = dt['LID'].unique() - dt = dt[['read_index', 'position', 'Reactivity']] - dt = dt.sort_values(by=['read_index', 'position']) - dt = (dt.pivot(index=["read_index"], columns="position", values="Reactivity") - .rename_axis(columns=None) - .reset_index()) + df_centroid = pd.DataFrame(columns=["LID", "cluster", "position", "centroid"]) + seqlen = int(df["position"].max()) + for i in clusters["cluster"].unique(): + dt = df.loc[df["cluster"] == i] + lids = dt["LID"].unique() + dt = dt[["read_index", "position", "Reactivity"]] + dt = dt.sort_values(by=["read_index", "position"]) + dt = ( + dt.pivot(index=["read_index"], columns="position", values="Reactivity") + .rename_axis(columns=None) + .reset_index() + ) dt.fillna(value=0, inplace=True) - dt.drop(columns=['read_index'], inplace=True) + dt.drop(columns=["read_index"], inplace=True) X = dt.to_numpy() # get kmeans for all reads in cluster km = KMeans(n_clusters=1, random_state=0, n_init="auto").fit(X) for lid in lids: - df_centroid = pd.concat([df_centroid, pd.DataFrame({'LID': lid, 'cluster': i, 'position': np.arange(0, seqlen+1), - 'centroid': np.array(km.cluster_centers_).flatten()})]) - df_centroid['method'] = 'kmeans' - df_centroid['contig'] = re.sub('\[|\]|\\s|\\n+|\'|\"', '', seq) + df_centroid = pd.concat( + [ + df_centroid, + pd.DataFrame( + { + "LID": lid, + "cluster": i, + "position": np.arange(0, seqlen + 1), + "centroid": np.array(km.cluster_centers_).flatten(), + } + ), + ] + ) + df_centroid["method"] = "kmeans" + df_centroid["contig"] = re.sub("\[|\]|\\s|\\n+|'|\"", "", seq) # save centroids - dbins.insert_centroids(df_centroid,seq_ids) + dbins.insert_centroids(df_centroid, seq_ids) return + ###### Heatmap of Reads Against Summed Reactivity Across different Metrics ###### def heatmap(df, seq, tit): - fig = plt.gcf() - plt.figure(figsize=(60, 60)) - ax = plt.gca() - df = df[['position', 'contig', 'read_id', 'Reactivity']] + assert save_path is not None, "save_path is not set — call set_save_path()" + + df = df[["position", "contig", "read_id", "Reactivity"]] result = df.pivot(index="position", columns="read_id", values="Reactivity") - seq_len = len(df['position'].unique()) - sns.heatmap(result, cmap='crest', - yticklabels=range(0,seq_len), ax=ax) - plt.yticks(fontsize=10) - plt.xticks(fontsize=10) - plt.xlabel("Read Number", fontsize=100, weight='bold') - plt.ylabel("Position", fontsize=100, weight='bold') - plt.title(seq + " Reactivity", fontsize=100, weight='bold') - img_path = save_path + seq + tit + '_heatmap.png' - plt.savefig(img_path, dpi=300) - #plt.showblock=False) - return img_path + seq_len = len(df["position"].unique()) + + img_path = save_path + seq + tit + "_heatmap.png" + + with matplotlib.rc_context({"font.size": 10}): + fig = Figure(figsize=(60, 60), dpi=300) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + sns.heatmap(result, cmap="crest", yticklabels=range(0, seq_len), ax=ax) + ax.set_xlabel("Read Number", fontsize=100, fontweight="bold") + ax.set_ylabel("Position", fontsize=100, fontweight="bold") + ax.set_title(f"{seq} Reactivity", fontsize=100, fontweight="bold") + ax.tick_params(axis="x", labelsize=10) + ax.tick_params(axis="y", labelsize=10) + fig.tight_layout() + fig.savefig(img_path, dpi=300, bbox_inches="tight") + + return str(img_path) def den_array(df, clust_num, plot=True): - tit = '_kmeans_' + str(df['LID'].unique()) - seq = str(df['contig'].unique()) - seq = re.sub('\[|\]|\\s|\\n+|\'|\"', '', seq) + tit = "_kmeans_" + str(df["LID"].unique()) + seq = str(df["contig"].unique()) + seq = re.sub("\[|\]|\\s|\\n+|'|\"", "", seq) distk, reads, read_list = distancematrix(df, seq=seq) - clusters, img_path_dend = get_clusters(X=distk, read_list=read_list, seq=seq, tit=tit, clust_num=clust_num) - #centroids + clusters, img_path_dend = get_clusters( + X=distk, read_list=read_list, seq=seq, tit=tit, clust_num=clust_num + ) + # centroids get_centroids(X=df, clusters=clusters, seq=seq) if plot: img_path_heatmap = heatmap(df, seq, tit) - #get ordered similarity from most to least similar for clusters + # get ordered similarity from most to least similar for clusters clustermap(distk, seq=seq, tit=tit) # correlation between reads img_path_corr = corrmatrix(seq=seq, data=reads, tit=tit) diff --git a/DashML/Landscape/Cluster/Cluster_mode.py b/DashML/Landscape/Cluster/Cluster_mode.py index 4d5fab56..a6591c77 100644 --- a/DashML/Landscape/Cluster/Cluster_mode.py +++ b/DashML/Landscape/Cluster/Cluster_mode.py @@ -5,8 +5,8 @@ import numpy as np import pandas as pd import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt +from matplotlib.figure import Figure +from matplotlib.backends.backend_agg import FigureCanvasAgg import seaborn as sns from scipy.cluster.hierarchy import dendrogram, linkage, cophenet, fcluster from sklearn.cluster import AgglomerativeClustering, KMeans @@ -17,41 +17,48 @@ import DashML.Database_fx.Insert_DB as dbins -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) -save_path=None +save_path = None -#reactivity format for RNAcofold -def scale_reactivities(reactivities): - min = reactivities.min() - max = reactivities.max() - smin = 0 - smax = 2 - reactivities = ((reactivities - min) / (max - min)) * (smax - smin) + smin +# reactivity format for RNAcofold +def scale_reactivities(reactivities): + try: + rmin = float(np.nanmin(reactivities)) + rmax = float(np.nanmax(reactivities)) + smin, smax = 0.0, 2.0 + + reactivities = ((reactivities - rmin) / (rmax - rmin)) * (smax - smin) + smin + except ValueError: # raised if `y` is empty. + pass return reactivities + # Create linkage matrix and then plot the dendrogram -def distancematrix(df, seq='hcv'): +def distancematrix(df, seq="hcv"): ##### Compute Distance Between Reads ##### # reactivity kmode hamming - dx = df.sort_values(by=['read_id', 'position']) - read_list = dx['read_id'].unique() - dx = (dx.pivot(index=["read_id"], columns="position", values="Predict") - .rename_axis(columns=None) - .reset_index()) + dx = df.sort_values(by=["read_id", "position"]) + read_list = dx["read_id"].unique() + dx = ( + dx.pivot(index=["read_id"], columns="position", values="Predict") + .rename_axis(columns=None) + .reset_index() + ) dx.fillna(value=0, inplace=True) - dx = dx.drop(columns=['read_id']) + dx = dx.drop(columns=["read_id"]) mm = dx.to_numpy() - #clustermap - ham = distance.cdist(mm, mm, 'hamming') - #ham = np.divide(ham,np.linalg.norm(ham)) - #print(ham) + # clustermap + ham = distance.cdist(mm, mm, "hamming") + # ham = np.divide(ham,np.linalg.norm(ham)) + # print(ham) return ham, mm, read_list + def plot_dendrogram(model, **kwargs): # Create linkage matrix and then plot the dendrogram @@ -77,122 +84,169 @@ def plot_dendrogram(model, **kwargs): def get_clusters(X, seq, read_list, tit, clust_num, plot=False): # Heirarchical Clustering - #perform clustering - model = AgglomerativeClustering(n_clusters=clust_num, linkage='complete', metric='precomputed', - compute_distances=True, compute_full_tree=True) - model = model.fit(X) - df = pd.DataFrame(columns=['read_id', 'cluster']) - df['read_id'] = read_list - df['cluster'] = model.labels_ + model = AgglomerativeClustering( + n_clusters=clust_num, + linkage="complete", + metric="precomputed", + compute_distances=True, + compute_full_tree=True, + ).fit(X) + + df = pd.DataFrame({"read_id": read_list, "cluster": model.labels_}) + + img_path = None if plot: - # plot the top three levels of the dendrogram - plt.figure(figsize=(10, 7)) - plt.title("Hamming Hierarchical Clustering Dendrogram") - plt.xlabel("Number of points in node (or index of point).") - ax = plt.gca() - plot_dendrogram(model, truncate_mode="level", p=5, ax=ax) + assert save_path is not None, "save_path is not set" img_path = save_path + seq + tit + "_dendrogram.png" - plt.savefig(img_path, dpi=600) - #plt.showblock=False) + + # OO Matplotlib (no pyplot) + with matplotlib.rc_context({"font.size": 9}): + fig = Figure(figsize=(10, 7), dpi=100) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + ax.set_title("Hamming Hierarchical Clustering Dendrogram") + ax.set_xlabel("Number of points in node (or index of point).") + plot_dendrogram(model, truncate_mode="level", p=5, ax=ax) + fig.tight_layout() + fig.savefig(img_path, dpi=600, bbox_inches="tight") + fig.clear() + del fig + return df, img_path + # hybrid correlation matrix using distance instance of correlation -def clustermap(X, seq='HCV', tit='Hamming'): +def clustermap(X, seq="HCV", tit="Hamming"): # cannot contain nan X = np.nan_to_num(X) - z_col = linkage(np.transpose(squareform(X)), method='complete') + + # Linkages from distance matrix + z_col = linkage(np.transpose(squareform(X)), method="complete") print(z_col.shape) - z_row = linkage(squareform(X), method='complete') + z_row = linkage(squareform(X), method="complete") print(z_row.shape) - g = sns.clustermap(data=X, row_linkage=z_row, col_linkage=z_col, figsize = (8,8)) - #g.ax_heatmap.set_xlabel("Cluster Labels", labelpad=1, fontweight='extra bold') - #g.ax_heatmap.set_ylabel("Clusters Labels", labelpad=1, fontweight='extra bold') - #plt.title(tit + " Inter-Cluster Distances ") - plt.ylabel("Distance (Min-Max)") - g.savefig(save_path + seq + tit + '_clustermap.png', dpi=600) - #plt.showblock=False) - #print(dir(g)) - #print(g.data2d) - return + assert save_path is not None, "save_path is not set" + out = save_path + seq + tit + "_clustermap.png" + + # Seaborn without pyplot; save via ClusterGrid + g = sns.clustermap(data=X, row_linkage=z_row, col_linkage=z_col, figsize=(8, 8)) + g.ax_heatmap.set_ylabel("Distance (Min-Max)") + g.savefig(out, dpi=600) + g.figure.clear() + return str(out) + + +def corrmatrix(seq="HCV", data=None, tit="Hamming"): + assert data is not None, "data must be provided" + assert save_path is not None, "save_path is not set" + img_path = save_path + seq + tit + "_read_corr_matrix.png" + with matplotlib.rc_context({"font.size": 9}): + fig = Figure(figsize=(8, 6), dpi=600) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + sns.heatmap(np.corrcoef(data), ax=ax) + ax.set_xlabel("Reads") + ax.set_ylabel("Reads") + ax.set_title("Hamming Read Correlations") + fig.tight_layout() + fig.savefig(img_path, dpi=600, bbox_inches="tight") + fig.clear() + del fig -def corrmatrix(seq='HCV', data=None, tit="Hamming"): - g = sns.heatmap(np.corrcoef(data)) - g.set_xlabel("Reads") - g.set_ylabel("Reads") - g.set_title( "Hamming Read Correlations") - img_path = save_path + seq + tit + '_read_corr_matrix.png' - g.get_figure().savefig(img_path, dpi=600) - #plt.showblock=False) - #print(np.corrcoef(data)) return img_path + #### kmode to get centroids #### def get_centroids(data, clusters, seq): # set cluster labels to reads - df = data.merge(clusters, on=['read_id'], how='left') - df['method'] = 'hamming' - seq_ids = df['LID'].unique().flatten() + df = data.merge(clusters, on=["read_id"], how="left") + df["method"] = "hamming" + seq_ids = df["LID"].unique().flatten() # save all clusters dbins.insert_clusters(df, seq_ids) - df_centroid = pd.DataFrame(columns=['LID', 'cluster', 'position', 'centroid']) - seqlen = int(df['position'].max()) - for i in clusters['cluster'].unique(): - dt = df.loc[df['cluster']==i] - lids = dt['LID'].unique() - dt = dt[['read_id', 'position', 'Predict']] - dt = dt.sort_values(by=['read_id', 'position']) - dt = (dt.pivot(index=["read_id"], columns="position", values="Predict") - .rename_axis(columns=None) - .reset_index()) + df_centroid = pd.DataFrame(columns=["LID", "cluster", "position", "centroid"]) + seqlen = int(df["position"].max()) + for i in clusters["cluster"].unique(): + dt = df.loc[df["cluster"] == i] + lids = dt["LID"].unique() + dt = dt[["read_id", "position", "Predict"]] + dt = dt.sort_values(by=["read_id", "position"]) + dt = ( + dt.pivot(index=["read_id"], columns="position", values="Predict") + .rename_axis(columns=None) + .reset_index() + ) dt.fillna(value=0, inplace=True) - dt.drop(columns=['read_id'], inplace=True) + dt.drop(columns=["read_id"], inplace=True) X = dt.to_numpy() # get kmeans for all reads in cluster - km = KModes(n_clusters=1, init='Huang', n_init=5, verbose=0).fit(X) - #unique, counts = np.unique(km.cluster_centroids_, return_counts=True) - #print(unique, counts) + km = KModes(n_clusters=1, init="Huang", n_init=5, verbose=0).fit(X) + # unique, counts = np.unique(km.cluster_centroids_, return_counts=True) + # print(unique, counts) for lid in lids: - df_centroid = pd.concat([df_centroid, pd.DataFrame({'LID': lid, 'cluster': i, 'position': np.arange(0, seqlen+1), - 'centroid': np.array(km.cluster_centroids_).flatten()})]) - df_centroid['method'] = 'hamming' - df_centroid['contig'] = re.sub('\[|\]|\\s|\\n+|\'|\"', '', seq) + df_centroid = pd.concat( + [ + df_centroid, + pd.DataFrame( + { + "LID": lid, + "cluster": i, + "position": np.arange(0, seqlen + 1), + "centroid": np.array(km.cluster_centroids_).flatten(), + } + ), + ] + ) + df_centroid["method"] = "hamming" + df_centroid["contig"] = re.sub("\[|\]|\\s|\\n+|'|\"", "", seq) # save centroids - dbins.insert_centroids(df_centroid,seq_ids) + dbins.insert_centroids(df_centroid, seq_ids) return ###### Heatmap of Reads Against Summed Reactivity Across different Metrics ###### def heatmap(df, seq, tit): - fig = plt.gcf() - plt.figure(figsize=(60, 60)) - ax = plt.gca() - df = df[['position', 'contig', 'read_id', 'Reactivity']] + assert save_path is not None, "save_path is not set" + img_path = save_path + seq + tit + "_heatmap.png" + + df = df[["position", "contig", "read_id", "Reactivity"]] result = df.pivot(index="position", columns="read_id", values="Reactivity") - seq_len = len(df['position'].unique()) - sns.heatmap(result, cmap='crest', - yticklabels=range(0,seq_len), ax=ax) - plt.yticks(fontsize=10) - plt.xticks(fontsize=10) - plt.xlabel("Read Number", fontsize=100, weight='bold') - plt.ylabel("Position", fontsize=100, weight='bold') - plt.title(seq + " Reactivity", fontsize=100, weight='bold') - img_path = save_path + seq + tit + '_heatmap.png' - plt.savefig(img_path, dpi=300) - #plt.showblock=False) - return img_path + seq_len = len(df["position"].unique()) + + with matplotlib.rc_context({"font.size": 10}): + fig = Figure(figsize=(60, 60), dpi=300) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + + sns.heatmap(result, cmap="crest", yticklabels=range(0, seq_len), ax=ax) + ax.tick_params(axis="x", labelsize=10) + ax.tick_params(axis="y", labelsize=10) + ax.set_xlabel("Read Number", fontsize=100, weight="bold") + ax.set_ylabel("Position", fontsize=100, weight="bold") + ax.set_title(f"{seq} Reactivity", fontsize=100, weight="bold") + + fig.tight_layout() + fig.savefig(img_path, dpi=300, bbox_inches="tight") + fig.clear() + del fig + + return str(img_path) + def den_array(df, clust_num, plot=True): - tit = '_hamming_' + str(df['LID'].unique()) - seq = str(df['contig'].unique()) - seq = re.sub('\[|\]|\\s|\\n+|\'|\"', '', seq) + tit = "_hamming_" + str(df["LID"].unique()) + seq = str(df["contig"].unique()) + seq = re.sub("\[|\]|\\s|\\n+|'|\"", "", seq) distk, reads, read_list = distancematrix(df, seq=seq) - clusters, img_path_dend = get_clusters(X=distk, read_list=read_list, seq=seq, tit=tit, plot=plot, clust_num=clust_num) - #centroids + clusters, img_path_dend = get_clusters( + X=distk, read_list=read_list, seq=seq, tit=tit, plot=plot, clust_num=clust_num + ) + # centroids get_centroids(data=df, clusters=clusters, seq=seq) if plot: img_path_heatmap = heatmap(df, seq, tit) - #get ordered similarity from most to least similar for clusters + # get ordered similarity from most to least similar for clusters clustermap(distk, seq=seq, tit=tit) # correlation between reads img_corr = corrmatrix(seq=seq, data=reads, tit=tit) diff --git a/DashML/Landscape/Cluster/Cluster_native.py b/DashML/Landscape/Cluster/Cluster_native.py index e35adf2f..dd299ed4 100644 --- a/DashML/Landscape/Cluster/Cluster_native.py +++ b/DashML/Landscape/Cluster/Cluster_native.py @@ -1,154 +1,194 @@ import sys, re import platform +from matplotlib.figure import Figure import numpy as np import pandas as pd import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt + +import matplotlib +from matplotlib.figure import Figure +from matplotlib.backends.backend_agg import FigureCanvasAgg import seaborn as sns from scipy.spatial import distance from sklearn.metrics.cluster import adjusted_mutual_info_score import DashML.Database_fx.Insert_DB as dbins import DashML.Database_fx.Select_DB as dbsel -matplotlib.use("Agg") -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) save_path = None -#reactivity format for RNAcofold -def scale_reactivities(reactivities): - min = reactivities.min() - max = reactivities.max() - smin = 0 - smax = 2 - reactivities = ((reactivities - min) / (max - min)) * (smax - smin) + smin +# reactivity format for RNAcofold +def scale_reactivities(reactivities): + try: + rmin = float(np.nanmin(reactivities)) + rmax = float(np.nanmax(reactivities)) + smin, smax = 0.0, 2.0 + + reactivities = ((reactivities - rmin) / (rmax - rmin)) * (smax - smin) + smin + except ValueError: # raised if `y` is empty. + pass return reactivities # Create distance matrices between native sequence and cluster centroids using hamming and kmeans + def distance_matrix(df_kmeans, df_hamming, native_lid, seq="HCV"): native_rx = get_native_rx(native_lid) ##### Compute Distance Between centroids and native ##### - #reactivity kmeans distance - dx = df_kmeans.sort_values(by=['LID','cluster', 'position']) - dx = (dx.pivot(index=["cluster"], columns="position", values="reactivity") - .rename_axis(columns=None) - .reset_index()) + # reactivity kmeans distance + dx = df_kmeans.sort_values(by=["LID", "cluster", "position"]) + dx = ( + dx.pivot(index=["cluster"], columns="position", values="reactivity") + .rename_axis(columns=None) + .reset_index() + ) dx.fillna(value=0, inplace=True) - dx.drop(columns=['cluster'], inplace=True) + dx.drop(columns=["cluster"], inplace=True) mm = dx.to_numpy() # kmeans distance matrix - km = distance.cdist(native_rx, mm, 'euclidean') + km = distance.cdist(native_rx, mm, "euclidean") k_corr = mm - #save centroid distances - df = df_kmeans.sort_values(by=['LID', 'cluster','position']) - df = df[['LID', 'cluster']].reset_index(drop=True) - df = df.groupby(by=['LID', 'cluster']).mean().reset_index() - df['distance'] = np.array(km).flatten() - df['method'] = 'kmeans' + # save centroid distances + df = df_kmeans.sort_values(by=["LID", "cluster", "position"]) + df = df[["LID", "cluster"]].reset_index(drop=True) + df = df.groupby(by=["LID", "cluster"]).mean().reset_index() + df["distance"] = np.array(km).flatten() + df["method"] = "kmeans" dbins.insert_centroid_distance(df) # reactivity kmode distance - dx = df_hamming.sort_values(by=['LID', 'cluster', 'position']) - dx['reactivity'] = np.where(dx['reactivity'] == -1, 1, 0) - dx = (dx.pivot(index=["cluster"], columns="position", values="reactivity") - .rename_axis(columns=None) - .reset_index()) + dx = df_hamming.sort_values(by=["LID", "cluster", "position"]) + dx["reactivity"] = np.where(dx["reactivity"] == -1, 1, 0) + dx = ( + dx.pivot(index=["cluster"], columns="position", values="reactivity") + .rename_axis(columns=None) + .reset_index() + ) dx.fillna(value=0, inplace=True) - dx.drop(columns=['cluster'], inplace=True) + dx.drop(columns=["cluster"], inplace=True) mm = dx.to_numpy() # hamming distance matrix - ham = distance.cdist(native_rx, mm, 'hamming') + ham = distance.cdist(native_rx, mm, "hamming") ham_corr = mm # save centroid distances - df = df_hamming.sort_values(by=['LID', 'cluster', 'position']) - df = df[['LID', 'cluster']].reset_index(drop=True) - df = df.groupby(by=['LID', 'cluster']).mean().reset_index() - df['distance'] = np.array(ham).flatten() - df['method'] = 'hamming' + df = df_hamming.sort_values(by=["LID", "cluster", "position"]) + df = df[["LID", "cluster"]].reset_index(drop=True) + df = df.groupby(by=["LID", "cluster"]).mean().reset_index() + df["distance"] = np.array(ham).flatten() + df["method"] = "hamming" dbins.insert_centroid_distance(df) - #correlation between clusters matrix - ham_map = distance.cdist(mm, mm, 'hamming') + # correlation between clusters matrix + ham_map = distance.cdist(mm, mm, "hamming") return k_corr, ham_corr, km, ham + # get reactivities for native structures using hamming binary reactivities # remove sections with no data for each cluster, only compare non-zero positions # TODO: get reactivities for native structures using native reactivities calculated in Native_Nanopore # remove sections with no data for each cluster, only compare non-zero positions def get_native_rx(lid): df_native = dbsel.select_structure(lid) - df_native = df_native[['base_type']] - df_native['base_type'] = np.where(df_native['base_type']=='S', 1, 0) + df_native = df_native[["base_type"]] + df_native["base_type"] = np.where(df_native["base_type"] == "S", 1, 0) mm = df_native.to_numpy() return mm.T -###### Correlation between Representatives Centroids ###### -def corrmatrixh(seq='HCV', data=None, tit="Hamming", clust_num=5): - plt.clf() - plt.figure(figsize=(20, 20)) - print(data.shape) - arr = np.zeros((clust_num, clust_num)) - for i in range(0, clust_num-1): - for j in range(0, clust_num-1): - arr[i,j] = adjusted_mutual_info_score(data[i,:], data[j,:]) - g = sns.heatmap(arr, annot=True) - g.set_xlabel("Clusters") - g.set_ylabel("Clusters") - g.set_title( tit + " Cluster Correlations (Adjusted Mutual Information)") - g.get_figure().savefig(save_path + seq +'_'+ tit + '_corr_matrix.png', dpi=600) - #plt.showblock=False) - #print(np.corrcoef(data)) - -def corrmatrixk(seq='HCV', data=None, tit="Kmeans", clust_num=5): - plt.clf() - plt.figure(figsize=(20, 20)) - arr = np.zeros((clust_num, clust_num)) + +def corrmatrixh(seq="HCV", data=None, tit="Hamming", clust_num=5): + assert data is not None, "data must be provided" + assert save_path is not None, "save_path is not set" + + arr = np.zeros((clust_num, clust_num), dtype=float) for i in range(0, clust_num - 1): for j in range(0, clust_num - 1): - arr[i, j] = adjusted_mutual_info_score(data[i, :], data[j, :]) - #g = sns.heatmap(np.nan_to_num(np.corrcoef(data)), annot=True) - g = sns.heatmap(arr, annot=True) - g.set_xlabel("Clusters") - g.set_ylabel("Clusters") - g.set_title( tit + " Cluster Correlations (Adjusted Mutual Information)") - img_path = save_path + seq +'_'+ tit + '_corr_matrix.png' - g.get_figure().savefig(img_path, dpi=600) - #plt.showblock=False) - return img_path - + arr[i, j] = adjusted_mutual_info_score( + np.asarray(data[i, :]), np.asarray(data[j, :]) + ) + + with matplotlib.rc_context({"font.size": 12}): + fig = Figure(figsize=(20, 20), dpi=600) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + g = sns.heatmap(arr, annot=True, ax=ax) + g.set_xlabel("Clusters") + g.set_ylabel("Clusters") + g.set_title(f"{tit} Cluster Correlations (Adjusted Mutual Information)") + + out = save_path + seq + "_" + tit + "_corr_matrix.png" + fig.tight_layout() + fig.savefig(out, dpi=600, bbox_inches="tight") + fig.clear() + del fig + + +def corrmatrixk(seq="HCV", data=None, tit="Kmeans", clust_num=5): + assert data is not None, "data must be provided" + assert save_path is not None, "save_path is not set" + + arr = np.zeros((clust_num, clust_num), dtype=float) + for i in range(0, clust_num - 1): + for j in range(0, clust_num - 1): + arr[i, j] = adjusted_mutual_info_score( + np.asarray(data[i, :]), np.asarray(data[j, :]) + ) + + with matplotlib.rc_context({"font.size": 12}): + fig = Figure(figsize=(20, 20), dpi=600) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + g = sns.heatmap(arr, annot=True, ax=ax) + g.set_xlabel("Clusters") + g.set_ylabel("Clusters") + g.set_title(f"{tit} Cluster Correlations (Adjusted Mutual Information)") + img_path = save_path + seq + "_" + tit + "_corr_matrix.png" + fig.tight_layout() + fig.savefig(img_path, dpi=600, bbox_inches="tight") + fig.clear() + del fig + return str(img_path) ###### Heatmap of Cluster Representatives Against native structure ###### def heatmap(X, seq="HCV", method="Kmeans"): - fig = plt.gcf() - plt.figure(figsize=(20, 8)) - plt.title( seq + " " + method + " Centroid Distance from Native Structure.\n") - ax = sns.heatmap(X, cmap='crest', annot=True, square=True) - ax.set(xlabel="Clusters", ylabel=seq) - plt.savefig(save_path + seq + '_' + method + '_heatmap.png', dpi=600) - #plt.show) + assert save_path is not None, "save_path is not set" + + with matplotlib.rc_context({"font.size": 12}): + fig = Figure(figsize=(20, 8), dpi=600) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + ax.set_title(f"{seq} {method} Centroid Distance from Native Structure.\n") + g = sns.heatmap(X, cmap="crest", annot=True, square=True, ax=ax) + g.set(xlabel="Clusters", ylabel=seq) + + out = save_path + seq + "_" + method + "_heatmap.png" + fig.tight_layout() + fig.savefig(out, dpi=600, bbox_inches="tight") + fig.clear() + del fig def den_array(lids, native_lid, plot=True): df = dbsel.select_centroids(lids) - seq = str(df['contig'].unique()) - seq = re.sub('\[|\]|\\s|\\n+|\'|\"', '', seq) - df_kmeans = df[df['method']=='kmeans'] - df_hamming = df[df['method']=='hamming'] - clust_num = len(df['cluster'].unique()) - - k_corr, ham_corr, dist_k, dist_ham = distance_matrix(df_kmeans, df_hamming, native_lid=native_lid, seq=seq) + seq = str(df["contig"].unique()) + seq = re.sub("\[|\]|\\s|\\n+|'|\"", "", seq) + df_kmeans = df[df["method"] == "kmeans"] + df_hamming = df[df["method"] == "hamming"] + clust_num = len(df["cluster"].unique()) + + k_corr, ham_corr, dist_k, dist_ham = distance_matrix( + df_kmeans, df_hamming, native_lid=native_lid, seq=seq + ) if plot: # heatmap of cluster distances from native sequences heatmap(dist_ham, method="Hamming") @@ -156,6 +196,8 @@ def den_array(lids, native_lid, plot=True): # todo: get ordered similarity from most to least similar for clusters # TODO: correlation between reads or clusters??? corrmatrixh(seq=seq, data=ham_corr, tit="Hamming", clust_num=clust_num) - img_corr_means = corrmatrixk(seq=seq, data=k_corr, tit="Kmeans", clust_num=clust_num) + img_corr_means = corrmatrixk( + seq=seq, data=k_corr, tit="Kmeans", clust_num=clust_num + ) return img_corr_means return None diff --git a/DashML/Predict/Gmm_Analysis.py b/DashML/Predict/Gmm_Analysis.py index e828423a..d1d8d76f 100644 --- a/DashML/Predict/Gmm_Analysis.py +++ b/DashML/Predict/Gmm_Analysis.py @@ -1,11 +1,6 @@ import os import sys, re import pandas as pd -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt -import itertools -import matplotlib as mpl import numpy as np from scipy import linalg from sklearn import mixture @@ -13,8 +8,6 @@ import DashML.Database_fx.Insert_DB as dbins - - data_path = None save_path = None @@ -24,15 +17,17 @@ dmso_sequences = None acim_sequences = None seq = None -THRESHOLD = .02 +THRESHOLD = 0.02 + def get_metric(dfr): #### predict modification based on percent modified read depth #### - mean = dfr['percent_modified'].mean() - dfr['Predict'] = np.where(dfr['percent_modified'] > mean, -1, 1) + mean = dfr["percent_modified"].mean() + dfr["Predict"] = np.where(dfr["percent_modified"] > mean, -1, 1) print(dfr.columns) return dfr + #### GMM for each position across all reads ### # predict clusters [mod vs unmod] for each position # plot clusters for a few positions, for paper @@ -45,62 +40,92 @@ def positional_gmm(): dmso = dmso_sequences acim = acim_sequences - d = dmso[['read_index', 'position', 'contig', 'event_level_mean', 'event_length']] - d = d.groupby(by=['read_index', 'position', 'contig']).mean().reset_index() - d = d.pivot(index='read_index', columns=['position'], values=['event_level_mean', 'event_length']) + d = dmso[["read_index", "position", "contig", "event_level_mean", "event_length"]] + d = d.groupby(by=["read_index", "position", "contig"]).mean().reset_index() + d = d.pivot( + index="read_index", + columns=["position"], + values=["event_level_mean", "event_length"], + ) d.fillna(0, inplace=True) - dfn = pd.DataFrame(columns=['LID', 'position', 'contig', 'read_depth', 'percent_modified', 'Predict']) - - for lid in acim['LID'].unique(): - a = acim[acim['LID']==lid] - a = a[['read_index', 'position', 'contig', 'event_level_mean', 'event_length']] - a = a.groupby(by=['read_index', 'position', 'contig']).mean().reset_index() - a = a.pivot(index='read_index', columns=['position'], values=['event_level_mean', 'event_length']) + dfn = pd.DataFrame( + columns=[ + "LID", + "position", + "contig", + "read_depth", + "percent_modified", + "Predict", + ] + ) + + for lid in acim["LID"].unique(): + a = acim[acim["LID"] == lid] + a = a[["read_index", "position", "contig", "event_level_mean", "event_length"]] + a = a.groupby(by=["read_index", "position", "contig"]).mean().reset_index() + a = a.pivot( + index="read_index", + columns=["position"], + values=["event_level_mean", "event_length"], + ) a.fillna(0, inplace=True) - #iterate over positions + # iterate over positions for i in a.columns.get_level_values(1).unique(): - #print(i) + # print(i) ### acim data for all reads at position i ### - b = pd.concat([a['event_level_mean'][i],a['event_length'][i]], axis=1) + b = pd.concat([a["event_level_mean"][i], a["event_length"][i]], axis=1) b = b.to_numpy() ### dmso data for all reads at position i training ### - c = pd.concat([d['event_level_mean'][i], d['event_length'][i]], axis=1) + c = pd.concat([d["event_level_mean"][i], d["event_length"][i]], axis=1) c = c.to_numpy() # training gaussian mixture model, separate into modified vs unmodified num_components = 2 - #{‘full’, ‘tied’, ‘diag’, ‘spherical’} - gmm = BayesianGaussianMixture(n_components = num_components, max_iter=10000, - random_state=1) - model = gmm.fit(c) #train on dmso + # {‘full’, ‘tied’, ‘diag’, ‘spherical’} + gmm = BayesianGaussianMixture( + n_components=num_components, max_iter=10000, random_state=1 + ) + model = gmm.fit(c) # train on dmso dmso_labels = model.predict(c) unique, counts = np.unique(dmso_labels, return_counts=True) - #print("dmso") + # print("dmso") dml = dict(zip(unique, counts)) - dmso_percent_modified = (dml.get(1) if dml.get(1) is not None else 0)/ len(dmso_labels) + dmso_percent_modified = (dml.get(1) if dml.get(1) is not None else 0) / len( + dmso_labels + ) - #n_components = gmm.covariances_ - #num_components = n_components - #print(num_components) + # n_components = gmm.covariances_ + # num_components = n_components + # print(num_components) - #gplot.plot_sigma_vector(gmm.means_, gmm.covariances_) - #gplot.plot_gaussian_mixture(gmm.means_, gmm.covariances_, gmm.weights_) + # gplot.plot_sigma_vector(gmm.means_, gmm.covariances_) + # gplot.plot_gaussian_mixture(gmm.means_, gmm.covariances_, gmm.weights_) - #predictions from gmm - #print("acim") + # predictions from gmm + # print("acim") labels = model.predict(b) unique, counts = np.unique(labels, return_counts=True) al = dict(zip(unique, counts)) - acim_percent_modified = (al.get(1) if al.get(1) is not None else 0) / len(labels) + acim_percent_modified = (al.get(1) if al.get(1) is not None else 0) / len( + labels + ) percent_modified = acim_percent_modified - dmso_percent_modified - #print("dmso_percent_modified: ", dmso_percent_modified) - #print("acim_percent_modified: ", acim_percent_modified) + # print("dmso_percent_modified: ", dmso_percent_modified) + # print("acim_percent_modified: ", acim_percent_modified) if percent_modified > THRESHOLD: predict = -1 - else: predict = 1 - dfn.loc[len(dfn.index)] = [lid, i, seq, len(labels), percent_modified, predict] + else: + predict = 1 + dfn.loc[len(dfn.index)] = [ + lid, + i, + seq, + len(labels), + percent_modified, + predict, + ] print(dfn.head()) print("based on mean") diff --git a/DashML/Predict/Gmm_Analysis2.py b/DashML/Predict/Gmm_Analysis2.py index c0351723..c14ed528 100644 --- a/DashML/Predict/Gmm_Analysis2.py +++ b/DashML/Predict/Gmm_Analysis2.py @@ -1,11 +1,7 @@ import os import sys, re import pandas as pd -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt import itertools -import matplotlib as mpl import numpy as np from scipy import linalg from sklearn import mixture @@ -13,21 +9,22 @@ import DashML.Database_fx.Insert_DB as dbins - #### dmso #### #### dmso #### dmso_sequences = None acim_sequences = None seq = None -THRESHOLD = .02 +THRESHOLD = 0.02 + def get_metric(dfr): #### predict modification based on percent modified read depth #### - mean = dfr['percent_modified'].mean() - dfr['Predict'] = np.where(dfr['percent_modified'] > mean, -1, 1) - #print(dfr.columns) + mean = dfr["percent_modified"].mean() + dfr["Predict"] = np.where(dfr["percent_modified"] > mean, -1, 1) + # print(dfr.columns) return dfr + #### GMM for each position across all reads ### # predict clusters [mod vs unmod] for each position # plot clusters for a few positions, for paper @@ -40,72 +37,102 @@ def positional_gmm(): dmso = dmso_sequences acim = acim_sequences - d = dmso[['read_index', 'position', 'contig', 'event_level_mean', 'event_length']] - d = d.groupby(by=['read_index', 'position', 'contig']).mean().reset_index() - #d = d.pivot(index='read_index', columns=['position'], values=['event_level_mean', 'event_length']) + d = dmso[["read_index", "position", "contig", "event_level_mean", "event_length"]] + d = d.groupby(by=["read_index", "position", "contig"]).mean().reset_index() + # d = d.pivot(index='read_index', columns=['position'], values=['event_level_mean', 'event_length']) d.fillna(0, inplace=True) - acim = acim[['LID','read_index', 'position', 'contig', 'event_level_mean', 'event_length']] - acim = acim.groupby(by=['read_index', 'position', 'contig']).mean().reset_index() + acim = acim[ + ["LID", "read_index", "position", "contig", "event_level_mean", "event_length"] + ] + acim = acim.groupby(by=["read_index", "position", "contig"]).mean().reset_index() acim.fillna(0, inplace=True) - dfn = pd.DataFrame(columns=['LID', 'position', 'contig', 'read_depth', 'percent_modified', 'Predict']) - - for lid in acim['LID'].unique(): - a = acim[acim['LID']==lid] - a = a.loc[a['event_level_mean']!=0] - - #iterate over positions - for pos in a['position'].unique(): - at = a.loc[a['position']==pos] + dfn = pd.DataFrame( + columns=[ + "LID", + "position", + "contig", + "read_depth", + "percent_modified", + "Predict", + ] + ) + + for lid in acim["LID"].unique(): + a = acim[acim["LID"] == lid] + a = a.loc[a["event_level_mean"] != 0] + + # iterate over positions + for pos in a["position"].unique(): + at = a.loc[a["position"] == pos] ### acim data for all reads at position i ### - b = pd.concat([at['event_level_mean'],at['event_length']], axis=1) + b = pd.concat([at["event_level_mean"], at["event_length"]], axis=1) b = b.to_numpy() ### dmso data for all reads at position i training ### - c = pd.concat([d.loc[d['position']==pos, 'event_level_mean'],d.loc[d['position']==pos,'event_length']], axis=1) - if len(c) < 2: #gmm requires at least 2 + c = pd.concat( + [ + d.loc[d["position"] == pos, "event_level_mean"], + d.loc[d["position"] == pos, "event_length"], + ], + axis=1, + ) + if len(c) < 2: # gmm requires at least 2 dfn.loc[len(dfn.index)] = [lid, pos, seq, len(c), 0, 0] else: c = c.to_numpy() # training gaussian mixture model, separate into modified vs unmodified num_components = 2 - #{‘full’, ‘tied’, ‘diag’, ‘spherical’} - gmm = BayesianGaussianMixture(n_components = num_components, max_iter=10000, - random_state=1) - model = gmm.fit(c) #train on dmso + # {‘full’, ‘tied’, ‘diag’, ‘spherical’} + gmm = BayesianGaussianMixture( + n_components=num_components, max_iter=10000, random_state=1 + ) + model = gmm.fit(c) # train on dmso dmso_labels = model.predict(c) - #print('c',dmso_labels) + # print('c',dmso_labels) unique, counts = np.unique(dmso_labels, return_counts=True) - #print("dmso") + # print("dmso") dml = dict(zip(unique, counts)) - #print(dml) - dmso_percent_modified = (dml.get(1) if dml.get(1) is not None else 0)/ len(dmso_labels) + # print(dml) + dmso_percent_modified = ( + dml.get(1) if dml.get(1) is not None else 0 + ) / len(dmso_labels) - #n_components = gmm.covariances_ - #num_components = n_components - #print(num_components) + # n_components = gmm.covariances_ + # num_components = n_components + # print(num_components) - #gplot.plot_sigma_vector(gmm.means_, gmm.covariances_) - #gplot.plot_gaussian_mixture(gmm.means_, gmm.covariances_, gmm.weights_) + # gplot.plot_sigma_vector(gmm.means_, gmm.covariances_) + # gplot.plot_gaussian_mixture(gmm.means_, gmm.covariances_, gmm.weights_) - #predictions from gmm - #print("acim") + # predictions from gmm + # print("acim") labels = model.predict(b) - #print('b',labels) + # print('b',labels) unique, counts = np.unique(labels, return_counts=True) al = dict(zip(unique, counts)) - #print('b', al) - acim_percent_modified = (al.get(1) if al.get(1) is not None else 0) / len(labels) + # print('b', al) + acim_percent_modified = ( + al.get(1) if al.get(1) is not None else 0 + ) / len(labels) percent_modified = acim_percent_modified - dmso_percent_modified - #print("dmso_percent_modified: ", dmso_percent_modified) - #print("acim_percent_modified: ", acim_percent_modified) + # print("dmso_percent_modified: ", dmso_percent_modified) + # print("acim_percent_modified: ", acim_percent_modified) if percent_modified > THRESHOLD: predict = -1 - else: predict = 1 - dfn.loc[len(dfn.index)] = [lid, pos, seq, len(labels), percent_modified, predict] - - #print(dfn.head()) - #print("based on mean") + else: + predict = 1 + dfn.loc[len(dfn.index)] = [ + lid, + pos, + seq, + len(labels), + percent_modified, + predict, + ] + + # print(dfn.head()) + # print("based on mean") dbins.insert_gmm(dfn) diff --git a/DashML/Predict/Lof.py b/DashML/Predict/Lof.py index 46835978..e8576fb2 100644 --- a/DashML/Predict/Lof.py +++ b/DashML/Predict/Lof.py @@ -2,16 +2,13 @@ import os.path import numpy as np import pandas as pd -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt from matplotlib.legend_handler import HandlerPathCollection from sklearn.neighbors import LocalOutlierFactor import DashML.Database_fx.Insert_DB as dbins -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) #### source #### dmso_sequences = None @@ -22,16 +19,22 @@ def get_metric(dfr): #### aggregate counts of Predict, read_depth, #### another decider based on the percentage of modified reads ### - dfr = dfr[['position', 'contig', 'Predict']] - dfr['Predict'] = dfr['Predict'].astype('category') - dfr = dfr.groupby(['position', 'contig', 'Predict'], observed=False).size().unstack(fill_value=0) + dfr = dfr[["position", "contig", "Predict"]] + dfr["Predict"] = dfr["Predict"].astype("category") + dfr = ( + dfr.groupby(["position", "contig", "Predict"], observed=False) + .size() + .unstack(fill_value=0) + ) dfr.reset_index(inplace=True) - dfr['read_depth'] = dfr[-1] + dfr[1] - dfr['percent_modified'] = dfr[-1] / dfr['read_depth'] + dfr["read_depth"] = dfr[-1] + dfr[1] + dfr["percent_modified"] = dfr[-1] / dfr["read_depth"] #### predict modification based on percent modified read depth #### - mean = dfr['percent_modified'].mean() - dfr['Predict'] = np.where(dfr['percent_modified'] > mean, -1, 1) + mean = dfr["percent_modified"].mean() + dfr["Predict"] = np.where(dfr["percent_modified"] > mean, -1, 1) return dfr + + def update_legend_marker_size(handle, orig): "Customize size of the legend marker" handle.update_from(orig) @@ -44,13 +47,13 @@ def get_novelty(): def novelty_signal(): ### train on all dmso #### - dftrain = dmso[['event_level_mean']] + dftrain = dmso[["event_level_mean"]] X_train = dftrain.to_numpy() - dfx = acim[['event_level_mean']] + dfx = acim[["event_level_mean"]] X = dfx.to_numpy() # fit the model for novelty detection (novelty=True) - clf = LocalOutlierFactor(n_neighbors=5, novelty=True, contamination=.4) + clf = LocalOutlierFactor(n_neighbors=5, novelty=True, contamination=0.4) clf.fit(X_train) # DO NOT use predict, decision_function and score_samples on X_train as this # would give wrong results but only on new unseen data (not used in X_train), @@ -59,24 +62,24 @@ def novelty_signal(): # n_errors = (y_pred_test != ground_truth).sum() X_scores = clf.negative_outlier_factor_ - #df = pd.DataFrame({"Shape_Map": ground_truth, "Lof_Novelty": y_pred_test, "BaseType":np.array(dft["BaseType"])}) - acim['predict_signal'] = y_pred_test - acim['varna_signal'] = np.where(y_pred_test == -1, 1, 0) - #dft.to_csv(save_path + seq + "_lof_signal.csv") - #local metric + # df = pd.DataFrame({"Shape_Map": ground_truth, "Lof_Novelty": y_pred_test, "BaseType":np.array(dft["BaseType"])}) + acim["predict_signal"] = y_pred_test + acim["varna_signal"] = np.where(y_pred_test == -1, 1, 0) + # dft.to_csv(save_path + seq + "_lof_signal.csv") + # local metric # dft = get_metric(dft) # sp = save_path + seq + "_lof_signal_metric.csv" # mx.get_Metric(dft, seq, sp) def novelty_dwell(): ### train on all dmso #### - dftrain = dmso[['event_length']] + dftrain = dmso[["event_length"]] X_train = dftrain.to_numpy() - dfx = acim[['event_length']] + dfx = acim[["event_length"]] X = dfx.to_numpy() # fit the model for novelty detection (novelty=True) - clf = LocalOutlierFactor(n_neighbors=5, novelty=True, contamination=.4) + clf = LocalOutlierFactor(n_neighbors=5, novelty=True, contamination=0.4) clf.fit(X_train) # DO NOT use predict, decision_function and score_samples on X_train as this # would give wrong results but only on new unseen data (not used in X_train), @@ -84,12 +87,22 @@ def novelty_dwell(): y_pred_test = clf.predict(X) # df = pd.DataFrame({"Shape_Map": ground_truth, "Lof_Novelty": y_pred_test, "BaseType":np.array(dft["BaseType"])}) - #dft = acim[['position', 'contig', 'read_index']] - acim['predict_dwell'] = y_pred_test - acim['varna_dwell'] = np.where(y_pred_test == -1, 1, 0) - + # dft = acim[['position', 'contig', 'read_index']] + acim["predict_dwell"] = y_pred_test + acim["varna_dwell"] = np.where(y_pred_test == -1, 1, 0) novelty_signal() novelty_dwell() - acim = acim[['LID', 'contig', 'read_index', 'position','predict_signal', 'predict_dwell', 'varna_dwell', 'varna_signal']] + acim = acim[ + [ + "LID", + "contig", + "read_index", + "position", + "predict_signal", + "predict_dwell", + "varna_dwell", + "varna_signal", + ] + ] dbins.insert_lof(acim) diff --git a/DashML/Predict/Peak_Analysis_BC.py b/DashML/Predict/Peak_Analysis_BC.py index ddfeec6b..3bb8a4cc 100644 --- a/DashML/Predict/Peak_Analysis_BC.py +++ b/DashML/Predict/Peak_Analysis_BC.py @@ -7,9 +7,6 @@ import math import scipy.signal from scipy import stats -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt from sklearn import preprocessing from sklearn import metrics from sklearn.metrics import RocCurveDisplay @@ -21,8 +18,8 @@ import DashML.Database_fx.Insert_DB as dbins -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) #### source #### @@ -33,17 +30,22 @@ def unit_vector_norm(x): # unit vector normalization - x = np.where(x != 0, ((x- np.min(x))/ (np.max(x)- np.min(x))), 0) - #x = stats.zscore(x) + x = np.where(x != 0, ((x - np.min(x)) / (np.max(x) - np.min(x))), 0) + # x = stats.zscore(x) return x + def ksm(df): bc_acim = df bc_dmso = dmso_sequences - c = ['basecall_reactivity'] - xr = bc_dmso[c].to_numpy().flatten() # unit_vector_norm(bc_dmso[c].to_numpy()).flatten() - yr = bc_acim[c].to_numpy().flatten() # unit_vector_norm(bc_acim[c].to_numpy()).flatten() + c = ["basecall_reactivity"] + xr = ( + bc_dmso[c].to_numpy().flatten() + ) # unit_vector_norm(bc_dmso[c].to_numpy()).flatten() + yr = ( + bc_acim[c].to_numpy().flatten() + ) # unit_vector_norm(bc_acim[c].to_numpy()).flatten() end = len(xr) ##### manual cdf to measure ks ##### # CDF @@ -52,27 +54,30 @@ def ksm(df): y1 = np.sort(yr) def ecdf(x, v): - res = np.searchsorted(x, v, side='right') / x.size + res = np.searchsorted(x, v, side="right") / x.size return res delta_ecdf = np.subtract(yr, xr) - log_dens = [] n = 0 width = 2 - for i in range(0,len(delta_ecdf)): + for i in range(0, len(delta_ecdf)): m = n n = n + width if n < len(delta_ecdf) and m < len(delta_ecdf): - kde = KernelDensity(kernel="gaussian", bandwidth="silverman").fit(delta_ecdf.reshape(-1, 1)[m:n]) - l= kde.score_samples(delta_ecdf.reshape(-1,1)[m:n]) + kde = KernelDensity(kernel="gaussian", bandwidth="silverman").fit( + delta_ecdf.reshape(-1, 1)[m:n] + ) + l = kde.score_samples(delta_ecdf.reshape(-1, 1)[m:n]) log_dens.append(l) else: - kde = KernelDensity(kernel="gaussian", bandwidth="silverman").fit(delta_ecdf.reshape(-1, 1)[m-width:]) + kde = KernelDensity(kernel="gaussian", bandwidth="silverman").fit( + delta_ecdf.reshape(-1, 1)[m - width :] + ) l = kde.score_samples(delta_ecdf.reshape(-1, 1)[m:]) log_dens.append(l) - break; + break tmp = [] for l in log_dens: @@ -81,10 +86,10 @@ def ecdf(x, v): p = find_peaks(delta_ecdf, plateau_size=[0, 10]) - #real values of delta_ecdf normalized + # real values of delta_ecdf normalized peaks_value = unit_vector_norm(delta_ecdf) - #set peaks in vector of all bases + # set peaks in vector of all bases s1 = np.ones(len(delta_ecdf)) s1[p[0]] = -1 # s2 = np.zeros(len(delta_ecdf)) @@ -92,40 +97,58 @@ def ecdf(x, v): # if i in p[0]: # s2[i] = n - #print(sequence) + # print(sequence) dfc = pd.DataFrame() - #print(p[1].get('peak_heights')) - #dfc['peak_weight'] = s2 - dfc['LID'] = df['LID'] - dfc['position'] = df['position'] - dfc['contig'] = contig - dfc['is_peak'] = s1 - dfc['peak_height'] = delta_ecdf - dfc['insertion'] = bc_dmso['insertion'].astype('float32').to_numpy() - bc_acim['insertion'].astype('float32').to_numpy() - dfc['mismatch'] = bc_dmso['mismatch'].astype('float32').to_numpy() - bc_acim['mismatch'].astype('float32').to_numpy() - dfc['deletion'] = bc_dmso['deletion'].astype('float32').to_numpy() - bc_acim['deletion'].astype('float32').to_numpy() - dfc['quality'] = bc_dmso['quality'].astype('float32').to_numpy() - bc_acim['quality'].astype('float32').to_numpy() - dfc['basecall_reactivity'] = bc_dmso['basecall_reactivity'].astype('float32').to_numpy() - bc_acim['basecall_reactivity'].astype('float32').to_numpy() - dfc['aligned_reads'] = bc_dmso['aligned_reads'].astype('float32').to_numpy() - bc_acim['aligned_reads'].astype('float32').to_numpy() - #dfc['is_peak2'] = np.where(dfc['peak_height'] > .008, 1, 0) - #dfc.to_csv(sequence + "_weightcompare.csv") - #print(dfc.head(100)) + # print(p[1].get('peak_heights')) + # dfc['peak_weight'] = s2 + dfc["LID"] = df["LID"] + dfc["position"] = df["position"] + dfc["contig"] = contig + dfc["is_peak"] = s1 + dfc["peak_height"] = delta_ecdf + dfc["insertion"] = ( + bc_dmso["insertion"].astype("float32").to_numpy() + - bc_acim["insertion"].astype("float32").to_numpy() + ) + dfc["mismatch"] = ( + bc_dmso["mismatch"].astype("float32").to_numpy() + - bc_acim["mismatch"].astype("float32").to_numpy() + ) + dfc["deletion"] = ( + bc_dmso["deletion"].astype("float32").to_numpy() + - bc_acim["deletion"].astype("float32").to_numpy() + ) + dfc["quality"] = ( + bc_dmso["quality"].astype("float32").to_numpy() + - bc_acim["quality"].astype("float32").to_numpy() + ) + dfc["basecall_reactivity"] = ( + bc_dmso["basecall_reactivity"].astype("float32").to_numpy() + - bc_acim["basecall_reactivity"].astype("float32").to_numpy() + ) + dfc["aligned_reads"] = ( + bc_dmso["aligned_reads"].astype("float32").to_numpy() + - bc_acim["aligned_reads"].astype("float32").to_numpy() + ) + # dfc['is_peak2'] = np.where(dfc['peak_height'] > .008, 1, 0) + # dfc.to_csv(sequence + "_weightcompare.csv") + # print(dfc.head(100)) return dfc + def get_bc_reactivity_peaks(): # basecall error file acim = acim_sequences - #print(acim.columns) + # print(acim.columns) dft = pd.DataFrame() - for lid in acim['LID'].unique(): - df = acim[acim['LID']==lid] + for lid in acim["LID"].unique(): + df = acim[acim["LID"] == lid] peak_weight = ksm(df=df) if peak_weight is not None: dft = pd.concat([dft, peak_weight], ignore_index=True) - - dft['Predict'] = np.where(dft['is_peak'] == -1, -1, 1) - dft['VARNA'] = np.where(dft['is_peak'] == -1, 1, 0) + dft["Predict"] = np.where(dft["is_peak"] == -1, -1, 1) + dft["VARNA"] = np.where(dft["is_peak"] == -1, 1, 0) dbins.insert_basecall(dft) return diff --git a/DashML/Predict/Predict_BPP.py b/DashML/Predict/Predict_BPP.py index fcc7e854..26b5335b 100644 --- a/DashML/Predict/Predict_BPP.py +++ b/DashML/Predict/Predict_BPP.py @@ -1,8 +1,5 @@ import sys, os import traceback -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt import pandas as pd import numpy as np from DashML.Predict import Predict_Fold as pbpp @@ -19,51 +16,53 @@ MAX_THREADS = 100 + +# reactivity format for RNAcofold def scale_reactivities(reactivities): - min = reactivities.min() - max = reactivities.max() - smin = 0 - smax = 2 + try: + rmin = float(np.nanmin(reactivities)) + rmax = float(np.nanmax(reactivities)) + smin, smax = 0.0, 2.0 - reactivities = ((reactivities - min) / (max - min)) * (smax - smin) + smin + reactivities = ((reactivities - rmin) / (rmax - rmin)) * (smax - smin) + smin + except ValueError: # raised if `y` is empty. + pass return reactivities + ### TODO incomplete def get_prob_pool_sep(df, lids, threshold, continue_reads=False): # get probabilities for each lid-lid interactions ids = [] ssids = [] - if continue_reads== True: + if continue_reads == True: dfr = dbsel.select_continued_reads(lids) - dfx = df.groupby(['LID', 'read_index']).size().reset_index() - dfx = dfx.merge(dfr, on=['LID', 'read_index'], how='left') - dfx['completed'].fillna(0, inplace=True) - dfx = dfx.loc[dfx['completed'] == 0] - dfx = dfx[['LID', 'read_index']] - ids = dfx.loc[:, 'LID':'read_index'].values.tolist() + dfx = df.groupby(["LID", "read_index"]).size().reset_index() + dfx = dfx.merge(dfr, on=["LID", "read_index"], how="left") + dfx["completed"].fillna(0, inplace=True) + dfx = dfx.loc[dfx["completed"] == 0] + dfx = dfx[["LID", "read_index"]] + ids = dfx.loc[:, "LID":"read_index"].values.tolist() del dfr del dfx else: - dfx = df.groupby(['LID', 'read_index']).size().reset_index() - dfx = dfx[['LID', 'read_index']] - ids = dfx.loc[:, 'LID':'read_index'].values.tolist() + dfx = df.groupby(["LID", "read_index"]).size().reset_index() + dfx = dfx[["LID", "read_index"]] + ids = dfx.loc[:, "LID":"read_index"].values.tolist() del dfx - def get_rx(lid, read): try: - dt = df.loc[(df['LID'] == int(lid)) & (df['read_index'] == int(read))] - reactivity = dt['Reactivity_score'].to_numpy().flatten() + dt = df.loc[(df["LID"] == int(lid)) & (df["read_index"] == int(read))] + reactivity = dt["Reactivity_score"].to_numpy().flatten() reactivity = scale_reactivities(reactivity) ssid = pbpp.get_probabilities(lid, lid, reactivity=reactivity, read=read) except Exception as e: raise Exception(str(e)) - print("rx exception",e) + print("rx exception", e) return - - def get_probpool(): with ThreadPoolExecutor(max_workers=MAX_THREADS) as executor: for lid, read in ids: @@ -72,51 +71,56 @@ def get_probpool(): get_probpool() # add bpp to rdf - dbins.insert_read_depth_full_update(lids, threshold, rx_threshold=.7) + dbins.insert_read_depth_full_update(lids, threshold, rx_threshold=0.7) return df + def get_prob_pool(df, lids, threshold, continue_reads=False): # get probabilities for each lid-lid interactions ids = [] ssids = [] df_ps = pd.DataFrame() - - if continue_reads== True: + if continue_reads == True: dfr = dbsel.select_continued_reads(lids) - dfx = df.groupby(['LID', 'read_index']).size().reset_index() - dfx = dfx.merge(dfr, on=['LID', 'read_index'], how='left') - dfx['completed'].fillna(0, inplace=True) - dfx = dfx.loc[dfx['completed'] == 0] - dfx = dfx[['LID', 'read_index']] - ids = dfx.loc[:, 'LID':'read_index'].values.tolist() + dfx = df.groupby(["LID", "read_index"]).size().reset_index() + dfx = dfx.merge(dfr, on=["LID", "read_index"], how="left") + dfx["completed"].fillna(0, inplace=True) + dfx = dfx.loc[dfx["completed"] == 0] + dfx = dfx[["LID", "read_index"]] + ids = dfx.loc[:, "LID":"read_index"].values.tolist() del dfr del dfx else: - dfx = df.groupby(['LID', 'read_index']).size().reset_index() - dfx = dfx[['LID', 'read_index']] - ids = dfx.loc[:, 'LID':'read_index'].values.tolist() + dfx = df.groupby(["LID", "read_index"]).size().reset_index() + dfx = dfx[["LID", "read_index"]] + ids = dfx.loc[:, "LID":"read_index"].values.tolist() del dfx - def get_rx(lid, read): - dt = df.loc[(df['LID'] == int(lid)) & (df['read_index'] == int(read))] - reactivity = dt['Reactivity_score'].to_numpy().flatten() + dt = df.loc[(df["LID"] == int(lid)) & (df["read_index"] == int(read))] + reactivity = dt["Reactivity_score"].to_numpy().flatten() reactivity = scale_reactivities(reactivity) ssid = pbpp.get_probabilities(lid, lid, reactivity=reactivity, read=read) df_ps = dbsel.select_max_structure_probabilities(ssid) if len(df_ps) > 0: - df_ps['position'] = df_ps['position'] - 1 - if 'base_pair_prob' in dt.columns: - dt = dt.drop(columns='base_pair_prob') - dt = dt.merge(df_ps, on=['LID', 'position', 'read_index'], how='left') - dt['completed'] = 1 + df_ps["position"] = df_ps["position"] - 1 + if "base_pair_prob" in dt.columns: + dt = dt.drop(columns="base_pair_prob") + dt = dt.merge(df_ps, on=["LID", "position", "read_index"], how="left") + dt["completed"] = 1 ### Adjust Predict Based on Threshold #### - dt['Predict'] = np.where( - ((dt['Predict'] == -1) & (dt['Reactivity_score'] < .7) & (dt['base_pair_prob'] > threshold)), 1, - dt['Predict']) - if 'base_pair_prob' in dt.columns: - dt['base_pair_prob'].fillna(0, inplace=True) + dt["Predict"] = np.where( + ( + (dt["Predict"] == -1) + & (dt["Reactivity_score"] < 0.7) + & (dt["base_pair_prob"] > threshold) + ), + 1, + dt["Predict"], + ) + if "base_pair_prob" in dt.columns: + dt["base_pair_prob"].fillna(0, inplace=True) try: dbins.insert_read_depth_full(dt) except Exception as e: @@ -133,79 +137,112 @@ def get_probpool(): return df -def get_predict_probabilities(df, lids, threshold, reactivity=False, continue_reads=False): - df.sort_values(by=['LID', 'position'], inplace=True) + +def get_predict_probabilities( + df, lids, threshold, reactivity=False, continue_reads=False +): + df.sort_values(by=["LID", "position"], inplace=True) # get probabilities for each lid-lid interactions df_ps = pd.DataFrame() if reactivity == True: - #read_depth_full - if 'Reactivity_score' in df.columns: #reactivity for each read read_depth_full - df = get_prob_pool(df, lids, threshold, continue_reads=continue_reads) - - #read_depth - elif 'RNAFold_Shape_Reactivity' in df.columns: # averaged reactivity read_depth, single value + # read_depth_full + if "Reactivity_score" in df.columns: # reactivity for each read read_depth_full + df = get_prob_pool(df, lids, threshold, continue_reads=continue_reads) + + # read_depth + elif ( + "RNAFold_Shape_Reactivity" in df.columns + ): # averaged reactivity read_depth, single value for lid in lids: - reactivity = df.loc[df['LID']==int(lid), 'RNAFold_Shape_Reactivity'].to_numpy().flatten() - ssid = pbpp.get_probabilities(lid, lid, reactivity=reactivity, read=-1) - df_ps = pd.concat([df_ps, dbsel.select_max_structure_probabilities(ssid)]) - df_ps['position'] = df_ps['position'] - 1 - df['read_index'] = -1 - df = df.merge(df_ps, on=['LID', 'position', 'read_index'], how='left') + reactivity = ( + df.loc[df["LID"] == int(lid), "RNAFold_Shape_Reactivity"] + .to_numpy() + .flatten() + ) + ssid = pbpp.get_probabilities(lid, lid, reactivity=reactivity, read=-1) + df_ps = pd.concat( + [df_ps, dbsel.select_max_structure_probabilities(ssid)] + ) + df_ps["position"] = df_ps["position"] - 1 + df["read_index"] = -1 + df = df.merge(df_ps, on=["LID", "position", "read_index"], how="left") ### Adjust Predict Based on Threshold #### - df['Predict'] = np.where(((df['Predict'] == -1) & (df['RNAFold_Shape_Reactivity']<.8) - & (df['base_pair_prob'] > threshold)), 1, df['Predict']) - df['base_pair_prob'].fillna(0, inplace=True) - - else: # no reactivity for unmodified and averaged if unmodified reactivity not calculated + df["Predict"] = np.where( + ( + (df["Predict"] == -1) + & (df["RNAFold_Shape_Reactivity"] < 0.8) + & (df["base_pair_prob"] > threshold) + ), + 1, + df["Predict"], + ) + df["base_pair_prob"].fillna(0, inplace=True) + + else: # no reactivity for unmodified and averaged if unmodified reactivity not calculated for lid in lids: ssid = pbpp.get_probabilities(lid, lid, reactivity=None, read=-2) df_ps = pd.concat([df_ps, dbsel.select_max_structure_probabilities(ssid)]) - df_ps['position'] = df_ps['position'] -1 - df['read_index'] = -2 - df = df.merge(df_ps, on=['LID', 'position', 'read_index'], how='left') + df_ps["position"] = df_ps["position"] - 1 + df["read_index"] = -2 + df = df.merge(df_ps, on=["LID", "position", "read_index"], how="left") print(df) ### Adjust Predict Based on Threshold #### - df['Predict'] = np.where( - ((df['Predict'] == -1) & (df['Reactivity_score'] < .7) & (df['base_pair_prob'] > threshold)), 1, - df['Predict']) - df['base_pair_prob'].fillna(0, inplace=True) + df["Predict"] = np.where( + ( + (df["Predict"] == -1) + & (df["Reactivity_score"] < 0.7) + & (df["base_pair_prob"] > threshold) + ), + 1, + df["Predict"], + ) + df["base_pair_prob"].fillna(0, inplace=True) return df -#separate run for deriving vienna probabilities for predicted lids -def get_predict_probability(lids, threshold=.95, continue_reads=False): + +# separate run for deriving vienna probabilities for predicted lids +def get_predict_probability(lids, threshold=0.95, continue_reads=False): try: if isinstance(lids, int): lids = str(lids) lids = [lids] df = dbsel.select_read_depth_full_all(lids) - df.sort_values(by=['LID', 'position'], inplace=True) + df.sort_values(by=["LID", "position"], inplace=True) # get probabilities for each lid-lid interactions df_ps = pd.DataFrame() # check if reactivities exist in read depth full tables, if not return - if 'Reactivity_score' in df.columns: # reactivity for each read read_depth_full + if "Reactivity_score" in df.columns: # reactivity for each read read_depth_full df = get_prob_pool_sep(df, lids, threshold, continue_reads=continue_reads) - else: # no reactivity for unmodified and averaged if unmodified reactivity not calculated + else: # no reactivity for unmodified and averaged if unmodified reactivity not calculated for lid in lids: ssid = pbpp.get_probabilities(lid, lid, reactivity=None, read=-2) if ssid is not None: - df_ps = pd.concat([df_ps, dbsel.select_max_structure_probabilities(ssid)]) - df_ps['position'] = df_ps['position'] -1 - df['read_index'] = -2 + df_ps = pd.concat( + [df_ps, dbsel.select_max_structure_probabilities(ssid)] + ) + df_ps["position"] = df_ps["position"] - 1 + df["read_index"] = -2 if len(df_ps) > 0: - df = df.merge(df_ps, on=['LID', 'position', 'read_index'], how='left') + df = df.merge(df_ps, on=["LID", "position", "read_index"], how="left") print(df) ### Adjust Predict Based on Threshold #### - df['Predict'] = np.where( - ((df['Predict'] == -1) & (df['Reactivity_score'] < .7) & (df['base_pair_prob'] > threshold)), 1, - df['Predict']) - df['base_pair_prob'].fillna(0, inplace=True) + df["Predict"] = np.where( + ( + (df["Predict"] == -1) + & (df["Reactivity_score"] < 0.7) + & (df["base_pair_prob"] > threshold) + ), + 1, + df["Predict"], + ) + df["base_pair_prob"].fillna(0, inplace=True) except Exception as e: raise Exception(str(e)) print("Probability Error ", e) diff --git a/DashML/Predict/Predicts.py b/DashML/Predict/Predicts.py index 76b01a63..dfb411bc 100644 --- a/DashML/Predict/Predicts.py +++ b/DashML/Predict/Predicts.py @@ -10,7 +10,7 @@ ### Note: Predict > 3 is noisier with more coverage, better for average predictions ### Predict > 4 stabilizes in the deconvolution offers improved predictions -#TODO: get all non reactive signle molecule base pair probabilities +# TODO: get all non reactive signle molecule base pair probabilities # and get intermolecular non reactive bp probabilities for caomparison #### TODO ## could use dmso reactivities ???, need to calculate dmso reactivities and integrate @@ -19,46 +19,63 @@ # TODO: optimize average and indiv reads, both should use maximal clusters for calculation ### TODO: predict based on cluster size -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) MAX_THREADS = None bpp.MAX_THREADS = MAX_THREADS -def scale_reactivities(reactivities): - min = reactivities.min() - max = reactivities.max() - smin = 0 - smax = 2 - reactivities = ((reactivities - min) / (max - min)) * (smax - smin) + smin +# reactivity format for RNAcofold +def scale_reactivities(reactivities): + try: + rmin = float(np.nanmin(reactivities)) + rmax = float(np.nanmax(reactivities)) + smin, smax = 0.0, 2.0 + + reactivities = ((reactivities - rmin) / (rmax - rmin)) * (smax - smin) + smin + except ValueError: # raised if `y` is empty. + pass return reactivities + def get_mods(lids, continue_reads=False, vienna=False): # get sequence data from db df = dbsel.select_predict(lids) df.fillna(value=0, inplace=True) - df = df.groupby(['LID', 'contig', 'read_index', 'position', 'Sequence'], as_index=False)[[ 'Predict_BC', - 'Predict_Signal', 'Predict_Dwell', 'Predict_Lofd', 'Predict_Lofs', - 'Predict_Gmm']].agg(lambda x: x.mode().iloc[0]).reset_index(drop=True) + df = ( + df.groupby( + ["LID", "contig", "read_index", "position", "Sequence"], as_index=False + )[ + [ + "Predict_BC", + "Predict_Signal", + "Predict_Dwell", + "Predict_Lofd", + "Predict_Lofs", + "Predict_Gmm", + ] + ] + .agg(lambda x: x.mode().iloc[0]) + .reset_index(drop=True) + ) if len(df) <= 0: raise Exception("No predictions found. Please run predictions.") return - # -1 is outlier, 1 is inlier - df['Predict_BC'] = np.where(df['Predict_BC'] == -1, 1, 0) - df['Predict_Signal'] = np.where(df['Predict_Signal'] == -1, 1, 0) - df['Predict_Dwell'] = np.where(df['Predict_Dwell'] == -1, 1, 0) - df.loc[(df['Predict_Dwell'] == 1) & (df['Sequence'] != 'G'), "Predict_Dwell"] = 3 - df.loc[(df['Predict_Dwell'] == 1) & (df['Sequence'] == 'G'), "Predict_Dwell"] = 2 - df.loc[df['Predict_Dwell'] == 0, "Predict_Dwell"] = -2 - df['Predict_Lofd'] = np.where(df['Predict_Lofd'] == -1, 1, 0) - df['Predict_Lofs'] = np.where(df['Predict_Lofs'] == -1, 1, 0) - df['Predict_Gmm'] = np.where(df['Predict_Gmm'] == -1, 1, 0) #1,1,0? + df["Predict_BC"] = np.where(df["Predict_BC"] == -1, 1, 0) + df["Predict_Signal"] = np.where(df["Predict_Signal"] == -1, 1, 0) + df["Predict_Dwell"] = np.where(df["Predict_Dwell"] == -1, 1, 0) + df.loc[(df["Predict_Dwell"] == 1) & (df["Sequence"] != "G"), "Predict_Dwell"] = 3 + df.loc[(df["Predict_Dwell"] == 1) & (df["Sequence"] == "G"), "Predict_Dwell"] = 2 + df.loc[df["Predict_Dwell"] == 0, "Predict_Dwell"] = -2 + df["Predict_Lofd"] = np.where(df["Predict_Lofd"] == -1, 1, 0) + df["Predict_Lofs"] = np.where(df["Predict_Lofs"] == -1, 1, 0) + df["Predict_Gmm"] = np.where(df["Predict_Gmm"] == -1, 1, 0) # 1,1,0? df.fillna(value=0, inplace=True) - #print(df.head()) + # print(df.head()) print(len(df)) # calculates average reactivity over all reads @@ -67,16 +84,22 @@ def get_mods(lids, continue_reads=False, vienna=False): # should use maximal clusters for this scoring metric # Least optimal method def mean_read(df): - df['Reactivity'] = df['Predict_BC'] + df['Predict_Signal'] + df['Predict_Dwell'] + \ - df['Predict_Lofs'] + df['Predict_Lofd'] + df['Predict_Gmm'] - #dfr = df[['position', 'contig', 'read_index', 'Sequence', 'Reactivity']] - df['VARNA'] = np.where(df['Reactivity'] >=6, 1, 0) - df['Predict'] = np.where(df['VARNA'] == 1, -1, 1) - #mx.get_Metric(df, seq + "_mean_") + df["Reactivity"] = ( + df["Predict_BC"] + + df["Predict_Signal"] + + df["Predict_Dwell"] + + df["Predict_Lofs"] + + df["Predict_Lofd"] + + df["Predict_Gmm"] + ) + # dfr = df[['position', 'contig', 'read_index', 'Sequence', 'Reactivity']] + df["VARNA"] = np.where(df["Reactivity"] >= 6, 1, 0) + df["Predict"] = np.where(df["VARNA"] == 1, -1, 1) + # mx.get_Metric(df, seq + "_mean_") # TODO add probabilities, opt dbins.insert_rx_full(df) - #mean_read(df) + # mean_read(df) def read_depth(df): print("Read_Depth") @@ -84,72 +107,84 @@ def read_depth(df): #### another decider based on the percentage of modified reads ### ##### Prediction per read, all reads - df['Reactivity'] = df['Predict_BC'] + df['Predict_Signal'] + df['Predict_Dwell'] + \ - df['Predict_Lofs'] + df['Predict_Lofd'] + df['Predict_Gmm'] - df['VARNA'] = np.where(df['Reactivity'] >= 6, 1, 0) - df['Predict'] = np.where(df['VARNA'] == 1, -1, 1) + df["Reactivity"] = ( + df["Predict_BC"] + + df["Predict_Signal"] + + df["Predict_Dwell"] + + df["Predict_Lofs"] + + df["Predict_Lofd"] + + df["Predict_Gmm"] + ) + df["VARNA"] = np.where(df["Reactivity"] >= 6, 1, 0) + df["Predict"] = np.where(df["VARNA"] == 1, -1, 1) # percent reactivities for shape in full set - df['Reactivity_score'] = df['Reactivity']/8 - #adjust read predict based on probabilities - if vienna==True: - if continue_reads==False: + df["Reactivity_score"] = df["Reactivity"] / 8 + # adjust read predict based on probabilities + if vienna == True: + if continue_reads == False: dbins.insert_read_depth_full_clear(df) - dfx = bpp.get_predict_probabilities(df, [lids], .95, reactivity=True, - continue_reads=False) + dfx = bpp.get_predict_probabilities( + df, [lids], 0.95, reactivity=True, continue_reads=False + ) else: - dfx = bpp.get_predict_probabilities(df, [lids], .95, reactivity=True, - continue_reads=True) + dfx = bpp.get_predict_probabilities( + df, [lids], 0.95, reactivity=True, continue_reads=True + ) else: - if continue_reads==False: - #delete + if continue_reads == False: + # delete dbins.insert_read_depth_full_clear(df) - #insert + # insert dbins.insert_read_depth_full(df) else: dbins.insert_read_depth_full(df) - ##### Prediction percent modified and base pairing probability, averaged - dfr = df[['LID','position', 'contig', 'read_index', 'Reactivity']] + dfr = df[["LID", "position", "contig", "read_index", "Reactivity"]] df_rd = dbsel.select_readdepth(lids) - dfr['Predict'] = df['Predict'].astype('category') - dfr = (dfr.groupby(['LID', 'position', 'contig', 'Predict'], observed=False).size(). - unstack(fill_value=0).reset_index()) - dfr = dfr.merge(df_rd, on=['LID', 'position'], how='left') - if -1 not in dfr.columns: #no modifications - dfr[-1]=0 - dfr['percent_modified'] = dfr[-1]/dfr['read_depth'] + dfr["Predict"] = df["Predict"].astype("category") + dfr = ( + dfr.groupby(["LID", "position", "contig", "Predict"], observed=False) + .size() + .unstack(fill_value=0) + .reset_index() + ) + dfr = dfr.merge(df_rd, on=["LID", "position"], how="left") + if -1 not in dfr.columns: # no modifications + dfr[-1] = 0 + dfr["percent_modified"] = dfr[-1] / dfr["read_depth"] # input for RNAfold shape reactivity - dfr['RNAFold_Shape_Reactivity'] = scale_reactivities(dfr['percent_modified']) - mean = dfr['percent_modified'].mean() - std = dfr['percent_modified'].std() - #print("mean ", mean) - + dfr["RNAFold_Shape_Reactivity"] = scale_reactivities(dfr["percent_modified"]) + mean = dfr["percent_modified"].mean() + std = dfr["percent_modified"].std() + # print("mean ", mean) #### Predict modification based on percent modified read depth #### #### Averaged over all transcripts based on read depth ### Most accurate correlates with cluster findings!!!! - #TODO: certainty = read depth/ tot reads, variance = unmod prediction/tot umod reads - #TODO gmm by read and check values - dfr['Predict'] = np.where(dfr['percent_modified'] > mean + std/3, -1, 1) + # TODO: certainty = read depth/ tot reads, variance = unmod prediction/tot umod reads + # TODO gmm by read and check values + dfr["Predict"] = np.where(dfr["percent_modified"] > mean + std / 3, -1, 1) ### adjust prediction with base pairing probabilities dfr.fillna(value=0, inplace=True) - dfr = bpp.get_predict_probabilities(dfr, [lids], .95, reactivity=True, - continue_reads=continue_reads) - dfr.rename(columns={-1: 'Out_num', 1: 'In_num'}, inplace=True) - dfr['VARNA'] = np.where(dfr['Predict'] == -1, 1, 0) - #print(dfr.head()) - #insert to db + dfr = bpp.get_predict_probabilities( + dfr, [lids], 0.95, reactivity=True, continue_reads=continue_reads + ) + dfr.rename(columns={-1: "Out_num", 1: "In_num"}, inplace=True) + dfr["VARNA"] = np.where(dfr["Predict"] == -1, 1, 0) + # print(dfr.head()) + # insert to db dbins.insert_read_depth(dfr) return read_depth(df) + # putative sequences designated in db # sequences = ['RNAse_P', "cen_3'utr", "cen_3'utr_complex", 'cen_FL', 'cen_FL_complex', # "ik2_3'utr_complex", 'ik2_FL_complex', 'T_thermophila', 'ik2_FL', 'HCV', # "ik2_3'utr", 'HSP70_HSPA1A'] ##### test -#get_mods(('43'),continue_reads=False) -#sys.exit(0) +# get_mods(('43'),continue_reads=False) +# sys.exit(0) diff --git a/DashML/Predict/run_predict.py b/DashML/Predict/run_predict.py index 2ffe378a..93d4d5b9 100644 --- a/DashML/Predict/run_predict.py +++ b/DashML/Predict/run_predict.py @@ -1,8 +1,8 @@ import sys, os import numpy as np +from matplotlib.figure import Figure +from matplotlib.backends.backend_agg import FigureCanvasAgg import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt import platform import pandas as pd from platformdirs import user_documents_path @@ -22,29 +22,38 @@ #### Paper GMM positional induced clusters reflect predictions, nice # TODO: add to Predict, + for indicate -pd.set_option('display.max_columns', None) -pd.set_option('display.max_rows', None) +pd.set_option("display.max_columns", None) +pd.set_option("display.max_rows", None) pd.options.mode.chained_assignment = None + def get_lids(seq, temp, type1, type2, complex): ### Source Data #### # TODO: Add modified modified case, generalizing when type1==type2 # generate new lid, duplicated Mod table for that contig with new id in stored proc # take into account multiple runs require diff ids # then run sub predictions normally - if type1==type2: + if type1 == type2: ### predicts on unmodified data, same data sets - acim_sequences = dbsel.select_unmod(seq, temp, type1=type1, type2=type2, complex=complex) - dmso_sequences = dbsel.select_unmod(seq, temp, type1=type1, type2=type1, complex=complex) + acim_sequences = dbsel.select_unmod( + seq, temp, type1=type1, type2=type2, complex=complex + ) + dmso_sequences = dbsel.select_unmod( + seq, temp, type1=type1, type2=type1, complex=complex + ) else: - acim_sequences = dbsel.select_mod(seq, temp, type1=type1, type2=type2, complex=complex) + acim_sequences = dbsel.select_mod( + seq, temp, type1=type1, type2=type2, complex=complex + ) dmso_sequences = dbsel.select_unmod(seq, temp, type1=type1, type2=type1) - unmod_seq_ids = dmso_sequences['LID'].unique() + unmod_seq_ids = dmso_sequences["LID"].unique() unmod_seq_ids = str(unmod_seq_ids).replace("[", "").replace("]", "") - mod_seq_ids = acim_sequences['LID'].unique() + mod_seq_ids = acim_sequences["LID"].unique() mod_seq_ids = str(mod_seq_ids).replace("[", "").replace("]", "") return mod_seq_ids + + ####Predict by Desc ##### def predict(seq, temp, type1, type2, complex): @@ -54,16 +63,22 @@ def predict(seq, temp, type1, type2, complex): # take into account multiple runs require diff ids # then run sub predictions normally - if type1==type2: + if type1 == type2: ### predicts on unmodified data, same data sets - acim_sequences = dbsel.select_unmod(seq, temp, type1=type1, type2=type2, complex=complex) - dmso_sequences = dbsel.select_unmod(seq, temp, type1=type1, type2=type1, complex=complex) + acim_sequences = dbsel.select_unmod( + seq, temp, type1=type1, type2=type2, complex=complex + ) + dmso_sequences = dbsel.select_unmod( + seq, temp, type1=type1, type2=type1, complex=complex + ) else: - acim_sequences = dbsel.select_mod(seq, temp, type1=type1, type2=type2, complex=complex) + acim_sequences = dbsel.select_mod( + seq, temp, type1=type1, type2=type2, complex=complex + ) dmso_sequences = dbsel.select_unmod(seq, temp, type1=type1, type2=type1) - unmod_seq_ids = dmso_sequences['LID'].unique() + unmod_seq_ids = dmso_sequences["LID"].unique() unmod_seq_ids = str(unmod_seq_ids).replace("[", "").replace("]", "") - mod_seq_ids = acim_sequences['LID'].unique() + mod_seq_ids = acim_sequences["LID"].unique() mod_seq_ids = str(mod_seq_ids).replace("[", "").replace("]", "") if (len(acim_sequences) <= 0) or (len(dmso_sequences) <= 0): @@ -83,21 +98,21 @@ def bc_analysis(): return bc.contig = seq bc.get_bc_reactivity_peaks() - bc_analysis() + bc_analysis() #### peak analysis def peak_analysis(): print("Running Peak Anaylysis for " + seq) - if type1==type2: + if type1 == type2: # TODO fix for mod mod peak.acim_sequences = dbsel.select_peaks_unmod(mod_seq_ids, unmod_seq_ids) else: peak.acim_sequences = dbsel.select_peaks(mod_seq_ids, unmod_seq_ids) peak.contig = seq peak.get_reactivity_peaks() - peak_analysis() + peak_analysis() ##### lof def lof_analysis(): @@ -106,8 +121,8 @@ def lof_analysis(): Lof.dmso_sequences = dmso_sequences Lof.seq = seq Lof.get_novelty() - lof_analysis() + lof_analysis() ##### gmm def gmm_analysis(): @@ -117,6 +132,7 @@ def gmm_analysis(): gmx.seq = seq gmx.seq_ids = mod_seq_ids gmx.positional_gmm() + gmm_analysis() return mod_seq_ids @@ -131,19 +147,19 @@ def predict_lids(unmod_lids, mod_lids): acim_sequences = dbsel.select_mod_lid(mod_lids) dmso_sequences = dbsel.select_unmod_lid(unmod_lids) - if (len(acim_sequences) <= 0) or (len(dmso_sequences) <= 0): - seq_err = ("No sequences found. Please load sequence data. modified {} unmodified {}". - format(len(acim_sequences), len(dmso_sequences))) - raise Exception (seq_err) + seq_err = "No sequences found. Please load sequence data. modified {} unmodified {}".format( + len(acim_sequences), len(dmso_sequences) + ) + raise Exception(seq_err) print("Error:", seq_err, file=sys.stderr) return None - unmod_seq_ids = dmso_sequences['LID'].unique() + unmod_seq_ids = dmso_sequences["LID"].unique() unmod_seq_ids = str(unmod_seq_ids).replace("[", "").replace("]", "") - mod_seq_ids = acim_sequences['LID'].unique() + mod_seq_ids = acim_sequences["LID"].unique() mod_seq_ids = str(mod_seq_ids).replace("[", "").replace("]", "") - seq = str(acim_sequences['contig'].unique()[0]).replace("[", "").replace("]", "") + seq = str(acim_sequences["contig"].unique()[0]).replace("[", "").replace("]", "") ### TODO: remove extra parameters from subtasks, and update subtask inserts @@ -158,8 +174,8 @@ def bc_analysis(): sys.exit(0) bc.contig = seq bc.get_bc_reactivity_peaks() - bc_analysis() + bc_analysis() #### peak analysis def peak_analysis(): @@ -167,8 +183,8 @@ def peak_analysis(): peak.acim_sequences = dbsel.select_peaks(mod_seq_ids, unmod_seq_ids) peak.contig = seq peak.get_reactivity_peaks() - peak_analysis() + peak_analysis() ##### lof def lof_analysis(): @@ -177,8 +193,8 @@ def lof_analysis(): Lof.dmso_sequences = dmso_sequences Lof.seq = seq Lof.get_novelty() - lof_analysis() + lof_analysis() ##### gmm def gmm_analysis(): @@ -188,10 +204,12 @@ def gmm_analysis(): gmx.seq = seq gmx.seq_ids = mod_seq_ids gmx.positional_gmm() + gmm_analysis() return mod_seq_ids + #### Predict Reactivity for All Reads and Averaged ##### def get_mods(lids, continue_reads, vienna=False): predict_mods.MAX_THREADS = 100 # vary by platform speed up computation @@ -199,105 +217,143 @@ def get_mods(lids, continue_reads, vienna=False): #### Get Vienna Probabilities for Exisitng Predicitons ##### -def base_prob(lids, threshold=.95, continue_reads=False, gui=None): +def base_prob(lids, threshold=0.95, continue_reads=False, gui=None): predict_probs.MAX_THREADS = 100 # vary by platform speed up computation - predict_probs.get_predict_probability(lids, threshold=threshold, continue_reads=continue_reads) + predict_probs.get_predict_probability( + lids, threshold=threshold, continue_reads=continue_reads + ) #### Get Graphs of Predictions for GUI ##### def get_graph_ave(lids): df = dbsel.select_read_depth_ave(lids) - df.rename(columns={'Rnafold_shape_reactivity': 'Reactivity'}, inplace=True) - seq_name = 'AverageReactivity_' + df['contig'].unique()[0] + '_' + str.replace(str(lids), ',', '-') - filtered_df = df[(df['position'] >= 0) & (df['position'] <= 100)] - ax = filtered_df.plot.bar(x='position', y='Reactivity', legend=False, figsize=(8, 2), width=0.8) - # Set the y-axis minimum - ax.set_ylim(0, 2) - # Set ticks only every 5th position - positions = df['position'].values - tick_indices = np.arange(0, len(positions), 5) # every 5th bar - ax.set_xticks(tick_indices) - ax.set_xticklabels([positions[i] for i in tick_indices], rotation=45, fontsize=6) - ax.tick_params(axis='y', labelsize=6) - ax.set_xlabel('Position', fontsize=8) - ax.set_ylabel('Reactivity', fontsize=8) - ax.set_title('Average Predicted Reactivity', fontsize=9) - subdir = output / "Figures" - subdir.mkdir(parents=True, exist_ok=True) - figname1 = subdir / f"{seq_name}.png" - figname1 = str(figname1.resolve()) - fig = ax.get_figure() - fig.tight_layout() - plt.savefig(figname1) - plt.close(fig) - ##plt.show) + df.rename(columns={"Rnafold_shape_reactivity": "Reactivity"}, inplace=True) + seq_name = ( + "AverageReactivity_" + + df["contig"].unique()[0] + + "_" + + str.replace(str(lids), ",", "-") + ) + filtered_df = df[(df["position"] >= 0) & (df["position"] <= 100)] + + # ---- First figure: Reactivity ---- + with matplotlib.rc_context({"font.size": 8}): + fig1 = Figure(figsize=(8, 2), dpi=100) + FigureCanvasAgg(fig1) + ax1 = fig1.add_subplot(111) + filtered_df.plot.bar( + x="position", + y="Reactivity", + legend=False, + width=0.8, + ax=ax1, + ) + ax1.set_ylim(0, 2) + positions = df["position"].values + tick_indices = np.arange(0, len(positions), 5) # every 5th bar + ax1.set_xticks(tick_indices) + ax1.set_xticklabels( + [positions[i] for i in tick_indices], rotation=45, fontsize=6 + ) + ax1.tick_params(axis="y", labelsize=6) + ax1.set_xlabel("Position", fontsize=8) + ax1.set_ylabel("Reactivity", fontsize=8) + ax1.set_title("Average Predicted Reactivity", fontsize=9) + + subdir = Path(output) / "Figures" + subdir.mkdir(parents=True, exist_ok=True) + figname1 = subdir / f"{seq_name}.png" + fig1.tight_layout() + fig1.savefig(figname1, bbox_inches="tight") + fig1.clear() + del fig1 # print current data subdir = output / "Data" subdir.mkdir(parents=True, exist_ok=True) - f_name = ('AvePredictedReactivity_' + - df['contig'].unique()[0] + '_' + str.replace(str(lids), ',', '-')) + '.csv' + f_name = ( + "AvePredictedReactivity_" + + df["contig"].unique()[0] + + "_" + + str.replace(str(lids), ",", "-") + ) + ".csv" file_path = subdir / f_name file_path = str(file_path.resolve()) - df[['LID', 'contig', 'position', 'Reactivity']].to_csv(file_path, index=False) + df[["LID", "contig", "position", "Reactivity"]].to_csv(file_path, index=False) # graph of base pairing probabilities - df = df.rename(columns={'Base_pair_prob': 'Base Pairing Prob.'}) - filtered_df = filtered_df.rename(columns={'Base_pair_prob': 'Base Pairing Prob.'}) - seq_name = 'AveMaxBasePairingProbability_' + df['contig'].unique()[0] + '_' + str.replace(str(lids), ',', '-') - ax = filtered_df.plot.bar( - x='position', y='Base Pairing Prob.', legend=False, - figsize=(8, 2), width=0.8 + df = df.rename(columns={"Base_pair_prob": "Base Pairing Prob."}) + filtered_df = filtered_df.rename(columns={"Base_pair_prob": "Base Pairing Prob."}) + seq_name = ( + "AveMaxBasePairingProbability_" + + df["contig"].unique()[0] + + "_" + + str.replace(str(lids), ",", "-") ) - # Set the y-axis minimum - ax.set_ylim(0, 1) - # Set ticks only every 5th position - positions = df['position'].values - tick_indices = np.arange(0, len(positions), 5) # every 5th bar - ax.set_xticks(tick_indices) - ax.set_xticklabels([positions[i] for i in tick_indices], rotation=45, fontsize=8) - # Set font size on y-axis ticks - ax.tick_params(axis='y', labelsize=8) - ax.set_xlabel('Position', fontsize=8) - ax.set_ylabel('Probability', fontsize=8) - ax.set_title('Ave. Max Base Pairing Prob. (ViennaRNA)', fontsize=9) - subdir = output / "Figures" - subdir.mkdir(parents=True, exist_ok=True) - figname2 = subdir / f"{seq_name}.png" - figname2 = str(figname2.resolve()) - fig = ax.get_figure() - fig.tight_layout() - plt.savefig(figname2) - plt.close(fig) - # #plt.show) + # ---- Second figure: Probability ---- + with matplotlib.rc_context({"font.size": 8}): + fig2 = Figure(figsize=(8, 2), dpi=100) + FigureCanvasAgg(fig2) + ax2 = fig2.add_subplot(111) + + filtered_df.plot.bar( + x="position", + y="Base Pairing Prob.", + legend=False, + width=0.8, + ax=ax2, + ) + ax2.set_ylim(0, 1) + positions = df["position"].values + tick_indices = np.arange(0, len(positions), 5) # every 5th bar + ax2.set_xticks(tick_indices) + ax2.set_xticklabels( + [positions[i] for i in tick_indices], rotation=45, fontsize=8 + ) + ax2.tick_params(axis="y", labelsize=8) + ax2.set_xlabel("Position", fontsize=8) + ax2.set_ylabel("Probability", fontsize=8) + ax2.set_title("Ave. Max Base Pairing Prob. (ViennaRNA)", fontsize=9) + + subdir = Path(output) / "Figures" + subdir.mkdir(parents=True, exist_ok=True) + figname2 = subdir / f"{seq_name}.png" + fig2.tight_layout() + fig2.savefig(figname2, bbox_inches="tight") + fig2.clear() + del fig2 # print current data subdir = output / "Data" subdir.mkdir(parents=True, exist_ok=True) - f_name = ('AveMaxBasePairingProbability_' + - df['contig'].unique()[0] + '_' + str.replace(str(lids), ',', '-')) + '.csv' + f_name = ( + "AveMaxBasePairingProbability_" + + df["contig"].unique()[0] + + "_" + + str.replace(str(lids), ",", "-") + ) + ".csv" file_path = subdir / f_name file_path = str(file_path.resolve()) - #print(df.columns) - df[['LID', 'contig', 'position', 'Base Pairing Prob.']].to_csv(file_path, index=False) - - return figname1, figname2 - + # print(df.columns) + df[["LID", "contig", "position", "Base Pairing Prob."]].to_csv( + file_path, index=False + ) + return str(figname1.resolve()), str(figname2.resolve()) #### Run Unmodified & Modified Data Separately #### Different Runs of same contig are grouped ### USER DEFINED: LID SELECTED BY DROP DOWN OF POPOPULATED BY DB #### -#seq = 'HCV' #contig in FAFSA/nanopolish -#temp = 37 #37, 42 other -#type1 = 'dmso' #unmodified condition, drop down -#type2 = 'acim' #modified condition, drop down -#complex = 0 #is contig in complex -threshold = .95 -#unmod_lids = 52 -#lids = 37 +# seq = 'HCV' #contig in FAFSA/nanopolish +# temp = 37 #37, 42 other +# type1 = 'dmso' #unmodified condition, drop down +# type2 = 'acim' #modified condition, drop down +# complex = 0 #is contig in complex +threshold = 0.95 +# unmod_lids = 52 +# lids = 37 ### Calculate modifications @@ -305,10 +361,10 @@ def get_graph_ave(lids): def run_predict(unmod_lids, lids, continue_reads=False, vienna=False): try: ### Get LIDs for modification parameters - #mod_lids = get_lids(seq=seq, temp=temp, type1=type1, type2=type2, complex=complex) + # mod_lids = get_lids(seq=seq, temp=temp, type1=type1, type2=type2, complex=complex) #### Run Sub Predictions - #predict(seq=seq, temp=temp, type1=type1, type2=type1, complex=complex) + # predict(seq=seq, temp=temp, type1=type1, type2=type1, complex=complex) predict_lids(unmod_lids, lids) ### Run Combined Predictions by Read and Averaged @@ -323,7 +379,7 @@ def run_predict(unmod_lids, lids, continue_reads=False, vienna=False): ### METRICS ### get metrics on averaged predictions ### - #metric_avg.get_metrics(mod_lids) + # metric_avg.get_metrics(mod_lids) ### get metrics per read ### # TODO add import @@ -335,8 +391,8 @@ def run_predict(unmod_lids, lids, continue_reads=False, vienna=False): print("Predictions Complete") -#run_predict(26,43, continue_reads=False, vienna=True) +# run_predict(26,43, continue_reads=False, vienna=True) ### Calculate modification probabilities over all predictions Vienna (optional) # run if done separetly (default) -#base_prob(lids, threshold, continue_reads=False) +# base_prob(lids, threshold, continue_reads=False)