Example makefile for data processing

Makefile

# Makefile for sequence analysis

D3 = 10229609_627E8AAXX_s_1.barcode.G1366_D3.q20l50.fq
D4 = 10229609_627E8AAXX_s_1.barcode.G1300_D4.q20l50.fq
T2 = 10229609_627E8AAXX_s_1.barcode.G1393_T2.q20l50.fq

DIPLOIDS = $(D3) $(D4)
SEQUENCES = $(DIPLOIDS) $(T2)

TEST = test.fq

.PHONY: all
all: called.vcf

.PHONY: test
test: $(TEST:.fq=.bam)

.PHONY: clean
clean: index-clean align-clean snp-clean

MAXDIFF = 0.04
# number of threads to use in alignment
THREADS = 16
# minimum score for bwasw alignment
THRESH = 37

# Gmax_109.fa from http://www.plantgdb.org/XGDB/phplib/download.php?GDB=Gm
GMAX = Gmax_109.fa
INDEX = $(GMAX).bwt

.PHONY: index
index: $(INDEX)

$(INDEX): $(GMAX)
	bwa index -a bwtsw $<

.PRECIOUS: $(INDEX)
.PRECIOUS: $(SEQUENCES:.fq=.sai)
.PRECIOUS: $(SEQUENCES:.fq=.sam)
.PRECIOUS: $(TEST:.fq=.sam)

%.sai: %.fq $(INDEX)
	bwa aln -t $(THREADS) -I $(GMAX) $< > $@

#%.sam: %.sai $(INDEX)
#	bwa samse -r "@RG\tID:$<\tPL:ILLUMINA" $(GMAX) $< $*.fq > $@

%.sam: %.fq $(INDEX_BWTSW)
	bwa bwasw -t $(THREADS) -T $(THRESH) $(GMAX) $< > $@

.PRECIOUS: $(SEQUENCES:.fq=.bam)
.PRECIOUS: $(TEST:.fq=.bam)

# mpileup expects sorted bam files, so sort them here
%.bam: %.sam
	samtools view -b -h -S -u $< | samtools sort - $*

%.fai: %.fa
	samtools faidx $<
	mv $<.fai $@

snps.bcf: $(DIPLOIDS:.fq=.bam) $(GMAX:.fq=.fai)
	samtools mpileup -DABSug -f $(GMAX:.fq=.fai) \
		$(DIPLOIDS:.fq=.bam) | tee snps.pileup > snps.bcf
# -C 50

called.vcf: snps.bcf
	bcftools view -s samples.samples -t pair -g snps.bcf > called.vcf

.PHONY: index-clean
index-clean:
	rm -f *.amb *.ann *.pac *.rpac *.bwt *.rbwt *.sa *.rsa *.fai

.PHONY: align-clean
align-clean:
	rm -f *.sam *.bam

.PHONY: snp-clean
snp-clean:
	rm -f *.bcf *.vcf *.pileup