4 `prot_finder.pl` is a script to search for query protein homologs in annotated bacterial genomes with **BLASTP**. A companion bash shell script pipeline is available, `prot_finder_pipe.sh`.
6 Included is also the script `prot_binary_matrix.pl` to create a binary presence/absence matrix (e.g. for [**iTOL**](http://itol.embl.de/)) from the `prot_finder.pl` output. Additionally, two downstream scripts are provided to wrangle these binary presence/absence matrix: `transpose_matrix.pl` to transpose a delimited TEXT matrix and `binary_group_stats.pl` to get overall presence/absence statistics for groups of columns in a delimited binary TEXT matrix (in the style of [`po2group_stats.pl`](/po2group_stats)).
89 The script `prot_finder.pl` is intended to search for homologous proteins in annotated bacterial genomes. For this purpose, a previous [**BLASTP**](http://blast.ncbi.nlm.nih.gov/Blast.cgi), either [legacy or plus](https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastDocs&DOC_TYPE=Download), needs to be run with query protein sequences against a **BLASTP** database of subject proteins (e.g. all proteins from several *Escherichia coli* genomes).
91 The script [`cds_extractor.pl`](/cds_extractor) (with options **-p -f**) can be used to create multi-FASTA protein files of all non-pseudo CDS from RichSeq genome files to create the needed subject **BLASTP** database. Present locus tags will be used as FASTA IDs, but see [`cds_extractor.pl`](/cds_extractor) for a description of the format. Query protein sequences for the **BLASTP** need a **unique** FASTA ID.
93 The **BLASTP** report file (option **-r**), the subject protein multi-FASTA file (option **-s**), and optionally the query protein (multi-)FASTA (option **-q**) file are then given to `prot_finder.pl`. Significant **BLASTP** subject hits are filtered according to the given cutoffs (options **-i**, **-cov_q**, and **-cov_s**) and the result is printed as an informative tab-separated result table to *STDOUT*. To apply global identity/coverage cutoffs to subject hits high-scoring pairs (HSPs) are tiled (see http://www.bioperl.org/wiki/HOWTO:Tiling and http://search.cpan.org/dist/BioPerl/Bio/Search/Hit/GenericHit.pm). Additionally, the subject protein sequences with significant query hits are written to result multi-FASTA files, named according to the respective query FASTA IDs (optionally including the query sequence with option **-q**).
95 Optionally, [**Clustal Omega**](http://www.clustal.org/omega/) can be called (option **-a** with optional **-p**) to create multiple alignments (FASTA format) for each of the resulting multi-FASTA files. These alignments can be used to calculate phylogenies e.g. with **RAxML** (http://sco.h-its.org/exelixis/software.html) or **MEGA** (http://www.megasoftware.net/).
97 Run the script [`cds_extractor.pl`](/cds_extractor) (with options **-p -f**) and the **BLASTP** manually or use the bash shell wrapper script `prot_finder_pipe.sh` (see below ['prot_finder_pipe bash script pipeline'](#prot_finder_pipe-bash-script-pipeline)) to execute the whole pipeline including `prot_finder.pl` (with optional option **-q**). For additional options of the pipeline shell script see below ['prot_finder_pipe.sh options'](#prot_finder_pipesh-options). Be aware that some options in `prot_finder_pipe.sh` corresponding to options in `prot_finder.pl` have different names (**-c** instead of **-cov_q**, **-k** instead of **-cov_s**, and **-o** instead of **-p**; also **-f** has a different meaning). If [`cds_extractor.pl`](/cds_extractor) is used in the pipeline (option **-f** of the shell script) the working folder has to contain the annotated bacterial genome subject files (in RichSeq format, e.g. EMBL or GENBANK format). Also, the Perl scripts [`cds_extractor.pl`](/cds_extractor) (only for `prot_finder_pipe.sh` option **-f**) and `prot_finder.pl`have to be either contained in the current working directory or installed in the global *PATH*. **BLASTP** (legacy and/or plus) and **Clustal Omega** binaries have to be installed in global *PATH*, or for **Clustal Omega** you can give the path to the binary with option **-o**. In the pipeline **BLASTP** is run with **disabled** query filtering, locally optimal Smith-Waterman alignments, and increasing the number of database sequences to show alignments to 500 for [**BioPerl**](http://www.bioperl.org) parsing (legacy: **-F F -s T -b 500**, plus: **-seg no -use_sw_tback -num_alignments 500**). The pipeline script ends with the *STDERR* message 'Pipeline finished!', if this is not the case have a look at the log files in the result directory for errors.
99 The resulting tab-separated table with significant **BLASTP** hits (from `prot_finder.pl` or `prot_finder_pipe.sh`) can be given to the script `prot_binary_matrix.pl`, either as *STDIN* or as a file, to create a presence/absence matrix of the results. See below ['prot_binary_matrix.pl options'](#prot_binary_matrixpl-options) for the `prot_binary_matrix.pl` options. By default a tab-delimited binary presence/absence matrix for query hits per subject organism will be printed to *STDOUT*. Use option **-t** to count all query hits per subject organism, not just the binary presence/absence. This presence/absence matrix can be given to the script `transpose_matrix.pl`, either as *STDIN* or as a file, to transpose the matrix, i.e. rows will become columns and columns rows. Actually, `transpose_matrix.pl` can be used to transpose any delimited TEXT matrix (see below ['transpose_matrix.pl options'](#transpose_matrixpl-options)).
101 The presence/absence matri(x|ces) can e.g. be loaded into the Interactive Tree Of Life website ([**iTOL**](http://itol.embl.de/)) to associate the data with a phylogenetic tree. [**iTOL**](http://itol.embl.de/) likes individual comma-separated input files, thus use `prot_binary_matrix.pl` options **-s -c** for this purpose. However, the organism names have to have identical names to the leaves of the phylogenetic tree, thus manual adaptation, e.g. in a spreadsheet software (like [**LibreOffice Calc**](https://www.libreoffice.org/discover/calc/)), might be needed. **Careful**, subject organisms without a significant **BLASTP** hit won't be included in the `prot_finder.pl` result table and hence can't be included by `prot_binary_matrix.pl`. If needed add them manually to the result matri(x|ces).
103 At last, script `binary_group_stats.pl` can be used to categorize columns of the binary presence/absence matrix from `prot_binary_matrix.pl` according to group affiliations. `binary_group_stats.pl` is based upon [`po2group_stats.pl`](/po2group_stats), which does the same thing for genomes in an ortholog/paralog output matrix from a [Proteinortho5](http://www.bioinf.uni-leipzig.de/Software/proteinortho/) calculation. Actually, `binary_group_stats.pl` can work with any delimited TEXT **binary** matrix (option **-i**). But, all fields of the binary matrix need to be filled with either a **0** indicating absence or a **1** indicating presence, i.e. all rows need to have the same number of columns. Use option **-d** to set the delimiter of the input matrix, default is set to tab-delimited/separated matrices.
105 Also, column headers in the first row and row headers in the first column are **mandatory** for the input binary matrix. Only alphanumeric (a-z, A-Z, 0-9), underscore (_), dash (-), and period (.) characters are allowed for the **column headers** and **group names** in the group file (option **-g**) to avoid downstream problems with the operating/file system. As a consequence, also no whitespaces are allowed in these! Additionally, **column headers**, **row headers**, and **group names** need to be **unique**.
107 The group affiliations of the columns are intended to get overall presence/absence statistics for groups of columns and not simply single columns of the matrix. Percentage inclusion (option **-cut_i**) and exclusion (option **-cut_e**) cutoffs can be set to define how strict the presence/absence of column groups within a row are defined. Of course groups can also hold only single column headers to get single column statistics. Group affiliations are defined in a mandatory **tab-delimited** group input file (option **-g**), including the column headers from the input binary matrix, with **minimal two** and **maximal four** groups.
109 See the *README.md* of [`po2group_stats.pl`](/po2group_stats) for an explanation of the logic behind the categorization and the resulting group binary matrix and venn diagram of `binary_group_stats.pl` (and of course its [options](#binary_group_statspl-options) and [output](#binary_group_statspl-output)).
195 Path to subject multi-FASTA protein sequence file (\*.faa) created with [`cds_extractor.pl`](/cds_extractor) (and its options **-p -f**), which was used to create the **BLASTP** database
209 Force output to an existing result folder, otherwise ask user to remove content of existing folder. Careful, files from a previous analysis might not be overwritten if different to current analysis.
213 Path to query (multi-)FASTA protein sequence file (\*.faa) with **unique** FASTA IDs, which was used as query in the **BLASTP**. Will include each query protein sequence in the respective multi-FASTA 'query-ID_hits.faa' result file.
257 File extension for files in the **current** working directory to use for [`cds_extractor.pl`](/cds_extractor) (e.g. 'embl' or 'gbk'); excludes shell script option **-s**
263 Path to subject protein multi-FASTA file (\*.faa) already created with [`cds_extractor.pl`](/cds_extractor) (and its options **-p -f**), will not run [`cds_extractor.pl`](/cds_extractor); excludes shell script option **-f**
309 Path to executable **Clustal Omega** binary if not in global *PATH*; requires shell script option **-a** (corresponds to [`prot_finder.pl` option](#optional-prot_finderpl-options) **-p**)
323 Separate presence/absence files for each query protein printed to the result directory [default without **-s** = *STDOUT* matrix for all query proteins combined]
339 Use the locus_tag **prefixes** in the subject_ID column of the `prot_finder.pl` output (instead of the subject_organism column) as organism IDs to associate query hits to organisms. The subject_ID column will include locus_tags if they're annotated for a genome (see the [`cds_extractor.pl`](/cds_extractor) format description). Useful if the [`cds_extractor.pl`](/cds_extractor) output doesn't include strain names for 'o=' in the FASTA IDs, because the prefix of a locus_tag should be unique for a genome (see http://www.ncbi.nlm.nih.gov/genbank/genomesubmit_annotation).
373 Tab-delimited file with group affiliation for the columns from the input binary matrix with **minimal two** and **maximal four** groups (easiest to create in a spreadsheet software and save in tab-separated format). **All** column headers from the input binary matrix need to be included. Column headers and group names can only include alphanumeric (a-z, A-Z, 0-9), underscore (_), dash (-), and period (.) characters (no whitespaces allowed either). Example format with two column headers in group A, three in group B and D, and one in group C:
396 Percentage inclusion cutoff for column presence counts in a group per row, has to be > 0 and <= 1. Cutoff will be rounded according to the column header number in each group and has to be > the rounded exclusion cutoff in this group. \[default = 0.9\]
408 Plot venn diagram from the group binary matrix (except 'unspecific' and 'underrepresented' categories) with function `venn` from **R** package **gplots**, requires option **-b**
416 Include singleton/column-specific rows for each column header in the output, activates also overall column header presence ('1') counts in final stats matrix for columns with singletons
424 Report overall presence counts for all column headers (appended to the final stats matrix), also those without singletons; will include all overall column header presence counts without option **-s**
436 Multi-FASTA file(s) of subject CDS protein sequences; will be removed with [`prot_finder_pipe.sh` option](#optional-prot_finder_pipesh-options) **-m**
444 The resulting tab-delimited output table with the significant subject **BLASTP** hits is printed to *STDOUT*. Redirect (e.g. to a file in the result directory, options **-d -f**) or pipe into another tool as needed (e.g. `prot_binary_matrix.pl`).
452 Multi-FASTA protein files of significant subject hits for each query protein (named after the respective query FASTA ID), optionally includes the respective query protein sequence (with option **-q**)
456 Index file of the subject protein file for fast sequence retrieval (can be deleted if no further **BLASTPs** are needed with these subject sequences); will be removed with [`prot_finder_pipe.sh` option](#optional-prot_finder_pipesh-options) **-m**
460 Lists all query sequence IDs without significant subject hits; with option **-b** includes also queries with significant hits but *without* a best blast hit for a subject
464 Optional log file of verbose **Clustal Omega** *STDOUT*/*STDERR* messages; will be removed with [`prot_finder_pipe.sh` option](#optional-prot_finder_pipesh-options) **-m**
476 In addition to the [`cds_extractor.pl` output](#cds_extractorpl-output) and the [`prot_finder.pl` output](#prot_finderpl-output) the pipeline also creates the following non-essential output files, which will be removed with [`prot_finder_pipe.sh` option](#optional-prot_finder_pipesh-options) **-m**:
579 For [citation](https://github.com/aleimba/bac-genomics-scripts#citation), [installation](https://github.com/aleimba/bac-genomics-scripts#installation-recommendations), and [license](https://github.com/aleimba/bac-genomics-scripts#license) information please see the repository main [*README.md*](https://github.com/aleimba/bac-genomics-scripts/blob/master/README.md).
598 * tab-delimited output table with filtered **BLASTP** hits now printed to *STDOUT* instead of file and can be used as *STDIN* for `prot_binary_matrix.pl`
599 * result files now created in result directory to unclutter working dir (new options **-dir_result** and **-force_dir**) and replaced subroutine 'file_exist' with 'empty_dir' (also no *STDOUT* message which files were created anymore)
607 * changes to subroutine 'split_fasta_header' to work more robustly with `split` instead of a regex and adapted to `cds_extractor.pl` v0.7+ output format
618 * new option **-best_hit** to include only best **BLASTP** hit (highest identity) for each subject locus_tag, if a subject protein has hits to several query proteins
621 * corrected bug in hash structure to store **BLASTP** hits, query acc/IDs (keys) and array reference of subject locus_tags (values), as several queries can have the same subject locus_tag as hit
637 * New option **-locus_tag** to use the locus_tag prefix in column subject_ID of the `prot_finder.pl` output as ID, which is in most cases the locus_tag, instead of the subject_organism; controls if the locus_tags follow NCBI standards
