#!/bin/bash

# Call run_alphafold.sh with the specified flags
mkdir -p "/mnt/myprotein/output/output_$FLOAT_JOB_ID"
RESULT="/mnt/myprotein/output/output_$FLOAT_JOB_ID"

exec 1>>"$RESULT/stdout_job_$FLOAT_JOB_ID" 2>>"$RESULT/stderr_job_$FLOAT_JOB_ID"

function do_MSA() {
  python /app/alphafold/run_alphafold_msa.py \
    --fasta_paths=/mnt/myprotein/${PROTEIN_NAME} \
    --uniref90_database_path=/mnt/mydb/uniref90/uniref90.fasta \
    --mgnify_database_path=/mnt/mydb/mgnify/mgy_clusters_2022_05.fa \
    --data_dir=/mnt/mydb/ \
    --template_mmcif_dir=/mnt/mydb/pdb_mmcif/mmcif_files \
    --obsolete_pdbs_path=/mnt/mydb/pdb_mmcif/obsolete.dat \
    --pdb70_database_path=/mnt/mydb/pdb70/pdb70 \
    --bfd_database_path=/mnt/mydb/bfd/bfd_metaclust_clu_complete_id30_c90_final_seq.sorted_opt \
    --uniref30_database_path=/mnt/mydb/uniref30/UniRef30_2021_03 \
    --output_dir=$RESULT \
    --max_template_date=2022-01-01 \
    --db_preset=full_dbs \
    --model_preset=monomer \
    --benchmark=False \
    --use_precomputed_msas=False \
    --num_multimer_predictions_per_model=5 \
    --models_to_relax='best' \
    --use_gpu_relax=False \
    --logtostderr
}

function do_inference() {
  python /app/alphafold/run_alphafold.py \
    --fasta_paths=/mnt/myprotein/${PROTEIN_NAME} \
    --uniref90_database_path=/mnt/mydb/uniref90/uniref90.fasta \
    --mgnify_database_path=/mnt/mydb/mgnify/mgy_clusters_2022_05.fa \
    --data_dir=/mnt/mydb/ \
    --template_mmcif_dir=/mnt/mydb/pdb_mmcif/mmcif_files \
    --obsolete_pdbs_path=/mnt/mydb/pdb_mmcif/obsolete.dat \
    --pdb70_database_path=/mnt/mydb/pdb70/pdb70 \
    --bfd_database_path=/mnt/mydb/bfd/bfd_metaclust_clu_complete_id30_c90_final_seq.sorted_opt \
    --uniref30_database_path=/mnt/mydb/uniref30/UniRef30_2021_03 \
    --output_dir=$RESULT \
    --max_template_date=2022-01-01 \
    --db_preset=full_dbs \
    --model_preset=monomer \
    --benchmark=False \
    --use_precomputed_msas=True \
    --num_multimer_predictions_per_model=5 \
    --models_to_relax='best' \
    --use_gpu_relax=True \
    --logtostderr
}

if [ "${DO_STEP1}" = "True" ]; then
  echo "msa"
  do_MSA
  yes | /opt/memverge/bin/float migrate --rerun --sync -f -e DO_STEP1=False --instType g4dn.2xlarge
else
  echo "inference"
  do_inference
fi

Input sequence file (in .fasta format) location relative to your storage mount point. 
For example, input/7fc7_96.fasta

Input sequence file (in .fasta format) location relative to your storage mount point. 
For example, input/7ktt_1142.fasta

#run_wgs_model_sention30_s3.sh
RESULT="/mnt/mys3"
sudo rm -rf "$RESULT/output/output_$FLOAT_JOB_ID"
sudo mkdir -p "$RESULT/output/output_$FLOAT_JOB_ID"
RESULT="$RESULT/output/output_$FLOAT_JOB_ID"

echo "after mkdir" 

set -exvuo pipefail
if [ -z ${FLOAT_LOG_PATH+x} ]; then
    FLOAT_LOG_PATH="."
fi

exec > $FLOAT_LOG_PATH/run.log 2>&1

cd /data
export PATH="/opt/aws/dist":$PATH

ref_build="b37"
s3_bucket_basename="s3://omics-us-west-2/reference"
ref_idx_files=(
  ""
  ".fai"
  ".amb"
  ".ann"
  ".bwt"
  ".pac"
  ".sa"
)

case "$ref_build" in
  hg38_alt|hg38_gatk)
    ref_base="hg38/Homo_sapiens_assembly38.fasta"
    has_alt=true
    VCFS=(
      "hg38/dbsnp_138.hg38.vcf.gz"
      "hg38/Homo_sapiens_assembly38.known_indels.vcf.gz"
      "hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"
    )
    ;;
  hg38|hg38_noalt|hs38)
    ref_base="hg38/hs38.fa"
    has_alt=false
    VCFS=(
      "hg38/dbsnp_138.hg38.vcf.gz"
      "hg38/Homo_sapiens_assembly38.known_indels.vcf.gz"
      "hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"
    )
    ;;
  b37_gatk)
    ref_base="b37/Homo_sapiens_assembly19.fasta"
    has_alt=false
    VCFS=(
      "b37/dbsnp_138.b37.vcf.gz"
      "b37/1000G_phase1.indels.b37.vcf.gz"
      "b37/Mills_and_1000G_gold_standard.indels.b37.vcf.gz"
    )
    ;;
  b37|hs37d5)
    ref_base="b37/hs37d5.fa"
    has_alt=false
    VCFS=(
      "b37/dbsnp_138.b37.vcf.gz"
      "b37/1000G_phase1.indels.b37.vcf.gz"
      "b37/Mills_and_1000G_gold_standard.indels.b37.vcf.gz"
    )
    ;;
  hg19|ucsc_hg19)
    ref_base="hg19/ucsc.hg19.fasta"
    has_alt=false
    VCFS=(
      "hg19/dbsnp_138.hg19.vcf.gz"
      "hg19/1000G_phase1.indels.hg19.sites.vcf.gz"
      "hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz"
    )
    ;;
  quickstart)
    ref_base="quickstart/ucsc.hg19_chr22.fasta"
    has_alt=false
    VCFS=(
      "quickstart/dbsnp_135.hg19_chr22.vcf.gz"
      "quickstart/1000G_phase1.snps.high_confidence.hg19_chr22.sites.vcf.gz"
      "quickstart/Mills_and_1000G_gold_standard.indels.hg19_chr22.sites.vcf.gz"
    )
    ;;
  *)
    echo "ERROR: unknown genome build"
    exit 1
    ;;
esac

for idx in "${ref_idx_files[@]}"; do
  aws s3 cp "$s3_bucket_basename/$ref_base$idx" . --quiet
done
if [[ "$has_alt" == true ]]; then
  aws s3 cp "$s3_bucket_basename/${ref_base}.alt" . --quiet
fi
for i in $(seq 1 "${#VCFS[@]}"); do
  i=$((i - 1))
  vcf="${VCFS[$i]}"
  aws s3 cp "$s3_bucket_basename/$vcf" . --quiet
  aws s3 cp "$s3_bucket_basename/${vcf}.tbi" . --quiet
done

echo $(date) "[Task] download ref data"
start=`date +%s`
aws s3 cp s3://omics-us-west-2/workflows-scrnaseq-reference/igenomes/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.fasta . --quiet
aws s3 cp s3://omics-us-west-2/workflows-scrnaseq-reference/igenomes/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.fasta.fai . --quiet
aws s3 cp s3://omics-us-west-2/workflows-scrnaseq-reference/igenomes/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.dict . --quiet
aws s3 cp s3://omics-us-west-2/workflows-scrnaseq-reference/igenomes/Homo_sapiens/GATK/GRCh37/Sequence/BWAIndex/human_g1k_v37_decoy.fasta.amb . --quiet
aws s3 cp s3://omics-us-west-2/workflows-scrnaseq-reference/igenomes/Homo_sapiens/GATK/GRCh37/Sequence/BWAIndex/human_g1k_v37_decoy.fasta.ann . --quiet
aws s3 cp s3://omics-us-west-2/workflows-scrnaseq-reference/igenomes/Homo_sapiens/GATK/GRCh37/Sequence/BWAIndex/human_g1k_v37_decoy.fasta.bwt . --quiet
aws s3 cp s3://omics-us-west-2/workflows-scrnaseq-reference/igenomes/Homo_sapiens/GATK/GRCh37/Sequence/BWAIndex/human_g1k_v37_decoy.fasta.pac . --quiet
aws s3 cp s3://omics-us-west-2/workflows-scrnaseq-reference/igenomes/Homo_sapiens/GATK/GRCh37/Sequence/BWAIndex/human_g1k_v37_decoy.fasta.sa . --quiet

end=`date +%s`
echo $(date) "[Task End] time used: " $((end-start)) "sec"

## download 1.fastq.gz and 2.fastq.gz
# wget https://omics-us-east-1.s3.amazonaws.com/sample-inputs/4414139/HG002.novaseq.pcr-free.32x.R1.fastq.gz --no-check-certificate -O 1.fastq.gz --quiet
# wget https://omics-us-east-1.s3.amazonaws.com/sample-inputs/4414139/HG002.novaseq.pcr-free.32x.R2.fastq.gz --no-check-certificate -O 2.fastq.gz --quiet

#wget https://storage.googleapis.com/brain-genomics-public/research/sequencing/fastq/novaseq/wgs_pcr_free/30x/HG003.novaseq.pcr-free.30x.R1.fastq.gz --no-check-certificate -O 1.fastq.gz --quiet
#wget https://storage.googleapis.com/brain-genomics-public/research/sequencing/fastq/novaseq/wgs_pcr_free/30x/HG003.novaseq.pcr-free.30x.R2.fastq.gz --no-check-certificate -O 2.fastq.gz --quiet

#wget https://storage.googleapis.com/brain-genomics-public/research/sequencing/fastq/novaseq/wgs_pcr_free/50x/HG007.novaseq.pcr-free.50x.R1.fastq.gz --no-check-certificate -O 1.fastq.gz --quiet
#wget https://storage.googleapis.com/brain-genomics-public/research/sequencing/fastq/novaseq/wgs_pcr_free/50x/HG007.novaseq.pcr-free.50x.R2.fastq.gz --no-check-certificate -O 2.fastq.gz --quiet

# *******************************************
# Script to perform DNAscope variant calling
# using a single sample with fastq files
# named 1.fastq.gz and 2.fastq.gz
# *******************************************

# Update with the fullpath location of your sample fastq
SM="sample" #sample name
RGID="rg_$SM" #read group ID
PL="ILLUMINA" #or other sequencing platform
FASTQ_FOLDER="/data"
FASTQ_1="/mnt/mys3/$FILE1"
FASTQ_2="/mnt/mys3/$FILE2" #If using Illumina paired data

# Update with the location to the DNAscope model file
# Model files can be found at, https://github.com/Sentieon/sentieon-models
DNASCOPE_MODEL=/models/DNAscopeIlluminaWGS2.0.bundle

# Update with the location of the reference data files
FASTA_DIR="/data"
FASTA="$FASTA_DIR/human_g1k_v37_decoy.fasta"
KNOWN_DBSNP="$FASTA_DIR/dbsnp_138.b37.vcf.gz"
KNOWN_INDELS="$FASTA_DIR/1000G_phase1.indels.b37.vcf.gz"
KNOWN_MILLS="$FASTA_DIR/Mills_and_1000G_gold_standard.indels.b37.vcf.gz"

cd /data
# Other settings
PCRFREE=true # The data was sequenced with a PCR-free library prep
NT=$(nproc) #number of threads to use in computation, set to number of cores in the server
START_DIR="$PWD" #Determine where the output files will be stored

# ************************************************************************************************************************************************************************

# ******************************************
# 0. Setup
# ******************************************
WORKDIR="$START_DIR/${SM}"
mkdir -p $WORKDIR
cd $WORKDIR

# ******************************************
# 1. Mapping reads with BWA-MEM, sorting
# ******************************************
#The results of this call are dependent on the number of threads used. To have number of threads independent results, add chunk size option -K 10000000

echo $(date) "start migrating"
#/opt/memverge/bin/float migrate --sync -f -c 36:36 -m 72 -j ${FLOAT_JOB_ID}
#/opt/memverge/bin/float migrate --sync -f -t c5a.8xlarge -j ${FLOAT_JOB_ID}
/opt/memverge/bin/float migrate --sync -f -t c5a.8xlarge -j ${FLOAT_JOB_ID}

NT=$(nproc)
echo $(date) "[Task] BWA-MEM, sorting"
start=`date +%s`
( sentieon bwa mem -R "@RG\tID:$RGID\tSM:$SM\tPL:$PL" \
    -t $NT -K 5000000 -x $DNASCOPE_MODEL/bwa.model \
    $FASTA $FASTQ_1 $FASTQ_2 || { echo -n 'BWA error'; exit 1; } ) | \
    sentieon util sort -r $FASTA -o sorted.bam -t $NT \
    --sam2bam -i - || { echo "Alignment failed"; exit 1; }

# Move output files to RESULT directory
cp sorted.bam sorted.bam.bai "$RESULT/"

end=`date +%s`
echo $(date) "[Task End] time used: " $((end-start)) "sec"

# ******************************************
# 2. Metrics
# ******************************************
echo $(date) "start migrating"
#/opt/memverge/bin/float migrate --sync -f -c 48:48 -m 96 -j ${FLOAT_JOB_ID}
#/opt/memverge/bin/float migrate --sync -f -t c5a.12xlarge -j ${FLOAT_JOB_ID}
/opt/memverge/bin/float migrate --sync -f -t c5a.12xlarge -j ${FLOAT_JOB_ID}
NT=$(nproc)
echo $(date) "[Task] Metrics"
start=`date +%s`
sentieon driver -r $FASTA -t $NT -i sorted.bam \
    --algo MeanQualityByCycle mq_metrics.txt --algo QualDistribution qd_metrics.txt \
    --algo GCBias --summary gc_summary.txt gc_metrics.txt --algo AlignmentStat \
    --adapter_seq '' aln_metrics.txt --algo InsertSizeMetricAlgo is_metrics.txt || \
    { echo "Metrics failed"; exit 1; }

sentieon plot GCBias -o gc-report.pdf gc_metrics.txt
sentieon plot QualDistribution -o qd-report.pdf qd_metrics.txt
sentieon plot MeanQualityByCycle -o mq-report.pdf mq_metrics.txt
sentieon plot InsertSizeMetricAlgo -o is-report.pdf is_metrics.txt

# Move output files to RESULT directory
cp mq_metrics.txt qd_metrics.txt gc_summary.txt gc_metrics.txt aln_metrics.txt is_metrics.txt gc-report.pdf qd-report.pdf mq-report.pdf is-report.pdf "$RESULT/"


end=`date +%s`
echo $(date) "[Task End] time used: " $((end-start)) "sec"

# ******************************************
# 3. Remove Duplicate Reads. It is possible
# to remove instead of mark duplicates
# by adding the --rmdup option in Dedup
# ******************************************
echo $(date) "[Task] dedup"
start=`date +%s`
sentieon driver -t $NT -i sorted.bam --algo LocusCollector \
    --fun score_info score.txt || { echo "LocusCollector failed"; exit 1; }

sentieon driver -t $NT -i sorted.bam --algo Dedup \
    --score_info score.txt --metrics dedup_metrics.txt deduped.bam || \
    { echo "Dedup failed"; exit 1; }

# Move output files to RESULT directory
cp dedup_metrics.txt deduped.bam "$RESULT/"

end=`date +%s`
echo $(date) "[Task End] time used: " $((end-start)) "sec"

# ******************************************
# 4a. DNAscope variant calling
# ******************************************
echo $(date) "[Task] DNAscope variant calling"
start=`date +%s`

indel_model_arg=""
if [ "$PCRFREE" = true ]; then
    indel_model_arg="--pcr_indel_model none"
fi
sentieon driver -r $FASTA -t $NT -i deduped.bam \
    --algo DNAscope $indel_model_arg --model $DNASCOPE_MODEL/dnascope.model \
    -d $KNOWN_DBSNP output-ds_tmp.vcf.gz || \
    { echo "DNAscope failed"; exit 1; }
# Move output files to RESULT directory
cp output-ds_tmp.vcf.gz "$RESULT/"

end=`date +%s`
echo $(date) "[Task End] time used: " $((end-start)) "sec"

# ******************************************
# 4b. Variant filtering and genotyping
# ******************************************
echo $(date) "[Task] Variant filtering and genotyping"
start=`date +%s`

sentieon driver -r $FASTA -t $NT \
    --algo DNAModelApply --model $DNASCOPE_MODEL/dnascope.model \
    -v output-ds_tmp.vcf.gz output-ds.vcf.gz || \
    { echo "DNAModelApply failed"; exit 1; }

# Move output files to RESULT directory
cp output-ds.vcf.gz "$RESULT/"

end=`date +%s`
echo $(date) "[Task End] time used: " $((end-start)) "sec"

Input R1 fastq file in .gz format in your s3 bucket. 
For example, input/HG002.novaseq.pcr-free.32x.R1.fastq.gz

Input R2 fastq file in .gz format in your s3 bucket. 
For example, input/HG002.novaseq.pcr-free.32x.R2.fastq.gz