#!/usr/bin/env bash

#
# This recipe finds low similarity regions in nCov, SARS and bat SARS
#

#
# Output files are in the bed directory.
#

set -uex

# Novel corona virus genome.
nCov=refs/NC_045512.fa

# SARS virus genome
SARS=refs/KT444582.fa

# Bat SARS virus genome
batSARS=refs/MG772933.fa

# Make a folder for the BED files.
mkdir -p bed

# Make a folder for the BAM files
mkdir -p bam

# Which regions of nCov are covered by SARS.
blastn -query $nCov -subject $SARS -outfmt "6 qacc sacc qcovhsp qcovs qstart qend sstart send" | sort -k5n > nCov-SARS-blast.txt

# Which regions of nCov are covered by batSARS.
blastn -query $nCov -subject $batSARS -outfmt "6 qacc sacc qcovhsp qcovs qstart qend sstart send" | sort -k5n > nCov-batSARS-blast.txt

# Covered regions relative to SARS.
cat nCov-SARS-blast.txt| cut -f 1,5,6 > bed/nCov-SARS-covered.bed

# Covered regions relative to batSARS.
cat nCov-batSARS-blast.txt| cut -f 1,5,6 > bed/nCov-batSARS-covered.bed

# Low similiarity regions relative to SARS.
bedtools complement -i bed/nCov-SARS-covered.bed -g ${nCov}.fai > bed/nCov-SARS-missing.bed

# Low similiarity regions relative to batSARS.
bedtools complement -i bed/nCov-batSARS-covered.bed -g ${nCov}.fai > bed/nCov-batSARS-missing.bed

# Align SARS vs nCov.
minimap2 -a $nCov $SARS | samtools sort > bam/nCov-SARS.bam;

# Align batSARS vs nCov.
minimap2 -a $nCov $batSARS | samtools sort > bam/nCov-batSARS.bam;

# Index all resulting BAM files.
find bam -name '*.bam' | parallel samtools index {}