multiple runs and maq

[Layla]

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

[zee]

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:

for file in `ls *sequence.txt`
do
maq fastq2bfq $file $file.bfq
maq map $file.map genome.bfa $file.bfq
done

[dakl]

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

[zee]

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.

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

Then for each batch run a loop

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.

[Layla]

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]

OK, it is a simple change to the ff:

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

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]

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

[jnfass]

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 ...

[dakl]

Hi all,

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

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.

#!/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
}