partition-before-pggb

#!/usr/bin/env bash

# exit when any command fails
set -eo pipefail


# wfmash's default values
SEGMENT_LENGTH=5000
MAP_PCT_ID=90
MASH_KMER=19
MASH_KMER_THRES=0.001

# wfmash's parameters
input_fasta=false
segment_length=$SEGMENT_LENGTH
block_length=false
map_pct_id=$MAP_PCT_ID
n_mappings=1
no_splits=false
sparse_map=false
mash_kmer=$MASH_KMER
mash_kmer_thres=$MASH_KMER_THRES
hg_filter_ani_diff=false
exclude_delim="#"

# seqwish's default values
MIN_MATCH_LENGTH=23
SPARSE_FACTOR=0
TRANSCLOSE_BATCH=10M

# seqwish's parameters
min_match_length=$MIN_MATCH_LENGTH
sparse_factor=$SPARSE_FACTOR
transclose_batch=$TRANSCLOSE_BATCH

# smoothxg's default values
MAX_PATH_JUMP=0
MAX_EDGE_JUMP=0
TARGET_POA_LENGTH=700,900,1100
POA_PADDING=0.001
PAD_MAX_DEPTH=100
CONSENSUS_PREFIX=Consensus_

# smoothxg's parameters
skip_normalization=false
n_haps=false
max_path_jump=$MAX_PATH_JUMP
max_edge_jump=$MAX_EDGE_JUMP
target_poa_length=$TARGET_POA_LENGTH
poa_params=false
poa_padding=$POA_PADDING
pad_max_depth=$PAD_MAX_DEPTH
run_abpoa=false
run_global_poa=false
write_maf=false
consensus_spec=false # Disabled due to PGGB's #133 and #182 issues
consensus_prefix=$CONSENSUS_PREFIX

# odgi's parameters
do_viz=true
do_layout=true
do_stats=false

# vg's parameter
vcf_spec=false

# multiqc's parameter
multiqc=false

# default values
THREADS=$(getconf _NPROCESSORS_ONLN 2>/dev/null || getconf NPROCESSORS_ONLN 2>/dev/null || echo 1)

# general parameters
output_dir=false
input_temp_dir=false
input_paf=false
resume=false
threads=$THREADS
poa_threads=0
keep_intermediate_files=false
compress=false
names_with_params=false
show_version=false
show_help=false

# not exposed parameters
no_merge_segments=false
block_ratio_min=0
reduce_redundancy=true

# function that parses metric suffix numerical arguments
parse_numeric() {
  local input=$1
  local value=${input%[KkMmGgTt]*}
  local suffix=${input#$value}
  if [[ -z $suffix ]]; then
    if [[ $value =~ ^[0-9]+(\.[0-9]+)?$ ]]; then
      echo $value
      return 0
    else
      echo "[pggb] ERROR: Invalid input" >&2
      return 1
    fi
  fi
  case $suffix in
    K|k) value=$(echo "$value * 1000" | bc) ;;
    M|m) value=$(echo "$value * 1000000" | bc) ;;
    G|g) value=$(echo "$value * 1000000000" | bc) ;;
    T|t) value=$(echo "$value * 1000000000000" | bc) ;;
    *) echo "[pggb] ERROR: Invalid suffix or unsupported suffix. Supported metric suffixes are k, K, m, M, g, G, t, T." >&2; return 1 ;;
  esac
  printf "%.0f" $value
}

if [ $# -eq 0 ]; then
    show_help=true
fi

# read the options
cmd=$0" "$@
TEMP=`getopt -o i:o:D:a:p:c:g:s:l:K:F:k:x:f:B:Xn:j:P:O:Me:t:T:vhASY:G:Q:d:I:R:NbrmZzV: --long input-fasta:,output-dir:,temp-dir:,input-paf:,map-pct-id:,n-mappings:,hg-filter-ani-diff:,segment-length:,block-length-min:,mash-kmer:,mash-kmer-thres:,min-match-length:,sparse-map:,sparse-factor:,transclose-batch:,skip-normalization,n-haplotypes:,path-jump-max:,subpath-min:,edge-jump-max:,threads:,poa-threads:,skip-viz,do-layout,help,no-merge-segments,stats,exclude-delim:,poa-length-target:,poa-params:,poa-padding:,run-abpoa,global-poa,write-maf,consensus-spec:,consensus-prefix:,pad-max-depth:,block-id-min:,block-ratio-min:,no-splits,resume,keep-temp-files,multiqc,compress,names-with-params,vcf-spec:,version -n 'pggb' -- "$@"`
eval set -- "$TEMP"

# extract options and their arguments into variables.
while true ; do
    case "$1" in
        -i|--input-fasta) input_fasta=$2 ; shift 2 ;;
        -s|--segment-length) segment_length=$(parse_numeric $2) ; shift 2 ;;
        -l|--block-length) block_length=$(parse_numeric $2) ; shift 2 ;;
        -p|--map-pct-id) map_pct_id=$2 ; shift 2 ;;
        -c|--n-mappings) n_mappings=$2 ; shift 2 ;;
        -g|--hg-filter-ani-diff) hg_filter_ani_diff=$2 ; shift 2 ;;
        -N|--no-splits) no_splits=true ; shift ;;
        -x|--sparse-map) sparse_map=$2 ; shift 2 ;;
        -K|--mash-kmer) mash_kmer=$2 ; shift 2 ;;
        -F|--mash-kmer-thres) mash_kmer_thres=$2 ; shift 2 ;;
        -Y|--exclude-delim) exclude_delim=$2 ; shift 2 ;;
        -k|--min-match-length) min_match_length=$2 ; shift 2 ;;
        -f|--sparse-factor) sparse_factor=$2 ; shift 2 ;;
        -B|--transclose-batch) transclose_batch=$(parse_numeric $2) ; shift 2 ;;
        -X|--skip-normalization) skip_normalization=true ; shift ;;
        -n|--n-haplotypes) n_haps=$2 ; shift 2 ;;
        -j|--path-jump-max) max_path_jump=$2 ; shift 2 ;;
        -e|--edge-jump-max) max_edge_jump=$2 ; shift 2 ;;
        -G|--poa-length-target) target_poa_length=$2 ; shift 2 ;;
        -P|--poa-params) poa_params=$2 ; shift 2 ;;
        -O|--poa-padding) poa_padding=$2 ; shift 2 ;;
        -d|--pad-max-depth) pad_max_depth=$2 ; shift 2 ;;
        -b|--run-abpoa) run_abpoa=true ; shift ;;
        -z|--global-poa) run_global_poa=true ; shift ;;
        -M|--write-maf) write_maf=true ; shift ;;
        #-C|--consensus-spec) consensus_spec=$2 ; shift 2 ;;
        -Q|--consensus-prefix) consensus_prefix=$2 ; shift 2 ;;
        -v|--skip-viz) do_viz=false ; do_layout=false; shift ;;
        -S|--stats) do_stats=true ; shift ;;
        -V|--vcf-spec) vcf_spec=$2 ; shift 2 ;;
        -m|--multiqc) multiqc=true ; shift ;;
        -o|--output-dir) output_dir=$2 ; shift 2 ;;
        -D|--temp-dir) input_temp_dir=$2 ; shift 2 ;;
        -a|--input-paf) input_paf=$2 ; shift 2 ;;
        -r|--resume) resume=true ; shift ;;
        -t|--threads) threads=$2 ; shift 2 ;;
        -T|--poa-threads) poa_threads=$2 ; shift 2 ;;
        -A|--keep-temp-files) keep_intermediate_files=true ; shift ;;
        -Z|--compress) compress=true ; shift ;;
        --names-with-params) names_with_params=true ; shift ;;
        --version) show_version=true ; shift ;;
        -h|--help) show_help=true ; shift ;;
        --) shift ; break ;;
        *) echo "$2" "Internal error!" ; exit 1 ;;
    esac
done

PGGB_RELEASE_VERSION="v0.7.2"
SCRIPT_DIR=$( cd -- "$(dirname -- "$(readlink -f "${BASH_SOURCE[0]}" )" )" &> /dev/null && pwd )
GIT_VERSION=$( cd "$SCRIPT_DIR"; git describe --always --tags --long 2> /dev/null || echo $PGGB_RELEASE_VERSION )

if [ "$show_version" == true ]; then
    echo "pggb $GIT_VERSION"
    cd - &> /dev/null
    exit
fi

# Check input parameters
if [ "$show_help" == true ]; then
    padding=$(printf %${#0}s) # prints as many spaces as the length of $0
    echo "usage: $0 -i <input-fasta> -o <output-dir> [options]"
    echo "options:"
    echo "   [wfmash]"
    echo "    -i, --input-fasta FILE      input FASTA/FASTQ file"
    echo "    -s, --segment-length N      segment length for mapping [default: "$SEGMENT_LENGTH"]"
    echo "    -l, --block-length N        minimum block length filter for mapping [default: 5*segment-length]"
    echo "    -p, --map-pct-id PCT        percent identity for mapping/alignment [default: "$MAP_PCT_ID"]"
    echo "    -c, --n-mappings N          number of mappings for each segment [default: 1]"
    echo "    -g, --hg-filter-ani-diff N  filter out mappings unlikely to be N% less than the best mapping [default: "$HG_FILTER_ANI_DIFF"%]"
    echo "    -N, --no-splits             disable splitting of input sequences during mapping [default: enabled]"
    echo "    -x, --sparse-map N          keep this fraction of mappings ('auto' for giant component heuristic) [default: 1.0]"
    echo "    -K, --mash-kmer N           kmer size for mapping [default: "$MASH_KMER"]"
    echo "    -F, --mash-kmer-thres N     ignore the top % most-frequent kmers [default: "$MASH_KMER_THRES"]"
    echo "    -Y, --exclude-delim C       skip mappings when the query and target have the same"
    echo "                                prefix before the last occurrence of the given character C [default: #, assuming PanSN-spec]"
    echo "   [seqwish]"
    echo "    -k, --min-match-len N       filter exact matches below this length [default: "$MIN_MATCH_LENGTH"]"
    echo "    -f, --sparse-factor N       keep this randomly selected fraction of input matches [default: no sparsification]"
    echo "    -B, --transclose-batch      number of bp to use for transitive closure batch [default: "$TRANSCLOSE_BATCH"]"
    echo "   [smoothxg]"
    echo "    -X, --skip-normalization    do not normalize the final graph [default: normalize the graph]"
    echo "    -n, --n-haplotypes N        number of haplotypes"
    echo "    -j, --path-jump-max         maximum path jump to include in block [default: "$MAX_PATH_JUMP"]"
    echo "    -e, --edge-jump-max N       maximum edge jump before breaking [default: "$MAX_EDGE_JUMP"]"
    echo "    -G, --poa-length-target N,M target sequence length for POA, one per pass [default: "$TARGET_POA_LENGTH"]"
    echo "    -P, --poa-params PARAMS     score parameters for POA in the form of match,mismatch,gap1,ext1,gap2,ext2"
    echo "                                may also be given as presets: asm5, asm10, asm15, asm20"
    echo "                                [default: 1,19,39,3,81,1 = asm5]"
    echo "    -O, --poa-padding N         pad each end of each sequence in POA with N*(mean_seq_len) bp [default: "$POA_PADDING"]"
    echo "    -d, --pad-max-depth N       depth/haplotype at which we don't pad the POA problem [default: "$PAD_MAX_DEPTH"]"
    echo "    -b, --run-abpoa             run abPOA [default: SPOA]"
    echo "    -z, --global-poa            run the POA in global mode [default: local mode]"
    echo "    -M, --write-maf             write MAF output representing merged POA blocks [default: off]"
    echo "    -Q, --consensus-prefix P    use this prefix for consensus path names [default: "$CONSENSUS_PREFIX"]"
    #echo "    -C, --consensus-spec SPEC   consensus graph specification: write consensus graphs to"
    #echo "                                BASENAME.cons_[spec].gfa; where each spec contains at least a min_len parameter"
    #echo "                                (which defines the length of divergences from consensus paths to preserve in the"
    #echo "                                output), optionally a file containing reference paths to preserve in the output,"
    #echo "                                a flag (y/n) indicating whether we should also use the POA consensus paths, a"
    #echo "                                minimum coverage of consensus paths to retain (min_cov), and a maximum allele"
    #echo "                                length (max_len, defaults to 1e6); implies -a; example:"
    #echo "                                cons,100,1000:refs1.txt:n,1000:refs2.txt:y:2.3:1000000,10000"
    #echo "                                [default: off]"
    echo "   [odgi]"
    echo "    -v, --skip-viz              don't render visualizations of the graph in 1D and 2D [default: make them]"
    echo "    -S, --stats                 generate statistics of the seqwish and smoothxg graph [default: off]"
    echo "   [vg]"
    echo "    -V, --vcf-spec SPEC         specify a set of VCFs to produce with SPEC = REF[:LEN][,REF[:LEN]]*"
    echo "                                the paths matching ^REF are used as a reference, while the sample haplotypes"
    echo "                                are derived from path names, assuming they match the PanSN; e.g. '-V chm13',"
    echo "                                a path named HG002#1#ctg would be assigned to sample HG002 phase 1."
    echo "                                If LEN is specified and greater than 0, the VCFs are decomposed, filtering "
    echo "                                sites whose max allele length is greater than LEN. [default: off]"
    echo "   [multiqc]"
    echo "    -m, --multiqc               generate MultiQC report of graphs' statistics and visualizations,"
    echo "                                automatically runs odgi stats [default: off]"
    echo "   [general]"
    echo "    -o, --output-dir PATH       output directory"
    echo "    -D, --temp-dir PATH         directory for temporary files"
    echo "    -a, --input-paf FILE        input PAF file; the wfmash alignment step is skipped"
    echo "    -r, --resume                do not overwrite existing outputs in the given directory"
    echo "                                [default: start pipeline from scratch]"
    echo "    -t, --threads N             number of compute threads to use in parallel steps [default: "$threads"]"
    echo "    -T, --poa-threads N         number of compute threads to use during POA (set lower if you OOM during smoothing)"
    echo "    -A, --keep-temp-files       keep intermediate graphs"
    echo "    -Z, --compress              compress alignment (.paf), graph (.gfa, .og), and MSA (.maf) outputs with pigz,"
    echo "                                and variant (.vcf) outputs with bgzip"
    echo "    --names-with-params         put parameter values in filenames, instead of hashes"
    echo "    --version                   display the version of pggb"
    echo "    -h, --help                  this text"
    echo
    echo "Use wfmash, seqwish, smoothxg, odgi, gfaffix, and vg to build, project and display a pangenome graph."
    exit
fi

# Function to check if a tool is in PATH
check_tool_availability() {
  local tool_name=$1
  local tool_command=$2

  if ! command -v $tool_command &> /dev/null; then
    # Instead of exiting immediately, add the tool name to the array of missing tools
    missing_tools+=("$tool_name")
  fi
}

# Check availability of each tool
check_tool_availability "wfmash" "wfmash"
check_tool_availability "seqwish" "seqwish"
check_tool_availability "smoothxg" "smoothxg"
check_tool_availability "odgi" "odgi"
check_tool_availability "gfaffix" "gfaffix"
if [ "$vcf_spec" != "false" ]; then
check_tool_availability "vg" "vg"
check_tool_availability "vcfbub" "vcfbub"
check_tool_availability "vcfwave" "vcfwave"
check_tool_availability "bcftools" "bcftools"
fi
if [[ $multiqc == true ]]; then
check_tool_availability "MultiQC" "multiqc"
fi

# If there are missing tools, report them and exit
if [ ${#missing_tools[@]} -ne 0 ]; then
  echo "[pggb] ERROR: the following tools are not installed or not in the PATH:"
  for tool in "${missing_tools[@]}"; do
    echo "- $tool"
  done
  exit 1
fi

if [ "$input_fasta" = "false" ] || [ "$output_dir" = "false" ]; then
    >&2 echo "[pggb] ERROR: mandatory argument: -i/--input-fasta and -o/--output-dir"
    exit
elif [ ! -f "${input_fasta}.fai" ]; then
    echo "[pggb] ERROR: index for $input_fasta does not exist. Please create it using 'samtools faidx $input_fasta'."
    exit 1
elif [ "$n_mappings" -lt 1 ]; then
    >&2 echo "[pggb] ERROR: -c/--n-mappings must be greater than or equal to 1"
    exit
fi

# Check Pangenome Sequence Naming (PanSN)
pansn_not_respected=false
while IFS= read -r line; do
  if [[ ! $line =~ ^([^#]+#)[0-9]+#[^#]+$ ]] && [[ "$pansn_not_respected" == "false" ]]; then
    pansn_not_respected=$line
    break
  fi
done < <(cut -f 1 "${input_fasta}.fai")
if [ "$pansn_not_respected" != "false" ]; then
  >&2 echo "[pggb] warning: there are sequence names (like '$pansn_not_respected') that do not match the Pangenome Sequence Naming (PanSN)."
  
  pansn_error=false
  if [ "$vcf_spec" != "false" ]; then
      pansn_error=true
      >&2 echo "[pggb] ERROR: -V/--vcf-spec cannot be used if the Pangenome Sequence Naming (PanSN) is not respected."
  fi
  if [ "$n_haps" = "false" ] || [ "$n_haps" -lt 1 ]; then
      pansn_error=true
      >&2 echo "[pggb] ERROR: -n/--n-haplotypes must be greater than or equal to 1 when the Pangenome Sequence Naming (PanSN) is not respected."
  fi

  if [ $pansn_error = true ]; then
    exit
  fi
fi

# Alignment
if [[ $input_paf == false ]]; then
  mapper=wfmash
  mapper_version=$(wfmash --version 2>&1)
  mapper_letter='W'
else
  mapper=external
  mapper_version=unknown
  mapper_letter='E'
fi

if [[ $block_length == false ]]; then
    block_length=$(echo "$segment_length * 5" | bc)
fi

if [[ $hg_filter_ani_diff == false ]]; then
    # Convert map_pct_id to a format that bc can handle if it's not already a decimal
    map_pct_id_decimal=$(printf "%.2f" "$map_pct_id")

    # Using bc for floating-point comparison
    if [[ $(echo "$map_pct_id_decimal >= 90" | bc) -eq 1 ]]; then
        hg_filter_ani_diff=30
    elif [[ $(echo "$map_pct_id_decimal >= 80 && $map_pct_id_decimal < 90" | bc) -eq 1 ]]; then
        hg_filter_ani_diff=10
    else
        hg_filter_ani_diff=0
    fi
fi

paf_spec=$mapper_letter-s$segment_length-l$block_length-p$map_pct_id-n$n_mappings-K$mash_kmer-F$mash_kmer_thres-x$sparse_map-g$hg_filter_ani_diff

split_cmd=""
if [[ $no_splits == true ]]; then
    split_cmd=-N
    paf_spec="$paf_spec"-N
fi

# Normalization ($n_haps is checked in this part of the script because it is also used for the 'auto' mapping sparsification)
if [[ $n_haps == false ]]; then
    n_haps=$(cut -f 1 "${input_fasta}.fai" | cut -f 1,2 -d '#' | sort | uniq | wc -l)
fi

sparse_map_cmd=""
if [[ $sparse_map == "auto" ]]; then
    # set sparse mapping using giant component heuristic
    # we keep 10x log(n)/n mappings
    # if this is < 1, otherwise we keep all
    n=$n_haps
    sparse_map_frac=$(echo "x=l($n)/$n * 10; if (x < 1) { x } else { 1 }"  | bc -l | cut -c -8)
    sparse_map_cmd="-x $sparse_map_frac"
elif [[ $sparse_map != false ]]; then
    sparse_map_cmd="-x $sparse_map"
fi

merge_cmd=""
if [[ $no_merge_segments == true ]]; then
    merge_cmd=-M
    paf_spec="$paf_spec"-M
fi

if [[ $exclude_delim != false ]]; then
    exclude_cmd="-Y "$exclude_delim
    paf_spec="$paf_spec"-Y
else
    exclude_cmd=-X
    paf_spec="$paf_spec"-X
fi

if [[ "$input_paf" == false ]]; then
    if [[ "$names_with_params" == true ]]; then
        prefix_paf="$input_fasta"."$paf_spec"
    else
        prefix_paf="$input_fasta".$(echo "$paf_spec" | sha256sum | head -c 7)
    fi
else
    prefix_paf="$input_paf"
fi


# Graph induction
if [[ "$names_with_params" == true ]]; then
    prefix_seqwish="$prefix_paf".k$min_match_length-f$sparse_factor-B$transclose_batch
else
    prefix_seqwish="$prefix_paf".$(echo k$min_match_length-f$sparse_factor-B$transclose_batch | sha256sum | head -c 7)
fi

# poa param suggestions from minimap2
# - asm5, --poa-params 1,19,39,3,81,1, ~0.1 divergence
# - asm10, --poa-params 1,9,16,2,41,1, ~1 divergence
# - asm20, --poa-params 1,4,6,2,26,1, ~5% divergence
# between asm10 and asm20 ~ 1,7,11,2,33,1
poa_params_cmd=""
if [[ $poa_params == false ]]; then
    poa_params_cmd="-p 1,19,39,3,81,1"
else
    if [[ $poa_params == "asm5" ]]; then
        poa_params_cmd="-p 1,19,39,3,81,1"
    elif [[ $poa_params == "asm10" ]]; then
        poa_params_cmd="-p 1,9,16,2,41,1"
    elif [[ $poa_params == "asm15" ]]; then
        poa_params_cmd="-p 1,7,11,2,33,1"
    elif [[ $poa_params == "asm20" ]]; then
        poa_params_cmd="-p 1,4,6,2,26,1"
    else
        poa_params_cmd="-p $poa_params"
    fi
fi

block_id_min=$(echo "scale=4; $map_pct_id / 100.0" | bc)
if [[ "$names_with_params" == true ]]; then
    prefix_smoothed="$prefix_seqwish".h$n_haps-G$target_poa_length-j$max_path_jump-e$max_edge_jump-d$pad_max_depth-I$block_id_min-R$block_ratio_min-p$poa_params-O$poa_padding.smooth
else
    prefix_smoothed="$prefix_seqwish".$(echo h$n_haps-G$target_poa_length-j$max_path_jump-e$max_edge_jump-d$pad_max_depth-I$block_id_min-R$block_ratio_min-p$poa_params-O$poa_padding | sha256sum | head -c 7).smooth
fi
prefix_smoothed_output="$prefix_smoothed"


fmt="%C\n%Us user %Ss system %P cpu %es total %MKb max memory"
timer=$(which time)


# Directories
if [ ! -e "$output_dir" ]; then
  mkdir "$output_dir"
fi
prefix_paf="$output_dir"/$(basename "$prefix_paf")
prefix_smoothed_output="$output_dir"/$(basename "$prefix_smoothed")

# If the temporary directory is not explicitly set, set it equal to the output directory
if [[ "$input_temp_dir" == false ]]; then
  temp_dir="$output_dir" # NOTE: the output directory always exists before this statement
else
  temp_dir="$input_temp_dir"
fi
temp_dir_was_created=false
if [ ! -e "$temp_dir" ]; then
  mkdir "$temp_dir"
  temp_dir_was_created=true
fi
prefix_mappings_paf="$temp_dir"/$(basename "$prefix_paf")
prefix_seqwish="$temp_dir"/$(basename "$prefix_seqwish")
prefix_smoothed="$temp_dir"/$(basename "$prefix_smoothed")

if [[ $poa_threads == 0 ]]; then
    poa_threads=$threads
fi

date=`date "+%m-%d-%Y_%H:%M:%S"`
log_file="$prefix_smoothed_output".$date.log
param_file="$prefix_smoothed_output".$date.params.yml

# write parameters to log_file:
echo -e "Starting pggb on $date\n" > "$log_file"
echo -e "Command: $cmd\n" >> "$log_file"
echo -e "PARAMETERS\n" >> "$log_file"
cat <<EOT | tee -a "$log_file" "$param_file" >/dev/null
general:
  input-fasta:        $input_fasta
  output-dir:         $output_dir
  temp-dir:           $temp_dir
  resume:             $resume
  compress:           $compress
  threads:            $threads
  poa_threads:        $poa_threads
pggb:
  version:            $GIT_VERSION
$mapper:
  version:            $mapper_version
  segment-length:     $segment_length
  block-length:       $block_length
  map-pct-id:         $map_pct_id
  n-mappings:         $n_mappings
  no-splits:          $no_splits
  sparse-map:         $sparse_map
  mash-kmer:          $mash_kmer
  mash-kmer-thres:    $mash_kmer_thres
  hg-filter-ani-diff: $hg_filter_ani_diff
  exclude-delim:      $exclude_delim
  no-merge-segments:  $no_merge_segments
seqwish:
  version:            $(seqwish --version 2>&1)
  min-match-len:      $min_match_length
  sparse-factor:      $sparse_factor
  transclose-batch:   $transclose_batch
smoothxg:
  version:            $(smoothxg --version 2>&1)
  skip-normalization: $skip_normalization
  n-haplotypes:       $n_haps
  path-jump-max:      $max_path_jump
  edge-jump-max:      $max_edge_jump
  poa-length-target:  $target_poa_length
  poa-params:         ${poa_params_cmd:3}
  poa_padding:        $poa_padding
  run_abpoa:          $run_abpoa
  run_global_poa:     $run_global_poa
  pad-max-depth:      $pad_max_depth
  write-maf:          $write_maf
  consensus-spec:     $consensus_spec
  consensus-prefix:   $consensus_prefix
  block-id-min:       $block_id_min
  block-ratio-min:    $block_ratio_min
odgi:
  version:            $(odgi version -v 2>&1)
  viz:                $do_viz
  layout:             $do_layout
  stats:              $do_stats
gfaffix:
  version:            v$(gfaffix --version | cut -f 2 -d ' ')
  reduce-redundancy:  $reduce_redundancy
vg:
  version:            $(vg version | head -n 1 | cut -f 3 -d ' ')
  deconstruct:        $vcf_spec
reporting:
  version:            v$(multiqc --version | cut -f 3 -d ' ')
  multiqc:            $multiqc
EOT

#-------------------------------------------------------------------------------
echo -e "\nRunning partitioning\n" >> "$log_file"
if [[ "$input_paf" == false ]]; then
  if [[ ! -s "$prefix_mappings_paf".mappings.$mapper.paf || $resume == false ]]; then
    ($timer -f "$fmt" wfmash \
        -s $segment_length \
        -l $block_length \
        -p $map_pct_id \
        -n $n_mappings \
        $split_cmd \
        $sparse_map_cmd \
        -k $mash_kmer \
        -H $mash_kmer_thres \
        $exclude_cmd \
        -t $threads \
        --tmp-base $temp_dir \
        $merge_cmd \
        "$input_fasta" \
        --lower-triangular \
        --hg-filter-ani-diff $hg_filter_ani_diff \
        --approx-map \
        > "$prefix_mappings_paf".mappings.$mapper.paf) 2> >(tee -a "$log_file")
    fi
fi

SCRIPT_DIR=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )

if [[ ! -s "$prefix_mappings_paf".mappings.$mapper.paf.edges.list.txt \
  || ! -s "$prefix_mappings_paf".mappings.$mapper.paf.edges.weights.txt \
  || ! -s "$prefix_mappings_paf".mappings.$mapper.paf.vertices.id2name.txt
  || $resume == false ]]; then
  ($timer -f "$fmt" python3 "$SCRIPT_DIR"/scripts/paf2net.py -p "$prefix_mappings_paf".mappings.$mapper.paf) 2> >(tee -a "$log_file")
fi

# Not perfect, but at least one community will exist
if [[ ! -s "$prefix_paf".community.0.txt  || $resume == false ]]; then
  ($timer -f "$fmt" python3 "$SCRIPT_DIR"/scripts/net2communities.py \
      -e "$prefix_mappings_paf".mappings.$mapper.paf.edges.list.txt \
      -w "$prefix_mappings_paf".mappings.$mapper.paf.edges.weights.txt \
      -n "$prefix_mappings_paf".mappings.$mapper.paf.vertices.id2name.txt \
      --accurate-detection --output-prefix "$prefix_paf") 2> >(tee -a "$log_file")
fi

ls "$prefix_paf".community.*.txt | while read community; do
  prefix_fasta="$output_dir"/$(basename "$community" .txt)
  if [[ $compress == true ]]; then
    if [[ ! -s "$prefix_fasta".fa.gz || $resume == false ]]; then
      ($timer -f "$fmt" samtools faidx "$input_fasta" $(cat "$community") | \
        bgzip -@ "$threads" -c > "$prefix_fasta".fa.gz) 2> >(tee -a "$log_file")
      ($timer -f "$fmt" samtools faidx "$prefix_fasta".fa.gz) 2> >(tee -a "$log_file")
    fi
  else
    if [[ ! -s "$prefix_fasta".fa || $resume == false ]]; then
      ($timer -f "$fmt" samtools faidx "$input_fasta" $(cat "$community") \
        > "$prefix_fasta".fa) 2> >(tee -a "$log_file")
      ($timer -f "$fmt" samtools faidx "$prefix_fasta".fa) 2> >(tee -a "$log_file")
    fi
  fi
done

# Command lines preparation
params=""
if [[ $no_splits != false ]]; then
  params="$params --no-splits"
fi
if [[ $sparse_map != false ]]; then
  params="$params -x $sparse_map"
fi
if [[ $exclude_delim != false ]]; then
  params="$params -Y \"$exclude_delim\""
fi
if [[ $skip_normalization != false ]]; then
  params="$params --skip-normalization"
fi
if [[ $run_abpoa != false ]]; then
  params="$params --run-abpoa"
fi
if [[ $run_global_poa != false ]]; then
  params="$params --global-poa"
fi
if [[ $write_maf != false ]]; then
  params="$params --write-maf"
fi
#consensus_spec=false # Disabled due to PGGB's #133 and #182 issues
if [[ $do_viz != true || $do_layout != true ]]; then
  params="$params --skip-viz"
fi
if [[ $run_abpoa != false ]]; then
  params="$params --do-stats"
fi
if [[ $vcf_spec != false ]]; then
  params="$params -V $vcf_spec"
fi
if [[ $multiqc != false ]]; then
  params="$params --multiqc"
fi
# the output directory will be set for each community
if [[ "$input_temp_dir" != false ]]; then
  # The temporary directory will be the same for each community
  params="$params --temp-dir $temp_dir"
#else when not specified, the temporary directory will be the output directory
fi

if [[ $input_paf != false ]]; then
  params="$params -a $input_paf"
fi
if [[ $resume != false ]]; then
  params="$params --resume"
fi
params="$params --threads $threads"
params="$params --poa-threads $poa_threads"
if [[ $keep_intermediate_files != false ]]; then
  params="$params --keep-temp-file"
fi
if [[ $compress != false ]]; then
  params="$params --compress"
fi



if [[ $compress == true ]]; then
  fasta_ext="fa.gz"
else
  fasta_ext="fa"
fi
echo ""
ls "$prefix_paf".community.*.$fasta_ext | while read community; do
cat <<EOT | tee -a "$log_file" "$param_file"
pggb -i $community \\
     -o $community.out \\
     -s $segment_length -l $block_length -p $map_pct_id -c $n_mappings -K $mash_kmer -F $mash_kmer_thres -g $hg_filter_ani_diff \\
     -k $min_match_length -f $sparse_factor -B $transclose_batch \\
     -n $n_haps -j $max_path_jump -e $max_edge_jump -G $target_poa_length $poa_params_cmd -O $poa_padding -d $pad_max_depth -Q $consensus_prefix \\
    $params
EOT
done
echo ""

if [[ $keep_intermediate_files == true ]]; then
  if [[ "$output_dir" != "$temp_dir" ]]; then
    # `|| true` to avoid `mv` fail if there are missing files to move
    mv -f "$prefix_mappings_paf".mappings.$mapper.paf "$output_dir" 2> /dev/null || true
    mv -f "$prefix_paf".community.*.txt "$output_dir" 2> /dev/null || true
    mv -f "$prefix_mappings_paf".mappings.$mapper.paf.*.txt "$output_dir" 2> /dev/null || true
  fi
else
  rm -f "$prefix_mappings_paf".mappings.$mapper.paf
  rm -f "$prefix_paf".community.*.txt
  rm -f "$prefix_mappings_paf".mappings.$mapper.paf.*.txt
fi

if [[ $temp_dir_was_created == true ]]; then
  rm -r "$temp_dir"
fi

if [[ $compress == true ]]; then
    # FASTA files have already been compressed
    if [[ $keep_intermediate_files == true ]]; then
      pigz -f -q -p $threads "$prefix_paf"*.paf -v
      ls "$output_dir"/*.txt | while read f; do pigz -f -q -p $threads "$f" -v; done
    fi
fi
#-------------------------------------------------------------------------------