inst/extcode/blastSeq.R

In WubingZhang/PhageR: Functions for phage display data analysis

#' Blastp
#'
#' @docType methods
#' @name blastSeq
#' @rdname blastSeq
#'
#' @param sequence A vector of sequences.
#' @param fa A fasta file in which sequences are applied to blast.
#' @param word_size Parameter used in blast, e.g. nchar(seq)/2.
#' @param gapopen Parameter used in blast.
#' @param gapextend Parameter used in blast.
#' @param db  Database to be search.
#' @param blast_dir Directory to blast API.
#'
#' @return A data frame with blast results.
#'
#' @author Wubing Zhang
#' @import Biostrings
#' @import rBLAST
#' @export
#'
blastSeq <- function(sequence = NULL, fa = NULL, word_size = 3, gapopen = 5, gapextend = 2,
                     db = "~/Jobs/Project/XunBaihui/_Data/6_PhageDisplay/TB_Proteosome/UP000001584_83332.fasta",
                     blast_dir = "/Applications/ncbi-blast-2.9.0+/bin/"){
  require(Biostrings)
  require(rBLAST)
  if(is.null(sequence)) peptide.seq <- readAAStringSet(fa)
  if(is.null(fa)){
    fa = tempfile("peptide", fileext = c(".fa"))
    names(sequence) = paste0(">", sequence)
    write.table(sequence, fa, quote = FALSE, col.names = FALSE, sep = "\n")
    peptide.seq <- readAAStringSet(fa)
  }
  if(is.null(sequence)&is.null(fa)) stop("No sequence is found !!!")

  Sys.setenv(PATH = paste(Sys.getenv("PATH"), blast_dir, sep= .Platform$path.sep))
  makeblastdb(db, dbtype = "prot")
  TB_Proteosome <- blast(db, type = "blastp")
  blastRes = predict(TB_Proteosome, peptide.seq,
                     BLAST_args = paste0("-word_size ", floor(word_size),
                                         " -evalue 20000 -matrix PAM30 -gapopen ", gapopen,
                                         " -gapextend ", gapextend, " -comp_based_stats 0"))
  return(blastRes)
}