From 1034af8877f255adf18c1dc3fd59d6295400d417 Mon Sep 17 00:00:00 2001 From: Katzi93_on_mac Date: Tue, 21 Apr 2026 21:17:07 -0400 Subject: [PATCH 1/6] created nupack wrapper and a way to plot the dot bracket notation for dna --- .../sec_struc_pred/Nupack_wrapper.py | 263 +++++++++++ .../simple_binding_example.toml | 16 + .../sec_struc_pred/default_outputs/.gitignore | 2 + .../sec_struc_pred/dot_brac_plotter.py | 427 ++++++++++++++++++ .../sec_struc_pred/scripts/analyse_strands.py | 46 ++ 5 files changed, 754 insertions(+) create mode 100644 crisscross_kit/sec_struc_pred/Nupack_wrapper.py create mode 100644 crisscross_kit/sec_struc_pred/default_configs/simple_binding_example.toml create mode 100644 crisscross_kit/sec_struc_pred/default_outputs/.gitignore create mode 100644 crisscross_kit/sec_struc_pred/dot_brac_plotter.py create mode 100644 crisscross_kit/sec_struc_pred/scripts/analyse_strands.py diff --git a/crisscross_kit/sec_struc_pred/Nupack_wrapper.py b/crisscross_kit/sec_struc_pred/Nupack_wrapper.py new file mode 100644 index 0000000..23fb9fa --- /dev/null +++ b/crisscross_kit/sec_struc_pred/Nupack_wrapper.py @@ -0,0 +1,263 @@ +import tomllib +from pathlib import Path + +from nupack import Model, SetSpec, Strand, Tube, tube_analysis + + +class NupackTubeConfig: + """Small config object for NUPACK nucleic acid tube analysis.""" + + def __init__( + self, + sequences, + material="dna", + celsius=37, + sodium=0.05, + magnesium=0.015, + max_complex_size=2, + compute=None, + options=None, + tube_name="tube_1", + strand_prefix="strand", + ): + self.sequences = sequences + self.material = material + self.celsius = celsius + self.sodium = sodium + self.magnesium = magnesium + self.max_complex_size = max_complex_size + self.compute = compute or ["pfunc", "mfe"] + self.options = options or {} + self.tube_name = tube_name + self.strand_prefix = strand_prefix + + def to_dict(self): + return { + "sequences": self.sequences, + "material": self.material, + "celsius": self.celsius, + "sodium": self.sodium, + "magnesium": self.magnesium, + "max_complex_size": self.max_complex_size, + "compute": self.compute, + "options": self.options, + "tube_name": self.tube_name, + "strand_prefix": self.strand_prefix, + } + + def write(self, path): + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text(self.to_toml()) + return path + + def to_toml(self): + sequences = ",\n ".join( + f'["{sequence}", {concentration_nM}]' + for sequence, concentration_nM in self.sequences + ) + compute = ", ".join(f'"{item}"' for item in self.compute) + options = "\n".join(f"{key} = {value}" for key, value in self.options.items()) + + text = f"""sequences = [ + {sequences} +] + +material = "{self.material}" +celsius = {self.celsius} +sodium = {self.sodium} +magnesium = {self.magnesium} +max_complex_size = {self.max_complex_size} +compute = [{compute}] +tube_name = "{self.tube_name}" +strand_prefix = "{self.strand_prefix}" + +[options] +{options} +""" + return text + + @classmethod + def from_file(cls, path): + with Path(path).open("rb") as f: + return cls(**tomllib.load(f)) + + +def _load_config(config): + if isinstance(config, NupackTubeConfig): + return config.to_dict() + if isinstance(config, (str, Path)): + return NupackTubeConfig.from_file(config).to_dict() + return dict(config) + + +def _complex_to_dict(complex_obj, complex_result, concentration_m, strand_sequences): + minimum_free_energy = None + if complex_result.mfe: + minimum_free_energy = { + "structure": str(complex_result.mfe[0].structure), + "energy_kcal_per_mol": float(complex_result.mfe[0].energy), + "stack_energy_kcal_per_mol": float(complex_result.mfe[0].stack_energy), + } + + return { + "name": complex_obj.name, + "strands": [strand.name for strand in complex_obj.strands], + "sequences": [strand_sequences[strand.name] for strand in complex_obj.strands], + "concentration_nM": float(concentration_m) * 1e9, + "total_free_energy": complex_result.free_energy, + "minimum_free_energy": minimum_free_energy, + } + + +def run_nupack_analysis(config): + config = _load_config(config) + + strands = {} + strand_sequences = {} + for i, (sequence, concentration_nM) in enumerate(config["sequences"], start=1): + strand_name = f"{config.get('strand_prefix', 'strand')}_{i}" + strand = Strand(sequence, name=strand_name) + strands[strand] = concentration_nM * 1e-9 + strand_sequences[strand_name] = sequence + + model = Model( + material=config.get("material", "dna"), + celsius=config.get("celsius", 37), + sodium=config.get("sodium", 0.05), + magnesium=config.get("magnesium", 0.015), + ) + tube = Tube( + strands=strands, + complexes=SetSpec(max_size=config.get("max_complex_size", 2)), + name=config.get("tube_name", "tube_1"), + ) + + analysis_result = tube_analysis( + tubes=[tube], + model=model, + compute=config.get("compute", ["pfunc", "mfe"]), + options=config.get("options", {}), + ) + + tube_result = analysis_result[tube] + complexes = {} + for complex_obj, concentration_m in tube_result.complex_concentrations.items(): + complex_result = analysis_result[complex_obj] + complexes[complex_obj.name] = _complex_to_dict( + complex_obj, + complex_result, + concentration_m, + strand_sequences, + ) + + return { + "config": config, + "direct_config": config, + "fraction_bases_unpaired": tube_result.fraction_bases_unpaired, + "complexes": complexes, + } + + +def _format_optional_number(value, digits=3): + if value is None: + return "NA" + return f"{value:.{digits}f}" + + +def _result_rows(result): + rows = [] + for name, data in result["complexes"].items(): + minimum_free_energy = data["minimum_free_energy"] + minimum_free_energy_value = None + minimum_free_energy_structure = None + if minimum_free_energy is not None: + minimum_free_energy_value = minimum_free_energy["energy_kcal_per_mol"] + minimum_free_energy_structure = minimum_free_energy["structure"] + + rows.append({ + "complex": name, + "sequences": " + ".join(data["sequences"]), + "concentration_nM": data["concentration_nM"], + "total_free_energy": data["total_free_energy"], + "minimum_free_energy": minimum_free_energy_value, + "minimum_free_energy_structure": minimum_free_energy_structure, + }) + return rows + + +def _summary_text(title, result, rows): + lines = [ + "", + title, + "=" * len(title), + f"Fraction bases unpaired: {_format_optional_number(result['fraction_bases_unpaired'])}", + "", + f"{'Complex':20s} " + f"{'Sequences':30s} " + f"{'Conc nM':>12s} " + f"{'Total free energy':>18s} " + f"{'Minimum free energy':>20s} " + f"{'Dot bracket':>25s}", + "-" * 138, + ] + + for row in rows: + lines.append( + f"{row['complex']:20s} " + f"{row['sequences']:30s} " + f"{row['concentration_nM']:12.3f} " + f"{_format_optional_number(row['total_free_energy']):>18s} " + f"{_format_optional_number(row['minimum_free_energy']):>20s} " + f"{str(row['minimum_free_energy_structure'] or 'NA'):>25s}" + ) + + return "\n".join(lines) + + +def _config_rows(config): + rows = [] + for key, value in config.items(): + if key == "sequences": + for i, (sequence, concentration_nM) in enumerate(value, start=1): + rows.append({"parameter": f"sequence_{i}", "value": sequence}) + rows.append({"parameter": f"sequence_{i}_concentration_nM", "value": concentration_nM}) + elif key == "options": + for option_key, option_value in value.items(): + rows.append({"parameter": f"options.{option_key}", "value": option_value}) + else: + rows.append({"parameter": key, "value": value}) + return rows + + +def _save_binding_summary(save, result_rows, config_rows): + import pandas as pd + + excel_path = Path(save) + excel_path.parent.mkdir(parents=True, exist_ok=True) + + with pd.ExcelWriter(excel_path) as writer: + pd.DataFrame(result_rows).to_excel(writer, sheet_name="results", index=False) + pd.DataFrame(config_rows).to_excel(writer, sheet_name="config", index=False) + + +def print_binding_summary(title, result, save=""): + rows = _result_rows(result) + text = _summary_text(title, result, rows) + print(text) + + if save: + config = result.get("direct_config", result["config"]) + output_path = Path(save) / f"{title}.xlsx" + _save_binding_summary(output_path, rows, _config_rows(config)) + + +if __name__ == "__main__": + test_config = NupackTubeConfig( + sequences=[ + ["GCGTATGC", 1000], + ["GCATACGC", 1000], + ], + ) + result = run_nupack_analysis(test_config) + print_binding_summary("NUPACK analysis example", result) diff --git a/crisscross_kit/sec_struc_pred/default_configs/simple_binding_example.toml b/crisscross_kit/sec_struc_pred/default_configs/simple_binding_example.toml new file mode 100644 index 0000000..ed1a833 --- /dev/null +++ b/crisscross_kit/sec_struc_pred/default_configs/simple_binding_example.toml @@ -0,0 +1,16 @@ +sequences = [ + ["GCGTATGC", 1000], + ["GCATACGC", 1000], + ["TTTTTTTT", 100] +] + +material = "dna" +celsius = 37 +sodium = 0.05 +magnesium = 0.015 +max_complex_size = 2 +compute = ["pfunc", "mfe"] +tube_name = "simple_binding_example" +strand_prefix = "seq" + +[options] diff --git a/crisscross_kit/sec_struc_pred/default_outputs/.gitignore b/crisscross_kit/sec_struc_pred/default_outputs/.gitignore new file mode 100644 index 0000000..d6b7ef3 --- /dev/null +++ b/crisscross_kit/sec_struc_pred/default_outputs/.gitignore @@ -0,0 +1,2 @@ +* +!.gitignore diff --git a/crisscross_kit/sec_struc_pred/dot_brac_plotter.py b/crisscross_kit/sec_struc_pred/dot_brac_plotter.py new file mode 100644 index 0000000..fe47f25 --- /dev/null +++ b/crisscross_kit/sec_struc_pred/dot_brac_plotter.py @@ -0,0 +1,427 @@ +from pathlib import Path +from math import cos, pi, sin +from random import Random +from xml.sax.saxutils import escape + +import numpy as np + + +BASE_SIMULATION_PARAMS = { + "backbone_distance": 34, + "pair_distance_factor": 1.5, + "start_radius": 95, + "reference_sequence_length": 16, + "relaxation_steps": 15000, + "step_size": 1.25, + "start_jitter": 8, + "backbone_spring": 0.135, + "pair_spring": 0.175, + "backbone_straightness": 0.268, + "centering": 0.002, + "neighbor_repulsion": 40, + "base_repulsion": 1500, + "scale": 0.5, + "row_padding": 1000, + "left_padding": 650, + "text_gap": 600, + "brownian_jitter": 0, + "text_line_width": 105, +} +DEFAULT_SIMULATION_PARAMS = dict(BASE_SIMULATION_PARAMS) + + +def set_simulation_params(**params): + unknown_params = set(params) - set(DEFAULT_SIMULATION_PARAMS) + if unknown_params: + raise KeyError(f"Unknown simulation params: {sorted(unknown_params)}") + DEFAULT_SIMULATION_PARAMS.update(params) + + +def set_simulation_param(name, value): + set_simulation_params(**{name: value}) + + +def reset_simulation_params(): + DEFAULT_SIMULATION_PARAMS.clear() + DEFAULT_SIMULATION_PARAMS.update(BASE_SIMULATION_PARAMS) + + +def get_simulation_params(): + return dict(DEFAULT_SIMULATION_PARAMS) + + +def parse_dot_bracket(dot_bracket): + pairs = [] + stack = [] + base_index = 0 + + for char in dot_bracket: + if char == "+": + continue + if char == "(": + stack.append(base_index) + elif char == ")": + pairs.append((stack.pop(), base_index)) + base_index += 1 + + return pairs + + +def _simulation_params(simulation_params=None): + params = dict(DEFAULT_SIMULATION_PARAMS) + if simulation_params: + params.update(simulation_params) + params["pair_distance"] = params["backbone_distance"] * params["pair_distance_factor"] + return params + + +def _initial_radius(sequence_length, params): + return params["start_radius"] * sequence_length / params["reference_sequence_length"] + + +def _complex_radius(complex_data, params): + sequence_length = sum(len(sequence) for sequence in complex_data["sequences"]) + return _initial_radius(sequence_length, params) + + +def _circle_positions(sequences, center_x, center_y, radius, params): + sequence = "".join(sequences) + positions = [] + labels = [] + n = len(sequence) + rng = Random(1) + + for i, base in enumerate(sequence): + angle = -pi / 2 + 2 * pi * i / n + x = center_x + radius * cos(angle) + rng.uniform(-params["start_jitter"], params["start_jitter"]) + y = center_y + radius * sin(angle) + rng.uniform(-params["start_jitter"], params["start_jitter"]) + positions.append((x, y)) + labels.append(base) + + return positions, labels + + +def _add_springs(forces, positions, edges, target_distance, strength): + if len(edges) == 0: + return + + left = edges[:, 0] + right = edges[:, 1] + delta = positions[right] - positions[left] + distance = np.linalg.norm(delta, axis=1) + distance[distance == 0] = 0.001 + force = strength * (distance - target_distance) + force_vectors = delta * (force / distance)[:, None] + + np.add.at(forces, left, force_vectors) + np.add.at(forces, right, -force_vectors) + + +def _backbone_terms(sequences): + edges = [] + triples = [] + start = 0 + + for sequence in sequences: + for i in range(start, start + len(sequence) - 1): + edges.append((i, i + 1)) + for i in range(start + 1, start + len(sequence) - 1): + triples.append((i - 1, i, i + 1)) + start += len(sequence) + + return np.array(edges, dtype=int), np.array(triples, dtype=int) + + +def _add_straightness(forces, positions, triples, strength): + if len(triples) == 0: + return + + left = triples[:, 0] + mid = triples[:, 1] + right = triples[:, 2] + target = (positions[left] + positions[right]) / 2 + force_vectors = (target - positions[mid]) * strength + + np.add.at(forces, mid, force_vectors) + np.add.at(forces, left, -force_vectors * 0.5) + np.add.at(forces, right, -force_vectors * 0.5) + + +def _repulsion_matrix(num_positions, backbone_edges, params): + repulsion = np.full((num_positions, num_positions), params["base_repulsion"], dtype=float) + np.fill_diagonal(repulsion, 0) + if len(backbone_edges): + left = backbone_edges[:, 0] + right = backbone_edges[:, 1] + repulsion[left, right] = params["neighbor_repulsion"] + repulsion[right, left] = params["neighbor_repulsion"] + return repulsion + + +def _add_repulsion(forces, positions, repulsion): + dx = positions[:, 0][None, :] - positions[:, 0][:, None] + dy = positions[:, 1][None, :] - positions[:, 1][:, None] + distance_squared = dx * dx + dy * dy + distance_squared[distance_squared == 0] = 0.001 + + force = repulsion / distance_squared + distance = np.sqrt(distance_squared) + fx = force * dx / distance + fy = force * dy / distance + + forces[:, 0] -= fx.sum(axis=1) + forces[:, 1] -= fy.sum(axis=1) + + +def _relax_positions(sequences, pairs, center_x, center_y, params): + sequence_length = sum(len(sequence) for sequence in sequences) + radius = _initial_radius(sequence_length, params) + positions, labels = _circle_positions(sequences, center_x, center_y, radius, params) + positions = np.array(positions, dtype=float) + backbone_edges, backbone_triples = _backbone_terms(sequences) + pair_edges = np.array(pairs, dtype=int) + repulsion = _repulsion_matrix(len(positions), backbone_edges, params) + rng = np.random.default_rng(2) + center = np.array([center_x, center_y]) + + for _ in range(params["relaxation_steps"]): + forces = np.zeros_like(positions) + + _add_springs(forces, positions, backbone_edges, params["backbone_distance"], params["backbone_spring"]) + _add_straightness(forces, positions, backbone_triples, params["backbone_straightness"]) + _add_springs(forces, positions, pair_edges, params["pair_distance"], params["pair_spring"]) + _add_repulsion(forces, positions, repulsion) + forces += (center - positions) * params["centering"] + forces += rng.uniform(-params["brownian_jitter"], params["brownian_jitter"], size=forces.shape) + + movement = np.clip(forces * params["step_size"], -4, 4) + positions += movement + + return positions.tolist(), labels + + +def _backbone_point_sets(sequences, positions): + point_sets = [] + start = 0 + for sequence in sequences: + end = start + len(sequence) + point_sets.append(positions[start:end]) + start = end + return point_sets + + +def _complex_metadata_lines(complex_data): + lines = [] + if "total_free_energy" in complex_data: + lines.append(f"Total free energy: {complex_data['total_free_energy']:.4f} kcal/mol") + if "minimum_free_energy" in complex_data: + lines.append(f"Minimum free energy: {complex_data['minimum_free_energy']:.4f} kcal/mol") + if "concentration_nM" in complex_data: + lines.append(f"Concentration: {complex_data['concentration_nM']:.4f} nM") + return lines + + +def _wrap_text(text, max_chars): + return [text[i:i + max_chars] for i in range(0, len(text), max_chars)] or [""] + + +def _complex_text_lines(complex_data, params): + lines = [complex_data["title"]] + for i, sequence in enumerate(complex_data["sequences"], start=1): + for line in _wrap_text(sequence, params["text_line_width"]): + lines.append(f"seq_{i}: {line}") + for line in _wrap_text(complex_data["dot_bracket"], params["text_line_width"]): + lines.append(f"dot: {line}") + lines.extend(_complex_metadata_lines(complex_data)) + return lines + + +def _longest_text_line(complexes, params): + return max(len(line) for complex_data in complexes for line in _complex_text_lines(complex_data, params)) + + +def _base_class(base): + return f"base-circle base-circle-{base.upper()}" + + +def _complex_svg(complex_data, center_x, center_y, radius, params): + title = complex_data["title"] + sequences = complex_data["sequences"] + dot_bracket = complex_data["dot_bracket"] + pairs = parse_dot_bracket(dot_bracket) + positions, labels = _relax_positions(sequences, pairs, center_x, center_y, params) + if len(labels) != len(dot_bracket.replace("+", "")): + raise ValueError(f"{title}: sequence length does not match dot-bracket length") + + text_x = center_x + radius + params["text_gap"] + text_y = center_y - 85 + svg = [] + for i, line in enumerate(_complex_text_lines(complex_data, params)): + class_name = "title" if i == 0 else "subtle" + svg.append(f'{escape(line)}') + + for points in _backbone_point_sets(sequences, positions): + path = " ".join( + f"{'M' if i == 0 else 'L'} {x:.1f} {y:.1f}" + for i, (x, y) in enumerate(points) + ) + svg.append(f'') + + for left_i, right_i in pairs: + x1, y1 = positions[left_i] + x2, y2 = positions[right_i] + svg.append(f'') + + for (x, base_y), base in zip(positions, labels): + svg.append(f'') + + for (x, base_y), base in zip(positions, labels): + svg.append(f'{escape(base)}') + + return svg + + +def _write_plot_output(svg_text, output_path, file_format): + output_path = Path(output_path) + output_path.parent.mkdir(parents=True, exist_ok=True) + + if file_format == "svg": + output_path.write_text(svg_text) + elif file_format == "pdf": + try: + import cairosvg + cairosvg.svg2pdf(bytestring=svg_text.encode(), write_to=str(output_path)) + except (ImportError, OSError) as exc: + raise RuntimeError( + "PDF output needs cairosvg and the Cairo system library. " + "In this conda environment, try: conda install -c conda-forge cairo cairosvg. " + "Or use file_format='svg'." + ) from exc + else: + raise ValueError("file_format must be 'svg' or 'pdf'") + + return output_path + + +def plot_dot_bracket_complexes(complexes, output_path, simulation_params=None, file_format="svg"): + params = _simulation_params(simulation_params) + row_padding = params["row_padding"] + max_radius = max(_complex_radius(complex_data, params) for complex_data in complexes) + row_heights = [2 * _complex_radius(complex_data, params) + row_padding for complex_data in complexes] + text_width = _longest_text_line(complexes, params) * 8 + view_width = int(max(900, 2 * max_radius + params["left_padding"] + params["text_gap"] + text_width + 80)) + width = int(view_width * params["scale"]) + height = int((sum(row_heights) + 60) * params["scale"]) + view_height = int(sum(row_heights) + 60) + + svg = [ + f'', + "", + ] + + y = 30 + for complex_data, row_height in zip(complexes, row_heights): + radius = _complex_radius(complex_data, params) + center_x = max_radius + params["left_padding"] + center_y = y + radius + row_padding / 2 + svg.extend(_complex_svg( + complex_data, + center_x, + center_y, + radius, + params, + )) + y += row_height + + svg.append("") + + return _write_plot_output("\n".join(svg), output_path, file_format) + + +def _safe_filename(title): + return "".join(char if char.isalnum() or char in "._- " else "_" for char in title).strip() + + +def _nupack_complex_to_plot_data(name, complex_data): + minimum_free_energy = complex_data.get("minimum_free_energy") + if not minimum_free_energy: + return None + + return { + "title": name, + "sequences": complex_data["sequences"], + "dot_bracket": minimum_free_energy["structure"], + "total_free_energy": complex_data.get("total_free_energy"), + "minimum_free_energy": minimum_free_energy.get("energy_kcal_per_mol"), + "concentration_nM": complex_data.get("concentration_nM"), + } + + +def plot_nupack_result(title, result, save="", simulation_params=None, file_format="svg"): + complexes = [] + for name, complex_data in result["complexes"].items(): + plot_data = _nupack_complex_to_plot_data(name, complex_data) + if plot_data: + complexes.append(plot_data) + + if not save: + save = Path(__file__).parent / "default_outputs" + + output_path = Path(save) / f"{_safe_filename(title)}.{file_format}" + return plot_dot_bracket_complexes( + complexes, + output_path, + simulation_params=simulation_params, + file_format=file_format, + ) + + +if __name__ == "__main__": + examples = [ + { + "title": "(seq_1)", + "sequences": ["CATATCCGCGTCGCTGCGCTCAGACCCACCACCACGCACC"], + "dot_bracket": ".......((...))((((................))))..", + "total_free_energy": -3.869, + "minimum_free_energy": -2.664, + "concentration_nM": 967.997, + }, + { + "title": "(seq_1+seq_1 long)", + "sequences": [ + "CATATCCGCGTCGCTGCGCTCAGACCCACCACCACGCACC", + "CATATCCGCGTCGCTGCGCTCAGACCCACCACCACGCACC", + ], + "dot_bracket": ".......(((.(((((((................))))..+.......))).)))((((................))))..", + "total_free_energy": -16.215, + "minimum_free_energy": -16.099, + "concentration_nM": 16.002, + }, + { + "title": "(seq_1+seq_1)", + "sequences": ["GCGTATGC", "GCGTATGC"], + "dot_bracket": "((.((.((+)).)).))", + }, + { + "title": "(seq_2+seq_3)", + "sequences": ["GCATACGC", "TTTTTTTT"], + "dot_bracket": "....(...+.)......", + }, + ] + + output = Path(__file__).parent / "default_outputs" / "test.svg" + plot_dot_bracket_complexes(examples, output) + print(f"Wrote {output}") diff --git a/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py b/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py new file mode 100644 index 0000000..5dad6f1 --- /dev/null +++ b/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py @@ -0,0 +1,46 @@ +from pathlib import Path +import sys + +sys.path.append(str(Path(__file__).resolve().parents[1])) + +from Nupack_wrapper import NupackTubeConfig, print_binding_summary, run_nupack_analysis +from dot_brac_plotter import plot_nupack_result, set_simulation_params + + +title = "test1" +standard_out = Path(__file__).resolve().parents[1] / "default_outputs" + +direct_config = NupackTubeConfig( + sequences=[ + ["TTTAAAATTTACTGGGCGCTGCAAGTTTTTTTTTTT", 1000], + ["ACTTGCAGGCCCAGTTTTTTTAAAAAAAAAAAAA", 1000], + ], + material="dna", + celsius=37, + sodium=0.05, + magnesium=0.015, + max_complex_size=2, + compute=["pfunc", "mfe"], + tube_name=title, + strand_prefix="seq", +) + +set_simulation_params(brownian_jitter=0.0) +set_simulation_params(relaxation_steps= 55000) + +direct_result = run_nupack_analysis(direct_config) +print_binding_summary(title, direct_result, save=standard_out) +plot_nupack_result( + title, + direct_result, + save=standard_out, + file_format="pdf", # file_format="svg", +) + +# config_path = Path(__file__).resolve().parents[1] / "default_configs" / "simple_binding_example.toml" +# direct_config.write(config_path) +# +# file_result = run_nupack_analysis(config_path) +# print_binding_summary("Same analysis loaded from TOML", file_result) +# print() +# print(f"Wrote TOML config to: {config_path}") From 310fe5e96accb00abd970e8856cf28ba5499099e Mon Sep 17 00:00:00 2001 From: Katzi93_on_mac Date: Wed, 22 Apr 2026 09:27:10 -0400 Subject: [PATCH 2/6] added now module to pyproject.toml. added init file for new model. added 5' and 3' makers to the plot --- crisscross_kit/pyproject.toml | 1 + crisscross_kit/sec_struc_pred/__init__.py | 1 + .../sec_struc_pred/dot_brac_plotter.py | 35 +++++++++++++++++++ .../sec_struc_pred/scripts/analyse_strands.py | 13 +++---- 4 files changed, 42 insertions(+), 8 deletions(-) create mode 100644 crisscross_kit/sec_struc_pred/__init__.py diff --git a/crisscross_kit/pyproject.toml b/crisscross_kit/pyproject.toml index e400b88..42ecca2 100644 --- a/crisscross_kit/pyproject.toml +++ b/crisscross_kit/pyproject.toml @@ -63,6 +63,7 @@ packages = [ "orthoseq_generator.streamlit_app", "orthoseq_generator.streamlit_app.tabs", "eqcorr2d", + "sec_struc_pred", ] include-package-data = true diff --git a/crisscross_kit/sec_struc_pred/__init__.py b/crisscross_kit/sec_struc_pred/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/crisscross_kit/sec_struc_pred/__init__.py @@ -0,0 +1 @@ + diff --git a/crisscross_kit/sec_struc_pred/dot_brac_plotter.py b/crisscross_kit/sec_struc_pred/dot_brac_plotter.py index fe47f25..44e8309 100644 --- a/crisscross_kit/sec_struc_pred/dot_brac_plotter.py +++ b/crisscross_kit/sec_struc_pred/dot_brac_plotter.py @@ -26,6 +26,7 @@ "text_gap": 600, "brownian_jitter": 0, "text_line_width": 105, + "terminal_circle_radius": 6, } DEFAULT_SIMULATION_PARAMS = dict(BASE_SIMULATION_PARAMS) @@ -210,6 +211,16 @@ def _backbone_point_sets(sequences, positions): return point_sets +def _strand_terminal_indices(sequences): + indices = [] + start = 0 + for sequence in sequences: + end = start + len(sequence) - 1 + indices.append((start, end)) + start = end + 1 + return indices + + def _complex_metadata_lines(complex_data): lines = [] if "total_free_energy" in complex_data: @@ -272,6 +283,25 @@ def _complex_svg(complex_data, center_x, center_y, radius, params): x2, y2 = positions[right_i] svg.append(f'') + terminal_bond_length = params["backbone_distance"] / 2 * 1.333 + for start_i, end_i in _strand_terminal_indices(sequences): + terminal_data = [ + (start_i, start_i + 1, "5'", "five-prime"), + (end_i, end_i - 1, "3'", "three-prime"), + ] + for index, neighbor_i, label, end_class in terminal_data: + x, y = positions[index] + neighbor_x, neighbor_y = positions[neighbor_i] + dx = x - neighbor_x + dy = y - neighbor_y + distance = (dx * dx + dy * dy) ** 0.5 or 1 + label_x = x + terminal_bond_length * dx / distance + label_y = y + terminal_bond_length * dy / distance + svg.append(f'') + terminal_radius = params["terminal_circle_radius"] + svg.append(f'') + svg.append(f'{label}') + for (x, base_y), base in zip(positions, labels): svg.append(f'') @@ -321,7 +351,12 @@ def plot_dot_bracket_complexes(complexes, output_path, simulation_params=None, f ".title { font-size: 18px; font-weight: 700; fill: #222; }", ".subtle { font-size: 13px; fill: #555; }", ".base { font-size: 13px; font-weight: 700; text-anchor: middle; fill: #111; }", + ".end-label { font-size: 9px; font-weight: 700; text-anchor: middle; fill: #333; }", ".base-circle { stroke: #222; stroke-width: 1.4; }", + ".terminal-circle { stroke: #666; stroke-width: 1.0; }", + ".terminal-five-prime { fill: #008b8b; }", + ".terminal-three-prime { fill: #ff7f11; }", + ".terminal-bond { stroke: #999; stroke-width: 1.5; stroke-linecap: round; }", ".base-circle-A { fill: #9bbcff; }", ".base-circle-C { fill: #91d7a7; }", ".base-circle-G { fill: #c1a7e8; }", diff --git a/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py b/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py index 5dad6f1..6c3f1e8 100644 --- a/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py +++ b/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py @@ -1,13 +1,10 @@ from pathlib import Path -import sys -sys.path.append(str(Path(__file__).resolve().parents[1])) +from sec_struc_pred.nupack_wrapper import NupackTubeConfig, print_binding_summary, run_nupack_analysis +from sec_struc_pred.dot_brac_plotter import plot_nupack_result, set_simulation_params -from Nupack_wrapper import NupackTubeConfig, print_binding_summary, run_nupack_analysis -from dot_brac_plotter import plot_nupack_result, set_simulation_params - -title = "test1" +title = "test12" standard_out = Path(__file__).resolve().parents[1] / "default_outputs" direct_config = NupackTubeConfig( @@ -24,9 +21,9 @@ tube_name=title, strand_prefix="seq", ) - +# parameters for the plotter. We use a relaxation algorithm. Its a toy version of a force field simulation set_simulation_params(brownian_jitter=0.0) -set_simulation_params(relaxation_steps= 55000) +set_simulation_params(relaxation_steps= 15000) direct_result = run_nupack_analysis(direct_config) print_binding_summary(title, direct_result, save=standard_out) From 190336e9469ffce143af7432d6b590c172eb455b Mon Sep 17 00:00:00 2001 From: Katzi93_on_mac Date: Wed, 22 Apr 2026 10:04:22 -0400 Subject: [PATCH 3/6] slight parameter adjustments for plotting --- crisscross_kit/sec_struc_pred/dot_brac_plotter.py | 10 +++++----- .../sec_struc_pred/scripts/analyse_strands.py | 6 ++++-- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/crisscross_kit/sec_struc_pred/dot_brac_plotter.py b/crisscross_kit/sec_struc_pred/dot_brac_plotter.py index 44e8309..b8ccd67 100644 --- a/crisscross_kit/sec_struc_pred/dot_brac_plotter.py +++ b/crisscross_kit/sec_struc_pred/dot_brac_plotter.py @@ -12,12 +12,12 @@ "start_radius": 95, "reference_sequence_length": 16, "relaxation_steps": 15000, - "step_size": 1.25, + "step_size": 1.00, "start_jitter": 8, - "backbone_spring": 0.135, - "pair_spring": 0.175, - "backbone_straightness": 0.268, - "centering": 0.002, + "backbone_spring": 0.035, + "pair_spring": 0.15, + "backbone_straightness": 0.148, + "centering": 0.0006, "neighbor_repulsion": 40, "base_repulsion": 1500, "scale": 0.5, diff --git a/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py b/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py index 6c3f1e8..dee5c0f 100644 --- a/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py +++ b/crisscross_kit/sec_struc_pred/scripts/analyse_strands.py @@ -9,8 +9,8 @@ direct_config = NupackTubeConfig( sequences=[ - ["TTTAAAATTTACTGGGCGCTGCAAGTTTTTTTTTTT", 1000], - ["ACTTGCAGGCCCAGTTTTTTTAAAAAAAAAAAAA", 1000], + ["TTTTTTTTTTTTAAAAAGGGGGCAGGCTTTTTTTTTTTTTTT", 1000], + ["AAAAAAAAAAAAAAAAAAAAAAAAAAACCCCTACTCGAAATTT", 1000], ], material="dna", celsius=37, @@ -24,6 +24,8 @@ # parameters for the plotter. We use a relaxation algorithm. Its a toy version of a force field simulation set_simulation_params(brownian_jitter=0.0) set_simulation_params(relaxation_steps= 15000) +set_simulation_params(centering=0.0006) +#set_simulation_params(reference_sequence_length=300)#"start_radius": 95,"reference_sequence_length": 16, direct_result = run_nupack_analysis(direct_config) print_binding_summary(title, direct_result, save=standard_out) From 5e1e3185f0d2e74b051022509ff7f2aaec64fcff Mon Sep 17 00:00:00 2001 From: Katzi93_on_mac Date: Wed, 22 Apr 2026 10:24:05 -0400 Subject: [PATCH 4/6] wrong file name ..... --- .../sec_struc_pred/{Nupack_wrapper.py => nupack_wrapper.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename crisscross_kit/sec_struc_pred/{Nupack_wrapper.py => nupack_wrapper.py} (100%) diff --git a/crisscross_kit/sec_struc_pred/Nupack_wrapper.py b/crisscross_kit/sec_struc_pred/nupack_wrapper.py similarity index 100% rename from crisscross_kit/sec_struc_pred/Nupack_wrapper.py rename to crisscross_kit/sec_struc_pred/nupack_wrapper.py From ab010c6055062c249a50c613f32b787a2518d66c Mon Sep 17 00:00:00 2001 From: matt Date: Wed, 22 Apr 2026 23:41:58 -0400 Subject: [PATCH 5/6] Re-formatted nupack wrapper to match repo spec and started regularizing graphics generated for the svg nupack plotters Also removed Cairo dependency --- .../sec_struc_pred/dot_brac_plotter.py | 253 ++++++++++++------ .../sec_struc_pred/nupack_wrapper.py | 190 ++++--------- 2 files changed, 227 insertions(+), 216 deletions(-) diff --git a/crisscross_kit/sec_struc_pred/dot_brac_plotter.py b/crisscross_kit/sec_struc_pred/dot_brac_plotter.py index b8ccd67..15911ab 100644 --- a/crisscross_kit/sec_struc_pred/dot_brac_plotter.py +++ b/crisscross_kit/sec_struc_pred/dot_brac_plotter.py @@ -20,10 +20,10 @@ "centering": 0.0006, "neighbor_repulsion": 40, "base_repulsion": 1500, - "scale": 0.5, - "row_padding": 1000, - "left_padding": 650, - "text_gap": 600, + "scale": 0.85, + "row_padding": 200, + "left_padding": 80, + "text_gap": 120, "brownian_jitter": 0, "text_line_width": 105, "terminal_circle_radius": 6, @@ -255,7 +255,12 @@ def _base_class(base): return f"base-circle base-circle-{base.upper()}" -def _complex_svg(complex_data, center_x, center_y, radius, params): +BASE_COLORS = {"A": "#9bbcff", "C": "#91d7a7", "G": "#c1a7e8", "T": "#d6b06a", "U": "#d6b06a"} +TERMINAL_COLORS = {"five-prime": "#008b8b", "three-prime": "#ff7f11"} + + +def _compute_complex_layout(complex_data, center_x, center_y, radius, params): + """Compute positions and drawing elements for a single complex (format-agnostic).""" title = complex_data["title"] sequences = complex_data["sequences"] dot_bracket = complex_data["dot_bracket"] @@ -264,94 +269,179 @@ def _complex_svg(complex_data, center_x, center_y, radius, params): if len(labels) != len(dot_bracket.replace("+", "")): raise ValueError(f"{title}: sequence length does not match dot-bracket length") - text_x = center_x + radius + params["text_gap"] - text_y = center_y - 85 - svg = [] - for i, line in enumerate(_complex_text_lines(complex_data, params)): - class_name = "title" if i == 0 else "subtle" - svg.append(f'{escape(line)}') - - for points in _backbone_point_sets(sequences, positions): - path = " ".join( - f"{'M' if i == 0 else 'L'} {x:.1f} {y:.1f}" - for i, (x, y) in enumerate(points) - ) - svg.append(f'') - - for left_i, right_i in pairs: - x1, y1 = positions[left_i] - x2, y2 = positions[right_i] - svg.append(f'') + text_lines = _complex_text_lines(complex_data, params) + backbone_sets = _backbone_point_sets(sequences, positions) terminal_bond_length = params["backbone_distance"] / 2 * 1.333 + terminals = [] for start_i, end_i in _strand_terminal_indices(sequences): - terminal_data = [ - (start_i, start_i + 1, "5'", "five-prime"), - (end_i, end_i - 1, "3'", "three-prime"), - ] - for index, neighbor_i, label, end_class in terminal_data: + for index, neighbor_i, label, end_class in [ + (start_i, start_i + 1, "5'", "five-prime"), (end_i, end_i - 1, "3'", "three-prime"), + ]: x, y = positions[index] - neighbor_x, neighbor_y = positions[neighbor_i] - dx = x - neighbor_x - dy = y - neighbor_y - distance = (dx * dx + dy * dy) ** 0.5 or 1 - label_x = x + terminal_bond_length * dx / distance - label_y = y + terminal_bond_length * dy / distance - svg.append(f'') - terminal_radius = params["terminal_circle_radius"] - svg.append(f'') - svg.append(f'{label}') - - for (x, base_y), base in zip(positions, labels): - svg.append(f'') - - for (x, base_y), base in zip(positions, labels): - svg.append(f'{escape(base)}') + nx, ny = positions[neighbor_i] + dx, dy = x - nx, y - ny + dist = (dx * dx + dy * dy) ** 0.5 or 1 + lx = x + terminal_bond_length * dx / dist + ly = y + terminal_bond_length * dy / dist + terminals.append({"x": x, "y": y, "lx": lx, "ly": ly, "label": label, "end_class": end_class, + "radius": params["terminal_circle_radius"]}) + + return {"positions": positions, "labels": labels, "pairs": pairs, "radius": radius, + "text_lines": text_lines, "backbone_sets": backbone_sets, "terminals": terminals} + + +def _complex_svg(layout, text_x, text_y, struct_offset_x): + """Render a complex layout to SVG elements. Text on left, structure shifted right by struct_offset_x.""" + svg = [] + for i, line in enumerate(layout["text_lines"]): + cls = "title" if i == 0 else "subtle" + svg.append(f'{escape(line)}') + + for points in layout["backbone_sets"]: + path = " ".join(f"{'M' if i == 0 else 'L'} {x + struct_offset_x:.1f} {y:.1f}" + for i, (x, y) in enumerate(points)) + svg.append(f'') + + for left_i, right_i in layout["pairs"]: + x1, y1 = layout["positions"][left_i] + x2, y2 = layout["positions"][right_i] + svg.append(f'') + + for t in layout["terminals"]: + svg.append(f'') + svg.append(f'') + svg.append(f'{t["label"]}') + + for (x, y), base in zip(layout["positions"], layout["labels"]): + svg.append(f'') + + for (x, y), base in zip(layout["positions"], layout["labels"]): + svg.append(f'{escape(base)}') return svg -def _write_plot_output(svg_text, output_path, file_format): - output_path = Path(output_path) - output_path.parent.mkdir(parents=True, exist_ok=True) +def _render_complexes_to_matplotlib(layouts, text_col_width, view_width, view_height, output_path, params): + """Render complex layouts directly to PDF using matplotlib (no system dependencies needed).""" + import matplotlib.pyplot as plt + from matplotlib.patches import Circle as MplCircle - if file_format == "svg": - output_path.write_text(svg_text) - elif file_format == "pdf": - try: - import cairosvg - cairosvg.svg2pdf(bytestring=svg_text.encode(), write_to=str(output_path)) - except (ImportError, OSError) as exc: - raise RuntimeError( - "PDF output needs cairosvg and the Cairo system library. " - "In this conda environment, try: conda install -c conda-forge cairo cairosvg. " - "Or use file_format='svg'." - ) from exc - else: - raise ValueError("file_format must be 'svg' or 'pdf'") + fig_w = view_width * params["scale"] / 72 + fig_h = view_height * params["scale"] / 72 + fig, ax = plt.subplots(figsize=(fig_w, fig_h)) + ax.set_xlim(0, view_width) + ax.set_ylim(view_height, 0) + ax.set_aspect("equal") + ax.axis("off") + + row_padding = params["row_padding"] + text_x = params["left_padding"] + struct_center_x = text_col_width + params["text_gap"] + + y = 20 + for layout in layouts: + radius = layout["radius"] + row_height = 2 * radius + row_padding + center_y = y + radius + row_padding / 2 + text_y = center_y - 50 + + for i, line in enumerate(layout["text_lines"]): + fontsize = 11 if i == 0 else 8 + weight = "bold" if i == 0 else "normal" + color = "#222" if i == 0 else "#555" + ax.text(text_x, text_y + i * 16, line, + fontsize=fontsize, fontweight=weight, color=color, fontfamily="sans-serif", va="top") + + struct_offset_x = struct_center_x + for points in layout["backbone_sets"]: + xs = [p[0] + struct_offset_x for p in points] + ys = [p[1] for p in points] + ax.plot(xs, ys, color="#999", linewidth=1.5, solid_capstyle="round", solid_joinstyle="round", zorder=1) + + for left_i, right_i in layout["pairs"]: + x1, y1 = layout["positions"][left_i] + x2, y2 = layout["positions"][right_i] + ax.plot([x1 + struct_offset_x, x2 + struct_offset_x], [y1, y2], + color="#b3263a", linewidth=2.4, alpha=0.9, zorder=2) + + for t in layout["terminals"]: + ax.plot([t["x"] + struct_offset_x, t["lx"] + struct_offset_x], [t["y"], t["ly"]], + color="#999", linewidth=0.75, zorder=3) + circ = MplCircle((t["lx"] + struct_offset_x, t["ly"]), t["radius"], + facecolor=TERMINAL_COLORS[t["end_class"]], edgecolor="#666", linewidth=0.5, zorder=4) + ax.add_patch(circ) + ax.text(t["lx"] + struct_offset_x, t["ly"], t["label"], fontsize=4, fontweight="bold", color="#333", + ha="center", va="center", fontfamily="sans-serif", zorder=5) + + for (x, bx_y), base in zip(layout["positions"], layout["labels"]): + circ = MplCircle((x + struct_offset_x, bx_y), 12, facecolor=BASE_COLORS.get(base.upper(), "#ccc"), + edgecolor="#222", linewidth=0.7, zorder=6) + ax.add_patch(circ) + + for (x, bx_y), base in zip(layout["positions"], layout["labels"]): + ax.text(x + struct_offset_x, bx_y, base, fontsize=8, fontweight="bold", color="#111", + ha="center", va="center", fontfamily="sans-serif", zorder=7) + y += row_height + + fig.savefig(str(output_path), format="pdf", bbox_inches="tight", pad_inches=0.1) + plt.close(fig) + + +def _write_svg(svg_text, output_path): + output_path = Path(output_path) + output_path.parent.mkdir(parents=True, exist_ok=True) + output_path.write_text(svg_text) return output_path def plot_dot_bracket_complexes(complexes, output_path, simulation_params=None, file_format="svg"): + if file_format not in ("svg", "pdf"): + raise ValueError("file_format must be 'svg' or 'pdf'") + params = _simulation_params(simulation_params) row_padding = params["row_padding"] - max_radius = max(_complex_radius(complex_data, params) for complex_data in complexes) - row_heights = [2 * _complex_radius(complex_data, params) + row_padding for complex_data in complexes] - text_width = _longest_text_line(complexes, params) * 8 - view_width = int(max(900, 2 * max_radius + params["left_padding"] + params["text_gap"] + text_width + 80)) - width = int(view_width * params["scale"]) - height = int((sum(row_heights) + 60) * params["scale"]) - view_height = int(sum(row_heights) + 60) + max_radius = max(_complex_radius(cd, params) for cd in complexes) + row_heights = [2 * _complex_radius(cd, params) + row_padding for cd in complexes] + text_col_width = _longest_text_line(complexes, params) * 10 + struct_col_width = 2 * max_radius + 80 + view_width = int(params["left_padding"] + text_col_width + params["text_gap"] + struct_col_width + 40) + view_height = int(sum(row_heights) + 40) + + layouts = [] + y = 20 + for complex_data, row_height in zip(complexes, row_heights): + radius = _complex_radius(complex_data, params) + center_x = 0 + center_y = y + radius + row_padding / 2 + layouts.append(_compute_complex_layout(complex_data, center_x, center_y, radius, params)) + y += row_height + + output_path = Path(output_path) + output_path.parent.mkdir(parents=True, exist_ok=True) + if file_format == "pdf": + _render_complexes_to_matplotlib(layouts, text_col_width, view_width, view_height, output_path, params) + return output_path + + text_x = params["left_padding"] + struct_offset_x = text_col_width + params["text_gap"] + width = int(view_width * params["scale"]) + height = int(view_height * params["scale"]) svg = [ - f'', + f'', "", ] - y = 30 - for complex_data, row_height in zip(complexes, row_heights): - radius = _complex_radius(complex_data, params) - center_x = max_radius + params["left_padding"] + y = 20 + for layout, row_height in zip(layouts, row_heights): + radius = layout["radius"] center_y = y + radius + row_padding / 2 - svg.extend(_complex_svg( - complex_data, - center_x, - center_y, - radius, - params, - )) + text_y = center_y - 50 + svg.extend(_complex_svg(layout, text_x, text_y, struct_offset_x)) y += row_height svg.append("") - - return _write_plot_output("\n".join(svg), output_path, file_format) + return _write_svg("\n".join(svg), output_path) def _safe_filename(title): diff --git a/crisscross_kit/sec_struc_pred/nupack_wrapper.py b/crisscross_kit/sec_struc_pred/nupack_wrapper.py index 23fb9fa..050f441 100644 --- a/crisscross_kit/sec_struc_pred/nupack_wrapper.py +++ b/crisscross_kit/sec_struc_pred/nupack_wrapper.py @@ -5,21 +5,13 @@ class NupackTubeConfig: - """Small config object for NUPACK nucleic acid tube analysis.""" - - def __init__( - self, - sequences, - material="dna", - celsius=37, - sodium=0.05, - magnesium=0.015, - max_complex_size=2, - compute=None, - options=None, - tube_name="tube_1", - strand_prefix="strand", - ): + """Small config object for NUPACK nucleic acid tube analysis. + + sequences: list of [sequence_str, concentration_nM] pairs, e.g. [["GCGTATGC", 1000], ["GCATACGC", 500]]. + """ + + def __init__(self, sequences, material="dna", celsius=37, sodium=0.05, magnesium=0.015, + max_complex_size=2, compute=None, options=None, tube_name="tube_1", strand_prefix="strand"): self.sequences = sequences self.material = material self.celsius = celsius @@ -33,16 +25,10 @@ def __init__( def to_dict(self): return { - "sequences": self.sequences, - "material": self.material, - "celsius": self.celsius, - "sodium": self.sodium, - "magnesium": self.magnesium, - "max_complex_size": self.max_complex_size, - "compute": self.compute, - "options": self.options, - "tube_name": self.tube_name, - "strand_prefix": self.strand_prefix, + "sequences": self.sequences, "material": self.material, "celsius": self.celsius, + "sodium": self.sodium, "magnesium": self.magnesium, "max_complex_size": self.max_complex_size, + "compute": self.compute, "options": self.options, + "tube_name": self.tube_name, "strand_prefix": self.strand_prefix, } def write(self, path): @@ -52,30 +38,16 @@ def write(self, path): return path def to_toml(self): - sequences = ",\n ".join( - f'["{sequence}", {concentration_nM}]' - for sequence, concentration_nM in self.sequences - ) + sequences = ",\n ".join(f'["{seq}", {conc}]' for seq, conc in self.sequences) compute = ", ".join(f'"{item}"' for item in self.compute) - options = "\n".join(f"{key} = {value}" for key, value in self.options.items()) - - text = f"""sequences = [ - {sequences} -] - -material = "{self.material}" -celsius = {self.celsius} -sodium = {self.sodium} -magnesium = {self.magnesium} -max_complex_size = {self.max_complex_size} -compute = [{compute}] -tube_name = "{self.tube_name}" -strand_prefix = "{self.strand_prefix}" - -[options] -{options} -""" - return text + options = "\n".join(f"{k} = {v}" for k, v in self.options.items()) + + return (f'sequences = [{sequences}]\n\n' + f'material = "{self.material}"\ncelsius = {self.celsius}\n' + f'sodium = {self.sodium}\nmagnesium = {self.magnesium}\n' + f'max_complex_size = {self.max_complex_size}\ncompute = [{compute}]\n' + f'tube_name = "{self.tube_name}"\nstrand_prefix = "{self.strand_prefix}"\n\n' + f'[options]\n{options}\n') @classmethod def from_file(cls, path): @@ -94,68 +66,50 @@ def _load_config(config): def _complex_to_dict(complex_obj, complex_result, concentration_m, strand_sequences): minimum_free_energy = None if complex_result.mfe: + mfe = complex_result.mfe[0] minimum_free_energy = { - "structure": str(complex_result.mfe[0].structure), - "energy_kcal_per_mol": float(complex_result.mfe[0].energy), - "stack_energy_kcal_per_mol": float(complex_result.mfe[0].stack_energy), + "structure": str(mfe.structure), "energy_kcal_per_mol": float(mfe.energy), + "stack_energy_kcal_per_mol": float(mfe.stack_energy), } return { "name": complex_obj.name, - "strands": [strand.name for strand in complex_obj.strands], - "sequences": [strand_sequences[strand.name] for strand in complex_obj.strands], + "strands": [s.name for s in complex_obj.strands], + "sequences": [strand_sequences[s.name] for s in complex_obj.strands], "concentration_nM": float(concentration_m) * 1e9, - "total_free_energy": complex_result.free_energy, - "minimum_free_energy": minimum_free_energy, + "total_free_energy": complex_result.free_energy, "minimum_free_energy": minimum_free_energy, } def run_nupack_analysis(config): + """Run a NUPACK tube analysis. config.sequences holds [sequence_str, concentration_nM] pairs.""" config = _load_config(config) strands = {} strand_sequences = {} - for i, (sequence, concentration_nM) in enumerate(config["sequences"], start=1): + for i, (sequence, concentration_nM) in enumerate(config["sequences"], start=1): # (seq, conc_nM) pairs strand_name = f"{config.get('strand_prefix', 'strand')}_{i}" strand = Strand(sequence, name=strand_name) strands[strand] = concentration_nM * 1e-9 strand_sequences[strand_name] = sequence - model = Model( - material=config.get("material", "dna"), - celsius=config.get("celsius", 37), - sodium=config.get("sodium", 0.05), - magnesium=config.get("magnesium", 0.015), - ) - tube = Tube( - strands=strands, - complexes=SetSpec(max_size=config.get("max_complex_size", 2)), - name=config.get("tube_name", "tube_1"), - ) - - analysis_result = tube_analysis( - tubes=[tube], - model=model, - compute=config.get("compute", ["pfunc", "mfe"]), - options=config.get("options", {}), - ) + model = Model(material=config.get("material", "dna"), celsius=config.get("celsius", 37), + sodium=config.get("sodium", 0.05), magnesium=config.get("magnesium", 0.015)) + tube = Tube(strands=strands, complexes=SetSpec(max_size=config.get("max_complex_size", 2)), + name=config.get("tube_name", "tube_1")) + + analysis_result = tube_analysis(tubes=[tube], model=model, compute=config.get("compute", ["pfunc", "mfe"]), + options=config.get("options", {})) tube_result = analysis_result[tube] complexes = {} for complex_obj, concentration_m in tube_result.complex_concentrations.items(): complex_result = analysis_result[complex_obj] - complexes[complex_obj.name] = _complex_to_dict( - complex_obj, - complex_result, - concentration_m, - strand_sequences, - ) + complexes[complex_obj.name] = _complex_to_dict(complex_obj, complex_result, concentration_m, strand_sequences) return { "config": config, - "direct_config": config, - "fraction_bases_unpaired": tube_result.fraction_bases_unpaired, - "complexes": complexes, + "fraction_bases_unpaired": tube_result.fraction_bases_unpaired, "complexes": complexes, } @@ -168,63 +122,44 @@ def _format_optional_number(value, digits=3): def _result_rows(result): rows = [] for name, data in result["complexes"].items(): - minimum_free_energy = data["minimum_free_energy"] - minimum_free_energy_value = None - minimum_free_energy_structure = None - if minimum_free_energy is not None: - minimum_free_energy_value = minimum_free_energy["energy_kcal_per_mol"] - minimum_free_energy_structure = minimum_free_energy["structure"] - + mfe = data["minimum_free_energy"] + mfe_value = mfe["energy_kcal_per_mol"] if mfe is not None else None + mfe_structure = mfe["structure"] if mfe is not None else None rows.append({ - "complex": name, - "sequences": " + ".join(data["sequences"]), - "concentration_nM": data["concentration_nM"], - "total_free_energy": data["total_free_energy"], - "minimum_free_energy": minimum_free_energy_value, - "minimum_free_energy_structure": minimum_free_energy_structure, + "complex": name, "sequences": " + ".join(data["sequences"]), + "concentration_nM": data["concentration_nM"], "total_free_energy": data["total_free_energy"], + "minimum_free_energy": mfe_value, "minimum_free_energy_structure": mfe_structure, }) return rows def _summary_text(title, result, rows): - lines = [ - "", - title, - "=" * len(title), - f"Fraction bases unpaired: {_format_optional_number(result['fraction_bases_unpaired'])}", - "", - f"{'Complex':20s} " - f"{'Sequences':30s} " - f"{'Conc nM':>12s} " - f"{'Total free energy':>18s} " - f"{'Minimum free energy':>20s} " - f"{'Dot bracket':>25s}", - "-" * 138, - ] + header = (f"{'Complex':20s} {'Sequences':30s} {'Conc nM':>12s} " + f"{'Total free energy':>18s} {'Minimum free energy':>20s} {'Dot bracket':>25s}") + lines = ["", title, "=" * len(title), + f"Fraction bases unpaired: {_format_optional_number(result['fraction_bases_unpaired'])}", + "", header, "-" * 138] for row in rows: - lines.append( - f"{row['complex']:20s} " - f"{row['sequences']:30s} " - f"{row['concentration_nM']:12.3f} " - f"{_format_optional_number(row['total_free_energy']):>18s} " - f"{_format_optional_number(row['minimum_free_energy']):>20s} " - f"{str(row['minimum_free_energy_structure'] or 'NA'):>25s}" - ) + lines.append(f"{row['complex']:20s} {row['sequences']:30s} {row['concentration_nM']:12.3f} " + f"{_format_optional_number(row['total_free_energy']):>18s} " + f"{_format_optional_number(row['minimum_free_energy']):>20s} " + f"{str(row['minimum_free_energy_structure'] or 'NA'):>25s}") return "\n".join(lines) def _config_rows(config): + """Flatten config dict into parameter/value rows. Sequences are unpacked as (sequence_str, concentration_nM).""" rows = [] for key, value in config.items(): if key == "sequences": - for i, (sequence, concentration_nM) in enumerate(value, start=1): + for i, (sequence, conc) in enumerate(value, start=1): # (seq, conc_nM) pairs rows.append({"parameter": f"sequence_{i}", "value": sequence}) - rows.append({"parameter": f"sequence_{i}_concentration_nM", "value": concentration_nM}) + rows.append({"parameter": f"sequence_{i}_concentration_nM", "value": conc}) elif key == "options": - for option_key, option_value in value.items(): - rows.append({"parameter": f"options.{option_key}", "value": option_value}) + for opt_key, opt_value in value.items(): + rows.append({"parameter": f"options.{opt_key}", "value": opt_value}) else: rows.append({"parameter": key, "value": value}) return rows @@ -232,10 +167,8 @@ def _config_rows(config): def _save_binding_summary(save, result_rows, config_rows): import pandas as pd - excel_path = Path(save) excel_path.parent.mkdir(parents=True, exist_ok=True) - with pd.ExcelWriter(excel_path) as writer: pd.DataFrame(result_rows).to_excel(writer, sheet_name="results", index=False) pd.DataFrame(config_rows).to_excel(writer, sheet_name="config", index=False) @@ -247,17 +180,12 @@ def print_binding_summary(title, result, save=""): print(text) if save: - config = result.get("direct_config", result["config"]) + config = result["config"] output_path = Path(save) / f"{title}.xlsx" _save_binding_summary(output_path, rows, _config_rows(config)) if __name__ == "__main__": - test_config = NupackTubeConfig( - sequences=[ - ["GCGTATGC", 1000], - ["GCATACGC", 1000], - ], - ) + test_config = NupackTubeConfig(sequences=[["GCGTATGC", 1000], ["GCATACGC", 1000]]) result = run_nupack_analysis(test_config) print_binding_summary("NUPACK analysis example", result) From b6b2eb44846cf7dd9b1d5f2e1909d497e53654d2 Mon Sep 17 00:00:00 2001 From: Katzi93_on_mac Date: Thu, 23 Apr 2026 09:22:18 -0400 Subject: [PATCH 6/6] remove spaces from input sequence automatically. script to analize bakteria aptamers --- .../sec_struc_pred/nupack_wrapper.py | 4 + .../bakteria_aptamer/bakteria_aptamer.py | 96 +++++++++++++++++++ 2 files changed, 100 insertions(+) create mode 100644 crisscross_kit/sec_struc_pred/scripts/bakteria_aptamer/bakteria_aptamer.py diff --git a/crisscross_kit/sec_struc_pred/nupack_wrapper.py b/crisscross_kit/sec_struc_pred/nupack_wrapper.py index 050f441..699c8ff 100644 --- a/crisscross_kit/sec_struc_pred/nupack_wrapper.py +++ b/crisscross_kit/sec_struc_pred/nupack_wrapper.py @@ -4,6 +4,10 @@ from nupack import Model, SetSpec, Strand, Tube, tube_analysis +def _normalize_sequence(sequence): + return "".join(str(sequence).split()) + + class NupackTubeConfig: """Small config object for NUPACK nucleic acid tube analysis. diff --git a/crisscross_kit/sec_struc_pred/scripts/bakteria_aptamer/bakteria_aptamer.py b/crisscross_kit/sec_struc_pred/scripts/bakteria_aptamer/bakteria_aptamer.py new file mode 100644 index 0000000..5d6e7be --- /dev/null +++ b/crisscross_kit/sec_struc_pred/scripts/bakteria_aptamer/bakteria_aptamer.py @@ -0,0 +1,96 @@ +from pathlib import Path + +from sec_struc_pred.nupack_wrapper import NupackTubeConfig, print_binding_summary, run_nupack_analysis +from sec_struc_pred.dot_brac_plotter import plot_nupack_result, set_simulation_params + + +def reverse_complement(sequence): + return sequence.translate(str.maketrans("ACGTacgt", "TGCAtgca"))[::-1] + + +NELSON_SOCKET = "AATTACATCTCTCTCCCATCA" +KATZEN_LINKER = "TTA TTA TTA TTA" +NELSON_BINDING_REGION = reverse_complement(NELSON_SOCKET) +#NELSON_SOCKET= "TTTTT"+NELSON_SOCKET +PRIMER_5 = "CGT ACG GAA TTC GCT AGC" +PRIMER_3 = "GGA TCC GAG CTC CAC GTG" + +APTAMER_CORES = { + "STC_03_37": "CATATCCGCGTCGCTGCGCTCAGACCCACCACCACGCACC", + "STC_06_37": "GGGCGGGGGTGCTGGGGGAATGGAGTGCTGCGTGCTGCGG", + "STC_12_37": "GGACCGCAGGTGCACTGGGCGACGTCTCTGGGTGTGGTGT", +} + +OUTPUT_DIR = Path(__file__).resolve().parents[2] / "default_outputs" / "bak_apta" + + +def sequence_variants(core_sequence): + return { + "01_core_only": { + "sequence": core_sequence, + "include_nelson_socket": False, + }, + "02_core_linker_nelson": { + "sequence": f"{core_sequence}{KATZEN_LINKER}{NELSON_BINDING_REGION}", + "include_nelson_socket": True, + }, + "03_primers_core_linker_nelson": { + "sequence": f"{PRIMER_5}{core_sequence}{PRIMER_3}{KATZEN_LINKER}{NELSON_BINDING_REGION}", + "include_nelson_socket": True, + }, + "04_primers_core": { + "sequence": f"{PRIMER_5}{core_sequence}{PRIMER_3}", + "include_nelson_socket": False, + }, + } + + +def run_single_analysis(title, sequence, output_dir, include_nelson_socket=False): + sequences = [[sequence, 1000]] + if include_nelson_socket: + sequences.append([NELSON_SOCKET, 1000]) + + config = NupackTubeConfig( + sequences=sequences, + material="dna", + celsius=25, + sodium=0.05, + magnesium=0.015, + max_complex_size=2 if include_nelson_socket else 1, + compute=["pfunc", "mfe"], + tube_name=title, + strand_prefix="seq", + ) + + result = run_nupack_analysis(config) + print_binding_summary(title, result, save=output_dir) + plot_nupack_result( + title, + result, + save=output_dir, + file_format="pdf", + ) + + +# Parameters for the plotter. We use a relaxation algorithm. Its a toy version of a force field simulation. +set_simulation_params(brownian_jitter=0.0) +set_simulation_params(relaxation_steps=15000) +set_simulation_params(backbone_straightness=0.01) +set_simulation_params(centering=0.0002) +# set_simulation_params(reference_sequence_length=300) # "start_radius": 95, "reference_sequence_length": 16 + + +if __name__ == "__main__": + OUTPUT_DIR.mkdir(parents=True, exist_ok=True) + + for aptamer_name, core_sequence in APTAMER_CORES.items(): + aptamer_output_dir = OUTPUT_DIR / aptamer_name + aptamer_output_dir.mkdir(parents=True, exist_ok=True) + + for variant_name, variant in sequence_variants(core_sequence).items(): + run_single_analysis( + f"{aptamer_name}_{variant_name}", + variant["sequence"], + aptamer_output_dir, + include_nelson_socket=variant["include_nelson_socket"], + )