SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics
Similar Threads
Thread Thread Starter Forum Replies Last Post
For MAQ: Is there a Tool to convert sanger-format fastq file to illumina-fotmat fastq byb121 Bioinformatics 6 12-20-2013 01:26 AM
MAQ-convert fasta to fastq rururara Bioinformatics 0 12-07-2011 11:06 PM
Convert FASTQ to QSEQ nguyendofx Bioinformatics 2 11-13-2011 08:03 PM
how to convert fastq to export or qseq format? feng Bioinformatics 3 06-15-2011 05:46 AM
From fastq to qseq? DarioC Bioinformatics 3 11-13-2010 02:09 PM

Reply
 
Thread Tools
Old 05-05-2009, 12:07 PM   #1
tgenahmet
Member
 
Location: Phoenix

Join Date: Apr 2009
Posts: 11
Default convert qseq to fastq that maq or bwa can use

Hi all,

I'm new here.

I have the Illumina pipeline 1.3 installed. I'd like to use BWA for aligning. I need to get from qseq files to phred fastq files.

Is there a tool that already exists to accomplish this?
tgenahmet is offline   Reply With Quote
Old 05-06-2009, 06:56 AM   #2
acnoll
Member
 
Location: Kansas City

Join Date: Mar 2008
Posts: 14
Default

Since bwa doesn't use quality values one can disregard the calibrations or scaling that Illumina implements as part of the pipeline. A simple awk command should get you filtered input from the qseq files for bwa.

cat s_1_1_*_qseq.txt | awk -F '\t' '{gsub(/\./,"N", $9); if ($11 > 0) print "@"$1$2$3$4$5$6$8"\n"$9"\n""+"$1$2$3$4$5$6$8"\n"$10}' > s_1_1_bwa.txt

Or you could use the sequence.txt files directly.
acnoll is offline   Reply With Quote
Old 05-08-2009, 02:09 AM   #3
chuck
Member
 
Location: china

Join Date: Apr 2009
Posts: 13
Default still need conversion of qseq into fastq

Hello,

I have been working with the previous Solexa format, with two separate files *seq.txt and *prb.txt, but just starting receiving a second round of data in the qseq format.

Converting them into the previous *seq.txt or a fasta file is not a problem but still I would like to run these files in maq to align them against my reference.

I am a bit confused as well about the issue of (@Phred+64) or (@solexa+64). After conversion, all of the values range from 2 to almost 40, so I assume these are the former.

But, I have never run maq either and it seems it needs (@Phred+33), correct? I can write something to convert them, if no one else has a script handy, but I want to make sure that I am doing it properly.

All the best,
Chuck
chuck is offline   Reply With Quote
Old 05-08-2009, 09:21 AM   #4
TylerBackman
Member
 
Location: Riverside, CA

Join Date: Oct 2008
Posts: 13
Default

This perl function I wrote will convert all of the qseq files for a given lane into proper Phred fastq files. Sorry the forum seems to be removing all of the whitespace... You can PM me for a text file if you'd like.

Use it like this:
#!/usr/bin/perl
use strict;
use Carp;
generate_sequence_list("bustard_path", "0", "8", "output.fastq"); # where "0" means single end and "8" means lane 8

**** begin code ****
sub generate_sequence_list {
# **** BEGIN CONFIG OPTIONS ****
my $bustard_path = $_[0];
my $pair = $_[1]; # 0=single end, 1=first pair, 2=second pair
my $lane = $_[2];
my $output_fastq_file = $_[3];
# **** END CONFIG OPTIONS ****

my $this_tile = 1;
my $qfilter = "";

open(OUTFASTAQFILE, "> $output_fastq_file");

if($pair > 0){
$pair = "_" . $pair . "_" ;
} else {
$pair = "_1_";
}

while(-r $bustard_path . "/s_" . $lane . $pair . sprintf("%04d", $this_tile) . "_qseq.txt"){
my $filename = $bustard_path . "/s_" . $lane . $pair . sprintf("%04d", $this_tile) . "_qseq.txt";
open(INFILE, "< $filename");
foreach my $thisline (<INFILE>) {
my @this_line = split("\t", $thisline);
croak("Error: invalid column number in $filename\n") unless(scalar(@this_line) == 11);
if($this_line[10] == 1) {
$qfilter = "Y";
} else {
$qfilter = "N";
}
# Convert quality scores
my $quality_string = $this_line[9];
my @quality_array = split(//, $quality_string);
my $phred_quality_string = "";
# convert each char to Phred quality score
foreach my $this_char (@quality_array){
my $phred_quality = ord($this_char) - 64; # convert illumina scaled phred char to phred quality score
my $phred_char = chr($phred_quality + 33); # convert phred quality score into phred char
$phred_quality_string = $phred_quality_string . $phred_char;
}

# replace "." gaps with N
$this_line[8] =~ s/\./N/g;

# output line
print OUTFASTAQFILE "@" . $this_line[2] . ":" . $this_line[3] . ":" . # output label line
$this_line[4] . ":" . $this_line[5] . ":" . $qfilter . "\n" .
$this_line[8] . "\n+\n" . # output sequence
$phred_quality_string . "\n"; # output quality string
}
close(INFILE);
$this_tile++;
}
$this_tile--;
croak("Error: 99 or less tiles in lane\n") unless($this_tile > 99);
print "\tFound $this_tile tiles in lane $lane.\n";

close(OUTFASTAQFILE);
}
**** end code ****
__________________
@1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
+
""""""""""""""""""""""""""""""""""""

Last edited by TylerBackman; 05-08-2009 at 09:34 AM.
TylerBackman is offline   Reply With Quote
Old 08-20-2010, 09:11 AM   #5
ashrafi_h
Junior Member
 
Location: Davis

Join Date: Jan 2010
Posts: 7
Default Usage

Quote:
Originally Posted by TylerBackman View Post
This perl function I wrote will convert all of the qseq files for a given lane into proper Phred fastq files. Sorry the forum seems to be removing all of the whitespace... You can PM me for a text file if you'd like.

Use it like this:
#!/usr/bin/perl
use strict;
use Carp;
generate_sequence_list("bustard_path", "0", "8", "output.fastq"); # where "0" means single end and "8" means lane 8

**** begin code ****
sub generate_sequence_list {
# **** BEGIN CONFIG OPTIONS ****
my $bustard_path = $_[0];
my $pair = $_[1]; # 0=single end, 1=first pair, 2=second pair
my $lane = $_[2];
my $output_fastq_file = $_[3];
# **** END CONFIG OPTIONS ****

my $this_tile = 1;
my $qfilter = "";

open(OUTFASTAQFILE, "> $output_fastq_file");

if($pair > 0){
$pair = "_" . $pair . "_" ;
} else {
$pair = "_1_";
}

while(-r $bustard_path . "/s_" . $lane . $pair . sprintf("%04d", $this_tile) . "_qseq.txt"){
my $filename = $bustard_path . "/s_" . $lane . $pair . sprintf("%04d", $this_tile) . "_qseq.txt";
open(INFILE, "< $filename");
foreach my $thisline (<INFILE>) {
my @this_line = split("\t", $thisline);
croak("Error: invalid column number in $filename\n") unless(scalar(@this_line) == 11);
if($this_line[10] == 1) {
$qfilter = "Y";
} else {
$qfilter = "N";
}
# Convert quality scores
my $quality_string = $this_line[9];
my @quality_array = split(//, $quality_string);
my $phred_quality_string = "";
# convert each char to Phred quality score
foreach my $this_char (@quality_array){
my $phred_quality = ord($this_char) - 64; # convert illumina scaled phred char to phred quality score
my $phred_char = chr($phred_quality + 33); # convert phred quality score into phred char
$phred_quality_string = $phred_quality_string . $phred_char;
}

# replace "." gaps with N
$this_line[8] =~ s/\./N/g;

# output line
print OUTFASTAQFILE "@" . $this_line[2] . ":" . $this_line[3] . ":" . # output label line
$this_line[4] . ":" . $this_line[5] . ":" . $qfilter . "\n" .
$this_line[8] . "\n+\n" . # output sequence
$phred_quality_string . "\n"; # output quality string
}
close(INFILE);
$this_tile++;
}
$this_tile--;
croak("Error: 99 or less tiles in lane\n") unless($this_tile > 99);
print "\tFound $this_tile tiles in lane $lane.\n";

close(OUTFASTAQFILE);
}
**** end code ****
Could you please add usage to the code or an example as per how to run the code. I put it the directory where all of my qseq files are. I changed the line and pair parameters according to my file names but the $bustard_path is not clear to me.

Thanks a lot.
ashrafi_h is offline   Reply With Quote
Old 08-20-2010, 11:35 AM   #6
TylerBackman
Member
 
Location: Riverside, CA

Join Date: Oct 2008
Posts: 13
Default

Quote:
Originally Posted by ashrafi_h View Post
Could you please add usage to the code or an example as per how to run the code. I put it the directory where all of my qseq files are. I changed the line and pair parameters according to my file names but the $bustard_path is not clear to me.

Thanks a lot.
"bustard_path" is the path on your filesystem where the base caller output is located. If you're running it from within the folder where your qseq files are, you can probably use "." to denote the current directory
__________________
@1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
+
""""""""""""""""""""""""""""""""""""
TylerBackman is offline   Reply With Quote
Old 08-20-2010, 12:15 PM   #7
ashrafi_h
Junior Member
 
Location: Davis

Join Date: Jan 2010
Posts: 7
Default

I tried to run the script within the folder that I have the qseq files. something like the following. Name of qseq files are also like
s_6_1_0074_qseq.txt = lane 6 read 1

>perl Perl_qseq_to_fastq.pl .

Results:
perl Perl_qseq_to_fastq.pl .
Error: 99 or less tiles in lane
at Perl_qseq_to_fastq.pl line 61
main::generate_sequence_list('bustard_path', 1, 6, 'output.fastq') called at Perl_qseq_to_fastq.pl line 4

or

>Perl_qseq_to_fastq.pl ..

or from outside of the folder, lets say from one folders outside of qseq file folder

perl Perl_qseq_to_fastq.pl /BaseCalls
Error: 99 or less tiles in lane
at Perl_qseq_to_fastq.pl line 61
main::generate_sequence_list('bustard_path', 1, 6, 'output.fastq') called at Perl_qseq_to_fastq.pl line 4
ashrafi_h is offline   Reply With Quote
Old 01-05-2011, 04:42 AM   #8
rgregor
Member
 
Location: Ljubljana

Join Date: Jun 2010
Posts: 11
Default

Can somebody explain what is the qseq quality control filter (column 11, either 0 or 1) and how is it obtained (Illumina)? Do you usually disregard sequences with qseq_filter = 0?

http://jumpgate.caltech.edu/wiki/QSeq, see column 11

tnx,
Gregor
rgregor is offline   Reply With Quote
Old 02-19-2011, 01:13 PM   #9
lingfung.tang
Member
 
Location: San Francisco

Join Date: Dec 2010
Posts: 10
Default

and for the python script, what's the trim and pass_qc_msg refer to?
lingfung.tang is offline   Reply With Quote
Reply

Tags
bwa, fastq, maq, qseq

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off



All times are GMT -8. The time now is 08:07 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2022, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO