Hi, I tend in cases like this to shirk wrappers and stay closer to the command line usage of tools such as blast. Below is a generic example of run NCBI blastn (part of blast+ package) under R, and post filter the results. The approach should work with minor changes if you have not upgraded to NCBI's blast+. ~Malcolm s<- ## The sequences to be blasted and their fasta 'deflines' c('>s1','CTTTGTCTTCCCCTATGTCATCCTCAATCTCTATGAAAGCAACAC' ,'>s2','CTATGTCATCCTCAATCTCTATGAAAGCAACACCGCTACCATAGA' ) blast.f6<- ## The fields you want back from blast. c.f. `blastn -help` c('qseqid', 'sseqid', 'pident', 'qcovs') blast.out<- ## The system call to blastn system2('blastn',c('-db',"'nt'" ,'-outfmt',sprintf('"6 %s"',paste(collapse=' ',blast.f6)) , '-perc_identity',"'.90'" ,'-num_threads', 15 # use 'em if you got 'em ! ) ,input=s ,stdout=TRUE ) blast.out.df<- ## parse blast.out as a table and assign names to it `names<-`(read.table(sep='\t',textConnection(b)),blast.f6) # Query the data frame head(blast.out.df[with(b.df,pident>=90 & qcovs>95),],3) qseqid sseqid pident qcovs 1 s1 gi|380719094|gb|JQ281544.1| 100 100 2 s1 gi|161015434|gb|EU143287.2| 100 100 3 s1 gi|161015430|gb|EU143282.2| 100 100 ~ Malcolm >-----Original Message----- >From: bioconductor-bounces at r-project.org [mailto:bioconductor- bounces at r-project.org] On Behalf Of Martin Morgan >Sent: Wednesday, April 23, 2014 12:12 PM >To: prabhakar ghorpade; bioconductor at r-project.org >Subject: Re: [BioC] Blast analysis of two sequences in R > >On 04/21/2014 03:51 AM, prabhakar ghorpade wrote: >> Hello, >> I have following sequences for which I want to BLAST and see result only for sequences showing > 95% Query coverage and >>90% identity? >> >> Sequences >>> 1 >> CTTTGTCTTCCCCTATGTCATCCTCAATCTCTATGAAAGCAACAC >>> 2 >> CTATGTCATCCTCAATCTCTATGAAAGCAACACCGCTACCATAGA > >The 'annotate' package has 'blastSequences'; I'm not sure that it's useful >enough for your purposes. In the 'devel' branch (see > > http://bioconductor.org/developers/how-to/useDevel/ > >it has been updated to be more responsive and to return richer data, e.g., > > df = blastSequences("CTTTGTCTTCCCCTATGTCATCCTCAATCTCTATGAAAGCAACAC", > timeout=40, as="data.frame") > > > head(df, 1) > Hit_num Hit_id >1 1 gi|380719094|gb|JQ281544.1| > Hit_def Hit_accession Hit_len Hsp_num >1 Expression vector pAV-UCSF, complete sequence JQ281544 11534 1 > Hsp_bit-score Hsp_score Hsp_evalue Hsp_query-from Hsp_query-to Hsp_hit-from >1 82.4379 90 1.2063e-13 1 45 2126 > Hsp_hit-to Hsp_query-frame Hsp_hit-frame Hsp_identity Hsp_positive Hsp_gaps >1 2170 1 1 45 45 0 > Hsp_align-len Hsp_qseq >1 45 CTTTGTCTTCCCCTATGTCATCCTCAATCTCTATGAAAGCAACAC > Hsp_hseq >1 CTTTGTCTTCCCCTATGTCATCCTCAATCTCTATGAAAGCAACAC > Hsp_midline >1 ||||||||||||||||||||||||||||||||||||||||||||| > > > >Hope that helps; would be happy to hear of other R solutions. > >Martin > >> >> >> Can you please Suggest how can I select them in R in NCBI BLAST so that I get sequences showing > 95% Query coverage and >90% >identity. Is there program in R to select them? I want to detect number of organism showing uniques results for given sequences. >> Thanks. >> >> >> Dr. Ghorpade Prabhakar B. >> PhD Scholar ( Veterinary Biochemistry), >> IVRI, >> Izatnagar, >> Bareilly, U.P., >> India >> [[alternative HTML version deleted]] >> >> >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor at r-project.org >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor >> > > >-- >Computational Biology / Fred Hutchinson Cancer Research Center >1100 Fairview Ave. N. >PO Box 19024 Seattle, WA 98109 > >Location: Arnold Building M1 B861 >Phone: (206) 667-2793 > >_______________________________________________ >Bioconductor mailing list >Bioconductor at r-project.org >https://stat.ethz.ch/mailman/listinfo/bioconductor >Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor