#!/bin/bash
DEST_FILE=""
PAIRED=""
STRATA=""
UNALIGNED=""
HALFPAIRED=""
OUTPUT=10
MAX_ALIGNMENTS=0
TMP_DIR="/tmp/$$"
SUFFIX=':[^:]*$'
PROGRAM_NAME=$0
function usage {
  echo "Usage:" $0 --dest-file filename input_file1 input_file2 ...
  echo '
This script merges read hits generated by running gmapper across several
databases. It works for both sam and shrimp output format.

Parameters:

<out-file1> <out-file2> ...

--dest-file <dest-file>
        This parameter is required. It specifies the file where all hits will be
        stored.

--paired
--unpaired
        These parameters specify whether the hits were generated by running
        gmapper in paired mode (with -p) or not.

--half-paired
	This can be used to report half-paired alignments in paired mode.

--tmp-dir <tmp-dir>
        Specifies a temporary directory to use. Defaults to /tmp/<PID>.

--suffix <suffix>
        Specifies an awk-style regexp matching the ending of a read tag that
        specifies whether the read is the first or second in a pair. This suffix
        is only used when neither --paired nor --unpaired is given, and only if
        shrimp output format is detected (and not sam). Defaults to ":[^:]*$".

--report <amount>
	Output maximum amount alignments per read. If more alignments then this
	reported, then report only the top hits.

--max-alignments <amount>
	If there are more then amount hits per read, suppress all output for this
	read.

--strata
	If strata is enabled only output alignments with highest score for a read

--unaligned
	Output aligned entries

Note:
        If neither --paired nor --unpaired is given, a script is run to detect
        the correct setting.
'
}

if [ -z "$SHRIMP_FOLDER" ]; then
  echo 1>&2 "warning: SHRIMP_FOLDER undefined; assuming ."
  SHRIMP_FOLDER="."
fi

if [ ! -d "$SHRIMP_FOLDER/utils/" ]; then
  echo 1>&2 "error: cannot find directory $SHRIMP_FOLDER/utils/"
  exit 1
fi

until [ -z "$1" ]; do
  case $1 in
    "--dest-file")
      shift
      DEST_FILE=$1; shift;;
    "--paired")
      shift
      PAIRED="--paired";;
    "--unaligned")
      shift
      UNALIGNED="--sam-unaligned";;
    "--half-paired")
      shift
      HALF_PAIRED="--half-paired";;
    "--unpaired")
      shift
      PAIRED="--unpaired";;
    "--report")
      shift
      OUTPUT=$1; shift;;
    "--max-alignments")
      shift
      MAX_ALIGNMENTS=$1; shift;;
    "--strata")
      shift
      STRATA="--strata";;
    "--suffix")
      shift
      SUFFIX=$1; shift;;
    "--tmp-dir")
      shift
      TMP_DIR=$1; shift;;
    *) break;;
  esac
done
if [ -z "$DEST_FILE" ]; then
  echo "ERROR: --dest-file was not specified"
  usage
  exit
fi
if [ -z "$PAIRED" ]; then
  echo "WARNING: --paired or --unpaired was not specified, trying to auto-detect"
fi


# make temporary directory
if [ -d $TMP_DIR ]; then
  echo 1>&2 "error: temporary directory [$TMP_DIR] exists"
  exit 1
fi
mkdir -p $TMP_DIR

# detect output format
if [ ! -r "$1" ]; then
  echo 1>&2 "error: file [$1] unreadable"
  exit 1
fi

IGNORE=`head -1 $1 | cut -c 1`
case "$IGNORE" in
"#" )
  FORMAT="shrimp"
  ;;
"@" )
  FORMAT="sam"
  ;;
* )
  echo 1>&2 "Unrecognized format in file [$1]"
  exit 1
  ;;
esac
echo 1>&2 "Detected format: $FORMAT"
#warning dont quit
#if [ $FORMAT = "sam" ]; then
#	grep -m 1 "@PG.*ID:gmapper.*CL" $1 | grep " --sam-unaligned"
#	if [ $? -eq 0 ]; then 
#		echo "Warning: SAM output was created using --sam-unaligned, merging may include excess entries!"
#	fi
#fi

# if pairing is unspecified, detect it with script
if [ -z "$PAIRED" ]; then
  if head -n 10000 $1 | awk -v SUFFIX=${SUFFIX} -f ${SHRIMP_FOLDER}/utils/test-unpaired-${FORMAT}.awk ; then
    PAIRED="--unpaired"
  else
    PAIRED="--paired"
  fi
  echo 1>&2 "Detected pairing: $PAIRED"
fi

# write headers
$SHRIMP_FOLDER/utils/mappings-merge-headers-same-qr-diff-db "$@" 2>/dev/null \
  > $DEST_FILE

# extract individual bodies and sort them
while [ "$#" -gt 0 ]; do
  if [ ! -r "$1" ]; then
    echo 1>&2 "error: file [$1] unreadable"
    exit 1
  fi

  if [ "$PAIRED" = "--paired" ]; then
    $SHRIMP_FOLDER/utils/mappings-extract-body $1 2>/dev/null \
      | $SHRIMP_FOLDER/utils/mappings-collapse-pairs - 2>/dev/null \
      | sort -T ${TMP_DIR} \
      > $TMP_DIR/$(basename $1).sort
  else
    $SHRIMP_FOLDER/utils/mappings-extract-body $1 2>/dev/null \
      | sort -T ${TMP_DIR} \
      > $TMP_DIR/$(basename $1).sort
  fi
  shift
done

# merge into destination
if [ "$PAIRED" = "--paired" ]; then
  sort -T ${TMP_DIR} -m $TMP_DIR/*.sort \
    | $SHRIMP_FOLDER/utils/mappings-expand-pairs - 2>/dev/null \
    >> ${TMP_DIR}/merged.tmp
else
  sort -T ${TMP_DIR} -m $TMP_DIR/*.sort \
    >> ${TMP_DIR}/merged.tmp
fi

#grab strata or -o options
echo $SHRIMP_FOLDER/utils/mappings-top-hits.py $STRATA $UNALIGNED $PAIRED $HALF_PAIRED --max-alignments $MAX_ALIGNMENTS --report $OUTPUT ${TMP_DIR}/merged.tmp
$SHRIMP_FOLDER/utils/mappings-top-hits.py $STRATA $UNALIGNED $PAIRED $HALF_PAIRED --max-alignments $MAX_ALIGNMENTS --report $OUTPUT ${TMP_DIR}/merged.tmp > $DEST_FILE

# clean up
rm -rf $TMP_DIR
