Hello,
I'm hoping that someone can offer guidance about filtering reads from
fastQ files based on quality score.
Using Shortread (or another program in Bioconductor), is it possible
to filter a fastQ file to remove reads that have an average quality
score below a certain value, and then output the remaining reads as a
fastQ file? If so, could one apply this filter to paired-end data?
For example, it would be helpful to remove a read pair when either or
both individual reads in that pair fall below the quality threshold.
I see the defined filters in Shortread (strandFilter, nFilter, etc),
but I'm not sure how to create a custom filter based on average
quality score. My modest start on the problem and a few lines from
our fastQ files are below.
Any hints would be very much appreciated.
Jen
library(ShortRead)
reads <- readFastq("~/", pattern="query.fastq")
# import reads from fastq file into ShortReadQ Object
# replace "~/" with your path and query.fastq with filename
#### ?? how to filter out (remove) reads if average quality score <25?
(for example) ####
#### ?? how to apply this filter to paired-end data ####
writeFastq(filteredReads, "output.fastq")
# export filtered reads to fastq file
A couple of sample lines from the two fastq files... We have about
7.5 million reads per file. ("ANNACG" notes the adaptor for this
library, which oddly had a consistent N in the sequence.)
paired end file 1
@HWW-T7400-009_1_120_1879_1134#ANNACG/1
TGTTGTAACTCTTTAATTACCCACNGAGNNNNATTCATAATAATTGCTNNNNNNNNNNANTTATATTAAAACTGATATCATACAAAAATAAATAATCATT
+
AB?B>A>@@babb...@babb>3=&561&&&&4...@bbbb>BB%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@HWW-T7400-009_1_120_1879_1950#ANNACG/1
ACTATTATGTCTATTCGAGAGCATNATCNNNNGTACTGGATTATCGATNNNNNNNNNNGNGACATCAATTTTAGTAAACAGTTTCCATTCCTTACTTTAT
+
ABB9BBABBABB<BABB?B8>A?8&9A;&&&&9=A>BB@@A%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
paired end file 2
@HWW-T7400-009_1_120_1879_1134#ANNACG/2
GCGCAGTATATAATTTCCATTGTTTTTCTCAAGTATAACATACTGTATGATATTGTCGGATTAAAATAAAAATAAAATATTAATGATTATTTATTTTTGT
+
?cccbbaba?bbcbcc=b?...@b87abbc?7>b...@5?@@4=BBBCBC9BBCB?AABB;=AA@;>B<5...@?
@?6-?BBBBB>=BC???B=-AB?9=9ABB
@HWW-T7400-009_1_120_1879_1950#ANNACG/2
TCCTCGACATGCACTGAATAATTCATTGTGGATAGGATCATATATTACACCTATTTCAATATTTTTGTTTTGTTGTACAGCAATAGATAAAGTAAGGAAT
+
BBA@@B4<3BB?7=?=<=B7ABBB>b...@a@A1;BBA?AA<<=@>==?9B?>B@>a...@=0=6*a@B:89)>?
5B?>B@>=...@bb>4B2A=-53>B*&8...@81@
_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
