#!/bin/bash
# 02 November 2009, G. Robertson

#=============================================================================
BWADIR=/home/pubseq/BioSw/bwa/bwa-0.5.4
#=============================================================================
#CHROMINFO=/projects/mapp_etc/Mappability/fastq/mm9/mm9.chromInfo.txt
CHROMINFO=/home/grobertson/bin/x86_64/mm9.chrom.sizes
#CHROMINFO=/home/grobertson/bin/x86_64/mm8.chrom.sizes
#GENOMEFA=/projects/mapp_etc/Mappability/fastq/mm9/36mer/bwa_ind/mm9_build37_mouse.fasta
# 'chr1' is '1'
GENOMEFA=/projects/wtsspipeline/resources/Mus_musculus/bfa_NCBI-37/all_mouse_no_chr.fasta
#=============================================================================



# MM0570
BAMDIR=/projects/analysis/analysis5/MM0570/61U0GAAXX_5/bwa
BAMFILE=61U0GAAXX_5_withJunctionsOnGenome_dupsFlagged.bam
FILENAME=61U0GAAXX_5_withJunctionsOnGenome_dupsFlagged
OUTDIR=/projects/remc_bigdata/MORGEN/MM0570

#=============================================================================
#BAMFILE=$FILENAME.bam
SORTEDBAMFILE=$FILENAME.sorted.bam
SORTEDBAIFILE=$FILENAME.sorted.bam.bai
PILEUPFILE=$FILENAME.Q10_m512.pileup
COLSPFILE=$PILEUPFILE.3_cols
FXSTEPFILE=$COLSPFILE.fixedStep.wigdata
FXSTEPCHRFILE=$COLSPFILE.fixedStep.chr.wigdata
BWFILE=$PILEUPFILE.bw
#=============================================================================

echo; echo "Sorting..."
echo $BAMFILE
echo $SORTEDBAMFILE
echo "samtools sort $BAMDIR/$BAMFILE sorted && mv sorted.bam $OUTDIR/$SORTEDBAMFILE"
samtools sort $BAMDIR/$BAMFILE sorted && mv sorted.bam $OUTDIR/$SORTEDBAMFILE

echo "done"; echo "Indexing..."
echo "samtools index $OUTDIR/$SORTEDBAMFILE"
samtools index $OUTDIR/$SORTEDBAMFILE

# Nina: pileup by default exludes duplicates and chastity-flagged reads
# -m INT    filtering reads with bits in INT [1796]
# to retain duplicates but exclude chastity-flagged reads, use: -m 512
# SNP-calling does not use chastity or duplicates. i.e. use default pileup settings
# WTSS coverage, exclude chastity but retain duplicates. use: -m 512

# MAPQ filter:
# -q INT   minimum mapping quality [0]
echo "done"; echo "Pileup...require MAPQ>10, exclude chastity-flagged reads"
echo "samtools view -q 10 $OUTDIR/$SORTEDBAMFILE | samtools pileup -S -m 512 -cf $GENOMEFA - > $OUTDIR/$PILEUPFILE"
samtools view -q 10 $OUTDIR/$SORTEDBAMFILE | samtools pileup -S -m 512 -cf $GENOMEFA - > $OUTDIR/$PILEUPFILE

echo "done"; echo "Pileup..."
echo "samtools pileup -m 512 -cf $GENOMEFA $OUTDIR/$SORTEDBAMFILE > $OUTDIR/$PILEUPFILE"
samtools pileup -m 512 -cf $GENOMEFA $OUTDIR/$SORTEDBAMFILE > $OUTDIR/$PILEUPFILE

echo "done"; echo "Write a minimal set of 3 pileup columns"
echo "awk '{ print \$1 \"\t\" \$2 \"\t\" \$8 > \"$OUTDIR/$COLSPFILE\" }' $OUTDIR/$PILEUPFILE"
awk '{ print $1 "\t" $2 "\t" $8 > cf }' cf=$OUTDIR/$COLSPFILE $OUTDIR/$PILEUPFILE

echo "done"; echo "Transformimg pileup into fixedStep WIG data, gzipping..."
less $OUTDIR/$COLSPFILE | perl -pe '($c, $start, $depth) = split; if ($c ne $lastC || $start != $lastStart+1) { print "fixedStep chrom=$c start=$start step=1 span=1\n";  }  $_ = "$depth\n"; ($lastC, $lastStart) = ($c, $start);' | gzip -c > $OUTDIR/$FXSTEPFILE.gz

#echo "change to 'chr1' from '1', if you need to"
#less $OUTDIR/$FXSTEPFILE.gz | sed 's/chrom=/chrom=chr/' > $OUTDIR/$FXSTEPCHRFILE
#gzip $OUTDIR/$FXSTEPCHRFILE

echo "done"; echo "running wigToBigWig... (~/bin/x86_64/wigToBigWig)"
echo "wigToBigWig -clip $OUTDIR/$FXSTEPCHRFILE.gz $CHROMINFO $OUTDIR/$BWFILE"
wigToBigWig -clip $OUTDIR/$FXSTEPCHRFILE.gz $CHROMINFO $OUTDIR/$BWFILE

echo; echo "done!"