#!/bin/bash
##################################################################################
#Andy Rampersaud, 01.27.16
#Adapted from tophat/paired_end_mapping scripts by Tisha Melia
##################################################################################
#Assumptions for this job to run correctly:
#1. This script is intended to map paired end RNA-Seq reads
#2. You are mapping to mm9 (ucsc) reference genome
#3. Your data is organized in the following way:
#You have a data set dir such as:
#wax-es/aramp10/G83_Samples
#Within this dir you have sample specific folders such as:
#G83_M1
#G83_M2
#G83_M3
#G83_M4
#Within each sample specific folder you have a "fastq" folder with:
#Files: *_R1_*.fastq.gz and *_R2_*.fastq.gz such as:
#Waxman-TP17_CGATGT_L007_R1_001.fastq.gz
#Waxman-TP17_CGATGT_L007_R2_001.fastq.gz
##################################################################################
#Fill in the following information:
##################################################################################
#Information about your data set
#As mentioned above, you should have a data set dir containing your sample specific folders:
Dataset_DIR=/projectnb/wax-es/aramp10/G83_Samples
##################################################################################
#Samples to process
#To facilitate processing of samples in parallel we can use a text file that lists the samples to analyze
#Note: this text file is still valid even if there is only one sample to process
#You need to have a "Sample_Labels" dir within your Dataset_DIR
#Within the Sample_Labels dir have a Sample_Labels.txt such that:
################################################
#The text file is formatted like the following:
#----------------------------------------------
#Sample_DIR	Sample_ID	Description
#Sample_Waxman-TP17	G83_M1	Male 8wk-pool 1
#Sample_Waxman-TP18	G83_M2	Male 8wk-pool 2
#Sample_Waxman-TP19	G83_M3	Female 8wk-pool 1
#Sample_Waxman-TP20	G83_M4	Female 8wk-pool 2
#----------------------------------------------
#The 1st column: The Sample_DIR name
#The 2nd column: Waxman Lab Sample_ID 
#The 3rd column: Sample's description 
################################################
#Sample_Labels_DIR=${Dataset_DIR}/Sample_Labels
##################################################################################
#Bowtie2Index_DIR
#Need this dir for the Bowtie2 Index
#Feel free to use my dir but it's better practice to have a copy in your own Dataset_DIR
#Bowtie2Index_DIR=/projectnb/wax-es/aramp10/Bowtie2
##################################################################################
#GTF_Files_DIR
#Need this dir that contains the various GTF files
#Feel free to use my dir but it's better practice to have a copy in your own Dataset_DIR
#GTF_Files_DIR=/projectnb/wax-es/aramp10/GTF_Files
##################################################################################
#STRANDEDNESS
#There are three types of library construction. Please select your type of library construction
#
#1.) fr-unstranded (Standard Illumina) -> Reads from the left-most end of the fragment (in transcript coordinates) map to the transcript strand, and the right-most end maps to the opposite strand.
#2.) fr-firststrand (dUTP, NSR, NNSR) -> Same as above except we enforce the rule that the right-most end of the fragment (in transcript coordinates) is the first sequenced (or only sequenced for single-end reads). Equivalently, it is assumed that only the strand generated during first strand synthesis is sequenced.
#3.) fr-secondstrand (Ligation, Standard SOLiD) -> Same as above except we enforce the rule that the left-most end of the fragment (in transcript coordinates) is the first sequenced (or only sequenced for single-end reads). Equivalently, it is assumed that only the strand generated during second strand synthesis is sequenced.
#
# NOTE: If your dataset comes from Spira lab, your library type is likely to be "fr-unstranded". 
#       If your dataset comes from MIT/Encode, your library type is likely to be "fr-firststrand".
#If RNA-Seq library was prepared in-house using the NEBNext Ultra Directional kit, use option “fr-firststrand”
#If your data is unstranded then use "fr-unstranded"
#---------------------------------------------------------------------------------
#Choose one:
#STRANDEDNESS="fr-unstranded"
#STRANDEDNESS="fr-firststrand"
#STRANDEDNESS="fr-secondstrand"
#---------------------------------------------------------------------------------
##################################################################################
#Set the distance between your read pair
#This is the option description:
#---------------------------------------------------------------------------------
#-r/--mate-inner-dist <int>:
#This is the expected (mean) inner distance between mate pairs. For, example, for paired end runs with fragments selected at 300bp, where each end is 50bp, you should set -r to be 200. The default is 50bp. 
#---------------------------------------------------------------------------------
#There should be a corresponding Bioanalyzer report and/or gel analysis to determine the fragment length
#DISTANCE_BT_READ_PAIR -> (fragment length bp) - (2)*(read length bp)
#As in the above example:
#DISTANCE_BT_READ_PAIR -> (300 bp) - (2)*(50 bp) = 200 bp
#---------------------------------------------------------------------------------
#This script will do the arithmetic to calculate your (--mate-inner-dist)
#Input the Bioanalyzer length (bp) (directly from Bioanalyzer tracings)
#BIOANALYZER_LEN=277
#Note: Lengths directly from Bioanalyzer tracings still include adaptors (ligated to both ends)
#Usually the adaptor length is 60bp (60bp on each end: 120bp total)
#Input the adaptor length (bp) for one end:
#ADAPTOR_LEN=60
#Input the read length (bp)
#Should have this from the Read_Length job
#READ_LEN=99
#---------------------------------------------------------------------------------
##################################################################################
##################################################################################
#Please: DO NOT EDIT CODE BELOW
##################################################################################
##################################################################################
#Gene Annotation file to use
#---------------------------------------------------------------------------------
#For read mapping purposes we want the full gene body information (use the RefSeq_GeneBody.gtf)
#Presence of splice junctions should be based on the full gene structure for any isoform of a gene symbol
#For read counting purposes we can use either 
#1. RefSeq_GeneBody.gtf
#2. Intron_Only_Regions.gtf
#3. Exon_Only_Regions.gtf
#---------------------------------------------------------------------------------
#ANNOTATION_FILE=RefSeq_GeneBody.gtf
#---------------------------------------------------------------------------------
##################################################################################
#G83 Data:
#Samples have 99bp paired end reads
#Confirmed read length from "Read_Length" job
#The average of the fragment lengths
#Average(159, 162, 170, 166) = 164 bp (without adaptors)
#(Information from ChIPSeq_samples_57.xlsx)
#Average(272, 279, 280, 275) = 277 bp (with adaptors)
#(Information from G83_mapping_stats.docx)
#r parameter: (277 bp bp) (2)(60 bp) - (2)(99 bp) = -41 bp
#---------------------------------------------------------------------------------
#Do the calculation:
#FRAGMENT_LEN=$(echo "scale=4;$BIOANALYZER_LEN-(2*$ADAPTOR_LEN)" | bc)
#DISTANCE_BT_READ_PAIR=$(echo "scale=4;$FRAGMENT_LEN-(2*$READ_LEN)" | bc)
##################################################################################
#Time hour limit
#On SCC a 12-hour runtime limit is enforced on all jobs, unless specified explicitly. 
#A runtime limit can be specified in the format "hh:mm:ss"
#Dont change the following time limit value unless you know that your job is going to go over 12 hrs 
#TIME_LIMIT="12:00:00"
##################################################################################