Script for running reciprocal best blast hit (RBH) analysis of one query genome against one or multiple target genomes.
-q : Path to query genome as amino acid multi-FASTA file (.faa).
-i : Path to input directory. Location of target genomes (.faa files).
-o : Path to output directory. Location for output files to be created.
-t : OPTIONAL: Number of threads (CPUs) to use in the BLAST search. Default = 1.
-m : OPTIONAL: Modify headers of target genomes. If selected, FASTA headers of genome files specified with "-i" will be modified, replacing blank spaces (" ") with underscores ("_"). This prevents BLAST from truncating headers at blank spaces when assigning protein sequence names.
-c : OPTIONAL. If selected, C. crescentus CB15 annotations found in "scripts/CC_CCNA_conversion_table.tsv" will be added to output table. This requires that the query used is C. crescentus CB15, and that the query file is formatted properly. See "Formatting the query", below.
bash make_RBH.sh -q GCF_000006905.1_Caulobacter_vibrioides_CB15.NH.faa -i RBH_input/ -o RBH_output/ -t 4 -m -c
In our use case, the fasta headers of the query genome (GCF_000006905.1) were formatted accordingly, using "@" to separate variables, and placing protein accession IDs in third position:
>assembly@taxon_name@protein_accession@annotation
An example FASTA header for the "GCF_000006905.1_Caulobacter_vibrioides_CB15.NH.faa" file used as query:
>GCF_000006905.1@Caulobacter_vibrioides_CB15@WP_004615089.1@MULTISPECIES: 30S ribosomal protein S18
The script assumes that the query headers are formatted accordingly.
When including the parameter "-c", C. crescentus CB15/NA1000 annotations from "CC_CCNA_conversion_table.tsv" are added to the output table by joining with the query sequence variable "protein_accession". This essentially requires the use of the CB15 genome as query.
"make_RBH.sh" needs to be called from a directory that contains the "scripts/" folder.
The shell script "runrblast_JH.sh" requires command-line BLAST to be installed (uses makeblastdb and blastp) and the R script "generate_RBH_summary_table.R" uses the R library tidyverse.
The shell scripts "runrblast_JH.sh" and ""findrblast_JH.sh" are based on the scripts from Elizabeth Brooks (2020-08-18), accessed 2023-12-04, available at: https://morphoscape.wordpress.com/2020/08/18/reciprocal-best-hits-blast-rbhb/