hervx

#!/usr/bin/env bash
set -euo pipefail


function usage() { cat << EOF
hervx: Main entry point to the HERVx WDL and Cromwell pipeline.

USAGE:
  hervx <MODE> [OPTIONS] -r1 <SRR123_1.fastq> -r2 <SRR1234_2.fastq> -o <outdir>

SYNOPSIS:
  HERVx calculates Human Endogenous Retrovirus (HERV) expression in paired-end
RNA-seq data. Characterization of HERV retrotranscriptome is difficult due to
uncertainty that arises in fragment assignment because of sequence similarity.

  The HERVx pipeline runs cutadapt to remove adapter sequences and to perform
quality-trimming, bowtie2 to align reads against the Human reference genome (hg38),
SAMtools to convert from SAM to BAM format  and to sort reads by name, and Telescope
to characterize Human Endogenous Retrovirus (HERV) expression.

  Telescope directly addresses uncertainty in fragment assignment by reassigning
ambiguously mapped fragments to the most probable source transcript as
determined within a Bayesian statistical model.

Required Positional Argument:
  [1] MODE  [Type: Str] Define the cromwell executor mode [Default: local].
                        Vaild mode options include: <local|slurm>
                         a) local: uses local singularity cromwell backend.
                             The local EXECUTOR will run serially on compute
                             instance. This is useful for testing, debugging,
                             or when a users does not have access to a high
                             performance computing environment.
                         b) slurm: uses slurm and singularity cromwell backend.
                             The slurm EXECUTOR will submit jobs to the cluster.
                             It is recommended running hervx in this mode.
Required Arguments:
  -r1, --read-1  [Type: File] Input R1 FastQ file.
  -r2, --read-2  [Type: File] Input R2 FastQ file.
  -o,  --outdir  [Type: Path] Path to output directory.

OPTIONS:
  -b, --basename [Type: Str] Basename of the sample. This is the name of the
                             sample without the PATH and the file extension.
                              WHERE:
                                file extension = MateInfo + .fastq + .gz [opt]
                              Given: /tmp/S25_WT_1.fastq /tmp/S25_WT_2.fastq
                                The base name would be "S25_rep1". This string is
                                used to deterministically resolve output filenames.
                                If not provided, the base name will be resolved
                                automatically by cleaning the provided "-r1" input
                                filename against a list of common extensions.
  -t, --threads  [Type: Int]  Number of threads (Default = 2).
  -g, --gtf      [Type: File] Input annotation file in GTF format.
                              NOTE: User must choose from one of the following
                                GTFs in the container ("/opt2/refs/"):
                                  a) HERV_rmsk.hg38.v2.transcripts.gtf (Default)
                                  b) HERV_rmsk.hg38.v2.genes.gtf
                                  c) L1Base.hg38.v1.transcripts.gtf
                                  d) retro.hg38.v1.transcripts.gtf
  -p, --prior    [Type: Int]   Prior on theta for Telescope (Default = 200000).
                                 Equivalent to adding N non-unique reads.
                                 NOTE: It is recommended to set this prior to a
                                 large value. This increases the penalty for
                                 non-unique reads and improves accuracy.
  -m, --max-iter [Type: Int]   Maximum number of iterations to test convergence
                                 of the EM algorithm (Default = 200).
  -h, --help     [Type: Bool]  Displays usage and help information.

Example:
  $ hervx local -r1 $(dirname  "$0")/tests/small_S25.R1.fastq.gz -r2 $(dirname  "$0")/tests/small_S25.R2.fastq.gz --outdir /data/$USER/
  $ hervx slurm -r1 $(dirname  "$0")/tests/small_S25.R1.fastq.gz -r2 $(dirname  "$0")/tests/small_S25.R2.fastq.gz --outdir /data/$USER/

Version:
  0.0.1
EOF
}


# Functions
function err() { cat <<< "$@" 1>&2; }
function fatal() { cat <<< "$@" 1>&2; usage; exit 1; }
function abspath() { readlink -e "$1"; }
function parser() {
  # Adds parsed command-line args to GLOBAL $Arguments associative array
  # + KEYS = short_cli_flag ("r1", "r2", "o", ...)
  # + VALUES = parsed_user_value ("WT_1.fastq" "WT_2.fastq", "/data/$USER", ...)
  # @INPUT "$@" = user command-line arguments
  # @CALLS check() to see if the user provided all the required arguments

  while [[ $# -gt 0 ]]; do
    key="$1"
    case $key in
      -h  | --help) usage && exit 0;;
      -r1 | --read-1)  provided "$key" "${2:-}"; Arguments["r1"]="$2"; shift; shift;;
      -r2 | --read-2)  provided "$key" "${2:-}"; Arguments["r2"]="$2"; shift; shift;;
      -o  | --outdir)  provided "$key" "${2:-}"; Arguments["o"]="$2"; shift; shift;;
      -t  | --threads) provided "$key" "${2:-}"; Arguments["t"]="$2"; shift; shift;;
      -g  | --gtf)     provided "$key" "${2:-}"; Arguments["g"]="$2"; shift; shift;;
      -p  | --prior)   provided "$key" "${2:-}"; Arguments["p"]="$2"; shift; shift;;
      -m  | --max-iter) provided "$key" "${2:-}"; Arguments["m"]="$2"; shift; shift;;
      -b  | --basename) provided "$key" "${2:-}"; Arguments["b"]="$2"; shift; shift;;
      -*  | --*) err "Error: Failed to parse unsupported argument: '${key}'."; usage && exit 1;;
      *) err "Error: Failed to parse unrecognized argument: '${key}'. Do any of your inputs have spaces?"; usage && exit 1;;
    esac
  done

  # Check for required args
  check
}


function provided() {
  # Checks to see if the argument's value exists
  # @INPUT $1 = name of user provided argument
  # @INPUT $2 = value of user provided argument
  # @CALLS fatal() if value is empty string or NULL

  if [[ -z "${2:-}" ]]; then
     fatal "Fatal: Failed to provide value to '${1}'!";
  fi
}


function clean(){
  # Finds the base name of the sample
  # @INPUT $1 = From optional basename argument
  # @RETURN $bname = cleaned base name (PATH and EXT removed)

  local bname=${1:-}
  local exts=("_1.fastq" "_2.fastq" ".R1.fastq" ".R2.fastq" "_R1.fastq" "_R2.fastq")

  if [[ -z "$bname" ]]; then
     bname="${Arguments[r1]}"  # Determine base name from R1 input
  fi

  # Remove PATH and .gz extension
  bname=$(basename $bname | sed 's/.gz$//g')

  # Clean remaining extensions (MateInfo + )
  for ext in "${exts[@]}"; do
    if [[ $bname == *${ext} ]]; then
      bname=$(echo "$bname" | sed "s@$ext\$@@")
      break # only remove one extension
    fi
  done

  echo "$bname"
}


function check(){
  # Checks to see if user provided required arguments
  # @INPUTS $Arguments = Global Associative Array
  # @CALLS fatal() if user did NOT provide all the $required args

  # List of required arguments
  local required=("r1" "r2" "o")
  #echo -e "Provided Required Inputs"
  for arg in "${required[@]}"; do
    value=${Arguments[${arg}]:-}
    if [[ -z "${value}" ]]; then
      fatal "Failed to provide all required args.. missing ${arg}"
    fi
  done
}


function initialize(){
  # Step 1. Initialize ouput directory, make outdir and copy over pipeline resources
  # @INPUT $1 = Read1 FastQ file
  # @INPUT $2 = Read2 FastQ file
  # @INPUT $3 = Pipeline outdir
  # @INPUT $4 = HERVx repository

  local r1=$(abspath "$1")
  local r2=$(abspath "$2")

  # Initialize pipeline output directory
  mkdir -p "$3" || fatal "Failed to create output directory '${3}', please check permissions before proceeding again."
  ln -sf "$r1" "$3"
  ln -sf "$r2" "$3"
  cp -R "$4"/wdl "$3"
  cp -R "$4"/config "$3"

}


function setup(){
  # Step 2. Setups pipeline for cromwell execution, creates new pipeline inputs JSON file
  # @INPUT $1 = Template inputs.json
  # @INPUT $2 = Pipeline outdir

  # Copy and/or over-write previous input JSON file template
  local template="${2}/config/inputs.json"
  local fullpath="$(abspath "${2}")"
  cp "$1" "${template}"

  # Update pointer to input FastQ files to use symlinks in Pipeline outdir
  Arguments[r1]="${fullpath}/$(basename "${Arguments[r1]}")"
  Arguments[r2]="${fullpath}/$(basename "${Arguments[r2]}")"

  # Create new inputs.json using the template
  for key in "${!Arguments[@]}"; do value=${Arguments["$key"]};
    sed -i "s|__${key}__|${value}|" "$template";
  done
}


function orchestrate(){
  # Step 3. Run WDL HERVx pipeline with specified cromwell executor backend
  # @INPUT $1 = Cromwell Executor backend, either local or slurm
  # @INPUT $2 = Pipeline outdir

  module load cromwell 2> /dev/null || fatal "Fail to load 'cromwell/52', not installed on target system."

  # Cromwell executor
  executor=${1:-local}

  # Goto Pipeline Ouput directory
  # Create a local singularity cache in output directory
  # cache can be used instead of re-pulling from DockerHub everytime
  cd "$2" && export SINGULARITY_CACHEDIR="${PWD}/.singularity"
  mkdir -p $SINGULARITY_CACHEDIR

  # unsetting XDG_RUNTIME_DIR to avoid some unsighly but harmless warnings
  unset XDG_RUNTIME_DIR

  # Run WDL workflow with specified cromwell executor backend
  case "$executor" in
    local) java -Dconfig.file=config/local-singularity.conf -jar ${CROMWELL_JAR} run -i config/inputs.json wdl/hervx.wdl
        ;;
    slurm) java -Dconfig.file=config/slurm-singularity.conf -jar ${CROMWELL_JAR} run -i config/inputs.json wdl/hervx.wdl
        ;;
    *)  echo "${executor} is not available." && \
        fatal "Failed to provide valid cromwell backend: ${executor}. Please use either local or slurm."
        ;;
  esac
}

function main(){
  # Parses args and runs initializes output directory, setups cromwell, and runs pipeline
  # @INPUT "$@" = command-line arguments
  # @CALLS parser(), initialize(), setup(), cromwell()

  if [ $# -eq 0 ]; then usage; exit 1; fi

  # Associative array to store parsed args
  declare -Ag Arguments

  # Positional Argument for Cromwell Executor Backend
  case $1 in
    local | slurm) Arguments["e"]="$1";;
    -h    | --help | help) usage && exit 0;;
    -*    | --*) err "Error: Failed to provide required positional argument: <local|slurm>."; usage && exit 1;;
    *) err "Error: Failed to provide valid positional argument. '${1}' is not supported. Valid options are local or slurm"; usage && exit 1;;
  esac

  # Parses remaining user provided command-line arguments
  parser "${@:2}" # Remove first item of list

  # Setting defaults for non-required arguments
  Arguments[b]="$(clean "${Arguments[b]:-}")"
  Arguments[o]="$(abspath "${Arguments[o]%/}")" # clean outdir path (remove trailing '/')
  Arguments[t]="${Arguments[t]:-2}"             # default is 2 threads
  Arguments[g]="${Arguments[g]:-HERV_rmsk.hg38.v2.transcripts.gtf}" # default = HERV_rmsk.hg38.v2.transcripts.gtf
  Arguments[p]="${Arguments[p]:-200000}"                            # default = 200000
  Arguments[m]="${Arguments[m]:-200}"                               # default = 200

  # Print pipeline metadata prior to running
  hervxrepo=$(abspath $(dirname  "$0"))
  echo -e "HERVx (version 0.0.1)\t$(date)\t${hervxrepo}"
  echo -e "Running pipeline with the following parameters:"
  for key in "${!Arguments[@]}"; do echo -e "\t-${key}\t${Arguments["$key"]}"; done

  # Step 1. Initialize output directory
  initialize "${Arguments[r1]}" "${Arguments[r2]}" "${Arguments[o]}" "${hervxrepo}"

  # Step 2. Setup pipeline for cromwell execution
  setup "${hervxrepo}/templates/inputs.json" "${Arguments[o]}"

  # Step 3. Run WDL pipeline with cromwell
  orchestrate "${Arguments[e]}" "${Arguments[o]}"
}


# Main: check usage, parse args, and run HERVx pipeline
main "$@"