Using the Bioinfo Cluster

A script file is simply a file, typically with a .sh extension, which is used to tell Torque how to run a job; this can be as simple as a couple lines. The file starts with the line "#!/bin/sh" telling Linux to use /bin/sh to run the following list of Linux commands which will setup and run the application. For a simple single processor job, the following job will run blast on a single node, using ONE processor. If you would have ran:

laurichj@bioinfo:~/$ blastall -p blastp -d /data/NCBI/blast/nr/nr -i ~/blast/proteins.fa -e 1e-19

on the command line, you can transform that into a Torque script as:

#!/bin/sh
blastall -p blastp -d /data/NCBI/blast/nr/nr -i ~/blast/proteins.fa -e 1e-19

This script might be written into a file named blast.sh. However, a more flexible version would be:

#!/bin/sh 

DIR=~/blast 
INPUT=proteins.fa
PROGRAM=blastp
DATABASE=/data/NCBI/blast/nr/nr
OPTIONS="-e 1e-10" 

cd $DIR; blastall -p $PROGRAM -d $DATABASE -i $INPUT $OPTIONS

	This tells the computer to use /bin/sh to run this script.
	By setting variables to point to each of the blast locations, we can simplify the next blast jobs, simply by changing the names. The DIR value is the working directory for the job to be ran, otherwise the job will be run you home directory. The INPUT, DATABASE and PROGRAM are passed as is to blastall. However, OPTIONS is in quotes to prevent the shell from parsing it.
	Finally, we execute blastall.

This script will run a blast on a free node. All output from the program is saved into files. You no longer have to find a free node, or wait for nodes to be free or even log into them. All you need to do is submit that file to torque, using qsub as follows:

laurichj@bioinfo:~/$ qsub blast.sh
123.bioinfo.ucr.edu

The job is now queued, and will be ran on the next free node. The second line is the output from qsub. The number, 123, is the Job ID Number which is used to track the progress of the job with qstat. qstat will list the current running and queued jobs.

laurichj@bioinfo:~/others/blast$ qstat
Job id           Name             User             Time Use S Queue
---------------- ---------------- ---------------- -------- - -----
123.bioinfo      blast.sh         laurichj                0 R batch

Once the job is no longer listed, it has finished. The files blast.sh.o123 and blast.sh.e123 will now be available. The blast.sh will be whatever the name of your job script is, and the number at the end will be the job id number. The o file is the normal output, and in the case of blast will be the file we are interested in. The e file is the error output. Both of these will always exist, but the error output will usually be empty.

If the BLAST database is larger than 2GB, then the file needs to be moved from bioinfo to the local node. For a single processor job, using cp over NFS is sufficient. By copying the file prior to running blast and removing it after, the job will run faster.

#!/bin/sh

DIR=~/blast
INPUT=proteins.fa
PROGRAM=blastp
SOURCE_DIR=/data/local/NCBI/blast/nr/
LOCAL_DIR=/scratch/$USER/
DATABASE=nr
OPTIONS="-e 1e-10"

mkdirheir $LOCAL_DIR
cd $DIR
/bin/cp -r $SOURCE_DIR/$DATABASE.* $LOCAL_DIR/ 
blastall -p $PROGRAM -d $LOCAL_DIR/$DATABASE -i $INPUT -a 2 $OPTIONS
rm -rf $LOCAL_DIR 

	This will copy all database files to the local node.
	This will remove the file afterwords.

Moving the files in will only increase performance if the file(s) is/are too large to fit into memory and the program scans through the file(s) more than once. If we were to be blasting a single sequence, blast would only scan the database once therefore, copying the database would be wasted effort.

Sometimes a program needs user input to run, or you are just running a single blast and don't want to bother with writing a script file. In this case, in place of a script file you can tell Torque to run an interactive job by giving it the -I (capital i) switch. An interactive job behaves in much the same way as a normal batched job in that it will mark the node(s) as used, wait for enough nodes to become free and so forth. The difference is that qsub will no longer immediatly return. In stead it will wait for your job to execute and then provide you a shell in which you can run commands. For example, to run a simple blast you could:

laurichj@bioinfo:~/$ qsub -I
qsub: waiting for job 124059.bioinfo.ucr.edu to start
qsub: job 124059.bioinfo.ucr.edu ready

laurichj@node01:~$ cd blast/
laurichj@node01:~$ blastall -p blastp -d /data/TIGR/arabidopsis/blast/ATH1.pep  -i seq.fa

The worst, but simplest, way to reuse shell scripts has already be hinted at: use multiple directories, one for each job, each with a slightly modified script. Using the previous BLAST script, it is easy to illustrate this idea. Imagine you wished to run a set of proteins against multiple databases. In this case we will use Arabidopsis thaliana and Oryza sativa. To do this you can simply make two directories: ~/blast/ath and ~/blast/osa. Then place a copy of the blast script in each directory, named blast.sh with the INPUT and DATABASE variables modified. However, for more than a few combinations of input and databases, this becomes increasingly tedious.

Rather, a simple evolution of this will allow you to use the same file for multiple jobs. By leaving the options unspecified you can make a script act as a generic schema for a job. For instance, blast is usually ran relatively similar each time. The following file would make a generic schema for blast:

#!/bin/sh
#PBS -v DB,IN

PROGRAM=blastp
OPTIONS="-e 1e-10"

cd $PBS_O_WORKDIR
blastall -p $PROGRAM -d $DB -i $IN $OPTIONS

This variant makes a couple changes. First, the script sets the PROGRAM and OPTIONS variables since they tend to be the same for many jobs. Secondly, the script will change directories not to a directory specified in the script, but rather to $PBS_O_WORKDIR which is the directory qsub was run in. Finally, blastall is run. As a general schema, this should be saved in an easy location such as ~/torque/blastp.sh where future schemes can go.

However, DB and IN were not specified in the script. Where do they come from? The answer lies in the second line. That is a directive telling qsub to send the DB and IN environmental variables to the remote host. These are set on the command line to qsub. However, they are environmental variable so they go before the command:

laurichj@bioinfo:~/blast/foo$ DB=/data/TIGR/arabidopsis/ATH1.pep IN=input.fa qsub ~/torque/blastp.sh

Which will run input.fa against the ATH database.

Many programs, such as blast run relativly quickly, but usually have a large number of input chunks. In the case of blast a chunk would be a single sequence. In general, a chunk is the smallest amount of input that can be ran independantly. In order to speed up the processing of a large number of blast, a little bit of shell scripting can be used. The idea is to split the data into a few smaller groups of chunks, then run these chunks as independant jobs. For example, a blast job of 1000 sequences can be speed up by running 5 jobs of 200 sequences.

The first step is to create the smaller input files. This is very dependant on the program, however almost any program that takes sequences as input can make use of seqsplit. This utility will take a sequence file and convert it into many smaller files. It can do this in one of two ways, first it can create files with a given number of sequences (with the last file possilby being smaller) or it can create a number of approximatly even files. To do this, you specify the input file, the number of sequences per file or the number of files, a prefix for the output files and a file format (if not fasta).

laurichj@bioinfo:~/blast/$ seqsplit -i proteins.fa -l 100 -p proteins-

This will create files of 100 sequences named proteins-01, proteins-02 and so on. In order to create 100 files, simply change the -l to -L, like so:

laurichj@bioinfo:~/blast/$ seqsplit -i proteins.fa -L 100 -p proteins-

That will create the files: proteins-01, proteins-02 to proteins-100. To make the file names uniform you can extend the length of the numeric suffix:

laurichj@bioinfo:~/blast/$ seqsplit -i proteins.fa -L 100 -s 3 -p proteins-

The additional -s 3 will make the files proteins-001, proteins-002 and so forth.

If each segment of the input has a constant number of lines or bytes, then split can be used to split the input. If the data is small enough, it can always be done by hand.

Now, we need to write a script file capable of taking any input we want. To do this, we are just going to use the blast.sh from the Reuse Examples, modified to only need the input files specified, and to operate on multiple files:

#!/bin/sh
#PBS -v IN

DB=/data/NCBI/blast/nr/nr
PROGRAM=blastp
OPTIONS="-e 1e-10"

cd $PBS_O_WORKDIR
for input in $IN; do
  blastall -p $PROGRAM -d $DB -i $input $OPTIONS
done

If only a few jobs need to be submitted, it will be simplist to do this by hand:

laurichj@bioinfo:~/blast/$ IN=proteins-01 qsub blast.sh

laurichj@bioinfo:~/blast/$ IN=proteins-02 qsub blast.sh

laurichj@bioinfo:~/blast/$ IN=proteins-03 qsub blast.sh

laurichj@bioinfo:~/blast/$ IN=proteins-04 qsub blast.sh

However, if there are dozens or hundreds of files, a rather complicated command line needs to be run. While this can be ran straight from the command line (write it all on one line), it will likely be easier to save it to a file, and make that file executable. First, we use find to create a list of filenames to run. Then we pipe that output into xargs to print a certain number of them per line, we finally send that into a small while loop in order to submit our jobs. The full command line will look like (formated for some readability):

find -name "proteins-*" | xargs -n 10 echo | while read IN; do
  export IN;
  qsub blast.sh;
done

That will list all files that start with "proteins-" and then print those out 10 to a line. For each line it will "export" the IN variable so that qsub can see it. qsub will pass that variable to the script. The result will be that a single job will be submited for every 10 input files.