setup

#!/bin/bash
trap 'echo $0: line $LINENO: exit code $?' ERR
set -e

add_dependency() {
    label=$1
    check_command=$2
    extra_path=

    while ! PATH=$new_path:$extra_path:$PATH eval $check_command; do
	echo "$label not found"
	extra_path=$(readlink -f "$(read -p "Enter $label path: " -e extra_path; eval echo "$extra_path")")
    done
    [ ! "$extra_path" ] || new_path=$new_path${new_path:+:}$extra_path
}


DIR=$( cd -P "$( dirname "${BASH_SOURCE[0]}" )"; pwd )
REF_NAME=${1:-}

[ -d "$DIR/settings" ] || mkdir settings
[ -d "$DIR/data" ] || mkdir data

GLOBAL_SETTINGS=$DIR/settings/global
if [ -r "$GLOBAL_SETTINGS" ]; then
    echo "Global settings file already exists: $GLOBAL_SETTINGS"
    echo "If you need to regenerate it, remove the existing file"
else
    #
    # first the paths
    #
    new_path=
    add_dependency "python2.6+" "which python >/dev/null 2>&1 && [ \"\$(python --version 2>&1 | cut -d '.' -f 2)\" -ge 6 ]"
    add_dependency "samtools" "which samtools >/dev/null 2>&1"
    add_dependency "bedtools" "which intersectBed >/dev/null 2>&1"
    add_dependency "RepeatMasker" "which RepeatMasker >/dev/null 2>&1"
    add_dependency "bowtie2" "which bowtie2 >/dev/null 2>&1"
    add_dependency "gmapper" "which gmapper-ls >/dev/null 2>&1"

    echo "[[ \"\$PATH\" =~ $DIR/bin ]] || PATH=$DIR/bin:$new_path${new_path:+:}\$PATH" >"$GLOBAL_SETTINGS"
    echo "[[ \"\$AWKPATH\" =~ $DIR/bin ]] || export AWKPATH=$DIR/bin\${AWKPATH:+:}\$AWKPATH" >>"$GLOBAL_SETTINGS"
#    echo "[[ \"\$PYTHONPATH\" =~ $DIR/bin ]] || PYTHONPATH=$DIR/bin:\$PYTHONPATH" >>"$GLOBAL_SETTINGS"

    #
    # next, alu consensus sequences
    #
    if [[ -z "$ALU_ALL_FA" || ! -r "$ALU_ALL_FA" ]]; then
	echo "Need to build list of Alu consensus sequences"
	ALU_ALL_FA=$DIR/data/alus.pos.fa
	ALU_NO_POLYA_FA=$DIR/data/alus.hidden-polya.fa

        REPEATMASKER_DIR=$(dirname "$(readlink -f "$(PATH=$new_path:$PATH which RepeatMasker)")")

	echo -n "Generating list of alu consensus sequences... "
	$REPEATMASKER_DIR/util/queryRepeatDatabase.pl -species homo -class SINE | 
	$DIR/bin/grep-fasta "^>Alu[JSY][a-z0-9]*#" | 
	sed 's/^\(.*\)#SINE\/Alu/\1/' | 
	sed 's/^\(.*\)RepbaseID:/\1/' | 
	sed 's/^\(.*\)ALU/\1/' | 
	awk '/^>/ {if (substr($2,1,3)=="Alu"&&$2!=substr($1,2)) $1=$1"/"substr($2,4); print $1} /^[^>]/' >$ALU_ALL_FA

	if [ -z "$(head -n 1 $ALU_ALL_FA)" ]; then
	    echo "oops"
	    cat <<EOF
Your RepeatMasker installation does not include the Repbase libraries,
which are needed to derive the Alu consensus sequences. The libraries
can be downloaded from http://www.girinst.org/server/RepBase/index.php
EOF
	    rm -f "$GLOBAL_SETTINGS"
	    exit 1
	else
	    echo "done"
	fi

	echo "Building list with hidden poly-A streches"
	$DIR/bin/hide-polya <$ALU_ALL_FA >$ALU_NO_POLYA_FA
    fi
    echo "export ALU_ALL_FA=\"$ALU_ALL_FA\"" >>"$GLOBAL_SETTINGS"
    echo "export ALU_NO_POLYA_FA=\"$ALU_NO_POLYA_FA\"" >>"$GLOBAL_SETTINGS"
fi


while [ -z "$REF_NAME" ]; do
    read -p "Enter reference name: " -e REF_NAME
done
REF_SETTINGS=$DIR/settings/ref.$REF_NAME
if [ -r "$REF_SETTINGS" ]; then
    echo "Reference settings file already exists: $REF_SETTINGS"
    echo "If you need to regenerate it, remove the existing file"
else
    . "$GLOBAL_SETTINGS"
    #
    # fasta file
    #
    while [[ -z "$REF_FA" || ! -r "$REF_FA" ]]; do
	read -p "Enter reference genome FASTA file: " -e REF_FA
	if [[ "$REF_FA" && ! -r "$REF_FA" ]]; then
	    echo "$REF_FA: not found"
	    REF_FA=
	fi
    done
    echo "export REF_FA=\"$REF_FA\"" >"$REF_SETTINGS"

    # link if not already there
    REF_FA_LINK=$DIR/data/$REF_NAME.fa
    [ -f "$REF_FA_LINK" -a ! -L "$REF_FA_LINK" ] || ln -sf "$REF_FA" "$REF_FA_LINK"

    #
    # fasta index
    #
    CHROMINFO=$REF_FA.fai
    CHROMINFO_LINK="$REF_FA_LINK".fai
    if [ -r "$CHROMINFO" ]; then
	[ -f "$CHROMINFO_LINK" -a ! -L "$CHROMINFO_LINK" ] || ln -sf "$CHROMINFO" "$CHROMINFO_LINK"
    else
	CHROMINFO=$CHROMINFO_LINK
	echo -n "Creating fasta index... "
	samtools faidx "$REF_FA_LINK"
	echo "done"
    fi
    echo "export CHROMINFO=\"$CHROMINFO\"" >>"$REF_SETTINGS"

    #
    # bowtie2 indexes
    #
    while [[ -z "$BOWTIE2_INDEX" || ! -r "$BOWTIE2_INDEX.1.bt2" ]]; do
	read -p "Enter bowtie2 index prefix: " -e BOWTIE2_INDEX
	if [[ "$BOWTIE2_INDEX" && ! -r "$BOWTIE2_INDEX.1.bt2" ]]; then
	    echo "$BOWTIE2_INDEX.1.bt2: not found"
	    BOWTIE2_INDEX=
	fi
    done
    missing=$(diff <(cut -f 1 $CHROMINFO | sort) <(bowtie2-inspect -n $BOWTIE2_INDEX | sort) | grep '^>' | cut -c 3-)
    [ ! "$missing" ] || { echo "Some chromosomes in the bowtie2 index [$BOWTIE2_INDEX] are missing from the fasta file [$REF_FA]:"; echo "$missing"; exit 1; }

    echo "export BOWTIE2_INDEX=\"$BOWTIE2_INDEX\"" >>"$REF_SETTINGS"
    BOWTIE2_INDEX_LINK=$DIR/data/$REF_NAME
    for suffix in .1.bt2 .2.bt2 .3.bt2 .4.bt2 .rev.1.bt2 .rev.2.bt2; do
	[ -f "$BOWTIE2_INDEX_LINK$suffix" -a ! -L "$BOWTIE2_INDEX_LINK$suffix" ] || ln -sf "$BOWTIE2_INDEX$suffix" "$BOWTIE2_INDEX_LINK$suffix"
    done

    #
    # reference alu annotations
    #
    while [[ -z "$ALUS_BED" || ! -r "$ALUS_BED" ]]; do
	ALUS_BED=$DIR/data/alus.$REF_NAME.bed.gz
	echo "Reference Alu annotations not found"
	while [[ -z "$CHROM_FA_OUT" ]]; do
	    echo "RepeatMasker output files not found"
	    read -p "Enter path to RepeatMasker output file(s) (can use wildcards): " -e CHROM_FA_OUT
	    CHROM_FA_OUT=$(eval echo "$CHROM_FA_OUT" | xargs -n 1 readlink -e)
	done
	echo -n "Generating list of reference alus in bed format... "
	grep -h SINE/Alu $CHROM_FA_OUT | 
	sed 's/ \+/\t/g' |
	sed 's/^\t//' |
	sed 's/\t$//' |
	awk 'BEGIN{OFS="\t"} {pct=int($2+$3+$4);if(pct>100)pct=100; if($9=="C"){strand="-"}else{strand="+"}if(substr($12,1,1)=="("){start=$13;end=$14}else{start=$12;end=$13} print $5, $6-1, $7, $10, 100-pct, strand, start, end}' | 
	gzip >$ALUS_BED
	echo "done"
    done
    echo "export ALUS_BED=\"$ALUS_BED\"" >>"$REF_SETTINGS"
fi