Makefile.angsd

#-*- makefile -*-

MAKEDIR = $(dir $(lastword $(MAKEFILE_LIST)))
ifeq ($(findstring Makefile.ngsvars,$(MAKEFILE_LIST)),)
include $(MAKEDIR)Makefile.ngsvars
endif
ifeq ($(findstring Makefile.general,$(MAKEFILE_LIST)),)
include $(MAKEDIR)Makefile.general
endif

# Angsd search path. Since some of the programs are located in the
# misc subfolder, use absolute search path. The user must set this
# variable in the makefile that includes this file, or else set the
# executable variables below.
ifndef ANGSD_PATH
ANGSD_PATH=.
endif

# Angsd executable
ifndef ANGSD
ANGSD=$(ANGSD_PATH)/angsd
endif
# emoptim executable
ifndef ANGSD_EMOPTIM
ANGSD_EMOPTIM=$(ANGSD_PATH)/misc/emOptim
endif
# emoptim2 executable
ifndef ANGSD_EMOPTIM2
ANGSD_EMOPTIM2=$(ANGSD_PATH)/misc/emOptim2
endif
# optimSFS executable
ifndef ANGSD_OPTIMSFS
ANGSD_OPTIMSFS=$(ANGSD_PATH)/misc/optimSFS
endif
# bgid executable
ifndef ANGSD_BGID
ANGSD_BGID=$(ANGSD_PATH)/misc/bgid
endif
# calcstat executable
ifndef ANGSD_CALCSTAT
ANGSD_CALCSTAT=$(ANGSD_PATH)/misc/calcStat
endif
ifndef ANGSD_OPTION_THREADS
ANGSD_OPTION_THREADS=$(THREADS)
endif
##############################
# Program options
##############################
# GL option - use GATK model by default
ifndef ANGSD_OPTION_GL
ANGSD_OPTION_GL=2
endif
# Maf option - use option 10 by default (EM unknown + EM known minor)
ifndef ANGSD_OPTION_MAF
ANGSD_OPTION_MAF=10
endif
# MajorMinor - 5 default (major is ancestral, minor inferred from Genotype Likelihoods)
#
# NB: if the ancestral is not in the inferred (major/minor) tuple, the
# site is discarded
ifndef ANGSD_OPTION_MAJORMINOR
ANGSD_OPTION_MAJORMINOR=5
endif
# Option realsfs
ifndef ANGSD_OPTION_REALSFS
ANGSD_OPTION_REALSFS=1
endif
# Ancestral file
# Set to empty here; will throw error
ifndef ANGSD_OPTION_ANC
ANGSD_OPTION_ANC=$(REF)
endif
ifndef ANGSD_OPTION_REF
ANGSD_OPTION_REF=$(REF)
endif

# Window sizes and steps for thetas calculations
ifndef ANGSD_OPTION_WINDOW
ANGSD_OPTION_WINDOW=50000
endif

ifndef ANGSD_OPTION_STEP
ANGSD_OPTION_STEP=50000
endif

ifndef ANGSD_OPTION_MAXITER
ANGSD_OPTION_MAXITER=100
endif

ifndef ANGSD_OPTION_NCHR
ANGSD_OPTION_NCHR=
endif

ifndef ANGSD_OPTION_NSITES
ANGSD_OPTION_NSITES=
endif

ifndef ANGSD_OPTION_CHR
ANGSD_OPTION_CHR=
endif

##############################
# Utility variables
# 
# These are not angsd options but are used by the makefile
##############################
# String used to find chromosome label
ifndef ANGSD_CHR_LABEL
ANGSD_CHR_LABEL=scaffold
endif
ifndef ANGSD_CHROMOSOMES_EXTRA
ANGSD_CHROMOSOMES_EXTRA=chrM
endif
ifndef ANGSD_CHROMOSOMES_WIDTH
ANGSD_CHROMOSOMES_WIDTH=
endif

# Utility for generating chromosome names; if ANGSD_CHR_ENUM_FROM and
# ANGSD_CHR_ENUM_TO is set, enumerate chromosomes
ifndef ANGSD_CHR_ENUM_FROM
ANGSD_CHR_ENUM_FROM=1
endif
ifndef ANGSD_CHR_ENUM_TO
ANGSD_CHR_ENUM_TO=1
endif
# Enumration of chromosomes
ifndef ANGSD_CHR_ENUM
ANGSD_CHR_ENUM=$(shell seq $(ANGSD_CHR_ENUM_FROM) $(ANGSD_CHR_ENUM_TO))
endif
ifndef ANGSD_CHROMOSOMES
ANGSD_CHROMOSOMES=$(foreach i,${ANGSD_CHR_ENUM},$(ANGSD_CHR_LABEL)_$(i)) $(ANGSD_CHROMOSOMES_EXTRA)
endif
ifndef ANGSD_POPULATIONS
ANGSD_POPULATIONS=
endif


# Precious files
.PRECIOUS: %.em.ml %.em2.ml %_full.list %.list

# Create SFS file
#
# Requirement: here a bam list text file (i.e. names and locations of
# bam files for population)
ifndef ANGSD_REALSFS_OPTIONS
ANGSD_REALSFS_OPTIONS=-realSFS 1 -GL $(ANGSD_OPTION_GL) -anc $(ANGSD_OPTION_ANC) -nThreads $(ANGSD_OPTION_THREADS)
endif
%.sfs: %.list
	$(ANGSD) $(ANGSD_REALSFS_OPTIONS) -bam $< -out $(@:.sfs=).tmp && mv $(@:.sfs=).tmp.sfs $@

# Several rules for generating ML-estimate of sfs. Need to set a sfs optimization program variable
# 1. emOptim2
ifndef ANGSD_EMOPTIM2_OPTIONS
ANGSD_EMOPTIM2_OPTIONS=-maxIter $(ANGSD_OPTION_MAXITER) -P $(ANGSD_OPTION_THREADS)
endif
ifneq ($(ANGSD_OPTION_NSITES),)
ANGSD_EMOPTIM2_OPTIONS+=-nSites $(ANGSD_OPTION_NSITES)
endif
ifdef ANGSD_OPTION_TOLE
ANGSD_EMOPTIM2_OPTIONS+=-tole $(ANGSD_OPTION_TOLE)
endif
ifdef ANGSD_EMOPTIM2_OPTION_FILTER
ANGSD_EMOPTIM2_OPTIONS+=-filter $(ANGSD_EMOPTIM2_OPTION_FILTER)
endif

%.sfs.em2.ml: %.sfs
	$(ANGSD_EMOPTIM2) $< $(ANGSD_OPTION_NCHR) $(ANGSD_EMOPTIM2_OPTIONS) > $@.exp
	$(AWK) '{for (i=1; i<NF; i++) printf("%f ",exp($$i)); printf("%f", exp($$i))}' $@.exp > $@.tmp && mv $@.tmp $@

# 2. emOptim
ifndef ANGSD_EMOPTIM_OPTIONS
ANGSD_EMOPTIM_OPTIONS=-maxIter $(ANGSD_OPTION_MAXITER)
endif
%.sfs.em.ml: %.sfs
	$(ANGSD_EMOPTIM) -binput $< -nChr $(ANGSD_OPTION_NCHR) $(ANGSD_EMOPTIM_OPTIONS) > $@.tmp && mv $@.tmp $@

# 3. optimSFS
ifdef ANGSD_OPTIMSFS_OPTIONS
ANGSD_OPTIMSFS_OPTIONS=
endif
%.sfs.ml: %.sfs
	$(ANGSD_OPTIMSFS) -binput $< -nChr $(ANGSD_OPTION_NCHR) -outnames $@.tmp && mv $@.tmp $@

# doThetas
# Generate thetas files
# NB: by default this is done chromosome by chromosome

# The recipe assumes files are named %_CHROMOSOME.thetas.gz, where
# CHROMOSOME can contain one separating underscore '_'
ifndef ANGSD_DOTHETAS_OPTIONS
ANGSD_DOTHETAS_OPTIONS=-realSFS 1 -GL $(ANGSD_OPTION_GL) -anc $(ANGSD_OPTION_ANC) -doThetas 1 -nThreads $(ANGSD_OPTION_THREADS) -r $(ANGSD_OPTION_CHR)
endif


# additional rules that rely on links having been set up
%.thetas.gz: %.sfs.em2.ml.ln
	$(if $(findstring $(ANGSD_CHR_LABEL),$(notdir $*)),$(eval ANGSD_OPTION_CHR=$(word 2,$(subst _, ,$(notdir $*)))_$(word 3,$(subst _, ,$(notdir $*)))),$(eval ANGSD_OPTION_CHR=$(word 2,$(subst _, ,$(notdir $*)))))
	$(ANGSD) $(ANGSD_DOTHETAS_OPTIONS) -pest $(<:.ln=) -bam $*.list -out $*.tmp && mv $*.tmp.thetas.gz $@ && mv $*.tmp.arg $*.arg

# Regular rules
%.thetas.gz: %.sfs.em2.ml
	$(if $(findstring $(ANGSD_CHR_LABEL),$(notdir $*)),$(eval ANGSD_OPTION_CHR=$(word 2,$(subst _, ,$(notdir $*)))_$(word 3,$(subst _, ,$(notdir $*)))),$(eval ANGSD_OPTION_CHR=$(word 2,$(subst _, ,$(notdir $*)))))
	$(ANGSD) $(ANGSD_DOTHETAS_OPTIONS) -pest $< -bam $*.list -out $*.tmp && mv $*.tmp.thetas.gz $@ && mv $*.tmp.arg $*.arg

%.thetas.gz: %.sfs.em.ml
	$(if $(findstring $(ANGSD_CHR_LABEL),$(notdir $*)),$(eval ANGSD_OPTION_CHR=$(word 2,$(subst _, ,$(notdir $*)))_$(word 3,$(subst _, ,$(notdir $*)))),$(eval ANGSD_OPTION_CHR=$(word 2,$(subst _, ,$(notdir $*)))))
	$(ANGSD) $(ANGSD_DOTHETAS_OPTIONS) -pest $< -bam $*.list -out $*.tmp && mv $*.tmp.thetas.gz $@ && mv $*.tmp.arg $*.arg

# doMafs
ifndef ANGSD_DOMAFS_OPTIONS
ANGSD_DOMAFS_OPTIONS=-anc $(ANGSD_OPTION_ANC) -doMajorMinor $(ANGSD_OPTION_MAJORMINOR) -doMaf $(ANGSD_OPTION_MAF) -GL $(ANGSD_OPTION_GL) -nThreads $(ANGSD_OPTION_THREADS)
endif
ifneq ($(ANGSD_OPTION_CHR),)
ANGSD_DOMAFS_OPTIONS+=-r $(ANGSD_OPTION_CHR)
endif
ifdef ANGSD_DOMAFS_OPTION_FILTER
ANGSD_DOMAFS_OPTIONS+=-filter $(ANGSD_DOMAFS_OPTION_FILTER)
endif

%.mafs.gz: %.list
	$(ANGSD) $(ANGSD_DOMAFS_OPTIONS) -bam $< -out $(@:.mafs.gz=).tmp && rename $(@:.mafs.gz=).tmp $(@:.mafs.gz=) $(@:.mafs.gz=).tmp*

# doCounts: Generate mpileup from show options
ifndef ANGSD_DOCOUNTS_OPTIONS
ANGSD_DOCOUNTS_OPTIONS=-doCounts 1 -dumpCounts 2 -anc $(ANGSD_OPTION_ANC)
endif
%.counts.gz: %.list
	$(ANGSD) $(ANGSD_DOCOUNTS_OPTIONS) -bam $< -out $(@:.sfs=.counts.gz).tmp && rename $(@:.counts.gz=).tmp $(@:.counts.gz=) $(@:.counts.gz=).tmp*

# add -show to get mpileup
%.mpileup: %.list
	$(ANGSD) $(ANGSD_DOCOUNTS_OPTIONS) -show 1 -bam $< -out $(@:.sfs=.mpileup).tmp && mv $(@:.mpileup=).tmp.mpileup $@ 

# bgid
%.thetas.gz.bin: %.thetas.gz
	$(ANGSD_BGID) make_bed $<

ifndef ANGSD_BGID_DOSTAT_OPTIONS
ANGSD_BGID_DOSTAT_OPTIONS=-step $(ANGSD_OPTION_STEP) -win $(ANGSD_OPTION_WINDOW) -nChr $(ANGSD_OPTION_NCHR)
endif
%.thetas.gz.bgid_pestPG: %.thetas.gz %.thetas.gz.bin
	$(ANGSD_BGID) do_stat $< $(ANGSD_BGID_DOSTAT_OPTIONS) && mv $*.thetas.gz.pestPG $@

# Bed file is zero based; hence subtract one from first coordinate
%.bed: %.mafs.gz
	$(ZCAT) $< | $(GREP) -v chromo | $(AWK) '{print $$1,$$2-1,$$2}' > $@.tmp && mv $@.tmp $@

%.pos.gz: %.bed
	$(AWK) '{print $$1, $$3}' $< | $(GZIP) -c > $@.tmp && mv $@.tmp $@

%.pos: %.bed
	$(AWK) '{print $$1, $$3}' $< > $@.tmp && mv $@.tmp $@


##############################
# ngsTools
##############################
ifndef NGSTOOLS_PATH
NGSTOOLS_PATH=.
endif
ifndef NGSTOOLS_NGSFST
NGSTOOLS_NGSFST=$(NGSTOOLS_PATH)/ngsFST
endif
ifndef NGSTOOLS_OPTION_NSITES
NGSTOOLS_SITES=
endif
ifndef NGSTOOLS_OPTION_BLOCK_SIZE
NGSTOOLS_OPTION_BLOCK_SIZE=
endif

##############################
# GEM
##############################
ifndef GEM_PATH
GEM_PATH=.
endif

# This is needed as GEM internally calls GEM programs
export PATH := $(GEM_PATH):$(PATH)

ifndef GEM_INDEXER_OPTIONS
GEM_INDEXER_OPTIONS=-T 7
endif
%.gem: %.fasta
	$(GEM_PATH)/gem-indexer $(GEM_INDEXER_OPTIONS) -i $< -o $*.tmp && mv $*.tmp.gem $@; mv $*.tmp.log $*.log

ifndef GEM_MAPPABILITY_OPTIONS
GEM_MAPPABILITY_OPTIONS=-l 100 -T 7
endif
%.mappability: %.gem
	$(GEM_PATH)/gem-mappability $(GEM_MAPPABILITY_OPTIONS) -I $< -o $*.tmp && mv $*.tmp.mappability $@

##############################
# ngstk utilities
# see https://github.com/percyfal/ngstk
##############################
ifndef NGSTK
NGSTK=ngstk
endif
# Naming convention: POPA_POPB_CHROMOSOME_FRACTION_COVERAGE_ngstk.txt
%_ngstk.txt:
	$(eval POPA=$(word 1,$(subst _, ,$*)))
	$(eval POPB=$(word 2,$(subst _, ,$*)))
	$(eval CHR=$(word 3,$(subst _, ,$*)))
	$(if $(findstring $(ANGSD_CHR_LABEL),$(notdir $*)),$(eval CHR=$(word 3,$(subst _, ,$(notdir $*)))_$(word 4,$(subst _, ,$(notdir $*)))),$(eval CHR=$(word 3,$(subst _, ,$(notdir $*)))))
	$(if $(findstring $(ANGSD_CHR_LABEL),$(notdir $*)),$(eval FRAC=$(subst F,,$(word 5,$(subst _, ,$(notdir $*))))),$(eval FRAC=$(subst F,,$(word 4,$(subst _, ,$(notdir $*))))))
	$(if $(findstring $(ANGSD_CHR_LABEL),$(notdir $*)),$(eval COV=$(subst C,,$(word 6,$(subst _, ,$(notdir $*))))),$(eval COV=$(subst C,,$(word 5,$(subst _, ,$(notdir $*))))))
	$(NGSTK) angsd $(POPA)_$(CHR) $(POPB)_$(CHR) -f $(FRAC) -c $(COV) > $@.tmp && mv $@.tmp $@

##############################
# Plotting utilities
##############################
ifndef PLOTTHETAS
PLOTTHETAS=$(MAKEDIR)scripts/plotThetas.R
endif
ifndef PLOTNGSTK
PLOTNGSTK=$(MAKEDIR)scripts/plotNgstk.R
endif

%.pdf: %.bgid_pestPG
	$(PLOTTHETAS) $< $@.tmp && mv $@.tmp $@

%_ngstk.pdf: %_ngstk.txt
	$(PLOTNGSTK) $< $@.tmp && mv $@.tmp $@

##############################
# settings
##############################
.PHONY: angsd-settings angsd-header angsd-general-options angsd-program-options utility-programs

print-%:
	@echo '$*=$($*)'

angsd-settings: angsd-header print-ANGSD_PATH print-ANGSD print-ANGSD_EMOPTIM print-ANGSD_EMOPTIM2 print-ANGSD_OPTIMSFS print-ANGSD_BGID print-ANGSD_CALCSTAT angsd-general-options print-ANGSD_OPTION_THREADS print-ANGSD_OPTION_GL print-ANGSD_OPTION_MAF print-ANGSD_OPTION_MAJORMINOR print-ANGSD_OPTION_REALSFS print-ANGSD_OPTION_ANC print-ANGSD_OPTION_REF print-ANGSD_OPTION_WINDOW print-ANGSD_OPTION_STEP print-ANGSD_OPTION_MAXITER print-ANGSD_OPTION_NCHR print-ANGSD_OPTION_NSITES print-ANGSD_OPTION_CHR angsd-program-options print-ANGSD_REALSFS_OPTIONS print-ANGSD_EMOPTIM2_OPTIONS print-ANGSD_EMOPTIM_OPTIONS print-ANGSD_OPTIMSFS_OPTIONS print-ANGSD_DOTHETAS_OPTIONS print-ANGSD_DOMAFS_OPTIONS print-ANGSD_DOCOUNTS_OPTIONS print-ANGSD_BGID_DOSTAT_OPTIONS angsd-utility-vars  print-ANGSD_POPULATIONS  print-ANGSD_CHR_ENUM_FROM print-ANGSD_CHR_ENUM_TO print-ANGSD_CHR_ENUM print-ANGSD_CHR_LABEL print-ANGSD_CHROMOSOMES_EXTRA utility-programs print-NGSTOOLS_PATH print-NGSTOOLS_NGSFST print-NGSTOOLS_OPTION_BLOCK_SIZE print-NGSTOOLS_OPTION_NSITES print-GEM_PATH print-GEM_INDEXER_OPTIONS print-GEM_MAPPABILITY_OPTIONS print-NGSTK

angsd-header:
	@echo -e "\nMakefile.angsd options"
	@echo "========================"

angsd-general-options:
	@echo -e "\nGeneral options"
	@echo "------------------------"

angsd-program-options:
	@echo -e "\nProgram options"
	@echo "------------------------"

angsd-utility-vars:
	@echo -e "\nAngsd utility variables"
	@echo "------------------------------"

utility-programs: 
	@echo -e "\nUtility programs"
	@echo "========================"

angsd-target-names-header: 
	@echo -e "\nAngsd targets"
	@echo "------------------------"

angsd-target-names: angsd-target-names-header print-THETASTARGETS print-THETASBINTARGETS