#! /bin/sh
#SBATCH --job-name=cram_oncocnv
#SBATCH --partition=debug
#SBATCH --cpus-per-task 2
#SBATCH --output=%x.out
#SBATCH --error=%x.err
#SBATCH --workdir=/projects/tab_archive2/comp
#SBATCH --mem=4G

> /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err

/gsc/software/linux-x86_64-centos6/samtools-1.9/bin/samtools view -@ 20 -T /projects/alignment_references/Homo_sapiens/hg19a/genome/fasta/hg19a.fa -O cram,store_md=1,store_nm=1 -o /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.cram /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam 2> /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err
chmod 664 /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.cram

if [ $? -eq 0 ] && [ ! -s /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err ]
then
    echo CRAM success
else
    echo CRAM failure
    exit $?
fi

/gsc/software/linux-x86_64-centos6/samtools-1.9/bin/samtools index /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.cram 2> /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err
chmod 664 /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.cram.crai

if [ $? -eq 0 ] && [ ! -s /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err ]
then
    echo CRAM index success
else
    echo CRAM index failure
    exit $?
fi


/gsc/software/linux-x86_64-centos6/samtools-1.9/bin/samtools view -@ 20 -T /projects/alignment_references/Homo_sapiens/hg19a/genome/fasta/hg19a.fa --input-fmt-option decode_md=0 /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.cram | cut -f 1,2,3,4,5,6,7,8,9,10,11 | md5sum > /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.cram.non_binary.md5sum 2> /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err
chmod 664 /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.cram.non_binary.md5sum

if [ $? -eq 0 ] && [ ! -s /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err ]
then
    echo CRAM checksum success
else
    echo CRAM checksum failure
    exit $?
fi


/gsc/software/linux-x86_64-centos6/sambamba-0.6.1/sambamba_v0.6.1 view -t 20 /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam  | cut -f 1,2,3,4,5,6,7,8,9,10,11 | md5sum > /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam.non_binary.md5sum 2> /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err
chmod 664 /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam.non_binary.md5sum

if [ $? -eq 0 ] && [ ! -s /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err ]
then
    echo BAM checksum success
else
    echo BAM checksum failure
    exit $?
fi


read -r cram_md5 < /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.cram.non_binary.md5sum
read -r bam_md5 < /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam.non_binary.md5sum

if [ "$cram_md5" = "$bam_md5" ]
then
    echo Content MD5 matches
    #rm -f /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam
    #rm -f /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam.bai
else
    echo Content MD5 mismatch


    /gsc/software/linux-x86_64-centos6/sambamba-0.6.1/sambamba_v0.6.1 view -t 20 -F "unmapped and (cigar != '' or mapping_quality > 0)" /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam | cut -f 1,2,3,4 > /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.nonspec_unmapped.txt 2> /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err
    chmod 664 /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.nonspec_unmapped.txt
    if [ $? -eq 0 ] && [ ! -s /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err ]
    then
    echo Nonspec unmapped reads success
    else
    echo Nonspec unmapped reads failure
    exit $?
    fi

    /gsc/software/linux-x86_64-centos6/sambamba-0.6.1/sambamba_v0.6.1 view -t 20 -F "not (unmapped and (cigar != '' or mapping_quality > 0))" /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam | cut -f 1,2,3,4,5,6,7,8,9,10,11 | md5sum > /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.filtered_reads.bam.non_binary.md5sum 2> /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err
    chmod 664 /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.filtered_reads.bam.non_binary.md5sum
    if [ $? -eq 0 ] && [ ! -s /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err ]
    then
    echo Filtered reads from BAM checksum success
    else
    echo Filtered reads from BAM checksum failure
    exit $?
    fi

    /gsc/software/linux-x86_64-centos6/samtools-1.9/bin/samtools view -@ 20 -T /projects/alignment_references/Homo_sapiens/hg19a/genome/fasta/hg19a.fa --input-fmt-option decode_md=0 /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.cram | grep -vf /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.nonspec_unmapped.txt | cut -f 1,2,3,4,5,6,7,8,9,10,11 | md5sum > /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.filtered_reads.cram.non_binary.md5sum 2> /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err
    chmod 664 /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.filtered_reads.cram.non_binary.md5sum
    if [ $? -eq 0 ] && [ ! -s /projects/tab_archive2/comp/gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/cram_oncocnv.last_command.err ]
    then
    echo Filtered reads CRAM checksum success
    else
    echo Filtered reads CRAM checksum failure
    exit $?
    fi

    read -r filtered_reads_cram_md5 < /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.filtered_reads.cram.non_binary.md5sum
    read -r filtered_reads_bam_md5 < /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.filtered_reads.bam.non_binary.md5sum

    if [ "$filtered_reads_cram_md5" = "$filtered_reads_bam_md5" ]
    then
    echo Filtered content MD5 matches
    #rm -f /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam
    #rm -f /gsc/www/bcgsc.ca/downloads/ybutterf/oncocnv/control.bam.bai
    else
    echo Filtered content MD5 mismatch
    fi

fi