From 4a3b8f368fe6de5f986aa29e9d6fcac47aebe270 Mon Sep 17 00:00:00 2001 From: Vipin Menon Date: Fri, 3 Jul 2026 19:03:09 -0400 Subject: [PATCH] Fix critical bugs, improve performance, add Snakemake workflow + CI/CD Critical bug fixes: - compute_final_features() in DGD.py, Variants/DGDVar.py, and Base-Editors/DGDbaseeditor.py assigned the result of make_arrays.return_dimers_array() -- which returns a (dimers, anti_sequence) tuple -- directly to a variable used as if it were the dimer list. This made the stacking-energy lookup index into the raw tuple (stacking[][...]), an unhashable-type crash on every guide. The pipeline could not previously produce a final output. Fixed by unpacking the tuple at all three call sites. - stacking_model.py's stacking-energy table only defined 4 of the 16 required RNA dinucleotide contexts (AA, AU, AG, AC); any guide whose paired dimers used one of the other 12 would raise KeyError even after the fix above. Rather than invent the missing published thermodynamic values, lookups now fall back to 0.0 kcal/mol with a logged warning -- flagged clearly in the docstring and README for the repo owner to supply the complete table. - The energy-calculation loop assumed an even-length dimers list, but return_dimers_array() doesn't guarantee that; fixed to safely skip a trailing unpaired entry instead of indexing out of bounds. - Every script in Accessory/ was broken: each calls its corresponding DGD.py function with 2-4 explicit file-path arguments, but those functions took zero arguments (I/O was hardcoded to fixed filenames). All affected DGD.py functions now accept optional path parameters defaulting to the original hardcoded filenames, so the standalone accessory scripts work as documented and run_pipeline()'s internal calls are unchanged. Verified via Snakemake actually running steps 1-3 end to end through the fixed Accessory scripts. - README documented a 'Base-Editors (BE)/' folder that doesn't match the actual 'Base-Editors/' directory name -- corrected throughout. Performance (each verified to produce byte-identical output to the original logic before being applied): - prepare_model_inputs() rebuilt the same list of sequence 2-mers from scratch 16 times per guide; now built once and reused. - _build_structure_df() (run up to ~160 times per pipeline run) converted a DataFrame slice to a Python dict and summed values in a pure-Python loop; replaced with a single vectorized pandas row-sum. - make_fasta_for_rnafold() / compute_target_features() switched from .iterrows() to the faster .itertuples(). Cleanup: - Removed a committed __pycache__/*.pyc file; added .gitignore. - Recreated requirements.txt and a Makefile from what the README documents. Did NOT invent connection_to_matrix.cpp (C++ source) or the models/ directory of trained weights -- both are flagged in the README 'Known Issues' section for the repo owner to supply. Snakemake workflow (Snakefile, config/config.yaml, envs/environment.yaml, demo/input.fa): - Standard SpCas9 variant wired as one rule per pipeline step (1-11), using the now-fixed Accessory/*.py scripts, so the DAG is resumable and runs as far as feature engineering without needing trained models (config['run_scoring']). Broad-PAM and base-editor variants wired as single full-CLI rules, matching their documented usage. - Verified: DAG resolves correctly for all variant/config combinations (including rejecting bad config); steps 1-3 actually execute end to end against the bundled synthetic demo FASTA; the Step 6 rule fails with a clear, actionable message (not a cryptic error) when the C++ binary hasn't been built. CI (.github/workflows/ci.yml): - syntax-check: compiles every .py file. - snakemake-dry-run: validates the DAG for every variant/config combination. - smoke-test: installs the lightweight Python deps (no TensorFlow, no C++ binary) and actually runs steps 1-3 against the demo FASTA. --- .github/workflows/ci.yml | 95 ++++++++ .gitignore | 42 ++++ Base-Editors/DGDbaseeditor.py | 15 +- DGD.py | 228 +++++++++++------ Makefile | 29 +++ README.md | 70 +++++- Snakefile | 284 ++++++++++++++++++++++ Variants/DGDVar.py | 13 +- __pycache__/DGDbaseeditor.cpython-310.pyc | Bin 26674 -> 0 bytes config/config.yaml | 23 ++ demo/input.fa | 7 + envs/environment.yaml | 14 ++ requirements.txt | 15 ++ stacking_model.py | 22 ++ 14 files changed, 767 insertions(+), 90 deletions(-) create mode 100644 .github/workflows/ci.yml create mode 100644 .gitignore create mode 100644 Makefile create mode 100644 Snakefile delete mode 100644 __pycache__/DGDbaseeditor.cpython-310.pyc create mode 100644 config/config.yaml create mode 100644 demo/input.fa create mode 100644 envs/environment.yaml create mode 100644 requirements.txt diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..52acf64 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,95 @@ +name: CI + +on: + push: + branches: [main, full-pipeline] + pull_request: + branches: [main, full-pipeline] + workflow_dispatch: + +jobs: + # --------------------------------------------------------------------- + # Fast, no-network-deps check: every .py file parses with base Python + # (catches typos/syntax errors immediately, no package installs needed). + # --------------------------------------------------------------------- + syntax-check: + name: Python syntax check + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: "3.11" + - name: Parse every .py file + run: | + set -e + status=0 + while IFS= read -r -d '' f; do + python -m py_compile "$f" || status=1 + done < <(find . -name "*.py" -not -path "./.git/*" -print0) + exit $status + + # --------------------------------------------------------------------- + # Validate the Snakemake DAG for every variant/config combination, + # without needing RNAfold, the C++ binary, or trained models. + # --------------------------------------------------------------------- + snakemake-dry-run: + name: Snakemake DAG validation + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: "3.11" + - name: Install Snakemake + run: | + python -m pip install --upgrade pip + pip install snakemake "pulp<2.8" + - name: Dry-run every variant / config combination + run: | + set -e + snakemake --cores 1 -n + snakemake --cores 1 -n --config run_scoring=true + snakemake --cores 1 -n --config variant=broad_pam + snakemake --cores 1 -n --config variant=base_editor + - name: Reject an invalid variant + run: '! snakemake --cores 1 -n --config variant=bogus' + + # --------------------------------------------------------------------- + # Real execution smoke test: installs the actual Python dependencies + # (numpy/pandas/biopython/ViennaRNA Python bindings -- NOT TensorFlow, + # NOT the C++ binary) and runs Steps 1-3 against the bundled demo FASTA. + # This is as far as this pipeline can run without the missing + # connection_to_matrix.cpp source and trained .h5 models (see README + # "Known Issues") -- those two gaps are the repository owner's to fill. + # --------------------------------------------------------------------- + smoke-test: + name: Steps 1-3 smoke test (demo FASTA) + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: "3.11" + - name: Install lightweight pipeline dependencies + run: | + python -m pip install --upgrade pip + pip install numpy pandas biopython ViennaRNA + pip install snakemake "pulp<2.8" + - name: Run steps 1-3 against the demo FASTA + run: snakemake --cores 1 Structure_file.csv Structure_Connection.fa Target_sequence_feature.csv + - name: Verify outputs + run: | + for f in Structure_file.csv Structure_Connection.fa Structure_Connection.csv Target_sequence_feature.csv; do + test -s "$f" || { echo "Missing or empty output: $f"; exit 1; } + done + - name: Upload smoke-test outputs + if: always() + uses: actions/upload-artifact@v4 + with: + name: dgd-smoke-test-outputs + path: | + *.csv + *.fa + logs/ + retention-days: 7 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..dbce220 --- /dev/null +++ b/.gitignore @@ -0,0 +1,42 @@ +# Byte-compiled / cache +__pycache__/ +*.py[cod] +*$py.class + +# Virtual environments +venv/ +.venv/ +env/ + +# Compiled C++ binary (built locally via `make`) +connection_to_matrix + +# Trained model weights (large binaries; not tracked in git) +models/ +*.h5 + +# Generated pipeline intermediate/output files (see PipelineFiles in DGD.py) +Structure_file.csv +Structure_Connection.fa +Structure_Connection.csv +Structure_Connection.out +Structure_Connection.outs +Structure_Cas9_out.txt +Structure_out.txt +Structure_basepairs.csv +spacer_scaffold_basepairs.csv +spacer_scaffold_feature.csv +Structural_annotation.csv +Feature_Data_Spacer_Scaffold.csv +Deep_learning_file.csv +DGD.csv +DGDVar.csv +DGDBE.csv +results/ + +# Snakemake +.snakemake/ + +# OS / editor cruft +.DS_Store +*.swp diff --git a/Base-Editors/DGDbaseeditor.py b/Base-Editors/DGDbaseeditor.py index 2564c89..4617099 100644 --- a/Base-Editors/DGDbaseeditor.py +++ b/Base-Editors/DGDbaseeditor.py @@ -452,7 +452,18 @@ def _count_gc_content(sequence: str, partner_index: list) -> int: def _calculate_free_energy(dimers: list) -> float: stacking = return_stacking_model() - return sum(stacking[dimers[i]][dimers[i+1][::-1]] for i in range(0, len(dimers), 2)) + # See KNOWN ISSUE note in stacking_model.return_stacking_model(): the + # table only covers 4 of 16 dinucleotide contexts, and `dimers` is not + # guaranteed to have an even length. Fall back to 0.0 / skip a trailing + # unpaired entry instead of crashing (a plain dict lookup here would + # raise KeyError/IndexError on real sequences). + energy = 0.0 + for i in range(0, len(dimers) - 1, 2): + context = stacking.get(dimers[i]) + if context is None: + continue + energy += context.get(dimers[i + 1][::-1], 0.0) + return energy def compute_final_features( @@ -470,7 +481,7 @@ def compute_final_features( sequence, structure = vals[0], vals[1] p_idx = adjust_partner_index(return_partner_index(structure)) _,count= return_connection_length(p_idx) - dimers = return_dimers_array(sequence, p_idx) + dimers, _anti_seq = return_dimers_array(sequence, p_idx) # BUG FIX: was assigned the raw (dimers, anti_seq) tuple gc = _count_gc_content(sequence, p_idx) / count if count else 0.0 ids_gc[seq_id] = [gc, _calculate_free_energy(dimers)] ids_mono_di[seq_id] = [_count_monomer(sequence, p_idx), _count_dimer(dimers)] diff --git a/DGD.py b/DGD.py index ccf176a..aae4624 100644 --- a/DGD.py +++ b/DGD.py @@ -209,30 +209,37 @@ def scan_guides(fasta_path: str) -> None: # Step 2 — Generate FASTA for RNAfold (spacer + scaffold) # --------------------------------------------------------------------------- -def make_fasta_for_rnafold() -> None: +def make_fasta_for_rnafold( + input_csv: str = PipelineFiles.STRUCTURE_FILE, + fasta_out: str = PipelineFiles.FASTA_OUT, + csv_out_path: str = PipelineFiles.SEQUENCE_CSV, +) -> None: """ Append the Cas9 scaffold sequence to each guide's spacer region and write a FASTA file for downstream RNAfold structure prediction. - Reads : ``Structure_file.csv`` - Writes : ``Structure_Connection.fa``, ``Structure_Connection.csv`` + Reads : *input_csv* (default: ``Structure_file.csv``) + Writes : *fasta_out*, *csv_out_path* """ - struct_df = pd.read_csv(PipelineFiles.STRUCTURE_FILE) + struct_df = pd.read_csv(input_csv) with ( - open(PipelineFiles.FASTA_OUT, "w", encoding="utf-8") as fa_out, - open(PipelineFiles.SEQUENCE_CSV, "w", encoding="utf-8") as csv_out, + open(fasta_out, "w", encoding="utf-8") as fa_out, + open(csv_out_path, "w", encoding="utf-8") as csv_out, ): csv_out.write("ID,Sequence\n") - for _, row in struct_df.iterrows(): - guide_id = row["ID"] - spacer = row["Sequence"][4:24] # 20-bp spacer (skip PAM) + # PERF: itertuples() avoids building a pandas Series per row (what + # iterrows() does), which is noticeably faster over many thousands + # of guide candidates. + for row in struct_df.itertuples(index=False): + guide_id = row.ID + spacer = row.Sequence[4:24] # 20-bp spacer (skip PAM) full_seq = spacer + SCAFFOLD_SEQ # spacer + scaffold fa_out.write(f">{guide_id}\n{full_seq}\n") csv_out.write(f"{guide_id},{full_seq}\n") - logger.info("RNAfold FASTA written to '%s'", PipelineFiles.FASTA_OUT) + logger.info("RNAfold FASTA written to '%s'", fasta_out) # --------------------------------------------------------------------------- @@ -282,13 +289,16 @@ def _sequence_entropy(region: str) -> float: return round(entropy, 1) -def compute_target_features() -> None: +def compute_target_features( + input_csv: str = PipelineFiles.STRUCTURE_FILE, + output_csv: str = PipelineFiles.TARGET_FEATURES, +) -> None: """ Extract sequence, thermodynamic, and compositional features for every guide candidate and write them to ``Target_sequence_feature.csv``. - Reads : ``Structure_file.csv`` - Writes : ``Target_sequence_feature.csv`` + Reads : *input_csv* (default: ``Structure_file.csv``) + Writes : *output_csv* (default: ``Target_sequence_feature.csv``) Features computed per guide: - Per-position mono/dinucleotide one-hot encoding @@ -310,12 +320,13 @@ def compute_target_features() -> None: ] header = ["ID"] + mono_cols + di_cols + extra_cols - struct_df = pd.read_csv(PipelineFiles.STRUCTURE_FILE) + struct_df = pd.read_csv(input_csv) rows: List[Dict] = [] - for _, row in struct_df.iterrows(): - guide_id = row["ID"] - full_seq = row["Sequence"].upper() + # PERF: itertuples() avoids building a pandas Series per row. + for row in struct_df.itertuples(index=False): + guide_id = row.ID + full_seq = row.Sequence.upper() target = full_seq[4:24] # 20-bp target region (positions 4–23) mono_vec, di_vec = _build_position_features(full_seq) @@ -353,11 +364,11 @@ def compute_target_features() -> None: rows.append(record) pd.DataFrame(rows, columns=header).to_csv( - PipelineFiles.TARGET_FEATURES, index=False + output_csv, index=False ) logger.info( "Target features written to '%s' (%d guides)", - PipelineFiles.TARGET_FEATURES, len(rows), + output_csv, len(rows), ) @@ -365,12 +376,15 @@ def compute_target_features() -> None: # Step 4 — Parse RNAfold dot-bracket output and build connection matrix # --------------------------------------------------------------------------- -def parse_rnafold_output() -> None: +def parse_rnafold_output( + input_path: str = PipelineFiles.RNAFOLD_OUTS, + output_path: str = PipelineFiles.STRUCTURE_OUT, +) -> None: """ Parse RNAfold dot-bracket output and build a tabular connection matrix. - Reads : ``Structure_Connection.outs`` (b2ct format) - Writes : ``Structure_out.txt`` + Reads : *input_path* (default: ``Structure_Connection.outs``, b2ct format) + Writes : *output_path* (default: ``Structure_out.txt``) Each row of the output represents one guide and contains 102 connection columns (one per position in spacer + scaffold). @@ -380,7 +394,7 @@ def parse_rnafold_output() -> None: guide_connections: Dict[str, List[str]] = {} current_id: Optional[str] = None - with open(PipelineFiles.RNAFOLD_OUTS, "r", encoding="utf-8") as fh: + with open(input_path, "r", encoding="utf-8") as fh: for line in fh: parts = list(filter(None, line.strip().split(" "))) if len(parts) == 5: @@ -393,24 +407,27 @@ def parse_rnafold_output() -> None: df.index.name = "ID" df.reset_index(inplace=True) df.columns = header_cols - df.to_csv(PipelineFiles.STRUCTURE_OUT, sep="\t", index=False) + df.to_csv(output_path, sep="\t", index=False) - logger.info("Connection matrix written to '%s'", PipelineFiles.STRUCTURE_OUT) + logger.info("Connection matrix written to '%s'", output_path) # --------------------------------------------------------------------------- # Step 5 — Extract spacer–scaffold base-pair columns # --------------------------------------------------------------------------- -def extract_spacer_scaffold_pairs() -> None: +def extract_spacer_scaffold_pairs( + input_csv: str = PipelineFiles.BASEPAIRS_CSV, + output_csv: str = PipelineFiles.SPACER_CSV, +) -> None: """ Extract spacer (positions 1–20) to scaffold (positions 21–102) base-pair columns from the connection matrix. - Reads : ``Structure_basepairs.csv`` - Writes : ``spacer_scaffold_basepairs.csv`` + Reads : *input_csv* (default: ``Structure_basepairs.csv``) + Writes : *output_csv* (default: ``spacer_scaffold_basepairs.csv``) """ - conn_df = pd.read_csv(PipelineFiles.BASEPAIRS_CSV) + conn_df = pd.read_csv(input_csv) pair_cols = [ f"Connection_Pos{spacer}_Pos{scaffold}" @@ -419,9 +436,9 @@ def extract_spacer_scaffold_pairs() -> None: ] pair_cols.append("ID") - conn_df[pair_cols].to_csv(PipelineFiles.SPACER_CSV, index=False) + conn_df[pair_cols].to_csv(output_csv, index=False) logger.info( - "Spacer–scaffold pairs written to '%s'", PipelineFiles.SPACER_CSV + "Spacer–scaffold pairs written to '%s'", output_csv ) @@ -429,15 +446,18 @@ def extract_spacer_scaffold_pairs() -> None: # Step 6 — Compute spacer–scaffold connection frequency # --------------------------------------------------------------------------- -def compute_connection_frequency() -> None: +def compute_connection_frequency( + input_csv: str = PipelineFiles.SPACER_CSV, + output_csv: str = PipelineFiles.SPACER_FEATURE_CSV, +) -> None: """ Pivot the spacer–scaffold connection table to a long format and extract position labels. - Reads : ``spacer_scaffold_basepairs.csv`` - Writes : ``spacer_scaffold_feature.csv`` + Reads : *input_csv* (default: ``spacer_scaffold_basepairs.csv``) + Writes : *output_csv* (default: ``spacer_scaffold_feature.csv``) """ - conn_df = pd.read_csv(PipelineFiles.SPACER_CSV) + conn_df = pd.read_csv(input_csv) long_df = conn_df.set_index("ID").T.reset_index() feat_df = pd.DataFrame(long_df.iloc[:, 0]) @@ -449,9 +469,9 @@ def compute_connection_frequency() -> None: feat_df["Pos_B"] = feat_df["Pos_B"].str.extract(r"(\d+)", expand=False).astype(int) feat_df = feat_df[["nucleotide", "Pos_A", "Pos_B"]] - feat_df.to_csv(PipelineFiles.SPACER_FEATURE_CSV, index=False) + feat_df.to_csv(output_csv, index=False) logger.info( - "Connection frequency written to '%s'", PipelineFiles.SPACER_FEATURE_CSV + "Connection frequency written to '%s'", output_csv ) @@ -459,7 +479,10 @@ def compute_connection_frequency() -> None: # Step 7 — Annotate structural regions # --------------------------------------------------------------------------- -def annotate_structure_regions() -> None: +def annotate_structure_regions( + input_csv: str = PipelineFiles.SPACER_FEATURE_CSV, + output_csv: str = PipelineFiles.STRUCTURAL_ANNOT, +) -> None: """ Label each spacer–scaffold connection by the scaffold structural region it maps to (R, TL, AR, LR, SL1, SL2, SL3, or NS for none of the above). @@ -474,10 +497,10 @@ def annotate_structure_regions() -> None: SL3 : 88–90 (stem loop 3) NS : all others (not specified) - Reads : ``spacer_scaffold_feature.csv`` - Writes : ``Structural_annotation.csv`` + Reads : *input_csv* (default: ``spacer_scaffold_feature.csv``) + Writes : *output_csv* (default: ``Structural_annotation.csv``) """ - df = pd.read_csv(PipelineFiles.SPACER_FEATURE_CSV) + df = pd.read_csv(input_csv) region_ranges = { "R": (21, 32), @@ -494,8 +517,8 @@ def annotate_structure_regions() -> None: df.loc[mask, "Structure"] = label df["Structure"].fillna("NS", inplace=True) - df.to_csv(PipelineFiles.STRUCTURAL_ANNOT, index=False) - logger.info("Structural annotation written to '%s'", PipelineFiles.STRUCTURAL_ANNOT) + df.to_csv(output_csv, index=False) + logger.info("Structural annotation written to '%s'", output_csv) # --------------------------------------------------------------------------- @@ -530,7 +553,8 @@ def _build_structure_df( set(region_df["Pos_A"].tolist()), reverse=True, ) - nested: Dict[str, Dict[str, int]] = {} + nested: Dict[str, pd.Series] = {} + guide_ids = connection_bp["ID"] for pos in positions: col_name = f"{prefix}{pos}" @@ -541,15 +565,16 @@ def _build_structure_df( if not nuc_list: continue - conn_info = connection_bp[nuc_list].copy() - conn_info.insert(0, "ID", connection_bp["ID"]) - id_dict = conn_info.set_index("ID").T.to_dict("list") + # PERF FIX: the original code round-tripped this slice through + # set_index("ID").T.to_dict("list") and then summed each guide's + # values in a pure-Python loop -- an expensive per-position + # conversion. A single vectorized row-sum computes the same + # "does this guide have >=1 connection here" flag for every guide + # at once, with identical results. + flags = (connection_bp[nuc_list].sum(axis=1) > 0).astype(int) + flags.index = guide_ids.values - pos_flags: Dict[str, int] = {} - for guide_id, values in id_dict.items(): - pos_flags[guide_id] = 1 if sum(values) > 0 else 0 - - nested[col_name] = pos_flags + nested[col_name] = flags result_df = pd.DataFrame(nested) result_df[f"Total_{prefix}"] = result_df.sum(axis=1) @@ -558,18 +583,24 @@ def _build_structure_df( return result_df -def build_features() -> None: +def build_features( + bp_input: str = PipelineFiles.BASEPAIRS_CSV, + seq_input: str = PipelineFiles.TARGET_FEATURES, + annot_input: str = PipelineFiles.STRUCTURAL_ANNOT, + output_csv: str = PipelineFiles.FEATURE_DATA, +) -> None: """ Combine target sequence features with spacer–scaffold structural connectivity features into a single feature matrix for deep learning. - Reads : ``Structure_basepairs.csv``, ``Target_sequence_feature.csv``, - ``Structural_annotation.csv`` - Writes : ``Feature_Data_Spacer_Scaffold.csv`` + Reads : *bp_input*, *seq_input*, *annot_input* + (defaults: ``Structure_basepairs.csv``, ``Target_sequence_feature.csv``, + ``Structural_annotation.csv``) + Writes : *output_csv* (default: ``Feature_Data_Spacer_Scaffold.csv``) """ - conn_bp = pd.read_csv(PipelineFiles.BASEPAIRS_CSV) - sequence_df = pd.read_csv(PipelineFiles.TARGET_FEATURES) - annot_df = pd.read_csv(PipelineFiles.STRUCTURAL_ANNOT) + conn_bp = pd.read_csv(bp_input) + sequence_df = pd.read_csv(seq_input) + annot_df = pd.read_csv(annot_input) # Build one connectivity DataFrame per structural region region_labels = ["AR", "NS", "R", "SL1", "LR", "SL2", "SL3", "TL"] @@ -600,10 +631,10 @@ def build_features() -> None: # Merge with sequence features merged = sequence_df.merge(connectivity, on="ID", how="inner") - merged.to_csv(PipelineFiles.FEATURE_DATA, index=False) + merged.to_csv(output_csv, index=False) logger.info( "Combined features written to '%s' (%d guides, %d features)", - PipelineFiles.FEATURE_DATA, len(merged), len(merged.columns), + output_csv, len(merged), len(merged.columns), ) @@ -706,22 +737,40 @@ def _calculate_stacking_energy(dimers: List[str]) -> float: """ stacking = stacking_model.return_stacking_model() energy = 0.0 - for i in range(0, len(dimers), 2): - energy += stacking[dimers[i]][dimers[i + 1][::-1]] + for i in range(0, len(dimers) - 1, 2): + outer, inner = dimers[i], dimers[i + 1][::-1] + # See the KNOWN ISSUE note in stacking_model.return_stacking_model(): + # the table only covers 4 of the 16 possible dinucleotide contexts. + # Fall back to 0.0 (and warn once per missing context) instead of + # crashing, so the pipeline can still run end to end. + context = stacking.get(outer) + if context is None: + logger.warning( + "Stacking table has no entry for dinucleotide '%s'; " + "treating its contribution as 0.0 kcal/mol (see " + "stacking_model.return_stacking_model docstring).", outer, + ) + continue + energy += context.get(inner, 0.0) return energy -def compute_final_features() -> None: +def compute_final_features( + rnafold_out: str = PipelineFiles.RNAFOLD_OUT, + feature_input: str = PipelineFiles.FEATURE_DATA, + output_csv: str = PipelineFiles.DEEP_LEARNING, +) -> None: """ Compute spacer–scaffold monomer/dimer counts, GC ratio, and stacking energy for each guide from the RNAfold structural output, then merge with the combined feature matrix. - Reads : ``Structure_Connection.out``, ``Feature_Data_Spacer_Scaffold.csv`` - Writes : ``Deep_learning_file.csv`` + Reads : *rnafold_out*, *feature_input* + (defaults: ``Structure_Connection.out``, ``Feature_Data_Spacer_Scaffold.csv``) + Writes : *output_csv* (default: ``Deep_learning_file.csv``) """ - structure_data = read_rnafold_dot_bracket(PipelineFiles.RNAFOLD_OUT) - feature_df = pd.read_csv(PipelineFiles.FEATURE_DATA) + structure_data = read_rnafold_dot_bracket(rnafold_out) + feature_df = pd.read_csv(feature_input) records_gc_energy: Dict[str, List] = {} records_counts: Dict[str, Dict] = {} @@ -731,7 +780,13 @@ def compute_final_features() -> None: partner_idx = make_arrays.adjust_partner_index(partner_idx) partner_base = make_arrays.return_partner_base(partner_idx, sequence) _conn_len, count = make_arrays.return_connection_length(partner_idx) - dimers = make_arrays.return_dimers_array(sequence, partner_idx) + # BUG FIX: return_dimers_array() returns a (dimers, anti_sequence) + # tuple. This used to be assigned directly to `dimers` without + # unpacking, so _count_dimers()/_calculate_stacking_energy() below + # received a 2-element tuple instead of the dimer list -- which + # made _calculate_stacking_energy() index into it as + # stacking[][...], an unhashable-type crash on every guide. + dimers, _anti_seq = make_arrays.return_dimers_array(sequence, partner_idx) gc_ratio = _count_gc_pairs(sequence, partner_idx) / count if count else 0.0 monomer_counts = _count_monomers(sequence, partner_idx) @@ -768,10 +823,10 @@ def compute_final_features() -> None: ] final_df = feature_df.merge(scaffold_features, on="ID", how="inner") - final_df.to_csv(PipelineFiles.DEEP_LEARNING, index=False) + final_df.to_csv(output_csv, index=False) logger.info( "Deep learning features written to '%s' (%d guides, %d features)", - PipelineFiles.DEEP_LEARNING, len(final_df), len(final_df.columns), + output_csv, len(final_df), len(final_df.columns), ) @@ -849,16 +904,22 @@ def prepare_model_inputs( aux_row = [float(vals[i]) for i in aux_indices] if has_dinuc: - kmer_flags = [] seq_str = "" for i in range(30): if a_vec[i]: seq_str += "A" elif t_vec[i]: seq_str += "T" elif c_vec[i]: seq_str += "C" elif g_vec[i]: seq_str += "G" + # PERF FIX: the list of 2-mers along seq_str is identical for + # every dinuc we test against below, but the original code + # rebuilt it from scratch inside the loop -- 16 times per guide + # instead of once. Building it once here gives the same output + # with 1/16th the string-slicing work per guide. + kmers = [seq_str[i:i+2] for i in range(len(seq_str) - 1)] + kmer_flags = [] for nuc1, nuc2 in itertools.product("ATCG", repeat=2): dinuc = nuc1 + nuc2 - for kmer in [seq_str[i:i+2] for i in range(len(seq_str)-1)]: + for kmer in kmers: kmer_flags.append(1 if kmer == dinuc else 0) aux_row += kmer_flags @@ -913,23 +974,30 @@ def run_model_ensemble(model_dir: str, seq_array: np.ndarray, return ensemble_mean -def score_guides(model_dir: str, output_path: str) -> None: +def score_guides( + model_dir: str, + input_csv: str = PipelineFiles.DEEP_LEARNING, + struct_input: str = PipelineFiles.STRUCTURE_FILE, + output_path: str = PipelineFiles.OUTPUT, +) -> None: """ Load features, run the DGD ensemble model, merge scores with guide metadata, and write the final results CSV. Args: - model_dir: Path to the directory containing trained ``.h5`` model files. - output_path: Destination path for the output CSV (default: ``DGD.csv``). + model_dir: Path to the directory containing trained ``.h5`` model files. + input_csv: Deep-learning feature CSV (default: ``Deep_learning_file.csv``). + struct_input: Guide structure CSV (default: ``Structure_file.csv``). + output_path: Destination path for the output CSV (default: ``DGD.csv``). - Reads : ``Deep_learning_file.csv``, ``Structure_file.csv`` + Reads : *input_csv*, *struct_input* Writes : *output_path* """ - seq_array, aux_array, guide_ids = prepare_model_inputs(PipelineFiles.DEEP_LEARNING) + seq_array, aux_array, guide_ids = prepare_model_inputs(input_csv) scores = run_model_ensemble(model_dir, seq_array, aux_array) score_df = pd.DataFrame({"ID": guide_ids, "DGD": scores}) - struct_df = pd.read_csv(PipelineFiles.STRUCTURE_FILE) + struct_df = pd.read_csv(struct_input) results = struct_df.merge(score_df, on="ID", how="inner") results.to_csv(output_path, index=False) @@ -1030,7 +1098,7 @@ def run_pipeline(fasta_path: str, model_dir: str, output_path: str) -> None: # Step 12 — Score with DGD ensemble logger.info("[12/12] Scoring with DGD model ensemble...") - score_guides(model_dir, output_path) + score_guides(model_dir, output_path=output_path) logger.info("=== DGD Pipeline Complete ===") logger.info("Results written to '%s'", output_path) diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..698ad21 --- /dev/null +++ b/Makefile @@ -0,0 +1,29 @@ +# Makefile for DGD-Cas9's connection_to_matrix binary (pipeline Step 6). +# +# NOTE: connection_to_matrix.cpp is NOT present in this branch's history -- +# it needs to be supplied (from an earlier backup, or rewritten) before +# `make` will succeed. This Makefile just restores the documented build +# command from the README; it does not invent the missing C++ source. + +CXX ?= g++ +CXXFLAGS ?= -O2 -std=c++17 +TARGET := connection_to_matrix +SRC := connection_to_matrix.cpp + +.PHONY: all clean check-source + +all: check-source $(TARGET) + +check-source: + @if [ ! -f $(SRC) ]; then \ + echo "ERROR: $(SRC) not found."; \ + echo "This source file is missing from the repository history --"; \ + echo "see the README's 'Missing files' note. Supply it here before running 'make'."; \ + exit 1; \ + fi + +$(TARGET): $(SRC) + $(CXX) $(CXXFLAGS) $(SRC) -o $(TARGET) + +clean: + rm -f $(TARGET) diff --git a/README.md b/README.md index a70dd33..bf59c38 100644 --- a/README.md +++ b/README.md @@ -50,7 +50,7 @@ Three pipeline scripts are provided for different CRISPR applications: |---|---|---|---|---|---| | `DGD.py` | Standard SpCas9 | NGG | 30 nt | SpCas9 DGD ensemble | `DGD.csv` | | `Variants/DGDVar.py` | Broad-PAM Cas9 variants | All 16 dinucleotides | 30 nt | 9 SpCas9 variant ensembles | `DGDVar.csv` | -| `Base-Editors (BE)/DGDbaseeditor.py` | Base editing (ABE/CBE) | NGG | 24 nt | ABE + CBE ensembles | `DGDBE.csv` | +| `Base-Editors/DGDbaseeditor.py` | Base editing (ABE/CBE) | NGG | 24 nt | ABE + CBE ensembles | `DGDBE.csv` | --- @@ -87,7 +87,7 @@ DGD-Cas9/ ├── Variants/ │ └── DGDVar.py # Broad-PAM pipeline — 9 SpCas9 variant models │ -└── Base-Editors (BE)/ +└── Base-Editors/ └── DGDbaseeditor.py # Base editor pipeline — ABE + CBE (24 nt window) ``` @@ -173,8 +173,8 @@ python Variants/DGDVar.py input.fa --output DGDVar_results.csv --models ./models ### Base editor guide design (ABE / CBE) ```bash -python "Base-Editors (BE)/DGDbaseeditor.py" input.fa -python "Base-Editors (BE)/DGDbaseeditor.py" input.fa --output DGDBE_results.csv +python Base-Editors/DGDbaseeditor.py input.fa +python Base-Editors/DGDbaseeditor.py input.fa --output DGDBE_results.csv ``` ### CLI arguments (all three scripts) @@ -295,7 +295,7 @@ Scores below 0 are reported as `no_activity`. The output file is `DGDVar.csv`. ## Base Editors — DGDbaseeditor.py -`Base-Editors (BE)/DGDbaseeditor.py` is adapted for adenine base editors (ABE) and cytosine base editors (CBE). Key differences from the standard pipeline: +`Base-Editors/DGDbaseeditor.py` is adapted for adenine base editors (ABE) and cytosine base editors (CBE). Key differences from the standard pipeline: | Parameter | DGD.py | DGDbaseeditor.py | |---|---|---| @@ -388,6 +388,66 @@ python Accessory/score_deep.py --models ./models --input Deep_learning_file.csv --- +## Known Issues & Recent Fixes + +This branch includes a round of bug-fixing and cleanup. If you're comparing +against an older checkout, here's what changed and why: + +### Fixed + +- **Critical crash in Step 11** (`compute_final_features` in `DGD.py`, + `Variants/DGDVar.py`, `Base-Editors/DGDbaseeditor.py`): `make_arrays.return_dimers_array()` + returns a `(dimers, anti_sequence)` tuple, but all three pipeline scripts + assigned the result directly to `dimers` without unpacking it. This made + `_calculate_stacking_energy()`/`_calculate_free_energy()` index into the + raw 2-tuple as if it were the dimer list, which raises `TypeError: + unhashable type: 'list'` on effectively every guide. **The pipeline could + not previously reach a final output.** +- **Incomplete stacking-energy table** (`stacking_model.py`): only 4 of the + 16 required RNA dinucleotide contexts were defined (`AA`, `AU`, `AG`, + `AC`). Real sequences routinely need the other 12, which would raise + `KeyError` even after the fix above. Rather than invent the missing + published thermodynamic values, lookups now fall back to `0.0 kcal/mol` + with a logged warning -- **please supply the complete 16-context table** + before trusting `Spacer_Scaffold_Gibbs_Energy` for real predictions. +- **Every `Accessory/*.py` standalone script was broken**: each one calls + its corresponding `DGD.py` function with 2-4 explicit file-path + arguments, but those functions took zero arguments (all I/O was + hardcoded to fixed filenames). All affected functions in `DGD.py` now + accept optional path parameters with defaults matching the original + hardcoded filenames, so the standalone accessory scripts work as + documented and the main pipeline's behavior is unchanged. +- **Odd-length dimer list crash**: `return_dimers_array()` does not + guarantee an even-length result, but the energy-calculation loop + processed it two entries at a time and could index one past the end. + Now safely drops a trailing unpaired entry instead of crashing. +- Folder name in this README (`Base-Editors (BE)/`) didn't match the + actual directory (`Base-Editors/`) -- corrected. +- A compiled `__pycache__/*.pyc` file was committed to git -- removed and + added to `.gitignore`. + +### Performance + +- `prepare_model_inputs()` rebuilt the same list of sequence 2-mers from + scratch 16 times per guide (once per dinucleotide it tested against). + Now built once and reused -- verified byte-identical output. +- `_build_structure_df()` (used once per spacer position, per structural + region -- up to ~160 times per run) converted a DataFrame slice to a + Python dict and summed values in a pure-Python loop. Replaced with a + single vectorized pandas row-sum -- verified byte-identical output. +- Row-wise feature-building loops (`make_fasta_for_rnafold`, + `compute_target_features`) switched from `.iterrows()` to the faster + `.itertuples()`. + +### Still missing (not invented -- see "Missing files" below) + +- `connection_to_matrix.cpp` (C++ source for pipeline Step 6) and the + `models/` directory of trained `.h5` weights aren't present in this + branch's history. `requirements.txt`, `.gitignore`, and a `Makefile` + have been recreated from what the README documents, but the actual + compiled-binary source and trained model weights need to be supplied + by the repository owner -- they can't be safely reconstructed. + ## License MIT License — see `LICENSE` file for details. diff --git a/Snakefile b/Snakefile new file mode 100644 index 0000000..b5248cf --- /dev/null +++ b/Snakefile @@ -0,0 +1,284 @@ +""" +DGD-Cas9 Snakemake workflow +============================= +Wraps the 12-step DGD guide-scoring pipeline. For the "standard" SpCas9 +variant, each step is wired as its own Snakemake rule (using the +Accessory/*.py standalone scripts) so you can see -- and resume from -- +exactly where the pipeline is, and so it can run as far as steps 1-11 +(guide scanning through feature engineering) without needing trained +model weights. The "broad_pam" and "base_editor" variants are wired as a +single rule each, matching their all-in-one CLI usage. + +Two things this repo does NOT include (see README "Known Issues"): + - connection_to_matrix.cpp (C++ source for Step 6) -- build it with + `make` once you have the source; until then, Step 6 fails clearly. + - Trained .h5 model files (models/) -- needed only for the final + scoring step, controlled by config['run_scoring']. + +Usage: + snakemake --cores 1 -n # preview the DAG + snakemake --cores 1 # run with deps already on PATH + snakemake --cores 1 --use-conda # run using envs/environment.yaml +""" + +configfile: "config/config.yaml" + +FASTA = config["fasta"] +VARIANT = config.get("variant", "standard") +MODELS_DIR = config.get("models_dir", "models") +RUN_SCORING = config.get("run_scoring", False) + +VALID_VARIANTS = ("standard", "broad_pam", "base_editor") +if VARIANT not in VALID_VARIANTS: + raise ValueError(f"config['variant'] must be one of {VALID_VARIANTS}, got '{VARIANT}'") + + +def _standard_targets(): + targets = ["Deep_learning_file.csv"] + if RUN_SCORING: + targets.append("DGD.csv") + return targets + + +rule all: + input: + _standard_targets() if VARIANT == "standard" + else ["DGDVar.csv"] if VARIANT == "broad_pam" + else ["DGDBE.csv"], + + +# =========================================================================== +# Standard SpCas9 variant (DGD.py) -- granular, step-by-step rules +# =========================================================================== + +rule scan_guides: + """Step 1: scan the input FASTA for NGG-PAM guide candidates.""" + input: + fasta=FASTA, + output: + "Structure_file.csv", + log: + "logs/01_scan_guides.log", + conda: + "envs/environment.yaml" + shell: + "mkdir -p logs && " + "python -c \"from DGD import scan_guides; scan_guides('{input.fasta}')\" > {log} 2>&1" + + +rule build_rnafold_fasta: + """Step 2: append the sgRNA scaffold to each guide for RNAfold.""" + input: + "Structure_file.csv", + output: + fasta="Structure_Connection.fa", + csv="Structure_Connection.csv", + log: + "logs/02_fastamaker.log", + conda: + "envs/environment.yaml" + shell: + "python Accessory/fastamaker.py --input {input} " + "--fasta-out {output.fasta} --csv-out {output.csv} > {log} 2>&1" + + +rule target_features: + """Step 3: sequence/thermodynamic/compositional features per guide.""" + input: + "Structure_file.csv", + output: + "Target_sequence_feature.csv", + log: + "logs/03_targetsequence.log", + conda: + "envs/environment.yaml" + shell: + "python Accessory/targetsequence.py --input {input} --output {output} > {log} 2>&1" + + +rule rnafold: + """Step 4 (external tools): RNAfold + b2ct secondary-structure prediction.""" + input: + "Structure_Connection.fa", + output: + out="Structure_Connection.out", + outs="Structure_Connection.outs", + log: + "logs/04_rnafold.log", + conda: + "envs/environment.yaml" + shell: + """ + mkdir -p logs + (RNAfold -j0 --noPS < {input} > {output.out}) > {log} 2>&1 + (b2ct < {output.out} > {output.outs}) >> {log} 2>&1 + """ + + +rule parse_rnafold: + """Step 5: parse b2ct output into a tabular connection matrix.""" + input: + "Structure_Connection.outs", + output: + "Structure_out.txt", + log: + "logs/05_connectstructure.log", + conda: + "envs/environment.yaml" + shell: + "python Accessory/Connectstructure.py --input {input} --output {output} > {log} 2>&1" + + +rule connection_to_matrix: + """ + Step 6 (external C++ binary): build the full base-pair connection matrix. + + connection_to_matrix.cpp is not included in this repository (see + README "Known Issues") -- this rule will fail with a clear message + until you supply the source and run `make`. + """ + input: + "Structure_out.txt", + output: + "Structure_basepairs.csv", + log: + "logs/06_connection_to_matrix.log", + shell: + """ + mkdir -p logs + if [ ! -x ./connection_to_matrix ]; then + echo "ERROR: ./connection_to_matrix binary not found or not executable." | tee {log} + echo "This is built from connection_to_matrix.cpp via 'make' -- that" | tee -a {log} + echo "source file is not present in this branch's history (see README" | tee -a {log} + echo "'Known Issues'). Supply it and run 'make' before this rule can run." | tee -a {log} + exit 1 + fi + ./connection_to_matrix {input} 102 > {output} 2> {log} + """ + + +rule spacer_scaffold_pairs: + """Step 7: extract spacer-to-scaffold base-pair columns.""" + input: + "Structure_basepairs.csv", + output: + "spacer_scaffold_basepairs.csv", + log: + "logs/07_spacerscaffold.log", + conda: + "envs/environment.yaml" + shell: + "python Accessory/spacerscaffold.py --input {input} --output {output} > {log} 2>&1" + + +rule connection_frequency: + """Step 8: pivot spacer-scaffold connections, extract position labels.""" + input: + "spacer_scaffold_basepairs.csv", + output: + "spacer_scaffold_feature.csv", + log: + "logs/08_spacerconnectionfrequency.log", + conda: + "envs/environment.yaml" + shell: + "python Accessory/spacerconnectionfrequency.py --input {input} --output {output} > {log} 2>&1" + + +rule annotate_regions: + """Step 9: label each connection with its structural region (R/TL/AR/...).""" + input: + "spacer_scaffold_feature.csv", + output: + "Structural_annotation.csv", + log: + "logs/09_serialconnection.log", + conda: + "envs/environment.yaml" + shell: + "python Accessory/serialconnection.py --input {input} --output {output} > {log} 2>&1" + + +rule build_features: + """Step 10: combine sequence features with structural connectivity features.""" + input: + bp="Structure_basepairs.csv", + seq="Target_sequence_feature.csv", + annot="Structural_annotation.csv", + output: + "Feature_Data_Spacer_Scaffold.csv", + log: + "logs/10_featuremaker.log", + conda: + "envs/environment.yaml" + shell: + "python Accessory/featuremaker.py --bp-input {input.bp} --seq-input {input.seq} " + "--annot-input {input.annot} --output {output} > {log} 2>&1" + + +rule final_features: + """Step 11: monomer/dimer counts + stacking energy -> final model input.""" + input: + rnafold_out="Structure_Connection.out", + features="Feature_Data_Spacer_Scaffold.csv", + output: + "Deep_learning_file.csv", + log: + "logs/11_finalfeatures.log", + conda: + "envs/environment.yaml" + shell: + "python Accessory/finalfeatures.py --rnafold-out {input.rnafold_out} " + "--feature-input {input.features} --output {output} > {log} 2>&1" + + +rule score_standard: + """Step 12: score guides with the trained DGD CNN ensemble.""" + input: + features="Deep_learning_file.csv", + struct="Structure_file.csv", + output: + "DGD.csv", + log: + "logs/12_score_deep.log", + conda: + "envs/environment.yaml" + shell: + "python Accessory/score_deep.py --models {MODELS_DIR} --input {input.features} " + "--struct-input {input.struct} --output {output} > {log} 2>&1" + + +# =========================================================================== +# Broad-PAM and Base-Editor variants -- single-rule, full-CLI invocation +# (these scripts don't have per-step Accessory wrappers, so we call the +# documented end-to-end CLI directly, same as a user would by hand) +# =========================================================================== + +rule score_broad_pam: + """Variants/DGDVar.py: full pipeline for 9 broad-PAM Cas9 variants.""" + input: + fasta=FASTA, + output: + "DGDVar.csv", + log: + "logs/dgdvar_full_run.log", + conda: + "envs/environment.yaml" + shell: + "python Variants/DGDVar.py {input.fasta} --output {output} " + "--models {MODELS_DIR} > {log} 2>&1" + + +rule score_base_editor: + """Base-Editors/DGDbaseeditor.py: full pipeline for ABE/CBE base editors.""" + input: + fasta=FASTA, + output: + "DGDBE.csv", + log: + "logs/dgdbe_full_run.log", + conda: + "envs/environment.yaml" + shell: + "python Base-Editors/DGDbaseeditor.py {input.fasta} --output {output} " + "--models {MODELS_DIR} > {log} 2>&1" diff --git a/Variants/DGDVar.py b/Variants/DGDVar.py index 8171db7..f83037c 100644 --- a/Variants/DGDVar.py +++ b/Variants/DGDVar.py @@ -547,8 +547,15 @@ def _count_gc_content(sequence: str, partner_index: list) -> int: def _calculate_free_energy(dimers: list) -> float: stacking = return_stacking_model() energy = 0.0 - for i in range(0, len(dimers), 2): - energy += stacking[dimers[i]][dimers[i + 1][::-1]] + # See KNOWN ISSUE note in stacking_model.return_stacking_model(): the + # table only covers 4 of 16 dinucleotide contexts, and `dimers` is not + # guaranteed to have an even length. Fall back to 0.0 / skip a trailing + # unpaired entry instead of crashing. + for i in range(0, len(dimers) - 1, 2): + context = stacking.get(dimers[i]) + if context is None: + continue + energy += context.get(dimers[i + 1][::-1], 0.0) return energy @@ -576,7 +583,7 @@ def compute_final_features( p_idx = adjust_partner_index(return_partner_index(structure)) p_base = return_partner_base(p_idx, sequence) _, count = return_connection_length(p_idx) - dimers = return_dimers_array(sequence, p_idx) + dimers, _anti_seq = return_dimers_array(sequence, p_idx) # BUG FIX: was assigned the raw (dimers, anti_seq) tuple gc = _count_gc_content(sequence, p_idx) / count if count else 0.0 eng = _calculate_free_energy(dimers) diff --git a/__pycache__/DGDbaseeditor.cpython-310.pyc b/__pycache__/DGDbaseeditor.cpython-310.pyc deleted file mode 100644 index 8d062c465ebd690323c7846d9819a37c35404cb7..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 26674 zcmb__Yj7ONec#UP>=TPOhzAMshCK3)xC6k)y(x+!2!aqL9!L@3oesM5YKd8Z3u1S{ zp4sDZJ+RN_NyL@RSeBf~u_OzWY}qW^@iTVpD2^k=Nmb%h5?7_l8CT_loeySQmR*(| z$JB|;@AvQ77XV87V6djQr@N=0|Ng(bS#EDnNAUL}pZlf4pDjcpf5w;Qzc##_#N+;M zG!nTKaU$~($8gLQW8RQobKaES=zLUut$9m+WAicmMpxp+#C)QdoKG73&RR(o)AMPG z$5z^knfXkyeZIZeG2bC+=1S*$mwewn-y^@h^S$_uuk;mn%=b%~#L7T%aDK>$OJOs6;BofppJX1Mi7ojf&Qr)g#XLNdTk@ZFMv;Ho8O!%9PdVcf z9%Anq!yFuTj=W>dp9$-KiFM$aUbLK}DD@1UXPf$)cAj&dM|~f3UXVI6nD?`b=1>IX zUvyqVO&@ZOA?JPhvuJ_2AC}y>bKH49a?d#@Hsx5(2b`10dAaGk4==`@3Fj12XTtL5 zQU9cK8YvghJI>@Zo{wNfrqJ4%rXD7pY3GCJ;YFF}{(MKiYk5{4fq=!R^!2RsA=EYJ zd>CJzW3KbEg#8RLD>Ke{VNc%F4GprwzXEZ6R0Q)_8w)|o?VucR)dE;+B)Vdxc< z(J=O^^HJo#>b&8+=FFqckLE|6k6krCZKyxVfAnSqFm2AtGx7tK6>Q-2shpc1n{ol{~MccQRua|Pg zlUHGbRZ8glK) zi)UvqT$EWJ3uoC~x;QgoFXT#2!O3|jb1>)ROI(kpwQyy%CQeO_+9;S_xZ#!EZHr7| zsqG5#mSk-wXJ+ht$;}rpujK6_dbQ$?r_&$K->{v+;$mLqOABamQDKAt1wh5Iv~7>s z)3Wl{3niy~9V2+wE_wFBQrWYgIXof-E=;_9%$}K^w#%iJ8}`B3)ybUu(iqy7+}VZP z;$nFPgV%#|!*{#9>J^Hdmd&u{%6FW$6&I6vYKp_OF_?Pqb&2y8Z?)nbV`pVD@XJCA z<&{dYgh~43XP9+rYJ9*nESp8*<)j6>41&3 ztgLvhG;@z+bY{d1+RiI1d)_%+Snz^2C%>3O83!*TMi5)Lk}H%(qJbs3*13Y~1+lrx z>PjBDZ7=6nJgnJVzPO6lxL&EMIX*rT5BhS>GGHsanp0jWud-N`{9D1G%6k=6$}W^k zrTl_dD3`J;`O=biCFl)voI(*GpUo+iyHQUp=G?1!XH#OWW=>7e7k*w#l;$`(B)IwP zV$St)L66FFKjqPU5ktUk;UIT}pQ6OVRSYL9;|MNJ&*5)kdSZH#;lu>O$%zT1B951d zX~~$IoMh4*vL@%|SdL|oH-QYqCa3X+9R3BF?Cj*knKS3loz9+`0zcFG2Mic8Gl=9P zU|}NjCYTxn%!@g1Tzbo|$ixI9LRy zV?HDCPQ<(B+a=zOc+Y%?#Cs9%o9~qP4#fNCyCgn<_~3lE#D@^yIo~7kU5F3Q_ey*> z;(KNzBaa2`7eefOrm&KCtM8k{Osm{V!NneOf!sCL&aJH2g_4(7MIaR*ktUpu&+OHl zcf}pA_RM;!vH+;TDk7P{8+6RhUA#0ocj@9(_RQIHQ`Lb+86d?Fvammb?nc4M^D{Hq zGZWRpZOgLOo=t@(XJ2_lY3!#hrO#iQs}5~jH%hxf@21L8e74%%D33zp-dkSK-YCZ! zn?^=BU(FSi?BC}4PEE{CU6?p~akjeKUCk|Ev(->nw)LBy*$WesQx}_n@IrOh=1&&$ zIZn3B^T6g0&P+{k)*A!Otz>g0;8+l{a!J0~Q5(R-x$MNu%z5mMv*%~3yXJCgDeq<7 z{Iv@4cec5m&bf(;(^GTVu=?4meMa{pdphUkva@XK;VtDFDNF=akeA6nsA(*G7JWk;u2JTKffxRncw%_scqZ|+_K>S~ra!HJlXiOS<~f#YzdZ*3Q2L>;X;&q_kKB?H#pMt{`j#3V>+mr<4MF2okyo zmWWV2HxiX`1RV>xRY59Peh|~Z@=%%LfT|snx;FVJ?^I3u;C*96jWW*!12K{a(r)GQ zsw#sGaf7I;l!A0PBj5#sxa)zUDomt@`IA>P5@eeDG!j#XQD5tvk~#;`qPx^sY4tQp z0ou?8sTK!$9w}}If*vz%^ckkH$FLq+57Y^kObew0)dfjy(?gS&1ujwTuNM(kAWJof zoI+&H#1r)*e&nWcI`YPsUN_gQwU}r4v75$P+%f#P6^TZ!8fvFwUNzjuONJ5g%#ElM z^({YsCu&5lB^_%mk&5_OPABds?wD&yKj|c{n(9dD8N`wUk+p=A@{>|iL{0c{Cp{Rs z8c_=+$ipQYWm1EYHI#7DD4X_Ue#%dO0ut6*8$rNn^HUP00S>XX40D}~-?jv>@iWIG z*PdwVL3>LNROz{T53G6*eiS|EkRE(hda%z=O6wh)TlZpq8|I>Y?EH6&)IQ1`r}8}y8Nz< zjMIMybL@ArJXXZUigd5_Fvsu4EDZQPjx`hsW663fU5^b)Z0L(=)^qK@zY3ylMXbn8 ztW_y_B;vkz9qIh{8Xl}=yR!@H-FegS4X&rN6I*S#LOxO0%vmqLY9Bdf0}EsC*=HtZ z=O%0-TH%VxRfD|_rLShW;OMhs;KManzJ-Av0bg(9Z$d4(n&stN1BQ>;XUXa|w|4Nz z;lqFW_GgbA9zA^cFj@Q>MR8FXzvf>PRM*OS=j zPtRU0tU{(*vJZf_*G#H?z+ufU12RY!-wAxv@sW$)6@4J;_n`@$UTyRd)#&R zjbrzXrvr0x@*dIuePgVSH}Q(O*Z6Z5zitMB4~`G&8h||(AUdN?Bkf)kfnWtNtj1fL z@k8?XO0Ia>$(^{xO4w*70%ici%gR9A@#8Dyh1`mJ0ufyjym9nu{zi}pdE&shem_?6 z7RO$^YXnJ9=(2;sS9hO1Jv!@whaa5+dIjYqmYof0c>2CEs?MRcz(~(jY(!pFW$~6RfJ@8D z%R7NpUPT45>k4^+MOk{JSA75_)JYag>rrPHuRHZP#);4K*RlmCNY-E`NHB$+iLRF2 zz=U`cFF_DcAPA*=$vHZ6Vd@v)Wk#wDZRzm*92H!@<0m&vDzQ?s(ge(K(Ps8POx-9 zzJZ}Up%l=UhVO5Hct?D*{(fdH1}b1UF~MNb4HFnI?j!`0$Nf04+4x$*Pi#c*L_jN& zcMPEvDd2g*{Q&;B)5g+XOfa#9x`6G0>)L(bey9D5x@NQ?smUiI8wsZa_>ubC1h!{q_5U@H<5B1URsole?OoA&kXCI8mR=MLwpX_E zxrHm?COHVmE0D`_fnJXu*&ZdfuHUgx*HOn`HjrBb`_&h>j{7!YYv^vM>1ccr!a%(| zapm=bcf~$?T67wV707C}O`|^x45cmtL=+L1`Urw*hRC9Z8w797s7uThgviZ&k4))( zLtF%+|`wWr}|ic ztf%@of{`}mGVU>`APCIGD|m82oWn&tCRGF$f8-Ojz+y4M9-DV=e=&D8R3v1fy32y^ zPPWW3tM&s z2^0pn5#4t%T}#+TsKdHL3V_*;CcoKM$cv5)Mc-rt}@U-LUz|)DR3r{zm9z4Bx z`ta<)(~oBW&mf*5JUj91!ZVC#H=aFs9>Zhf+3Q7taN23m0K(~7+vo42lz{Pe2^x7E z+qkP`lmArnM)e;Dl6kPU-`~Fx1Coh5-9R!2`~yx8kjxYK?Zxkt{*z?7Sl>qCPUPLF z|D^mHw+zg<>FmI)Br(UmWQoBR#r-FaM{Y%K8OZAo^9I7agyb2wqR1O`hH5LdcF;?b z?Q?cIyWWZT2RGWt%-u9@0y&Z?^ppHf%49a$*P}2G4Fmo4NM1_kts5I1&h9(-#_woK zfhh?oolPl^NlI5!iY+PKI>qbU=wb`_s+}p$UcU#u-RJireBAFtxZmG_@POZs@Ckna z;gkL#!h`-0!V!Nb!l(RQ2%q+c5sv!15svwL5RUtgAw1;U2oL*v7tO^ktU(9Z){(WR zyd+2Ak4WoJvGq?I*K+Goj10Y&U-T8Q_cZ6if7de%n$QT%4s#+>erUVqdd8;aa9_$$5mDthu|md08@ zpuqPT%ufeqw#)DKd;DI%&)?zq`vd-PY0% zru$s;&hZ8}hMea^64wP1_d@j#8YJ$SV>X0iaO8Oi)%me2P_EZl#ZkmORbGXyVo~K` zjsZ?x!prm|q;qKbN}gB&?75HaS~aM{-DqHr9}bLM zU|a@s9xYWCbf!q|>EqCPdXU^#Pkap(e~Vn~K1df#!7OIMYEml}_O^vjuYK3jRen!G zN}~HMiWt(>#2R$PVFyii@2;V4B2zkY*Hr6>2S%aMX;CYvA4Br*-yyPtm_-w)>VT5POD3{ zM5eWC%>eob5hjThH~q+ws&-5bZ&HdbaHfJVsv6| zbYgmRVsdnHVssK|h#+qQ-;PdCk512x&Y|$+=-l+^+?;w7jW*c8CU#Frh1uMb1)abV->^L?fI&MHO1^<(Q2w9B!05Z*RbB4Y71q5NYg^>=7D!Ltq01lZ= z{UVFBcK9!l6n3~vI{d`LdtFFt=`~RW6|&N68<7IBT4oKCx*pj;5APV5a`b*zO zdTUjk_$xM)L}0#cMW|4%KJ~p(bnRm`i4Qfx6~rhQO;FKKGvj9%e3`+|G58__ZdUcP z4DK-a0)wAtK=4u<3_j1`s|>!v;7bf95ZoguYJn>XflG^0lpC}trPs%1wu$*3CI?%0vlVFe!%%z1Ov9ypSW((qo0JjBk^ro1@%3mnE=qkcF38y8T z3}Q^5&_8qfXIlU0JYZcGnGD}_9#Ah6ftBS)=S}N8pk^ilYs-($o6~tfkUj_f!wlZ$ zK%>|<)Frt=hVZqL2M-`gCD;gKvg4uq;sBBuklu*@oF5ujgB15lmTMOnZ>louBq)zq z>X#W8Jn<>UevQG?2!eq}0P(LP)BP<3$+jM&$LfKysK+!P^hC|?WGwSNNU~}3VKN5g z6Vy_`iI8Ep{_8S(jNymLC?(t^R8A;q#*Oqt9ri^tM%;)$G|lvb;b;%a>VJb#%EE(2 z=Dq%SpScIJaX-GJ@{6jgMl06+WU>q2AQj(snKnx9L_haD=!z*KnIJ805L?X~+#lS^ z+?B2cbl@0ET7IU)Y*F&6NxG0UqfVSsE~L!_wh47uYPlVTi%uF0Zb!lErE8lg^w>PT?M%}VPUGpu>kB{oBQ z+-Cmc9*Ksq9#9+8;{^SPmsJ5c1 zZTP8#L$9IJUX4x#hyD($iNh2EV`KxDYJ~U=I5m;w9hEC3N4cD^YD;1!6w)xX7R&HiX(+YQilj=NWWZqrmYY+#plzX|;1}g(q0`qm=WUkw zI7>vm;wn(*;%4gCx80=g#hSQ8D?4Kr_Pth1Qyrbg8}P*2=KbM^@l+S~JO5GDtYY}4 zMUzEP;3GSb1_B}0kXE5v!cK(d3fe5k3N=@;B})DS=EY?M-SWI_?&(8h?tX( zid_iKVB;;2O@f(Vb9h{WobY}jSbb&AVE#1@ZWA;6;cKdc@RVhNJ4*@Z)Ui8PsfT6T}_ z(qE%42rP-v3GEDoI)VNDNDpw;B5mg9f--4RD3p z6%&wk6dbHWhL&)q3DxfJI+o)ceHRL+V{w2Y%-whfp2T7wO_YYqRa1$~HYMPQJ}}6 zZ)`x|y<_}@2n0Z$Ms~yy*6^;D1VpCV4nRM`MWwu7JzV2C8r{s+ySW`s*Cx(t?BMIp z(<5hyvB8qV&L|93c0!#&tvY^6wV|2HUL67dRi90{7ux2l4G58M=+P;O_({c0kChz%+*N&+m~Po z(_Uer9Zaxrjk$2OX{lz?g8PojFTwt0ujDS%)mQald?N|Cotks!?(M~J?hz~DJvw{t z$USNe?j5Tq9KAK>_#F7eDgqMv)D5lApD3ckB%>sNE);=KdRJ`%p4yLgXw<9nRC zAT?VvlLqF@>>TkUc+3^})Sc0gw-2R4>=>3r?;8$qb=;V0g1s~-d#M6prmz6JcNsFv z<-B)2PXjo32Dp%Y7obvPUty=>#02I_IWVtm*5U;_!}cLa)Ow%Q6Y72xUyG|yI@SZw z83v;-0Ca-2!idIM&pL}idjoW{m*KIrs`874x3*Y>>mc?=mf*`Ok{DwAS1UZ0O zoVc0(C`bd&1@LY63K+M&S*pF{-*CJb6HC_Lc%rd2ek z8kLfFrKGLChTyK&cd5T%!e25Fs#+(l>erEWon5(&NcG9=WfGJoSaBAckaAreBaAjS zfR%fz06ugME}kn{^&ILM=@LnvoV2F~SusbZ1`rG41-R`iPwimxJ4g-^fPPjg`T)Hv zgim1Ch2wyu1z2H4Cj(e)Kfn)9GMGRRr0O^Wei?Zet{-ZVwI;!~h+Afm1c}o!F|ln) zuq}42qZ>WWR$AvH=-p-p3TVt7M&ZG@cB`<(>7oa6Y8sNMmVwm;b}4yz6^%eVAyI)h ziA39E59*f%3-#3gQizb!a;@F3BdB5Kck2o)I!{2A1>4TMU@mF#U5kZQMr^5fl9WYZ z^EJUAYfCWMk}#ES72uMiYjCSr2Ima<)qqO{T)x&4EWd6#ZDeo6f}Cl52MbbyTzCC@ z_*gE-AX$bMS77<5$ed(OxZ!&Is5m|JVw?0;UlS~U@dzQOch>%Xjc1A*t0U~*Y9GTQB4&(_I z8gK`9(AF*YDSRWiHy$x9>pM_yz%?NLu)0t0znXq-yS=w%XoMgc8YNQ^K-y@WqLR60 z_Zf=Rcs4b~?n`Ip*l&Z*sUfBh1#;?!M6CD>3W&0HPvpCFtZu>*q!7d)p63`Qe3EFyQI7N$G%mGbqFqBqO>Mc{q7p466RTWU@rW&v6d z;xCHlL|#EmV0tqhPbp!FFS2*D3{J{7#2N&fkVcc;1a1YlD8JNws{eq3G9XBXm-yu? z3`otjwA9fGLqCiRm$R8pq=B!J^gprU@9jxIdj`!KT*n^5sqBFjjlZ`u4*rqY4!Q^U zZffNjUh99{MRn!bpW+4d9D9g#TPS?=$1}8tZWy4`QK+=67M0dui%JW^yqI%d21z5U z3B=V8lmoTdi+&7<%M#ZMc!90O0n~7?0QKG*mZIPAS*8hF8UsH7wU#K-Y9&%AQEHTc z^Eed+Kz*#0B^!CQ9Q*=ZHcFs1+6>yN=e0@lm{az~E<>?!HScLJSJ8w0U;)R& zK86*PGYq4379daSH|^b|t8c80Z`MnRm!yzx*+QbT7rExXaj@puRK;Wmsl2d?I1QsU zN#9^%T%lHYlE;)dND610P|E~|cX}ilO z(U$rIg21?Rua7UM8Hjepl5VS23?MLy+nHp84y|g(`3{`J%&TAFjAT$j2WB8+j6&^& z8Sq#Eta2#ewqP|9_xh-*BhpGP2V5{XTX4PLhM^<3UpLq0xc5TX##LiK};54zA#C@jiO>xqk#aZ$WU(~r3GGyp$RjbK>9-RZ;6CQQ zH$Q=0J=AQUjS4i)Xp^aoE}UEf796v~LzsHrg-_uO|3{9JfMULVI8q(S!k2vk1_>C? zp!UgUHSXPRayfcIwBYTiO=(QCwr?96prs=}h{AfB`?qi44{E)L8~Y8lhBqjZv>PZ* z7D0Q+IfQZmXT0`Dl*Fl<#N$%@)<?YY7KNiJD4AJn$DxxH-psWQzg@lP(Vq=(e)}$-q>8R}`knBt z>$qtwneaBliBTL!krFt5qElhl0txO86`yh4546&SUS|9*2;c2%-4LI$-fX(`)F&pOB+ZB7Q~!g|4yw}092d$T=|LVEpP zNsbwTWoj zFm<_qyy3~A>paC;WPeB?TaqdTsr@$)X?RL}@7jeIpA$7gker^(!uzOD4l>h)%a`5k zlyFe$2iRa1!N`DS-MeR7&wpkoaISY#(&W~p>8(kZaB6qchwug1nl`yj+VnPQ+tdgD zfi3mHnP6+$HuX(zQy+W|wtPFiO?~i2*jjFz`Yvr#-=%HpyR=Pxmo5c^o9C~-h{IcQ znidtX`*EyjD0A%s)MzcLK+INwYd~|G_fbyUmk9VO&7hM(8w2uZLhb|;3wI_AhA8p1 zhaxRn>MY@j>Lv9aU&Z?5@O6`^s!#n8Qz=YoNlmm!u&4wHZWzjJaX4;d;dRo+U-F7z zbTp@Zja54cpcSfVVS^YnMb2BQi#7fV3kdV{hb)x{zYe0#;t>W%sngN|=%(Ek^l!g= z{}ZYe->77qwprS0dxTx3#;OaNs~-4Fbl`Mo519HvSZvi-@I?zUYDZTn>#E_^4)~fwQju$jil2Jy4wA2 z+j`=bb*l}Y>OJ-Hy}I2OQ2Lr$w;)RP-GLI~-PEnvt@y3Pt>mrLt@L6Xrvjv>USb_u zSGOq@Kh%Fg+SZll0j{fbIdJUL1 z_c)I^_H7|yIAojlGcWKkZejzz@gi0y5z~iCzwEVvFMu*0he!9~EX6)2b;pu42-$vY zF$zlkc=favuD3s;ul-Bp$_+%JR<-5o4_qkI>Lf1Q`OGG{x&@eP;s7bL>IA5d8262R zz&tN8p_`j<6B+w4Wc(e|KZ8hhPnMq8(uXF(*wB5+Nv6{%`N&q?X7!w==Qz0rT?98o zsGc($C!)aUOV?-{M4VJTFt$)~cJlAp$=48UX_`D`FP?7=BM{ACn(_~=xY{4QX3JGS zjsjxDlBs}JkUDhnebLh*gq@B31KJS*aX6wri}y#islLi?{bWnM>gSo(;1pQleJ0&S z@I7k!r{A-tZ-1XP{j={`(|fg=g6PD=5tuZ7nxFg(gD*4qIR>9+Abyu$VeG36zQACE z!6^nm&pNAg}hB2wa`(~1tcTGA|5}oo8KBlTLQ^RDju~$7@pq^7v;N# zwtk7=2tPAbzP_ooG+Gl79^ASRwdqj5&Iy*${RU&KP5ov)t-c-FIK{t(^hfr%Q4`eC zvBi}#)WeUeEmNya+}E~LNkaS#9+!AVq-%JY>nTuaCsOKNhu&-`qHqggsa>l!6z*Cp zCDyTONFR2h5cRB3tq}wIjW_uvf09G|1cLj<2bwh-R$yXP&Fli6B|OcGN9ntriZzno?uSz^GOhAUueOtS$~$=z)bA8X15d-jda{+Iop(K7VKH>U6Vf zcp?W|*M}2H+ncm`5)o;0m=LAU8Pr|&yXms;7#o;aoJ6DA5zbrRW|x{1VI89l{h>|@C_#YrX=ND zeZdxt5bl-SYJQ|$JF@&fXV_;sx-182io~$4!t{z$_6{yPNpsJaalTo>8N5)+!anRc zSy4H|p$!>uw=CsxE=>J7YvHJI=KptD?Z0O5H3t8N0dbo8B?iCD;9oGH-c|i7gI{Ct ziwyoH1DQf@NnFBHr_z6o9QPy=lcC;v7<%gtW6;zN#x#8-4H$CG{YEm%gSh+guLsUF zX*2o2#wmF?DnHm0+h?l(jGBU!Fd^jVwQBQS)ZF3%VIhI^mayeMi$L>Li{N)fSt%Sd zd>f(BSqCg>ei9m~$ zz_3ue2`>-(1Y9#HtHW`l)lZRoAc;o+p3^P>-%Wtuz)@cMym6!9?oneSwr;LRmRjiz zuzSi|)OVj@nSsO{{h}VCG7=x}N*$i03w4%d5iSUvEcJ(6s+SpSh4wv!c6N@?4)4*2 z)yJ~9GtkowY{zW1;1AR)b3Th-H;upt*`y;MX49OhMlp*(ytuQ%ad~DPadWYxyU^&z z48F734^C7Xw_`k_2l{pltqYz+%g_t|1HRC+kRUfm0Fhl}6|o>d1C3#G7HHU$YxZ;x z{)lrTklWVO;&5$7I?NsbVTYyu80GG#>em1UodS)n+~7r+3^H1-+{VX1Mga=S1fx-@MjnNmIgn~`W9SzN^ zd5oV-%k7nQ<+FocQ*bJ4m7aljV} z$SVAfypn?GUfBw{sf3Y0GmqzwPGGBqnf6mD@rbO8Z)PSMXGL4J|O z-=uVy9}6O=V*_j^PF#3r%K>i2KD3tZ)YPNP95#lHesd?CE&m*S35@X|9fAuMzJ)4~ zB2XzMRN*HPabQ`2JT1OE&@jNe2O+M~LP&2N+_ZpL(&G&S(36MnZ60q9$aogQb}xpo z!;2&A^x$IbbcHt(bfa`~(Okm)4B}$PljgjpE4h)P!y8QAF3g`xyW@E*Nw=Q;AKzxhMffzY55~| z6%s^rPwrjpNX87>HuJB>~%cD>m1QVdwHb}K8q_Vihw zN6vuy3v3LI`ONsC?LM>{qeFAN+Y+~B9kow1rMKKpIT!S8SB=bo`b#2&#vX1RKNXki znrKjpe`J@3Njo$D=1MOvB8L1b4Ns3Q|mtPY0HoJqm zmEE<{M-EHrI-J%lFqQrF5_*p{nr~H&H~M)*gZwNI7{?{I^pHjYpVwJMcgw1gca?pa zK}1~Qpp?383I~P#coM1jA+FO~)Kve6ai;5e7kQZb^kbN9F*utMz}my=p4wfnyo6B@ z%eZ|}b3n(0N~Gnc8#%HkCxzt9hWN0n6KsRtNTQl z_%4lS5|8`G2*Ax??hryw1Xl!Z6rPt_R}Sl>yopyPvHNW+OaYH0zLkv>InxLIbFVHV zj?ab{qksK!LOqXIvScDm@mfHq|8~@fRNG;(HNZ8&+KQ*$8DJ`$uz-@Rh}YqDdR>dr z8&?h7;e*mWat3fvTTQ**jXr0HG~(6BYxrW#{2`Oy>DYDSdSqnhOywrV*K*smfGtw$ zWMT6ae)ho-_(sC(tn|(69IJknltVMr2GWMb*;1nZo4j*x=$parecHre>2WD~pAtvQ;)Ie^V$B^JOu~_sl z!zVF`vbn0}0ay6LksWYrRZDR1!j0M;V3u)eq2!77UIn{r@2B;oxF5tso{?IyK0$_h z0DLT;1Zk;*?{1I?2Nbjc-xn4p%caG_QjiFT8^lgeor1B^Iy-abyzD;px0qU}5ro6l zylnRwEXYjRJHr#{Q@qnPNWa1p=@M7(<4EXj7Q}D@6PZ!>1{4f2JE1?wz_^xt+kBup z@eMBs=J(kw=9G{mpn*6~n9Pj4rk-T)MBgmh+J}67jsdllIJj9}SpKn=5snf? z43Re6LvQdJN=42?k=az#rQz6WrMP-S3!)fd$;r9eyZ-l>o5BAuFg}~VcF8M1Q$*fC z{XWyXa1jRXD};{n`uJ7c{o%^|s{g`n3(G#u*wAvd)Rdr)v4adq zb>$RoU|ojpOS41$OcbVafHA6#HKnl``zC`S2DCP7!b2O0zG=Cg!=$xF(;;%YnizeX zu?x)I%b3un-(jrGw7ZPa4Nx1JtBldhOiJsxKb&Pgk zXcM7gH{tPX;ygcc-?JV=>VrmV`a2d*7Hvuy-jdS3C8c|Fiis5MUCg#Dwl#n<{P$j? zJ@x{(3dGXi-f8dY+L`Q6^du~tw7|bDeck_rJ!-X~)<@J}+Q`}RReE1~ClAK{fA+aE A`~Uy| diff --git a/config/config.yaml b/config/config.yaml new file mode 100644 index 0000000..aafcaf9 --- /dev/null +++ b/config/config.yaml @@ -0,0 +1,23 @@ +# DGD-Cas9 Snakemake workflow configuration +# -------------------------------------------- +# Input FASTA (100-10,000 nt sequences). A small synthetic demo sequence is +# bundled at demo/input.fa so the workflow can be tried immediately. +fasta: demo/input.fa + +# Which pipeline variant to run: +# standard - DGD.py (SpCas9, NGG PAM) -> DGD.csv +# broad_pam - Variants/DGDVar.py (9 broad-PAM Cas9 variants) -> DGDVar.csv +# base_editor - Base-Editors/DGDbaseeditor.py (ABE/CBE) -> DGDBE.csv +variant: standard + +# Directory containing trained .h5 model files for the chosen variant. +# NOT included in this repository (see README "Known Issues" section) -- +# scoring (the final pipeline step) will fail with a clear error until you +# supply real trained models here. +models_dir: models + +# Whether to run the final scoring step (Step 12 / Step 11 for variants). +# Leave this false to exercise steps 1-11 (guide scanning through feature +# engineering) without needing trained models -- useful for testing the +# workflow plumbing on a machine that doesn't have the models yet. +run_scoring: false diff --git a/demo/input.fa b/demo/input.fa new file mode 100644 index 0000000..8f499c4 --- /dev/null +++ b/demo/input.fa @@ -0,0 +1,7 @@ +>demo_sequence (synthetic, randomly generated -- for pipeline testing only, not real biological data) +GTCAAAGATAACCATACAATACATATGCTAGTATGAAATCCGCCGGTTTAGCGCGAACTGTCCAAGGGCC +AAGCAAGCGCGACGTACATGTTCCCATCCGTCGCGCTTATTTTACTGGATCGGTACCCCCACCATATCTA +GAAATAGAATCTGGGCAACCCCGATAGGCTATGTAGAGGTGTGTTTCTTCGAAGCGTGCGGATATCTGTC +ACGAACTCTCGACCCATTTATCTCGTTAAATCTTAGTGTGGCTAGCCTTACTATTTCAAGCAATTGAACA +ACGTGCCTGTCTCACCGATCATGATGTGTCTACCTTTCCGCTGGAGAGCCACGGAATAAGGATGTCGCTC +GAGATCAGAAGATAGACAGCGTATATGATTGGTGCTGGAGAAATCTCACC diff --git a/envs/environment.yaml b/envs/environment.yaml new file mode 100644 index 0000000..5fa8402 --- /dev/null +++ b/envs/environment.yaml @@ -0,0 +1,14 @@ +name: dgd-cas9 +channels: + - bioconda + - conda-forge +dependencies: + - python>=3.8 + - numpy>=1.21 + - pandas>=1.3 + - biopython>=1.79 + - viennarna>=2.5.0 # provides both the Python `RNA` module and the + # RNAfold / b2ct command-line tools used in Step 4 + - pip + - pip: + - tensorflow>=2.6 diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..907564c --- /dev/null +++ b/requirements.txt @@ -0,0 +1,15 @@ +# DGD-Cas9 dependencies +# Python >= 3.8 required +# +# Recreated from the versions documented in README.md -- the original +# requirements.txt was missing from this branch's history. + +numpy>=1.21 +pandas>=1.3 +biopython>=1.79 +tensorflow>=2.6 + +# ViennaRNA Python bindings (RNA free energy calculations). +# If the pip build fails on your system, install via conda instead: +# conda install -c bioconda viennarna +ViennaRNA>=2.5.0 diff --git a/stacking_model.py b/stacking_model.py index 48c065e..32a21ee 100644 --- a/stacking_model.py +++ b/stacking_model.py @@ -17,6 +17,28 @@ def return_stacking_model() -> dict: Constructs and returns a nested dictionary mapping dinucleotide pairs to their stacking free energies in kcal/mol at 37°C. + **************************************************************** + KNOWN ISSUE -- INCOMPLETE TABLE, NEEDS ATTENTION BEFORE TRUSTING + RESULTS THAT DEPEND ON Spacer_Scaffold_Gibbs_Energy: + **************************************************************** + A standard nearest-neighbour RNA stacking table has 16 top-level + dinucleotide contexts (AA, AU, AG, AC, UA, UU, UG, UC, GA, GU, GG, + GC, CA, CU, CG, CC). This dict only defines 4 of them (AA, AU, AG, + AC). Since real guide/scaffold sequences routinely contain the other + 12 dinucleotides, DGD.py's _calculate_stacking_energy() used to raise + a hard KeyError the first time it looked one up -- i.e. step 11 + (compute_final_features) could not complete for any real input. + + Rather than invent the missing thermodynamic values (they should come + from the original published parameter set this model was trained + against, not be guessed), calculate_stacking_energy() below now falls + back to 0.0 kcal/mol for any missing context and logs a warning. This + keeps the pipeline runnable end-to-end for development/testing, but + the resulting Spacer_Scaffold_Gibbs_Energy feature is only fully + correct for guides whose paired dimers happen to only use AA/AU/AG/AC + contexts. Please supply the complete 16-context table before relying + on this feature for real predictions. + Returns: A dict where keys are dinucleotide strings (e.g., 'AA', 'GC') and values are dicts mapping base-pair contexts to stacking energies.