Motif Analysis Pipeline

Unified pipeline for discovering transcription factor motifs that contribute to the regulatory activity of CRE tiles, as described in the paper's Motif analysis section.

Overview

For each CRE and fragment, the pipeline:

1. Extracts all 8–14 nt k-mers from each sequence in the tile.
2. Filters k-mers by requiring both presence and absence in at least 20% of sequences.
3. Tests association between k-mer presence and expression (log2 RNA/DNA) using the Brunner-Munzel test (p < 1e-3).
4. Matches significant k-mers to known TF motifs using TomTom (p < 0.01) against both JASPAR CORE 2022 and HOCOMOCO v12 CORE databases.
5. Cross-references the corresponding TFs with expression data from two iNeuron RNA-seq datasets (Cui et al. and Song et al.) and the Allen Brain Atlas Brainspan dataset to prioritize biologically relevant TFs.

Directory Structure

motif_analysis/
├── motif_analysis.py           # Unified analysis script
├── run_motif_analysis.sh             # Runner script (all CREs or single CRE/fragment)
├── README.md          # This file
├── data/              # Input data (see "Input Data Requirements" below)
│   ├── mpra_sequences.txt                   # MPRA sequence data
│   ├── JASPAR_CORE_2022.meme                # JASPAR motif database
│   ├── HOCOMOCO_v12_CORE.meme               # HOCOMOCO v12 CORE motif database
│   ├── HOCOMOCO_v12_CORE_annotation.jsonl   # HOCOMOCO annotation file
│   ├── brainspan_rows_metadata.csv          # Brainspan gene metadata
│   ├── brainspan_expression_matrix.csv      # Brainspan expression matrix
│   ├── rnaseq_cui_rep1.tsv                  # Cui et al. replicate 1
│   ├── rnaseq_cui_rep2.tsv                  # Cui et al. replicate 2
│   ├── rnaseq_song_rep1.tsv                 # Song et al. replicate 1
│   └── rnaseq_song_rep2.tsv                 # Song et al. replicate 2
└── output/            # Results (created automatically per CRE)
    └── <cre>/
        ├── <cre>frag<N>_results.csv         # Final results table
        ├── plots_<cre>frag<N>.pdf           # KDE plots (present vs absent)
        ├── motifs.txt / motifs.meme         # Intermediate motif files
        ├── combined_genes.txt               # Candidate gene list
        ├── tomtom_out_<cre>frag<N>_jaspar/  # TomTom JASPAR output
        └── tomtom_out_<cre>frag<N>_h12core/ # TomTom HOCOMOCO output

Usage

Run all CREs, all fragments (1–11)

./run_motif_analysis.sh

Run a single CRE (all fragments)

./run_motif_analysis.sh cre6

Run a single CRE + fragment

./run_motif_analysis.sh cre6 3

Run directly with Python

python3 motif_analysis.py <cre> <frag>
# e.g.
python3 motif_analysis.py creA 7

Supported CRE names: cre2, cre6, creA, mecp2prom

Environment

Dependencies

• Python 3
• pandas, numpy, scipy, matplotlib, seaborn
• beautifulsoup4
• mygene
• MEME Suite — the tomtom binary must be on your PATH (or set the path manually in motif_analysis.py)

Input Data Requirements

To run this pipeline with a different experiment, the following files must be provided in the data/ directory. Look for the ### USER CONFIG sections in motif_analysis.py for exactly where each file is referenced.

1. MPRA Sequences (mpra_sequences.txt)

Tab-separated file with at minimum the following columns:

Column	Description
mutation_ref_seq	DNA sequence of the mutant tile
log2	log2(RNA/DNA) expression value
group_id	Identifier linking sequences to CREs/fragments (e.g., cre6frag1)

The group_id column is used to filter data for a specific CRE and fragment. Its format must be <cre_name>frag<fragment_number>.

2. Motif Databases (MEME format)

• JASPAR_CORE_2022.meme — JASPAR CORE motif database in MEME format. Download from JASPAR.
• HOCOMOCO_v12_CORE.meme — HOCOMOCO v12 CORE motif database in MEME format. Download from HOCOMOCO.
• HOCOMOCO_v12_CORE_annotation.jsonl — HOCOMOCO annotation file (JSONL format). Each line is a JSON object containing name and masterlist_info.species.HUMAN.gene_symbol.

These can be replaced with other motif databases (in MEME format), but the annotation file and gene name extraction logic may need adjustment.

3. RNA-seq Expression Data (TSV)

Two RNA-seq datasets, each with 2 replicates. All files are tab-separated with at minimum:

Column	Description
gene_id	Gene identifier (see note below)
TPM	Transcripts per million

Important: The script expects Dataset 1 (Cui et al.) to use gene symbols as gene_id (e.g., NRF1), and Dataset 2 (Song et al.) to use Ensembl gene IDs (e.g., ENSG00000169057). If your datasets use a different ID format, update the lookup logic in motif_analysis.py (see ### USER CONFIG sections, specifically REFERENCE_GENE_SYMBOL, REFERENCE_GENE_ENSEMBL, and the exp_datasets loading code).

4. Brainspan Expression Data

Download from the Allen Brain Atlas BrainSpan (RNA-Seq Gencode v10 summarized to genes):

• brainspan_rows_metadata.csv — CSV with columns: gene_symbol, row_num (and optionally gene_id, ensembl_gene_id, entrez_id). This corresponds to rows_metadata.csv in the BrainSpan download.
• brainspan_expression_matrix.csv — CSV with genes as rows and brain regions/timepoints as columns. This corresponds to expression_matrix.csv in the BrainSpan download.

Note: These files are not included in the repository due to their size. Download them from BrainSpan and place them in the data/ directory before running the pipeline.

5. TomTom (MEME Suite)

The script calls TomTom as an external binary. By default it is auto-detected from your PATH. If TomTom is not on your PATH, set the path manually in motif_analysis.py (see ### USER CONFIG):

tomtom_path = "/path/to/meme/bin/tomtom"

6. Analysis Parameters

All key thresholds are defined as named constants at the top of motif_analysis.py (see ### USER CONFIG: Analysis parameters):

K_MIN = 8                           # Minimum k-mer length
K_MAX = 14                          # Maximum k-mer length
PRESENCE_ABSENCE_THRESHOLD = 0.2    # Min fraction present/absent (20%)
BM_PVALUE_THRESHOLD = 1e-3          # Brunner-Munzel significance threshold
TOMTOM_PVALUE_THRESHOLD = 0.01      # TomTom motif match p-value threshold
BRAINSPAN_PERCENTILE = 65           # Percentile cutoff for Brainspan expression filter

Output Format

The main output file (<cre>frag<N>_results.csv) contains one row per identified TF with the following columns:

Column	Description
Gene	TF gene symbol
K-mer	Matched k-mer sequence
Starting position	Average position of the k-mer in the tile
Percentage (absence/presence)	Minority percentage of k-mer presence/absence
Brunner-Munzel p-value	Statistical significance of expression difference
Mean difference (abs)	Absolute mean log2 difference (present vs absent)
Median difference (abs)	Absolute median log2 difference
TOMTOM p-value	TomTom motif match p-value
TOMTOM E-value	TomTom E-value
TOMTOM q-value	TomTom q-value
TPM (Cui et al_1)	TPM in Cui et al. replicate 1
TPM (Cui et al_2)	TPM in Cui et al. replicate 2
TPM (Song et al_1)	TPM in Song et al. replicate 1
TPM (Song et al_2)	TPM in Song et al. replicate 2
TPM (Brainspan)	Mean TPM from Brainspan