#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

Sentieon DNAscope Germline SNV/INDEL Variant Calling

Highlight

Detail