Workflow for identifying and classifying homologous gene/protein sequences

public 1yr ago Version: v3.0 0 bookmarks

View Workflow

Help improve this workflow!

This workflow has been published but could be further improved with some additional meta data:

Keyword(s) in category input

You can help improve this workflow by suggesting the addition or removal of keywords, suggest changes and report issues, or request to become a maintainer of the Workflow .

Introduction

A nalysis of MO lecular E volution with BA tch E ntry (AMOEBAE) is a customizable bioinformatics workflow for identifying homologues (and potential orthologues) of genes of interest among a mid-size sampling of genomes. This workflow is designed to be run on high-performance computing (HPC) clusters and is executed via the SnakeMake workflow management system. Code for steps in this workflow is written primarily in Python3 , relying heavily on the Biopython library, and apply bioinformatics packages including BLAST+ and HMMER3 to input data files. AMOEBAE is open-source, and all dependencies are freely available. Lael D. Barlow is the author.

AMOEBAE is useful for certain mid-scale comparative genomics studies that might otherwise require time-intensive and repetitive manual/visual manipulation of data. Webservices such as those provided by NCBI and EMBL-EBI provide a means to readily investigate the evolution of one or a few genes via similarity searching, and large-scale analysis workflows such as OrthoMCL and OrthoFinder attempt to rapidly perform orthology prediction for all genes among several genomes. AMOEBAE addresses analyses which are too cumbersome to be performed via webservices or simple scripts and yet require a level of detail and flexibility not offered by large-scale analysis workflows. AMOEBAE is useful for analyzing the distribution of homologues of up to approximately 30 genes/proteins among a sampling of no more than approximately 100 eukaryotic genomes, especially when follow-up with custom phylogenetic analysis is planned.

AMOEBAE serves this purpose by providing several unique features. The core functionality of AMOEBAE is to run sequence similarity searches with multiple algorithms, multiple queries, and multiple databases simultaneously and to allow highly customizable implementation of reciprocal-best-hit search strategies. The output includes detailed summaries of results in the form of a spreadsheet and presence/absence plots. A particular advantage of AMOEBAE compared to other workflows is its functionality for parsing results of TBLASTN (which searches nucleotide sequences with peptide sequence queries) search results. This allows rapid identification of High-scoring Segment Pair (HSP) clusters at separate gene loci, automatic checking of those loci against information in genome annotation files, and systematic use of the Exonerate package where possible for obtaining exon predictions. In addition, AMOEBAE provides many options which can be tailored to the specific genes/proteins being analyzed. Despite the complexity of this workflow, analyses performed using AMOEBAE can be reproduced via SnakeMake .

The output files include a plot of the number of identified homologues (potential orthologues) of several genes across several genomes, as well as a spreadsheet in CSV format providing a detailed summary of search results.

Here's a diagram of the steps in the overall workflow:

Here's an example coulson plot output by the workflow:

Protocol

See the workflow protocol for instructions and guidelines for running the AMOEBAE workflow.

Questions and bug reporting

Please use the issue tracker on the GitHub webpage to report any problems you encounter while using AMOEBAE.

How to cite AMOEBAE

Please cite the AMOEBAE GitHub repository (or alternative permanent repositories if relevant).

Acknowledgments

AMOEBAE was initially developed at the Dacks Laboratory at the University of Alberta, and was supported by National Sciences and Engineering Council of Canada (NSERC) Discovery grants RES0021028, RES0043758, and RES0046091 awarded to Joel B. Dacks, as well as an NSERC Postgraduate Scholarship-Doctoral awarded to Lael D. Barlow.

We acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC).

Cette recherche a été financée par le Conseil de recherches en sciences naturelles et en génie du Canada (CRSNG).

Also, help with testing AMOEBAE has been kindly provided by numerous members of the Dacks laboratory.

License

Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at

http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

Code Snippets

DOCTEX=amoebae_commands
pdflatex --shell-escape $DOCTEX.tex &&\
bibtex $DOCTEX &&\
makeglossaries $DOCTEX &&\
pdflatex --shell-escape $DOCTEX.tex &&\
pdflatex --shell-escape $DOCTEX.tex

Shell From line 2 of documentation/build_pdf.sh

run:
    for rule in workflow.rules:
        print(rule.name)
        print(rule.docstring)

SnakeMake From line 106 of workflow/Snakefile

shell:
    """
	python3 {input.script} {input.csv_file} results/temp_databases
    """

SnakeMake From line 123 of workflow/Snakefile

shell:
    """
	python3 {input.script} {input.csv_file} \
            results/temp_single_seq_queries \
            results/temp_multi_seq_queries
    """

SnakeMake From line 143 of workflow/Snakefile

shell:
    'sc=$( python3 {input.split_char_script} ' + database_csv + ' {input.input_file} ) && '
    'sp=$( python3 {input.split_pos_script} ' + database_csv + ' {input.input_file} ) && '
    'python3 {input.amoebae_script} '
     'add_to_dbs '
         '{input.input_file} '
         '{input.data_dir} '
         '--split_char $sc '
         '--split_pos $sp '

SnakeMake From line 176 of workflow/Snakefile

shell:
    'sc=$( python3 {input.split_char_script} ' + database_csv + ' {input.input_file} ) && '
    'sp=$( python3 {input.split_pos_script} ' + database_csv + ' {input.input_file} ) && '
    'python3 {input.amoebae_script} '
        'add_to_dbs '
            '{input.input_file} '
            '{input.data_dir} '
                 '--split_char $sc '
                 '--split_pos $sp '

SnakeMake From line 212 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
    'add_to_queries '
        '{input.input_file} '
        '{input.data_dir} '

SnakeMake From line 243 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
    'add_to_queries '
        '{input.input_file} '
        '{input.data_dir} '

SnakeMake From line 271 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
        'align_fa '
            '{input.input_file} '
            '--output_format fasta '

SnakeMake From line 296 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
       'add_to_queries '
           '{input.input_file} '
           '{input.data_dir} '

SnakeMake From line 323 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
       'add_to_queries '
           '{input.input_file} '
           '{input.data_dir} '

SnakeMake From line 351 of workflow/Snakefile

run:
    # Get list of .faa and .fna database files (excluding .gff3 files).
    fa_database_names_all = [x for x in database_names_all if \
                             x.endswith('.faa') or x.endswith('.fna')]

    # Get alignment names from multifasta names.
    alignment_fasta_query_names = [x.rsplit('.', 1)[0] + '.afaa' for x in \
    multi_fasta_query_names_all]

    # Write lists to files.
    with open(output.dbs, 'w') as dbsfh, open(output.queries, 'w') as queriesfh:
        # Write database name list to file.
        dbsfh.write('\n'.join(fa_database_names_all))
        # Write query name list to file.
        queriesfh.write('\n'.join(single_fasta_query_names_all + alignment_fasta_query_names))

    # Write basic pre-filled coulson plot organization CSV file.
    with open(output.coulson_plot_org, 'w') as coulson_plot_org:
        # Make a non-redundant list of query names.
        nonredun_query_names = list(set(single_fasta_query_names_all + \
                               alignment_fasta_query_names))

        # Get a list of just the query titles from the query names.
        nonredun_query_titles = list(set([x.split('_', 1)[0].rsplit('.', 1)[0] for x in nonredun_query_names]))
        nonredun_query_titles.sort()

        # Write lines, just using the first character of the query name as the
        # category name.
        for i in nonredun_query_titles:
            coulson_plot_org.write(i[0] + ',' + i + '\n')

SnakeMake From line 373 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} ' + \
        'get_redun_hits ' + \
            'results/Ref_seqs_1 ' + \
            'results/AMOEBAE_Data ' + \
            '--db_list_file {input.db_list_file} ' + \
            '--query_list_file {input.query_list_file} ' + \
            '--csv_file ' + \
            'results/Ref_seqs_1/0_ref_seqs.csv ' + \
            '--predict_redun_hit_selection ' + \
            '--max_number_of_hits_to_summarize 200 ' + \
    ' && cp results/Ref_seqs_1/0_ref_seqs_auto_predictions.csv ' + \
           'results/Ref_seqs_1_auto_predictions.csv'

SnakeMake From line 430 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} ' + \
        'setup_fwd_srch ' + \
            'results ' + \
            '{input.query_list_file} ' + \
            '{input.db_list_file} ' + \
            '--outdir {output.outdir} ' + \
    '&& ' + \
    'python3 {input.amoebae_script} ' + \
        'run_fwd_srch ' + \
            '{output.outdir} ' + \
            'results/AMOEBAE_Data '

SnakeMake From line 468 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
        'sum_fwd_srch ' + \
            '{input.fwd_srch_out} ' + \
            '{output.csv_out_path} ' + \
            '{input.main_data_dir} ' + \
            '--max_hits_to_sum 5 ' + \
            '--max_gap_between_tblastn_hsps 1000 '

SnakeMake From line 501 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
        'setup_rev_srch ' + \
            'results ' + \
            '{input.fwd_srch_csv} ' + \
            '{input.reverse_srch_db_file} ' + \
            '{input.main_data_dir} ' + \
            '--outdir  {output.out_dir} ' + \
            '--aasubseq ' + \
    '&& ' + \
    'python3 {input.amoebae_script} '
        'run_rev_srch ' + \
            '{output.out_dir} ' + \
            '{input.main_data_dir} '

SnakeMake From line 531 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
        'sum_rev_srch ' + \
            '{input.fwd_srch_csv} ' + \
            '{input.rev_srch_dir} ' + \
            '{output.csv_out_path} ' + \
            '{input.main_data_dir} ' + \
            '--redun_hit_csv {input.ref_seq_csv} ' + \
            '--aasubseq ' + \
            '--min_evaldiff 2 '

SnakeMake From line 568 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
        'interp_srchs ' + \
            '{input.rev_srch_csv} ' + \
            '--out_csv_path {output.interp_csv} '

SnakeMake From line 597 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
        'find_redun_seqs ' + \
            '{input.interp_csv} ' + \
            '{input.main_data_dir} ' + \
            '--out_csv_path {output.interp_csv_with_ali_col} ' + \
            '--add_ali_col ' + \
    '&& ' + \
    'python3 {input.amoebae_script} '
        'find_redun_seqs ' + \
            '{output.interp_csv_with_ali_col} ' + \
            '{input.main_data_dir} ' + \
            '--out_csv_path {output.nonredun_seq_csv} ' + \
            '--min_length 55'

SnakeMake From line 628 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
        'plot ' + \
            '{input.nonredun_seq_csv} ' + \
            '--complex_info {input.coulson_org_csv} ' + \
            '--row_order {input.row_order} ' + \
            '--out_pdf {output.main_plot_pdf} '

SnakeMake From line 669 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
        'csv_to_fasta ' + \
            '{input.nonredun_seq_csv} ' + \
            '--abbrev ' + \
            '--split_by_query_title '

SnakeMake From line 696 of workflow/Snakefile

shell:
    'python3 {input.amoebae_script} '
        'csv_to_fasta ' + \
            '{input.fwd_hit_seq_csv} ' + \
            '--abbrev ' + \
            '--split_by_query_title ' + \
            '--all_hits '

SnakeMake From line 722 of workflow/Snakefile

shell:
    'mkdir -p {output.ali_dir} && '
    'for X in {input.fasta_dir}/*.fa; '
    'do '
    'python3 {input.amoebae_script} ' + \
        'align_fa ' + \
            '{input.fasta_dir}/$( basename $X ) ' + \
            '--output_dir_path {output.ali_dir} ' + \
            '--output_format fasta '
    '; done '

SnakeMake From line 750 of workflow/Snakefile

shell: 
    'bash documentation/build_pdf.sh'

SnakeMake From line 772 of workflow/Snakefile

shell: 
    'snakemake --cores 1 -p --rulegraph | dot -Tpdf > {output} && '
    'snakemake --cores 1 -p --rulegraph | dot -Tpng > images/example_workflow_diagram.png '