RegEvol

RegEvol is a computational framework for detecting directional selection in regulatory sequences through phenotypic predictions and phenotype-to-fitness models. This ongoing project aims to identify signatures of positive selection acting on transcription factor binding sites across mammals and drosophila, using both statistical and evolutionary approaches.

---

Project Overview

The pipeline includes the following key steps:

Regulatory Data Processing

• TF ChIP-seq Peak Calling: Extraction and processing of transcription factor binding sites from publicly available ChIP-seq datasets in mammals and drosophila.
• Homologous Sequence Retrieval: Identification of orthologous regulatory regions using whole-genome alignments from the Zoonomia Consortium.

Predictive Modelling

• gkm-SVM Training: Training of gapped k-mer support vector machine (gkm-SVM) models on species-specific regulatory sequences to predict binding affinity changes.

Positive Selection Detection

• Permutation Test: Application of the method from Liu & Robinson-Rechavi, 2020 to detect extreme change in binding affinity.
• RegEvol Test: A novel maximum likelihood-based test to infer evolutionary regimes of regulatory regions from predicted binding affinity shifts.

Comparative and Functional Analyses

• Simulations: Sequence evolution simulations under controlled parameters to validate the performance and robustness of selection detection methods.
• Conservation Metrics: Comparison of selection signals with evolutionary conservation scores (phastCons and phyloP).
• Functional Enrichment: Gene Ontology enrichment analysis of genes associated with positively selected regulatory regions.
• Divergence vs. Polymorphism: Integration of interspecies divergence and intraspecies polymorphism data to interpret evolutionary dynamics.

---

Repository structure

├── 1_chipseq/                    # Scripts to run nf-core ChIP-seq pipeline
│   ├── run.nextflow.sh           # Main wrapper script
│   └── logs/                     # Logs generated by pipeline runs
│
├── 2_selection_tests/            # Scripts and Snakemake workflow for RegEvol and Permutation tests 
│   ├── Snakemake/                # Snakefile and workflow rules
│   ├── envs/                     # Conda environment YAML files
│   └── scripts/                  # Python/R scripts used by the workflow
│
├── 3_analysis/                   # Scripts for downstream analyses and figures
│
├── config/                       # Global configuration files
│
├── docs/                         # ChiP-seq samples description and download files
│
└── README.md                     # Project documentation

---

Requirements

Before running the pipelines, ensure the following tools are installed on your system:

• Conda / Mamba – for managing environments and dependencies.
• Singularity – required to use Docker/Singularity containers in Snakemake rules.

This workflow also depends on:

• Snakemake : to run the RegEvol pipeline (version 7.24.0 was used in development).
• Nextflow : to run the ChIP-seq peaks calling pipeline (version 22.10 was used in development).

Option 1: Use your existing installation

If you already have Snakemake and Nextflow installed on your system, you can skip creating the bundled environment and run the workflow directly.

Option 2: Create the provided environment

If you don’t, you can create the environment with:

mamba env create -f config/workflows.yaml

---

Usage

1. Run the ChIP calling pipeline (Nextflow)

The workflow uses the nf-core/chipseq (2.1.0) Nextflow pipeline to process ChIP-seq data.

Required inputs

The pipeline requires the following three input files:

1. Reference genome sequence (FASTA format)
   → Located in data/genome_sequences/<species>/.

2. Sample sheet (CSV format)
   → Describes the FASTQ files, sample groups, and metadata (see the nf-core/chipseq sample sheet specification).
   → Located in docs/ChIP-seq/<species>/<sample>_samples_input.csv.

3. Annotations (GTF or GFF file)
   → Located in data/genome_sequences/<species>/.

Data download

If FASTQ files are not already present in data/ChIP-seq/<species>/<sample>/, the pipeline will automatically download them using the list of URLs provided in: docs/ChIP-seq/<sp>/DL_<sample>.txt

Run

./1_chipseq/run.nextflow.sh --sp <species> --sample <sample> [options]

Required arguments:

• --sp → species name (must match those in config/params.sh)
• --sample → sample name (used for result subfolders)

Optional arguments (can be passed in any order using --option value):

• --threads → number of threads to allocate (default: 1)
• --peaksType → "Narrow" for transcription factors, "Broad" for histone marks (default: "Narrow")
• --system → Execution mode: local or SLURM (default: "local")
• --container → Container type: "conda" or "singularity" (default: "conda")
• --resume → "resume" to resume a previous run, "false" to start fresh (default: "false")
• --skip → "true" to skip SPP and MultiQC steps, "false" to run all steps (default: "false")
• --help → display usage and exit

Example:

./1_chipseq/run.nextflow.sh --sp human --sample Wilson --threads 16 --system SLURM --container singularity

2. Run the RegEvol pipeline (Snakemake)

Required inputs

The pipeline requires the following three input files:

1. Whole Genome Alignment (maf.gz format)
   → Expected file is located and named as data/genome_alignments/<species>/ancestral_alignment.maf.gz
   → Contains sequences of the species of interest aligned with an ancestral state. → The name of the Ancestral state is defined by anc_name and can be adapted in config/params.sh.

2. Peaks coordinates (BED format)
   → Located in results/peaks/<species>/<sample>/ the file should be named <TF>.peaks.bed or <TF>.peaks_UCSC_names.bed (see below for chromosome names) → This file is automatically generated if the above ChIP calling pipeline was used

3. BSgenome object (tar.gz format)
   → Used by 2_selection_tests/scripts/generate_negative_sequences.R to obtain a matching negative set of sequences for model training.
   → If BSgenome is publicly available for the species of interest, it will be automatically downloaded; otherwise, it will look in: data/BSgenome/
   → Home-made BSgenomes used in this study are available in Zenodo: https://zenodo.org/records/18379807
   → Guidelines to build a new BSgenome can be found in 2_selection_tests/scripts/utils/build_BSgenome/

Global and Snakemake configurations

Global RegEvol parameters for species, ancestral sequences, chromosome, and annotation source are defined in config/params.sh Snakemake-specific parameters for subworkflows and selection tests are defined in /2_selection_tests/Snakemake/config/

Chromosome names

Chromosome convention names should match between the peaks BED file and the assembly used for the whole-genome alignment (i.e "chr1" for UCSC vs "1" for Ensembl). The convention used in the genome alignment file should be indicated for each species in /2_selection_tests/Snakemake/config/. Currently, using the Zoonomia alignment, the default suffix argument is set to "_UCSC_names" for human, mouse, caroli, spretus and rat. And to "" for all the other species.

The suffix is supposed to match the provided peaks BED file: results/peaks/<species>/<sample>/<TF>.peaks<suffix>.bed. If your BED file doesn't follow the correct convention, you could rename the chromosome and adjust the BED filename yourself. The current pipeline is also able to automatically convert chromosome names from Ensembl to UCSC convention if you have both reference genomes in /data/genome_sequences/{sp}/. In this case, the genome file names should match the arguments UCSC_Assembly and Ensembl_Assembly in /2_selection_tests/Snakemake/config/.
The current default for human is: UCSC_Assembly="hg38.fa.gz" and Ensembl_Assembly="Homo_sapiens.GRCh38.dna_sm.toplevel.fa.gz".

Run

cd ./2_selection_tests/Snakemake/
./run.Snakemake.sh --sp <species> --sample <sample> --baseDir <path_to_RegEvol_folder> [options]

Required arguments:

• --sp → species name (must match those in config/params.sh)
• --sample → sample name (used for result subfolders)

Optional arguments (can be passed in any order using --option value):

• --threads → number of threads to allocate (default: 1)
• --system → execution mode: local or SLURM (default: local)
• --baseDir → path to base directory (default: three levels up from Snakefile)
• --extend_config → add any RegEvol-specific parameters defined in /2_selection_tests/Snakemake/config/ (default: '')
• --dryRun → add to run Snakemake in dry-run mode
• --unlock → add to run Snakemake with the --unlock argument
• --help → display usage and exit

Example:

cd ./2_selection_tests/Snakemake/
./run.snakemake.sh --sp human --sample Wilson --threads 10 --dryRun --system SLURM --extend_config TF_source='bed'