# Accession number for the genome.
GCF = GCF_000848505.1

# The name of the genome
GENOME_NAME = ebola-1976

# The referencefasta file for the genome
GENOME_FASTA = fasta/${GENOME_NAME}.fa

# The reference GFF file.
GENOME_GFF = fasta/${GENOME_NAME}.gff

# The source genome in not identical to the reference.
SOURCE_GENOME = fasta/${GENOME_NAME}.source.fa

# Sequencing reads.
R1 = fastq/read_1.fq
R2 = fastq/read_2.fq

BAM = bam/${GENOME_NAME}.bam

# Simulation parameters

# Error rate for the reads.
ERROR_RATE = 0

# Mutation rate for the reads.
MUTATION_RATE = 0

# Indel rate for the reads.
INDEL_RATE = 0

# Number of reads.
N_READS = 2000

# DNA fragment size.
FRAG_SIZE = 500

# Standard deviation for the reads.
STD_DEV = 10

# Length of the reads.
LENGTH = 100

# Seed for the random number generator.
SEED = 0

# 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 --

help:
	@echo "# Read the source code"

# Download the genome
${GENOME_FASTA} ${GENOME_GFF} ${SOURCE_GENOME} :
	# Download the genome
	datasets download genome accession ${GCF} --include genome,gff3,gtf

	# Create the unpack directory
	mkdir -p ncbi 

	# Unpack the data into the ncbi directory
	unzip -o -j -d ncbi ncbi_dataset.zip

	# Create the reference directory
	mkdir -p $(dir ${GENOME_FASTA})

	# Copy the genome to a simpler name
	cp -f ncbi/*.fna ${GENOME_FASTA}

	# Initially assume that thereal genome is the same as the reference genome
	cp -f ${GENOME_FASTA} ${SOURCE_GENOME}

	# Copy the gff file
	cp -f ncbi/*.gff ${GENOME_GFF}

	# Print the statistics of the genome
	seqkit stats ${GENOME_FASTA}

	# Compute bwa index for the reference genome
	bwa index ${GENOME_FASTA}

genome: ${GENOME_FASTA} ${GENOME_GFF} ${SOURCE_GENOME}
	ls -l ${GENOME_FASTA} ${GENOME_GFF} ${SOURCE_GENOME}

align: ${BAM}
	ls -l ${BAM}

${BAM}: $(GENOME_FASTA) ${R1} ${R2}
	mkdir -p $(dir ${BAM})
	bwa mem ${GENOME_FASTA} ${R1} ${R2} | samtools sort --write-index -o $@

# Simulate the reads from the source genome.
fastq: $(SOURCE_GENOME) 
	mkdir -p $(dir ${R1})
	wgsim -1 ${LENGTH} -2 ${LENGTH} -d ${FRAG_SIZE} -e ${ERROR_RATE} -r ${MUTATION_RATE} -R ${INDEL_RATE} \
		-s ${STD_DEV} -N ${N_READS} -S ${SEED} ${SOURCE_GENOME} ${R1} ${R2} > mutations.txt
	seqkit stats ${R1} ${R2}

all: fastq align
	ls -l ${BAM}

clean:
	rm -rf ${R1} ${R2} ${BAM}

# Visualize with IGV 
igv: fastq align ${BAM}
	echo "new" | ${NC}
	echo "genome ${DIR}/$(GENOME_FASTA)" | ${NC}
	echo "load ${DIR}/$(BAM)" | ${NC}

.PHONY: all clean igv fastq align genome