[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