hi Martin, thanks for your quick answer. I can't get round to the fact that using srdistance() will remove the whole read, instead of just editing the adapter off and keeping the remaining fragment of the read, which it is useful information. That is the difference I see between using srdistance() and the approach using pairwiseAlignment(). The second, as it is shown in the Biostrings vignette using DNAStringSets, will only edit the adapter sequence from the read. I may be wrong but the problem I find is that I cannot work back from a DNAStringSet to a ShortReadQ object that will include quality calls for the remaining bases of the edited reads, and thus ending up with a full clean ShortReadQ object that can be transformed back to fastq format.
Regarding appending two ShortReadQ objects, the idea came up early in this thread when I run out of computer memory and I was suggested to split the original data so that I could do the adapter search in small pieces. I would then have to put them back together. Also, another situation where needed is merging two flowcell lanes (or more) coming from the same sample when you are analyzing large genomes and do not have enough coverage on one single lane. I guess this could be done outside R using cat. Dave > To: [email protected] > CC: [email protected]; [email protected] > Subject: Re: [Bioc-sig-seq] adapter removal > From: [email protected] > Date: Wed, 14 Jan 2009 06:03:26 -0800 > > Hi Dave -- > > "David A.G" <[email protected]> writes: > > > Hi Martin and list, > > going back to the adapter topic, I tried the srdistance() function and it > > works really quick! > > yeah it's great; Herve Pages & Patrick Aboyoun deserve credit for this. > > > I have a question regarding the edit distance. After having used tables() I > > find that the second top element is CTGGACTTTGTAGGATACCCTCGCTTTCCTGCTC with > > 890 occurrences. The edit distance for this sequence is of course 0. But if > > I > > try only part of it as an adapter, CTGGACTTTGTAGGATA, the lowest edit > > distance > > is 17 and there are 14758 reads with this sequence. > > Would you just use 17 as the distance cutoff? > > If I wanted to find just those 890 occurrences in an object 'fq' with > the adapter at the start, I'd use srdistance on a DNAStringSet created > from just those nucleotides, along the lines of > > > d <- srdistance(DNAStringSet(sread(fq), 1, 17), "CTGGACTTTGTAGGATA") > > fg[d != 0] > > It's not really clear that the adapters will always be exactly at the > beginning, etc., so this might just be a starting point. Creating a > new DNAStringSet is not particularly memory- or CPU-intensive, so > creating it 'on the fly' as in the above example should be > reasonable. I think Patrick is working on a one-stop-shopping adapter > removal function in Biostrings. > > > Also,is it possible to merge two ShortReadQ objects, or, having subset an > > ShortReadQ object in two by indexing, how do I put them back together? > > The ability to 'append' two ShortReadQ objects was added just > yesterday > > > append(fq, fq) > > It should be available in the 'devel' branch using biocLite after > about 1pm Seattle time, all being well; look for ShortRead version >= > 1.1.33. It's possible to do this with the released version of > ShortRead, and I can provide off-line hints if you'd like, but I'd > recommend working with the development version of R and these > packages. > > I'm curious to know the use case here; the need to append hasn't come > up in my own work. > > Martin > > > Cheers, > > Dave > >> To: [email protected] > >> CC: [email protected]; [email protected] > >> Subject: Re: [Bioc-sig-seq] adapter removal > >> From: [email protected] > >> Date: Fri, 9 Jan 2009 19:19:18 -0800 > >> > >> "David A.G" <[email protected]> writes: > >> > >> > Jim, thanks for pointing that out, I hadn´t faced that issue yet. The > > "reads" object I created only had the sequences and not the qualities, and > > after editing the adapters it is interesting to only have the quality calls > > for the remaining bases. > >> > > >> > Tyler and Martin, thanks for the information on writing back to fastq. > >> > >> There is a now an initial version of writeFastq in the 'devel' version > >> of ShortRead. It is memory and time efficient; try > >> > >> > example(writeFastq) > >> > ?writeFastq > >> > >> please let me know if there are problems, or other missing parts of > >> ShortRead. > >> > >> Martin > >> > >> > Dave > >> > > >> >> Date: Thu, 8 Jan 2009 15:12:01 -0800 > >> >> From: [email protected] > >> >> To: [email protected] > >> >> Subject: Re: [Bioc-sig-seq] adapter removal > >> >> > >> >> Here's another "quick and dirty" ShortReadQ-fastq exporter. It's very > > slow... > >> >> > >> >> exportFastaQ <- function(shortReadQObject, fileName) { > >> >> ## dump ShortReadQ object data to FASTQ format > >> >> entries <- length(shortReadQObject) > >> >> i <- 1 > >> >> while (i <= entries){ > >> >> cat(append("@", as.character(id(shortReadQObject)[[i]])), sep="", > >> >> file=fileName, append=TRUE) > >> >> cat("\n", file=fileName, append=TRUE) > >> >> cat(as.character(sread(shortReadQObject)[i]), sep="", > >> >> file=fileName, append=TRUE) > >> >> cat("\n", file=fileName, append=TRUE) > >> >> cat(append("+", as.character(id(shortReadQObject)[[i]])), sep="", > >> >> file=fileName, append=TRUE) > >> >> cat("\n", file=fileName, append=TRUE) > >> >> cat(as.character(quality(shortReadQObject)[[i]]), sep="", > >> >> file=fileName, append=TRUE) > >> >> cat("\n", file=fileName, append=TRUE) > >> >> i <- i + 1 > >> >> } > >> >> } > >> >> > >> >> Sincerely, > >> >> Tyler William H Backman > >> >> Bioinformatics Analyst > >> >> Institute for Integrative Genome Biology (IIGB) > >> >> Department of Botany and Plant Sciences > >> >> E-mail: [email protected] > >> >> 1007 Noel T. Keen Hall > >> >> University of California > >> >> Riverside, CA 92521 > >> >> > >> >> > Hi Jim -- > >> >> > > >> >> > "James W. MacDonald" <[email protected]> writes: > >> >> > > >> >> >> As a follow-up question, once you have removed the adapter sequences > >> >> >> is it then possible to write the ShortReadQ object back out in a > >> >> >> fastq > >> >> >> format for aligning with e.g., bowtie? I know I can use > >> >> >> write.XstringSet() to output in fasta format, but I would like to > >> >> >> feed > >> >> >> bowtie the quality scores as well. > >> >> > > >> >> > Actually a 'writeFastq' is on my 'do this week' list. In the mean time > >> >> > my best (untested) suggestion is along the lines of > >> >> > > >> >> > n <- length(fq) > >> >> > lines <- character(n) > >> >> > lines[(seq_len(n)) * 4 - 3] <- > >> >> > paste("@", as.character(id(fq)), sep="") > >> >> > lines[(seq_len(n)) * 4 - 2] <- > >> >> > as.character(sread(rfq1)) > >> >> > lines[(seq_len(n)) * 4 - 1] <- > >> >> > paste("+", as.character(id(fq)), sep="") > >> >> > lines[(seq_len(n)) * 4] <- > >> >> > as.character(quality(quality(fq))) > >> >> > writeLines(lines, "my.fastq") > >> >> > > >> >> > This will be unnecessarily memory intensive. > >> >> > > >> >> >> Or is the speed reasonably comparable to use the PDict functions in > >> >> >> Biostrings to do the aligning? > >> >> > > >> >> > Bowtie is very fast. One potential 'gotcha' is the orientation of > >> >> > reads throughout the workflow -- Eland _sequence.txt files are the > >> >> > reads as they come from the machine; MAQ (and I think Bowtie) aligned > >> >> > files report reads that align to the '-' strand in their reverse > >> >> > complement sequence, with qualities reversed to match the reverse > >> >> > compelement read, and with the alignment position recorded at the > >> >> > 'left-most' (rather than 5', say) end. Not that it will be a problem > >> >> > in the current use case, just something to be aware of as you move > >> >> > between tools. > >> >> > > >> >> > Martin > >> >> > > >> >> >> Best, > >> >> >> > >> >> >> Jim > >> >> >> > >> >> >> > >> >> >> > >> >> >> Martin Morgan wrote: > >> >> >>> "David A.G" <[email protected]> writes: > >> >> >>> > >> >> >>>> Dear list, > >> >> >>>> > >> >> >>>> I have some experience with Bioconductor but am newbie to this list > >> >> >>>> and to NGS. I am trying to remove some adapters from my solexa > >> >> >>>> s_N_sequence.txt file using Biostrings and ShortRead packages and > > the > >> >> >>>> vignettes. I managed to read in the text file and got to save the > >> >> >>>> reads as follows > >> >> >>>> > >> >> >>>> fqpattern <- "s_4_sequence.txt" > >> >> >>>> f4 <- file.path(analysisPath(sp), fqpattern) > >> >> >>>> fq4 <- readFastq(sp, fqpattern) > >> >> >>>> reads <- sread(fq4) #"reads" contains more than 4 million 34-length > >> >> >>>> fragments > >> >> >>>> > >> >> >>>> Having the following adapter sequence: > >> >> >>>> > >> >> >>>> adapter <- DNAString("ACGGATTGTTCAGT") > >> >> >>>> > >> >> >>>> I tried to mimic the example in the Biostring vignette as follows: > >> >> >>>> > >> >> >>>> > >> >> >>>> myAdapterAligns <- pairwiseAlignment(reads, adapter, type = > > "overlap") > >> >> >>>> > >> >> >>>> but after more than two hours the process is still running. > >> >> >>>> > >> >> >>>> I am running R 2.8.0 on a 64bit linux machine (Kubuntu 2.6.24) with > >> >> >>>> 4Gb RAM, and I only have some 30Mb free RAM left. I found a thread > > on > >> >> >>>> adapter removal but does not clear things much to me, since as far > > as > >> >> >>>> I understood the option mentioned in the thread is not appropriate > >> >> >>>> (quote :(though apparently this is not entirely satisfactory, see > > the > >> >> >>>> second entry!)). > >> >> >>>> > >> >> >>>> Is this just a memory issue or am I doing something wrong? Shall I > >> >> >>>> leave the process to run for longer? > >> >> >>> As a bit of a follow-up, a reasonable strategy is to figure out how > >> >> >>> long / how much memory the calculation takes on smaller chunks of > >> >> >>> data. > >> >> >>> Also, here I 'clean' (remove reads with 'N's) and then calculate the > >> >> >>> srdistance to your adapter on a lane with 5978132 reads in about a > >> >> >>> minute > >> >> >>> > >> >> >>>> fq <- readFastq(sp, "s_1_sequence.txt") > >> >> >>>> fq > >> >> >>> class: ShortReadQ > >> >> >>> length: 5978132 reads; width: 53 cycles > >> >> >>>> system.time(dist <- srdistance(clean(fq), "ACGGATTGTTCAGT")) > >> >> >>> user system elapsed 60.960 0.000 60.961 I might then > >> >> >>> tabulate the distances > >> >> >>> > >> >> >>>> table(dist[[1]]) > >> >> >>> and subset fq based on some criterion > >> >> >>> fq[dist[[1]] > 4] > >> >> >>> There's a second interface to this, too, using filters, e.g., > >> >> >>> fq[srdistanceFilter("ACGGATTGTTCAGT", 4)(fq)] > >> >> >>> Martin > >> >> >>> > >> >> >>>> TIA for your help, > >> >> >>>> > >> >> >>>> Dave > >> >> >>>> > >> >> >>>> _________________________________________________________________ > >> >> >>>> Show them the way! Add maps and directions to your party invites. > >> >> >>>> > >> >> >>>> [[alternative HTML version deleted]] > >> >> >>>> > >> >> >>>> _______________________________________________ > >> >> >>>> Bioc-sig-sequencing mailing list > >> >> >>>> [email protected] > >> >> >>>> https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing > >> >> >>> > >> >> >> > >> >> >> -- > >> >> >> James W. MacDonald, M.S. > >> >> >> Biostatistician > >> >> >> Hildebrandt Lab > >> >> >> 8220D MSRB III > >> >> >> 1150 W. Medical Center Drive > >> >> >> Ann Arbor MI 48109-5646 > >> >> >> 734-936-8662 > >> >> > > >> >> > -- > >> >> > Martin Morgan > >> >> > Computational Biology / Fred Hutchinson Cancer Research Center > >> >> > 1100 Fairview Ave. N. > >> >> > PO Box 19024 Seattle, WA 98109 > >> >> > > >> >> > Location: Arnold Building M2 B169 > >> >> > Phone: (206) 667-2793 > >> >> > > >> >> > _______________________________________________ > >> >> > Bioc-sig-sequencing mailing list > >> >> > [email protected] > >> >> > https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing > >> >> > > >> >> > > >> >> > >> >> _______________________________________________ > >> >> Bioc-sig-sequencing mailing list > >> >> [email protected] > >> >> https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing > >> > > >> > _________________________________________________________________ > >> > Show them the way! Add maps and directions to your party invites. > >> > > >> > [[alternative HTML version deleted]] > >> > > >> > _______________________________________________ > >> > Bioc-sig-sequencing mailing list > >> > [email protected] > >> > https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing > >> > >> -- > >> Martin Morgan > >> Computational Biology / Fred Hutchinson Cancer Research Center > >> 1100 Fairview Ave. N. > >> PO Box 19024 Seattle, WA 98109 > >> > >> Location: Arnold Building M2 B169 > >> Phone: (206) 667-2793 > > > > ------------------------------------------------------------------------------ > > > > Get news, entertainment and everything you care about at Live.com. [[Check > > it > > out!]] > > -- > Martin Morgan > Computational Biology / Fred Hutchinson Cancer Research Center > 1100 Fairview Ave. N. > PO Box 19024 Seattle, WA 98109 > > Location: Arnold Building M2 B169 > Phone: (206) 667-2793 _________________________________________________________________ [[alternative HTML version deleted]]
_______________________________________________ Bioc-sig-sequencing mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
