INPUT="$1"   # input file
OUTPUT="$2"  # output file
REF="$3"     # reference genome file

# Make sure that it doesn't stop if the exit code is != 0
set +e

# Run picard
picard ValidateSamFile \
  INPUT="$INPUT" \
  OUTPUT="$OUTPUT" \
  REFERENCE_SEQUENCE="$REF" \
  MODE=SUMMARY

# Issue exit code
exitcode=$?
if [ $exitcode -eq 1 ]
then
    exit 1
else
    exit 0
fi

shell:
  "Rscript --vanilla workflow/scripts/cleanPhenotypes.R"

shell:
  "bash workflow/scripts/getReference.sh {params.release} data/external/reference/ >{log} 2>&1"

shell:
  "Rscript --vanilla workflow/scripts/data_processing_founders/getAccessionGenotypes.R --outdir data/external/founder_genotypes "
  ">{log} 2>&1"

shell:
  "Rscript --vanilla workflow/scripts/data_processing_founders/getMagicMosaics.R"
  ">{log} 2>&1"

shell:
  "Rscript --vanilla workflow/scripts/data_processing_founders/imputeMagicGeno.R"
  ">{log} 2>&1"

  shell:
    "(bwa mem -M -t {threads} "
    "{params.extra} "
	  "{input.ref} {input.reads} "
    " | "
    "samtools sort -T data/intermediate/mapped/.samtools_{wildcards.name}.{wildcards.unit} "
    "-o {output} - ) "
    "2>{log}"

shell:
  "(samtools merge -@ $(({threads} - 1)) {output.bam} {input}; "
  "samtools index -@ $(({threads} - 1)) {output.bam} {output.bai}) "
  ">{log} 2>&1"

shell:
  "(picard MarkDuplicates "
  "INPUT={input} "
  "OUTPUT={output.bam} "
  "METRICS_FILE={output.metrics} "
  "REMOVE_DUPLICATES=true; "

  shell:
    "gatk3 -T IndelRealigner "
    "-R {input.ref} "
    "--baq CALCULATE_AS_NECESSARY "
    "-I {input.bam} "
    "-o {output.bam} "
    "-targetIntervals {input.intervals} "
    "-noTags "
    ">{log} 2>&1; "
    "touch {output.bai}" # this is just to ensure the .bai has newer date than .bam


rule picard_ValidateSamFile:
  input:
    "data/intermediate/realigned/{name}.bam"

shell:
  "bash workflow/scripts/ValidateSamFile.sh {input} {output} {params.ref} >{log} 2>&1"

shell:
  '''
  END=$(cat {input.ref}.fai | awk -F'\\t' '$1=={wildcards.chrom}' | cut -f2)
  harp like \
    --bam {input.bam} \
    --refseq {input.ref} \
    --region "{wildcards.chrom}:1-$END" \
    --snps {input.snps} \
    --out logs/harp_like/{wildcards.name}_chrom{wildcards.chrom} \
    --stem data/intermediate/harp_like/{wildcards.name}_chrom{wildcards.chrom} \
    >{log} 2>&1
  '''

shell:
  '''
  END=$(cat {input.ref}.fai | awk -F'\\t' '$1=={wildcards.chrom}' | cut -f2)
  echo "Starting harp"
  harp freq \
    --hlk {input.hlk} \
    --region "{wildcards.chrom}:1-$END" \
    --window_width {params.winsize} \
    --window_step {params.winstep} \
    --stem {params.stem} \
    >{log} 2>&1
  echo "Finished harp"
  rm -r {params.stem}.output
  '''

shell:
  '''
  Rscript --vanilla workflow/scripts/compilePoolseqFreqs.R \
    --indir data/intermediate/harp_freq/ \
    --outdir data/processed/poolseq/ \
    --window {wildcards.window} \
    >{log} 2>&1
  '''

shell:
  '''
  Rscript --vanilla workflow/scripts/simulations/compileSimulations.R >{log} 2>&1
  '''