____ ____ ____ ____ _ ____
U| _"\ u U /"___|u U /"___| U /"___|u U /"\ u U| _"\ u
\| |_) |/ \| | _ / \| | u \| | _ / \/ _ \/ \| |_) |/
| __/ | |_| | | |/__ | |_| | / ___ \ | __/
|_| \____| \____| \____| /_/ \_\ |_|
||>>_ _)(|_ _// \\ _)(|_ \\ >> ||>>_
(__)__) (__)__) (__)(__) (__)__) (__) (__) (__)__)
PGCGAP is a pipeline for prokaryotic comparative genomics analysis. It can take the pair-end reads, Oxford reads or PacBio reads as input. In addition to genome assembly, gene prediction and annotation, it can also get common comparative genomics analysis results such as phylogenetic trees of single-core proteins and core SNPs, pan-genome, whole-genome Average Nucleotide Identity (ANI), orthogroups and orthologs, COG annotations, substitutions (snps) and insertions/deletions (indels) with only one line of commands.
The software was tested successfully on Windows WSL, Linux x64 platform and macOS. Because this software relies on a large number of other softwares, so it is recommended to install with Bioconda. Because PGCGAP relies on both Gubbins and Orthofinder, which are developed in different versions of python, Gubbins must be installed separately. Once Orthofinder was upgraded to python 3, PGCGAP can be installed with only one command.
Step1: Install Gubbins
If the system is installed with python 3 version of miniconda, Gubbins can be installed directly via conda.
$conda install gubbins
If the python 2 version of miniconda is installed on the system, users need to create a new python 3 environment to install Gubbins. And the Gubbins installation directory need to be added to the environment variable.
#Create an environment called gubbins
$conda create -n gubbins
#Activate the gubbins environment
$conda activate gubbins
#Installation of Gubbins
$conda install gubbins
View the installation path of gubbins and then add this path to the environment variable.
$whereis gubbins
Exit the gubbins environment
$conda deactivate
Step2: Install PGCGAP
$conda create -n pgcgap
$conda activate pgcgap
$conda install pgcgap
$conda deactivate
Step3: Setup COG database (Users should execute this after first installation of pgcgap)
$conda activate pgcgap
$pgcgap --setup-COGdb
$conda deactivate
$pgcgap --help
$pgcgap [functions] [options]
$pgcgap --setup-COGdb
Functions:
[–All] Perform Assemble, Annotate, CoreTree, Pan, OrthoF, ANI and pCOG functions with one command
[–Assemble] Assemble reads into contigs
[–Annotate] Genome annotation
[–CoreTree] Construct single-core proteins tree and SNPs tree of single core genes
[–Pan] Run "roary" pan genome pipeline with gff3 files, and construct a phylogenetic tree with the sing-copy core proteins called by roary
[–OrthoF] Identify orthologous protein sequence families with "OrthoFinder"
[–ANI] Compute whole-genome Average Nucleotide Identity ( ANI )
[–MASH] Genome and metagenome similarity estimation using MinHash
[–pCOG] Run COG annotation for each strain (*.faa), and generate a table containing the relative abundance of each flag for all strains
[–VAR] Rapid haploid variant calling and core genome alignment with "Snippy"
Global Options:
[–strain_num (INT)] [Required by "–All", "–CoreTree", "–Pan", "–VAR" and "–COG"] The total number of strains used for analysis, not including the reference genome
[–ReadsPath (PATH)] [Required by "–All", "–Assemble" and "–VAR"] Reads of all strains as file paths ( Default ./Reads/Illumina )
[–scafPath (PATH)] [Required by "–All", "–Annotate" and "–MASH"] Path for contigs/scaffolds (Default "Results/Assembles/Scaf/Illumina")
[–AAsPath (PATH)] [Required by "–All", "–CoreTree", "–Pan", "–OrthoF" and "–pCOG"] Amino acids of all strains as fasta file paths, ( Default "./Results/Annotations/AAs" )
[–reads1 (STRING)] [Required by "–All", "–Assemble" and "–VAR"] The suffix name of reads 1 ( for example: if the name of reads 1 is "YBT-1520_L1_I050.R1.clean.fastq.gz", "YBT-1520" is the strain same, so the suffix name should be ".R1.clean.fastq.gz")
[–reads2 (STRING)] [Required by "–All", "–Assemble" and "–VAR"] The suffix name of reads 2( for example: if the name of reads 2 is "YBT-1520_2.fq", the suffix name should be "_2.fq" )
[–Scaf_suffix (STRING)] [Required by "–All", "–Annotate" "MASH" and "–ANI"] The suffix of scaffolds or genomes. Here, "-8.fa" for Illumina data, ".contigs.fasta" for PacBio data and Oxford data. Users can also fill in other suffixes according to the actual situation (Default -8.fa)
[–codon (INT)] [Required by "–All", "–Annotate", "–CoreTree" and "–Pan"] Translation table ( Default 11 )
[–suffix_len (INT)] [Required by "–All", "–Assemble" and "–VAR"] (Strongly recommended) The suffix length of the reads, that is the length of your reads name minus the length of your strain name. For example the –suffix_len of "YBT-1520_L1_I050.R1.clean.fastq.gz" is 26 ( "YBT-1520" is the strain name ) ( Default 0 )
[–logs (STRING)] Name of the log file ( Default Logs.txt )
[–threads (INT)] Number of threads to be used ( Default 4 )
Local Options:
–Assemble
[–platform (STRING)] Sequencing Platform, "illumina", "pacbio" and "oxford" available (Default illumina)
[–kmmer (INT)] [Required] k-mer size for genome assembly of Illumina data ( Default 81 )
[–genomeSize (FLOAT)] [Required] An estimate of the size of the genome. Common suffices are allowed, for example, 3.7m or 2.8g. Needed by PacBio data and Oxford data (Default Unset)
–Annotate
[–genus (STRING)] Genus name of your strain ( Default "NA" )
[–species (STRING)] Species name of your strain ( Default "NA")
–CoreTree
[–CDsPath (PATH)] [Required] CDs of all strains as fasta file paths ( Default "./Results/Annotations/CDs" ), if set to "NO", the SNPs of single-copy core genes will not be called
[-c (FLOAT)] Sequence identity threshold, ( Default 0.5)
[-n (INT)] Word_length, -n 2 for thresholds 0.4-0.5, -n 3 for thresholds 0.5-0.6, -n 4 for thresholds 0.6-0.7, -n 5 for thresholds 0.7-1.0 ( Default 2 )
[-G (INT)] Use global (set to 1) or local (set to 0) sequence identity, ( Default 0 )
[-t (INT)] Tolerance for redundance ( Default 0 )
[-aL (FLOAT)] Alignment coverage for the longer sequence. If set to 0.9, the alignment must covers 90% of the sequence ( Default 0.5 )
[-aS (FLOAT)] Alignment coverage for the shorter sequence. If set to 0.9, the alignment must covers 90% of the sequence ( Default 0.7 )
[-g (INT)] If set to 0, a sequence is clustered to the first cluster that meet the threshold (fast cluster). If set to 1, the program will cluster it into the most similar cluster that meet the threshold (accurate but slow mode, Default 1)
[-d (INT)] length of description in .clstr file. if set to 0, it takes the fasta defline and stops at first space ( Default 0 )
–Pan
–OrthoF
–ANI
[–queryL (FILE)] [Required] The file containing paths to query genomes, one per line ( Default scaf.list )
[–refL (FILE)] [Required] The file containing paths to reference genomes, one per line. ( Default scaf.list )
[–ANIO (FILE)] The name of output file ( Default "Results/ANI/ANIs" )
–VAR
[–refgbk (FILE)] [Required] The full path and name of reference genome in GENBANK format ( recommended ), fasta format is also OK. For example: "/mnt/g/test/ref.gbk"
[–qualtype (STRING)] [Required] Type of quality values (solexa (CASAVA < 1.3), illumina (CASAVA 1.3 to 1.7), sanger (which is CASAVA >= 1.8)). ( Default sanger )
[–qual (INT)] Threshold for trimming based on average quality in a window. ( Default 20 )
[–length (INT)] Threshold to keep a read based on length after trimming. ( Default 20 )
[–mincov (INT)] The minimum number of reads covering a site to be considered ( Default 10 )
[–minfrac (FLOAT)] The minimum proportion of those reads which must differ from the reference ( Default 0.9 )
[–minqual (INT)] The minimum VCF variant call "quality" ( Default 100 )
[–ram (INT)] Try and keep RAM under this many GB ( Default 8 )
[–tree_builder (STRING)] Application to use for tree building [raxml|fasttree|hybrid] ( Default fasttree)
[–iterations (INT)] Maximum No. of iterations for gubbins ( Default 5 )
Paths of external programs
Not needed if they were in the environment variables path. Users can check with the "–check-external-programs" option for the essential programs.
[–abyss-bin (PATH)] Path to abyss binary file. Default tries if abyss is in PATH;
[–canu-bin (PATH)] Path to canu binary file. Default tries if canu is in PATH;
[–prodigal-bin (PATH)] Path to prodigal binary file. Default tries if prodigal is in PATH;
[–prokka-bin (PATH)] Path to prokka binary file. Default tries if prokka is in PATH;
[–cd-hit-bin (PATH)] Path to cd-hit binary file. Default tries if cd-hit is in PATH;
[–mafft-bin (PATH)] Path to mafft binary file. Default tries if mafft is in PATH;
[–fasttree-bin (PATH)] Path to the fasttree binary file. Default tries if fasttree is in PATH;
[–pal2nal-bin (PATH)] Path to the pal2nal.pl binary file. Default tries if pal2nal.pl is in PATH;
[–snp-sites-bin (PATH)] Path to the snp-sites binary file. Default tries if snp-sites is in PATH;
[–roary-bin (PATH)] Path to the roary binary file. Default tries if roary is in PATH;
[–orthofinder-bin (PATH)] Path to the orthofinder binary file. Default tries if orthofinder is in PATH;
[–fastANI-bin (PATH)] Path to the fastANI binary file. Default tries if fastANI is in PATH;
[–gubbins-bin (PATH)] Path to the run_gubbins.py binary file. Default tries if run_gubbins.py is in PATH;
[–snippy-bin (PATH)] Path to the snippy binary file. Default tries if snippy is in PATH;
[–sickle-bin (PATH)] Path to the sickle-trim binary file. Default tries if sickle is in PATH.
[–mash-bin (PATH)] Path to the mash binary file. Default tries if mash is in PATH.
Setup COG database
Check the required external programs (It is strongly recommended that this step be performed after the installation of pcgp):
$pgcgap --check-external-programs
Example 1: Perform all functions, take the Escherichia coli as an example, total 6 strains for analysis.
Notice: For the sake of flexibility, The "VAR" function needs to be added additionally.
$pgcgap --All --platform illumina --ReadsPath Reads/Illumina --reads1 _1.fastq.gz --reads2 _2.fastq.gz --suffix_len 11 --kmmer 81 --genus Escherichia --species "Escherichia coli" --codon 11 --strain_num 6 --threads 4 --VAR --refgbk /mnt/h/PGCGAP_Examples/Reads/MG1655.gbff --qualtype sanger
Example 2: Genome assembly.
Illumina reads assembly
In this dataset, the naming format of the genome is "strain_1.fastq.gz" and "strain_2.fastq.gz". The string after the strain name is "_1.fastq.gz", and its length is 11, so "–suffix_len" is set to 11.
pgcgap --Assemble --platform illumina --ReadsPath Reads/Illumina --reads1 _1.fastq.gz --reads2 _2.fastq.gz --kmmer 81 --threads 4 --suffix_len 11
Oxford reads assembly
Oxford nanopore only produces one reads file, so only the parameter of "–reads1" needs to be set, where the value is ".fasta". "–genomeSize" is the estimated genome size, and users can check the genome size of similar strains in NCBI database for reference. The parameter was set to "4.8m" here. The suffix of the reads file here is ".fasta" and its length is 6, so "–suffix_len" is set to 6.
$pgcgap --Assemble --platform oxford --ReadsPath Reads/Oxford --reads1 .fasta --genomeSize 4.8m --threads 4 --suffix_len 6
PacBio also produces only one reads file "pacbio.fastq", the parameter settings are similar to Oxford. The strain name is "pacbio" with the suffix ".fastq" and the suffix length is 6, so "–suffix_len" is set to 6.
$pgcgap --Assemble --platform pacbio --ReadsPath Reads/PacBio --reads1 .fastq --genomeSize 4.8m --threads 4 --suffix_len 6
Example 3: Constructing single-copy core protein tree and core SNPs tree
$pgcgap --CoreTree --CDsPath Results/Annotations/CDs --AAsPath Results/Annotations/AAs --codon 11 --strain_num 6 --threads 4
Example 4: Constructing single-copy core protein tree only.
$pgcgap --CoreTree --CDsPath NO --AAsPath Results/Annotations/AAs --codon 11 --strain_num 6 --threads 4
Example 5: Conduct pan-genome analysis and construct a phylogenetic tree of single-copy core proteins called by roary.
$pgcgap --Pan --codon 11 --strain_num 6 --threads 4 --GffPath Results/Annotations/GFF --PanTree --AAsPath Results/Annotations/AAs
Example 6: Inference of orthologous gene groups.
$pgcgap --OrthoF --threads 4 --AAsPath Results/Annotations/AAs
Example 7: Compute whole-genome Average Nucleotide Identity (ANI).
$pgcgap --ANI --threads 4 --queryL scaf.list --refL scaf.list --ANIO Results/ANI/ANIs --Scaf_suffix .fa
Example 8: Genome and metagenome similarity estimation using MinHash
$pgcgap --MASH --scafPath <PATH> --Scaf_suffix <STRING>
Example 9: Run COG annotation for each strain.
$pgcgap --pCOG --threads 4 --strain_num 6 --AAsPath Results/Annotations/AAs
Example 10: Variants calling and phylogenetic tree construction based on reference genome.
$pgcgap --VAR --threads 4 --refgbk /mnt/h/PGCGAP_Examples/Reads/MG1655.gbff --ReadsPath Reads/Illumina --reads1 _1.fastq.gz --reads2 _2.fastq.gz --suffix_len 11 --strain_num 6 --qualtype sanger --PanTree
The directory where the PGCGAP software runs.
Pair-end reads of all strains in a directory or PacBio reads or Oxford nanopore reads (Default: ./Reads/Illumina/ under the working directory).
Genomes files (complete or draft) in a directory (Default: Results/Assembles/Scaf/Illumina under the working directory).
QUERY_LIST and REFERENCE_LIST files containing full paths to genomes, one per line (default: scaf.list under the working directory). If the "–Assemble" function was run first, the list file will be generated automatically.
Amino acids file (With ".faa" as the suffix) and nucleotide (With ".ffn" as the suffix) file of each strain placed into two directories (default: "./Results/Annotations/AAs/" and "./Results/Annotations/CDs/"). The ".faa" and ".ffn" files of same strain should have the same prefix name. The name of protein IDs and gene IDs should be started with the strain name. The "Prokka" software was suggested to generate the input files. If the "–Annotate" function was run first, the files will be generated automatically. If the "–CDsPath" was set to "NO", the nucleotide files will not be needed.
Genomes files (complete or draft) in a directory (Default: Results/Assembles/Scaf/Illumina under the working directory).
A set of protein sequence files (one per species) in FASTA format under a directory (default: "./Results/Annotations/AAs/"). If the "–Annotate" function was run first, the files will be generated automatically.
GFF3 files (With ".gff" as the suffix) of each strain placed into a directory. They must contain the nucleotide sequence at the end of the file. All GFF3 files created by Prokka are valid (default: ./Results/Annotations/GFF/). protein sequence files (one per species) in FASTA format under another directory were also needed (default: "./Results/Annotations/AAs/"). If the "–Annotate" function was run first, the files will be generated automatically.
Amino acids file (With ".faa" as the suffix) of each strain placed into a directory (default: ./Results/Annotations/AAs/). If the "–Annotate" function was run first, the files will be generated automatically.
Pair-end reads of all strains in a directory (default: ./Reads/Over/ under the working directory).
The full path of reference genome in fasta format or GenBank format (must be provided).
Results/Assembles/Illumina/*_assembly
Directories contain Illumina assembly files and information of each strain.
Results/Assembles/PacBio
Directories contain PacBio assembly files and information of each strain.
Results/Assembles/Illumina/Oxford
Directories contain Oxford nanopore assembly files and information of each strain.
Results/Assembles/Scaf/Illumina
Directory contains Illumina contigs/scaffolds of all strains.
Results/Assembles/Scaf/Oxford
Directory contains Oxford nanopore contigs/scaffolds of all strains.
Results/Assembles/Scaf/PacBio
Directory contains PacBio contigs/scaffolds of all strains.
Results/Annotations/*_annotation
directories contain annotation files of each strain.
Results/Annotations/AAs
Directory contain amino acids sequences of all strains.
Results/Annotations/CDs
Directory contain nucleotide sequences of all strains.
Results/Annotations/GFF
Directory contain the master annotation of all strains in GFF3 format.
Results/ANI/ANIs
The file contains comparation information of genome pairs. The document is composed of five columns, each of which represents query genome, reference genome, ANI value, count of bidirectional fragment mappings, total query fragments.
Results/ANI/ANIs.matrix
file with identity values arranged in a phylip-formatted lower triangular matrix.
Results/ANI/ANIs.heatmap
An ANI matrix of all strains.
Results/ANI/ANI_matrix.pdf
The heatmap plot of "ANIs.heatmap".
Results/CoreTrees/faa/ALL.core.protein.fasta
Concatenated and aligned sequences file of single-copy core proteins.
Results/CoreTrees/faa2ffn/ALL.core.nucl.fasta
Concatenated and aligned sequences file of single-copy core genes.
Results/CoreTrees/faa2ffn/ALL.core.snp.fasta
Core SNPs of single-copy core genes in fasta format.
Results/CoreTrees/ALL.core.protein.nwk
The phylogenetic tree file of single-copy proteins for all strains.
Results/CoreTrees/faa2ffn/ALL.core.snp.nwk
The phylogenetic tree file of SNPs of single-copy core genes for all strains.
Results/CoreTrees/"Other_files"
Intermediate directories and files.
Results/MASH/MASH
Pairwise distance between pair genomes, each column represents Reference-ID, Query-ID, Mash-distance, P-value, and Matching-hashes, respectively.
Results/MASH/MASH2
Pairwise similarity between pair genomes, each column represents Reference-ID, Query-ID, similarity, P-value, and Matching-hashes, respectively.
Results/MASH/MASH.heatmap
A similarity matrix of all genomes.
Results/MASH/MASH.heatmap
A heat map plot of "MASH_matrix.pdf".
Results/PanGenome/Pangenome_Pie.pdf
A 3D pie chart and a fan chart of the breakdown of genes and the number of isolate they are present in.
Results/PanGenome/pangenome_frequency.pdf
A graph with the frequency of genes versus the number of genomes.
Results/PanGenome/Pangenome_matrix.pdf
A figure showing the tree compared to a matrix with the presence and absence of core and accessory genes.
Results/PanGenome/Core/Roary.core.protein.nwk
A phylogenetic tree based on single-copy core proteins called by roary software.
Results/PanGenome/Other_files
see roary outputs.
*.COG.xml, *.2gi.table, *.2id.table, *.2Sid.table
Intermediate files.
*.2Scog.table
The super COG table of each strain.
*.2Scog.table.pdf
A plot of super COG table in pdf format.
All_flags_relative_abundances.table A table containing the relative abundance of each flag for all strains.
Results/Variants/directory-named-in-strains
directories containing substitutions (snps) and insertions/deletions (indels) of each strain. See Snippy outputs for detail.
Results/Variants/Core
The directory containing Core SNP phylogeny files.
PGCGAP is free software, licensed under GPLv3.
Please report any issues to the issues page or email us at liaochenlanruo@webmail.hzau.edu.cn.
If you use this software please cite: (Please keep an eye on it as it will be noted soon!)
Check the log file named in "strain_name.log" under Results/Variants/
$conda install java-jdk=8.0.112
Click here for more solutions.
When running prokka of Assemble founction, this error could happened, the error message shows as following:
Error: A JNI error has occurred, please check your installation and try again
Exception in thread "main" java.lang.UnsupportedClassVersionError: minced has been compiled by a more recent version of the Java Runtime (class file version 55.0), this version of the Java Runtime only recognizes class file versions up to 52.0
at java.lang.ClassLoader.defineClass1(Native Method)
at java.lang.ClassLoader.defineClass(ClassLoader.java:763)
at java.security.SecureClassLoader.defineClass(SecureClassLoader.java:142)
at java.net.URLClassLoader.defineClass(URLClassLoader.java:468)
at java.net.URLClassLoader.access$100(URLClassLoader.java:74)
at java.net.URLClassLoader$1.run(URLClassLoader.java:369)
at java.net.URLClassLoader$1.run(URLClassLoader.java:363)
at java.security.AccessController.doPrivileged(Native Method)
at java.net.URLClassLoader.findClass(URLClassLoader.java:362)
at java.lang.ClassLoader.loadClass(ClassLoader.java:424)
at sun.misc.Launcher$AppClassLoader.loadClass(Launcher.java:349)
at java.lang.ClassLoader.loadClass(ClassLoader.java:357)
at sun.launcher.LauncherHelper.checkAndLoadMain(LauncherHelper.java:495)
[01:09:40] Could not determine version of minced - please install version 2.0 or higher
Users can downgrade the minced to version 0.3 to solve this problem.
$conda install minced=0.3
Click here for detail informations.
This error may happen when running function "VAR" on macOS. It is an error of openssl. Users can solve this problem as following:
#Firstly, install brew if have not installed before
$ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
#Install openssl with brew
$brew install openssl
#Create the soft link for libraries
$ln -s /usr/local/opt/openssl/lib/libcrypto.1.0.0.dylib /usr/local/lib/
$ln -s /usr/local/opt/openssl/lib/libssl.1.0.0.dylib /usr/local/lib/
Click here for more informations
This warning may happen when running function "Pan". It is a warning of Roary software. The content of line 61 is "require Encode::ConfigLocal;". Users can ignore the warning. Click here for details.
V1.0.3
V1.0.4
V1.0.5
V1.0.6
V1.0.7
V1.0.8
V1.0.9