> e.coli <- readAAStringSet(filepath = "uniprot-proteome_UP000000625.fasta")
> summary(e.coli)
     Length       Class        Mode 
       4306 AAStringSet          S4 
> some.peptide <- AAString("DYWRALQNRIREGHVEDVYAYRRRQ")
> summary(some.peptide)
  Length    Class     Mode 
      25 AAString       S4 
> x <- matchPDict(pdict = some.peptide, subject = e.coli)
Error in matchPDict(pdict = some.peptide, subject = e.coli) : 
  please use vmatchPDict() when 'subject' is an XStringSet object (multiple sequence)

> x <- vmatchPDict(pdict = some.peptide, subject = e.coli)
Error in .local(pdict, subject, max.mismatch, min.mismatch, with.indels,  : 
  vmatchPDict() is not ready yet, sorry

What am I doing wrong? I would like to match a list of peptide sequences (here as AAStringSet) vs. a list of Proteins (also AAStringSet).

Greetings,

Tobi

Hi Tobi,

You have only 1 peptide so why use matchPDict/vmatchPDict? Functions in the matchPDict family are meant to be used on a set of short motifs. If you only have 1 motif, use functions from the matchPattern family instead. More precisely, if you have 1 sequence in your subject, use matchPattern (but first subset your subject with e.coli[[1]] if it's an AAStringSet object). And if you have more than 1 sequence in your subject, use vmatchPattern. 

H.

I'm not aware of the notion of iterator in R in general.

If you only have 10-20k queries, you can simply use a for loop. The overhead of looping over an AAStringSet object of length 20k is not much, only 6 sec on my machine:

library(Biostrings)
library(hgu95av2probe)
peptides <- translate(DNAStringSet(hgu95av2probe)[1:20000])
system.time( for (i in seq_along(peptides)) p <- peptides[[i]] )
#   user  system elapsed 
#  6.233   0.000   6.245 

H.

Following Martin Morgans talk yesterday at the European bioconductor meeting in Basel I also implemented a nicer looking version of the above:

> foo <- function(x){vmatchPattern(pattern = x, subject = H.sapiens)}
> system.time(
+   l <- BiocGenerics::lapply(peptides, foo)
+ )
   user  system elapsed 
  4.349   0.044   4.453 

@Martin: I hope this the the "poetry" style of code you were mentioning! ;-)

...but I still have the feeling that this solution might be "correct and robust" but in no way "efficient". Another problem when writing those code was that the lapply() function to handle the S4 AAStringSet object was hidden in the BiocGenerics package and not loaded by default once you load Biostrings. What is the reason for this? It took me quite a while to find this out and it looks completely unnecessary to me. Why not put all functions that handle instances of a certain class in the package that def. the class???

> e.coli <- readAAStringSet(filepath = "uniprot-proteome_UP000000625.fasta")
> summary(e.coli)
     Length       Class        Mode
       4306 AAStringSet          S4
> peptides <- AAStringSet(c("DYWRALQNRIREGHVEDVYAYRRRQ", "QRWQSVLARDPNADGEFVFAV"))
> summary(peptides)
     Length       Class        Mode
          2 AAStringSet          S4
> x <- vmatchPDict(pdict = peptides, subject = e.coli)
Error in .local(pdict, subject, max.mismatch, min.mismatch, with.indels,  :
  vmatchPDict() is not ready yet, sorry

Finally, I want to match a set of AAStrings vs. another set of AAStrings. I simply used the trivial case of n =1 above. How to do this? If vmatchPict() is not ready yet, how can I loop over a set of AAStrings (queries) efficiently? I have 10-20k queries. The core of that loop could be:

vmatchPattern(pattern = peptides[[i]], subject = e.coli)

Is there something like an iterator for the AAStringSet class?

A downstream problem. I tried to use the index returned by vmatchPattern to the org. object and get this strange error message:

> extractAllMatches(subject = H.sapiens, mindex = l[[1]])
Error in extractAllMatches(subject = H.sapiens, mindex = l[[1]]) : 
  'subject' must be an XString or MaskedXString object
> class(H.sapiens)
[1] "AAStringSet"
attr(,"package")
[1] "Biostrings"

How can that be? T thought AAStrings is derived from XString?

A related question: If AhoCorasickTrie is soooo much faster for exact string matching why not use it under the hood when calling:

vmatchPattern(..., algorithm = "AhoCorasick-exact")

???