SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to merge multiple sequencing runs vinay052003 Bioinformatics 4 01-31-2012 04:34 AM
Using multiple MIDs in titanium sequence runs JurgenP 454 Pyrosequencing 8 07-23-2010 08:24 AM
Maq multiple hits file Lesley Bioinformatics 1 12-09-2009 04:12 AM
maq multiple hits information? baohua100 Bioinformatics 15 08-13-2009 02:31 PM
Maq multiple map output (-H) torrey Bioinformatics 0 03-06-2009 11:34 AM

Reply
 
Thread Tools
Old 05-05-2009, 04:52 AM   #1
Layla
Member
 
Location: London

Join Date: Sep 2008
Posts: 58
Default multiple runs and maq

How do you do an alignment (./ maq map command) in maq when you have run more than one lane on a Solexa machine on the same sample?
Can you simply append the reads,

./maq map n.read1.bfq in.read2.bfq in.read3.bfq etc etc?


Cheers
L
Layla is offline   Reply With Quote
Old 05-05-2009, 07:12 AM   #2
zee
NGS specialist
 
Location: Malaysia

Join Date: Apr 2008
Posts: 249
Default

Layla,

If you have tried that command by now you know it would never work.

At most 2 bfq files can be given and these are assumed to contain paired-end reads. If one file is given then just single-lane.

So, if you have all your *sequence.txt files use a for loop on these astq files:
Quote:

for file in `ls *sequence.txt`
do
maq fastq2bfq $file $file.bfq
maq map $file.map genome.bfa $file.bfq
done
zee is offline   Reply With Quote
Old 05-05-2009, 07:24 AM   #3
dakl
Member
 
Location: sweden

Join Date: May 2009
Posts: 15
Default

Zee,

Would it be possible to parallelize this over several CPU cores in a simple manner? Kind of like PBS for cluster jobs, but locally.

cheers
D
dakl is offline   Reply With Quote
Old 05-05-2009, 07:37 AM   #4
zee
NGS specialist
 
Location: Malaysia

Join Date: Apr 2008
Posts: 249
Default

Dakl,

I use novoalign for most of my mutli-core jobs, but it should be possible to do something similar with maq.

If you have a large database to search you might run into some problems. I would split all my files into batches of N-1, N= no. CPUs - I like to keep one free for system, IO,etc.

Quote:

ls *.fastq | split -l <N-1> BATCH


Then for each batch run a loop

Quote:
for file in `cat BATCH...`
do
maq fastq2bfq $file $file.bfq
maq map $file.map genome.bfa $file.bfq
done &
And dont forget the "&" which places each loop in the background.
zee is offline   Reply With Quote
Old 05-13-2009, 08:35 AM   #5
Layla
Member
 
Location: London

Join Date: Sep 2008
Posts: 58
Default

Hi zee,

Thankyou for your reply but I have been receiving a bizarre error:Assertion failed: (fp_bfa), function ma_match, file match.cc, line 516

for file in `ls *.bfq`
do
./maq map $file.map genome.bfa $file.bfq
done

From the paired end experiment I have a total of 3 pairs stored in one folder:
s_1.bfq s_2.bfq
t_1.bfq t_2.bfq
u_1.bfq u_2.bfq

I didnt understand how/where the ./maq map loop looks at s_1.bfq s_2.bfq and then t_1.bfq t_2.bfq, u_1.bfq u_2.bfq.

Thank you for your help

L
Layla is offline   Reply With Quote
Old 05-13-2009, 08:42 AM   #6
zee
NGS specialist
 
Location: Malaysia

Join Date: Apr 2008
Posts: 249
Default

OK, it is a simple change to the ff:

Code:
for base in `echo s t u`; do
  ./maq map $base.map genome.bfa $base"_1.bfq" $base"_2.bfq"
done

Quote:
Originally Posted by Layla View Post
Hi zee,

Thankyou for your reply but I have been receiving a bizarre error:Assertion failed: (fp_bfa), function ma_match, file match.cc, line 516

for file in `ls *.bfq`
do
./maq map $file.map genome.bfa $file.bfq
done

From the paired end experiment I have a total of 3 pairs stored in one folder:
s_1.bfq s_2.bfq
t_1.bfq t_2.bfq
u_1.bfq u_2.bfq

I didnt understand how/where the ./maq map loop looks at s_1.bfq s_2.bfq and then t_1.bfq t_2.bfq, u_1.bfq u_2.bfq.

Thank you for your help

L
zee is offline   Reply With Quote
Old 05-13-2009, 09:23 AM   #7
Layla
Member
 
Location: London

Join Date: Sep 2008
Posts: 58
Default

Thanx Zee. Since multiple runs are carried out to increase the number of reads, why is it that a separate .map file is being created for each pair (so a total of 3)? Is the purpose not to merge all the pairs and generate a single .map file to increase genome coverage?

Actually whilst writing, is this where I can use the ./maq merge command?

Cheers
L
Layla is offline   Reply With Quote
Old 05-13-2009, 11:34 AM   #8
jnfass
Member
 
Location: Davis, CA

Join Date: Aug 2008
Posts: 88
Default

You got it ... after you've generated all the maps, use maq merge to combine them into one map, from which you can generate a pileup, consensus, etc ...
jnfass is offline   Reply With Quote
Old 05-14-2009, 01:18 AM   #9
dakl
Member
 
Location: sweden

Join Date: May 2009
Posts: 15
Default

Hi all,

Since Maq is optimized for ~2M reads as input, I managed to do the following:

Code:
time maq fastq2bfq -n 2000000 ../50a_fastq.single.fastq 50a
to create several bfq-files containing the reads, and then use the perl module Parallel::ForkManager to fork the process. See the script below for details.

Code:
#!/usr/bin/perl -w

use strict;
use Parallel::ForkManager;

my $pm = new Parallel::ForkManager(4); # number of parallel processes is 4
while(<>){
        chomp;

        # Forks and returns the pid for the child:
        my $pid = $pm->start and next; 
        
        qx/ maq match -c $_.map ~\/hg18\/hg18.bfa $_/; 
        
    $pm->finish; # Terminates the child process
}
dakl is offline   Reply With Quote
Reply

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 05:21 PM.


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