#
# Variant calling
#
# Uses the NCBI datasets tool to download genome and annotation files.
#
# Uses bwa to align and then bcftools to call the variants.
#

#
# The NCBI accession number
#
ACC = GCF_000848505.1

# The name used to name the files.
NAME = ebola-mayinga-1976

# The zip archive containing the NCBI genome and annotation files.
NCBI_ZIP = ncbi/$(ACC).zip

# Genome fasta file.
FASTA = fasta/$(NAME).fasta

# GFF annotation file.
GFF = gff/$(NAME).gff

# GTF annotation file.
GTF = gtf/$(NAME).gtf

# SRR accession number
SRR=SRR1553425

# The sample name for the SRR file. 
SAMPLE = EM100

# Number of reads to download
N = 20000

# If it is a paired end sample, the read names.
R1 = fastq/$(SRR)_1.fastq
R2 = fastq/$(SRR)_2.fastq

# BAM file derived from the sample name
BAM = bam/$(SAMPLE).bam

# BWA index sentinel file.
INDEX = $(FASTA).bwt

# The output VCF file.
VCF=vcf/$(SAMPLE).vcf.gz

# Annotation parameters for the mpileup command
ANN=-d 100 --annotate 'INFO/AD,FORMAT/DP,FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/SP'

# Call parameters for the call command
CALL=--ploidy 2 --annotate 'FORMAT/GQ' 

# Current directory
DIR = $(shell pwd)

# Netcat command to send commands to IGV
NC = nc localhost 60151

# --- All definitions above this line. All code below this line ---

# Print the usage information.
usage:
	@echo "#"
	@echo "# --- Variant calling inputs ---"
	@echo "#"
	@echo "# ACC=$(ACC)"
	@echo "# NAME=$(NAME)"
	@echo "#"
	@echo "# SRR=$(SRR)"
	@echo "# SAMPLE=$(SAMPLE)"
	@echo "# N=$(N)"
	@echo "#"
	@echo "# --- Generated file names ---"
	@echo "#"
	@echo "# FASTA=$(FASTA)"
	@echo "# GFF=$(GFF)"
	@echo "#"
	@echo "# R1=$(R1)"
	@echo "# R2=$(R2)"
	@echo "#"
	@echo "# --- Output files ---"
	@echo "#"
	@echo "# BAM=$(BAM)"
	@echo "# BAM.bai=$(BAM).bai"
	@echo "# VCF=$(VCF)"
	@echo "#"
	@echo "# make all"
	@echo "#"

# The directory containing the NCBI zip file.
NCBI_DIR = $(dir $(NCBI_ZIP))

# Download the NCBI zip file.
$(NCBI_ZIP):
	mkdir -p $(NCBI_DIR)
	datasets download genome accession $(ACC) --filename $(NCBI_ZIP) --include genome,gff3,gtf
	unzip -n -d $(NCBI_DIR) $(NCBI_ZIP)

# Get the NCBI zip file
zip: $(NCBI_ZIP)
	ls -lh $(NCBI_ZIP)

# Extract the genome and annotation files from the NCBI zip file.
$(FASTA) $(GFF) $(GTF): $(NCBI_ZIP)
	# Create the directory for the FASTA file.
	mkdir -p $(dir $(FASTA))
	cp $(NCBI_DIR)/ncbi_dataset/data/GCF_000848505.1/*.fna $(FASTA)
	
	# Create the directory for the GFF file.
	mkdir -p $(dir $(GFF))
	cp $(NCBI_DIR)/ncbi_dataset/data/GCF_000848505.1/genomic.gff $(GFF)

	# Create the directory for the GTF file.
	mkdir -p $(dir $(GTF))
	cp $(NCBI_DIR)/ncbi_dataset/data/GCF_000848505.1/genomic.gtf $(GTF)

genome: $(FASTA) $(GFF) $(GTF)
	ls -lh $(FASTA)
	ls -lh $(GFF)
	ls -lh $(GTF)

# Download the fastq files.
fastq: $(R1) $(R2)

$(R1) $(R2):
	mkdir -p $(dir $@)
	fastq-dump -X $(N) --split-files -O $(dir $@) $(SRR)

# Index the genome.
index: $(INDEX)

$(INDEX): $(FASTA)
	bwa index $<

# Align the fastq files to the genome.
align: $(BAM)
	ls -lh $(BAM)

$(BAM): $(INDEX) $(R1) $(R2)
	mkdir -p $(dir $@)
	bwa mem -t 4 $(FASTA) $(R1) $(R2) | samtools sort -o $@ 
	samtools index $@
	samtools flagstat $@

# Call the variants.
call: $(VCF)
	ls -lh $(VCF)

# Create the VCF file
$(VCF): $(BAM) $(FASTA)
	# Create the directory for the VCF file.
	mkdir -p $(dir $@)

	# Run the variant calling process.
	bcftools mpileup ${ANN} -O u -f $(FASTA) $(BAM) | \
		bcftools call $(CALL) -mv -O u | \
		bcftools norm -f $(FASTA) -d all -O u | \
		bcftools sort --write-index -O z -o $(VCF)

# Run all steps.
all: genome call 

# Visualize with IGV 
igv: fastq align ${BAM} $(GFF)
	echo "new" | ${NC}
	echo "genome ${DIR}/$(FASTA)" | ${NC}
	echo "load ${DIR}/$(BAM)" | ${NC}
	echo "load ${DIR}/$(VCF)" | ${NC}
	echo "load ${DIR}/$(GFF)" | ${NC}
	
# Remove some of the generated files.
clean:
	rm -rf $(R1) $(R2) $(BAM) $(BAM).bai $(VCF) $(VCF).tbi

# Clean up all subdirectories.
realclean:
	rm -rf fasta gff gtf ncbi bam vcf fastq

.PHONY: usage all genome fastq index align clean call igv