functional-enrichment

# Functional Enrichment (HOMER + R)

## Overview

- **Validate input**: Accept BED/peak files with genomic coordinates or gene lists; check format and genome assembly.
- **Map regions to genes**: Convert regions to a unique gene set using HOMER `annotatePeaks.pl`.
- **Run GO enrichment**: Use HOMER `findGO.pl` (or `annotatePeaks.pl -go`) for BP/MF/CC.
- **Run KEGG enrichment**: Use HOMER `findGO.pl -kegg` (or `annotatePeaks.pl -kegg`).
- **Collect outputs**: Save tidy tables for downstream plotting and a compact summary of top terms.
- **Visualize in R**: Create barplots and dotplots (GO/KEGG) with `ggplot2` from standardized outputs.
- **QC & troubleshooting**: Provide checks for genome mismatch, chromosome naming, and low-signal inputs.

## Inputs & Outputs

### Inputs (choose one):
#### Option 1: Input is a genomic region file (BED/narrowPeak/broadPeak)
Genomic region formats supported:
- **BED files**: Standard genomic interval format
- **narrowPeak**: narrow peak format
- **broadPeak**: broad peak format

#### Option 2: Input is a gene list (txt)
- `gene_list.txt` with one official gene symbol per line (no header). And an optional `gene_list_background.txt` with one official gene symbol per line (no header).

### Outputs (directory layout):
```bash
${sample}_functional_enrichment/
    results/
      ${sample}.anno_genomic_features.txt
      ${sample}.anno_genomic_features_stats.txt
      biological_process.txt
      cellular_component.txt  
      molecular_function.txt  

      kegg.txt                
      biocyc.txt              
      chromosome.txt  
      cosmic.txt
      interactions.txt  
      interpro.txt
      gene3d.txt
      pathwayInteractionDB.txt
      pfam.txt
      prints.txt    
      prosite.txt   
      reactome.txt
      smpdb.txt
      wikipathways.txt

      gwas.txt          
      lipidmaps.txt           
      msigdb.txt                
      smart.txt

    tables/
      ${sample}.gene_list.txt
      go_bp.tsv
      go_mf.tsv
      go_cc.tsv
      kegg.tsv
    logs/
      ${sample}.anno_genomic_features.log # if genome region file is provided
      findGO.log
```


## Decision Tree


### Step 0 — Gather Required Information from the User

Before calling any tool, **ask the user**:

1. Sample name (`sample`): used as prefix and for the output directory `${sample}_functional_enrichment`.
2. Genome assembly (`genome`): e.g. `hg38`, `mm10`, `danRer11`.  
   - **Never** guess or auto-detect.

---

### Step 1: Initialize Project

1. Make director for this project:

Call:
- `mcp__project-init-tools__project_init`

with:
- `sample`: the user-provided sample name
- `task`: de_novo_motif_discovery

The tool will:
- Create `${sample}_functional_enrichment` directory.
- Get the full path of the `${sample}_functional_enrichment` directory, which will be used as `${proj_dir}`.

---


### Step 2: Prepare genome file for homer

Call:
- `mcp__homer-tools__check_genome_installation`

With:
- `genome`: the user-provided genome assembly, e.g. `hg38`, `mm10`, `danRer11`

The tool will:
- Check if the genome is installed in HOMER.
- If not, install the genome.


---


### Step 3 (Optional): Standardize chromosome names for BED files

This step is optional. Only perform this step if the input file is a BED file. If the input file is a gene list, skip this step.

From `1` format to `chr1` format
From `MT` format to `chrM` format

Call:
- `mcp__file-format-tools__standardize_bed_chrom_names`

with:
- `input_bed`: the user-provided BED file
- `output_bed`: the path to save the standardized BED file

The tool will:
- Standardize the chromosome names in the BED file.
- Return the path of the standardized BED file.


---

### Step 4 (Optional): Convert gene ID to gene symbol

This step is optional. Only perform this step if the input file is a gene list file. If the input file is a BED file, skip this step.

Call:
- `mcp__mygene-tools__convert_gene_ids_mygene`

With:
- `input_ids_file`: the user-provided gene list file. May end with `.txt`.
- `scopes`: the source ID type for mygene (e.g., 'ensembl.gene', 'symbol', 'entrezgene', 'uniprot', or a comma-separated list).
- `fields`: the comma-separated target fields to retrieve from mygene (e.g., 'symbol,ensembl.gene,uniprot,entrezgene').
- `species`: the species for mygene (e.g., 'human', 'mouse', 'zebrafish', or NCBI taxon ID like '9606').
- `out_file`: the path to save the converted gene list file. In this skill, it is the full path of the `${sample}_functional_enrichment` directory returned by `mcp__project-init-tools__project_init`
- `batch_size`: the batch size for mygene.querymany (default 1000).

The tool will:
- Convert the gene ID to gene symbol.
- Return the path of the converted gene list file.

---


### Step 5: GO enrichment analysis

#### Option 1: from genomic regions file

Only if the input file is a BED file. If the input file is a gene list, call tools in Option 2.

1. annotate the genomic regions using Homer's `annotatePeaks.pl` with `-go` option. If user also provides a background genome region file, like a control peak file, also call this tool for the background genome region file. Use a different `${sample}` as the sample name for the background sample.

Call:
`mcp__homer-tools__annotate_genomic_features`

With:
- `sample`: the user-provided sample name
- `proj_dir`: directory to save the genomic feature annotation results. In this skill, it is the full path of the `${sample}_functional_enrichment` directory returned by `mcp__project-init-tools__project_init`
- `regions_bed`: the user-provided regions file in BED format. May end with `.bed`, `.narrowPeak`, `.broadPeak`, etc.
- `genome`: the user-provided genome assembly, e.g. `hg38`, `mm10`, `danRer11`
- `ann`: "custom homer annotation file (created by assignGenomeAnnotation.pl), (default: None).
- `size_given`: keep original region sizes (default: True)
- `cpg`: include CpG information (default: False)
- `go`: `True` to perform GO enrichment analysis.

The tool will:
- Annotate the genomic regions using Homer's `annotatePeaks.pl`.
- Return the path of the annotated regions file under `${proj_dir}/results/` directory, and the path to the log file under `${proj_dir}/logs/` directory.
    - `${proj_dir}/results/${sample}.anno_genomic_features.txt`
    - `${proj_dir}/results/${sample}.anno_genomic_features_stats.txt`
    - `${proj_dir}/logs/${sample}.anno_genomic_features.log`

---

2. (optional) extract the genes from the annotated regions file if neccessary for future analysis or the target gene list is requested by user. If not requested, skip this step.

Call:
`mcp__file-format-tools__extract_gene_list`

With:
- `sample`: the user-provided sample name
- `proj_dir`: directory to save the genomic feature annotation results. In this skill, it is the full path of the `${sample}_functional_enrichment` directory returned by `mcp__project-init-tools__project_init`

The tool will:
- Extract the genes from the annotated regions file.
- Return the path of the gene list file under `${proj_dir}/tables/` directory.
    - `${proj_dir}/tables/${sample}.gene_list.txt`


---

#### Option 2: from gene list file

Only if the input file is a gene list file. If the input file is a BED file, call tools in Option 1.

Call:
`mcp__homer-tools__gene_function_enrichment`

With:
- `sample`: the user-provided sample name
- `proj_dir`: directory to save the GO & KEGG enrichment results. In this skill, it is the full path of the `${sample}_functional_enrichment` directory returned by `mcp__project-init-tools__project_init`
- `gene_list_file`: the user-provided gene list file. May end with `.txt`.
- `organism`: the user-provided organism name, e.g. `human`, `mouse`, `zebrafish`, etc.
- `background_gene_list_file`: the user-provided background gene list file. May end with `.txt`. If not provided, set this parameter to `None`.

The tool will:
- Find the GO enrichment for the gene list.
- Return the path of the GO & KEGG enrichment results under `${proj_dir}/results/` directory.
    - `${proj_dir}/results/biological_process.txt`
    - `${proj_dir}/results/kegg.txt`
    - ... other GO and KEGG enrichment results files.
- Return the path of the log file under `${proj_dir}/logs/` directory.
    - `${proj_dir}/logs/${sample}.find_go_and_kegg_enrichment.log`


---




> **Alternative direct from BED**  
> `annotatePeaks.pl peaks.bed hg38 -go results/{run}/tables/go_dir -genomeOntology`  
> `annotatePeaks.pl peaks.bed hg38 -kegg results/{run}/tables/kegg_dir`


## Notes & Best Practices

- **Genome & naming**: Ensure the HOMER genome key matches the species; chromosome naming must be consistent (`chr1` vs `1`).
- **BED format**: Tab-delimited, ≥3 columns, 0-based coordinates, no header.
- **Multiple testing**: Prefer FDR (BH) if provided; otherwise fallback to P-value.
- **Background set**: `-bg` helps reduce bias; choose a reasonable universe (e.g., all expressed or all accessible regions → genes).
- **Direct-from-BED**: `annotatePeaks.pl -go/-kegg` is convenient; the gene-list route yields uniform TSVs for plotting.

## Troubleshooting

- **Many NAs after annotation**: Check genome version, chromosome naming, BED formatting, and headers.
- **Empty/weak enrichment**: Ensure sufficient genes (suggest ≥50), verify species of symbols, tune thresholds or background.
- **Column name drift**: HOMER versions may differ; adjust R column mappings if needed.