HI ,

I have a problem getting quality scores for the following large file.

My first step is to check the quality  to generate some stats for each fastq file, however it seems that i cannot allocate enough memory even though i have enough memory (10Gb).

> fq
class: ShortReadQ
length: 24377237 reads; width: 100 cycles

qual <- as(quality(fq), "matrix")

Error in asMethod(object) : long vectors not supported yet: memory.c:3324

How can i get around this to get the quality scores ?

Ok tried I have tried the following:

f <- FastqStreamer(file, 100000) #chunks of 100000
allqual <- list()
i=1
while (length(fq <- yield(f)) ) {
  chunkqual <- as(quality(fq), "matrix")
  allqual[[i]] <- chunkqual
  i = i+1
}
close(f)
qual <- do.call(rbind,allqual)

#quantile
qt <- apply(qual,2,quantile)

It crashes inside the loop. Maybe there is a better way to do more memory efficiently ?

> sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=ca_AD.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=ca_AD.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=ca_AD.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets
[8] methods   base     

other attached packages:
 [1] ShortRead_1.24.0        GenomicAlignments_1.2.2 Rsamtools_1.18.3       
 [4] GenomicRanges_1.18.4    GenomeInfoDb_1.2.4      Biostrings_2.34.1      
 [7] XVector_0.6.0           IRanges_2.0.1           S4Vectors_0.4.0        
[10] BiocParallel_1.0.3      BiocGenerics_0.12.1    

loaded via a namespace (and not attached):
 [1] base64enc_0.1-2     BatchJobs_1.6       BBmisc_1.9         
 [4] Biobase_2.26.0      bitops_1.0-6        brew_1.0-6         
 [7] checkmate_1.5.2     codetools_0.2-11    DBI_0.3.1          
[10] digest_0.6.8        fail_1.2            foreach_1.4.2      
[13] grid_3.1.2          hwriter_1.3.2       iterators_1.0.7    
[16] lattice_0.20-30     latticeExtra_0.6-26 RColorBrewer_1.1-2
[19] RSQLite_1.0.0       sendmailR_1.2-1     stringr_0.6.2      
[22] tools_3.1.2         zlibbioc_1.12.0    

Depending upon what you're interested in doing with the quality scores, you could try using narrow() to look at a reduced number of cycles.  This code will get you the quality matrix for the first 20 bases from each read:

qual <- as( quality( narrow( fq, start = 1, width = 20 ) ), "matrix" )

You can then summarise each column or whatever, and then repeat the process with another set of cycles.

Use FastqStreamer() to iterate (process your quality assessment in chunks, say of 1M reads per chunk) or FastqSampler() (to read a random subset), as illustrated on the help pages,e.g., ?FastqStreamer.

Thanks for the update.

What i´m trying to achieve here is ( so I might be wrong the way i do it) is to compute the quantile at each cycle across all reads. This is working pretty wall with bacterial genomes using 16s amplicons . Now working with higher datasets (metagenomics samples doing shotgun sequencing i get much more data)

So for read1 i will have the score for each cycle (i.g 1..100)

for read i will also get the score for each cycle (again 1..100)

at the end i want a matrix (basically what the as(quality(fq),"matrix") would return. With this matrix i compute the quantile per cycle and generate a boxplot that show the quality per cycle.

See attached plot. Obviously i could just use fastQC from Babraham bioinformatics but that would not be that funny

https://www.dropbox.com/s/h3v6c564gequxe3/Screen%20Shot%202015-06-02%20at%2022.46.42.png?dl=0

crashes in my case are problems allocating memory. SInce i don´t need the sequences but just the scores i thought i could just get the quality scores.

Thanks for answers.

I will follow your advice and work with 1  cycle at a time.

Best,

david