Entering edit mode
Hi all,
I'm trying to read in a fasta sequence, extract the "gene sequences"
and
write these out to a fasta file. I can read the sequences with
read.DNAStringSet(), obtain an XStringViews object with Views(), but
I'm
having trouble knowing how to obtain the reverse complement sequence
for
the genes on the "-" strand. I can get them with a reverseComplement()
of the XStringViews object, but I can't overwrite the elements of this
object. So my solution involves dealing separately with all the genes
on
the "+" strand and those on the "-" strand. Is there an easier way?
An example:
file <- system.file("extdata", "someORF.fa", package="Biostrings")
x <- read.DNAStringSet(file, "fasta")[[1]]
names <- c("a","b","c")
starts <- c(10,384,947)
ends <- starts+20
strands <- c("+","-","+")
myViews <- Views(x, start=starts, end=ends)
names(myViews) <- names
revViews <- reverseComplement(myViews[strands=="-"])
posViews <- myViews[strands=="+"]
Many thanks,
Cei
sessionInfo()
R version 2.9.0 (2009-04-17)
i386-apple-darwin8.11.1
locale:
en_GB.UTF-8/en_GB.UTF-8/C/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] tools stats graphics grDevices datasets utils
methods
[8] base
other attached packages:
[1] Biostrings_2.12.1 IRanges_1.2.2 Biobase_2.4.1