Martin,
  I was unaware of the trimEnds() function, it worked wonderfully without
having to use your modified trimEnds().  I imagine it is because the
cluster has a lot of memory available.  I ran it a second time using the
modification, but it was much slower.

Thanks again!


On Fri, May 17, 2013 at 12:20 PM, Martin Morgan <mtmorgan@fhcrc.org> wrote:

> On 05/17/2013 09:19 AM, Sam McInturf wrote:
>
>> Hello everyone,
>>    I am trying to do some quality trimming on some RNA seq reads
>> (Illumina,
>> 100 bp, single end).  I have been following the pipeline from the UCR
>> manuals<http://manuals.**bioinformatics.ucr.edu/home/**
>> ht-seq#TOC-Trimming-Low-**Quality-3-Tails-from-R<http://manuals.bioinformatics.ucr.edu/home/ht-seq#TOC-Trimming-Low-Quality-3-Tails-from-R>
>> >
>>
>> but
>> I have run into an error I can't seem to resolve.
>>
>> library(ShortRead)
>> reads <- readFastq("myFile.fastq")
>>
>> qualityCutoff <- 10 # remove read tails with quality lower than this
>>
>> seqs <- sread(reads) # get sequence list
>>
>> qual <- PhredQuality(quality(quality(**reads))) # get quality score list
>> as
>> PhredQuality
>>
>>  length(qual)                 #my length is positive
>>>
>> [1] 39145018
>>
>>  myqual_mat <- matrix(charToRaw(as.character(**unlist(qual))),
>>>
>> nrow=length(qual), byrow=TRUE)
>>
>
> I'm guessing that what's happening is that as.character(unlist(qual)) is
> making a very large vector sum(width(qual)), larger than can be represented
> in (your) version of R (the output of sessionInfo() would be helpful...).
>
> But I don't think you need to go this route;  have you looked at
> ShortRead::trimEnds & friends? For instance, after running
>
>   exmample(trimEnds)
>
> I have an object 'rfq' with some qualities
>
> > quality(rfq)
> class: SFastqQuality
> quality:
>   A BStringSet instance of length 256
>       width seq
>   [1]    36 ]]]]]]]]]]]]Y]Y]]]]]]]]]]]]**VCHVMPLAS
>   [2]    36 ]]]]]]]]]]]]Y]]]]]]]]]YPV]T][**PZPICCK
>   [3]    36 ]]]]Y]]]]]V]T]]]]]T]]]]]V]**TMJEUXEFLA
>   [4]    36 ]]]]]]]]]]]]]]]]]]]]]]T]]]]**RJRZTQLOA
>   [5]    36 ]]]]]]]]]]]]]]]]]T]]]]]]]]]]**MJUJVLSS
>   ...   ... ...
> [252]    36 ]]]]]]]]]]]Y]]]]]NY]]]]Y]VR]**MJQNSAOC
> [253]    36 ]]]]]]]]]]]]]]]]]]]]]]]Y]]]**VTVVRVMSM
> [254]    36 ]]]]Y]Y]]]]]]]OYY]]]Y]]]**YYVVTSZUOOHH
> [255]    36 ]]]]]]]]]]]]]]]T]]]]]]]]]V]T[**OVXEJSJ
> [256]    36 ]]]]]]]]]]]Y]]T]T]]]]**TRYVMEVVRSRHHNH
>
> And I can simultaneously trim reads and sequences with
>
>     trimEnds(rfq, "V")
>
> which you can see in action as
>
> > quality(trimEnds(rfq, "V"))
> class: SFastqQuality
> quality:
>   A BStringSet instance of length 256
>       width seq
>   [1]    27 ]]]]]]]]]]]]Y]Y]]]]]]]]]]]]
>   [2]    31 ]]]]]]]]]]]]Y]]]]]]]]]YPV]T][**PZ
>   [3]    32 ]]]]Y]]]]]V]T]]]]]T]]]]]V]**TMJEUX
>   [4]    31 ]]]]]]]]]]]]]]]]]]]]]]T]]]]**RJRZ
>   [5]    28 ]]]]]]]]]]]]]]]]]T]]]]]]]]]]
>   ...   ... ...
> [252]    28 ]]]]]]]]]]]Y]]]]]NY]]]]Y]VR]
> [253]    27 ]]]]]]]]]]]]]]]]]]]]]]]Y]]]
> [254]    31 ]]]]Y]Y]]]]]]]OYY]]]Y]]]**YYVVTSZ
> [255]    32 ]]]]]]]]]]]]]]]T]]]]]]]]]V]T[**OVX
> [256]    24 ]]]]]]]]]]]Y]]T]T]]]]TRY
>
> The following might be helpful to transform a large file --
>
> setMethod(trimEnds, "character",
>     function(object, a, left = TRUE, right = TRUE,
>              relation = c("<=", "=="), ...,
>              destination, yieldSize=1000000, ranges = FALSE)
> {
>     if (missing('destination'))
>         stop("'destination' required")
>     strm <- FastqStreamer(object, yieldSize)
>     tot <- nNuc <- nTrimNuc <- 0L
>     while (length(fq <- yield(strm))) {
>         tot <- tot + length(fq)
>         nNuc <- nNuc + sum(width(fq))
>         fq <- trimEnds(fq, a, left, right, relation, ..., ranges=ranges)
>         nTrimNuc <- nTrimNuc + sum(width(fq))
>         writeFastq(fq, destination, "a")
>     }
>     list(destination=destination,
>          c(TotalReads=tot, TotalNucleotides=nNuc,
>            TrimmedNucleotides=nNuc - nTrimNuc))
> })
>
> provide the first argument to trimEnds as a simple character vector, and
> provide a 'destination' file name. For example
>
>   ## just a convenient path to a fastq file
>   fl <- file.path(analysisPath(sp), "s_1_sequence.txt")
>   trimEnds(fl, "V", destination=tempfile())
>
> A variation of this will be added to the 'devel' version of ShortRead, so
> if there are suggestions for improvements please let me know.
>
> Martin
>
>
>
>
>  Error in .Call(.NAME, ..., PACKAGE = PACKAGE) :
>>    negative length vectors are not allowed
>>
>> as(qual, matrix) gives the same error, and as.matrix(qual) does as well,
>> not sure how those two commands are different, but tried anyway.
>>
>> If I run the exact same commands, instead of the full 39145018 reads, I
>> use
>> the first 100, it works like a charm.  I am working on a linux cluster, if
>> that is important
>>
>> Does anyone have an idea?
>>
>> Cheers!
>>
>>
>
> --
> Computational Biology / Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N.
> PO Box 19024 Seattle, WA 98109
>
> Location: Arnold Building M1 B861
> Phone: (206) 667-2793
>



-- 
Sam McInturf

	[[alternative HTML version deleted]]

