websoft/bin/run_large_data.pl

   1 #!/usr/bin/perl
   2
   3 =head1 NAME
   4
   5 run_large_data.pl - script to run a large dataset against Jpred using the cluster
   6
   7 =cut
   8
   9 use strict;
  10 use warnings;
  11
  12 use Getopt::Long;
  13 use Pod::Usage;
  14 use File::Basename;
  15
  16 # temporary. Should put this module somewhere more global.
  17 use lib '/homes/ccole/lib';
  18 use Cluster::ArrayJob 0.3;
  19
  20 my $file;
  21 my $scriptName = "jpredLarge.pl";
  22 my $title = 'jpredLarge';
  23 my $VERBOSE = 0;
  24 my $DEBUG = 0;
  25 my $help;
  26 my $man;
  27
  28 GetOptions (
  29    'in=s'      => \$file,
  30    'title=s'   => \$title,
  31    'verbose!'  => \$VERBOSE,
  32    'debug!'    => \$DEBUG,
  33    'man'       => \$man,
  34    'help|?'    => \$help,
  35 ) or pod2usage();
  36
  37 pod2usage(-verbose => 2) if ($man);
  38 pod2usage(-verbose => 1) if ($help);
  39 pod2usage(-msg => 'Please supply a valid filename.') if (!$file or !-e $file);
  40
  41 print "Reading fasta file...\n";
  42 my $total = split_fasta($file);
  43 print "Created $total fasta files\n";
  44
  45 write_script($scriptName);
  46
  47
  48 ## set up Array Job and submit perl script
  49 ## Just use all the defaults.
  50 print "Submitting $total Jpred jobs to the cluster as an array job...\n" if $VERBOSE;
  51 my $sgeArray = Cluster::ArrayJob->new();
  52 $sgeArray->taskRange("1-$total");
  53 $sgeArray->queue('64bit-pri.q');
  54 $sgeArray->setCWD();
  55 $sgeArray->jobname($title);
  56 $sgeArray->setJobShare(0);
  57 $sgeArray->setPriority(-100);
  58 $sgeArray->setResourceRequest( 'ram' => '6G' );
  59 $sgeArray->setEnv( 'PERL5LIB' => '/homes/www-jpred/live/lib' );
  60 $sgeArray->submit($scriptName);
  61
  62 print "Checking status of array job...\n" if $VERBOSE;
  63 while (1) {
  64    my $status = $sgeArray->getStatus() or die "ERROR - unable to get SGE job status: ", $sgeArray->error();
  65    if ($status eq '-1') {
  66       print "Job has finished\n" if $VERBOSE;
  67       last;
  68    } elsif ($status eq '1') {
  69       die "ERROR - unable to get SGE job status: ", $sgeArray->error();
  70    } else {
  71       my $out = '';
  72       foreach my $k (sort keys %$status) {
  73          $out .= " $k:".$status->{$k};
  74       }
  75       print "Job status: $out\n" if $VERBOSE;
  76    }
  77    sleep(30);
  78 }
  79
  80 ## Finally, check that all jobs completed successfully.
  81 my $jobID = $sgeArray->jobid();
  82 my @files = glob "$title.e$jobID.*";
  83
  84 my $errors = 0;
  85 ## check that the correct number of SGE files exist - should be the same as the number of tasks
  86 my $nFiles = scalar @files;
  87 if ($nFiles != $total) {
  88    warn "Warning - found $nFiles SGE output files where $total expected\n";
  89    ++$errors;
  90 }
  91
  92 ## check that the right number of Jnet outputs were generated
  93 my @jnets = glob "*.jnet";
  94 my $nJnets = scalar @jnets;
  95 if ($total != $nJnets) {
  96    warn "Warning - found $nJnets Jnet predictions where $total expected\n";
  97 }
  98 print "All Jpred searches completed!\n";
  99 exit;
 100
 101
 102 ## write script with code from __DATA__ below
 103 sub write_script {
 104    my $file = shift;
 105
 106    open(my $OUT, ">$file") or die "ERROR - unable to open file '$file': $!\n";
 107    while (<DATA>) {
 108       print $OUT $_;
 109    }
 110    close($OUT);
 111    close(DATA);
 112 }
 113
 114 ## short function to split a Fasta file
 115 ## into one file per sequence
 116 sub split_fasta {
 117    my $file = shift;
 118
 119    open(my $FAS, "<", $file) or die "ERROR - unable to open '$file': ${!}\nDied";
 120    my $num = 0;
 121    my $OUT;
 122    while(<$FAS>) {
 123       ## check first line has appriate start
 124       if ($. == 1) {
 125          die "ERROR - this file does not appear to be a Fasta file\n" unless (/^>/);
 126       }
 127       ## for each new record open a new file and close the preceeding one
 128       if (/^>/) {
 129          ++$num;
 130          my $out = "$num.fasta";
 131          close($OUT) if ($OUT); # needed for first one; can't close if nothing's open
 132          open($OUT, ">", $out) or die "ERROR - unable to open '$out' for write: ${!}\nDied";
 133       }
 134       print $OUT $_;
 135
 136    }
 137    close($OUT);
 138    return($num)
 139 }
 140
 141 =head1 SYNOPSIS
 142
 143 run_large_data.pl --in <file> [--verbose] [--debug] [--man] [--help]
 144
 145 =head1 DESCRIPTION
 146
 147 Running a large set of sequences against Jpred can be a but of a pain as Jpred can be very memory intensive and can take a while.
 148
 149 Therefore, use of the cluster is best, but managing the jobs can be a bit fiddly. So, this script does it all for you!
 150
 151 Just provide a Fasta file with all the sequences you want to run and the script will submit each of them to the cluster separately, monitor their progress and report if there have been any failures. Job done!
 152
 153 =head1 OPTIONS
 154
 155 =over 5
 156
 157 =item B<--in>
 158
 159 Input file (fasta format).
 160
 161 =item B<--verbose|--no-verbose>
 162
 163 Toggle verbosity. [default:none]
 164
 165 =item B<--debug|--no-debug>
 166
 167 Toggle debugging output. [default:none]
 168
 169 =item B<--help>
 170
 171 Brief help.
 172
 173 =item B<--man>
 174
 175 Full manpage of program.
 176
 177 =back
 178
 179 =head1 AUTHOR
 180
 181 Chris Cole <christian@cole.name>
 182
 183 =cut
 184
 185 __DATA__
 186 #!/usr/bin/perl
 187
 188 use strict;
 189 use warnings;
 190
 191 my $pwd = $ENV{PWD};  # get current directory
 192 my $dir = $ENV{TMPDIR};  # get SGE tmpdir
 193 my $task = $ENV{'SGE_TASK_ID'};  # get SGE array task ID
 194 die "ERROR - not in SGE array job context. Please submit as an array job.\n" if (!$task or $task eq 'undefined');
 195
 196 die "ERROR - file '$task.fasta' not found at '$pwd'. Check path.\n" unless (-e "$task.fasta");
 197 if (-s "$task.jnet") {
 198    print "A Jnet prediction already exists for '$task.fasta'. Skipping...\n";
 199    exit();
 200 }
 201
 202 my $cmd = "cd $dir && /homes/www-jpred/live/jpred/jpred --sequence $pwd/$task.fasta --output $task";
 203 print "Running CMD: $cmd\n";
 204 system($cmd) == 0 or die "ERROR - system() died\n";
 205 if (-s "$dir/$task.jnet") {
 206    print "Jnet successful! Copying files from cluster to $pwd...";
 207    system("cp $dir/$task.* $pwd/") ==0 or die "ERROR - system() died\n";
 208    print "Done\n";
 209 } else {
 210    #my $out = `ls -ltr $dir`;
 211    #print "ls -ltr:\n$out\n";
 212    print "No Jnet prediction found for '$task'. Something failed...\n"
 213 }