From 7d1343e153d792d354cf3a9340d1bf1ebae94c7b Mon Sep 17 00:00:00 2001 From: promptadmin Date: Tue, 3 Nov 2026 17:02:06 +0000 Subject: [PATCH 1/3] [upstream-sync] skills/tooluniverse-image-analysis/references/segmentation.md from mims-harvard/ToolUniverse@aafb6fb6 [prompt] --- .../references/segmentation.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/upstream/mims-harvard-ToolUniverse/skills/tooluniverse-image-analysis/references/segmentation.md b/upstream/mims-harvard-ToolUniverse/skills/tooluniverse-image-analysis/references/segmentation.md index 488f059d..e3994735 100644 --- a/upstream/mims-harvard-ToolUniverse/skills/tooluniverse-image-analysis/references/segmentation.md +++ b/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}") ``` -- 2.54.0 From 7cab31de90c016d92f52e54fce464127f58d915c Mon Sep 17 00:00:00 2001 From: promptadmin Date: Tue, 3 Nov 2026 17:02:09 +0000 Subject: [PATCH 2/3] [upstream-sync] skills/tooluniverse-image-analysis/references/statistical_analysis.md from mims-harvard/ToolUniverse@aafb6fb6 [prompt] --- .../references/statistical_analysis.md | 148 ++++++++++-------- 1 file changed, 82 insertions(+), 66 deletions(-) diff --git a/upstream/mims-harvard-ToolUniverse/skills/tooluniverse-image-analysis/references/statistical_analysis.md b/upstream/mims-harvard-ToolUniverse/skills/tooluniverse-image-analysis/references/statistical_analysis.md index 8061e0c8..39cac5b8 100644 --- a/upstream/mims-harvard-ToolUniverse/skills/tooluniverse-image-analysis/references/statistical_analysis.md +++ b/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( = "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. -- 2.54.0 From e55512404e65d6ff666dec394ed80bd5e715af30 Mon Sep 17 00:00:00 2001 From: promptadmin Date: Tue, 3 Nov 2026 17:02:12 +0000 Subject: [PATCH 3/3] [upstream-sync] skills/tooluniverse-rnaseq-deseq2/SKILL.md from mims-harvard/ToolUniverse@aafb6fb6 [catalogue] --- .../tooluniverse-rnaseq-deseq2/SKILL.md | 38 +++++++++---------- 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/upstream/mims-harvard-ToolUniverse/catalogue/skills/tooluniverse-rnaseq-deseq2/SKILL.md b/upstream/mims-harvard-ToolUniverse/catalogue/skills/tooluniverse-rnaseq-deseq2/SKILL.md index 3c7ca7ce..446102b8 100644 --- a/upstream/mims-harvard-ToolUniverse/catalogue/skills/tooluniverse-rnaseq-deseq2/SKILL.md +++ b/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 /counts.csv \ --metadata /meta.csv \ --design "~sex" --contrast "sex,M,F" \ - --subset-col celltype --subset-values "CD4,CD8" \ + --subset-col celltype --subset-values "," \ --min-row-sum 10 --shrink apeglm \ - --report-genes FAM138A \ + --report-genes "" \ --workdir /tmp/deseq2_run ``` Output highlights (parseable): ``` -# CONTRAST sex_M_vs_F: n=37496 n_tested=26591 +# CONTRAST : n= n_tested= # 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 []: 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 -- 2.54.0