Hi Everyone at Bioconductor,

I've been doing work with large multiple sequence alignments in which I find the polymorphic sites, keeping them and their base positions, and rejecting all the non-polymorphic sites. I usually use a package like ape for this - which implements the DNAbin class. This time the MSA I have is simply too large and can't be represented as a DNAbin class and R gives an error about not supporting Long vectors of Raw. I've heard and read that the sequence types in Bioconductor are good because they don't store the entire sequence in memory and are built to handle very large sequences efficiently. This will be my first use of Bioconductor. I've currently identified the polymorphic sites using the following code:

library(Biostrings)
origMAlign <- readDNAStringSet(filepath = "~/Dropbox/MySequences.fas", format = "fasta")
consensusM <- consensusMatrix(origMAlign)[1:4,]
polymorphic <- which(colSums(consensusM != 0) > 1)
for(i in 1:length(origMAlign))
  origMAlign[[i]] <- origMAlign[[i]][polymorphic]

Is this the most efficient way to do this with the Biostring classes - being R, use of a loop feels instinctively wrong! 

Thanks,

Ben Ward - UEA & TSL.

Hi Ben,

A much faster way is to use an IntegerList subscript to subset origMAlign:

i <- rep.int(IntegerList(polymorphic), length(origMAlign))
origMAlign[i]

In the S4Vectors/IRanges/GenomicRanges/Biostrings infrastructure (and for S4 objects that rely on that infrastructure) we've extended the capabilities of the [ operator so you can:

1. Use an Rle (logical-, integer-, numeric-, character-, or factor-Rle) or Ranges subscript to subset a Vector derivative.
2. Use a list-like object (i.e. list or List) to subset a List derivative 'x'. In that case, the i-th list element in the subscript is used to subset the i-th list element in 'x'.

Hope this helps,

H.