[Bioc-sig-seq] mixing 2 fastQ files

[Martin Morgan]

On 02/05/2010 03:04 AM, Ramzi TEMANNI wrote:
> Hi,
> I've found a working solution but not sure whether it's the best one, but
> could be useful for those who have similar interest :
> rfq.sf1=append(rfq.sf1[*-n.sf.ids*],rfq.phix1[*n.ph.ids*])

This might also work

subsetAssignShortReadQ <-
    function(x, i, value)
{
    qcls <- class(quality(x))
    if (qcls != class(quality(value)))
        stop("'quality' classes differ")

    sread <- sread(x)
    sread[i] <- sread(value)
    quality <- quality(quality(x))
    quality[i] <- quality(quality(value))

    ShortReadQ(sread, do.call(qcls, list(quality)))
}

example(readFastq)
from <- sample(length(rfq), 100)
to <- sample(length(rfq), 100)
subsetAssignShortReadQ(rfq, to, rfq[from])

Martin

> Regards,
> Ramzi
> 
> ----------------------------------------------------------------
> 
> 
> On Fri, Feb 5, 2010 at 11:53 AM, Ramzi TEMANNI <ramzi.temanni at gmail.com>wrote:
> 
>> Hi,
>> I'm trying to generate simulated  fastq data by mixing read from original
>> sample with Phix reads
>> So i load my sample data and phix data
>> library(ShortRead)
>> rfq.sf1<- readFastq(dir, pattern=sample.f1)
>> rfq.phix1 <- readFastq(dir, pattern=phix.f2)
>> specify the percentage of reads in the original file to be replaced by phix
>> reads and select ramdomly the reads to be replaced in sample and the reads
>> to be selected in phix
>> nbreads=fx.ratio/100*length(rfq.sf1)
>> n.sf.ids=sample(1:length(rfq.sf1),nbreads)
>> n.ph.ids=sample(1:length(rfq.phix1),nbreads)
>> but when I call
>> rfq.sf1[n.sf.ids]=rfq.phix1[n.ph.ids]
>> I get an error
>> Error in rfq.sf1[n.sf.ids] = rfq.phix1[n.ph.ids] :
>>   object of type 'S4' is not subsettable
>> I tried then
>> sread(rfq.sf1[n.sf.ids])=sread(rfq.phix1[n.ph.ids])
>> quality(rfq.sf1[n.sf.ids])=quality(rfq.phix1[n.ph.ids])
>> And i get an error too
>> Error in sread(rfq.sf1[n.sf.ids]) = sread(rfq.phix1[n.ph.ids]) :
>>   could not find function "sread<-"
>> Error in quality(rfq.sf1[n.sf.ids]) = quality(rfq.phix1[n.ph.ids]) :
>>   could not find function "quality<-"
>> what will be the method to use to replace theses reads ?
>> Thanks in advance for suggestions.
>>
>>
>> Kind regards
>> Ramzi
>> ----------------------------------------------------------------
>>
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> Bioc-sig-sequencing mailing list
> Bioc-sig-sequencing at r-project.org
> 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 M1 B861
Phone: (206) 667-2793