#!/bin/bash -l
##################################################################################
#Andy Rampersaud, 05.03.16
#Adapted from tophat/paired_end_mapping scripts by Tisha Melia
#This script is called by TopHat_Paired_End.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 TopHat_Paired_End.sh
##################################################################################
#checking the command line arg
#-ne : "is not equal to"
if [ $# -ne 8 ] ; then
      echo "Need 8 arguments for the qsub command:"
      echo "qsub -N ${Job_Name}'_'${Sample_ID} -P waxmanlab -l h_rt=${TIME_LIMIT} TopHat_Paired_End.qsub ${Sample_ID} ${Dataset_DIR} ${Sample_Labels_DIR} ${Bowtie2Index_DIR} ${GTF_Files_DIR} ${ANNOTATION_FILE} ${STRANDEDNESS} ${DISTANCE_BT_READ_PAIR}"
      exit 0
fi
#process the command line arguments
Sample_ID=$1
Dataset_DIR=$2
Sample_Labels_DIR=$3
Bowtie2Index_DIR=$4
GTF_Files_DIR=$5
ANNOTATION_FILE=$6
STRANDEDNESS=$7
DISTANCE_BT_READ_PAIR=$8

#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.

#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 "Bowtie2Index_DIR:"
echo ${Bowtie2Index_DIR}
echo "GTF_Files_DIR:"
echo ${GTF_Files_DIR}
echo "ANNOTATION_FILE:"
echo ${ANNOTATION_FILE}
echo "STRANDEDNESS:"
echo ${STRANDEDNESS}
echo "DISTANCE_BT_READ_PAIR:"
echo ${DISTANCE_BT_READ_PAIR}
echo "-----------------------"
echo "End of variable list"
echo "-----------------------"

#You can use option "-pe omp N", where N is the number of cores ( number of course on our system can be any integer between 1 and 16):
#$ -pe omp 16

# 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 boost
#module load boost/1.54.0
#If you're using tophat/2.1.1 (need boost/1.58.0) 
module load boost/1.58.0
#--------------------------------------
#module avail -t 2>&1 | grep -i samtools
module load samtools/samtools-0.1.19_gnu446
#--------------------------------------
#module avail -t 2>&1 | grep -i bowtie2
module load bowtie2/2.2.2
#--------------------------------------
#module avail -t 2>&1 | grep -i cufflinks
module load cufflinks/2.2.1
#--------------------------------------
#module avail -t 2>&1 | grep -i tophat
#module load tophat/2.0.13
#Newer version:
#tophat/2.1.1
module load tophat/2.1.1
#--------------------------------------
##################################################################################
#module help boost/1.54.0
#----------- Module Specific Help for 'boost/1.54.0' ---------------
#Boost provides free peer-reviewed portable C++ source libraries.
#http://www.boost.org/

#module help samtools/samtools-0.1.19_gnu446
#----------- Module Specific Help for 'samtools/samtools-0.1.19_gnu446' ---------------------------
#sets the environment for samtools (0.1.19) built using GNU Compilers
#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.
#http://samtools.sourceforge.net/

#module help bowtie2/2.2.2
#----------- Module Specific Help for 'bowtie2/2.2.2' --------------
#bowtie2 2.2.2: Ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences.
#Sets the environment for Bowtie2 2.2.2
#Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences. It is particularly good at aligning reads of about 50 up to 100s or 1,000s of characters and particularly good at aligning to relatively long (e.g. mammalian) genomes. Bowtie 2 indexes the genome with an FM Index to keep its memory footprint small: for the human genome, its memory footprint is typically around 3.2 GB. Bowtie 2 supports gapped, local, and paired-end alignment modes.

#module help cufflinks/2.2.1
#----------- Module Specific Help for 'cufflinks/2.2.1' ------------
#cufflinks 2.2.1 Transcript assembly, differential expression, and differential regulation for RNA-Seq.
#Cufflinks assembles transcripts, estimates their abundances, and tests for differential expression and regulation in RNA-Seq samples. It accepts aligned RNA-Seq reads and assembles the alignments into a parsimonious set of transcripts. Cufflinks then estimates the relative abundances of these transcripts based on how many reads support each one, taking into account biases in library preparation protocols.
#For more information on cufflinks, please see https://cole-trapnell-lab.github.io/cufflinks/

#module help tophat/2.0.13
#----------- Module Specific Help for 'tophat/2.0.13' --------------
#tophat 2.0.13 Fast splice junction mapper for RNA-Seq reads.
#TopHat is a fast splice junction mapper for RNA-Seq reads. It aligns RNA-Seq reads to mammalian-sized genomes using the ultra high-throughput short read aligner Bowtie, and then analyzes the mapping results to identify splice junctions between exons.
#For more information on tophat, please see http://ccb.jhu.edu/software/tophat
#--------------------------------------
#module help tophat/2.1.1
#----------- Module Specific Help for 'tophat/2.1.1' ---------------

#tophat 2.1.1 Spliced read mapper for RNA-Seq

#TopHat is a fast splice junction mapper for RNA-Seq reads. It aligns RNA-Seq reads to mammalian-sized genomes using the ultra high-throughput short read aligner Bowtie, and then analyzes the mapping results to identify splice junctions between exons.

#Note: This application requres either bowtie or bowtie2
#For more information on tophat, please see http://ccb.jhu.edu/software/tophat/index.shtml
##################################################################################

# copy user input data files to scratch
#Copy over the *R1*.fastq.gz file(s):
cp ${Dataset_DIR}/${Sample_ID}/fastq/*_1.fastq.gz .
#Copy over the *R2*.fastq.gz file(s):
cp ${Dataset_DIR}/${Sample_ID}/fastq/*_2.fastq.gz .
#Need the Bowtie2 index:
cp -r ${Bowtie2Index_DIR}/Bowtie2Index .
#Copy RefSeq_GeneBody.gtf
cp ${GTF_Files_DIR}/RefSeq_GeneBody.gtf .
#-----------------------------------------------------
#As described by TopHat2:
#<reads1_1[,...,readsN_1]> 	A comma-separated list of files containing reads in FASTQ or FASTA format. When running TopHat with paired-end reads, this should be the *_1 ("left") set of files.
#<[reads1_2,...readsN_2]> 	A comma-separated list of files containing reads in FASTA or FASTA format. Only used when running TopHat with paired end reads, and contains the *_2 ("right") set of files. The *_2 files MUST appear in the same order as the *_1 files. 
#-----------------------------------------------------
#We need a comma separated list for the *_1 ("left") set of files:
READ1_FILE_NAME=*_1.fastq.gz	
for i in ${READ1_FILE_NAME}; 
	do
	READ1=`echo $READ1$i,`; 
done
Left_List=${READ1%?}
echo
echo "Here is the Left_List:"
echo ${Left_List}
echo
#-----------------------------------------------------
#We need a comma separated list for the **_2 ("right") set of files:
READ2_FILE_NAME=*_2.fastq.gz	
for i in ${READ2_FILE_NAME}; 
	do
	READ2=`echo $READ2$i,`; 
done
Right_List=${READ2%?}
echo
echo "Here is the Right_List:"
echo ${Right_List}
echo
#-----------------------------------------------------
#Set the standard deviation for the ${DISTANCE_BT_READ_PAIR}
#--mate-std-dev <int>
#The standard deviation for the distribution on inner distances between mate pairs. The default is 20bp.  
#Looks like we use a wide standard deviation:
PARAM_STD=150
#-----------------------------------------------------
#Make output dir:
STORAGE_DIR=${Dataset_DIR}/${Sample_ID}/fastq
OUTPUT_DIR=$TMPDIR/tophat2
##############################
if [[ ! -d $OUTPUT_DIR ]]; 
then
mkdir $OUTPUT_DIR
fi
##############################
echo
echo 'List files in the scratch directory:'
echo

ls -alh

echo
echo 'Starting to run my commands'
echo

#echo
#echo 'Unzip files:'
#echo

#time gzip -d *.gz

#echo
#echo 'Finished unzipping'
#echo

# run my commands.
echo
echo 'Starting tophat'
echo
##################################################################################
#tophat --version
#TopHat v2.0.13
#tophat
#tophat: 
#TopHat maps short sequences from spliced transcripts to whole genomes.

#Usage:
#    tophat [options] <bowtie_index> <reads1[,reads2,...]> [reads1[,reads2,...]] \
#                                    [quals1,[quals2,...]] [quals1[,quals2,...]]
#-o/--output-dir                <string>    [ default: ./tophat_out         ]
#-G/--GTF                       <filename>  (GTF/GFF with known transcripts)
#-p/--num-threads               <int>       [ default: 1                   ]
#--library-type                 <string>    (fr-unstranded, fr-firststrand,
#                                                fr-secondstrand)
##################################################################################
#tophat --version
#TopHat v2.1.1
#tophat
#tophat: 
#TopHat maps short sequences from spliced transcripts to whole genomes.

#Usage:
#    tophat [options] <bowtie_index> <reads1[,reads2,...]> [reads1[,reads2,...]] \
#                                    [quals1,[quals2,...]] [quals1[,quals2,...]]

#Options:
#    -v/--version
#    -o/--output-dir                <string>    [ default: ./tophat_out         ]
#    --bowtie1                                  [ default: bowtie2              ]
#    -N/--read-mismatches           <int>       [ default: 2                    ]
#    --read-gap-length              <int>       [ default: 2                    ]
#    --read-edit-dist               <int>       [ default: 2                    ]
#    --read-realign-edit-dist       <int>       [ default: "read-edit-dist" + 1 ]
#    -a/--min-anchor                <int>       [ default: 8                    ]
#    -m/--splice-mismatches         <0-2>       [ default: 0                    ]
#    -i/--min-intron-length         <int>       [ default: 50                   ]
#    -I/--max-intron-length         <int>       [ default: 500000               ]
#    -g/--max-multihits             <int>       [ default: 20                   ]
#    --suppress-hits
#    -x/--transcriptome-max-hits    <int>       [ default: 60                   ]
#    -M/--prefilter-multihits                   ( for -G/--GTF option, enable
#                                                 an initial bowtie search
#                                                 against the genome )
#    --max-insertion-length         <int>       [ default: 3                    ]
#    --max-deletion-length          <int>       [ default: 3                    ]
#    --solexa-quals
#    --solexa1.3-quals                          (same as phred64-quals)
#    --phred64-quals                            (same as solexa1.3-quals)
#    -Q/--quals
#    --integer-quals
#    -C/--color                                 (Solid - color space)
#    --color-out
#    --library-type                 <string>    (fr-unstranded, fr-firststrand,
#                                                fr-secondstrand)
#    -p/--num-threads               <int>       [ default: 1                   ]
#    -R/--resume                    <out_dir>   ( try to resume execution )
#    -G/--GTF                       <filename>  (GTF/GFF with known transcripts)
#    --transcriptome-index          <bwtidx>    (transcriptome bowtie index)
#    -T/--transcriptome-only                    (map only to the transcriptome)
#    -j/--raw-juncs                 <filename>
#    --insertions                   <filename>
#    --deletions                    <filename>
#    -r/--mate-inner-dist           <int>       [ default: 50                  ]
#    --mate-std-dev                 <int>       [ default: 20                  ]
#    --no-novel-juncs
#    --no-novel-indels
#    --no-gtf-juncs
#    --no-coverage-search
#    --coverage-search
#    --microexon-search
#    --keep-tmp
#    --tmp-dir                      <dirname>   [ default: <output_dir>/tmp ]
#    -z/--zpacker                   <program>   [ default: gzip             ]
#    -X/--unmapped-fifo                         [use mkfifo to compress more temporary
#                                                 files for color space reads]

#Advanced Options:
#    --report-secondary-alignments
#    --no-discordant
#    --no-mixed

#    --segment-mismatches           <int>       [ default: 2                ]
#    --segment-length               <int>       [ default: 25               ]

#    --bowtie-n                                 [ default: bowtie -v        ]
#    --min-coverage-intron          <int>       [ default: 50               ]
#    --max-coverage-intron          <int>       [ default: 20000            ]
#    --min-segment-intron           <int>       [ default: 50               ]
#    --max-segment-intron           <int>       [ default: 500000           ]
#    --no-sort-bam                              (Output BAM is not coordinate-sorted)
#    --no-convert-bam                           (Do not output bam format.
#                                                Output is <output_dir>/accepted_hits.sam)
#    --keep-fasta-order
#    --allow-partial-mapping

#Bowtie2 related options:
#  Preset options in --end-to-end mode (local alignment is not used in TopHat2)
#    --b2-very-fast
#    --b2-fast
#    --b2-sensitive
#    --b2-very-sensitive

#  Alignment options
#    --b2-N                         <int>       [ default: 0                ]
#    --b2-L                         <int>       [ default: 20               ]
#    --b2-i                         <func>      [ default: S,1,1.25         ]
#    --b2-n-ceil                    <func>      [ default: L,0,0.15         ]
#    --b2-gbar                      <int>       [ default: 4                ]

#  Scoring options
#    --b2-mp                        <int>,<int> [ default: 6,2              ]
#    --b2-np                        <int>       [ default: 1                ]
#    --b2-rdg                       <int>,<int> [ default: 5,3              ]
#    --b2-rfg                       <int>,<int> [ default: 5,3              ]
#    --b2-score-min                 <func>      [ default: L,-0.6,-0.6      ]

#  Effort options
#    --b2-D                         <int>       [ default: 15               ]
#    --b2-R                         <int>       [ default: 2                ]

#Fusion related options:
#    --fusion-search
#    --fusion-anchor-length         <int>       [ default: 20               ]
#    --fusion-min-dist              <int>       [ default: 10000000         ]
#    --fusion-read-mismatches       <int>       [ default: 2                ]
#    --fusion-multireads            <int>       [ default: 2                ]
#    --fusion-multipairs            <int>       [ default: 2                ]
#    --fusion-ignore-chromosomes    <list>      [ e.g, <chrM,chrX>          ]

#    --fusion-do-not-resolve-conflicts          [this is for test purposes  ]

#SAM Header Options (for embedding sequencing run metadata in output):
#    --rg-id                        <string>    (read group ID)
#    --rg-sample                    <string>    (sample ID)
#    --rg-library                   <string>    (library ID)
#    --rg-description               <string>    (descriptive string, no tabs allowed)
#    --rg-platform-unit             <string>    (e.g Illumina lane ID)
#    --rg-center                    <string>    (sequencing center name)
#    --rg-date                      <string>    (ISO 8601 date of the sequencing run)
#    --rg-platform                  <string>    (Sequencing platform descriptor)

#    for detailed help see http://ccb.jhu.edu/software/tophat/manual.shtml
##################################################################################
#------------------------------------------------------------------------------
#Jobs are taking too long to complete
#Need the option: --no-coverage-search
#------------------------------------------------------------------------------
echo "tophat --no-coverage-search -o ${OUTPUT_DIR} -G ${ANNOTATION_FILE} -p $NSLOTS --library-type ${STRANDEDNESS} -r ${DISTANCE_BT_READ_PAIR} --mate-std-dev ${PARAM_STD} Bowtie2Index/genome ${Left_List} ${Right_List}"
tophat --no-coverage-search -o ${OUTPUT_DIR} -G ${ANNOTATION_FILE} -p $NSLOTS --library-type ${STRANDEDNESS} -r ${DISTANCE_BT_READ_PAIR} --mate-std-dev ${PARAM_STD} Bowtie2Index/genome ${Left_List} ${Right_List}
##################################################################################
echo
echo 'Ending tophat'
echo

echo
echo 'Starting samtools commands'
echo
#----------------------------------------------------------
#Site for explaining flag definitions:
#http://broadinstitute.github.io/picard/explain-flags.html
#http://davetang.org/muse/2014/03/06/understanding-bam-flags/
#0x100 = not primary alignment
#samtools view options:
#-f INT   required flag, 0 for unset [0]
#-F INT   filtering flag, 0 for unset [0]
#----------------------------------------------------------
#Require that the read record is not primary alignment:
samtools view -f 0x100 -b ${OUTPUT_DIR}/accepted_hits.bam > ${OUTPUT_DIR}/non_primary.bam
#Filter out all the read records that are not primary alignment
#Otherwords: extract all primary alignment reads
samtools view -F 0x100 -b ${OUTPUT_DIR}/accepted_hits.bam > ${OUTPUT_DIR}/primary.bam
#Generate index for a bam file
samtools index ${OUTPUT_DIR}/accepted_hits.bam
samtools index ${OUTPUT_DIR}/primary.bam
#Collect statistics for a bam file
samtools flagstat ${OUTPUT_DIR}/non_primary.bam > ${OUTPUT_DIR}/statistics_for_nonprimary_reads.txt
samtools flagstat ${OUTPUT_DIR}/primary.bam > ${OUTPUT_DIR}/statistics_for_primary_reads.txt
samtools flagstat ${OUTPUT_DIR}/accepted_hits.bam > ${OUTPUT_DIR}/statistics_for_all_accepted_reads.txt
samtools flagstat ${OUTPUT_DIR}/unmapped.bam > ${OUTPUT_DIR}/statistics_for_unmapped_reads.txt
#numunique_hits.txt contains the mapping statistics of all mappable reads, ***excluding*** multi-mappable reads (identical to b)). This is not applicable for paired-end reads.
samtools view ${OUTPUT_DIR}/accepted_hits.bam | cut -f 1 | sort | uniq | wc -l > numunique_hits.txt
#We now want to get the uniquely mapped reads:
#------------------------------------------------------------------------------
cd ${OUTPUT_DIR}
#http://www.researchgate.net/post/How_can_I_get_uniquely_mapped_reads_from_Tophatv2012
#We can filter the *_primary.bam
#First I need to get the header lines
#------------------------------------------------------------------------------
#Easier alternative: samtools view with (-H) option:
#-H       print header only (no alignments)
samtools view -H 'primary.bam' > Header.txt
#samtools view (without option) will print alignments (no header) :
samtools view 'primary.bam' > temp1.sam
#Run grep command to extract the uniquely mapped reads:
#http://onetipperday.blogspot.com/2012/07/deeply-understanding-sam-tags.html
#NH is "Number of reported alignments that contains the query in the current record", which should be >1 for multiple-mapped reads.
grep -w "NH:i:1" temp1.sam > temp2.sam
#Add header:
cat Header.txt temp2.sam > temp3.sam
#Convert SAM to BAM:
samtools view -bS temp3.sam > temp4.bam
#Sort the BAM file:
#Usage:   samtools sort [options] <in.bam> <out.prefix>
#Options: -n        sort by read name
#Without any options it will sort by coordinates
samtools sort temp4.bam 'primary_unique'
#http://sourceforge.net/p/samtools/mailman/message/31318733/
#The samtools index command expects the input file to be coordinate sorted
#Generate index for a bam file
samtools index 'primary_unique'.bam
#Collect statistics for a bam file
samtools flagstat 'primary_unique'.bam > statistics_for_'primary_unique'_reads.txt
#Remove temp files:
rm Header.txt
rm temp*.sam
rm temp*.bam
cd ..
#------------------------------------------------------------------------------
echo
echo 'Ending samtools commands'
echo
##################################################################################
#Output file descriptions:
#1) accepted_hits.bam -> The mapped reads for your sample. If your reads are mappable to several places, this file will have up to 20 of those.
#3) Several mapping statistics:
#a) statistics_for_all_accepted_reads.txt contains the mapping statistics of all mappable reads, which include multi-mappable reads.
#b) statistics_for_primary_reads.txt contains the mapping statistics of all mappable reads, ***excluding*** multi-mappable reads.
#c) statistics_for_nonprimary_reads.txt contains the mapping statistics for multi-mappable reads.
#d) statistics_for_unmapped_reads.txt contains the mapping statistics for unmapped rads.
#e) numunique_hits.txt contains the mapping statistics of all mappable reads, ***excluding*** multi-mappable reads (identical to b)). This is not applicable for paired-end reads.
#---------------------------------------------------------------------------------
#primary.bam -> The mapped reads that only contain the best mappable place for each of your read that was able to be mapped.
#junctions.bed. A UCSC BED track of junctions reported by TopHat. Each junction consists of two connected BED blocks, where each block is as long as the maximal overhang of any read spanning the junction. The score is the number of alignments spanning the junction.
#insertions.bed and deletions.bed. UCSC BED tracks of insertions and deletions reported by TopHat. 
#---------------------------------------------------------------------------------
#primary_unique.bam -> The subset of the primary.bam where each read only has 1 reported alignment
#In other words, these are the uniquely mapped reads
#The uniquely mapped reads is the read set we want to use for both counting and visualization
#---------------------------------------------------------------------------------
##################################################################################
#remove stuff that we dont need
#For DESeq, we need the primary_unique.bam
rm ${OUTPUT_DIR}/non_primary.bam
rm ${OUTPUT_DIR}/accepted_hits.bam
rm ${OUTPUT_DIR}/accepted_hits.bam.bai
rm ${OUTPUT_DIR}/primary.bam
rm ${OUTPUT_DIR}/primary.bam.bai
rm ${OUTPUT_DIR}/unmapped.bam

#Rename the output files
Output_List=${OUTPUT_DIR}/*
for output_file in $Output_List
do
output_file_name=`basename $output_file`
#Add the ${Sample_ID} as a file name prefix
mv ${OUTPUT_DIR}/$output_file_name ${OUTPUT_DIR}/${Sample_ID}'_'$output_file_name
done

#---------------------------------------------------------------------------------
echo
echo 'Starting spliced_read_counts'
echo
#cd into the ${OUTPUT_DIR}
cd ${OUTPUT_DIR}
######################
OUTPUT_FILE=$OUTPUT_DIR/${Sample_ID}'_spliced_read_counts.txt'
######################
if [ -f $OUTPUT_FILE ]
then 
rm $OUTPUT_FILE
else
touch $OUTPUT_FILE
fi
######################
#It would be useful to know the number of mapped reads that have skipped regions from the reference genome
#This information is in the CIGAR string
#http://zenfractal.com/2013/06/19/playing-with-matches/
#Aligned reads in a SAM or BAM file typically have a Compact Idiosyncratic Gapped Alignment Report (CIGAR) string that expresses how the read is mapped to the reference genome.
#The CIGAR string is in the 6th column:
#https://samtools.github.io/hts-specs/SAMv1.pdf
#M stands for: alignment match (can be a sequence match or mismatch)
#N stands for: skipped region from the reference
SKIPPED_READS=$(samtools view ${Sample_ID}'_primary_unique.bam' | awk '($6 ~ /N/)' | cut -f6 | wc -l)
#Need a regular expression to specifically find these spliced reads:
#Bunch of numbers then "M", bunch of numbers then "N", bunch of numbers then "M":
SPLICED_READS=$(samtools view ${Sample_ID}'_primary_unique.bam' | awk '($6 ~ /[0-9]+[M][0-9]+[N][0-9]+[M]/)' | cut -f6 | wc -l )
#Print header to output file:
echo 'SAMPLE_ID' $'\t''skipped_region_READS' $'\t''splice_junction_READS' >> $OUTPUT_FILE
#Print to output file
echo $Sample_ID $'\t'$SKIPPED_READS $'\t'$SPLICED_READS >> $OUTPUT_FILE
#cd out of ${OUTPUT_DIR}
cd ..
echo
echo 'Ending spliced_read_counts'
echo
#---------------------------------------------------------------------------------

# copy the output files to users storage dir
cp *.txt ${OUTPUT_DIR}
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 "=========================================================="