Which package to run blast
1
0
Entering edit mode
biomiha ▴ 20
@biomiha-11346
Last seen 4 months ago
UK/Cambridge

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

blast • 8.9k views
ADD COMMENT
1
Entering edit mode

I think annotate is the only one of these packages that is a Bioconductor package? For me example(blastSequences) performs several queries that complete without timing out; does it help to increase the timeout argument to something larger? Can you share the specific example of what you are trying to query? As far as I can tell from the blastSequence code, this is entirely on the latency of the blast server.

ADD REPLY
0
Entering edit mode

Thanks Martin. For some reason it was constantly timing out yesterday but it seems perfectly fine today. Definitely the best package in my opinion.

ADD REPLY
0
Entering edit mode

I'm trying to use blastSequences but I always get the same error

>example(blastSequences)

blstSq> ## x can be an entrez gene ID
blstSq> blastSequences(17702, timeout=40, as="data.frame")
Error: failed to load external entity "http://www.ncbi.nlm.nih.gov/blast/Blast.cgi?QUERY=17702&DATABASE=nr&HITLIST_SIZE=10&FILTER=L&EXPECT=10&PROGRAM=blastn&CMD=Put"

Any idea on what is missing?

Thanks,

Xavier

ADD REPLY
0
Entering edit mode

Are you using a current version of the package? Please post the output of the command sessionInfo(). NCBI is expecting an https:// url

ADD REPLY
1
Entering edit mode
Malcolm Cook ★ 1.6k
@malcolm-cook-6293
Last seen 3 months ago
United States

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)
ADD COMMENT
0
Entering edit mode

Thanks Malcolm. I'll certainly have a look at your solution if I want to run a local db blast but I just wanted an interface with the NCBI hosted databases.

ADD REPLY
1
Entering edit mode

I have improved the example to:

  • provide the -remote option (to get job to run at NCBI on their database)
  • omit data.table example (which was extraneous and incorrect)
  • demonstrate use of -entrez_query option, which hints at the power of filtering results using entrez queries.
ADD REPLY
0
Entering edit mode

Hi Malcolm,

This still requires installing a local instance of the blast program, which I think is not necessary if one only wants to use the NCBI databases.

Thanks,

Miha

 

ADD REPLY

Login before adding your answer.

Traffic: 549 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6