#!/bin/bash
##################################################################################
# Andy Rampersaud, 02.23.16
#This script would be used to summarize TopHat2 mapping statistics 
#Way to run script:
#Usage: 
#./TopHat_Paired_End_Summary.sh
##################################################################################
#---------------------------------------------------------------------------------
#Source the setup file to initialize variables
source ../00_Setup_Pipeline/01_Pipeline_Setup.sh
#---------------------------------------------------------------------------------
#Check that each variable prints a value to the terminal:
echo "-----------------------"
echo "Start of variable list:"
echo "-----------------------"
echo "Dataset_DIR:"
echo ${Dataset_DIR}
echo "Sample_Labels_DIR:"
echo ${Sample_Labels_DIR}
echo "Bowtie2Index_DIR:"
echo ${Bowtie2Index_DIR}
echo "GTF_Files_DIR:"
echo ${GTF_Files_DIR}
echo "ANNOTATION_FILE:"
echo ${ANNOTATION_FILE}
echo "STRANDEDNESS:"
echo ${STRANDEDNESS}
echo "TIME_LIMIT:"
echo ${TIME_LIMIT}
echo "-----------------------"
echo "End of variable list"
echo "-----------------------"
#---------------------------------------------------------------------------------
##################################################################################
#---------------------------------------------------------------------------------
#Retrieve the job name for this step:
Current_DIR=$(pwd)
#Extract the folder name:
DIR_name=`basename ${Current_DIR}`
#---------------------------------------------------------------------------------
OUTPUT_DIR=${Dataset_DIR}/Scripts/${DIR_name}/Job_Summary
###############################
if [ ! -d ${OUTPUT_DIR} ]; then
mkdir -p ${OUTPUT_DIR}
else
#Remove dir:
rm -r ${OUTPUT_DIR}
#Make new dir:
mkdir -p ${OUTPUT_DIR}
fi
###############################
OUTPUT_FILE=$OUTPUT_DIR/TopHat2_Stats.txt
######################
if [ -f $OUTPUT_FILE ]
then 
rm $OUTPUT_FILE
else
touch $OUTPUT_FILE
fi
######################
OUTPUT_FILE_BestMapped=$OUTPUT_DIR/TopHat2_Stats_BestMapped.txt
######################
if [ -f $OUTPUT_FILE_BestMapped ]
then 
rm $OUTPUT_FILE_BestMapped
else
touch $OUTPUT_FILE_BestMapped
fi
######################
OUTPUT_FILE_UniqueReads=$OUTPUT_DIR/TopHat2_Stats_UniqueReads.txt
######################
if [ -f $OUTPUT_FILE_UniqueReads ]
then 
rm $OUTPUT_FILE_UniqueReads
else
touch $OUTPUT_FILE_UniqueReads
fi
######################
OUTPUT_FILE_SpliceReads=$OUTPUT_DIR/TopHat2_Stats_SpliceReads.txt
######################
if [ -f $OUTPUT_FILE_SpliceReads ]
then 
rm $OUTPUT_FILE_SpliceReads
else
touch $OUTPUT_FILE_SpliceReads
fi
######################
cd ${Current_DIR}
#-----------------------------------------------
#Minor issue with jobs using multiple cores
#Empty job output log files (*.pe* and *.po*) are created
#Remove them if they exist:
echo 'Removing  *.pe* files...'
#Need the 2>/dev/null to supress the "No such file or directory"
count=`ls -1 *.pe* 2>/dev/null  | wc -l`
if [ $count != 0 ]
then 
rm *.pe*
fi 
echo 'Removing  *.po* files...'
#Need the 2>/dev/null to supress the "No such file or directory"
count=`ls -1 *.po* 2>/dev/null  | wc -l`
if [ $count != 0 ]
then 
rm *.po*
fi 
#-----------------------------------------------
################################################
#A text file (Sample_Labels.txt) is needed to run this script
SCRIPT_DIR=$(pwd)
cp $Sample_Labels_DIR/Sample_Labels.txt $SCRIPT_DIR
################################################
#The text file is formatted like the following:
#----------------------------------------------
#Sample_ID	Description	
#G110_M1	STAT5_High	
#G110_M2	STAT5_High	
#----------------------------------------------
#The 1st column: Waxman Lab Sample_ID 
#The 2nd column: Sample's description 
################################################
#Print header to output file:
echo 'Read_Mate' $'\t''SAMPLE_ID' $'\t''DESCRIPTION' $'\t''TOTAL_SEQUENCED_READS' $'\t''MULTIMAPPED_READS_ONLY' $'\t''MAPPED_READS_(BEST_LOCATION_MAPPED)' $'\t''BEST_LOCATION_MAPPED_RATIO' >> $OUTPUT_FILE
################################################
#Print header to output file:
echo 'SAMPLE_ID' $'\t''DESCRIPTION' $'\t''Aligned_Pairs' $'\t''Concordant_pair_alignment_rate' $'\t''TOTAL_SEQUENCED_READS' $'\t''MAPPED_READS_(BEST_LOCATION_MAPPED)' $'\t''MAPPED_READS_(BEST_LOCATION_MAPPED)_RATIO' >> $OUTPUT_FILE_BestMapped
################################################
#Print header to output file:
echo 'SAMPLE_ID' $'\t''DESCRIPTION' $'\t''Aligned_Pairs' $'\t''Concordant_pair_alignment_rate' $'\t''TOTAL_SEQUENCED_READS' $'\t''uniquely_mapped_reads' $'\t''uniquely_mapped_reads_RATIO' >> $OUTPUT_FILE_UniqueReads
################################################
#Print header to output file:
echo 'SAMPLE_ID' $'\t''DESCRIPTION' $'\t''Aligned_Pairs' $'\t''Concordant_pair_alignment_rate' $'\t''TOTAL_SEQUENCED_READS' $'\t''splice_junction_READS' $'\t''SPLICED_READ_RATIO' >> $OUTPUT_FILE_SpliceReads
################################################
#Text file has a header line to ignore:
tail -n +2 Sample_Labels.txt > Sample_Labels.temp
#Use a while loop to run jobs
while IFS=$'\t' read -r -a myArray
do
#---------------------------
##Check that text file is read in properly:
#echo 'Sample_DIR:'
Sample_DIR=${myArray[0]}
#echo 'Sample_ID:'
Sample_ID=${myArray[1]}
#echo $Sample_ID
#echo 'Description:'
Description=${myArray[2]}
#echo $Description
#---------------------------
echo $Sample_ID
#Need to cd to sample specific TopHat2 folder
cd ${Dataset_DIR}/$Sample_ID/fastq/tophat2
#---------------------------------------------------------------------------------
#Summary for $OUTPUT_FILE
#---------------------------------------------------------------------------------
#Get counts from *_align_summary.txt file
TOTAL_SEQUENCED_READS_Left=$(grep 'Input' $Sample_ID'_align_summary.txt'  | awk '{print $3}' | awk 'NR==1 {print $1}')
TOTAL_SEQUENCED_READS_Right=$(grep 'Input' $Sample_ID'_align_summary.txt'  | awk '{print $3}' | awk 'NR==2 {print $1}')
#Get counts from *_align_summary.txt file
MULTIMAPPED_READS_ONLY_Left=$(grep 'of these:' $Sample_ID'_align_summary.txt'  | awk '{print $3}' | awk 'NR==1 {print $1}')
MULTIMAPPED_READS_ONLY_Right=$(grep 'of these:' $Sample_ID'_align_summary.txt'  | awk '{print $3}' | awk 'NR==2 {print $1}')
#Get counts from *_align_summary.txt file
MAPPED_READS_BEST_LOCATION_MAPPED_Left=$(grep 'Mapped'  $Sample_ID'_align_summary.txt'  | awk '{print $3}' | awk 'NR==1 {print $1}')
MAPPED_READS_BEST_LOCATION_MAPPED_Right=$(grep 'Mapped'  $Sample_ID'_align_summary.txt'  | awk '{print $3}' | awk 'NR==2 {print $1}')
#Calculate percentages
BEST_LOCATION_MAPPED_RATIO_Left=$(echo "scale=4;$MAPPED_READS_BEST_LOCATION_MAPPED_Left/$TOTAL_SEQUENCED_READS_Left" | bc)
BEST_LOCATION_MAPPED_RATIO_Right=$(echo "scale=4;$MAPPED_READS_BEST_LOCATION_MAPPED_Right/$TOTAL_SEQUENCED_READS_Right" | bc)
#Print to output file
echo 'Left reads:' $'\t'$Sample_ID $'\t'$Description $'\t'$TOTAL_SEQUENCED_READS_Left $'\t'$MULTIMAPPED_READS_ONLY_Left $'\t'$MAPPED_READS_BEST_LOCATION_MAPPED_Left $'\t'$BEST_LOCATION_MAPPED_RATIO_Left >> $OUTPUT_FILE
echo 'Right reads:' $'\t'$Sample_ID $'\t'$Description $'\t'$TOTAL_SEQUENCED_READS_Right $'\t'$MULTIMAPPED_READS_ONLY_Right $'\t'$MAPPED_READS_BEST_LOCATION_MAPPED_Right $'\t'$BEST_LOCATION_MAPPED_RATIO_Right >> $OUTPUT_FILE
#---------------------------------------------------------------------------------
#Summary for $OUTPUT_FILE_BestMapped
#---------------------------------------------------------------------------------
#----------------------------------------------
#Need to get the number of Aligned pairs:
Aligned_Pairs=$(grep 'Aligned pairs:'  $Sample_ID'_align_summary.txt'  | awk '{print $3}')
#Get the "concordant pair alignment rate"
Align_Rate=$(grep 'concordant pair alignment rate.'  $Sample_ID'_align_summary.txt'  | awk '{print $1}')
#----------------------------------------------
#Need the TOTAL_SEQUENCED_READS = (left reads + right reads)
TOTAL_SEQUENCED_READS=$(echo "scale=4;$TOTAL_SEQUENCED_READS_Left + $TOTAL_SEQUENCED_READS_Right" | bc)
#Need the *_statistics_for_primary_reads.txt
BEST_Mapped_Reads=$(grep 'total' $Sample_ID'_statistics_for_primary_reads.txt'  | awk '{print $1}')
#Calculate percentages
BEST_Mapped_Reads_RATIO=$(echo "scale=4;$BEST_Mapped_Reads/$TOTAL_SEQUENCED_READS" | bc)
#Print to output file
echo $Sample_ID $'\t'$Description $'\t'$Aligned_Pairs $'\t'$Align_Rate $'\t'$TOTAL_SEQUENCED_READS $'\t'$BEST_Mapped_Reads $'\t'$BEST_Mapped_Reads_RATIO >> $OUTPUT_FILE_BestMapped
#---------------------------------------------------------------------------------
#Summary for $OUTPUT_FILE_UniqueReads
#---------------------------------------------------------------------------------
#Need to summarize the uniquely mapped reads:
uniquely_mapped_reads=$(grep 'total' $Sample_ID'_statistics_for_primary_unique_reads.txt'  | awk '{print $1}')
#Calculate percentage:
uniquely_mapped_reads_RATIO=$(echo "scale=4;$uniquely_mapped_reads/$TOTAL_SEQUENCED_READS" | bc)
#Print to output file
echo $Sample_ID $'\t'$Description $'\t'$Aligned_Pairs $'\t'$Align_Rate $'\t'$TOTAL_SEQUENCED_READS  $'\t'$uniquely_mapped_reads $'\t'$uniquely_mapped_reads_RATIO >> $OUTPUT_FILE_UniqueReads
#---------------------------------------------------------------------------------
#Summary for $OUTPUT_FILE_SpliceReads
#---------------------------------------------------------------------------------
#Print 2nd line of the *_spliced_read_counts.txt:
#Reminder: the following counts are from the single BAM file
#SKIPPED_READS=$(awk 'NR == 2' $Sample_ID'_spliced_read_counts.txt' | awk '{print $2}')
SPLICED_READS=$(awk 'NR == 2' $Sample_ID'_spliced_read_counts.txt' | awk '{print $3}')
#----------------------------------------------
#Calculate percentage:
SPLICED_READS_RATIO=$(echo "scale=4;$SPLICED_READS/$TOTAL_SEQUENCED_READS" | bc)
#Print to output file
echo $Sample_ID $'\t'$Description $'\t'$Aligned_Pairs $'\t'$Align_Rate $'\t'$TOTAL_SEQUENCED_READS  $'\t'$SPLICED_READS $'\t'$SPLICED_READS_RATIO >> $OUTPUT_FILE_SpliceReads
#----------------------------------------------
done < Sample_Labels.temp
cd ${Current_DIR}
#Remove the temp file:
rm Sample_Labels.temp
echo '#--------------------------------------------------------------------------'
echo 'Check out '${OUTPUT_DIR}
#echo 'Check out '$OUTPUT_FILE
#echo 'Check out '$OUTPUT_FILE_BestMapped
#echo 'Check out '$OUTPUT_FILE_UniqueReads
#echo 'Check out '$OUTPUT_FILE_SpliceReads
echo '#--------------------------------------------------------------------------'
##################################################################################