What to do next depends on your specific research question and the purpose of your analysis. Regarding Biostrings, for learning what can be done and how to use it the best place to look would be the package web page here. Take a look at the vignettes, starting possibly with the one named "Biostrings Quick Overview".

Thank you for you reply.

However I am still struggling.

the blastSequences() does not work for me.

estimated response time 33 seconds
elapsed time 35 seconds
timeout after 61 seconds; wait another 60 seconds? [y/n]

(I have pressed "y" many times and increased the timeout time to 320 sec)

A problem which I think is related to large query sequences - as pointed at by K. Keenan. which is why he has written an alternative function blastSeqKK(), adding the argument of attempts, this one (from here: http://rstudio-pubs-static.s3.amazonaws.com/12097_1352791b169f423f910d93222a4c2d85.html):

# function definition blastSeqKK <- function (x, database = "nr", hitListSize = "10", filter = "L", expect = "10", program = "blastn", attempts = 10) { baseUrl <- "http://www.ncbi.nlm.nih.gov/blast/Blast.cgi" query <- paste("QUERY=", as.character(x), "&DATABASE=", database, "&HITLIST_SIZE=", hitListSize, "&FILTER=", filter, "&EXPECT=", expect, "&PROGRAM=", program, sep = "") url0 <- sprintf("%s?%s&CMD=Put", baseUrl, query) results <- tempfile() Sys.sleep(5) require(XML) post <- htmlTreeParse(url0, useInternalNodes = TRUE) x <- post[["string(//comment()[contains(., \"QBlastInfoBegin\")])"]] rid <- sub(".*RID = ([[:alnum:]]+).*", "\\1", x) rtoe <- as.integer(sub(".*RTOE = ([[:digit:]]+).*", "\\1", x)) url1 <- sprintf("%s?RID=%s&FORMAT_TYPE=XML&CMD=Get", baseUrl, rid) Sys.sleep(rtoe) .tryParseResult <- function(url, attempts){ for (i in 1:(attempts+1)) { result <- tryCatch({ xmlTreeParse(url, useInternalNodes=TRUE, error = xmlErrorCumulator(immediate=FALSE)) }, error=function(err) NULL) if (!is.null(result)) return(result) Sys.sleep(10) } stop(paste("no results after ", attempts, " attempts; please try again later", sep = "")) } result <- .tryParseResult(url1, attempts) qseq <- xpathApply(result, "//Hsp_qseq", xmlValue) hseq <- xpathApply(result, "//Hsp_hseq", xmlValue) require(Biostrings) res <- list() for (i in seq_len(length(qseq))) { res[i] <- DNAMultipleAlignment(c(hseq[[i]], qseq[[i]]), rowmask = as(IRanges(), "NormalIRanges"), colmask = as(IRanges(), "NormalIRanges")) } res }

However, as I understand this there is no longer the choice/argument of getting the output as a "data.fame", like in the original function:

blastSequences(x, database, hitListSize, filter, expect, program,
      timeout=40, as=c("DNAMultipleAlignment", "data.frame", "XML"))

So this is where I am stuck...

I tried to simply (naïve) add the argument:

# function definition
blastSeqKK <- function (x, database = "nr", hitListSize = "10", 
                        filter = "L", expect = "10", program = "blastn",timeout=40,
 as=c("DNAMultipleAlignment", "data.frame", "XML"),
                        attempts = 10) 

But it was "ignored" and the output when running:

res <- lapply(mySeq, function(x){
  # collapse seq into string
  seqCollapse <- paste(toupper(as.character(x)), collapse = "")
  # run blast
  blastRes <- blastSeqKK(x = seqCollapse, hitListSize = 2, timeout=60, 
    attempts = 30, as=c("data.frame") )
  return(blastRes)

Was still this:

> head(res)
[[1]]
[[1]][[1]]
DNAMultipleAlignment with 2 rows and 324 columns
     aln
[1] CTGAATCTGCGAACGGCTCCGCACACCAGTTGCAAA...GATGGCAGCGTAAGGGACATGCTATGGTTACAACGG
[2] CTGAATCTGCGTACGGCTCCGCACACCAGTTGCAAA...GATGGCAGCGTAAGGGACATGCTATGG-TACAACGG

[[1]][[2]]
DNAMultipleAlignment with 2 rows and 324 columns
     aln
[1] CTGAATCTGCGAACGGCTCCGCACACCAGTTGCAAA...GATGGCAGCGTAAGGGACATGCTATGGTTACAACGG
[2] CTGAATCTGCGTACGGCTCCGCACACCAGTTGCAAA...GATGGCAGCGTAAGGGACATGCTATGG-TACAACGG

[[2]]
[[2]][[1]]
DNAMultipleAlignment with 2 rows and 307 columns
     aln

Sooo, I don't know how to proceed right now, so again any suggestions would be much appreciated.

/Mathilde

Thank you so much for the detailed answer.

Your first advice worked, and I can get the "data.fame" output - perfect :)

However the function blastSeqKK still only works when the fasta file contains only 2 sequences - more, and it gives the error of "no result after xx attempts" I have set the attempts up to 50.

I also tried your 2nd advice and let it run for 4h before I stopped it - this was a file containing only 10 sequences.
I think I will try to leave it overnight, just to see if it eventually comes to and end - However I don't understand why it is struggling so? - A manual blast in NCBI's webpage takes only a few sec for a file with all (42) seq. Hhm...