Filtering raw sequencing reads

Question

I have to filter the raw sequencing reads based on the following criteria:

1. Remove reads containing adapters

2. Remove reads containing N > 10% (N represents base that could not be determined)

3. Remove reads where the Qscore (Quality value) of over 50% bases of the read is <= 5

The reads are paired-end fastq files. I have been able to do the first part by using trimmomatic:

trimmomatic \
  PE in_1.fq.gz in_2.fq.gz \
  pair_out_1 unpair_out_1 pair_out_2 unpair_out_2 \
  ILLUMINACLIP:adapter.fa:2:30:5

The other 2 steps are a bit more challenging since I have not found a tool that perform exactly what is required. I could write some custom code of course but I would like to use the existing tools if possible. No need to reinvent the wheel.

Do you have any suggestions on how to proceed?

Answer 1

I think you can come close to this with cutadapt:

cutadapt \
  -a ${FWD_ADAPTER} -A ${REV_ADAPTER} \         ## adapter
  --max-n 0.1 \                                 ## do not allow > 10% Ns
  -q 5 \                                        ## remove bases with basequal < 5
  --minimum-length (readlength/2) \             ## minimum length after trimming
  -o read1_trimmed.fq -p read2_trimmed.fq \     ## outputs
  read1.fq read2.fq                             ## inputs

This will remove reads with more than 10% of bases being N and remove bases with base quality below 5. If you then set --minimum-length to half the read length, this should also remove reads where more than half of the bases had a base quality below 5 and therefore got trimmed away. Poor base quality in Illumina reads typically happens towards the 3' end so I think this little workaround is a sufficient proxy to achieve your goal. A base quality of 5 is very low in any case, so if sequencing went normal the typical quality should be notably higher. I would be surprised if this third filter is really necessary and will eliminate a notable number of reads.