ntJoin

#!/usr/bin/make -rRf

# ntJoin: Scaffold assemblies using reference assemblies and minimizer graphs
# Written by Lauren Coombe @lcoombe
# ntJoin v1.1.0

# Input files
target=None
references=None
reference_config=None

# List of edge weights
reference_weights=None
target_weight=1

# Prepare targets
ifneq ($(reference_config), None)
references_list=$(shell cat $(reference_config) | awk -F"," '{print $$1}')
reference_weights_list=$(shell cat $(reference_config) | awk -F"," '{print $$2}')
else
references_list=$(references)
reference_weights_list=$(reference_weights)
endif

min_ref_targets=$(addsuffix .k$(k).w$(w).tsv, $(references_list))
fai_ref_targets=$(addsuffix .fai, $(references_list))

# Path to minimizer assemble code
assemble_path=$(shell dirname $(realpath $(MAKEFILE_LIST)))

# Window size
w=1000

# Kmer size
k=32

# Run overlap logic?
overlap=True

# Window size for overlap detection
overlap_w=10

# Kmer size for overlap detection
overlap_k=15

# Number of threads
t=4

# Number of threads for ntjoin_assemble.py (Change with caution -- multi-processing increases memory)
assemble_t=1

# Minimum edge weight
n=1

# Minimum gap size between scaffolds
g=20

# Gap size between trimmed overlapping segments (used if overlap=True)
overlap_g?=$g

# Maximum gap size between scaffolds (0 if no maximum)
G=0

# Whether Mann-Kendall test will be run to orient contigs (True) or not (False)
mkt=False

# Output AGP file?
agp=False

# The proportion of increasing/decreasing minimizer positions needed to orient a contig
m=90

# Reference genome (for analysis)
ref=None

# Run Quast in large mode (1) or not (0)
large=0

# Cut the input contigs and putative misassemblies?
no_cut=False

# Log time for each step?
time=False

# Prefix for output files
prefix=out.k$(k).w$(w).n$(n)

SHELL=bash -e -o pipefail
ifeq ($(shell zsh -e -o pipefail -c 'true' 2>/dev/null; echo $$?), 0)
# Set pipefail to ensure that all commands of a pipe succeed.
SHELL=zsh -e -o pipefail
# Report run time and memory usage with zsh.
export REPORTTIME=1
export TIMEFMT=time user=%U system=%S elapsed=%E cpu=%P memory=%M job=%J
endif

# Record run time and memory usage in a file using GNU time
ifeq ($(time),True)
ifneq ($(shell command -v gtime),)
log_time=command gtime -v -o $@.time
else
log_time=command time -v -o $@.time
endif
else
log_time=
endif

# Compress in parallel
ifneq ($(shell command -v pigz),)
gzip=pigz -p$t -f
else
gzip=gzip -f
endif

# Set-up PYTHON PATH
PYTHONPATH=$(assemble_path)/src/btllib/install/lib/btllib/python

help:
	@echo ""
	@echo "ntJoin: Scaffolding assemblies using reference assemblies and minimizer graphs"
	@echo "ntJoin v1.1.0"
	@echo "Usage: ntJoin assemble target=<target scaffolds> references='List of reference assemblies' reference_weights='List of weights per reference assembly'"
	@echo ""
	@echo "Options:"
	@echo "target			Target assembly to be scaffolded in fasta format"
	@echo "references		List of reference files (separated by a space, in fasta format)"
	@echo "target_weight		Weight of target assembly [1]"
	@echo "reference_weights	List of weights of reference assemblies"
	@echo "prefix			Prefix of intermediate output files [out.k<k>.w<w>.n<n>]"
	@echo "t			Number of threads [4]"
	@echo "k			K-mer size for minimizers [32]"
	@echo "w			Window size for minimizers (bp) [1000]"
	@echo "n			Minimum graph edge weight [1]"
	@echo "g			Minimum gap size (bp) [20]"
	@echo "G			Maximum gap size (bp) (0 if no maximum) [0]"
	@echo "m			Minimum percentage of increasing/decreasing minimizer positions to orient contig [90]"
	@echo "mkt			If True, use Mann-Kendall Test to predict contig orientation (computationally-intensive, overrides 'm') [False]"
	@echo "agp			If True, output AGP file describing output scaffolds [False]"
	@echo "no_cut		    	If True, will not cut contigs at putative misassemblies [False]"
	@echo "overlap			If True, attempts to detect and trim overlaps between joined sequences [True]"
	@echo "time		    	If True, will log the time for each step [False]"
	@echo "reference_config	Config file with reference assemblies and reference weights as comma-separated values (See README for example)"
	@echo "			 This is optional, and will override the 'references' and 'reference_weights' variables if specified"
	@echo ""
	@echo "Notes: "
	@echo "	- Ensure the lists of reference assemblies and weights are in the same order, and that both are space-separated"
	@echo "	- Ensure all assembly files are in the current working directory"
	@echo ""

assemble: check_params \
	$(min_ref_targets) \
	$(target).k$(k).w$(w).tsv \
	$(fai_ref_targets) \
	$(target).fai \
	$(target).k$(k).w$(w).n$(n).assigned.scaffolds.fa \
	$(target).k$(k).w$(w).n$(n).unassigned.scaffolds.fa \
	$(target).k$(k).w$(w).n$(n).all.scaffolds.fa

analysis: check_params_analysis \
	$(addsuffix .bam.bai, $(references_list)) \
	$(target).bam.bai \
	$(target).k$(k).w$(w).n$(n).all.scaffolds.fa.bam.bai

all: check_params assemble analysis

check_params:
ifeq ($(references), None)
ifeq ($(reference_config), None)
	$(error ERROR: Must set references)
endif
endif

ifeq ($(reference_weights), None)
ifeq ($(reference_config), None)
	$(error ERROR: Must set reference_weights)
endif
endif

ifeq ($(target), None)
	$(error ERROR: Must set target)
endif

check_params_analysis:
ifeq ($(ref), None)
	$(error must set ref)
endif

version:
	@echo "ntJoin v1.1.0"
	@echo "Written by Lauren Coombe (lcoombe@bcgsc.ca)"

.PHONY: help all version analysis assemble check_params jupiter
.DELETE_ON_ERROR: $(prefix).n$(n).mx.dot
.SECONDARY:

%.k$(k).w$(w).tsv: %
	$(log_time) $(assemble_path)/src/indexlr --long --pos -k $(k) -w $(w) -t $(t) $< > $@

%.fai: %
	$(log_time) samtools faidx $<

%.fa.gz: %.fa
	$(log_time) $(gzip) $<

# Set ntJoin options
ifeq ($(mkt), True)
    ntjoin_assembly_opts += --mkt
endif
ifeq ($(no_cut), True)
    ntjoin_assembly_opts += --no_cut
endif
ifeq ($(agp), True)
    ntjoin_assembly_opts += --agp
endif
ifeq ($(overlap), True)
    ntjoin_assembly_opts += --overlap --overlap_gap $(overlap_g) --btllib_t $(t)
    ntjoin_assembly_opts += --overlap_k $(overlap_k) --overlap_w $(overlap_w)
endif
	
$(target).k$(k).w$(w).n$(n).assigned.scaffolds.fa: $(target).k$(k).w$(w).tsv $(min_ref_targets)
	PYTHONPATH=$(PYTHONPATH) $(log_time)  $(assemble_path)/bin/ntjoin_assemble.py -p $(prefix) -n $(n) -s $< -l $(target_weight) \
	-r "$(reference_weights_list)" -k $(k) -g $(g) -G $(G) -t $(assemble_t) $(ntjoin_assembly_opts) $(min_ref_targets)

$(target).k$(k).w$(w).n$(n).unassigned.scaffolds.fa: $(target).k$(k).w$(w).n$(n).assigned.scaffolds.fa
	touch $@

%.all.scaffolds.fa: %.assigned.scaffolds.fa %.unassigned.scaffolds.fa
	$(log_time) cat $^ > $@

%.fa.bam: %.fa
	minimap2 -a -x asm5 -r100000 -t $(t) $(ref) $< |samtools view -b |samtools sort -o $@

%.bam.bai: %.bam
	samtools index $<

# Run QUAST
quast_$(prefix)/report.tsv: $(references_list) $(target) $(target).k$(k).w$(w).n$(n).all.scaffolds.fa
ifeq ($(large), 1)
	quast -t $(t) -o quast_$(prefix) -r $(ref) --fast --scaffold-gap-max-size 100000 --split-scaffolds \
	--large $^
else
	quast -t $(t) -o quast_$(prefix) -r $(ref) --fast --scaffold-gap-max-size 100000 --split-scaffolds \
	$^
endif