tooluniverse-epigenomics
Genomics and epigenomics analysis: DNA methylation (CpG, 5mC, 5hmC, bisulfite, RRBS), m6A RNA modification (MeRIP-seq), ChIP-seq peaks, ATAC-seq accessibility, histone modifications, chromatin state, multi-omics integration. Combines pandas/scipy/pysam computation with ToolUniverse annotation tools.
By mims-harvard · 407 installs
npx skills add mims-harvard/tooluniverse --skill tooluniverse-epigenomics
Source repository · Upstream listing
Genomics and Epigenomics Data Processing
⚠️ TOP OF MIND RULE: long format methylation CSV — count ROWS, not unique positions
When the input is a long format methylation CSV (one row per (sample, CpG position)
e.g. columns Pos, Chromosome, MethylationPercentage ), "how many sites are
removed when filtering" almost always means rows removed , NOT unique position
removals. The two answers differ by a factor of ≈ n samples .
Question phrasing What it means
"how many sites are removed when filtering …" rows removed (= samples × positions failing the filter)
"how many unique CpG sites pass filter" unique positions (dedupe by Pos then filter)
❌ WRONG: df.drop duplicates(["Pos"]).query("MethylationPercentage<10 or 90") then len(filtered) → counts unique positions (typically 100–1500)
✅ RIGHT: df.query("MethylationPercentage<10 or MethylationPercentage 90") then len(df) len(filtered) → counts rows (typically 10k–30k)
If your answer is < 2000 when the data has 1000+ positions × 20+ samples, you
deduplicated too early. Re read the question's noun before reporting.
RULE ZERO — Check for pre computed results FIRST
Before following any instruction below, scan the data folder for:
executed.ipynb → read with tu run read executed notebook '{"data folder":"<path ","search":"<keyword "}' and cite its cell outputs as the authoritative answer
Pre computed result files (CSV/TSV with names like results , deseq , enrich , stats , simplified.csv ) → read directly and report the requested value
Canonical analysis scripts ( analysis.R , run .py , find .R , .Rmd ) → execute as is and read the output
Only follow this skill's re analysis recipe below if none of the above exist. Re running from raw data produces different numbers than the published answer and is much slower (often 5 10× turn count).
Production ready skill combining Python computation (pandas, scipy, numpy, pysam, statsmodels) with ToolUniverse annotation tools for epigenomics analysis.
LOOK UP, DON'T GUESS
When uncertain about any scientific fact, SEARCH databases first.
When to Use
Methylation data, ChIP seq peaks, ATAC seq, multi omics integration, genome wide epigenomic statistics. Keywords: methylation, CpG, ChIP seq, ATAC seq, histone, chromatin, epigenetic.
NOT for : RNA seq DEG, variant calling, gene enrichment, protein structure.
Key Principles
1. Data first Load/inspect before analysis
2. Question driven Extract specific numeric answer
3. Coordinate system awareness Track genome build (hg19/hg38/mm10), chr prefix
4. Statistical rigor FDR correction, effect size filtering
5. CpG identification Parse Illumina probe IDs, genomic coordinates
PRIMARY SCRIPT — methylation density.py (use FIRST for CpG density questions)
For long format methylation CSVs ( Pos, Chromosome, MethylationPercentage )
paired with chromosome length CSVs, ALWAYS run the bundled script before
hand rolling pandas. It deterministically computes every common metric in one
pass and avoids the rows vs sites pitfall that produces silently wrong answers.
The full JSON output contains every metric. Pick the one that matches the
question's wording (NOT a similar looking one):
Question phrasing Script field
"how many sites are removed when filtering …" rows removed
"how many unique CpG sites pass filter" unique pos after filter
"genome wide AVERAGE chromosomal density" density avg per chr
"density on chromosome X" density chromosome (pass chromosome X )
"total density across the genome" density total over genome
The two density numbers ( density avg per chr vs density total over genome )
typically differ by ~2× because CpGs are not uniformly distributed across
chromosomes; reporting one when the question asks for the other is the most
common failure mode here.
For "sites removed" questions, the long format CSV has multiple rows per CpG
position (one per sample), so rows removed is in the tens of thousands while
unique pos removed is in the hundreds. Match the granularity to the question.
Distinguish "rows" vs "unique sites" — methylation CSVs are usually long format
CpG methylation CSVs typically have ONE ROW PER (sample × CpG site) — so len(df) n unique sites . Before computing anything, decide which axis the question is asking about:
Question phrasing Axis Operation
"how many sites are removed when filtering" sample rows filter then count rows; do NOT dedupe by Pos . The CSV is in long format; "sites" here is row shaped. Subtract len(df filtered) from len(df) .
"how many unique CpG sites pass filter" unique positions dedupe by position (or Pos column), then filter
"genome wide average chromosomal density" per chromosome density MEAN of per chromosome densities: (n unique per chr / chr length).mean() . NOT total unique / total genome — that gives a different answer (typically ≈ ½ of the per chr mean for unevenly distributed CpGs).
"density on chromosome X" single chromosome unique positions on X / length(X). Be careful which species — check the question text for "Zebra Finch" vs "Jackdaw".
"chi square for uniform distribution across chromosomes" unique positions per chromosome filter rows first, then dedupe by (Chromosome, Pos) , then count per chromosome unique positions for chi square against expected = chr length / total length × n unique filtered
Sanity check : if your filtered count is two orders of magnitude smaller than the GT range, you likely deduped when the question wanted row level counts (or vice versa). Re run with the other axis and compare.
For the chi square uniformity test: expected counts = chromosome length / total genome length × n unique sites . The chi square statistic depends on the count granularity (rows vs unique sites) — a row level chi square gives a much higher chi square than a unique position chi square because n is larger.
Precedence : when an executed.ipynb exists, read its filtering code verbatim — df[(df.MethylationPercentage 90) (df.MethylationPercentage < 10)] (no dedup) and df.drop duplicates('Pos') (with dedup) yield wildly different counts on the same dataset.
Workflow
Phase 0: Question Parsing
Identify data files, specific statistic, thresholds, genome build. Categorize by keywords.
See ANALYSIS PROCEDURES.md for decision tree.
Phase 1: Methylation Processing
Load beta/M value matrix (CSV/TSV/parquet/HDF5)
Filter by variance, missing rate, probe type, chromosome, CpG island relation
Differential methylation: T test/Wilcoxon between groups + FDR
Age related CpG: Pearson/Spearman correlation + FDR
Chromosome density: CpG count / chromosome length
Phase 2: ChIP seq Peak Analysis
Load BED/narrowPeak/broadPeak, normalize chromosomes
Peak stats, annotation to genes, overlap analysis (Jaccard)
Phase 3: ATAC seq
NFR detection (<150bp peaks), region classification
Phase 4: Multi Omics Integration
Methylation expression correlation per probe gene (Pearson/Spearman + FDR)
ChIP seq + expression: promoter peaks vs expression levels
Phase 5: Clinical Data
Missing data analysis across modalities, complete case identification
Phase 6: ToolUniverse Annotation
ENCODE tools :
ENCODE search rnaseq experiments : assay type ("total RNA seq" default; fall back to "polyA plus RNA seq"), biosample , limit
ENCODE search histone experiments : target (e.g., "H3K27ac"), cell type / tissue / biosample , limit
GEO tools : GEO search rnaseq datasets , GEO search atacseq datasets both accept limit or max results
GTEx tools :
GTEx get median gene expression : gene symbol (NOT Ensembl ID)
GTEx query eqtl : gene symbol , tissue id (case sensitive exact, e.g., "Whole Blood" )
Other : ensembl lookup gene (requires species='homo sapiens' ), ensembl get regulatory features (NO "chr" prefix), SCREEN get regulatory elements , ChIPAtlas (requires operation param), SRA search experiments (library strategy: "ChIP Seq"/"Bisulfite Seq"/"ATAC seq")
Phase 7: Genome Wide Statistics
Global mean/median beta, probe variance, chromosome density, DMP counts.
See CODE REFERENCE.md for full implementations.
Common Patterns
Pattern Key Steps
Differential methylation Filter probes → groups → t test → FDR → threshold
Age related CpG density Correlate with age → FDR → map to chr → density ratio
Multi omics missing data Extract IDs → intersect → check NaN → complete case count
ChIP seq annotation Load peaks → annotate genes → classify regions
Methylation expression Align samples → correlate → FDR → anti correlations
GTEx Tissue IDs
Whole Blood, Liver, Lung, Breast Mammary Tissue, Brain Cortex, Heart Left Ventricle, Kidney Cortex, Thyroid, Adipose Subcutaneous, Muscle Skeletal
Evidence Grading
Grade Criteria
Strong padj < 0.01 AND abs(delta beta) = 0.2, replicated
Moderate padj < 0.05 AND abs(delta beta) = 0.1
Weak padj < 0.05 but delta beta < 0.1
Insufficient padj = 0.05 or no replication
Delta beta = 0.2 = strong effect. ChIP seq: q < 0.01, FE = 2 for confidence. ATAC seq NFR < 150bp = active regulatory. Always apply BH FDR. Verify genome build consistency.
Limitations
No pybedtools/pyBigWig: pure Python intervals
Illumina centric (450K/EPIC); uses t test/Wilcoxon (not limma)
No peak calling (assumes pre called)
API rate limits: ~20 genes per batch
Reference Files
CODE REFERENCE.md , TOOLS REFERENCE.md , ANALYSIS PROCEDURES.md , QUICK START.md