Find Homologous Proteins Between Clusters by Pairwise Alignment

Compares the proteins of a query cluster with those of the other clusters. It generates all protein combinations between the query cluster and each other cluster, aligns each pair with a global (Needleman-Wunsch) pairwise alignment from pwalign::pairwiseAlignment(), calculates sequence identity and similarity, and keeps the best hit of each gene above a minimum identity threshold.

Usage

protein_blast(
  data,
  query,
  id = "protein_id",
  start = "start",
  end = "end",
  cluster = "cluster",
  genes = NULL,
  identity = 30,
  parallel = TRUE
)

Arguments

Argument	Description
data	A dataframe or a character vector specifying the path to .gbk files. When a character vector is provided, it is interpreted as file paths to .gbk files which are then read and processed. The dataframe must contain columns for unique protein identifiers, cluster identifiers, protein sequences, and the start and end positions of each gene.
query	The name of the query cluster the other clusters are compared with.
id	The name of the column that contains the gene identifiers. Defaults to "protein_id".
start	The name of the column specifying the start positions of genes. Defaults to "start".
end	The name of the column specifying the end positions of genes. Defaults to "end".
cluster	The name of the column specifying the cluster names. Defaults to "cluster".
genes	An optional vector of gene identifiers to include in the analysis. Defaults to NULL.
identity	Minimum percentage identity of the global alignment for a hit to be kept. Defaults to 30.
parallel	Deprecated and ignored. The alignments are now fast enough without parallel processing. Kept so that existing code keeps working.

Value

A modified version of the input data dataframe, including additional columns BlastP (the protein of the query cluster that is the best hit; the name is kept for compatibility), identity, similarity and score (the synteny score of the cluster). Columns of data with these names, e.g. from an earlier call, are replaced.

Details

Despite its name, the function does not run BLAST: there is no seed search, local alignment, substitution matrix or E-value. Identity and similarity are percentages of the full length of the global alignment, gaps included, so proteins of very different length, or that share only one domain, get lower values than they would get from BLASTP. Set identity with this in mind.

Note

This function relies on the Biostrings and pwalign packages from Bioconductor for sequence alignment. Ensure these packages are installed.

References

Needleman, S. B. and Wunsch, C. D. (1970). A general method applicable to the search for similarities in the amino acid sequence of two proteins. Journal of Molecular Biology, 48(3), 443-453. doi:10.1016/0022-2836(70)90057-4

Examples

if (requireNamespace("Biostrings", quietly = TRUE) &&
    requireNamespace("pwalign", quietly = TRUE)) {
  genes <- data.frame(
    protein_id = c("A1", "A2", "B1", "B2"),
    cluster = c("A", "A", "B", "B"),
    start = c(1, 200, 1, 200),
    end = c(150, 350, 150, 350),
    translation = c("MKTAYIAKQRQISFVKSHFSRQ", "MSDNELLKRLGEAVLAHPTTG",
                    "MKTAYIAKQRQISFVKSHFSRQ", "MLPGGRTSWVQEHCDKAYNIF")
  )
  protein_blast(genes, query = "A", identity = 30)
}