Hello,

This might be a daft question but I'd be interested to know which R package to use to blast query a DNA sequence of around 500bp against the "blastn" database at NCBI. I've got a data frame of antigen sequences in R that instead of copy/pasting each into the blast website, I'd like to query through R.

There are about 5-6 different R packages that offer a blast sequence function of some type, e.g.

blast{BoSSA}

blastSequences{annotate}

blastSeq{hoardeR}

blast(rBLAST)

Then there's the rentrez package hosted by rOpenSci.org that allows searching of Entrez databases, BLAST being one of them.

I've tried them all and for some reason blastSequences{annotate} always seems to time out, blast{BoSSA} needs the input to be a DNAbin class and just returns NAs in the output, blastSeq{hoardeR} requests an email address as input and an XML as output and blast(rBLAST) seems to be an R wrapper for running BLAST+ locally.

What are other people's experiences with running a blast query from R? I imagine this would be quite a common task but there's not much information out there. Honestly, I would just like something akin to the qblast function from the Bio.Blast.NCBIWWW Biopython module.

Many thanks,

Miha

I recommend not using an R package.  Rather, make system calls to execute blast jobs remotely at NCBI, and get tabular results back in an R data.frame

Here is an example:

seqs <-
  ## The sequences to be blasted and their fasta 'deflines' as an character vector.  
  c('>s1','CTTTGTCTTCCCCTATGTCATCCTCAATCTCTATGAAAGCAACAC'
    ,'>s2','CTATGTCATCCTCAATCTCTATGAAAGCAACACCGCTACCATAGA')

blast.f6 <-
  ## The fields you want back from blast.  c.f. `blastn -help` for a full list of available fields.
  c('qseqid', 'sseqid', 'pident', 'qcovs')

blast.out <-
  ## Run the job, converting output to an R data.frame
  system2('blastn',c('-db',"'nt'"
                    ,'-outfmt',sprintf('"6 %s"',paste(collapse=' ',blast.f6))
                    ,'-perc_identity',"'.90'"
                    ,'-entrez_query','"Aichi virus 1[ORGANISM]"'
                    ,'-remote'
                     )
         ,input=seqs
         ,stdout=TRUE                   # capture the output in a character vector
          )

blast.out.df <-
  ## parse blast.out as a table and assign names to it
  `names<-`(read.table(quote=""
                      ,sep='\t'
                      ,textConnection(blast.out)
                       )
           ,blast.f6)