How To Run Bfast In Cluster Computing

Question

[Venky]

Hi,

I am trying to run Bfast tool in a PBS/Cluter computing but it has to be written with a script to run the entire protocol of bfast. Please help me in submitting the job using PBS scripts and bfast.

Thank you

Answer 1

[Michael]

It's a bit tricky to script bfast but very rewarding because its workflow is quite long. One idea would be to build a wrapper script around the indexing and alignment steps and then call these in a PBS script. I tried perl wrappers some time ago you can maybe copy the command sequences into your script, r just use the ideas or use the perl wrapper directly, but I don't know much about PBS scripts.

So here's the build index script without any warranty, I remember I translated that from a python script, check with the bfast documentation for the choice of the masks, creates a bfast path:

#!/usr/bin/env perl

use strict;
use warnings;
use File::Path;
use File::Basename;

my $dir_name = "bfast";
my $window_size = 14;
my @bfast_nt_masks = (
   "1111111111111111111111",
   "1111101110111010100101011011111",
   "1011110101101001011000011010001111111",
   "10111001101001100100111101010001011111",
   "11111011011101111011111111",
   "111111100101001000101111101110111",
   "11110101110010100010101101010111111",
   "111101101011011001100000101101001011101",
   "1111011010001000110101100101100110100111",
   "1111010010110110101110010110111011",
    );

my $in = $ARGV[0];
print($in,"\n");

my ($name,$path,$suffix) = fileparse($in);

my $outpath = File::Spec->catfile($path, $dir_name);
mkpath($outpath);
chdir($outpath) or die "Could not change to $outpath $!";
system( sprintf("bfast fasta2brg -f %s -A 0", $in) );
my $i = 1;
foreach my $mask (@bfast_nt_masks) {
  my $cmd = sprintf ("bfast index -d 1 -n 4 -f %s -A 0 -m %s -w %s -i %s",
             ($in, $mask, $window_size, $i++));
  print "running $cmd :\n";
  system($cmd);
}

And the alignment script for multiple files, creates a path bfastout:

#!/usr/bin/env perl

use strict; 
use warnings;

use File::Path;
use File::Basename;

my $CMD = `which bfast`;
chomp $CMD;
my $PARS_MATCH = "-A0  -n4 -i1";
my $PARS_POST =  "-a0 -O0"; 


my $genome = shift @ARGV;


mkpath("./bfastout");

foreach my $infile (@ARGV) {
  runbfast($infile);
}
print STDERR "done\n";

sub runbfast{
  my $infile = shift;
  my $outfile = "./bfastout/".(basename($infile)).".bfast";
  my $cmdline1 = join(" ", $CMD, "match", $PARS_MATCH, '-f', $genome, "-r", $infile, ">" , "$outfile-tmp1" ); 
  my $cmdline2 = join ( " ", $CMD, "localalign", "-m", "$outfile-tmp1",  '-f', $genome, ">" , "$outfile-tmp2" );   
  my $cmdline3 = join ( " ", $CMD, "postprocess","-i", "$outfile-tmp2", '-f', $genome, $PARS_POST, ">", $outfile);
  print ($cmdline1,"\n");
  system ($cmdline1);
  print ($cmdline2,"\n");
  system ($cmdline2);
  print ($cmdline3,"\n");
  system ($cmdline3);


}

__END__

[Venky]

Thanks a lot for your kind help.

Answer 2

[brentp]

We have done this by splitting the reads (as recommended by bfast and probably any alignment software). From there, for each set of reads you use a scripting language to create a shell script that does:

#PBS params
#PBS otherparams
bfast match [params] $input.fastq > $output.bmf
bfast localalign [params] $output.bmf > $output.baf
bfast postprocess [params] $output.baf > $output.sam
samtools stuff [sort, index bam] > $output.sorted.bam

Then you just submit each of those jobs to qsub. you can even use a template and just fill in the filenames.

[Venky]

The jobs for each step has to be given separately in qsub or together in a batch script??

[brentp]

together in a batch script seems the simplest.

[Venky]

Can you show me a piece of sample please??

Answer 3

[Sophia]

The bfast package itself contains actually a script named bfast.submit.pl that is supposed to serve this purpose...