Gene EXpansion and CONtraction analysis pipeline

# EXCON (v2.3.1) A Nextflow pipeline for gene family **EX**pansion and **CON**traction analysis across multiple species using CAFE5. Given a set of genome assemblies and annotations, EXCON builds orthogroups with OrthoFinder, fits and compares multiple CAFE models to identify gene families evolving at significantly different rates, and automatically selects the best-fitting model for downstream analysis. Optionally, GO enrichment analysis can be run on expanded and contracted gene families, and on genes grouped by chromosome. It works with any set of species that have a genome (fasta) and annotation (gff) file. (minimum of 5 species ideally up to around 30). Maximum 100 species (normally). ## Overview The general pipeline logic is as follows: * Downloads genome and annotation files from NCBI `[NCBIGENOMEDOWNLOAD]`, or you provide your own. * Unzips the files, if necessary `[GUNZIP]` * Sanitises GFF annotations and extracts longest isoform `[AGAT_SPKEEPLONGESTISOFORM]`. * Gets the protein sequences `[GFFREAD]`. * Renames the genes to gene name (as some will be isoform name) `RENAME_FASTA`. * Finds orthologous genes across species `[ORTHOFINDER_CAFE]`, or accepts a pre-computed tree and orthogroups to skip this step (see `--input_tree` / `--input_orthogroups`). * Rescales OrthoFinder branch lengths and converts to an ultrametric tree for CAFE `[RESCALE_TREE]`, `[CAFE_PREP]`. * Prepares gene count input, estimates the error model, and builds an ultrametric tree `[CAFE_PREP]`. * Runs CAFE5 with k=1 to k=`cafe_max_k` (default 6) rate categories in parallel `[CAFE_RUN_K]`. * Compares all k runs by AIC and selects the best k `[CAFE_SELECT_K]`. * Re-runs CAFE5 at the best k with Poisson birth-death (`-p`) `[CAFE_RUN_BEST]`. * Compares uniform vs Poisson model at the best k by likelihood, selects the winner `[CAFE_MODEL_COMPARE]`. * Plots gene family expansions and contractions for the best model `[CAFE_PLOT]`. ### Optional — GO enrichment (`--run_eggnog` or `--predownloaded_gofiles`) * GO annotation can be run in two ways: - **Run EggNOG-mapper** (`--run_eggnog`): assigns GO terms from scratch using the EggNOG database (~45 GB). Provide `--eggnog_data_dir` to reuse a pre-downloaded copy and avoid re-downloading on every run. - **Supply your own GO files** (`--predownloaded_gofiles`): if you already have gene-to-GO mappings, point to a directory of `{species_id}.go.txt` files and skip EggNOG entirely (see [Parameters](#go-annotation-optional) below). * Optionally downloads the EggNOG-mapper database `[EGGNOG_DOWNLOAD]`. Default is on, unless you provide it. * Optionally assigns GO terms to genes using EggNOG-mapper `[EGGNOGMAPPER]`, or reads GO terms from user-supplied files. * Optionally prepares GO gene lists from the best CAFE model results `[CAFE_GO_PREP]`. * Optionally runs GO enrichment in parallel, one job per species/node and direction (expansion/contraction) `[CAFE_GO_RUN]`. ### Optional — chromosome GO enrichment (`--chromo_go`, requires GO annotation) * Optionally plots GO enrichment of genes by chromosome `[CHROMO_GO]`. * Optionally summarizes GO enrichment by chromosome `[SUMMARIZE_CHROMO_GO]`. ### Optional — genome quality statistics (`--stats`) * Optionally describes genome assembly and annotation: - `[BUSCO_BUSCO]`: Completeness of the genome compared to expected gene set. - `[QUAST]`: Assembly contiguity statistics (N50 etc). - `[AGAT_SPSTATISTICS]`: Gene, exon, and intron statistics. ## Installation Nextflow pipelines require a few prerequisites. There is further documentation on the nf-core webpage [here](https://nf-co.re/docs/usage/installation), about how to install Nextflow. ### Prerequisites - [Docker](https://docs.docker.com/engine/install/) or [Singularity](https://docs.sylabs.io/guides/3.11/admin-guide/installation.html). - [Java](https://www.java.com/en/download/help/download_options.html) and [openJDK](https://openjdk.org/install/) >= 8 (**Please Note:** When installing Java versions are `1.VERSION` so `Java 8` is `Java 1.8`). - [Nextflow](https://www.nextflow.io/docs/latest/getstarted.html) >= `v25.10.0`. - When running nextflow with this pipeline, ideally run `NXF_VER=25.10.0` beforehand, to ensure functionality on this version. ### Install To install the pipeline please use the following commands but replace VERSION with a [release](https://github.com/Eco-Flow/excon/releases). `wget https://github.com/Eco-Flow/excon/archive/refs/tags/VERSION.tar.gz -O - | tar -xvf -` or `curl -L https://github.com/Eco-Flow/excon/archive/refs/tags/VERSION.tar.gz --output - | tar -xvf -` This will produce a directory in the current directory called `excon-VERSION` which contains the pipeline. ## Inputs ### Required * `--input /path/to/csv/file` - A singular csv file as input in one of the two formats stated below. This csv can take 2 forms: * A 2 field csv where each row is a unique species name followed by a Refseq genome reference ID (**NOT** a Genbank reference ID) i.e. `data/input_small-s3.csv`. The pipeline will download the relevant genome fasta file and annotation gff3 (or gff augustus) file. * A 3 field csv where each row is a unique species name, followed by an absolute path to a genome fasta file, followed by an absolute path to an annotation gff3 (or gff augustus) file. Input can be gzipped (.gz) or not. **Please Note:** The genome has to be chromosome level not contig level. 2 fields (Name,Refseq_ID): ``` Drosophila_yakuba,GCF_016746365.2 Drosophila_simulans,GCF_016746395.2 Drosophila_santomea,GCF_016746245.2 ``` 3 fields (Name,genome.fna,annotation.gff): ``` Drosophila_yakuba,data/Drosophila_yakuba/genome.fna.gz,data/Drosophila_yakuba/genomic.gff.gz Drosophila_simulans,data/Drosophila_simulans/genome.fna.gz,data/Drosophila_simulans/genomic.gff.gz Drosophila_santomea,data/Drosophila_santomea/genome.fna.gz,data/Drosophila_santomea/genomic.gff.gz ``` > **Note:** Genomes should be chromosome-level, not contig-level. RefSeq IDs must be used (not GenBank IDs). ## Parameters ### Core options | Parameter | Description | Default | |-----------|-------------|---------| | `--input` | Path to input CSV file | Required unless using `--input_tree` + `--input_orthogroups` without `--run_eggnog` or `--stats` | | `--outdir` | Output directory | `results` | | `--groups` | NCBI taxonomy group for genome download (e.g. `insects`, `bacteria`) | `insects` | | `--help` | Display help message | `false` | | `--custom_config` | Path to a custom Nextflow config file | `null` | ### Quality statistics (optional) | Parameter | Description | Default | |-----------|-------------|---------| | `--stats` | Run BUSCO, QUAST and AGAT statistics on genomes | `null` | | `--busco_lineage` | BUSCO lineage database (e.g. `insecta_odb10`) | `null` | | `--busco_mode` | BUSCO mode (`genome`, `proteins`, `transcriptome`) | `null` | | `--busco_lineages_path` | Path to local BUSCO lineage databases | `null` | | `--busco_config` | Path to BUSCO config file | `null` | ### OrthoFinder options (optional) | Parameter | Description | Default | |-----------|-------------|---------| | `--orthofinder_v2` | Use OrthoFinder v2.5.5 instead of v3.1.4. v2 produces Hierarchical Orthogroups (`N0.tsv`) which have lower copy-number variance and are better suited to CAFE5. v3 produces flat orthogroups (`Orthogroups.tsv`). Recommended for large datasets (>30 species) or when CAFE5 fails to converge. | `false` | | `--orthofinder_method` | Gene tree inference method: `msa` or `dendroblast` | `msa` | | `--orthofinder_search` | Sequence search program: `diamond`, `blast`, or `mmseqs2` | `diamond` | | `--orthofinder_msa_prog` | MSA program (requires `--orthofinder_method msa`): `mafft` or `muscle` | `mafft` | | `--orthofinder_tree` | Tree inference method (requires `--orthofinder_method msa`): `fasttree`, `raxml`, `raxml-ng`, or `iqtree` | `fasttree` | > **Note:** `-A` and `-T` are only valid when `-M msa` is set. If you set `--orthofinder_msa_prog` or `--orthofinder_tree` without `--orthofinder_method msa`, OrthoFinder will error. ### CAFE gene family evolution | Parameter | Description | Default | |-----------|-------------|---------| | `--skip_cafe` | Skip CAFE analysis | `null` | | `--cafe_max_k` | Maximum number of k rate categories to test (runs k=1 through k=N in parallel) | `6` | | `--cafe_max_differential` | Maximum gene count differential for CAFE filtering on retry | `50` | | `--tree_scale_factor` | Factor to multiply all OrthoFinder branch lengths by before `chronos()` converts the tree to a time tree for CAFE5. Lower values can cause numerical issues. | `1000` | | `--input_tree` | Path to a pre-computed rooted species tree (Newick format) — skips OrthoFinder when used with `--input_orthogroups` | `null` | | `--input_orthogroups` | Path to a pre-computed `Orthogroups.tsv` from a previous OrthoFinder run — skips OrthoFinder when used with `--input_tree` | `null` | > **Skipping OrthoFinder:** OrthoFinder is the slowest step in the pipeline. If you have already run it > (the results are in `results/orthofinder_cafe/ortho_cafe/`), you can reuse the outputs. > The orthogroups file to pass depends on which version of OrthoFinder was used: > > **OrthoFinder v2 (`--orthofinder_v2`)** — use `N0.tsv` (Hierarchical Orthogroups): > ``` > --input_tree results/orthofinder_cafe/ortho_cafe/Species_Tree/SpeciesTree_rooted_node_labels.txt \ > --input_orthogroups results/orthofinder_cafe/ortho_cafe/Phylogenetic_Hierarchical_Orthogroups/N0.tsv > ``` > > **OrthoFinder v3 (default)** — use `Orthogroups.tsv` (flat orthogroups): > ``` > --input_tree results/orthofinder_cafe/ortho_cafe/Species_Tree/SpeciesTree_rooted_node_labels.txt \ > --input_orthogroups results/orthofinder_cafe/ortho_cafe/Orthogroups/Orthogroups.tsv > ``` > Both `--input_tree` and `--input_orthogroups` must be supplied together. If either is omitted, OrthoFinder runs normally. > **No