[kmcarr]

My question is related, but not quite the same, which is: I need to get a sub-string of the sequence and the corresponding quality score for each of the entries, the file format untouched. My Illumina reads consists of 101bp long. I want remove the first 10bp and last 25bp then only the middle 66bp left.

The reason I need do this is my DNA sequence is methylated and the first 10 and last 25bp seem not having good quality in general so that I want get rid of them. My challenge is to remove both quality score and sequence correspondingly.

use Bio::SeqIO; use Bio::Seq::Quality;

$seqio = Bio::SeqIO->new('-format'=>'fastq' , '-file'=>'some.fasq');

my $out_fastq = Bio::SeqIO->new( -format => 'fastq', '-file'=> 'subset.fastq');

while((my $line = $seqio->next_seq() ) { 
# keep the id 
# substring the sequence (from 11 to 66) 
# substring the quality score (from 11 to 66) $out_fastq->write_seq($line); }

Or any tools there to do my job?

Thanks a lot!

Yifang

Yifang,

First rule of coding, never write your own when a tool already exist to do what you want. The FASTX-Toolkit has a utility to do exactly what you want, namely the fastx_trimmer. You give an input file (fasta or fastq), the postion of the first and last base you wish to keep (in your case 11 and 76) and it will produce a trimmed (bases and quality) file.

[yifangt]

Thanks a lot kmcarr! I was searching for it, just missed the inside. I will give it a try.

Best!
Yifang

[gpcr]

awk '{print ; getline } {print substr($0, 11, 76) ; getline; print ; getline ; print substr($0, 11, 76) }' input.fastq

[yifangt]

This is a cool script. Thanks gpcr! Yifang