Hello WDL ( Workflow Description Language )

This is a quick note about my first WDL workflow (Workflow Description Language) https://software.broadinstitute.org/wdl/.

As a Makefile, my workflow would be the following one:

NAME?=world
$(NAME)_sed.txt : $(NAME).txt
 sed 's/Hello/Goodbye/' $< > $@
$(NAME).txt:
 echo "Hello $(NAME)" > $@

Executed as:

$ make NAME=WORLD

echo "Hello WORLD" > WORLD.txt
sed 's/Hello/Goodbye/' WORLD.txt > WORLD_sed.txt

As a WDL, the workflow is described as below:

/**

print "Hello (message)" into a file

*/
task echotask {
 String outfile
 String ? message


 command {
  echo "Hello ${default="world" message}" > ${outfile}
  }
  
 output {
  File out = "${outfile}"
  }
}

/**

replace Hello with goodbye

*/
task sedtask {
 String outfile
 File infile

 
 command {
        sed 's/Hello/Goodbye/' ${infile}  > ${outfile}
  }
 
 output {
       File out = "${outfile}"
  }
}


workflow basicworkflow {
String message

call echotask {
 input: outfile="${message}.txt", message = "${message}"
 }
call sedtask {
 input: infile = echotask.out , outfile  = "${message}_sed.txt"
 }
}

I've downloaded cromwell-0.22.jar and wdltool-0.4.jar from github.


Generate a JSON describing the input parameters:

$ java -jar wdltool-0.4.jar inputs test.wdl  > inputs.json
$ cat inputs.json
{
  "basicworkflow.message": "String"
}

Customize the basicworkflow.message in "inputs.json":

$ cat inputs.json
{
  "basicworkflow.message": "WORLD"
}

Execute the workflow:

$ java -jar cromwell-0.22.jar run test.wdl inputs.json

[2016-10-27 12:36:29,69] [info] Slf4jLogger started
[2016-10-27 12:36:29,73] [info] RUN sub-command
[2016-10-27 12:36:29,73] [info]   WDL file: /home/lindenb/tmp/WDL/test.wdl
[2016-10-27 12:36:29,73] [info]   Inputs: /home/lindenb/tmp/WDL/inputs.json
[2016-10-27 12:36:29,76] [info] SingleWorkflowRunnerActor: Submitting workflow
[2016-10-27 12:36:30,43] [info] Running with database db.url = jdbc:hsqldb:mem:a402b714-c86f-414a-b03f-8e6ecb533e33;shutdown=false;hsqldb.tx=mvcc
[2016-10-27 12:36:37,67] [info] Metadata summary refreshing every 2 seconds.
[2016-10-27 12:36:38,20] [info] Workflow 172e1061-82fc-4240-9e80-e08aaafd6f3f submitted.
[2016-10-27 12:36:38,20] [info] SingleWorkflowRunnerActor: Workflow submitted 172e1061-82fc-4240-9e80-e08aaafd6f3f
[2016-10-27 12:36:38,38] [info] 1 new workflows fetched
[2016-10-27 12:36:38,38] [info] WorkflowManagerActor Starting workflow 172e1061-82fc-4240-9e80-e08aaafd6f3f
[2016-10-27 12:36:38,39] [info] WorkflowManagerActor Successfully started WorkflowActor-172e1061-82fc-4240-9e80-e08aaafd6f3f
[2016-10-27 12:36:38,39] [info] Retrieved 1 workflows from the WorkflowStoreActor
[2016-10-27 12:36:38,76] [info] MaterializeWorkflowDescriptorActor [172e1061]: Call-to-Backend assignments: basicworkflow.echotask -> Local, basicworkflow.sedtask -> Local
[2016-10-27 12:36:38,96] [info] WorkflowExecutionActor-172e1061-82fc-4240-9e80-e08aaafd6f3f [172e1061]: Starting calls: basicworkflow.echotask:NA:1
[2016-10-27 12:36:39,10] [info] SharedFileSystemAsyncJobExecutionActor [172e1061basicworkflow.echotask:NA:1]: echo "Hello WORLD" > WORLD.txt
[2016-10-27 12:36:39,10] [info] SharedFileSystemAsyncJobExecutionActor [172e1061basicworkflow.echotask:NA:1]: executing: /bin/bash /home/lindenb/tmp/WDL/cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-echotask/execution/script
[2016-10-27 12:36:39,11] [info] SharedFileSystemAsyncJobExecutionActor [172e1061basicworkflow.echotask:NA:1]: command: "/bin/bash" "/home/lindenb/tmp/WDL/cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-echotask/execution/script.submit"
[2016-10-27 12:36:39,14] [info] SharedFileSystemAsyncJobExecutionActor [172e1061basicworkflow.echotask:NA:1]: job id: 6056
[2016-10-27 12:36:39,20] [info] WorkflowExecutionActor-172e1061-82fc-4240-9e80-e08aaafd6f3f [172e1061]: Starting calls: basicworkflow.sedtask:NA:1
[2016-10-27 12:36:39,27] [info] SharedFileSystemAsyncJobExecutionActor [172e1061basicworkflow.sedtask:NA:1]: sed 's/Hello/Goodbye/' /home/lindenb/tmp/WDL/cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-sedtask/inputs/home/lindenb/tmp/WDL/cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-echotask/execution/WORLD.txt  > WORLD_sed.txt
[2016-10-27 12:36:39,27] [info] SharedFileSystemAsyncJobExecutionActor [172e1061basicworkflow.sedtask:NA:1]: executing: /bin/bash /home/lindenb/tmp/WDL/cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-sedtask/execution/script
[2016-10-27 12:36:39,27] [info] SharedFileSystemAsyncJobExecutionActor [172e1061basicworkflow.sedtask:NA:1]: command: "/bin/bash" "/home/lindenb/tmp/WDL/cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-sedtask/execution/script.submit"
[2016-10-27 12:36:39,28] [info] SharedFileSystemAsyncJobExecutionActor [172e1061basicworkflow.sedtask:NA:1]: job id: 6063
[2016-10-27 12:36:40,57] [info] WorkflowExecutionActor-172e1061-82fc-4240-9e80-e08aaafd6f3f [172e1061]: Workflow complete. Final Outputs: 
{
  "basicworkflow.echotask.out": "/home/lindenb/tmp/WDL/cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-echotask/execution/WORLD.txt",
  "basicworkflow.sedtask.out": "/home/lindenb/tmp/WDL/cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-sedtask/execution/WORLD_sed.txt"
}
[2016-10-27 12:36:40,61] [info] WorkflowManagerActor WorkflowActor-172e1061-82fc-4240-9e80-e08aaafd6f3f is in a terminal state: WorkflowSucceededState
{
  "outputs": {
    "basicworkflow.sedtask.out": "/home/lindenb/tmp/WDL/cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-sedtask/execution/WORLD_sed.txt",
    "basicworkflow.echotask.out": "/home/lindenb/tmp/WDL/cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-echotask/execution/WORLD.txt"
  },
  "id": "172e1061-82fc-4240-9e80-e08aaafd6f3f"
}
[2016-10-27 12:36:43,43] [info] SingleWorkflowRunnerActor workflow finished with status 'Succeeded'.

Check the final output:

$ cat ./cromwell-executions/basicworkflow/172e1061-82fc-4240-9e80-e08aaafd6f3f/call-sedtask/execution/WORLD_sed.txt
Goodbye WORLD

That's it,

Pierre

Divide-and-conquer in a #Makefile : recursivity and #parallelism.

This post is my notebook about implementing a divide-and-conquer strategy in GNU make.
Say you have a list of 'N' VCFs files. You want to create a list of:

• common SNPs in vcf1 and vcf2
• common SNPs in vcf3 and the previous list
• common SNPs in vcf4 and the previous list
• (...)
• common SNPs in vcfN and the previous list

Yes, I know I can do this using:grep -v '^#' f.vcf|cut -f 1,2,4,5 | sort | uniq

Using a linear Makefile it could look like:

list2: vcf1 vcf2
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@
list3: vcf3 list2
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@
list4: vcf4 list3
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@
list5: vcf5 list4
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@
(...)

We can speed-up the workflow using the parallel option of make -j (number-of-parallel-jobs) and using a divide-and-conquer strategy. Here, the targets 'list1_2' and 'list3_4' can be processed independently in parallel.

list1_2: vcf1 vcf2
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@

list3_4: vcf3 vcf4
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@

list1_4: list1_2 list3_4
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@

By using the internal 'Make' functions $(eval), $(call), $(shell) we can define a recursive method "recursive". This method takes two arguments which are the 0-based indexes of a VCF in the list of VCFs. Here is the Makefile:
Running the makefile:

$ make 
gunzip -c Sample18.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target0_1
gunzip -c Sample13.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target1_2
LC_ALL=C comm -12 target0_1 target1_2 > target0_2
gunzip -c Sample1.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target2_3
gunzip -c Sample19.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target3_4
gunzip -c Sample12.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target4_5
LC_ALL=C comm -12 target3_4 target4_5 > target3_5
LC_ALL=C comm -12 target2_3 target3_5 > target2_5
LC_ALL=C comm -12 target0_2 target2_5 > target0_5
gunzip -c Sample17.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target5_6
gunzip -c Sample16.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target6_7
LC_ALL=C comm -12 target5_6 target6_7 > target5_7
gunzip -c Sample9.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target7_8
gunzip -c Sample15.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target8_9
gunzip -c Sample5.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target9_10
LC_ALL=C comm -12 target8_9 target9_10 > target8_10
LC_ALL=C comm -12 target7_8 target8_10 > target7_10
LC_ALL=C comm -12 target5_7 target7_10 > target5_10
LC_ALL=C comm -12 target0_5 target5_10 > target0_10
gunzip -c Sample14.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target10_11
gunzip -c Sample3.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target11_12
LC_ALL=C comm -12 target10_11 target11_12 > target10_12
gunzip -c Sample11.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target12_13
gunzip -c Sample2.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target13_14
gunzip -c Sample6.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target14_15
LC_ALL=C comm -12 target13_14 target14_15 > target13_15
LC_ALL=C comm -12 target12_13 target13_15 > target12_15
LC_ALL=C comm -12 target10_12 target12_15 > target10_15
gunzip -c Sample20.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target15_16
gunzip -c Sample10.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target16_17
LC_ALL=C comm -12 target15_16 target16_17 > target15_17
gunzip -c Sample4.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target17_18
gunzip -c Sample8.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target18_19
gunzip -c Sample7.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target19_20
LC_ALL=C comm -12 target18_19 target19_20 > target18_20
LC_ALL=C comm -12 target17_18 target18_20 > target17_20
LC_ALL=C comm -12 target15_17 target17_20 > target15_20
LC_ALL=C comm -12 target10_15 target15_20 > target10_20
LC_ALL=C comm -12 target0_10 target10_20 > target0_20

and here is the generated workflow (drawn with make2graph ).
:
That's it
Pierre.

Parallelizing #RStats using #make

In the current post, I'll show how to use R as the main SHELL of GNU-Make instead of using a classical linux shell like 'bash'. Why would you do this ?

• awesomeness
• Make-based workflow management
• Make-based execution with --jobs. GNU make knows how to execute several recipes at once. Normally, make will execute only one recipe at a time, waiting for it to finish before executing the next. However, the '-j' or '--jobs' option tells make to execute many recipes simultaneously.

The following recipe has been tested with GNU-Make 4.0 and I'm not sure it would world with '<=3.81'.

The only problem is that R doesn't accept a multiline-argument on the command line (see http://stackoverflow.com/questions/21442674) so I created a wrapper 'mockR' that save the argument '-e "code"' into a file and pipe it into R:

(Edit1: A comment from madscientist : Re your script; you can save yourself some wear-and-tear on your disk and avoid the need for temp files and cleanup by just piping the input directly: echo "$R" | R --vanilla --no-readline --quiet . Just a thought. ")

(Edit2: the exit value of 'R' should also be returned by 'mockR'.)

This file is set as executable:

$ chmod u+x ./mockR

In the makefile, we tell 'make' to use 'mockR' instead of '/usr/bin/sh':

SHELL = ./mockR

The R code will be passed to 'mockR' using the argument '-e "code"'

.SHELLFLAGS= -e

We also set 'ONESHELL': "If .ONESHELL is mentioned as a target, then when a target is built all lines of the recipe will be given to a single invocation of the shell rather than each line being invoked separately"

.ONESHELL:

Example 1

We download the table 'knownGene' from the UCSC and we plot a pdf file 'countExons=f(txStart)'. Please, note that the targets are created using some R statements, NOT bash statements:

Now Invoke make

Example 2

Using a the eval and the call function we can make the previous 'Makefile' applicable for all the chromosomes:

Now Invoke make USING TRHEE PARALLEL JOBS

You can now watch the final pdf files:




That's it,
Pierre

GNU Make: saving the versions of the tools using 'order-only-prerequisites' : my notebook

Rule 3 of "Ten Simple Rules for Reproducible Computational Research". is
:

Archive the Exact Versions of All External Programs Used

.
I work with Makefile-based workflows: how can I save the version of each software used when I invoke 'make', whatever the target is ? A naive solution is to add a dependency to each target. For example, the following makefile takes a simple SAM file, convert it to BAM, sort and index. For each target, I've added a dependency named "dump_params" that append the version of samtools to a file "config.txt".

But that solution doesn't work because make re-builds all targets even if the top target already exists.

$ make

date << config.txt && \
 echo -n "Samtools " << config.txt && \
 samtools  2<&1 | grep Version << config.txt
samtools view -Sb samtools-0.1.18/examples/toy.sam < unsorted.bam
[samopen] SAM header is present: 2 sequences.
samtools sort unsorted.bam sorted
samtools index sorted.bam


$ make

date << config.txt && \
 echo -n "Samtools " << config.txt && \
 samtools  2<&1 | grep Version << config.txt
samtools view -Sb samtools-0.1.18/examples/toy.sam < unsorted.bam
[samopen] SAM header is present: 2 sequences.
samtools sort unsorted.bam sorted
samtools index sorted.bam

The solution I got via Stackoverflow is to use a order-only-prerequisites: "Order-only prerequisites can be specified by placing a pipe symbol (|) in the prerequisites list: any prerequisites to the left of the pipe symbol are normal; any prerequisites to the right are order-only... (...) Note that if you declare the same file to be both a normal and an order-only prerequisite, the normal prerequisite takes precedence (...)". The makefile with the 'order-only-prerequisites' is now:




And that works ! the final target is generated only once, but the file 'config.txt' is always generated.

$ make
date << config.txt && \
 echo -n "Samtools " << config.txt && \
 samtools  2<&1 | grep Version << config.txt
samtools view -Sb samtools-0.1.18/examples/toy.sam < unsorted.bam
[samopen] SAM header is present: 2 sequences.
samtools sort unsorted.bam sorted
samtools index sorted.bam

$ make
date << config.txt && \
 echo -n "Samtools " << config.txt && \
 samtools  2<&1 | grep Version << config.txt

$ make
date << config.txt && \
 echo -n "Samtools " << config.txt && \
 samtools  2<&1 | grep Version << config.txt

That's it,
Pierre

Update :another solution

Citing MadScientist's answer on stackoverflow : Another option is to use immediately expanded shell functions, like:

__dummy := $(shell echo "Makefile was run." >> config.txt)

Since it's immediately expanded the shell script will be invoked once, as the makefile is read in. There's no need to define a dump_params target or include it as a prerequisite. This is more old-school, but has the advantage that it will run for every invocation of make, without having to go through and ensure every target has the proper order-only prerequisite defined.

visualizing the dependencies in a Makefile

Update 2014: I wrote a C version at https://github.com/lindenb/makefile2graph.

I've just coded a tool to visualize the dependencies in a Makefile. The java source code is available on github at : https://github.com/lindenb/jsandbox/blob/master/src/sandbox/MakeGraphDependencies.java. This simple tool parses the ouput of

make -dq

( here option '-d' is 'Print lots of debugging information' and '-q' is 'Run no commands') and prints a graphiz-dot file.

Example

Below is a simple NGS workflow:

%.bam.bai : %.bam
 
file.vcf:  merged.bam.bai ref.fa
merged.bam : sorted1.bam sorted2.bam
sorted1.bam: lane1_1.fastq  lane1_2.fastq ref.fa
sorted2.bam: lane2_1.fastq  lane2_2.fastq ref.fa

Invoking the program:

make -d --dry-run | java -jar makegraphdependencies.jar

generates the following graphiz-dot file:

digraph G {
n9[label="sorted2.bam" ];
n3[label="merged.bam.bai" ];
n10[label="lane2_1.fastq" ];
n11[label="lane2_2.fastq" ];
n2[label="file.vcf" ];
n4[label="merged.bam" ];
n6[label="lane1_1.fastq" ];
n8[label="ref.fa" ];
n7[label="lane1_2.fastq" ];
n0[label="[ROOT]" ];
n5[label="sorted1.bam" ];
n1[label="Makefile" ];
n10->n9;
n11->n9;
n8->n9;
n4->n3;
n3->n2;
n8->n2;
n9->n4;
n5->n4;
n2->n0;
n1->n0;
n6->n5;
n8->n5;
n7->n5;
}

The result: (here using the google chart API for Graphviz)

That's it,
Pierre

Next Generation Sequencing, GNU-Make and .INTERMEDIATE

I gave a crash course about NGS to a few colleagues today. For my demonstration I wrote a simple Makefile. Basically, it downloads a subset of the human chromosome 22, indexes it with bwa, generates a set of fastqs with wgsim, align the fastqs, generates the *.sai, the *.sam, the *.bam, sorts the bam and calls the SNPs with mpileup.

SAMDIR=/usr/local/package/samtools-0.1.18
SAMTOOLS=$(SAMDIR)/samtools
BCFTOOLS=$(SAMDIR)/bcftools/bcftools
BWA=/usr/local/package/bwa-0.6.1/bwa

%.bam : %.sam
        $(SAMTOOLS) view -o $@ -b -S -T chr22.fa $<
%.bam.bai : %.bam
        $(SAMTOOLS) index $<

variations.vcf: variations.bcf
        $(BCFTOOLS) view $< > $@

variations.bcf : align.sorted.bam align.sorted.bam.bai
            $(SAMTOOLS) mpileup -uf chr22.fa $< | $(BCFTOOLS) view -bvcg - > $@


align.sam : random_1.sai random_2.sai  
        $(BWA) sampe chr22.fa $^ random_1.fq.gz random_2.fq.gz > $@

chr22.fa:
        curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr22.fa.gz" |\
        gunzip -c | grep -v NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN > $@


random_1.sai :  random_1.fq.gz chr22.fa.bwt
        $(BWA) aln -f $@ chr22.fa $<

random_2.sai :  random_2.fq.gz chr22.fa.bwt
        $(BWA) aln -f $@ chr22.fa $<

random_1.fq.gz random_2.fq.gz : chr22.fa
        $(SAMDIR)/misc/wgsim -N 1000000 $< random_1.fq random_2.fq > wgsim.output
        gzip --best random_1.fq random_2.fq

chr22.fa.bwt : chr22.fa
        $(BWA) index -a bwtsw $<

chr22.fa.fai :  chr22.fa
        $(SAMTOOLS) faidx $< 

align.sorted.bam : align.bam
        $(SAMTOOLS) sort $< align.sorted


clean:
        rm -f chr22.* *.bam *.vcf *.bcf *.sai *.gz *.fq *.bai  wgsim.output *.sam

I was asked about the weakness of "Make":
"The problem here, is that you need to generate a bunch of (possibly huge) files that won't be of use later: the *.sam, the unsorted *.bam, the *.sai, etc... if one of those files is removed everything will be re-compiled even if the final sorted.bam and the VCF exist."

I asked StackOverflow ( http://stackoverflow.com/questions/12199237 ) if a solution to fix this problem would exist and I received a nice solution from Eldar Abusalimov:

Use "intermediate files" feature of GNU Make:

Intermediate files are remade using their rules just like all other files. But intermediate files are treated differently in two ways.

The first difference is what happens if the intermediate file does not exist. If an ordinary file b does not exist, and make considers a target that depends on b, it invariably creates b and then updates the target from b. But if b is an intermediate file, then make can leave well enough alone. It won't bother updating b, or the ultimate target, unless some prerequisite of b is newer than that target or there is some other reason to update that target.

The second difference is that if make does create b in order to update something else, it deletes b later on after it is no longer needed. Therefore, an intermediate file which did not exist before make also does not exist after make. make reports the deletion to you by printing a rm -f command showing which file it is deleting.

Ordinarily, a file cannot be intermediate if it is mentioned in the makefile as a target or prerequisite. However, you can explicitly mark a file as intermediate by listing it as a prerequisite of the special target .INTERMEDIATE. This takes effect even if the file is mentioned explicitly in some other way.

You can prevent automatic deletion of an intermediate file by marking it as a secondary file. To do this, list it as a prerequisite of the special target .SECONDARY. When a file is secondary, make will not create the file merely because it does not already exist, but make does not automatically delete the file. Marking a file as secondary also marks it as intermediate.

So, adding the following line to the Makefile should be enough:

I've added the simple line below to my Makefile:

.INTERMEDIATE : align.sam random_1.sai random_2.sai align.bam variations.bcf

and I've invoked make:

[lindenb@srv-clc-02 20120830.demongs]$ make
curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr22.fa.gz" |\
    gunzip -c | grep -v N | head -n 100000 > chr22.fa
/usr/local/package/samtools-0.1.18/misc/wgsim -N 1000000 chr22.fa random_1.fq random_2.fq > wgsim.output
[wgsim_core] calculating the total length of the reference sequence...
[wgsim_core] 1 sequences, total length: 4999950
gzip -f --best random_1.fq random_2.fq
/usr/local/package/bwa-0.6.1/bwa index -a bwtsw chr22.fa
[bwa_index] Pack FASTA... 0.12 sec
[bwa_index] Construct BWT for the packed sequence...
[BWTIncCreate] textLength=9999900, availableWord=12703528
[bwt_gen] Finished constructing BWT in 5 iterations.
[bwa_index] 4.02 seconds elapse.
[bwa_index] Update BWT... 0.05 sec
[bwa_index] Pack forward-only FASTA... 0.08 sec
[bwa_index] Construct SA from BWT and Occ... 1.46 sec
[main] Version: 0.6.1-r104
[main] CMD: /usr/local/package/bwa-0.6.1/bwa index -a bwtsw chr22.fa
[main] Real time: 5.740 sec; CPU: 5.740 sec
/usr/local/package/bwa-0.6.1/bwa aln -f random_1.sai chr22.fa random_1.fq.gz
[bwa_aln] 17bp reads: max_diff = 2
[bwa_aln] 38bp reads: max_diff = 3
[bwa_aln] 64bp reads: max_diff = 4
[bwa_aln] 93bp reads: max_diff = 5
[bwa_aln] 124bp reads: max_diff = 6
[bwa_aln] 157bp reads: max_diff = 7
[bwa_aln] 190bp reads: max_diff = 8
[bwa_aln] 225bp reads: max_diff = 9
[bwa_aln_core] calculate SA coordinate... 26.07 sec
[bwa_aln_core] write to the disk... 0.08 sec
[bwa_aln_core] 262144 sequences have been processed.
[bwa_aln_core] calculate SA coordinate... 26.51 sec
[bwa_aln_core] write to the disk... 0.06 sec
[bwa_aln_core] 524288 sequences have been processed.
[bwa_aln_core] calculate SA coordinate... 26.80 sec
[bwa_aln_core] write to the disk... 0.08 sec
[bwa_aln_core] 786432 sequences have been processed.
[bwa_aln_core] calculate SA coordinate... 21.77 sec
[bwa_aln_core] write to the disk... 0.05 sec
[bwa_aln_core] 1000000 sequences have been processed.
[main] Version: 0.6.1-r104
[main] CMD: /usr/local/package/bwa-0.6.1/bwa aln -f random_1.sai chr22.fa random_1.fq.gz
[main] Real time: 104.702 sec; CPU: 104.675 sec
/usr/local/package/bwa-0.6.1/bwa aln -f random_2.sai chr22.fa random_2.fq.gz
[bwa_aln] 17bp reads: max_diff = 2
[bwa_aln] 38bp reads: max_diff = 3
[bwa_aln] 64bp reads: max_diff = 4
[bwa_aln] 93bp reads: max_diff = 5
[bwa_aln] 124bp reads: max_diff = 6
[bwa_aln] 157bp reads: max_diff = 7
[bwa_aln] 190bp reads: max_diff = 8
[bwa_aln] 225bp reads: max_diff = 9
[bwa_aln_core] calculate SA coordinate... 28.40 sec
[bwa_aln_core] write to the disk... 0.09 sec
[bwa_aln_core] 262144 sequences have been processed.
[bwa_aln_core] calculate SA coordinate... 28.94 sec
[bwa_aln_core] write to the disk... 0.07 sec
[bwa_aln_core] 524288 sequences have been processed.
[bwa_aln_core] calculate SA coordinate... 29.18 sec
[bwa_aln_core] write to the disk... 0.08 sec
[bwa_aln_core] 786432 sequences have been processed.
[bwa_aln_core] calculate SA coordinate... 23.07 sec
[bwa_aln_core] write to the disk... 0.05 sec
[bwa_aln_core] 1000000 sequences have been processed.
[main] Version: 0.6.1-r104
[main] CMD: /usr/local/package/bwa-0.6.1/bwa aln -f random_2.sai chr22.fa random_2.fq.gz
[main] Real time: 113.270 sec; CPU: 113.233 sec
/usr/local/package/bwa-0.6.1/bwa sampe chr22.fa random_1.sai random_2.sai random_1.fq.gz random_2.fq.gz > align.sam
[bwa_sai2sam_pe_core] convert to sequence coordinate...
[infer_isize] (25, 50, 75) percentile: (466, 500, 534)
[infer_isize] low and high boundaries: 330 and 670 for estimating avg and std
[infer_isize] inferred external isize from 220114 pairs: 499.899 +/- 49.771
[infer_isize] skewness: -0.006; kurtosis: -0.083; ap_prior: 1.00e-05
[infer_isize] inferred maximum insert size: 795 (5.93 sigma)
[bwa_sai2sam_pe_core] time elapses: 7.34 sec
[bwa_sai2sam_pe_core] changing coordinates of 6560 alignments.
[bwa_sai2sam_pe_core] align unmapped mate...
[bwa_paired_sw] 16796 out of 16796 Q17 singletons are mated.
[bwa_paired_sw] 27 out of 27 Q17 discordant pairs are fixed.
[bwa_sai2sam_pe_core] time elapses: 4.83 sec
[bwa_sai2sam_pe_core] refine gapped alignments... 0.97 sec
[bwa_sai2sam_pe_core] print alignments... 2.46 sec
[bwa_sai2sam_pe_core] 262144 sequences have been processed.
[bwa_sai2sam_pe_core] convert to sequence coordinate...
[infer_isize] (25, 50, 75) percentile: (466, 500, 534)
[infer_isize] low and high boundaries: 330 and 670 for estimating avg and std
[infer_isize] inferred external isize from 219840 pairs: 499.874 +/- 49.751
[infer_isize] skewness: 0.001; kurtosis: -0.072; ap_prior: 1.00e-05
[infer_isize] inferred maximum insert size: 795 (5.93 sigma)
[bwa_sai2sam_pe_core] time elapses: 7.36 sec
[bwa_sai2sam_pe_core] changing coordinates of 6668 alignments.
[bwa_sai2sam_pe_core] align unmapped mate...
[bwa_paired_sw] 16833 out of 16834 Q17 singletons are mated.
[bwa_paired_sw] 38 out of 38 Q17 discordant pairs are fixed.
[bwa_sai2sam_pe_core] time elapses: 4.78 sec
[bwa_sai2sam_pe_core] refine gapped alignments... 0.98 sec
[bwa_sai2sam_pe_core] print alignments... 2.55 sec
[bwa_sai2sam_pe_core] 524288 sequences have been processed.
[bwa_sai2sam_pe_core] convert to sequence coordinate...
[infer_isize] (25, 50, 75) percentile: (466, 500, 534)
[infer_isize] low and high boundaries: 330 and 670 for estimating avg and std
[infer_isize] inferred external isize from 220140 pairs: 500.062 +/- 49.780
[infer_isize] skewness: -0.000; kurtosis: -0.075; ap_prior: 1.00e-05
[infer_isize] inferred maximum insert size: 795 (5.93 sigma)
[bwa_sai2sam_pe_core] time elapses: 7.76 sec
[bwa_sai2sam_pe_core] changing coordinates of 6522 alignments.
[bwa_sai2sam_pe_core] align unmapped mate...
[bwa_paired_sw] 16804 out of 16806 Q17 singletons are mated.
[bwa_paired_sw] 33 out of 33 Q17 discordant pairs are fixed.
[bwa_sai2sam_pe_core] time elapses: 4.79 sec
[bwa_sai2sam_pe_core] refine gapped alignments... 1.07 sec
[bwa_sai2sam_pe_core] print alignments... 2.57 sec
[bwa_sai2sam_pe_core] 786432 sequences have been processed.
[bwa_sai2sam_pe_core] convert to sequence coordinate...
[infer_isize] (25, 50, 75) percentile: (466, 500, 534)
[infer_isize] low and high boundaries: 330 and 670 for estimating avg and std
[infer_isize] inferred external isize from 179161 pairs: 499.910 +/- 49.860
[infer_isize] skewness: -0.001; kurtosis: -0.075; ap_prior: 1.00e-05
[infer_isize] inferred maximum insert size: 796 (5.93 sigma)
[bwa_sai2sam_pe_core] time elapses: 6.22 sec
[bwa_sai2sam_pe_core] changing coordinates of 5351 alignments.
[bwa_sai2sam_pe_core] align unmapped mate...
[bwa_paired_sw] 13904 out of 13905 Q17 singletons are mated.
[bwa_paired_sw] 27 out of 27 Q17 discordant pairs are fixed.
[bwa_sai2sam_pe_core] time elapses: 3.97 sec
[bwa_sai2sam_pe_core] refine gapped alignments... 0.86 sec
[bwa_sai2sam_pe_core] print alignments... 2.10 sec
[bwa_sai2sam_pe_core] 1000000 sequences have been processed.
[main] Version: 0.6.1-r104
[main] CMD: /usr/local/package/bwa-0.6.1/bwa sampe chr22.fa random_1.sai random_2.sai random_1.fq.gz random_2.fq.gz
[main] Real time: 70.067 sec; CPU: 69.984 sec
/usr/local/package/samtools-0.1.18/samtools view -o align.bam -b -S -T chr22.fa align.sam
[samopen] SAM header is present: 1 sequences.
/usr/local/package/samtools-0.1.18/samtools sort align.bam align.sorted
/usr/local/package/samtools-0.1.18/samtools index align.sorted.bam
/usr/local/package/samtools-0.1.18/samtools mpileup -uf chr22.fa align.sorted.bam | /usr/local/package/samtools-0.1.18/bcftools/bcftools view -bvcg - > variations.bcf
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
[bcfview] 100000 sites processed.
[afs] 0:99762.176 1:154.256 2:83.568
[bcfview] 200000 sites processed.
[afs] 0:99790.512 1:146.516 2:62.972
[bcfview] 300000 sites processed.
[afs] 0:99750.990 1:143.270 2:105.740
[bcfview] 400000 sites processed.
[afs] 0:99764.146 1:156.783 2:79.071
[bcfview] 500000 sites processed.
[afs] 0:99749.220 1:177.647 2:73.133
[bcfview] 600000 sites processed.
[afs] 0:99758.419 1:160.146 2:81.435
[bcfview] 700000 sites processed.
[afs] 0:99769.451 1:155.968 2:74.581
[bcfview] 800000 sites processed.
[afs] 0:99761.914 1:149.704 2:88.382
[bcfview] 900000 sites processed.
[afs] 0:99746.200 1:158.466 2:95.334
[bcfview] 1000000 sites processed.
[afs] 0:99732.973 1:177.167 2:89.860
[bcfview] 1100000 sites processed.
[afs] 0:99815.394 1:109.322 2:75.284
[bcfview] 1200000 sites processed.
[afs] 0:99762.584 1:160.194 2:77.222
[bcfview] 1300000 sites processed.
[afs] 0:99748.201 1:155.610 2:96.189
[bcfview] 1400000 sites processed.
[afs] 0:99739.710 1:170.437 2:89.853
[bcfview] 1500000 sites processed.
[afs] 0:99729.049 1:176.956 2:93.995
[bcfview] 1600000 sites processed.
[afs] 0:99747.410 1:163.200 2:89.390
[bcfview] 1700000 sites processed.
[afs] 0:99751.841 1:171.932 2:76.228
[bcfview] 1800000 sites processed.
[afs] 0:99800.303 1:145.575 2:54.122
[bcfview] 1900000 sites processed.
[afs] 0:99828.697 1:126.069 2:45.234
[bcfview] 2000000 sites processed.
[afs] 0:99730.074 1:158.242 2:111.684
[afs] 0:87928.796 1:106.271 2:97.933
/usr/local/package/samtools-0.1.18/bcftools/bcftools view variations.bcf > variations.vcf
rm random_1.sai variations.bcf random_2.sai align.bam align.sam

here is the last line of the output:

rm random_1.sai variations.bcf random_2.sai align.bam align.sam

Makefile has removed the intermediate files
The working directory only contain the final files:

$ ls -la
total 168020
drwxr-xr-x  2 lindenb users    4096 Aug 30 17:46 .
drwxr-xr-x 12 lindenb users    4096 Aug 30 14:18 ..
-rw-r--r--  1 lindenb users 81537078 Aug 30 17:43 align.sorted.bam
-rw-r--r--  1 lindenb users    15096 Aug 30 17:43 align.sorted.bam.bai
-rw-r--r--  1 lindenb users  5099956 Aug 30 17:36 chr22.fa
-rw-r--r--  1 lindenb users      181 Aug 30 17:37 chr22.fa.amb
-rw-r--r--  1 lindenb users      41 Aug 30 17:37 chr22.fa.ann
-rw-r--r--  1 lindenb users  5000048 Aug 30 17:37 chr22.fa.bwt
-rw-r--r--  1 lindenb users      22 Aug 30 17:42 chr22.fa.fai
-rw-r--r--  1 lindenb users  1249989 Aug 30 17:37 chr22.fa.pac
-rw-r--r--  1 lindenb users  2500024 Aug 30 17:37 chr22.fa.sa
-rw-r--r--  1 lindenb users    1344 Aug 30 17:35 Makefile
-rw-r--r--  1 lindenb users 37831888 Aug 30 17:36 random_1.fq.gz
-rw-r--r--  1 lindenb users 37836209 Aug 30 17:36 random_2.fq.gz
-rw-r--r--  1 lindenb users  611760 Aug 30 17:46 variations.vcf
-rw-r--r--  1 lindenb users  101626 Aug 30 17:36 wgsim.output

That's it,

Pierre