efg-memsAligner

#!/bin/bash
set -e
set -o pipefail
thisfolder=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd ) # https://stackoverflow.com/questions/59895/how-do-i-get-the-directory-where-a-bash-script-is-located-from-within-the-script

# executables' absolute paths/commands (make sure they work!)
graphaligner=$thisfolder/tools/GraphAligner/bin/GraphAligner
efgmems=$thisfolder/tools/efg-mems/efg-mems
efggafsplitter=$thisfolder/tools/efg-gaf-splitter/efg-gaf-splitter
seqtk=seqtk

# default params
workingfolder="."
extendoptions="--max-cluster-extend 5 -b 10" # GraphAligner default extend options

print_help()
{
   echo "Pipeline to align long reads to indexable Elastic Founder Graphs, based on MEMs"
   echo "Currently, only 1 compute thread is supported, and all intermediate seeds"
   echo " are saved in the working folder"
   echo "	-h --help:  show this screen"
   echo "	-g graph.gfa: semi-repeat-free EFG in xGFA format"
   echo "	-f reads.fastq: reads in FASTQ format"
   echo "	-a alignmentsout.gaf: output alignments in GAF format"
   echo "	-w path: working folder for output and temporary files"
}

# https://stackoverflow.com/questions/12022592/how-can-i-use-long-options-with-the-bash-getopts-builtin
for arg in "$@"; do
	shift
	case "$arg" in
		'--help')                set -- "$@" '-h'   ;;
		*)                       set -- "$@" "$arg" ;;
	esac
done

while getopts "hg:f:a:w:" option; do
	case $option in
		h) # display help
			print_help
			exit;;
		g) # graph
			argg=true
			graph="$OPTARG" ;;
		f) # fastq reads
			argf=true
			reads="$OPTARG" ;;
		a) # output
			arga=true
			alignmentsout="$OPTARG" ;;
		w) # working folder
			argw=true
			workingfolder="$OPTARG" ;;
		\?) # invalid option
			echo "Error: Invalid option"
			print_help
			exit;;
	esac
done

if [[ "$argg" != true ]] || [[ "$argf" != true ]] ; then
	print_help
	exit
fi

# move to working folder
if [[ "$argw" = true ]] ; then
	workingfolder="${workingfolder%/}"
else
	workingfolder="."
fi

# find efg-mems seeds
$efgmems -a "ACGTXN#0" -k 20 --indexing --bdbwt \
	-o "$workingfolder/efgmems_seeds_pre_split.gaf" \
	<(cat <(awk 'NR % 4 == 1 || NR % 4 == 2' $reads | sed 's/^@/>/g' | cut -d' ' -f1) <(awk 'NR % 4 == 1 || NR % 4 == 2' $reads | sed 's/^@/>/g' | cut -d' ' -f1 | sed 's/^>/>rev_/g' | $seqtk seq -r) | tr "N" "X") \
	$graph

# split seeds and reverse if needed
$efggafsplitter \
	$graph \
	"$workingfolder/efgmems_seeds_pre_split.gaf" \
	> "$workingfolder/efgmems_seeds.gaf"

# GraphAligner extend
$graphaligner $extendoptions  \
	-t 1 \
	-g $graph \
	-f $reads \
	--realign "$workingfolder/efgmems_seeds.gaf" \
	-a $alignmentsout