promptadmin/drug-discovery-prompts
1
0
Fork 0
Code Issues Pull Requests 3 Packages Projects Releases Wiki Activity

[Upstream sync] mims-harvard/ToolUniverse (github) — 0 added, 3 modified #66

Open
promptadmin wants to merge 3 commits from upstream-sync/tooluniverse-20261103-aafb6f-tgab into main
pull from: upstream-sync/tooluniverse-20261103-aafb6f-tgab
merge into: :main
Conversation 0 Commits 3 Files Changed 3
+105 -91
3 changed files with 105 additions and 91 deletions
Download Patch File Download Diff File Expand all files Collapse all files
upstream/mims-harvard-ToolUniverse/catalogue/skills/tooluniverse-rnaseq-deseq2/SKILL.md
+18 -20
View File
@@ -1,8 +1,8 @@
--- ---
lineage_type: import lineage_type: import
upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-rnaseq-deseq2/SKILL.md upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/aafb6fb6/skills/tooluniverse-rnaseq-deseq2/SKILL.md
upstream_sha: e2520a96 upstream_sha: aafb6fb6
imported_at: 2026-06-26 imported_at: 2026-11-03
prompt_class: catalogue prompt_class: catalogue
upstream_changes: accepted upstream_changes: accepted
name: tooluniverse-rnaseq-deseq2 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. 3-way Venn region sizes and percentage-of-X interpretations.
```bash ```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 \ python scripts/r_deseq2_wrapper.py \
--counts <data-folder>/counts.csv \ --counts <data-folder>/counts.csv \
--metadata <data-folder>/meta.csv \ --metadata <data-folder>/meta.csv \
--design "~sex" --contrast "sex,M,F" \ --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 \ --min-row-sum 10 --shrink apeglm \
--report-genes FAM138A \ --report-genes "<GENE_SYMBOL>" \
--workdir /tmp/deseq2_run --workdir /tmp/deseq2_run
``` ```
Output highlights (parseable): 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_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=... # 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: 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 \ python scripts/r_deseq2_wrapper.py \
--counts .../raw_counts.csv \ --counts .../raw_counts.csv \
--metadata .../experiment_metadata.csv \ --metadata .../experiment_metadata.csv \
--design "~Replicate + Strain + Media" \ --design "~replicate + group + media" \
--multi-contrast "Strain,97,1;Strain,98,1;Strain,99,1" \ --multi-contrast "group,g1,ctrl;group,g2,ctrl;group,g3,ctrl" \
--exclude-samples "resub-5,resub-10,resub-33" \ --exclude-samples "sample-7,sample-12" \
--lfc-thr 1.5 --padj-thr 0.05 --basemean-thr 0 \ --lfc-thr 1.5 --padj-thr 0.05 --basemean-thr 0 \
--workdir /tmp/strain_venn --workdir /tmp/strain_venn
``` ```
@@ -130,11 +130,11 @@ this AFTER `r_deseq2_wrapper.py` if you need to explore the
```bash ```bash
python scripts/multi_strain_venn.py \ python scripts/multi_strain_venn.py \
--deg-csv "JBX97=/tmp/strain_venn/res_unshrunk_Strain_97_vs_1.csv" \ --deg-csv "g1=/tmp/venn_run/res_unshrunk_group_g1_vs_ctrl.csv" \
--deg-csv "JBX98=/tmp/strain_venn/res_unshrunk_Strain_98_vs_1.csv" \ --deg-csv "g2=/tmp/venn_run/res_unshrunk_group_g2_vs_ctrl.csv" \
--deg-csv "JBX99=/tmp/strain_venn/res_unshrunk_Strain_99_vs_1.csv" \ --deg-csv "g3=/tmp/venn_run/res_unshrunk_group_g3_vs_ctrl.csv" \
--padj-thr 0.05 --lfc-thr 1.5 \ --padj-thr 0.05 --lfc-thr 1.5 \
--target-set "JBX97,JBX99" --target-set "g1,g3"
``` ```
Output emits `# PCT |target∩ - others| / |...|` lines for four 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. - **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. - **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). - **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. - **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. - **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%`). - **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) ### 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 ### Step 3: Run PyDESeq2
@@ -391,13 +391,12 @@ See [troubleshooting.md](references/troubleshooting.md) for full debugging guide
## Known Limitations ## Known Limitations
- **PyDESeq2 vs R DESeq2**: Numerical differences exist for very low dispersion genes (<1e-05). For exact R reproducibility, use rpy2. - **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 ## Reference Files
- [question_parsing.md](references/question_parsing.md) - Extract parameters from questions - [question_parsing.md](references/question_parsing.md) - Extract parameters from questions
- [data_loading.md](references/data_loading.md) - Data loading and validation - [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 - [pydeseq2_workflow.md](references/pydeseq2_workflow.md) - Complete PyDESeq2 code examples
- [result_filtering.md](references/result_filtering.md) - Advanced filtering and extraction - [result_filtering.md](references/result_filtering.md) - Advanced filtering and extraction
- [dispersion_analysis.md](references/dispersion_analysis.md) - Dispersion diagnostics - [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 - [output_formatting.md](references/output_formatting.md) - Format answers correctly
- [worked_examples.md](references/worked_examples.md) - All 10 question patterns - [worked_examples.md](references/worked_examples.md) - All 10 question patterns
- [troubleshooting.md](references/troubleshooting.md) - Common issues and debugging - [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 - [edger_limma_voom.md](references/edger_limma_voom.md) - edgeR / limma-voom DE routes: command sequences, contracts, column crosswalk
## Utility Scripts ## Utility Scripts
upstream/mims-harvard-ToolUniverse/skills/tooluniverse-image-analysis/references/segmentation.md
+5 -5
View File
@@ -2,9 +2,9 @@
title: "Colony and Object Segmentation" title: "Colony and Object Segmentation"
task: "" task: ""
lineage_type: import lineage_type: import
upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-image-analysis/references/segmentation.md upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/aafb6fb6/skills/tooluniverse-image-analysis/references/segmentation.md
upstream_sha: e2520a96 upstream_sha: aafb6fb6
imported_at: 2026-06-26 imported_at: 2026-11-03
prompt_class: prompt prompt_class: prompt
upstream_changes: accepted upstream_changes: accepted
author: upstream author: upstream
@@ -485,8 +485,8 @@ def colony_morphometry_analysis(df, genotype_col='Genotype',
# Example usage (colony morphometry) # Example usage (colony morphometry)
df = pd.read_csv("Swarm_1.csv") df = pd.read_csv("colony_measurements.csv")
result = colony_morphometry_analysis(df, 'Genotype', 'Area', 'Circularity') result = colony_morphometry_analysis(df, 'genotype', 'area', 'circularity')
print(f"Genotype with largest area: {result['max_area_genotype']}") print(f"Genotype with largest area: {result['max_area_genotype']}")
print(f"Mean circularity: {result['max_area_circularity']:.4f}") print(f"Mean circularity: {result['max_area_circularity']:.4f}")
``` ```
upstream/mims-harvard-ToolUniverse/skills/tooluniverse-image-analysis/references/statistical_analysis.md
+82 -66
View File
@@ -2,9 +2,9 @@
title: "Statistical Analysis for Imaging Data" title: "Statistical Analysis for Imaging Data"
task: "" task: ""
lineage_type: import lineage_type: import
upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-image-analysis/references/statistical_analysis.md upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/aafb6fb6/skills/tooluniverse-image-analysis/references/statistical_analysis.md
upstream_sha: e2520a96 upstream_sha: aafb6fb6
imported_at: 2026-06-26 imported_at: 2026-11-03
prompt_class: prompt prompt_class: prompt
upstream_changes: accepted upstream_changes: accepted
author: upstream author: upstream
@@ -155,7 +155,7 @@ def mann_whitney_test(group1, group2):
### Dunnett's Test ### 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 ```python
from scipy import stats from scipy import stats
@@ -294,11 +294,11 @@ def extract_anova_interaction(anova_table, factor1, factor2):
return f_stat, p_val return f_stat, p_val
# Example usage # 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) print(anova_result)
# Extract interaction # 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}") 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 ## 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 ```python
def prepare_ratio_data(df, strain_col='StrainNumber', ratio_col='Ratio', def ratio_to_frequency(df, ratio_col, group_col=None, exclude_groups=None):
area_col='Area', exclude_strains=None): """Convert "a:b" composition strings to the fraction of the first component.
"""Prepare co-culture ratio data for regression analysis.
Converts ratio strings (e.g., "3:1") to frequency fractions.
Filters out pure strains if requested.
Args: Args:
df: DataFrame with swarming data df: DataFrame with one row per sample
strain_col: Column with strain identifiers ratio_col: Column holding ratio strings like "3:1"
ratio_col: Column with ratio strings group_col: Optional column identifying the sample group
area_col: Column with area measurements exclude_groups: Optional group values to drop before fitting -- for
exclude_strains: List of strain IDs to exclude (e.g., pure strains) 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() result = df.copy()
if exclude_strains: if exclude_groups:
result = result[~result[strain_col].isin(exclude_strains)] if group_col is None:
raise ValueError("exclude_groups needs group_col")
result = result[~result[group_col].isin(exclude_groups)]
# Parse ratio into frequency parts = result[ratio_col].str.split(':', expand=True).astype(int)
# Ratio format: "rhlI_D:lasI_D" (e.g., "3:1" means 3 parts rhlI, 1 part lasI) result['part_a'] = parts[0]
ratio_parts = result[ratio_col].str.split(':', expand=True).astype(int) result['part_b'] = parts[1]
result['rhlI_D'] = ratio_parts[0] result['frequency_a'] = result['part_a'] / (result['part_a'] + result['part_b'])
result['lasI_D'] = ratio_parts[1]
result['Frequency_rhlI'] = result['rhlI_D'] / (result['rhlI_D'] + result['lasI_D'])
return result return result
``` ```
@@ -525,7 +527,9 @@ def fit_polynomial_model(df, x_col, y_col, degree=2):
### Natural Spline Regression ### 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 ```python
from patsy import dmatrix 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 ```python
def format_answer(value, question_type): def round_for_report(value, convention):
"""Format answer for reporting. """Round a statistic for reporting under a named convention.
Args: Args:
value: Numeric value to format value: Numeric value to format
question_type: Type of answer expected convention: Which rounding convention to apply
Returns: Formatted value Returns: Formatted value
""" """
if question_type == "nearest_thousand": if convention == "nearest_thousand":
return int(round(value, -3)) return int(round(value, -3))
elif question_type == "percentage_int": elif convention == "percentage_int":
return int(round(value)) return int(round(value))
elif question_type == "percentage_2dec": elif convention == "percentage_2dec":
return round(value, 2) return round(value, 2)
elif question_type == "cohen_d": elif convention == "cohen_d":
return round(value, 3) return round(value, 3)
elif question_type == "statistic_3dec": elif convention == "statistic_3dec":
return round(value, 3) return round(value, 3)
elif question_type == "sample_size": elif convention == "sample_size":
return int(np.ceil(value)) return int(np.ceil(value))
elif question_type == "count": elif convention == "count":
return int(value) return int(value)
elif question_type == "r_squared": elif convention == "r_squared":
return round(value, 2) return round(value, 2)
elif question_type == "p_value": elif convention == "p_value":
if value < 0.0001: if value < 0.0001:
return f"{value:.2e}" return f"{value:.2e}"
else: 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 ```python
import pandas as pd import pandas as pd
import numpy as np import numpy as np
# Load data df = pd.read_csv("measurements.csv")
df = pd.read_csv("Swarm_2.csv")
# Prepare ratio data # If composition is recorded as a ratio string, convert it to a fraction first.
df_coculture = prepare_ratio_data( # 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, df,
strain_col='StrainNumber', ratio_col='composition',
ratio_col='Ratio', group_col='sample_group',
area_col='Area', exclude_groups=single_component_controls,
exclude_strains=['1', '98'] # Exclude pure strains
) )
# Fit models x, y = 'frequency_a', 'response'
quadratic = fit_polynomial_model(df_coculture, 'Frequency_rhlI', 'Area', degree=2)
cubic = fit_polynomial_model(df_coculture, 'Frequency_rhlI', 'Area', degree=3) models = {
spline = fit_natural_spline_model(df_coculture, 'Frequency_rhlI', 'Area', spline_df=4) '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) comparison = compare_regression_models(models)
print(comparison) print(comparison)
# Best model by R-squared best_name = comparison.iloc[0]['model']
best_model_name = comparison.iloc[0]['model'] best = models[best_name]
best_model = models[best_model_name]
# Extract answers print(f"Best model: {best_name} (R-squared {best['r_squared']:.4f})")
print(f"Best model: {best_model_name}") print(f"Peak at x={best['peak_x']:.4f}, y={best['peak_value']:.1f}")
print(f"R-squared: {best_model['r_squared']:.4f}") print(f"95% CI for the peak: [{best['peak_ci_lower']:.1f}, {best['peak_ci_upper']:.1f}]")
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}]")
``` ```
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.