longReadValidate

#!/bin/bash

############################################################
#  Program:
#  Author :
############################################################


## BEGIN SCRIPT
usage()
{
    cat << EOF

usage: $0 OPTIONS

OPTIONS:
    -h      Show this message
    -i      BEDPE input file name
    -s      slop (default: 0)
    -m      Moleculo file
    -p      PacBio file
    -k      keep temporary files

EOF
}

# Show usage when there are no arguments.
if test -z "$1"
then
    usage
    exit
fi

BIN_DIR=$( dirname $0 )
PB=
MO=
INPUT=
SLOP=0
KEEP=0

# Check options passed in.
while getopts "h m:p:i:s:k" OPTION
do
    case $OPTION in
        h)
            usage
            exit 1
            ;;
        m)
            MO=$OPTARG
            ;;
        p)
            PB=$OPTARG
            ;;
        i)
            INPUT=$OPTARG
            ;;
	s)
	    SLOP=$OPTARG
	    ;;
	k)
	    KEEP=1
	    ;;
        ?)
            usage
            exit
            ;;
    esac
done

if [[ ! -f $PB ]] || [[ ! -f $MO ]]
then
    usage
    echo -e "Error: -p and -m are required\n"
    exit 1
fi

# calculate the number of columns in the input file
NCOL=`head -n 1 $INPUT | awk '{ print NF }'`

pairToPair -type both -is -slop $SLOP -a $PB -b $INPUT \
    | awk '! array[$0]++' \
    > $INPUT.p.slop$SLOP.tmp

cat $INPUT.p.slop$SLOP.tmp | sort -k 25 > $INPUT.p.slop$SLOP.tmp.sort

cat $INPUT.p.slop$SLOP.tmp.sort \
    | $BIN_DIR/filter_pb.py \
    > $INPUT.p.slop$SLOP.tmp

pairToPair -type both -is -slop $SLOP -a $MO -b $INPUT \
    | awk '! array[$0]++' \
    > $INPUT.m.slop$SLOP.tmp

cat $INPUT.p.slop$SLOP.tmp \
    | awk '$18==$29' \
    | cut -f 19- \
    | sort -k7,7n \
    | groupBy -g 1,2,3,4,5,6 -c 1 -o count -full \
    | $BIN_DIR/zjoin_empty -r -a $INPUT -b stdin -1 7 -2 7 \
    | cut -f -$NCOL,$(($NCOL+$NCOL+1)) \
    | awk '{ if ($NF=="NA") { $NF=0 } print $0 }' OFS="\t" \
    | $BIN_DIR/zjoin_empty -r -1 7 -2 7 -a stdin -b <(cat $INPUT.m.slop$SLOP.tmp | awk '$18==$29' | cut -f 19- | sort -k7,7n | groupBy -g 1,2,3,4,5,6 -c 1 -o count -full) \
    | cut -f -$(($NCOL+1)),$(($NCOL+$NCOL+2)) \
    | awk '{ if ($NF=="NA") { $NF=0 } print $0 }' OFS="\t"

    # > $INPUT.slop$SLOP.val

if [[ "$KEEP" -eq 0 ]]
then
    rm $INPUT.p.slop$SLOP.tmp $INPUT.m.slop$SLOP.tmp $INPUT.p.slop$SLOP.tmp.sort
fi