[Bioperl-l] fastq splitter
Sean O'Keeffe
limericksean at gmail.com
Wed Feb 29 16:30:05 UTC 2012
Hi Chris,
Here's the perldoc for fastq - it does seem to indicate that the optional
descriptor (+) must match the first header. (See DESCRIPTION).
Maybe this version doesn't include the updated code you mention which Peter
Cock has worked on.
Sean.
======================================
Bio::SeqIO::fastq(3) User Contributed Perl Documentation
Bio::SeqIO::fastq(3)
NAME
Bio::SeqIO::fastq - fastq sequence input/output stream
SYNOPSIS
################## pertains to FASTQ parsing only
##################
# grabs the FASTQ parser, specifies the Illumina variant
my $in = Bio::SeqIO->new(-format => ’fastq-illumina’,
-file => ’mydata.fq’);
# simple ’fastq’ format defaults to ’sanger’ variant
my $out = Bio::SeqIO->new(-format => ’fastq’,
-file => ’>mydata.fq’);
# $seq is a Bio::Seq::Quality object
while (my $seq = $in->next_seq) {
$out->write_seq($seq); # convert Illumina 1.3 to Sanger format
}
# for 4x faster parsing, one can do something like this for raw
data
use Bio::Seq::Quality;
# $data is a hash reference containing all arguments to be passed
to
# the Bio::Seq::Quality constructor
while (my $data = $in->next_dataset) {
# process $data, such as trim, etc
my $seq = Bio::Seq::Quality->new(%$data);
# for now, write_seq only accepts Bio::Seq::Quality, but may
be modified
# to allow raw hash references for speed
$out->write_seq($data);
}
DESCRIPTION
This object can transform Bio::Seq and Bio::Seq::Quality objects to
and from FASTQ flat file databases.
FASTQ is a file format used frequently at the Sanger Centre and in
next-gen sequencing to bundle a FASTA sequence and its quality data. A
typical FASTQ entry takes the from:
@HCDPQ1D0501
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT.....
+HCDPQ1D0501
!’’*((((***+))%%%++)(%%%%).1***-+*’’))**55CCF>>>>>>CCCCCCC65.....
where:
@ = descriptor, followed by one or more sequence lines
+ = optional descriptor (if present, must match first one),
followed by one or
more qual lines
FASTQ and Bio::Seq::Quality mapping
FASTQ files have sequence and quality data on single line or
multiple lines, and the quality values are single-byte encoded. Data are
mapped very simply to Bio::Seq::Quality instances:
Data Bio::Seq::Quality
method
------------------------------------------------------------------------
first non-whitespace chars in descriptor id^
descriptor line desc^
sequence lines seq
quality qual*
FASTQ variant namespace
^ first nonwhitespace chars are id(), everything else after (to
end of line)
is in desc()
* Converted to PHRED quality scores where applicable (’solexa’)
FASTQ variants
This parser supports all variants of FASTQ, including Illumina v 1.0
and 1.3:
variant note
-----------------------------------------------------------
sanger original
solexa Solexa, Inc. (2004), aka Illumina 1.0
illumina Illumina 1.3
The variant can be specified by passing by either passing the
additional -variant parameter to the constructor:
my $in = Bio::SeqIO->new(-format => ’fastq’,
-variant => ’solexa’,
-file => ’mysol.fq’);
or by passing the format and variant together (Bio::SeqIO will now
handle this and convert it accordingly to the proper argument):
my $in = Bio::SeqIO->new(-format => ’fastq-solexa’,
-file => ’mysol.fq’);
Variants can be converted back and forth from one another; however,
due to the difference in scaling for solexa quality reads, converting from
’illumina’ or ’sanger’ FASTQ to solexa is not recommended.
....
============================================
On 28 February 2012 21:40, Fields, Christopher J <cjfields at illinois.edu>wrote:
> That should work. Can you send the output of 'perldoc Bio::SeqIO::fastq'?
> That should indicate what is being called.
>
> chris
>
> On Feb 28, 2012, at 5:50 PM, Sean O'Keeffe wrote:
>
> > $ perl -MBio::Root::Version -e 'print $Bio::Root::Version::VERSION,"\n"'
> > 1.006001
> >
> > Isn't that 1.6.1 - does it need upgrading ?
> >
> > On 28 February 2012 18:36, Sean O'Keeffe <limericksean at gmail.com> wrote:
> > Could be. I'll check.
> >
> > On 28 February 2012 17:17, Fields, Christopher J <cjfields at illinois.edu>
> wrote:
> > That's a bit odd. Are you using an old version of the FASTQ parser? It
> was revised a while ago, prior to the v1.6.1 release (the error matches one
> in the older parser)
> >
> > chris
> >
> > On Feb 28, 2012, at 4:01 PM, Sean O'Keeffe wrote:
> >
> > > Hi Chris,
> > > Unfortunately the read pairs are not consecutive. It seems they are
> cat'd together.
> > > I could use split -l on the line number that they're glued together I
> guess.
> > > If this is an overnight job for a bunch of files, I can wait so don't
> mind using the module if it worked.
> > >
> > > Someone pointed out I need to switch $seqin->desc to $inseq->desc.
> > > However, now it spits out fasta output instead of fastq and returns a
> bunch of warnings: Seq/Qual descriptions don't match; using sequence
> description
> > >
> > > Hmm.
> > >
> > > On 28 February 2012 16:50, Fields, Christopher J <
> cjfields at illinois.edu> wrote:
> > > Sean,
> > >
> > > If you trust the data enough, in that:
> > >
> > > 1) each record is 4 lines,
> > > 2) mate pairs are consecutive in the file, and
> > > 3) that read 1 always preceeds read 2 in the pair,
> > >
> > > then I would simply iterate through 4 lines at a time and dump to the
> two separate files, maybe using a flip-flop or simple record count and
> modulus switch. You can always run a check on the header with a regex if
> you don't trust it completely.
> > >
> > > Just from the sanity point-of-view, unless you're doing a lot of
> validation I wouldn't use Bio::SeqIO::fastq, unless you have some time on
> your hands and a relatively low number of seqs (it's notoriously slow at
> the moment).
> > >
> > > chris
> > >
> > > On Feb 28, 2012, at 3:11 PM, Sean O'Keeffe wrote:
> > >
> > > > Hi,
> > > > I'm trying to write a quick script to separate one large PE fastq
> file into
> > > > 2 separate files, one for each mate pair
> > > >
> > > > The file is of the format (mate1)
> > > > @HWI-ST156:445:C0EDLACXX:4:1101:1496:1039 1:N:0:ATCACG
> > > > CTGCTGGTAGTGCCCAAAGACCTCGAATACAATGGGCTTGGTTTTGATGT
> > > > +
> > > > BCCFFFFEHHHHHJJJJJHIIJIJJIIGIJJJJJJJIJJJI?FHJJIIJA
> > > >
> > > > && (mate2)
> > > >
> > > > @HWI-ST156:445:C0EDLACXX:4:2308:20877:199811 2:Y:0:ATCACG
> > > > TCATAAAAATAACAAAACCACCACCCCATACAAACTCTACTCATCTCCAC
> > > > +
> > > > ##################################################
> > > >
> > > >
> > > > My idea is to separate using a regex such that / 1:/ would be the
> first
> > > > mate pair and / 2:/ would go in the second mate file.
> > > > I implemented the code below but each output file is empty. Can
> someone
> > > > spot my error?
> > > >
> > > > Thanks,
> > > > Sean.
> > > >
> > > > my $infile = shift;
> > > > my $outfile1 = $infile."_1";
> > > > my $outfile2 = $infile."_2";
> > > >
> > > > my $seqin = Bio::SeqIO->new(
> > > > -file => "<$infile",
> > > > -format => "fastq",
> > > > );
> > > > my $seqout1 = Bio::SeqIO->new(
> > > > -file => ">$outfile1",
> > > > -format => "fastq",
> > > > );
> > > >
> > > > my $seqout2 = Bio::SeqIO->new(
> > > > -file => ">$outfile2",
> > > > -format => "fastq",
> > > > );
> > > > while (my $inseq = $seqin->next_seq) {
> > > > if ($seqin->desc =~ / 1:/){
> > > > $seqout1->write_seq($inseq);
> > > > } else {
> > > > $seqout2->write_seq($inseq);
> > > > }
> > > > }
> > > > _______________________________________________
> > > > Bioperl-l mailing list
> > > > Bioperl-l at lists.open-bio.org
> > > > http://lists.open-bio.org/mailman/listinfo/bioperl-l
> > >
> > >
> >
> >
> >
>
>
More information about the Bioperl-l
mailing list