#!/bin/bash

# Example exhaustive run of the MMSEQ pipeline (see http://bgx.org.uk/software/mmseq.html).
# You will need at least
# - a genome fasta
# - a transcript fasta
# - a GFF file for the transcripts in the above fasta
# - a set of paired fastq files

# You might also want
# - pre-built Bowtie index files for the genome fasta

# Note: where indicated, groups of commands may safely be run in parallel to save time

### SET THESE VARIABLES
# Place reference files (fasta, gff, ebwt) in $BOWTIE_INDEXES/
export BOWTIE_INDEXES=refs/ # trailing slash is necessary
GENOME_FASTA=Homo_sapiens.GRCh37.64.dna.toplevel.ref.fa
GENOME_EBWT=`basename $GENOME_FASTA .fa`
TRANSCRIPT_FASTA=Homo_sapiens.GRCh37.64.ref_transcripts.fa
GFF_FILE=Homo_sapiens.GRCh37.64.ref.gff

FQ_BASENAMES=( 3125_7 3125_2 ) # excludes _1.fastq and _2.fastq suffixes
FQDIR=fqs
READ_LENGTHS=( 37 37 )
EXPECTED_ISIZES=( 192 189 )
DEVIATION_THRESHOLDS=( 50 50 ) # e.g. bp from mean to reach 1.sd
MAX_ISIZES=( 400 400 ) # twice expected isizes should be OK
THREADS=16 # CPU cores
CHUNKMBS=128

# directories in which to save output files
TOPHATOUT=thout
PH2TEMP=ph2temp
PH2OUT=ph2out
RESULTS=results

### DIAGNOSTIC CHECKS
# check all necessary binaries and scripts are available
CMDS=( mmseq bam2hits samtools bcftools vcfutils.pl bowtie bowtie-build tophat ruby haploref.rb java )
for c in ${CMDS[@]}; do
  type -P $c &> /dev/null || { echo "Error: can't find $c. Make sure it's in your PATH." >&2; exit 1; }
done

# check that PolyHap2.jar is in CLASSPATH
PH=`java dataFormat.ImportSNPs 2>&1 | grep NoClassDefFoundError`
if [ ${#PH} -gt 0 ]; then
  echo "Error: make sure PolyHap2.jar is in your CLASSPATH."
  exit 1;
fi

# check these three arrays have the same length (= no. of samples)
if [ ! ${#FQ_BASENAMES[@]} -eq ${#EXPECTED_ISIZES[@]} ] ||
  [ ! ${#FQ_BASENAMES[@]} -eq ${#DEVIATION_THRESHOLDS[@]} ] ||
  [ ! ${#FQ_BASENAMES[@]} -eq ${#READ_LENGTHS[@]} ]; then
  echo "Error: FQS_BASENAMES, EXPECTED_ISIZES, DEVIATION_THRESHOLDS have different no. of elements."
  exit 1;
fi

for(( i=0; i < ${#EXPECTED_ISIZES[@]}; i++ )); do
  INNER_DIST[i]=$(( ${EXPECTED_ISIZES[i]} - 2*${READ_LENGTHS[i]} ))
done

# check the files and directories exist
[ -d $BOWTIE_INDEXES ] || { echo "Error: can't find bowtie indexes directory $BOWTIE_INDEXES." >&2; exit 1; } 
[ -s $BOWTIE_INDEXES/$GENOME_FASTA ] || { echo "Error: can't find non-empty genome FASTA $GENOME_FASTA in $BOWTIE_INDEXES." >&2; exit 1; }
[ -s $BOWTIE_INDEXES/$TRANSCRIPT_FASTA ] || { echo "Error: can't find non-empty transcript FASTA $TRANSCRIPT_FASTA in $BOWTIE_INDEXES." >&2; exit 1; }
[ -s $BOWTIE_INDEXES/$GFF_FILE ] || { echo "Error: can't find non-empty GFF file $GFF_FILE in $BOWTIE_INDEXES." >&2; exit 1; } 
for (( i=0; i < ${#FQ_BASENAMES[@]}; i++ )); do
  [ -s $FQDIR/${FQ_BASENAMES[i]}_1.fastq ] || { echo "Error: can't find non-empty FASTQ file $FQDIR/${FQ_BASENAMES[i]}_1.fastq." >&2; exit 1; } 
  [ -s $FQDIR/${FQ_BASENAMES[i]}_2.fastq ] || { echo "Error: can't find non-empty FASTQ file $FQDIR/${FQ_BASENAMES[i]}_2.fastq." >&2; exit 1; }
done

mkdir -p $TOPHATOUT
mkdir -p $RESULTS

### START EXHAUSTIVE PIPELINE
echo ""
echo "%%%%%%%%%%%% STEP B1 %%%%%%%%%%%%"
echo ""
if [ ! -s $BOWTIE_INDEXES/$GENOME_EBWT.1.ebwt ]; then
  cd $BOWTIE_INDEXES
  bowtie-build -f $GENOME_FASTA $GENOME_EBWT
  cd - > /dev/null
else
  echo "Skipping bowtie-build: $BOWTIE_INDEXES/$GENOME_EBWT.1.ebwt exists."
fi

echo ""
echo "%%%%%%%%%%%% STEP B2 %%%%%%%%%%%%"
echo ""
# these iterations may be run in parallel 
for (( i=0; i < ${#FQ_BASENAMES[@]}; i++ )); do
  if [ -s ${TOPHATOUT}/${FQ_BASENAMES[i]}.tophat/accepted_hits.bam ]; then
    echo "Skipping tophat: ${TOPHATOUT}/${FQ_BASENAMES[i]}.tophat/accepted_hits.sam exists."
  else
    tophat -G $BOWTIE_INDEXES/$GFF_FILE --no-novel-juncs --min-isoform-fraction 0.0 --min-anchor-length 3 -r ${INNER_DIST[i]} \
      -p $THREADS -o ${TOPHATOUT}/${FQ_BASENAMES[i]}.tophat $GENOME_EBWT $FQDIR/${FQ_BASENAMES[i]}_1.fastq $FQDIR/${FQ_BASENAMES[i]}_2.fastq
  fi
done
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP B3 %%%%%%%%%%%%"
echo ""
# these iterations may be run in parallel 
for (( i=0; i < ${#FQ_BASENAMES[@]}; i++ )); do
  cd ${TOPHATOUT}/${FQ_BASENAMES[i]}.tophat
  if [ -s accepted_hits.sorted.bam ]; then
    echo "Skipping sorting BAM: ${TOPHATOUT}/${FQ_BASENAMES[i]}.tophat/accepted_hits.sorted.bam exists."
  else
    echo "Sorting ${TOPHATOUT}/${FQ_BASENAMES[i]}.tophat/accepted_hits.bam..."
    samtools sort accepted_hits.bam accepted_hits.sorted
  fi
  cd - > /dev/null
done
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP B4 %%%%%%%%%%%%"
echo ""

if [ -s ${TOPHATOUT}/var.flt.vcf ]; then
  echo "Skipping variant calling: VCF file var.flt.vcf exists in ${TOPHATOUT}/."
else
  echo "Calling variants using sorted BAM files in ${TOPHATOUT}/*.tophat..."
  cd $TOPHATOUT
  SBAMS=( `find . -wholename "*.tophat/accepted_hits.sorted.bam"` )
  samtools mpileup -ugf $BOWTIE_INDEXES/$GENOME_FASTA ${SBAMS[@]} | bcftools view -bvcg - > var.raw.bcf
  bcftools view var.raw.bcf | vcfutils.pl varFilter > var.flt.vcf
  cd - > /dev/null
fi
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP B5 %%%%%%%%%%%%"
echo ""
if [ -s $PH2TEMP/phase.sh ]; then
  echo "Skipping importing SNPs: phase.sh exists in $PH2TEMP/."
else
  java dataFormat.ImportSNPs $BOWTIE_INDEXES/$GFF_FILE $TOPHATOUT/var.flt.vcf $PH2TEMP
fi
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP B6 %%%%%%%%%%%%"
echo ""
OF=`ls $PH2TEMP | grep .out | wc -l`
if [ $OF -gt 0 ]; then
  echo "Skipping phasing: $OF .out files exist in $PH2TEMP/."
else
  # the commands in the phase.sh script may be run in parallel 
  ./$PH2TEMP/phase.sh
fi
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP B7 %%%%%%%%%%%%"
echo ""
if [ -s $PH2OUT/SNPinfo ]; then
  echo "Skipping AB to ATGC conversion: $PH2OUT/SNPinfo exists."
else
  java dataFormat.ConvertABtoATGC $BOWTIE_INDEXES/$GFF_FILE $PH2TEMP $PH2OUT
fi
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP B8 %%%%%%%%%%%%"
echo ""
cd $PH2OUT
# these iterations may be run in parallel 
for (( i=0; i < ${#FQ_BASENAMES[@]}; i++)); do
  if [ -s  $BOWTIE_INDEXES/${FQ_BASENAMES[i]}.fa ]; then
    echo "Skipping creating custom transcript fasta: $BOWTIE_INDEXES/${FQ_BASENAMES[i]}.fa exists."
  else
    haploref.rb $BOWTIE_INDEXES/$TRANSCRIPT_FASTA $BOWTIE_INDEXES/$GFF_FILE SNPinfo ${FQ_BASENAMES[i]}.hap > $BOWTIE_INDEXES/${FQ_BASENAMES[i]}.fa
  fi
done
cd - > /dev/null
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP A1 %%%%%%%%%%%%"
echo ""
# these iterations may be run in parallel 
for (( i=0; i < ${#FQ_BASENAMES[@]}; i++)); do
  if [ ! -s $BOWTIE_INDEXES/${FQ_BASENAMES[i]}.1.ebwt ]; then
    cd $BOWTIE_INDEXES
    bowtie-build -f ${FQ_BASENAMES[i]}.fa ${FQ_BASENAMES[i]}
    cd - > /dev/null
  else
    echo "Skipping creating custom transcript bowtie indexes: $BOWTIE_INDEXES/${FQ_BASENAMES[i]}.*.ebwt exist."
  fi
done
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP A2 %%%%%%%%%%%%"
echo ""
# these iterations may be run in parallel 
for (( i=0; i < ${#FQ_BASENAMES[@]}; i++)); do
  if [ ! -s $RESULTS/${FQ_BASENAMES[i]}.sam ]; then
    bowtie -a --best --strata -S -X ${MAX_ISIZES[i]} -p $THREADS --chunkmbs $CHUNKMBS --quiet ${FQ_BASENAMES[i]} \
      -1 $FQDIR/${FQ_BASENAMES[i]}_1.fastq -2 $FQDIR/${FQ_BASENAMES[i]}_2.fastq $RESULTS/${FQ_BASENAMES[i]}.sam
  else
    echo "Skipping alignment to custom reference: $RESULTS/${FQ_BASENAMES[i]}.sam exists."
  fi
done
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP A3 %%%%%%%%%%%%"
echo ""
# these iterations may be run in parallel 
for (( i=0; i < ${#FQ_BASENAMES[@]}; i++)); do
  if [ ! -s $RESULTS/${FQ_BASENAMES[i]}.sorted.bam ]; then
    samtools view -bS $RESULTS/${FQ_BASENAMES[i]}.sam > $RESULTS/${FQ_BASENAMES[i]}.bam
    samtools sort -n $RESULTS/${FQ_BASENAMES[i]}.bam $RESULTS/${FQ_BASENAMES[i]}.sorted
  else
    echo "Skipping converting SAM to query-name sorted BAM: $RESULTS/${FQ_BASENAMES[i]}.sorted.bam exists."
  fi
done
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP A4 %%%%%%%%%%%%"
echo ""
# these iterations may be run in parallel 
for (( i=0; i < ${#FQ_BASENAMES[@]}; i++)); do
  if [ ! -s $RESULTS/${FQ_BASENAMES[i]}.hits ]; then
    bam2hits -i ${EXPECTED_ISIZES[i]} ${DEVIATION_THRESHOLDS[i]} \
      $BOWTIE_INDEXES/${FQ_BASENAMES[i]}.fa $RESULTS/${FQ_BASENAMES[i]}.sorted.bam > $RESULTS/${FQ_BASENAMES[i]}.hits
  else
    echo "Skipping creating hits file: $RESULTS/${FQ_BASENAMES[i]}.hits exists."
  fi
done
[ $? -eq 0 ] || { exit; }

echo ""
echo "%%%%%%%%%%%% STEP A5 %%%%%%%%%%%%"
echo ""
# these iterations may be run in parallel 
for (( i=0; i < ${#FQ_BASENAMES[@]}; i++)); do
  if [ ! -s $RESULTS/${FQ_BASENAMES[i]}.mmseq ]; then
    mmseq $RESULTS/${FQ_BASENAMES[i]}.hits $RESULTS/${FQ_BASENAMES[i]}
  else
    echo "Skipping running mmseq: $RESULTS/${FQ_BASENAMES[i]}.mmseq exists."
  fi
done


