[BioC] shortread quality

David martin vilanew at gmail.com
Wed Jun 6 17:00:07 CEST 2012


Hi,
I'm reading a fastq file from the solexa sequencer.
I would like to know how many reads have a phred score (>Q29).  The 
thing is that i get the densities so i don't really know how many reads 
from the total pass that filter. It's probaly easy for you so any hint 
would be helpful

library("ShortRead")
fqpattern <- "1102sdd_SN148_A_s_3_seq_GJH-85.txt"

path = getwd()
sp <- SolexaPath(path,dataPath=path,analysisPath=path)

# Read fastq File and save report
fq <- readFastq(sp, fqpattern)
qaSummary <- qa(fq,fqpattern)
save(qaSummary, file=file.path("./", paste(fqpattern,".rda",sep="" )))
report(qaSummary,dest="report")

#Quality

idx = which(qaSummary[["readQualityScore"]]["quality"] > 29)
a = cbind( qaSummary[["readQualityScore"]][idx,"quality"] , 
qaSummary[["readQualityScore"]][idx,"density"])
a #reads with a quality >Q29

#How to get the total number ? or percent compared to the total number 
of reads ?

thanks



More information about the Bioconductor mailing list