| name | tooluniverse-epigenomics |
| description | Production-ready genomics and epigenomics data processing for BixBench questions. Handles methylation array analysis (CpG filtering, differential methylation, age-related CpG detection, chromosome-level density), ChIP-seq peak analysis (peak calling, motif enrichment, coverage stats), ATAC-seq chromatin accessibility, multi-omics integration (expression + methylation correlation), and genome-wide statistics. Pure Python computation (pandas, scipy, numpy, pysam, statsmodels) plus ToolUniverse annotation tools (Ensembl, ENCODE, SCREEN, JASPAR, ReMap, RegulomeDB, ChIPAtlas). Supports BED, BigWig, methylation beta-value matrices, Illumina manifest files, and multi-sample clinical data. Use when processing methylation data, ChIP-seq peaks, ATAC-seq signals, or answering questions about CpG sites, differential methylation, chromatin accessibility, histone marks, or epigenomic statistics. |
Genomics and Epigenomics Data Processing
Production-ready computational skill for processing and analyzing epigenomics data. Combines local Python computation (pandas, scipy, numpy, pysam, statsmodels) with ToolUniverse annotation tools for regulatory context. Designed to solve BixBench-style questions about methylation, ChIP-seq, ATAC-seq, and multi-omics integration.
When to Use This Skill
Triggers:
- User provides methylation data (beta-value matrices, Illumina arrays) and asks about CpG sites
- Questions about differential methylation analysis
- Age-related CpG detection or epigenetic clock questions
- Chromosome-level methylation density or statistics
- ChIP-seq peak files (BED format) with analysis questions
- ATAC-seq chromatin accessibility questions
- Multi-omics integration (expression + methylation, expression + ChIP-seq)
- Genome-wide epigenomic statistics
- Questions mentioning "methylation", "CpG", "ChIP-seq", "ATAC-seq", "histone", "chromatin", "epigenetic"
- Questions about missing data across clinical/genomic/epigenomic modalities
- Regulatory element annotation for processed epigenomic data
Example Questions This Skill Solves:
- "How many patients have no missing data for vital status, gene expression, and methylation data?"
- "What is the ratio of filtered age-related CpG density between chromosomes?"
- "What is the genome-wide average chromosomal density of unique age-related CpGs per base pair?"
- "How many CpG sites show significant differential methylation (padj < 0.05)?"
- "What is the Pearson correlation between methylation and expression for gene X?"
- "How many ChIP-seq peaks overlap with promoter regions?"
- "What fraction of ATAC-seq peaks are in enhancer regions?"
- "Which chromosome has the highest density of hypermethylated CpGs?"
- "Filter CpG sites by variance > threshold and map to nearest genes"
- "What is the average beta value difference between tumor and normal for chromosome 17?"
NOT for (use other skills instead):
- Gene regulation lookup without data files -> Use existing epigenomics annotation pattern
- RNA-seq differential expression -> Use
tooluniverse-rnaseq-deseq2
- Variant calling/annotation from VCF -> Use
tooluniverse-variant-analysis
- Gene enrichment analysis -> Use
tooluniverse-gene-enrichment
- Protein structure analysis -> Use
tooluniverse-protein-structure-retrieval
Required Python Packages
import pandas as pd
import numpy as np
from scipy import stats
import statsmodels.stats.multitest as mt
import pysam
import gseapy
from tooluniverse import ToolUniverse
KEY PRINCIPLES
- Data-first approach - Load and inspect data files BEFORE any analysis
- Question-driven - Parse what the user is actually asking and extract the specific numeric answer
- File format detection - Automatically detect methylation arrays, BED files, BigWig, clinical data
- Coordinate system awareness - Track genome build (hg19, hg38, mm10), handle chr prefix differences
- Statistical rigor - Proper multiple testing correction, effect size filtering, sample size awareness
- Missing data handling - Explicitly report and handle NaN/missing values
- Chromosome normalization - Always normalize chromosome names (chr1 vs 1, chrX vs X)
- CpG site identification - Parse Illumina probe IDs (cg/ch probes), genomic coordinates
- Report-first - Create output file first, populate progressively
- English-first queries - Use English in all tool calls
Complete Workflow
Phase 0: Question Parsing and Data Discovery
CRITICAL FIRST STEP: Before writing ANY code, parse the question to identify what is being asked and what data files are available.
0.1 Discover Available Data Files
import os
import glob
data_dir = "."
all_files = glob.glob(os.path.join(data_dir, "**/*"), recursive=True)
methylation_files = [f for f in all_files if any(x in f.lower() for x in
['methyl', 'beta', 'cpg', 'illumina', '450k', '850k', 'epic', 'mval'])]
chipseq_files = [f for f in all_files if any(x in f.lower() for x in
['chip', 'peak', 'narrowpeak', 'broadpeak', 'histone'])]
atacseq_files = [f for f in all_files if any(x in f.lower() for x in
['atac', 'accessibility', 'openChromatin', 'dnase'])]
bed_files = [f for f in all_files if f.endswith(('.bed', '.bed.gz', '.narrowPeak', '.broadPeak'))]
bigwig_files = [f f all_files f.endswith((, , ))]
clinical_files = [f f all_files (x f.lower() x
[, , , , , ])]
expression_files = [f f all_files (x f.lower() x
[, , , , , ])]
manifest_files = [f f all_files (x f.lower() x
[, , , ])]
category, files [
(, methylation_files),
(, chipseq_files),
(, atacseq_files),
(, bed_files),
(, bigwig_files),
(, clinical_files),
(, expression_files),
(, manifest_files),
]:
files:
()
0.2 Parse Question Parameters
Extract these from the question:
| Parameter | Default | Example Question Text |
|---|
| Significance threshold | 0.05 | "padj < 0.05", "FDR < 0.01" |
| Beta difference threshold | 0 | " |
| Variance filter | None | "variance > 0.01", "top 5000 most variable" |
| Chromosome filter | All | "chromosome 17", "autosomes only" |
| Genome build | hg38 | "hg19", "GRCh37", "mm10" |
| CpG type filter | All | "cg probes only", "exclude ch probes" |
| Region filter | None | "promoter", "gene body", "intergenic" |
| Missing data handling | Report | "complete cases", "no missing data" |
| Specific comparison | Infer | "tumor vs normal", "old vs young" |
| Specific statistic | Infer | "density", "ratio", "count", "average" |
0.3 Decision Tree
Q: What type of epigenomics data?
METHYLATION -> Phase 1 (Methylation Processing)
CHIP-SEQ -> Phase 2 (ChIP-seq Processing)
ATAC-SEQ -> Phase 3 (ATAC-seq Processing)
MULTI-OMICS -> Phase 4 (Integration)
CLINICAL -> Phase 5 (Clinical Integration)
ANNOTATION -> Phase 6 (ToolUniverse Annotation)
Q: Is this a genome-wide statistics question?
YES -> Focus on chromosome-level aggregation (Phase 7)
NO -> Focus on site/region-level analysis
Phase 1: Methylation Data Processing
1.1 Load Methylation Data
import pandas as pd
import numpy as np
def load_methylation_data(file_path, **kwargs):
"""Load methylation beta-value or M-value matrix.
Expected format:
- Rows: CpG probes (cg00000029, cg00000108, ...)
- Columns: Samples (TCGA-XX-XXXX, ...)
- Values: Beta values (0-1) or M-values (log2 ratio)
"""
ext = os.path.splitext(file_path)[1].lower()
if ext in ['.csv']:
df = pd.read_csv(file_path, index_col=0, **kwargs)
elif ext in ['.tsv', '.txt']:
df = pd.read_csv(file_path, sep='\t', index_col=0, **kwargs)
elif ext in ['.parquet']:
df = pd.read_parquet(file_path, **kwargs)
elif ext in ['.h5', '.hdf5']:
df = pd.read_hdf(file_path, **kwargs)
else:
try:
df = pd.read_csv(file_path, sep='\t', index_col=0, **kwargs)
except Exception:
df = pd.read_csv(file_path, index_col=0, **kwargs)
return df
def detect_methylation_type(df):
"""Detect if data is beta values (0-1) or M-values (unbounded)."""
sample_vals = df.iloc[:1000, :5].values.flatten()
sample_vals = sample_vals[~np.isnan(sample_vals)]
if sample_vals.() >= sample_vals.() <= :
:
():
beta = np.clip(beta, , - )
np.log2(beta / ( - beta))
():
**mvalue / (**mvalue + )
1.2 Load Methylation Manifest / Probe Annotation
def load_probe_annotation(manifest_path):
"""Load Illumina methylation array manifest.
Common columns: IlmnID, Name, CHR, MAPINFO (position), Strand,
UCSC_RefGene_Name, UCSC_RefGene_Group, Relation_to_UCSC_CpG_Island
"""
for skiprows in [0, 7, 8]:
try:
manifest = pd.read_csv(manifest_path, skiprows=skiprows,
low_memory=False)
if 'CHR' in manifest.columns or 'chr' in manifest.columns:
break
if 'Name' in manifest.columns or 'IlmnID' in manifest.columns:
break
except Exception:
continue
col_map = {}
for col in manifest.columns:
lower = col.lower()
if lower in ['chr', 'chromosome']:
col_map[col] = 'chr'
elif lower in ['mapinfo', 'position', 'pos', 'start']:
col_map[col] = 'position'
elif lower in [, , , ]:
col_map[col] =
lower lower:
col_map[col] =
lower:
col_map[col] =
lower lower:
col_map[col] =
manifest = manifest.rename(columns=col_map)
manifest
():
chrom pd.isna(chrom):
chrom = (chrom).strip()
chrom.startswith():
chrom = + chrom
chrom
():
hg38 = {
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
}
hg19 = {
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
}
mm10 = {
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
: , : , : ,
}
genomes = {: hg38, : hg19, : mm10}
genomes.get(genome, hg38)
1.3 CpG Site Filtering
def filter_cpg_probes(df, manifest=None, filters=None):
"""Filter CpG probes based on various criteria.
Args:
df: Methylation matrix (probes x samples)
manifest: Probe annotation DataFrame
filters: dict with keys:
- 'variance_threshold': float, minimum variance across samples
- 'mean_beta_range': tuple (min, max), filter probes with extreme mean beta
- 'missing_threshold': float (0-1), max fraction of NaN allowed per probe
- 'chromosomes': list, keep only these chromosomes
- 'exclude_sex_chr': bool, remove chrX and chrY
- 'probe_type': 'cg' or 'ch', keep only one type
- 'cpg_island': str ('Island', 'Shore', 'Shelf', 'OpenSea')
- 'gene_group': str ('TSS200', 'TSS1500', 'Body', '1stExon', etc.)
- 'top_n_variable': int, keep top N most variable probes
"""
if filters is None:
filters = {}
probe_mask = pd.Series(True, index=df.index)
if 'probe_type' in filters:
ptype = filters['probe_type']
probe_mask &= df.index.str.startswith(ptype)
if 'missing_threshold' in filters:
threshold = filters['missing_threshold']
missing_frac = df.isna().mean(axis=1)
probe_mask &= missing_frac <= threshold
if 'variance_threshold' in filters:
var_threshold = filters['variance_threshold']
probe_var = df.var(axis=1, skipna=True)
probe_mask &= probe_var >= var_threshold
if 'mean_beta_range' filters:
min_beta, max_beta = filters[]
probe_mean = df.mean(axis=, skipna=)
probe_mask &= (probe_mean >= min_beta) & (probe_mean <= max_beta)
filters:
n = filters[]
probe_var = df.var(axis=, skipna=)
top_probes = probe_var.nlargest(n).index
probe_mask &= df.index.isin(top_probes)
manifest (manifest) > :
probe_id_col = manifest.columns manifest.columns[]
manifest_indexed = manifest.set_index(probe_id_col) probe_id_col manifest.columns manifest
filters manifest_indexed.columns:
valid_chr = [normalize_chromosome(c) c filters[]]
chr_probes = manifest_indexed[
manifest_indexed[].apply(normalize_chromosome).isin(valid_chr)
].index
probe_mask &= df.index.isin(chr_probes)
filters.get(, ) manifest_indexed.columns:
sex_chr = [, , , ]
nonsex_probes = manifest_indexed[
~manifest_indexed[].apply(normalize_chromosome).isin([, ])
].index
probe_mask &= df.index.isin(nonsex_probes)
filters manifest_indexed.columns:
relation = filters[]
island_probes = manifest_indexed[
manifest_indexed[]..contains(relation, na=, =)
].index
probe_mask &= df.index.isin(island_probes)
filters manifest_indexed.columns:
group = filters[]
group_probes = manifest_indexed[
manifest_indexed[]..contains(group, na=, =)
].index
probe_mask &= df.index.isin(group_probes)
filtered_df = df[probe_mask]
filtered_df
1.4 Differential Methylation Analysis
from scipy import stats
import statsmodels.stats.multitest as mt
def differential_methylation(beta_df, group1_samples, group2_samples,
test='ttest', correction='fdr_bh', alpha=0.05):
"""Perform differential methylation analysis between two groups.
Args:
beta_df: Beta-value matrix (probes x samples)
group1_samples: list of sample IDs for group 1
group2_samples: list of sample IDs for group 2
test: 'ttest', 'wilcoxon', or 'ks' (Kolmogorov-Smirnov)
correction: multiple testing correction method
alpha: significance threshold
Returns:
DataFrame with columns: mean_g1, mean_g2, delta_beta, pvalue, padj
"""
g1 = beta_df[group1_samples]
g2 = beta_df[group2_samples]
results = []
for probe in beta_df.index:
vals1 = g1.loc[probe].dropna().values
vals2 = g2.loc[probe].dropna().values
if len(vals1) < 2 or len(vals2) < 2:
results.append({
'probe': probe, 'mean_g1': np.nan, 'mean_g2': np.nan,
'delta_beta': np.nan, 'pvalue': np.nan
})
continue
mean1 = np.nanmean(vals1)
mean2 = np.nanmean(vals2)
delta = mean2 - mean1
if test == 'ttest':
stat, pval = stats.ttest_ind(vals1, vals2, equal_var=False)
elif test == 'wilcoxon':
stat, pval = stats.mannwhitneyu(vals1, vals2, alternative='two-sided')
elif test == 'ks':
stat, pval = stats.ks_2samp(vals1, vals2)
else:
stat, pval = stats.ttest_ind(vals1, vals2, equal_var=)
results.append({
: probe, : mean1, : mean2,
: delta, : pval
})
result_df = pd.DataFrame(results).set_index()
valid_pvals = result_df[].dropna()
(valid_pvals) > :
reject, padj, _, _ = mt.multipletests(valid_pvals.values, alpha=alpha, method=correction)
result_df.loc[valid_pvals.index, ] = padj
:
result_df[] = np.nan
result_df
():
dmps = dm_results[
(dm_results[] < alpha) &
(dm_results[].() >= delta_beta_threshold)
].copy()
dmps[] = np.where(dmps[] > , , )
dmps.sort_values()
1.5 Age-Related CpG Analysis
def identify_age_related_cpgs(beta_df, ages, method='correlation',
correction='fdr_bh', alpha=0.05):
"""Identify CpG sites associated with age.
Args:
beta_df: Beta-value matrix (probes x samples)
ages: Series or array of ages corresponding to samples
method: 'correlation' (Pearson/Spearman) or 'regression'
correction: multiple testing method
alpha: significance threshold
Returns:
DataFrame with correlation, p-value, adjusted p-value
"""
results = []
for probe in beta_df.index:
vals = beta_df.loc[probe].values
mask = ~np.isnan(vals) & ~np.isnan(ages.values if hasattr(ages, 'values') else ages)
if sum(mask) < 5:
results.append({'probe': probe, 'correlation': np.nan,
'pvalue': np.nan})
continue
if method == 'correlation':
corr, pval = stats.pearsonr(ages[mask] if hasattr(ages, '__getitem__') else
np.array(ages)[mask], vals[mask])
elif method == 'spearman':
corr, pval = stats.spearmanr(ages[mask] if hasattr(ages, '__getitem__') else
np.array(ages)[mask], vals[mask])
else:
corr, pval = stats.pearsonr(ages[mask] if hasattr(ages, '__getitem__') else
np.array(ages)[mask], vals[mask])
results.append({: probe, : corr, : pval})
result_df = pd.DataFrame(results).set_index()
valid_pvals = result_df[].dropna()
(valid_pvals) > :
reject, padj, _, _ = mt.multipletests(valid_pvals.values, alpha=alpha, method=correction)
result_df.loc[valid_pvals.index, ] = padj
:
result_df[] = np.nan
result_df
1.6 Chromosome-Level Methylation Statistics
def chromosome_cpg_density(cpg_probes, manifest, genome='hg38'):
"""Calculate CpG density per chromosome.
Args:
cpg_probes: list/Index of CpG probe IDs
manifest: probe annotation with chr and position columns
genome: genome build for chromosome lengths
Returns:
DataFrame with chr, n_cpgs, chr_length, density (CpGs per bp)
"""
chr_lengths = get_chromosome_lengths(genome)
probe_id_col = 'probe_id' if 'probe_id' in manifest.columns else manifest.columns[0]
if probe_id_col in manifest.columns:
probe_chr = manifest.set_index(probe_id_col)
else:
probe_chr = manifest
if 'chr' in probe_chr.columns:
chr_col = 'chr'
elif 'CHR' in probe_chr.columns:
chr_col = 'CHR'
else:
raise ValueError("No chromosome column found in manifest")
probe_chrs = probe_chr.loc[probe_chr.index.isin(cpg_probes), chr_col]
probe_chrs = probe_chrs.apply(normalize_chromosome)
chr_counts = probe_chrs.value_counts()
results = []
for chrom, count in chr_counts.items():
if chrom in chr_lengths:
length = chr_lengths[chrom]
density = count / length
results.append({
'chr': chrom,
'n_cpgs': count,
'chr_length': length,
'density_per_bp': density,
: density * ,
})
pd.DataFrame(results).sort_values(,
key= x: x..replace(, ).replace({: , : }).astype())
():
total_cpgs = density_df[].()
total_length = density_df[].()
total_cpgs / total_length
():
chr1 = normalize_chromosome(chr1)
chr2 = normalize_chromosome(chr2)
d1 = density_df[density_df[] == chr1][].values[]
d2 = density_df[density_df[] == chr2][].values[]
d1 / d2
Phase 2: ChIP-seq Peak Analysis
2.1 Load BED/Peak Files
def load_bed_file(file_path, format='bed'):
"""Load BED format file (standard BED, narrowPeak, broadPeak).
Standard BED: chrom, start, end, name, score, strand
narrowPeak: + signalValue, pValue, qValue, peak
broadPeak: + signalValue, pValue, qValue
"""
if format == 'narrowPeak' or file_path.endswith('.narrowPeak'):
names = ['chrom', 'start', 'end', 'name', 'score', 'strand',
'signalValue', 'pValue', 'qValue', 'peak']
elif format == 'broadPeak' or file_path.endswith('.broadPeak'):
names = ['chrom', 'start', 'end', 'name', 'score', 'strand',
'signalValue', 'pValue', 'qValue']
else:
with open(file_path, 'r') as f:
first_line = f.readline().strip()
while first_line.startswith('#') or first_line.startswith('track') or first_line.startswith('browser'):
first_line = f.readline().strip()
n_cols = (first_line.split())
bed_col_names = [, , , , , ,
, , , ,
, ]
names = bed_col_names[:n_cols]
df = pd.read_csv(file_path, sep=, header=, names=names,
comment=, low_memory=)
df = df[~df[].astype()..startswith((, ))]
df[] = df[].apply(normalize_chromosome)
df[] = pd.to_numeric(df[], errors=)
df[] = pd.to_numeric(df[], errors=)
df
():
peaks_df = peaks_df.copy()
peaks_df[] = peaks_df[] - peaks_df[]
stats_dict = {
: (peaks_df),
: peaks_df[].mean(),
: peaks_df[].median(),
: peaks_df[].(),
: peaks_df[].value_counts().to_dict(),
}
peaks_df.columns:
stats_dict[] = peaks_df[].mean()
stats_dict[] = peaks_df[].median()
peaks_df.columns:
stats_dict[] = peaks_df[].mean()
stats_dict
2.2 Peak Annotation
def annotate_peaks_to_genes(peaks_df, gene_annotation=None,
tss_upstream=2000, tss_downstream=500):
"""Annotate peaks to nearest gene / genomic feature.
Args:
peaks_df: BED DataFrame
gene_annotation: DataFrame with gene coordinates (chr, start, end, gene_name, strand)
tss_upstream: bp upstream of TSS to define promoter
tss_downstream: bp downstream of TSS to define promoter
Returns:
DataFrame with peak annotations
"""
if gene_annotation is None:
return peaks_df
annotated = peaks_df.copy()
annotations = []
for _, peak in peaks_df.iterrows():
peak_chr = peak['chrom']
peak_mid = (peak['start'] + peak['end']) // 2
chr_genes = gene_annotation[gene_annotation['chr'] == peak_chr]
if len(chr_genes) == 0:
annotations.append({
'nearest_gene': 'intergenic',
'distance_to_tss': np.nan,
'feature': 'intergenic'
})
continue
tss_positions = chr_genes.apply(
lambda g: g['start'] if g.get('strand', '+') == '+' else g['end'],
axis=1
)
distances = (peak_mid - tss_positions).()
nearest_idx = distances.idxmin()
nearest_gene = chr_genes.loc[nearest_idx]
distance = distances.loc[nearest_idx]
tss = tss_positions.loc[nearest_idx]
(peak_mid - tss) <= tss_upstream:
feature =
peak[] >= nearest_gene[] peak[] <= nearest_gene[]:
feature =
(peak_mid - tss) <= :
feature =
:
feature =
annotations.append({
: nearest_gene.get(, nearest_gene.name),
: (distance),
: feature
})
ann_df = pd.DataFrame(annotations, index=peaks_df.index)
pd.concat([peaks_df, ann_df], axis=)
():
annotated_peaks.columns:
{: (annotated_peaks)}
annotated_peaks[].value_counts().to_dict()
2.3 Peak Overlap Analysis
def find_overlaps(peaks_a, peaks_b, min_overlap=1):
"""Find overlapping peaks between two BED DataFrames.
Uses a simple interval overlap approach (pure Python, no pybedtools).
Args:
peaks_a: BED DataFrame (query)
peaks_b: BED DataFrame (subject)
min_overlap: minimum overlap in bp
Returns:
DataFrame of overlapping pairs
"""
overlaps = []
for chrom in peaks_a['chrom'].unique():
a_chr = peaks_a[peaks_a['chrom'] == chrom].sort_values('start')
b_chr = peaks_b[peaks_b['chrom'] == chrom].sort_values('start')
if len(b_chr) == 0:
continue
for _, a_peak in a_chr.iterrows():
for _, b_peak in b_chr.iterrows():
if b_peak['start'] >= a_peak['end']:
break
if b_peak['end'] <= a_peak['start']:
continue
overlap_start = max(a_peak['start'], b_peak['start'])
overlap_end = min(a_peak['end'], b_peak['end'])
overlap_bp = overlap_end - overlap_start
if overlap_bp >= min_overlap:
overlaps.append({
'a_chrom': chrom,
: a_peak[],
: a_peak[],
: b_peak[],
: b_peak[],
: overlap_bp,
})
pd.DataFrame(overlaps) overlaps pd.DataFrame()
():
chr_lengths = get_chromosome_lengths(genome)
total_genome = (chr_lengths.values())
coverage_a = (peaks_a[] - peaks_a[]).()
coverage_b = (peaks_b[] - peaks_b[]).()
overlaps = find_overlaps(peaks_a, peaks_b)
(overlaps) == :
intersection = overlaps[].()
union = coverage_a + coverage_b - intersection
intersection / union union >
Phase 3: ATAC-seq Analysis
3.1 ATAC-seq Peak Processing
def load_atac_peaks(file_path):
"""Load ATAC-seq peak file (typically narrowPeak format)."""
return load_bed_file(file_path, format='narrowPeak')
def atac_peak_statistics(peaks_df):
"""ATAC-seq specific statistics.
ATAC-seq peaks represent open chromatin regions.
"""
basic_stats = peak_statistics(peaks_df)
peaks_df = peaks_df.copy()
peaks_df['length'] = peaks_df['end'] - peaks_df['start']
nfr_peaks = peaks_df[peaks_df['length'] < 150]
nucleosome_peaks = peaks_df[peaks_df['length'] >= 150]
basic_stats['nfr_peaks'] = len(nfr_peaks)
basic_stats['nucleosome_peaks'] = len(nucleosome_peaks)
basic_stats['nfr_fraction'] = len(nfr_peaks) / len(peaks_df) if len(peaks_df) > 0 else 0
return basic_stats
def chromatin_accessibility_by_region(peaks_df, gene_annotation=None):
"""Calculate chromatin accessibility distribution across genomic regions."""
annotated = annotate_peaks_to_genes(peaks_df, gene_annotation)
regions = classify_peak_regions(annotated)
total = sum(regions.values())
region_fractions = {k: v / total for k, v in regions.items()}
{
: regions,
: region_fractions,
: total,
}
Phase 4: Multi-Omics Integration
4.1 Expression-Methylation Correlation
def correlate_methylation_expression(beta_df, expression_df, probe_gene_map,
method='pearson', correction='fdr_bh'):
"""Correlate methylation levels with gene expression.
Args:
beta_df: Methylation matrix (probes x samples)
expression_df: Expression matrix (genes x samples)
probe_gene_map: dict or Series mapping probe IDs to gene symbols
method: 'pearson' or 'spearman'
correction: multiple testing correction
Returns:
DataFrame with correlation, p-value per probe-gene pair
"""
common_samples = list(set(beta_df.columns) & set(expression_df.columns))
if len(common_samples) < 5:
raise ValueError(f"Not enough common samples: {len(common_samples)}")
beta_aligned = beta_df[common_samples]
expr_aligned = expression_df[common_samples]
results = []
for probe, gene in probe_gene_map.items():
if probe not in beta_aligned.index or gene not in expr_aligned.index:
continue
meth_vals = beta_aligned.loc[probe].values
expr_vals = expr_aligned.loc[gene].values
mask = ~np.isnan(meth_vals) & ~np.isnan(expr_vals)
if sum(mask) < 5:
continue
if method == 'pearson':
corr, pval = stats.pearsonr(meth_vals[mask], expr_vals[mask])
else:
corr, pval = stats.spearmanr(meth_vals[mask], expr_vals[mask])
results.append({
'probe': probe,
'gene': gene,
: corr,
: pval,
: (mask),
})
result_df = pd.DataFrame(results)
(result_df) > :
valid_pvals = result_df[].dropna()
(valid_pvals) > :
reject, padj, _, _ = mt.multipletests(valid_pvals.values, method=correction)
result_df.loc[valid_pvals.index, ] = padj
result_df
4.2 ChIP-seq + Expression Integration
def integrate_chipseq_expression(peaks_df, expression_df, gene_annotation,
tss_window=5000):
"""Integrate ChIP-seq peaks with gene expression.
Args:
peaks_df: ChIP-seq peaks (BED)
expression_df: Gene expression (genes x samples)
gene_annotation: Gene coordinates
tss_window: window around TSS for promoter peaks
Returns:
DataFrame with genes having promoter peaks and their expression
"""
annotated = annotate_peaks_to_genes(peaks_df, gene_annotation,
tss_upstream=tss_window)
promoter_peaks = annotated[annotated['feature'] == 'promoter']
peak_genes = promoter_peaks['nearest_gene'].unique()
common_genes = [g for g in peak_genes if g in expression_df.index]
result = pd.DataFrame({
'gene': common_genes,
'has_promoter_peak': True,
'mean_expression': [expression_df.loc[g].mean() for g in common_genes],
})
return result
Phase 5: Clinical Data Integration
5.1 Missing Data Analysis
def missing_data_analysis(clinical_df=None, expression_df=None,
methylation_df=None, sample_id_col=None):
"""Analyze missing data across multiple omics modalities.
For BixBench questions like:
"How many patients have no missing data for vital status, gene expression, and methylation?"
Args:
clinical_df: Clinical data (patients x variables)
expression_df: Expression matrix (genes x samples)
methylation_df: Methylation matrix (probes x samples)
sample_id_col: Column name for sample/patient IDs in clinical data
Returns:
dict with completeness statistics
"""
results = {}
clinical_samples = set()
if clinical_df is not None:
if sample_id_col and sample_id_col in clinical_df.columns:
clinical_samples = set(clinical_df[sample_id_col].dropna())
else:
clinical_samples = set(clinical_df.index)
results['clinical_samples'] = len(clinical_samples)
expression_samples = set()
if expression_df is not None:
expression_samples = set(expression_df.columns)
results['expression_samples'] = len(expression_samples)
methylation_samples = set()
if methylation_df is not None:
methylation_samples = set(methylation_df.columns)
results['methylation_samples'] = (methylation_samples)
all_sets = []
clinical_samples:
all_sets.append(clinical_samples)
expression_samples:
all_sets.append(expression_samples)
methylation_samples:
all_sets.append(methylation_samples)
(all_sets) > :
complete_samples = .intersection(*all_sets)
results[] = (complete_samples)
results[] = (complete_samples)
:
results[] =
clinical_df :
col clinical_df.columns:
n_missing = clinical_df[col].isna().()
n_total = (clinical_df)
results[] = n_missing
results[] = n_total - n_missing
results
():
sample_sets = []
name, df data_frames.items():
df :
variables name variables:
var variables[name]:
var df.columns:
complete = (df[df[var].notna()].index)
sample_sets.append(complete)
var df.index:
complete = (df.columns[df.loc[var].notna()])
sample_sets.append(complete)
:
sample_sets.append((df.columns))
sample_sets:
()
.intersection(*sample_sets)
Phase 6: ToolUniverse Annotation Integration
Use ToolUniverse tools for biological annotation of epigenomic findings.
6.1 Gene-Level Annotation
from tooluniverse import ToolUniverse
tu = ToolUniverse()
tu.load_tools()
def annotate_genes_with_tooluniverse(gene_list, tu):
"""Annotate a list of genes using ToolUniverse tools.
Uses:
- Ensembl for gene coordinates and cross-references
- SCREEN for regulatory elements near gene
- ChIPAtlas for ChIP-seq experiments
"""
annotations = {}
for gene in gene_list[:20]:
annotation = {'gene': gene}
try:
ens = tu.tools.ensembl_lookup_gene(id=gene, species='homo_sapiens')
if isinstance(ens, dict):
data = ens.get('data', ens)
annotation['ensembl_id'] = data.get('id', 'N/A')
annotation['chr'] = data.get('seq_region_name', 'N/A')
annotation['start'] = data.get('start', 'N/A')
annotation['end'] = data.get('end', 'N/A')
annotation['biotype'] = data.get('biotype', 'N/A')
except Exception:
pass
try:
screen = tu.tools.SCREEN_get_regulatory_elements(
gene_name=gene, element_type="enhancer", limit=
)
screen :
annotation[] =
Exception:
annotations[gene] = annotation
pd.DataFrame.from_dict(annotations, orient=)
6.2 ChIPAtlas Integration
def query_chipatlas_experiments(antigen, genome='hg38', cell_type=None, tu=None):
"""Query ChIPAtlas for available ChIP-seq experiments.
Args:
antigen: TF or histone mark name (e.g., 'H3K27ac', 'CTCF')
genome: genome build
cell_type: optional cell type filter
Returns:
ChIPAtlas experiment metadata
"""
if tu is None:
from tooluniverse import ToolUniverse
tu = ToolUniverse()
tu.load_tools()
params = {
'operation': 'get_experiment_list',
'genome': genome,
'antigen': antigen,
'limit': 50,
}
if cell_type:
params['cell_type'] = cell_type
return tu.tools.ChIPAtlas_get_experiments(**params)
6.3 Ensembl Regulatory Feature Annotation
def annotate_regions_with_ensembl(regions, species='human', tu=None):
"""Annotate genomic regions with Ensembl regulatory features.
Args:
regions: list of (chr, start, end) tuples
species: Ensembl species name
Returns:
dict of region -> regulatory features
"""
if tu is None:
from tooluniverse import ToolUniverse
tu = ToolUniverse()
tu.load_tools()
annotations = {}
for chrom, start, end in regions[:10]:
ens_chrom = chrom.replace('chr', '') if chrom.startswith('chr') else chrom
region_str = f"{ens_chrom}:{start}-{end}"
try:
result = tu.tools.ensembl_get_regulatory_features(
region=region_str, feature="regulatory", species=species
)
annotations[(chrom, start, end)] = result
except Exception as e:
annotations[(chrom, start, end)] = {'error': str(e)}
return annotations
Phase 7: Genome-Wide Statistics
7.1 Comprehensive Genome Statistics
def genome_wide_methylation_stats(beta_df, manifest=None, genome='hg38'):
"""Calculate comprehensive genome-wide methylation statistics.
Args:
beta_df: Methylation matrix (probes x samples)
manifest: Probe annotation
genome: genome build
Returns:
dict with genome-wide statistics
"""
stats_result = {
'total_probes': len(beta_df),
'total_samples': beta_df.shape[1],
'global_mean_beta': float(beta_df.mean().mean()),
'global_median_beta': float(beta_df.median().median()),
'global_std_beta': float(beta_df.values[~np.isnan(beta_df.values)].std()),
'missing_fraction': float(beta_df.isna().mean().mean()),
}
stats_result['sample_means'] = beta_df.mean(axis=0).describe().to_dict()
probe_var = beta_df.var(axis=1, skipna=True)
stats_result['probe_variance'] = {
'mean': float(probe_var.mean()),
'median': float(probe_var.median()),
'max': float(probe_var.max()),
}
stats_result['high_variance_probes'] = int((probe_var > 0.01).sum())
if manifest is not None:
density_df = chromosome_cpg_density(beta_df.index.tolist(), manifest, genome)
stats_result[] = density_df.to_dict()
stats_result[] = genome_wide_average_density(density_df)
stats_result
():
sig = dm_results[dm_results[] < alpha]
hyper = sig[sig[] > ]
hypo = sig[sig[] < ]
{
: (dm_results),
: (sig),
: (hyper),
: (hypo),
: (sig) / (dm_results) (dm_results) > ,
: (sig[].mean()) (sig) > ,
: (sig[].().()) (sig) > ,
}
ToolUniverse Tool Parameter Reference
Regulatory Annotation Tools
| Tool | Parameters | Returns |
|---|
ensembl_lookup_gene | id: str, species: str (REQUIRED) | {status, data: {id, display_name, seq_region_name, start, end, strand, biotype}, url} |
ensembl_get_regulatory_features | region: str (no "chr"), feature: str, species: str | {status, data: [...features...]} |
ensembl_get_overlap_features | region: str, feature: str, species: str | Gene/transcript overlap data |
SCREEN_get_regulatory_elements | gene_name: str, element_type: str, limit: int | cCREs (enhancers, promoters, insulators) |
ReMap_get_transcription_factor_binding | gene_name: str, cell_type: str, limit: int | TF binding sites |
RegulomeDB_query_variant | rsid: str | {status, data, url} regulatory score |
jaspar_search_matrices | search: str, collection: str, species: str | {count, results: [...matrices...]} |
ENCODE_search_experiments | assay_title: str, target: str, organism: str, limit: int | Experiment metadata |
ChIPAtlas_get_experiments | operation: str (REQUIRED: "get_experiment_list"), genome: str, antigen: str, cell_type: str, limit: int | Experiment list |
ChIPAtlas_search_datasets |
Gene Annotation Tools
| Tool | Parameters | Returns |
|---|
MyGene_query_genes | query: str | {hits: [{_id, symbol, ensembl, ...}]} |
MyGene_batch_query | gene_ids: list[str], fields: str | {results: [{query, symbol, ...}]} |
HGNC_get_gene_info | symbol: str | Gene symbol, aliases, IDs |
GO_get_annotations_for_gene | gene_id: str | GO annotations |
CRITICAL Tool Notes
- ensembl_lookup_gene: REQUIRES
species='homo_sapiens' parameter
- ensembl_get_regulatory_features: Region format is "17:start-end" (NO "chr" prefix)
- ChIPAtlas tools: ALL require
operation parameter (SOAP-style)
- FourDN tools: ALL require
operation parameter (SOAP-style)
- SCREEN: Returns JSON-LD format with
@context, @graph keys
Response Format Notes
- Methylation data: Typically stored as probes (rows) x samples (columns), beta values 0-1
- BED files: Tab-separated, 0-based half-open coordinates
- narrowPeak: 10-column BED extension with signalValue, pValue, qValue, peak
- Illumina manifests: Contains probe ID, chromosome, position, gene annotation
- Clinical data: Patient/sample-centric with clinical variables as columns
Fallback Strategies
| Scenario | Primary | Fallback |
|---|
| No manifest file | Load from data dir | Build minimal from Ensembl lookup |
| No pybedtools | Pure Python overlap | pandas-based interval intersection |
| No pyBigWig | Skip BigWig analysis | Use pre-computed summary tables |
| Missing clinical data | Report missing | Use available samples only |
| Low sample count | Parametric test | Use non-parametric (Wilcoxon) |
| Large dataset (>500K probes) | Full analysis | Sample or chunk-based processing |
Common Use Patterns
Pattern 1: Methylation Array Analysis
Input: Beta-value matrix + manifest + clinical data
Question: "How many CpGs are differentially methylated?"
Flow:
1. Load beta matrix, manifest, clinical data
2. Filter CpG probes (cg only, remove sex chr, variance filter)
3. Define groups from clinical data
4. Run differential_methylation()
5. Apply thresholds (padj < 0.05, |delta_beta| > 0.2)
6. Report count and direction (hyper/hypo)
Pattern 2: Age-Related CpG Density
Input: Beta-value matrix + manifest + ages
Question: "What is the density ratio of age-related CpGs between chr1 and chr2?"
Flow:
1. Load beta matrix and ages from clinical data
2. Run identify_age_related_cpgs()
3. Filter significant age-related CpGs
4. Map to chromosomes using manifest
5. Calculate chromosome_cpg_density()
6. Compute ratio between specified chromosomes
Pattern 3: Multi-Omics Missing Data
Input: Clinical + expression + methylation data files
Question: "How many patients have complete data for all modalities?"
Flow:
1. Load all data files
2. Extract sample IDs from each
3. Find intersection (common samples)
4. Check for NaN/missing within clinical variables
5. Report complete cases count
Pattern 4: ChIP-seq Peak Annotation
Input: BED/narrowPeak file
Question: "What fraction of peaks are in promoter regions?"
Flow:
1. Load BED file with load_bed_file()
2. Load or fetch gene annotation (Ensembl)
3. Run annotate_peaks_to_genes()
4. Classify regions with classify_peak_regions()
5. Calculate fraction in promoters
Pattern 5: Methylation-Expression Integration
Input: Beta matrix + expression matrix + probe-gene mapping
Question: "What is the correlation between methylation and expression?"
Flow:
1. Load both matrices
2. Build probe-gene map from manifest
3. Align samples across datasets
4. Run correlate_methylation_expression()
5. Report significant anti-correlations
Edge Cases
Missing Probe Annotation
When no manifest/annotation file is available:
- Extract chromosome from probe ID naming patterns if possible
- Use ToolUniverse Ensembl tools to build minimal annotation
- Report limitation: "chromosome mapping unavailable for X probes"
Mixed Genome Builds
When data uses different builds:
- Detect build from context (data README, file names, known coordinates)
- Use appropriate chromosome lengths for density calculations
- Do NOT mix hg19 and hg38 coordinates
Very Large Datasets
For datasets with >500K CpG sites:
- Use chunked processing for differential methylation
- Pre-filter by variance before statistical testing
- Use vectorized operations (avoid row-by-row loops where possible)
Sample ID Mismatches
Clinical and molecular data may use different ID formats:
- TCGA: barcode (TCGA-XX-XXXX-01A) vs patient ID (TCGA-XX-XXXX)
- Try truncating or matching partial IDs
- Report number of matched/unmatched samples
Limitations
- No native pybedtools: Uses pure Python interval operations (slower for very large BED files)
- No native pyBigWig: Cannot read BigWig files directly without package
- No R bridge: Does not use methylKit, ChIPseeker, or DiffBind
- Illumina-centric: Methylation functions designed for 450K/EPIC arrays
- Statistical simplicity: Uses t-test/Wilcoxon for differential methylation (not limma/bumphunter)
- No peak calling: Assumes peaks are pre-called; does not run MACS2 or similar
- API rate limits: ToolUniverse annotation limited to ~20 genes per batch
Summary
Genomics & Epigenomics Data Processing Skill provides:
- Methylation analysis - Beta-value processing, CpG filtering, differential methylation, age-related CpGs, chromosome density
- ChIP-seq analysis - Peak loading (BED/narrowPeak), peak annotation, overlap analysis, statistics
- ATAC-seq analysis - Chromatin accessibility peaks, NFR detection, region classification
- Multi-omics integration - Methylation-expression correlation, ChIP-seq + expression
- Clinical integration - Missing data analysis across modalities, complete case identification
- Genome-wide statistics - Chromosome-level density, genome-wide averages, density ratios
- ToolUniverse annotation - Ensembl, SCREEN, ChIPAtlas, JASPAR, ENCODE for biological context
Core packages: pandas, numpy, scipy, statsmodels, pysam, gseapy
ToolUniverse tools: 25+ tools across Ensembl, SCREEN, ENCODE, ChIPAtlas, JASPAR, ReMap, RegulomeDB, 4DN
Best for: BixBench-style quantitative questions about methylation data, ChIP-seq peaks, chromatin accessibility, and multi-omics integration
