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[<list>][...], 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 8d062c4..0000000
Binary files a/__pycache__/DGDbaseeditor.cpython-310.pyc and /dev/null differ
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.