#!/bin/bash -l
##################################################################################
#Andy Rampersaud, 08.07.17
#This script is called by setup_UCSC_BigWig.sh
##################################################################################
# Specify which shell to use
#$ -S /bin/bash
# Run on the current working directory
#$ -cwd

# Join standard output and error to a single file  
#$ -j n
# change to y if you want a single qlog file 
##################################################################################
#Initialize variables from UCSC_BigWig.sh
##################################################################################
#checking the command line arg
#-ne : "is not equal to"
if [ $# -ne 7 ] ; then
      echo "Need 7 arguments for the qsub command:"
      echo "qsub -N ${Job_Name}'_'${Sample_ID} -P waxmanlab -l h_rt=${TIME_LIMIT} UCSC_BigWig.qsub ${Sample_ID} ${Dataset_DIR} ${Sample_Labels_DIR} ${STRAND_RULE} ${BU_User} ${VM_DIR_UCSC} ${SCRIPT_DIR}"
      exit 0
fi

#http://www.ibm.com/developerworks/library/l-bash-parameters/
#Note: If you have more than 9 parameters, you cannot use $10 to refer to the tenth one. You must first either process or save the first parameter ($1), then use the shift command to drop parameter 1 and move all remaining parameters down 1, so that $10 becomes $9 and so on.

#http://unix.stackexchange.com/questions/104420/how-to-use-command-line-arguments-in-a-shell-script
#If you need access more than 9 command line arguments, you can use the shift command. Example: shift 2 renames $3 to $1, $4 to $2 etc.

#process the command line arguments
Sample_ID=$1
Dataset_DIR=$2
Sample_Labels_DIR=$3
STRAND_RULE=$4
BU_User=$5
VM_DIR_UCSC=$6
SCRIPT_DIR=$7

#Print variables (make sure they appear correctly):
echo "-----------------------"
echo "Start of variable list:"
echo "-----------------------"
echo "Sample_ID:"
echo ${Sample_ID}
echo "Dataset_DIR:"
echo ${Dataset_DIR}
echo "Sample_Labels_DIR:"
echo ${Sample_Labels_DIR}
echo "STRAND_RULE:"
echo ${STRAND_RULE}
echo "BU_User:"
echo ${BU_User}
echo "VM_DIR_UCSC:"
echo ${VM_DIR_UCSC}
echo "SCRIPT_DIR:"
echo ${SCRIPT_DIR}
echo "-----------------------"
echo "End of variable list"
echo "-----------------------"

# Now let's keep track of some information just in case anything goes wrong
echo "=========================================================="
#Use to calculate job time:
#Start_Time in seconds
Start_Time=$(date +"%s")
echo "Starting on : $(date)"
echo "Running on node : $(hostname)"
echo "Current directory : $(pwd)"
echo "Current job ID : $JOB_ID"
echo "Current job name : $JOB_NAME"
echo "Task index number : $SGE_TASK_ID"
echo "Parameter for multiple cores : $NSLOTS"
echo "=========================================================="

# Go to local scratch directory
echo
echo 'Change dir to scratch directory'
echo
cd ${TMPDIR}
echo
echo 'Print scratch directory location:'
echo
echo $TMPDIR
#--------------------------------------
echo
echo 'Loading required modules...'
echo
#Make sure the shebang line = #!/bin/bash -l
#Need the -l option to load modules
#Search for latest program installed:
#module avail -t 2>&1 | grep -i samtools
module load samtools/samtools-0.1.19_gnu446
#module avail -t 2>&1 | grep -i python
module load python2.7
#module avail -t 2>&1 | grep -i R-2.15.1
module load R/R-2.15.1_gnu-4.4.6   
#module avail -t 2>&1 | grep -i RSeQC                
module load rseqc/2.4
#--------------------------------------
#module help samtools/samtools-0.1.19_gnu446
#---------------------------------------------------------------------------------
#----------- Module Specific Help for 'samtools/samtools-0.1.19_gnu446' ---------------------------
#SAMtools 0.1.19_gnu446: Suite of programs for interacting with high-throughput sequencing data.

#SAM (Sequence Alignment/Map) format is a generic format for storing
#large nucleotide sequence alignments. SAM Tools provide various utilities
#for manipulating alignments in the SAM format, including sorting, merging,
#indexing and generating alignments in a per-position format.
#---------------------------------------------------------------------------------
#module help python2.7
#----------- Module Specific Help for 'python2.7/Python-2.7.3_gnu446' ---------------------------
#python2.7 Python-2.7.3_gnu446: A general purpose, interpretive programming language.

#For more information on python2.7, please see http://www.python.org
#---------------------------------------------------------------------------------
#module help R/R-2.15.1_gnu-4.4.6   
#----------- Module Specific Help for 'R/R-2.15.1_gnu-4.4.6' -------
#R 2.15.1  is a system (language) for statistical computation and graphics

#<<R consists of a language plus run-time enviroment with graphics, a debugger, access to certain system functions and the ability to run programs stored in script files>>

#For more information on R, visit http://cran.r-project.org/
#---------------------------------------------------------------------------------
#module help rseqc/2.4
#----------- Module Specific Help for 'rseqc/2.4' ------------------
#rseqc 2.4 An RNA-seq Quality Control Package
#RSeQC package provides a number of useful modules that can comprehensively evaluate high throughput sequence data especially RNA-seq data. Some basic modules quickly inspect sequence quality, nucleotide composition bias, PCR bias and GC bias, while RNA-seq specific modules evaluate sequencing saturation, mapped reads distribution, coverage uniformity, strand specificity, etc.
#For more information on rseqc, please see http://rseqc.sourceforge.net/
#---------------------------------------------------------------------------------
# copy user input data files to scratch
cp ${Dataset_DIR}/${Sample_ID}/fastq/tophat2/${Sample_ID}'_primary_unique.bam' .
#Initialize INPUT_BAM:
INPUT_BAM=${Sample_ID}'_primary_unique.bam'

echo
echo 'INPUT_BAM:'
echo ${INPUT_BAM}
echo

#Copy UCSC_BigWig required files:
cp ${SCRIPT_DIR}/Chrom_Sizes/mm9.chrom.sizes .
#Initialize variable:
CHROMSIZE_FILE=mm9.chrom.sizes
#cp ${SCRIPT_DIR}/UCSC_Tools/wigToBigWig .
cp -r ${SCRIPT_DIR}/UCSC_Tools .

STORAGE_DIR=${Dataset_DIR}/${Sample_ID}/fastq/tophat2
#Create UCSC_BigWig output folder to store files:
OUTPUT_DIR=${TMPDIR}/'UCSC_BigWig'
###############################
if [ ! -d $OUTPUT_DIR ]; then
mkdir $OUTPUT_DIR
fi
###############################

echo
echo 'List files in scratch directory:'
echo
ls -alh

echo
echo 'Starting to run my commands'
echo

echo
echo 'Starting samtools sort and index'
echo

#samtools sort only takes file name without extension
samtools sort ${INPUT_BAM} ${Sample_ID}'_sorted'
#samtools sort only takes file name with extension
samtools index ${Sample_ID}'_sorted'.bam
#Initilize variable:
INPUT_FILE=${Sample_ID}'_sorted.bam'

echo
echo 'Finished samtools sort and index'
echo

echo
echo 'Starting bam2wig.py command'
echo
#--------------------------------------------------------------------------------
#http://rseqc.sourceforge.net/#infer-experiment-py
#1. bam2wig.py converts all types of RNA-seq data from BAM format into wiggle format.
#2. If UCSC wigToBigWig tool was found, output wiggle file will be converted into bigwig format automatically.
#3. To use bam2wig.py, BAM file must be sorted and indexed properly using SAMtools. Below example shows how to sort and index BAM file using samTools
#--------------------------------------------------------------------------------
#bam2wig.py --help
#Usage: bam2wig.py [options]

#Options:
#  --version             show program's version number and exit
#  -h, --help            show this help message and exit
#  -i INPUT_FILE, --input-file=INPUT_FILE
#                        Alignment file in BAM format. BAM file must be sorted
#                        and indexed using samTools. .bam and .bai files should
#                        be placed in the same directory. HowTo:
#                        http://genome.ucsc.edu/goldenPath/help/bam.html
#  -s CHROMSIZE, --chromSize=CHROMSIZE
#                        Chromosome size file. Tab or space separated text file
#                        with 2 columns: first column is chromosome name,
#                        second column is size of the chromosome. Chromosome
#                        name (such as "chr1") should be consistent between
#                        this file and BAM file.
#  -o OUTPUT_PREFIX, --out-prefix=OUTPUT_PREFIX
#                        Prefix of output wiggle files(s). One wiggle file will
#                        be generated for non-strand specific data, two wiggle
#                        files ("Prefix_Forward.wig" and "Prefix_Reverse.wig")
#                        will be generated for strand specific RNA-seq data.
#  -t TOTAL_WIGSUM, --wigsum=TOTAL_WIGSUM
#                        Specified wigsum. Eg: 1,000,000,000 equals to coverage
#                        of 10 million 100nt reads. Ignore this option to
#                        disable normalization
#  -u, --skip-multi-hits
#                        Presence this option render the program to skip non-
#                        unique hit reads.
#  -d STRAND_RULE, --strand=STRAND_RULE
#                        How read(s) were stranded during sequencing. For
#                        example: --strand='1++,1--,2+-,2-+' means that this is
#                        a pair-end, strand-specific RNA-seq, and the strand
#                        rule is: read1 mapped to '+' => parental gene on '+';
#                        read1 mapped to '-' => parental gene on '-'; read2
#                        mapped to '+' => parental gene on '-'; read2 mapped to
#                        '-' => parental gene on '+'.  If you are not sure
#                        about the strand rule, run 'infer_experiment.py'
#                        default=none (Not a strand specific RNA-seq data).
#  -q MAP_QUAL, --mapq=MAP_QUAL
#                        Minimum mapping quality for an alignment to be called
#                        "uniquely mapped". default=30
#--------------------------------------------------------------------------------
#Notes about command options:
TOTAL_WIGSUM=1000000000
#The number you want your wigfiles to be normalized to. Please put 1000000000. For the normalization, all wig files are normalized to a random, fixed number that I chose (1000000000). So, all wig files that are created with this script are comparable. The normalization is basically enforcing that the area under the wig-file curves are constant for all samples (i.e. the "drawing area" for each sample are fixed, it depends on each sample where to distribute its drawing area).
#--------------------------------------------------------------------------------
#For some reason "none" is not recognized by bam2wig.py
#Need an if statement to check if ${STRAND_RULE} equals "none"
#------------------------------------------
if [ "${STRAND_RULE}" == "none" ];
then
#If the data is unstranded, then don't specify --strand option (default=none)
bam2wig_options="--input-file=$INPUT_FILE --chromSize=$CHROMSIZE_FILE --out-prefix=${Sample_ID} --wigsum=${TOTAL_WIGSUM}"
else
bam2wig_options="--input-file=$INPUT_FILE --chromSize=$CHROMSIZE_FILE --out-prefix=${Sample_ID} --wigsum=${TOTAL_WIGSUM} --strand=${STRAND_RULE}"
fi
#------------------------------------------
echo 'Printing command:'
echo "bam2wig.py "${bam2wig_options}
echo
#Run command:
bam2wig.py ${bam2wig_options}
#--------------------------------------------------------------------------------

echo
echo 'Finished bam2wig.py command'
echo

echo
echo 'Sample from wiggle files:'
echo

head -10 *.wig

#Need some processing of the wig files before wigToBigWig conversion
#Previously saw error with "-" in the 2nd column of wig file (cause an error with wigToBigWig conversion)
#Move the wigToBigWig to another folder (copy the whole UCSC_Tools/)
#Need to by-pass the bam2wig.py finding the wigToBigWig script

wig_list=*.wig
for wig_file in $wig_list
do
wig_file_name=${wig_file%\.wig}
echo
echo '#--------------------------------'
echo $wig_file
echo '#--------------------------------'
echo 'Line count where 2nd column does not equal "-":'
awk '$2!="-" {print $0}' $wig_file | wc -l 
echo 'Line count where 2nd column is only a negative sign:'
awk '$2=="-" {print $0}' $wig_file | wc -l 
echo 'Copy '$wig_file' to UCSC_Tools folder'
cp $wig_file ./UCSC_Tools
echo 'Copy mm9.chrom.sizes to UCSC_Tools folder'
cp mm9.chrom.sizes ./UCSC_Tools
echo 'Start wigToBigWig'
cd UCSC_Tools
#Change permissions to run program:
chmod 700 wigToBigWig
./wigToBigWig $wig_file mm9.chrom.sizes $OUTPUT_DIR/$wig_file_name.bw
echo 'End wigToBigWig'
cd ..
done

echo
echo 'List of files in UCSC_Tools folder:'
echo
ls -al ./UCSC_Tools/*

#As mentioned above:
#2. If UCSC wigToBigWig tool was found, output wiggle file will be converted into bigwig format automatically.
#For some reason when bam2wig.py calls the wigToBigWig tool it may return an error
#Better to call the wigToBigWig tool directly (as I do above)

#--------------------------------------------------------------------------------
#./wigToBigWig 
#wigToBigWig v 4 - Convert ascii format wig file (in fixedStep, variableStep
#or bedGraph format) to binary big wig format.
#usage:
#   wigToBigWig in.wig chrom.sizes out.bw
#Where in.wig is in one of the ascii wiggle formats, but not including track lines
#and chrom.sizes is two column: <chromosome name> <size in bases>
#and out.bw is the output indexed big wig file.
#Use the script: fetchChromSizes to obtain the actual chrom.sizes information
#from UCSC, please do not make up a chrom sizes from your own information.
#options:
#   -blockSize=N - Number of items to bundle in r-tree.  Default 256
#   -itemsPerSlot=N - Number of data points bundled at lowest level. Default 1024
#   -clip - If set just issue warning messages rather than dying if wig
#                  file contains items off end of chromosome.
#   -unc - If set, do not use compression.
#--------------------------------------------------------------------------------

#Copy bigWig files to $OUTPUT_DIR
#cp ${Sample_ID}.Forward.bw $OUTPUT_DIR
#cp ${Sample_ID}.Reverse.bw $OUTPUT_DIR

echo
echo 'Copy output to storage dir'
echo
cp -r $OUTPUT_DIR $STORAGE_DIR

echo
echo "List files in scratch"
echo
ls -alh

echo "=========================================================="
echo "Finished on : $(date)"
#Use to calculate job time:
#End_Time in seconds
End_Time=$(date +"%s")
diff=$(($End_Time-$Start_Time))
echo "$(($diff / 3600)) hours, $((($diff / 60) % 60)) minutes and $(($diff % 60)) seconds elapsed."
echo "=========================================================="