PyamilySeq

PyamilySeq is a Python tool for clustering gene sequences into groups based on sequence similarity identified by tools such as CD-HIT, BLAST, DIAMOND or MMseqs2. This work is an extension of the gene family / pangenome tool developed for the StORF-Reporter publication in NAR (https://doi.org/10.1093/nar/gkad814).

Features

• End-to-End: PyamilySeq can take a directory of GFF+FASTA files, run CD-HIT for clustering and process the results.
• Clustering input: Supports input from CD-HIT formatted files as well as CSV and TSV node-edge lists (MMseqs2 and -outfmt 6 from BLAST/DIAMOND).
• Reclustering: Allows for the addition of new sequences post-initial clustering - Ensures continuity of contemporary clustering results and highlights impact of novel gene predictions.
• 'Genus Mode': Unlike other 'pangenome' tools, PyamilySeq can identify gene groups found across multiple genera as unique entities (see below).
• Output: Generates a 'Roary/Panaroo' formatted presence-absence CSV formatted file for downstream analysis.
  • User-define species-/genus-wide gene groups - User has control over grouping parameters (core = 99/95% or min 6 genera etc).
  • Aligns representative sequences using MAFFT.
  • Output concatenated aligned sequences for tree building.
  • Optionally output sequences of each separate identified gene group.
  • Group-Splitter tool to split multi-copy gene groups.
  • Numerous additional tools to assist in the pre- and post-processing of data.

Installation

PyamilySeq probably requires Python 3.6 or higher and the levenshtein library (https://pypi.org/project/Levenshtein/) -
If levenshtein is not available, a Python implementation is utilised which is significantly slower.

Install using pip:

pip install PyamilySeq [optionally add -U]

PyamilySeq is currently still under active development so expect 'regular' updates with bugfixes and new features.
To update to the newest version add '-U' to end of the pip install command.

Example usage: Below are two examples of running PyamilySeq in its two main modes.

usage: PyamilySeq.py [-h] {Full,Partial} ...

PyamilySeq v1.3.4: A tool for gene clustering and analysis.

positional arguments:
  {Full,Partial}  Choose a mode: 'Full' or 'Partial'.
    Full          Full mode: PyamilySeq to cluster with CD-HIT and process output.
    Partial       Partial mode: PyamilySeq to process pre-clustered data.

options:
  -h, --help      show this help message and exit

'Full Mode': Will conduct clustering of sequences with CD-HIT as part of PyamilySeq run

PyamilySeq Full -output_dir .../PyamilySeq_10_AA_90_80_Full_GFFs -input_type combined -input_dir .../genomes/ -name_split_gff _combined.gff3

'Partial Mode': Will process the output of a sequence clustering from MMseqs, BLAST, DIAMOND etc.

PyamilySeq Partial -clustering_format CD-HIT -cluster_file .../all_10_combined_pep_CD-HIT_90_80.clstr  -original_fasta .../all_10_combined_pep.fasta -output_dir .../PyamilySeq_10_AA_90_80_Partial -write_groups 99 -align 

Note: using a '-clustering_format' other than the default CD-HIT, requires input to be two in two columns as below (Same format as MMseqs2 tsv and BLAST outfmt 6) - Genome name and sequence name are separated by '|'.

Escherichia_coli_110957|ENSB_lL-zIKt-gh0oSno	Escherichia_coli_110957|ENSB_lL-zIKt-gh0oSno
Escherichia_coli_110957|ENSB_lL-zIKt-gh0oSno	Escherichia_coli_113290|ENSB_2fj4rJ8e8Z9PNdX
Escherichia_coli_110957|ENSB_lL-zIKt-gh0oSno	Escherichia_coli_b185|ENSB_G_PVe28-ej8q-3S
Escherichia_coli_110957|ENSB_TIZS9kbTvShDvyX	Escherichia_coli_110957|ENSB_TIZS9kbTvShDvyX

Example output:

Running PyamilySeq v1.3.4
Calculating Groups
Number of Genomes: 10
Gene Groups
First_core_99: 2994
First_core_95: 0
First_core_15: 3266
First_core_0: 5466
Total Number of First Gene Groups (Including Singletons): 11726
Outputting gene_presence_absence file
Outputting gene group FASTA files
Combined FASTA file saved to: ../combined_group_sequences_dna.fasta
Processing gene group alignment
Thank you for using PyamilySeq -- A detailed user manual can be found at https://github.com/NickJD/PyamilySeq
Please report any issues to: https://github.com/NickJD/PyamilySeq/issues

Reclustering:

Reclustering can be used to see where additional sequences/genes lay in relation to a contemporary pangenome/gene grouping.

PyamilySeq Partial -clustering_format CD-HIT -cluster_file .../all_10_combined_pep_CD-HIT_90_80.clstr -reclustered .../all_10_combined_pep_CD-HIT_90_80_AND_StORFs_CD-HIT_90_80.clstr -original_fasta .../all_10_combined_pep_AND_StORFs.fasta -output_dir .../PyamilySeq_10_AA_90_80_Partial_Reclustered_StORFs -write_groups 99 -align 

As can be seen below, the additional sequences recovered by the StORF-Reporter annotation tool have 'extended' contemporary or created entirely new gene groups. 'First' corresponds to the groups identified from the first clustering round and 'Second' for the second. In 'reclustering' mode, First_core_# groups are unaffected thus retaining the initial grouping information.

Number of Genomes: 10
Gene Groups
First_core_99: 2994
First_core_95: 0
First_core_15: 3266
First_core_0: 5466
extended_core_99: 3
extended_core_95: 0
extended_core_15: 49
extended_core_0: 0
combined_core_99: 0
combined_core_95: 0
combined_core_15: 3
combined_core_0: 0
Second_core_99: 0
Second_core_95: 0
Second_core_15: 20
Second_core_0: 39
only_Second_core_99: 768
only_Second_core_95: 0
only_Second_core_15: 4472
only_Second_core_0: 8395
Total Number of First Gene Groups (Including Singletons): 11726
Total Number of Second Gene Groups (Including Singletons): 25359
Total Number of First Gene Groups That Had Additional Second Sequences But Not New Genomes: 5

PyamilySeq is separated into two main 'run modes', Full and Partial. They each have their own set of required and optional arguments.

PyamilySeq - Full Menu:

usage: PyamilySeq.py Full [-h] -output_dir OUTPUT_DIR -input_type {separate,combined,fasta} [-input_dir INPUT_DIR] [-input_fasta INPUT_FASTA] [-name_split_gff NAME_SPLIT_GFF] [-name_split_fasta NAME_SPLIT_FASTA] [-seq_type {AA,DNA}] [-gene_ident GENE_IDENT] [-c PIDENT] [-s LEN_DIFF] [-fast_mode]
                          [-group_mode {Species,Genus}] [-species_groups SPECIES_GROUPS] [-genus_groups GENUS_GROUPS] [-write_groups WRITE_GROUPS] [-write_individual_groups] [-align] [-align_aa] [-no_gpa] [-M MEM] [-T THREADS] [-verbose] [-v]

options:
  -h, --help            show this help message and exit
  -output_dir OUTPUT_DIR
                        Directory for all output files.
  -input_type {separate,combined,fasta}
                        Type of input files: 'separate' for matching FASTA and GFF files, 'combined' for GFF+FASTA, or 'fasta' for a prepared FASTA file.
  -input_dir INPUT_DIR  Directory containing GFF/FASTA files - Use with -input_type separate/combined.
  -input_fasta INPUT_FASTA
                        Input FASTA file - Use with - input_type fasta.
  -name_split_gff NAME_SPLIT_GFF
                        Substring to split filenames and extract genome names for gff files (e.g., '_combined.gff3') - Use with -input_type separate/combined.
  -name_split_fasta NAME_SPLIT_FASTA
                        Substring to split filenames and extract genome names for fasta files if named differently to paired gff files (e.g., '_dna.fasta') - Use with -input_type separate/combined.
  -seq_type {AA,DNA}
                        Clustering mode: 'DNA' or 'AA'.
  -gene_ident GENE_IDENT
                        Gene identifiers to extract sequences (e.g., 'CDS, tRNA').
  -c PIDENT             Sequence identity threshold for clustering (default: 0.90) - CD-HIT parameter '-c'.
  -s LEN_DIFF           Length difference threshold for clustering (default: 0.80) - CD-HIT parameter '-s'.
  -fast_mode            Enable fast mode for CD-HIT (not recommended) - CD-HIT parameter '-g'.
  -group_mode {Species,Genus}
                        Grouping mode: 'Species' or 'Genus'.
  -species_groups SPECIES_GROUPS
                        Gene groupings for 'Species' mode (default: '99,95,15').
  -genus_groups GENUS_GROUPS
                        Gene groupings for 'Genus' mode (default: '1-10').
  -write_groups WRITE_GROUPS
                        Output gene groups as a single FASTA file (specify levels: e.g., '-w 99,95'). - triggers '-wig'.
  -write_individual_groups
                        Output individual FASTA files for each group.
  -align                Align and concatenate sequences for 'core' groups specified with '-w'.
  -align_aa             Align sequences as amino acids.
  -no_gpa               Skip creation of gene_presence_absence.csv.
  -M MEM                Memory allocation for clustering (MB) - CD-HIT parameter '-M'.
  -T THREADS            Number of threads for clustering/alignment - CD-HIT parameter '-T' | MAFFT parameter '--thread'.
  -verbose              Print verbose output.
  -v, --version         Print version number and exit.

PyamilySeq - Partial Menu:

usage: PyamilySeq.py Partial [-h] -clustering_format {CD-HIT,MMseqs,BLAST} -cluster_file CLUSTER_FILE -original_fasta ORIGINAL_FASTA -output_dir OUTPUT_DIR [-reclustered RECLUSTERED] [-seq_tag SEQUENCE_TAG] [-group_mode {Species,Genus}] [-species_groups SPECIES_GROUPS] [-genus_groups GENUS_GROUPS]
                             [-write_groups WRITE_GROUPS] [-write_individual_groups] [-align] [-align_aa] [-no_gpa] [-M MEM] [-T THREADS] [-verbose] [-v]

options:
  -h, --help            show this help message and exit
  -clustering_format {CD-HIT,MMseqs,BLAST}
                        Clustering format used: CD-HIT, MMseqs2, or BLAST.
  -cluster_file CLUSTER_FILE
                        Cluster file containing pre-clustered groups from CD-HIT, MMseqs, BLAST etc.
  -original_fasta ORIGINAL_FASTA
                        FASTA file