This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How To Run Bfast In Cluster Computing

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

clustering next-gen sequencing

3 answers

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__

Thanks a lot for your kind help.

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.

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

together in a batch script seems the simplest.

Can you show me a piece of sample please??

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

Log in to answer this question.