How to select names with wildcard from AAStringSet (biostrings)

Question

[Benn]

Maybe a simple question, hopefully I get a simple answer.

I have imported a fasta file of proteins in R with biostrings.

> mySeqs <- readAAStringSet(fastaFile)   # from package Biostrings
> mySeqs
  A AAStringSet instance of length 318
      width seq                                             names               
  [1]    98 QVQLVQSGAEVKKPGASVKVSC...STAYMELRSLRSDDTAVYYCAR IGHV1-18*01
  [2]    92 QVQLVQSGAEVKKPGASVKVSC...TTDTSTSTAYMELRSLRSDDTA IGHV1-18*02
  [3]    98 QVQLVQSGAEVKKPGASVKVSC...STAYMELRSLRSDDMAVYYCAR IGHV1-18*03
  [4]    98 QVQLVQSGAEVKKPGASVKVSC...STAYMELRSLRSDDTAVYYCAR IGHV1-18*04
  [5]    98 QVQLVQSGAEVKKPGASVKVSC...STAYMELSRLRSDDTVVYYCAR IGHV1-2*01
  ...   ... ...
[314]    91 QVQLVQSGSELKKPGASVKVSC...FSLDTSVSTAYLQISTLKAEDT IGHV7-4-1*03
[315]    98 QVQLVQSGSELKKPGASVKVSC...SMAYLQISSLKAEDTAVYYCAR IGHV7-4-1*04
[316]    98 QVQLVQSGSELKKPGASVKVSC...SMAYLQISSLKAEDTAVCYCAR IGHV7-4-1*05
[317]    98 QVQLVQSGHEVKQPGASVKVSC...STAYLQISSLKAEDMAMYYCAR IGHV7-81*01
[318]    98 EAQLTESGGDLVHLEGPLRLSC...YMLYMQMISLRTQNMAAFNCAG IGHV8-51-1*02

So how can I select only sequences with a certain name, lets say all names that contain "IGHV1". I have tried some grep, but that does not work on these objects.

Thanks in advance.

Answer 1

[AK]

Try:

mySeqs[grepl("IGHV1", mySeqs@ranges@NAMES)]

Example:

> mySeqs
  A AAStringSet instance of length 10
     width seq                                                                                names               
 [1]    44 QVQLVQSGAEVKKPGASVKVSCSTAYMELRSLRSDDTAVYYCAR                                       IGHV1-18*01
 [2]    44 QVQLVQSGAEVKKPGASVKVSCTTDTSTSTAYMELRSLRSDDTA                                       IGHV1-18*02
 [3]    44 QVQLVQSGAEVKKPGASVKVSCSTAYMELRSLRSDDMAVYYCAR                                       IGHV1-18*03
 [4]    44 QVQLVQSGAEVKKPGASVKVSCSTAYMELRSLRSDDTAVYYCAR                                       IGHV1-18*04
 [5]    44 QVQLVQSGAEVKKPGASVKVSCSTAYMELSRLRSDDTVVYYCAR                                       IGHV1-2*01
 [6]    44 QVQLVQSGSELKKPGASVKVSCFSLDTSVSTAYLQISTLKAEDT                                       IGHV7-4-1*03
 [7]    44 QVQLVQSGSELKKPGASVKVSCSMAYLQISSLKAEDTAVYYCAR                                       IGHV7-4-1*04
 [8]    44 QVQLVQSGSELKKPGASVKVSCSMAYLQISSLKAEDTAVCYCAR                                       IGHV7-4-1*05
 [9]    44 QVQLVQSGHEVKQPGASVKVSCSTAYLQISSLKAEDMAMYYCAR                                       IGHV7-81*01
[10]    44 EAQLTESGGDLVHLEGPLRLSCYMLYMQMISLRTQNMAAFNCAG                                       IGHV8-51-1*02
> mySeqs[grepl("IGHV1", mySeqs@ranges@NAMES)]
  A AAStringSet instance of length 5
    width seq                                                                                 names               
[1]    44 QVQLVQSGAEVKKPGASVKVSCSTAYMELRSLRSDDTAVYYCAR                                        IGHV1-18*01
[2]    44 QVQLVQSGAEVKKPGASVKVSCTTDTSTSTAYMELRSLRSDDTA                                        IGHV1-18*02
[3]    44 QVQLVQSGAEVKKPGASVKVSCSTAYMELRSLRSDDMAVYYCAR                                        IGHV1-18*03
[4]    44 QVQLVQSGAEVKKPGASVKVSCSTAYMELRSLRSDDTAVYYCAR                                        IGHV1-18*04
[5]    44 QVQLVQSGAEVKKPGASVKVSCSTAYMELSRLRSDDTVVYYCAR                                        IGHV1-2*01
> 

[ATpoint]

Simpler than mine, +1. I always forget than Biostrings follows the same logic as IRanges/GenomicRanges when it comes to subsetting. You almost never need any apply based method :)

[Benn]

Great solution, thanks!

Answer 2

[ATpoint]

With a dummy dataset:

require(Biostrings)

example.dna <- DNAStringSet(c(`IGHV1-18*01`="ACGGTTAGT", 
                              `IGHV1-18*02`="GTCGGGATA", 
                              `IGHV2-18*03`="GGGTCTAAA"))

example.aa <- translate(example.dna)

## Make a named list of AAStrings, then combine to AAStringSet:
out.aa <- AAStringSet(x = sapply(grep("IGHV1", names(example.aa)), 
                                 function(x) {
                                   tmp <- list(example.aa[[x]])
                                   names(tmp) <- names(example.aa)[x]
                                   return(tmp)
                                  })
)

[Benn]

I'll go for the one-liner, but thanks anyways!