jeffersonfparil/compare_genomes

A comparative genomics workflow using Nextflow, conda, Julia and R

compare_genomes: a comparative genomics workflow to streamline the analysis of evolutionary divergence across eukaryotic genomes

Comparative genomics papers seldomly fully disclose their comparative genomics workflow. This hinders reproducibility and transferability of novel methods. To address this issue, we use Nextflow and Conda to improve the transferability of our lab group's comparative genomics workflow which has been used to assess the reference genome assembly of Lolium rigidum and a new genome assembly of Salvia hispanica.

Plead to community: looking for co-maintainers

Hello fellow scholars, I am seeking a fellow co-maintainer of this project. I have been having less and less time to devote on this open-source project since I left academia, and would like to keep this project thriving. Please contact me through email to discuss further. Cheers!

Installation

1. Download compare_genomes repository

git clone https://github.com/jeffersonfparil/compare_genomes.git

2. Install and initialise conda

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
sh ./Miniconda3-latest-Linux-x86_64.sh

3. Import and activate compare_genomes conda environment

conda env create -n compare_genomes --file compare_genomes/compare_genomes.yml
conda activate compare_genomes

Quickstart

Run the example pipeline:

cd compare_genomes
nano config/params.config # replace dir = '/data/TEST' with a valid path in your machine
nano config/process.config # modify cpus and memory to match the configuration of your machine
chmod +x run.sh
time ./run.sh

Setup parameters

Set up the parameters for the comparative genomics analysis. You will find the following files in compare_genomes/config folder (Note: you will need to include at least 3 species for gene family contraction and expansion analysis to work).

1. urls.txt: web links or absolute paths to the genome sequence, genome annotation, coding DNA sequence, and amino acid sequences for at least 3 species you wish to include in the analyses.

   • Formatted as headerless, two-columned, comma-separated file
   • Column 1: filename of the genome sequence, genome annotation, coding DNA sequence, and amino acid sequences (Note: the species names and extension names should be the consistent across these files, e.g. *.fna for the genomes, *.gff for the annotations, *.cds for the coding DNA sequences, and *.faa for the amino acid sequences)
   • Column 2: URL (uniform resource locator) of the zipped (*.gz or *.zip) or unzipped files for download. Alternatively, this can be the absolute path of pre-downloaded zipped (*.gz or *.zip) or unzipped files.

2. dates.txt: pairwise divergence times between the species you wish to include in the analyses (e.g. look up dvergence times from http://timetree.org/).

   • Formatted as headerless, two-columned, tab-delimited file
   • Column 1: two species separated by a comma with the same names used in the urls.txt
   • Column 2: time in million years, e.g. -160 for 160 million years ago

3. comparisons_4DTv.txt: list of species and pairs of species you wish to include in the estimation of transversion rates among 4-fold degenerate sites (4DTv). This statistic is used as a biological clock, where more mutations at the third codon position means more time has passed between a pair of sequences within and among species (Note: uses 2-copy genes for the estimation of 4DTv).

   • Formatted as headerless, one-columned file
   • Column 1: Species and pairs of species which should be the same names as in urls.txt, and species pairs are written for example as "Zea_mays X Oryza_sativa".

4. venn_species_max_5.txt: list of at most 5 species you wish to include in the Venn diagram comparing the differences and commonalities of orthogroups between species. It is currently not possible to fit more than 5 species in the Venn diagram because of the limitations of the plotting package used.

   • Formatted as headerless, one-columned file
   • Column 1: Species names which should be the same names as in urls.txt.

5. genes.txt: links to the gene sequences you wish to assess expansion/contraction, and non-synonymous/synonymous mutation (Ka/Ks) rates between pairs of sequences within and among species.

   • Formatted as headerless, three-columned, comma-separated file
   • Column 1: non-critical information: phenotype name or some identifying name
   • Column 2: non-critical information: species name which may be any species even ones not included in the analyses
   • Column 3: URL of the genes for download

6. params.config: configuration file listing the variables specific to the analyses you wish to perform.

   • dir: output directory.
   • species_of_interest: the focal species of interest which should be the same as in urls.txt, dates.txt, and genes.txt.
   • species_of_interest_panther_HMM_for_gene_names_url: URL to the specific Panther HMM database to extract gene names from, preferrably from the same species which will be used for gene ontology (GO) term enrichment analysis. See the version 17.0 list at http://data.pantherdb.org/ftp/sequence_classifications/17.0/PANTHER_Sequence_Classification_files (note the http as opposed to https).
   • panther_hmm_database_location: URL or absolute path to the PANTHER HMM database. The default is the URL to version 17 i.e. http://data.pantherdb.org/ftp/panther_library/17.0/PANTHER17.0_hmmscoring.tgz). If you have downloaded the database before then use the absolute path of your pre-downloaded and unzipped database directory.
   • panther_hmm_classifications_location: URL to the PANTHER HMM classifications file containing the GO terms associated with each HMM class. The default is version 7, i.e. http://data.pantherdb.org/ftp/hmm_classifications/17.0/PANTHER17.0_HMM_classifications
   • urls: location of urls.txt.
   • dates: location of dates.txt.
   • comparisons_4DTv: location of comparisons_4DTv.txt.
   • venn_species_max_5: location of venn_species_max_5.txt.
   • genes: location of genes.txt.
   • cafe5_n_gamma_cats: number of the gamma values (parameter of the substitution model) to use for the assessment of significant gene family expansion and contraction using CAFE5. If this is equal to one, then we use the substitution model without the gamma function.
   • cafe5_pvalue: significance threshold of gene family expansion and contraction.
   • go_term_enrich_genome_id: genome ID for the species specified in species_of_interest_panther_HMM_for_gene_names_url or any species you wish to use. Find the appropriate taxon ID from here.
   • go_term_enrich_annotation_id: code for the gene ontology level you wish to use, e.g. GO:0008150 for "Biological Process". See the list of GO codes here.
   • go_term_enrich_test: GO term enrichment test to perform which can be either FISHER (Fisher's Exact Test) or BINOMIAL" (binomial distribution test).
   • go_term_enrich_correction: multiple testing correction which can be NONE, FDR (False discovery rate), or BONFERRONI (Bonferroni correction).
   • go_term_enrich_ngenes_per_test: number of randomly sampled genes to include in each GO term enrichment analysis.
   • go_term_enrich_ntests: number GO term enrichment analyses to perform.

7. process.config: configuration file listing the computing resource allocation available to you. Assign the number of cpus and memory capacity to use for low and high resources intensive tasks:

   • LOW_MEM_LOW_CPU
   • HIGH_MEM_HIGH_CPU

Run

Once the parameters for your specific comparative genomics analysis are set up, you can run the entire workflow in one go, or you may wish to run each module separately.

1. Run the entire workflow (6-hour runtime using 32 cores @ 2GHz)

cd compare_genomes
chmod +x run.sh
time ./run.sh

2. Run each module individually but this needs to be in order as each module depends on the output of the previous one. Run these if you wish to troubleshoot and/or modify the workflow for analyses not implemented here, e.g. add a GO term enrichment analysis for significantly contracted gene families.

cd compare_genomes
nextflow run modules/setup.nf                               -c config/params.config ### SETUP
nextflow run modules/orthofinder.nf                         -c config/params.config ### DETERMINE ORTHOGROUPS
nextflow run modules/single_gene_orthogroups_tree.nf        -c config/params.config ### BUILD A PHYLOGENETIC TREE USING SINGLE-COPY GENE FAMILIES
nextflow run modules/single_gene_orthogroups_4DTv.nf        -c config/params.config ### COMPUTE THE TRANSVERSION RATES AMONNG 4-FOLD DEGENERATE SITES (4DTv)
nextflow run modules/gene_family_contraction_expansion.nf   -c config/params.config ### ASSESS SIGNIFICANT GENE FAMILY CONTRACTION AND EXPANSION (This may take days depending on your machine, number of species, and proteome sizes)
nextflow run modules/GO_enrichment.nf                       -c config/params.config ### GENE ONTOLOGY TERM ENRICHMENT ANALYSIS OF SIGNIFICANTLY EXPANDED GENE FAMILIES
nextflow run modules/assess_WGD.nf                          -c config/params.config ### ASSESS WHOLE GENOME DUPLICATION EVENTS
nextflow run modules/plot_tree_conex_venn_4DTv.nf           -c config/params.config ### PLOT THE PHYLOGENETIC TREE, CONTRACTION/EXPANSION, GENE SETS VENN DIAGRAM AND 4DTv
nextflow run modules/assess_specific_genes.nf               -c config/params.config ### ASSESS CONTRACTION/EXPANSION AND NON-SYNONYMOUS TO SYNONYMOUS NUCLEOTIDE SUBSTITION RATIOS

### Print the most recent error log
cat $(ls -r work/*/*/.command.err | tail -n1)

Workflow and example output plot

The main output file of the workflow is the summary figure in scalable vector graphics (svg) format. This is illustrated in the right panel of the figure above. Other important output files are the:

• OrthoFinder results found in ${dir}/ORTHOGROUPS/OrthoFinder/,
• GO annotations for each orthogroup in ${dir}/ORTHOGROUPS/orthogroups_gene_counts_families_go.out,
• gene family expansion/contraction results in ${dir}/CAFE_results/,
• 4DTv results in plain text files (${dir}/*.4DTv),
• GO term analysis results in plain text files (${dir}/*.goout), and
• results of specific gene analyses:
  • plain text gene family expansion/contraction (${dir}/SPECIFIC_GENES/*.conex),
  • Ka/Ks plots across 15 base-pair windows: (${dir}/SPECIFIC_GENES/*kaks.svg), and
  • Significant Ka/Ks peaks Ka/Ks plots: (${dir}/SPECIFIC_GENES/*PEAKS.csv).

The workflow generates other intermediate files, the purposes of which can be understood by reading the nextflow modules (${compare_genomes_path}/modules/*nf).

Troubleshooting

In cases where IQ-TREE 2 fails to build a time tree (${dir}/ORTHOGROUPS_SINGLE_GENE.NT.timetree.nex), review the estimated divergence times between species in dates.txt and remove divergence time estimates with low confidence based on the number of literature reports and dating studies. Empty files generated by the workflow are indicative of steps that did not completely finish successfully.

While we are working on emitting more informative error messages, please look for the most recent folder and log/error files in compare_genomes/work/<most_recent_folders>/<hash_signature>/. The files are hidden, i.e. .command.sh, .command.log and .command.err.

Bioinformatics tools

• OrthoFinder: Emms, David M., and Steven Kelly. “OrthoFinder: Phylogenetic Orthology Inference for Comparative Genomics.” Genome Biology 20, no. 1 (November 14, 2019): 238. https://doi.org/10.1186/s13059-019-1832-y.was
• HMMER: Mistry, Jaina, Robert D. Finn, Sean R. Eddy, Alex Bateman, and Marco Punta. “Challenges in Homology Search: HMMER3 and Convergent Evolution of Coiled-Coil Regions.” Nucleic Acids Research 41, no. 12 (July 1, 2013): e121. https://doi.org/10.1093/nar/gkt263.
• PantherHMM gene family models: Mi, Huaiyu, Anushya Muruganujan, Dustin Ebert, Xiaosong Huang, and Paul D Thomas. “PANTHER Version 14: More Genomes, a New PANTHER GO-Slim and Improvements in Enrichment Analysis Tools.” Nucleic Acids Research 47, no. D1 (January 8, 2019): D419–26. https://doi.org/10.1093/nar/gky1038.
• MACSE: Ranwez, Vincent, Sébastien Harispe, Frédéric Delsuc, and Emmanuel J. P. Douzery. “MACSE: Multiple Alignment of Coding SEquences Accounting for Frameshifts and Stop Codons.” PLOS ONE 6, no. 9 (September 16, 2011): e22594.
• IQ-TREE: Minh, Bui Quang, Heiko A Schmidt, Olga Chernomor, Dominik Schrempf, Michael D Woodhams, Arndt von Haeseler, and Robert Lanfear. “IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era.” Molecular Biology and Evolution 37, no. 5 (May 1, 2020): 1530–34. https://doi.org/10.1093/molbev/msaa015.
• ModelFinder (module of IQ-TREE): Kalyaanamoorthy, Subha, Minh, Bui Quang, Wong, Thomas KF, von Haeseler, Arndt, and Jermiin, Lars S. ”ModelFinder: Fast model selection for accurate phylogenetic estimates.” Nature Methods, 14 (2017):587–589. https://doi.org/10.1038/nmeth.4285
• CAFE v5: De Bie, Tijl, Nello Cristianini, Jeffery P. Demuth, and Matthew W. Hahn. “CAFE: A Computational Tool for the Study of Gene Family Evolution.” Bioinformatics 22, no. 10 (May 15, 2006): 1269–71. https://doi.org/10.1093/bioinformatics/btl097.
• Gene ontology (GO) enrichment analysis: Mi, Huaiyu, Anushya Muruganujan, Dustin Ebert, Xiaosong Huang, and Paul D Thomas. “PANTHER Version 14: More Genomes, a New PANTHER GO-Slim and Improvements in Enrichment Analysis Tools.” Nucleic Acids Research 47, no. D1 (January 8, 2019): D419–26. https://doi.org/10.1093/nar/gky1038.
• KaKs_calculator v2: Wang, Da-Peng, Hao-Lei Wan, Song Zhang, and Jun Yu. “γ-MYN: A New Algorithm for Estimating Ka and Ks with Consideration of Variable Substitution Rates.” Biology Direct 4, no. 1 (June 16, 2009): 20. https://doi.org/10.1186/1745-6150-4-20.

Citation

Paril, Jefferson, Tannaz Zare, and Alexandre Fournier-Level. “Compare_Genomes: A Comparative Genomics Workflow to Streamline the Analysis of Evolutionary Divergence Across Eukaryotic Genomes.” Current Protocols 3, no. 8 (2023): e876. https://doi.org/10.1002/cpz1.876.

compare_genomes's Issues

Can't generate venn diagram with R script

Hi,

All the intermediate files have been created, but I'm still getting an R error to generate the venn diagram:

`Attaching package: VennDiagram
The following object is masked from package:ape:

rotate

Error in x[[i]] : index out of bounds
Calls: venn.diagram -> unique`

Any suggestions?
Thanks so much!!!

Performance assessment on polyploid genomes

Allow for PANTHER databases to be linked instead of downloaded

Add a parameter here, and modify the setup module accordingly.

Suitable for no model organisms?

Hi, I still need to go through the details. I just installed and started the Test run. We have sequenced and assembled some ciliate genomes that have never been reported. I am a newbie in the field of evolutionary genomics. Can this pipeline be adapted to those non-model organisms?

I really appreciate any help you can provide.

assess_specific_genes.nf fails to open ORTHOGROUPS_SINGLE_GENE.NT.treefile

Hi, I am facing the below mention issue. I have tried to go through the module assess_specific_genes.nf and unable to resolve this on my own and would request your help. I am trying this with the test data provided in github.

"ORTHOGROUPS_SINGLE_GENE.NT.treefile" is generated. but still I am getting the below error.

Define the location of the results of OrthoFinder run, i.e. the most recent output folder.
Extract orthologs per gene (Outputs: {GENE}.ortho)
Move blast output into a seperate directory to clean-up a bit.
Extract gene families (Output: {GENE}.orthocounts)
Infer gene family expansion and contraction (Outpus: {GENE}.conex)
Failed to open ORTHOGROUPS_SINGLE_GENE.NT.treefile

grep: GSTU13_CAFE_results/Gamma_clade_results.txt: No such file or directory
Clean-up.
Output:
(1/3) {GENE}.ortho
(2/3) {GENE}.conex
(3/3) BLASTOUT/

Error executing process > 'PLOT'

Hello!
I tried to run compare_genomes on several species, got an error

Error in file(file, "r") : cannot open the connection Calls: read.nexus -> scan -> file In addition: Warning message: In file(file, "r") : cannot open file 'ORTHOGROUPS_SINGLE_GENE.NT.timetree.nex': No such file or directory Execution halted

Also, most of the output files (e.g. expanded_orthogroup) are empty. Is there something wrong with my files?

P.S. Panther17 is no longer available, will it be better to change from v17 to v18 everywhere?

Thank you.

assess_specific_genes is taking days

Hello,

I am running the last workflow in the pipeline to assess four specific genes. My analysis included 10 plant species. the problem is that, this step is been running for days and it seems to be stuck at the KAKS analysis, see the following:

``
executor > local (6)
[3d/1f0900] process > DOWNLOAD_GENE_AASEQ [100%] 1 of 1 ✔
[db/177126] process > IDENTIFY_GENES_FROM_ORTHOGROUPS [100%] 1 of 1 ✔
[53/7659f6] process > GENE_EXPANSION_CONTRACTION [100%] 1 of 1 ✔
[ad/25dfde] process > EXTRACT_GENE_CDS [100%] 1 of 1 ✔
[14/f4c9be] process > ALIGN_GENE_CDS [100%] 1 of 1 ✔
[91/6cddca] process > KAKS [ 0%] 0 of 1

``
Is this normal?

cannot open file 'ORTHOGROUPS_SINGLE_GENE.NT.timetree.nex': No such file or directory

Hi,

I'm running compare_genomes pipeline and am having this error:

Error in file(file, "r") : cannot open the connection
Calls: read.nexus -> scan -> file
In addition: Warning message:
In file(file, "r") :
cannot open file 'ORTHOGROUPS_SINGLE_GENE.NT.timetree.nex': No such file or directory
Execution halted

I tried doing some things that you have suggested on other issue posts, but nothing solved it.
I changed the Panther database to the 18.0, and also related to the divergence times, I used the estimates that are available on timetree, so I think it should be fine?
This error keep happening and I don't know what else could be done.

If you could please help me, it would be very much appreciated.

Cheers,
Jacqueline

cannot open file 'CAFE_results/Gamma_clade_results.txt': No such file or directory

I noticed that there was an error when performing the mission GENE_EXPANSION_CONTRACTION. I changed line 134 in the file assess_specific_genes.nf. Change $TREE to an absolute path or copy it to the SPECIFIC_GENES folder to keep 176 lines of code working properly. But I am not sure if this is a bug.

Hope to get your reply.

Make use of Nextflow's IO more effectively

And in doing so will enable the effective use of Nextflow's resume capabilities

Rerunning from orthogroups appends species names to protein sequence names which conflicts with downstream module expectations of

Fix it!

RE: 4DTv and WGD

Hello Jefferson,

My analysis includes my species of interest plus six other species. In my "comparisons_4DTv.txt", I listed single species as well as pairs of species, as follows:

``
mySpecies
SpeciesA
SpeciesB
SpeciesC
mySpecies X SpeciesA
mySpecies X SpeciesB
mySpecies X SpeciesC
mySpecies X SpeciesD
mySpecies X SpeciesE
mySpecies X SpeciesF

Species A, and B are the closest species to mySpeices in that order. Species C, D, and E are from a different family. Species F is included as a monocot outgroup. 

My questions:
1) In the 4DTv density plot, I only see  plots for the comparisons, but not the single species that i listed. Is this expected? 

2) All plots have non-zero peaks at different 4DTv accumulation rates. I am a bit unclear about the meaning of it. Does it mean that there was one WGD event but the different 4DTv peaks are just related to the divergence time between the two species? From the WF script, I see that it uses the pairwise paralogs and pairwise orthologs of the compared species. For instance, I assume, the density plot for "mySpecies X Species A" include the 4DTV rates from the pairwise paralogs from mySpecies and SpeciesA plus the pairwise Orthologs between mySpecies and SpeciesA. And so on. that's why the plots have different peaks, but all refer to a single WGD events.  Is this accurate understanding? 

2) In the main output directory, there is a "*.4DTv" file for every species in the analysis. When I made a density plot for these files using the 4th column (assumed to be 4DTv rates) of each file, the plot is totally different from the plot generated by the pipeline.  the files have the statistics for only pairwise genes, which is what the WF uses I assume. so, why the plots are different? in fact, my density plot of the 4DTv of mySpecies (attached), has only a zero-peak, which is different from what the pipeline produces in the comparisons.



![image](https://github.com/jeffersonfparil/compare_genomes/assets/9359920/ea3ddb6d-0e05-45f8-a9f8-e875ce456e81)



I am sure I am missing many things here, I really appreciate it if you could point them out. 
Many thanks.

Error on test dataset

Hi,

So, I did not have access to PANTHER 17 files, and so I did modify files to PANTHER 18. It is possible this is the cause somehow, but I am skeptical because of where the error occurs.

The pipeline completes Orthofinder, as far as I can tell, but fails after that. I see this error message in command.err :

Traceback (most recent call last):
  File "/work/soghigian_lab/apps/conda/envs/compare_genomes/bin/orthofinder", line 7, in <module>
    main(args)
  File "/work/soghigian_lab/apps/conda/envs/compare_genomes/bin/scripts_of/__main__.py", line 1778, in main
    GetOrthologues(speciesInfoObj, options, prog_caller)
  File "/work/soghigian_lab/apps/conda/envs/compare_genomes/bin/scripts_of/__main__.py", line 1540, in GetOrthologues
    orthologues.OrthologuesWorkflow(speciesInfoObj.speciesToUse, 
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/work/soghigian_lab/apps/conda/envs/compare_genomes/bin/scripts_of/orthologues.py", line 999, in OrthologuesWorkflow
    spTreeFN_ids, qSTAG = db.RunAnalysis(userSpeciesTree == None)
                          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/work/soghigian_lab/apps/conda/envs/compare_genomes/bin/scripts_of/orthologues.py", line 481, in RunAnalysis
    stag.Run_ForOrthoFinder(files.FileHandler.GetOGsTreeDir(), files.FileHandler.GetWorkingDirectory_Write(), self.ogSet.seqsInfo.speciesToUse, spTreeFN_ids)
  File "/work/soghigian_lab/apps/conda/envs/compare_genomes/bin/scripts_of/stag.py", line 262, in Run_ForOrthoFinder
    ProcessTrees(dir_in, dir_matrices, dir_trees_out, gene_to_species, qVerbose=False, qForOF=True)
  File "/work/soghigian_lab/apps/conda/envs/compare_genomes/bin/scripts_of/stag.py", line 191, in ProcessTrees
    t = tree.Tree(fn)
        ^^^^^^^^^^^^^
  File "/work/soghigian_lab/apps/conda/envs/compare_genomes/bin/scripts_of/tree.py", line 221, in __init__
    read_newick(newick, root_node = self, format=format)
  File "/work/soghigian_lab/apps/conda/envs/compare_genomes/bin/scripts_of/newick.py", line 208, in read_newick
    nw = open(newick, 'rU').read()
         ^^^^^^^^^^^^^^^^^^
ValueError: invalid mode: 'rU'

HyPhy integration for more efficient detection of signatures of selection on gene sequences

Improve documentation and debug Venn diagram plotting

• Make sure underscore or space requirements in the species names are properly documented
• Add helpful hints to fix conflicts with dates and genotype distances in iqtree's tree building step
• The Venn diagram can accommodate a maximum of 5 species, add a config file specifying the species the user wants to include in the diagram

Diagnosing pipeline failure

Hi there!

Thanks again for the great tool. I've been trying it on my personal data and I think it is finishing the OrthoFinder step now. I am not sure exactly where it fails, but I think it may be failing to complete aligning single copy genes. On the macse alignment step I noticed the error log had a number of messages like this:

Exception in thread "main" java.lang.OutOfMemoryError: Java heap space at align.ProfileAligner.<init>(ProfileAligner.java:120) at align.CodingMSA.<init>(CodingMSA.java:64) at align.CodingMSA.buildAlignmentReliable(CodingMSA.java:633) at align.CodingMSA.run(CodingMSA.java:659) at utils.MacseMain.main(MacseMain.java:426)

I have tried various configurations for memory and processors (last run was 3 processors + 120 gb of memory)...

This is a set of 5 insect genomes...

Process `PLOT` terminated with an error exit status (1)

Dear compare_genomes developers,

I'd like to seek for your help on the possible resolution of an error encountered while running the pipeline particularly at the "plot_tree_conex_venn_4DTv.nf". The log file says "cannot open file 'ORTHOGROUPS_SINGLE_GENE.NT.timetree.nex': No such file or directory" which could have likely led to "Process PLOT terminated with an error exit status (1)". Indeed, this *nex file is missing in the output. I was wondering if by any chance you would know how to address this issue. It appears that besides this step, all other steps were executed successfully with output files. Looking forward to your response! Thanks a lot!