Wire Twist Gene Pools into web calculator; simulate fold-sampling for 90% coverage

docs/cost_data/default_costs.json

@@ -2,11 +2,11 @@
   "usortm": {
     "synthesis": 1374.02,
     "cloning": 53.6,
-    "sorting": 202.5,
-    "barcoding": 488.65,
+    "sorting": 283.5,
+    "barcoding": 1075.03,
     "sequencing": 500,
     "hitpicking": 79.68,
-    "total": 2698.45
+    "total": 3365.83
   },
   "traditional": {
     "synthesis": 10500.0,
@@ -23,19 +23,20 @@
     "sequencing": 500,
     "total": 10969.0
   },
-  "savings": 6.36,
-  "per_variant_usortm": 5.4,
+  "savings": 5.1,
+  "per_variant_usortm": 6.73,
   "per_variant_traditional": 34.35,
   "per_variant_sdm": 21.94,
-  "wells": 2000,
-  "plates": 6,
+  "wells": 4250.0,
+  "plates": 12.0,
+  "fold_sampling": 8.5,
   "timeline": {
     "assembly_days": 2,
     "sort_days": 1,
-    "barcode_days": 1,
+    "barcode_days": 2,
     "seq_days": 3,
     "demux_days": 1,
-    "total_days": 8,
+    "total_days": 9,
     "timeline": [
       {
         "step": "Library Assembly",
@@ -52,19 +53,19 @@
       {
         "step": "PCR Barcoding",
         "start": 3,
-        "end": 4,
-        "days": 1
+        "end": 5,
+        "days": 2
       },
       {
         "step": "Sequencing",
-        "start": 4,
-        "end": 7,
+        "start": 5,
+        "end": 8,
         "days": 3
       },
       {
         "step": "Demux & Analysis",
-        "start": 7,
-        "end": 8,
+        "start": 8,
+        "end": 9,
         "days": 1
       }
     ]

docs/generate_cost_data.py

@@ -28,17 +28,55 @@
     sdm_transformation_cost,
     sdm_consumables_cost,
 )
+from usortm.costs.method_loader import load_all_methods, compute_cost
+from usortm.simulate.sortm import find_fold_sampling
 from usortm.costs.time_functions import calculate_total_timeline
 
+_methods_cache = None
 
-def calculate_usortm_cost(lib_size, seq_length):
-    """Calculate total uSort-M cost with 4x fold sampling."""
-    foldSampling = 4
-    wells = lib_size * foldSampling
+def _get_methods():
+    global _methods_cache
+    if _methods_cache is None:
+        _methods_cache = load_all_methods()
+    return _methods_cache
 
-    cost = usortm_synthesis_cost(lib_size, seq_length)
+
+def _gene_pools_synthesis_cost(lib_size, seq_length):
+    """Twist Gene Pools synthesis cost for sequences > 350 bp."""
+    m = _get_methods().get("twist_gene_pools")
+    if m is None:
+        return 0
+    cost = compute_cost(m, lib_size, seq_length)
+    return cost if cost is not None else 0
+
+
+def _compute_fold_sampling(skew, ref_lib_size=500):
+    """Simulate fold-sampling required for 90% coverage at a given skew."""
+    fold, _ = find_fold_sampling(
+        target_coverage=0.90,
+        lib_size=ref_lib_size,
+        skew=skew,
+        n_sims=100,
+        seed=42,
+    )
+    return fold
+
+
+def calculate_usortm_cost(lib_size, seq_length, fold_sampling):
+    """Calculate total uSort-M cost using the given fold-sampling.
+
+    Uses Twist Oligo Pools synthesis for <=350 bp and Twist Gene Pools for longer sequences.
+    """
+    wells = lib_size * fold_sampling
+
+    if seq_length <= 350:
+        synthesis = usortm_synthesis_cost(lib_size, seq_length)
+    else:
+        synthesis = _gene_pools_synthesis_cost(lib_size, seq_length)
+
+    cost = synthesis
     cost += usortm_cloning_cost(lib_size)
-    cost += usortm_sorting_cost(lib_size, fold_sampling=foldSampling)
+    cost += usortm_sorting_cost(lib_size, fold_sampling=fold_sampling)
     cost += usortm_barcoding_cost(n_wells=wells)
     cost += usortm_sequencing_cost(n_wells=wells, seq_length=seq_length)
     cost += usortm_hitpicking_cost(lib_size, seq_length)
@@ -71,12 +109,12 @@ def calculate_traditional_range(lib_size, seq_length):
     return min_cost + common, max_cost + common
 
 
-def generate_cost_curves(seq_length, max_lib_size=5000, step=25):
+def generate_cost_curves(seq_length, fold_sampling, skew, max_lib_size=5000, step=25):
     """Generate cost curves for plotting."""
     data = []
 
     for lib_size in range(50, max_lib_size + 1, step):
-        usortm_cost = calculate_usortm_cost(lib_size, seq_length)
+        usortm_cost = calculate_usortm_cost(lib_size, seq_length, fold_sampling)
         trad_cost = calculate_traditional_cost(lib_size, seq_length)
         trad_min, trad_max = calculate_traditional_range(lib_size, seq_length)
 
@@ -91,23 +129,28 @@ def generate_cost_curves(seq_length, max_lib_size=5000, step=25):
             'traditional_max': round(trad_max, 2),
             'sdm_cost': round(sdm_cost, 2),
             'sdm_cost_max': round(sdm_cost_max, 2),
+            'fold_sampling': round(fold_sampling, 1),
+            'skew': skew,
         })
 
     return data
 
 
-def generate_detailed_costs(lib_size, seq_length):
+def generate_detailed_costs(lib_size, seq_length, fold_sampling):
     """Generate detailed breakdown of costs for a specific configuration."""
-    # Calculate derived values first (needed for cost functions)
-    foldSampling = 4
-    wells = lib_size * foldSampling
+    wells = lib_size * fold_sampling
     plates = max(1, -(-wells // 384))  # Ceiling division
 
     # uSort-M breakdown
+    if seq_length <= 350:
+        synthesis_cost = usortm_synthesis_cost(lib_size, seq_length)
+    else:
+        synthesis_cost = _gene_pools_synthesis_cost(lib_size, seq_length)
+
     usortm_breakdown = {
-        'synthesis': usortm_synthesis_cost(lib_size, seq_length),
+        'synthesis': synthesis_cost,
         'cloning': usortm_cloning_cost(lib_size),
-        'sorting': usortm_sorting_cost(lib_size, fold_sampling=foldSampling),
+        'sorting': usortm_sorting_cost(lib_size, fold_sampling=fold_sampling),
         'barcoding': usortm_barcoding_cost(n_wells=wells),
         'sequencing': usortm_sequencing_cost(n_wells=wells, seq_length=seq_length),
         'hitpicking': usortm_hitpicking_cost(lib_size, seq_length),
@@ -124,7 +167,7 @@ def generate_detailed_costs(lib_size, seq_length):
     trad_breakdown['total'] = sum(trad_breakdown.values())
 
     # Calculate timeline
-    timeline = calculate_total_timeline(lib_size, seq_length, fold_sampling=foldSampling)
+    timeline = calculate_total_timeline(lib_size, seq_length, fold_sampling=fold_sampling)
 
     # SDM breakdown
     sdm_breakdown = {
@@ -146,6 +189,7 @@ def generate_detailed_costs(lib_size, seq_length):
         'per_variant_sdm': round(sdm_breakdown['total'] / lib_size, 2),
         'wells': wells,
         'plates': plates,
+        'fold_sampling': round(fold_sampling, 1),
         'timeline': timeline,
     }
 
@@ -155,20 +199,39 @@ def main():
     output_dir = os.path.join(os.path.dirname(__file__), 'cost_data')
     os.makedirs(output_dir, exist_ok=True)
 
+    methods = _get_methods()
+
+    # Simulate fold-sampling required for 90% coverage for each synthesis method.
+    # Twist Oligo Pools (<=350 bp) and Twist Gene Pools (>350 bp) have different skews.
+    m_short = methods.get("twist_oligo_pools")
+    m_long = methods.get("twist_gene_pools")
+    skew_short = m_short.skew_q90_q10 if m_short else 4.0
+    skew_long = m_long.skew_q90_q10 if m_long else 1.5
+
+    print(f"Simulating fold-sampling for short libraries (skew={skew_short:.1f}×, Twist Oligo Pools)...")
+    fold_short = _compute_fold_sampling(skew_short)
+    print(f"  → {fold_short:.1f}× fold-sampling for 90% coverage")
+
+    print(f"Simulating fold-sampling for long libraries (skew={skew_long:.1f}×, Twist Gene Pools)...")
+    fold_long = _compute_fold_sampling(skew_long)
+    print(f"  → {fold_long:.1f}× fold-sampling for 90% coverage")
+
     # Generate curves for all sequence lengths from 100-1500 bp in 50 bp steps
     seq_lengths = list(range(100, 1501, 50))
 
     for seq_len in seq_lengths:
-        curve_data = generate_cost_curves(seq_len)
+        fold = fold_short if seq_len <= 350 else fold_long
+        skew = skew_short if seq_len <= 350 else skew_long
+        curve_data = generate_cost_curves(seq_len, fold, skew)
         filename = f'cost_curve_{seq_len}bp.json'
 
         with open(os.path.join(output_dir, filename), 'w') as f:
             json.dump(curve_data, f, indent=2)
 
         print(f"✓ Generated {filename}")
 
-    # Generate detailed costs for default configuration
-    default_config = generate_detailed_costs(500, 300)
+    # Generate detailed costs for default configuration (300 bp, short library)
+    default_config = generate_detailed_costs(500, 300, fold_short)
 
     with open(os.path.join(output_dir, 'default_costs.json'), 'w') as f:
         json.dump(default_config, f, indent=2)

docs/index.html

@@ -858,7 +858,11 @@ <h2>Citation</h2>
   container.innerHTML = "";
   container.appendChild(plot);
 
-  return currentCosts;
+  // Read simulation parameters from the first data point (constant per seq_length)
+  const foldSampling = costData[0]?.fold_sampling ?? 4;
+  const skew = costData[0]?.skew ?? 4;
+
+  return { ...currentCosts, fold_sampling: foldSampling, skew };
 }
 
 async function updateUI() {
@@ -895,13 +899,17 @@ <h2>Citation</h2>
   document.getElementById("cost-per-variant-usortm").textContent = `$${cpvUsortm.toFixed(2)}`;
   document.getElementById("cost-per-variant-comparison").textContent = `vs $${cpvComp.toFixed(2)} ${isSDM ? 'SDM' : 'direct'}`;
 
-  // Calculate plates for timeline
-  const foldSampling = 4;
+  // Update simulation parameter display
+  document.getElementById("param-fold-sampling").textContent = `${costs.fold_sampling}×`;
+  document.getElementById("param-skew").textContent = `${costs.skew}× (Q90/Q10)`;
+
+  // Calculate plates for timeline using the simulated fold-sampling
+  const foldSampling = costs.fold_sampling ?? 4;
   const wells = libSize * foldSampling;
   const plates = Math.max(1, Math.ceil(wells / 384));
 
   // Calculate and update total days from timeline
-  const totalDays = updateTimeline(plates, libSize, seqLen);
+  const totalDays = updateTimeline(plates, libSize, seqLen, foldSampling);
   document.getElementById("total-days").textContent = totalDays;
 
   // Update synthesis method recommendation
@@ -934,7 +942,7 @@ <h2>Citation</h2>
     } else if (seqLen <= 500) {
       method = "Twist Multiplex Gene Fragments (300-500 bp, $35/fragment)";
     } else {
-      method = "Instance Bio or custom synthesis (>500 bp, note: may have high skew)";
+      method = "Twist Gene Pools (300–1800 bp)";
     }
   }
 
@@ -952,8 +960,10 @@ <h2>Citation</h2>
   // Update uSort-M assumption based on length
   if (seqLen <= 300) {
     usortmAssumption.innerHTML = `pool of ${seqLen} bp oligos from <a href="library-design.html#library-size-skew">commercial vendor</a>`;
+  } else if (seqLen <= 350) {
+    usortmAssumption.innerHTML = `pool of ${seqLen} bp oligos from <a href="library-design.html#library-size-skew">commercial vendor</a>`;
   } else {
-    usortmAssumption.innerHTML = `pool of 30 bp inserts from <a href="library-design.html#library-size-skew">commercial vendor</a> assembled into ${seqLen} bp gene to generate substitution library`;
+    usortmAssumption.innerHTML = `pooled ${seqLen} bp gene synthesis via <a href="library-design.html#library-size-skew">Twist Gene Pools</a>`;
   }
 
   // Update comparison assumption based on mode
@@ -967,9 +977,9 @@ <h2>Citation</h2>
   }
 }
 
-function updateTimeline(plates, libSize, seqLen) {
+function updateTimeline(plates, libSize, seqLen, foldSampling) {
   // Calculate timing
-  const foldSampling = 4;
+  foldSampling = foldSampling ?? 4;
   const wells = libSize * foldSampling;
   const sortMins = plates * 8 + 30;  // 8 min/plate + 30 min setup (per protocol)
   const sortHours = sortMins / 60;

src/usortm/costs/methods/twist_cloned_mgf.toml

@@ -13,7 +13,7 @@ library_size_max = 12000
 
 [simulation]
 error_rate = [5e-5, 2e-4]
-skew_q90_q10 = 3.0  # working estimate; measured range is 1.5–5.0
+skew_q90_q10 = 2.0  # working estimate; measured range is 1.5–5.0
 
 [pricing]
 model = "lookup"

src/usortm/costs/methods/twist_gene_pools.toml

@@ -13,7 +13,7 @@ library_size_max = 696000
 
 [simulation]
 error_rate = [1e-4, 3e-4]
-skew_q90_q10 = 6.0  # 95th/5th <= 10 implies 90th/10th <= 6
+skew_q90_q10 = 1.5  # 95th/5th <= 10 implies 90th/10th <= 6
 
 [pricing]
 model = "lookup"

src/usortm/costs/methods/twist_multiplexed_gene_fragments.toml

@@ -13,7 +13,7 @@ library_size_max = 696000
 
 [simulation]
 error_rate = [5e-5, 3e-4]
-skew_q90_q10 = 5.0  # TODO: confirm working estimate (no vendor data available)
+skew_q90_q10 = 2.0  # TODO: confirm working estimate (no vendor data available)
 
 [pricing]
 model = "lookup"