upstream/mims-harvard-ToolUniverse/catalogue/skills/tooluniverse-rnaseq-deseq2/SKILL.md

@@ -1,8 +1,8 @@
 ---
 lineage_type: import
-upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-rnaseq-deseq2/SKILL.md
-upstream_sha: e2520a96
-imported_at: 2026-06-26
+upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/aafb6fb6/skills/tooluniverse-rnaseq-deseq2/SKILL.md
+upstream_sha: aafb6fb6
+imported_at: 2026-11-03
 prompt_class: catalogue
 upstream_changes: accepted
 name: tooluniverse-rnaseq-deseq2
@@ -36,23 +36,23 @@ genes can use the unshrunken value. For multi-contrast runs it auto-emits
 3-way Venn region sizes and percentage-of-X interpretations.
 
 ```bash
-# Single-factor sex DE on a CD4/CD8 subset, with FAM138A LFC
+# Single-factor DE on a cell-type subset, reporting one gene's LFC
 python scripts/r_deseq2_wrapper.py \
     --counts <data-folder>/counts.csv \
     --metadata <data-folder>/meta.csv \
     --design "~sex" --contrast "sex,M,F" \
-    --subset-col celltype --subset-values "CD4,CD8" \
+    --subset-col celltype --subset-values "<type-1>,<type-2>" \
     --min-row-sum 10 --shrink apeglm \
-    --report-genes FAM138A \
+    --report-genes "<GENE_SYMBOL>" \
     --workdir /tmp/deseq2_run
 ```
 
 Output highlights (parseable):
 ```
-# CONTRAST sex_M_vs_F: n=37496 n_tested=26591
+# CONTRAST <contrast>: n=<genes-in> n_tested=<genes-after-filtering>
 # SIG_sex_M_vs_F_unshrunk_strict (padj<0.05 AND |LFC|>0.5 AND baseMean>10): n=...
 # SIG_sex_M_vs_F_shrunk_padjlfc (padj<0.05 AND |LFC|>0.5, NO baseMean): n=...
-# GENE FAM138A [sex_M_vs_F]: baseMean=... unshrunkLFC=... shrunkLFC=... padj=...
+# GENE <GENE_SYMBOL> [<contrast>]: baseMean=... unshrunkLFC=... shrunkLFC=... padj=...
 ```
 
 For a multi-strain Venn run with notebook-style outlier exclusion:
@@ -60,9 +60,9 @@ For a multi-strain Venn run with notebook-style outlier exclusion:
 python scripts/r_deseq2_wrapper.py \
     --counts .../raw_counts.csv \
     --metadata .../experiment_metadata.csv \
-    --design "~Replicate + Strain + Media" \
-    --multi-contrast "Strain,97,1;Strain,98,1;Strain,99,1" \
-    --exclude-samples "resub-5,resub-10,resub-33" \
+    --design "~replicate + group + media" \
+    --multi-contrast "group,g1,ctrl;group,g2,ctrl;group,g3,ctrl" \
+    --exclude-samples "sample-7,sample-12" \
     --lfc-thr 1.5 --padj-thr 0.05 --basemean-thr 0 \
     --workdir /tmp/strain_venn
 ```
@@ -130,11 +130,11 @@ this AFTER `r_deseq2_wrapper.py` if you need to explore the
 
 ```bash
 python scripts/multi_strain_venn.py \
-    --deg-csv "JBX97=/tmp/strain_venn/res_unshrunk_Strain_97_vs_1.csv" \
-    --deg-csv "JBX98=/tmp/strain_venn/res_unshrunk_Strain_98_vs_1.csv" \
-    --deg-csv "JBX99=/tmp/strain_venn/res_unshrunk_Strain_99_vs_1.csv" \
+    --deg-csv "g1=/tmp/venn_run/res_unshrunk_group_g1_vs_ctrl.csv" \
+    --deg-csv "g2=/tmp/venn_run/res_unshrunk_group_g2_vs_ctrl.csv" \
+    --deg-csv "g3=/tmp/venn_run/res_unshrunk_group_g3_vs_ctrl.csv" \
     --padj-thr 0.05 --lfc-thr 1.5 \
-    --target-set "JBX97,JBX99"
+    --target-set "g1,g3"
 ```
 
 Output emits `# PCT |target∩ - others| / |...|` lines for four
@@ -259,7 +259,7 @@ DESeq2 assumes that most genes are NOT differentially expressed — this is its
 - **LFC shrinkage and individual-gene queries**: When the analysis pipeline uses LFC shrinkage but the question asks for a SPECIFIC gene's log2 fold change (especially a low-baseMean gene like a lncRNA with baseMean<10), the natural answer is the UNSHRUNKEN LFC from the standard DESeq2 results table. Shrinkage is designed to pull noisy low-baseMean estimates toward zero — reporting the shrunken value of a low-baseMean gene gives ≈0 which doesn't represent the gene's actual differential expression. Report unshrunken (raw `results()` output) for individual-gene queries; report shrunken values only when the question is about visualization, ranking, or aggregate comparisons.
 - **Venn-diagram percentages — union denominator**: When a question asks "what percentage of genes are DE in X" or "% of genes uniquely/jointly DE in conditions A/B/C", and the analysis is multi-condition with a Venn diagram, the natural denominator is the UNION of all DE sets across the conditions, NOT the total tested gene count. Published Venn diagrams report percentages as `|set| / |union|`. Apply this whenever the question references multiple-condition DE comparisons or Venn-style overlaps — `|union|` is materially smaller than total-tested and gives a different number.
 - **Gene-length vs expression correlation — match the notebook's transform**: When the question asks for the Pearson correlation between gene length and gene expression, the answer depends critically on whether expression is log-transformed first. Raw expression spans ~5 orders of magnitude and the raw correlation can differ substantially from the log-transformed correlation (e.g., raw ≈ near-zero, log ≈ moderate-positive). Pattern: a per-cell-type correlation question (e.g., "among CD8 cells") typically expects the RAW Pearson r (small values like 0.02–0.08, matching the notebook table). A pooled-across-samples question that just says "protein-coding genes" without restricting to one cell type typically expects the LOG10–LOG10 Pearson r (moderate values like 0.30–0.40). Read the executed notebook first; if it only shows per-cell-type raw r values, the pooled "protein-coding only" answer is almost always log10(length) vs log10(mean expression), NOT raw. Run `gene_length_correlation.py`, then pick `raw_pearson_r` for the cell-type subset (cell-type-restricted question) or `log10_both_pearson_r` for ALL_SAMPLES (pooled "protein-coding only" question).
-- **"NOT in any other strain/condition" — enumerate ALL strains/conditions in metadata**: When a question asks "% of genes DE in A AND DE in B AND **NOT in any other strain**", "other" refers to the FULL set of strains in the metadata file — not just the obvious counterpart. List all unique values in the relevant metadata column (e.g., `Strain`, `Genotype`, `Condition`) and exclude the gene if it is DE in ANY of them outside the specified set. Restricting "other strain" to only the most-recently-mentioned counterpart inflates the count by including genes that are co-DE in the unmentioned strains. Common error: question mentions strains JBX97/JBX98/JBX99 (relative to JBX1) and the agent excludes only JBX98, missing additional strains in the design (e.g., a fourth strain or media condition with its own DE set).
+- **"NOT in any other strain/condition" — enumerate ALL strains/conditions in metadata**: When a question asks "% of genes DE in A AND DE in B AND **NOT in any other strain**", "other" refers to the FULL set of strains in the metadata file — not just the obvious counterpart. List all unique values in the relevant metadata column (e.g., `Strain`, `Genotype`, `Condition`) and exclude the gene if it is DE in ANY of them outside the specified set. Restricting "other strain" to only the most-recently-mentioned counterpart inflates the count by including genes that are co-DE in the unmentioned strains. Common error: a question names two or three groups and the agent treats "other" as only the one named alongside them, missing the remaining levels in the design (a further group, or a second factor such as media with its own DE set).
 - **R DESeq2's `apeglm` shrinkage rounds the same as the notebook's pydeseq2**: For most CD4/CD8-style sex-DE questions, both libraries give the same DEG count to within ±2 genes when you match the filter EXACTLY. If the agent's count differs by >5%, the most likely cause is the `baseMean>10` filter being silently applied when the published notebook had it commented out. Use `r_deseq2_wrapper.py` and read both `_strict` (with baseMean) and `_padjlfc` (without baseMean) — the published notebook's `len(sigs)` typically matches `_padjlfc` (the filter without baseMean), even when the question text mentions a baseMean threshold.
 - **PCA orientation pitfalls**: "samples-as-rows" usually means transposing a CSV that loaded with rows=genes. But some published outputs were computed without the transpose (rows=genes), giving very different PC1 percentages. If your strict samples-as-rows answer disagrees with the published value, ALSO run `pca_variance.py` with both orientations and look for a `cum_PC1` value that matches.
 - **Report units that match the question literally**: When the question asks "What percentage of...", report the value as a percentage (e.g., `10.6%`), not as a count or raw fraction. When it asks "How many...", report a count. When it asks "What is the p-value...", report the p-value, not the count of significant items. The grader treats unit/type mismatch as wrong even if the underlying computation was correct. Common errors: returning the count `91` for a percentage question (should be `91/N×100`), returning `0.05` (the threshold) when asked for a count, returning a fraction `0.05` when asked for a percentage (should be `5%`).
@@ -307,7 +307,7 @@ Load counts + metadata, ensure samples-as-rows/genes-as-columns, verify integer
 
 ### Step 2.5: Inspect Metadata (REQUIRED)
 
-List ALL metadata columns and levels. Categorize as biological interest vs batch/block. Build design formula with covariates first, factor of interest last. See [design_formula_guide.md](references/design_formula_guide.md).
+List ALL metadata columns and levels. Categorize as biological interest vs batch/block. Build design formula with covariates first, factor of interest last.
 
 ### Step 3: Run PyDESeq2
 
@@ -391,13 +391,12 @@ See [troubleshooting.md](references/troubleshooting.md) for full debugging guide
 ## Known Limitations
 
 - **PyDESeq2 vs R DESeq2**: Numerical differences exist for very low dispersion genes (<1e-05). For exact R reproducibility, use rpy2.
-- **gseapy vs R clusterProfiler**: Results may differ. See [r_clusterprofiler_guide.md](references/r_clusterprofiler_guide.md).
+- **gseapy vs R clusterProfiler**: Results may differ.
 
 ## Reference Files
 
 - [question_parsing.md](references/question_parsing.md) - Extract parameters from questions
 - [data_loading.md](references/data_loading.md) - Data loading and validation
-- [design_formula_guide.md](references/design_formula_guide.md) - Multi-factor design decision tree
 - [pydeseq2_workflow.md](references/pydeseq2_workflow.md) - Complete PyDESeq2 code examples
 - [result_filtering.md](references/result_filtering.md) - Advanced filtering and extraction
 - [dispersion_analysis.md](references/dispersion_analysis.md) - Dispersion diagnostics
@@ -405,7 +404,6 @@ See [troubleshooting.md](references/troubleshooting.md) for full debugging guide
 - [output_formatting.md](references/output_formatting.md) - Format answers correctly
 - [worked_examples.md](references/worked_examples.md) - All 10 question patterns
 - [troubleshooting.md](references/troubleshooting.md) - Common issues and debugging
-- [r_clusterprofiler_guide.md](references/r_clusterprofiler_guide.md) - R clusterProfiler via rpy2
 - [edger_limma_voom.md](references/edger_limma_voom.md) - edgeR / limma-voom DE routes: command sequences, contracts, column crosswalk
 
 ## Utility Scripts

upstream/mims-harvard-ToolUniverse/skills/tooluniverse-image-analysis/references/segmentation.md

@@ -2,9 +2,9 @@
 title: "Colony and Object Segmentation"
 task: ""
 lineage_type: import
-upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-image-analysis/references/segmentation.md
-upstream_sha: e2520a96
-imported_at: 2026-06-26
+upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/aafb6fb6/skills/tooluniverse-image-analysis/references/segmentation.md
+upstream_sha: aafb6fb6
+imported_at: 2026-11-03
 prompt_class: prompt
 upstream_changes: accepted
 author: upstream
@@ -485,8 +485,8 @@ def colony_morphometry_analysis(df, genotype_col='Genotype',
 
 
 # Example usage (colony morphometry)
-df = pd.read_csv("Swarm_1.csv")
-result = colony_morphometry_analysis(df, 'Genotype', 'Area', 'Circularity')
+df = pd.read_csv("colony_measurements.csv")
+result = colony_morphometry_analysis(df, 'genotype', 'area', 'circularity')
 print(f"Genotype with largest area: {result['max_area_genotype']}")
 print(f"Mean circularity: {result['max_area_circularity']:.4f}")
 ```

upstream/mims-harvard-ToolUniverse/skills/tooluniverse-image-analysis/references/statistical_analysis.md

@@ -2,9 +2,9 @@
 title: "Statistical Analysis for Imaging Data"
 task: ""
 lineage_type: import
-upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-image-analysis/references/statistical_analysis.md
-upstream_sha: e2520a96
-imported_at: 2026-06-26
+upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/aafb6fb6/skills/tooluniverse-image-analysis/references/statistical_analysis.md
+upstream_sha: aafb6fb6
+imported_at: 2026-11-03
 prompt_class: prompt
 upstream_changes: accepted
 author: upstream
@@ -155,7 +155,7 @@ def mann_whitney_test(group1, group2):
 
 ### Dunnett's Test
 
-**CRITICAL**: Dunnett's test compares each treatment group to a single control group. R uses `multcomp::glht()` with `mcp(Strain_Ratio = "Dunnett")`. Python equivalent uses scipy.
+**CRITICAL**: Dunnett's test compares each treatment group to a single control group. R uses `multcomp::glht()` with `mcp(<group_col> = "Dunnett")`. Python equivalent uses scipy.
 
 ```python
 from scipy import stats
@@ -294,11 +294,11 @@ def extract_anova_interaction(anova_table, factor1, factor2):
     return f_stat, p_val
 
 # Example usage
-anova_result = two_way_anova(df, 'NeuN_count', 'Condition', 'Sex')
+anova_result = two_way_anova(df, 'cell_count', 'treatment', 'batch')
 print(anova_result)
 
 # Extract interaction
-f_stat, p_val = extract_anova_interaction(anova_result, 'Condition', 'Sex')
+f_stat, p_val = extract_anova_interaction(anova_result, 'treatment', 'batch')
 print(f"Interaction F={f_stat:.3f}, p={p_val:.4f}")
 ```
 
@@ -404,36 +404,38 @@ def calculate_achieved_power(n_per_group, effect_size, alpha=0.05, alternative='
 
 ## Regression Modeling
 
-### Data Preparation for Co-culture Ratios
+### Turning Mixture Ratios Into a Continuous Predictor
+
+A mixture experiment often records composition as a ratio string such as
+`"3:1"`. Regression needs a number, so convert the ratio to the fraction of one
+component. Pass the column names your own data uses; nothing here assumes a
+particular dataset.
 
 ```python
-def prepare_ratio_data(df, strain_col='StrainNumber', ratio_col='Ratio',
-                       area_col='Area', exclude_strains=None):
-    """Prepare co-culture ratio data for regression analysis.
-
-    Converts ratio strings (e.g., "3:1") to frequency fractions.
-    Filters out pure strains if requested.
+def ratio_to_frequency(df, ratio_col, group_col=None, exclude_groups=None):
+    """Convert "a:b" composition strings to the fraction of the first component.
 
     Args:
-        df: DataFrame with swarming data
-        strain_col: Column with strain identifiers
-        ratio_col: Column with ratio strings
-        area_col: Column with area measurements
-        exclude_strains: List of strain IDs to exclude (e.g., pure strains)
+        df: DataFrame with one row per sample
+        ratio_col: Column holding ratio strings like "3:1"
+        group_col: Optional column identifying the sample group
+        exclude_groups: Optional group values to drop before fitting -- for
+            example single-component controls, which have no meaningful
+            mixture fraction
 
-    Returns: DataFrame with Frequency_rhlI column added
+    Returns: DataFrame with `part_a`, `part_b` and `frequency_a` added
     """
     result = df.copy()
 
-    if exclude_strains:
-        result = result[~result[strain_col].isin(exclude_strains)]
+    if exclude_groups:
+        if group_col is None:
+            raise ValueError("exclude_groups needs group_col")
+        result = result[~result[group_col].isin(exclude_groups)]
 
-    # Parse ratio into frequency
-    # Ratio format: "rhlI_D:lasI_D" (e.g., "3:1" means 3 parts rhlI, 1 part lasI)
-    ratio_parts = result[ratio_col].str.split(':', expand=True).astype(int)
-    result['rhlI_D'] = ratio_parts[0]
-    result['lasI_D'] = ratio_parts[1]
-    result['Frequency_rhlI'] = result['rhlI_D'] / (result['rhlI_D'] + result['lasI_D'])
+    parts = result[ratio_col].str.split(':', expand=True).astype(int)
+    result['part_a'] = parts[0]
+    result['part_b'] = parts[1]
+    result['frequency_a'] = result['part_a'] / (result['part_a'] + result['part_b'])
 
     return result
 ```
@@ -525,7 +527,9 @@ def fit_polynomial_model(df, x_col, y_col, degree=2):
 
 ### Natural Spline Regression
 
-**CRITICAL**: This must match R's `lm(Area ~ ns(Frequency_rhlI, df=4))`.
+**CRITICAL**: `patsy`'s default spline basis does not reproduce R's
+`ns()`. To match `lm(y ~ ns(x, df=k))`, use explicit quantile knots as
+below; the fitted values differ otherwise.
 
 ```python
 from patsy import dmatrix
@@ -650,37 +654,41 @@ def compare_regression_models(models_dict):
 
 ---
 
-## Answer Extraction Patterns
+## Reporting Numeric Results
 
-### Answer formatting
+### Rounding conventions
+
+Round once, at the point of reporting, and say which convention you used. The
+helper below is a place to keep that decision rather than scattering `round()`
+calls through an analysis.
 
 ```python
-def format_answer(value, question_type):
-    """Format answer for reporting.
+def round_for_report(value, convention):
+    """Round a statistic for reporting under a named convention.
 
     Args:
         value: Numeric value to format
-        question_type: Type of answer expected
+        convention: Which rounding convention to apply
 
     Returns: Formatted value
     """
-    if question_type == "nearest_thousand":
+    if convention == "nearest_thousand":
         return int(round(value, -3))
-    elif question_type == "percentage_int":
+    elif convention == "percentage_int":
         return int(round(value))
-    elif question_type == "percentage_2dec":
+    elif convention == "percentage_2dec":
         return round(value, 2)
-    elif question_type == "cohen_d":
+    elif convention == "cohen_d":
         return round(value, 3)
-    elif question_type == "statistic_3dec":
+    elif convention == "statistic_3dec":
         return round(value, 3)
-    elif question_type == "sample_size":
+    elif convention == "sample_size":
         return int(np.ceil(value))
-    elif question_type == "count":
+    elif convention == "count":
         return int(value)
-    elif question_type == "r_squared":
+    elif convention == "r_squared":
         return round(value, 2)
-    elif question_type == "p_value":
+    elif convention == "p_value":
         if value < 0.0001:
             return f"{value:.2e}"
         else:
@@ -691,42 +699,50 @@ def format_answer(value, question_type):
 
 ---
 
-## Complete Example: Spline Regression Workflow
+## Worked Example: Choosing Between Curve Shapes
+
+A measurement that rises and then falls across a continuous predictor can be
+fitted several ways. Fit the candidates, compare them on the same data, and
+report the peak from whichever wins -- rather than assuming a shape.
+
+Substitute your own file and column names; the variables below are placeholders.
 
 ```python
 import pandas as pd
 import numpy as np
 
-# Load data
-df = pd.read_csv("Swarm_2.csv")
+df = pd.read_csv("measurements.csv")
 
-# Prepare ratio data
-df_coculture = prepare_ratio_data(
+# If composition is recorded as a ratio string, convert it to a fraction first.
+# Drop single-component controls, which have no mixture fraction to place them
+# on the x axis.
+single_component_controls = ['control-a', 'control-b']
+
+fitted = ratio_to_frequency(
     df,
-    strain_col='StrainNumber',
-    ratio_col='Ratio',
-    area_col='Area',
-    exclude_strains=['1', '98']  # Exclude pure strains
+    ratio_col='composition',
+    group_col='sample_group',
+    exclude_groups=single_component_controls,
 )
 
-# Fit models
-quadratic = fit_polynomial_model(df_coculture, 'Frequency_rhlI', 'Area', degree=2)
-cubic = fit_polynomial_model(df_coculture, 'Frequency_rhlI', 'Area', degree=3)
-spline = fit_natural_spline_model(df_coculture, 'Frequency_rhlI', 'Area', spline_df=4)
+x, y = 'frequency_a', 'response'
+
+models = {
+    'quadratic': fit_polynomial_model(fitted, x, y, degree=2),
+    'cubic': fit_polynomial_model(fitted, x, y, degree=3),
+    'spline': fit_natural_spline_model(fitted, x, y, spline_df=4),
+}
 
-# Compare models
-models = {'quadratic': quadratic, 'cubic': cubic, 'spline': spline}
 comparison = compare_regression_models(models)
 print(comparison)
 
-# Best model by R-squared
-best_model_name = comparison.iloc[0]['model']
-best_model = models[best_model_name]
+best_name = comparison.iloc[0]['model']
+best = models[best_name]
 
-# Extract answers
-print(f"Best model: {best_model_name}")
-print(f"R-squared: {best_model['r_squared']:.4f}")
-print(f"Peak frequency: {best_model['peak_x']:.4f}")
-print(f"Peak area: {best_model['peak_value']:.1f}")
-print(f"95% CI: [{best_model['peak_ci_lower']:.1f}, {best_model['peak_ci_upper']:.1f}]")
+print(f"Best model: {best_name}  (R-squared {best['r_squared']:.4f})")
+print(f"Peak at x={best['peak_x']:.4f}, y={best['peak_value']:.1f}")
+print(f"95% CI for the peak: [{best['peak_ci_lower']:.1f}, {best['peak_ci_upper']:.1f}]")
 ```
+
+Report which model you chose and why. A peak read off a shape the data did not
+support is a different number from a peak read off the best-fitting one.