#!/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
#=============================================================================
# MM0472
#BAMDIR=/projects/wtsspipeline/dev/analysis/MM0472/new_pipeline
#FILENAME=MM0472_7_lanes_merged

# http://www.bcgsc.ca/data/sbs/viewer/illumina/libraries/MM0566
# Heart Atrioventricular cushions, E11.5; 36,770,600 reads
#BAMDIR=/projects/analysis/analysis5/MM0566/61FB9AAXX_8/bwa
#BAMFILE=61FB9AAXX_8_withJunctionsOnGenome_dupsFlagged.bam
#FILENAME=61FB9AAXX_8_withJunctionsOnGenome_dupsFlagged

# http://www.bcgsc.ca/data/sbs/viewer/illumina/libraries/MM0564
# Heart Atrioventricular cushions, E12.5; 37,439,456 reads
# /home/aldente/private/Projects/MORGEN/MM0564/AnalyzedData/MM05641./Solexa/Data/current/BaseCalls/BWA_2010-04-25/s_7_paired.dup.sorted.bam
#BAMDIR=/projects/analysis/analysis5/MM0564/61FB9AAXX_7/
#FILENAME=s_7_paired.dup

# 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

# MM0571
#BAMDIR=/projects/analysis/analysis5/MM0571/61U0GAAXX_6/bwa
#BAMFILE=61U0GAAXX_6_withJunctionsOnGenome_dupsFlagged.bam
#FILENAME=61U0GAAXX_6_withJunctionsOnGenome_dupsFlagged
#OUTDIR=/projects/remc_bigdata/MORGEN/MM0571

# MM0581
#BAMDIR=/projects/analysis/analysis5/MM0581/61U0GAAXX_7/bwa
#BAMFILE=61U0GAAXX_7_withJunctionsOnGenome_dupsFlagged.bam
#FILENAME=61U0GAAXX_7_withJunctionsOnGenome_dupsFlagged
# 22 October 2010
BAMDIR=/projects/analysis/analysis5/MM0581/meta_bwa
BAMFILE=MM0581_7_lanes_withJunctionsOnGenome_dupsFlagged.bam
FILENAME=MM0581_7_lanes_withJunctionsOnGenome_dupsFlagged
OUTDIR=/projects/remc_bigdata/MORGEN/MM0581_WTSS_definitive-endo_E8.5

# MM0510
#BAMDIR=/projects/wtsspipeline/dev/analysis/MM0510/new_pipeline
#BAMFILE=MM0510_3153RAAXX_8_merged.bam
#=============================================================================
#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 -f $GENOMEFA - > $OUTDIR/$PILEUPFILE"
#samtools view -q 10 $OUTDIR/$SORTEDBAMFILE | samtools pileup -S -m 512 -f $GENOMEFA - > $OUTDIR/$PILEUPFILE

#echo "done"; echo "Write a minimal set of 3 pileup columns"
#echo "awk '{ print \$1 \"\t\" \$2 \"\t\" \$4 > \"$OUTDIR/$COLSPFILE\" }' $OUTDIR/$PILEUPFILE"
#awk '{ print $1 "\t" $2 "\t" $4 > 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'"
#less $OUTDIR/$FXSTEPFILE.gz | sed 's/chrom=/chrom=chr/' > $OUTDIR/$FXSTEPCHRFILE
#echo "gzipping..."
#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

# HTTPDIR=/gsc/www/bcgsc.ca/downloads/genereg/wtss_assy/
# cp MM0566_61FB9AAXX_8_withJunctionsOnGenome_dupsFlagged.pileup.bw /gsc/www/bcgsc.ca/downloads/genereg/wtss_assy/MM0566_61FB9AAXX_8_withJunctionsOnGenome_dupsFlagged.pileup.bw

#echo; echo "done!"