#!/bin/bash
#$ -m e
#$ -pe smp 4
#$ -l mem_free=36G,h_vmem=38G
#$ -cwd


### IMPORTANT ###

# reads should be gzipped fastq files
# name of the file should have the following format:
# samplename_readpair_replicate this could look like this
# total_R1_A.fastq.gz and total_R2_A.fastq.gz
# the third argument is "" or "QRNA" this depends on if the library was prepared with the qRNA kit
# usage: pairedRNAseqPipeline.sh  total_R1_A.fastq.gz total_R2_A.fastq.gz QRNA
# esample usage on cluster: qsub -V -e total_A.e -o total_A.o pairedRNAseqPipeline.sh  total_R1_A.fastq.gz total_R2_A.fastq.gz QRNA


# TO DO:
# QRNA read counting
# Move log files to analysis folder
# Make all relative paths
# Create summaries



THRESHOLD="5" # Minimal number for junction reads
strandSPthresh="0.75" # Minimal strand-specificity threshold
STAR="/data/ohler/Neel/bioTools/STAR-STAR_2.4.0j/bin/Linux_x86_64/STAR" # STAR Executable
GTF="/data/ohler/Neel/Genomes/human/hg19/annotation/v19/gencode.v19.chr_patch_hapl_scaff.annotation.ERCC.gtf" # GTF File
CPS="/data/ohler/Neel/Genomes/human/hg19/annotation/v19/gencode.v19.cps" # CPS FILE
INDEX="/data/ohler/Neel/Genomes/human/hg19/STARindex/100ERCC" # STAR Index
FASTA="/data/ohler/Neel/Genomes/human/hg19/genomes/v19/GRCh37.p13.genome.fa" # Genome
CHRSIZE="/data/ohler/Neel/metLab/data/accessory/gencodev19.lengths"
SCRIPTS="/data/ohler/Neel/bioTools/" #Scripts
#QC="/data/ohler/Neel/bioTools/usr/local/bin/" # RSeqQC scripts
#QC="~/.guix-profile/bin/" # RSeqQC scripts
QC="" # RSeqQC scripts
QChouse="/data/ohler/Neel/bioTools/RSeQC-2.4/accessory/hg19.HouseKeepingGenes.bed" # RSeqQC files
QCgene="/data/ohler/Neel/bioTools/RSeQC-2.4/accessory/hg19_GENCODE_GENE_V19_comprehensive.bed" # RSeqQC files
rRNA_bwt="/data/ohler/Neel/Genomes/human/hg19/rRNA/bowtie/v19rRNA" # rRNA bowtie index
BEDTOOLS=""
strictPARS="
--runThreadN 4 \
--outFilterType BySJout --outFilterMultimapNmax 20 --outFilterMismatchNmax 999 \
--outFilterMismatchNoverLmax 0.04 --outFilterIntronMotifs RemoveNoncanonicalUnannotated \
--alignIntronMin 20 --alignIntronMax 1000000 \
--alignMatesGapMax 10000 \
--alignSJoverhangMin 8 --alignSJDBoverhangMin 1 \
--sjdbScore 1 \
--chimSegmentMin 15 --chimJunctionOverhangMin 15 \
--outSJfilterReads Unique --outSJfilterCountUniqueMin 5 2 2 2 \
--outReadsUnmapped Fastx \
--outSAMtype BAM SortedByCoordinate --outSAMmode Full \
--outSAMattributes NH HI AS NM MD --outSAMprimaryFlag OneBestScore --outSAMstrandField intronMotif \
--outWigType bedGraph "
#--quantMode TranscriptomeSAM  \


RSEM="/data/ohler/Lorenzo/bins/rsem_1_2_11/"
bowPATH="/data/ohler/Neel/bioTools/bowtie/"
rsemINDEX="/data/ohler/Neel/Genomes/human/hg19/RSEMercc/genome"
pri_rsemINDEX="/data/ohler/Neel/Genomes/human/hg19/RSEMercc/pri"
#rsemINDEX="/data/ohler/Neel/Genomes/human/hg19/RSEMercc/genome"
rsemPARS=" -p 4 --no-bam-output --estimate-rspd -q --paired-end "
    


### DO NOT MODIFY BELOW ####
[ $# -eq 0 ] && { echo "No fastq files provided"; exit 1; }

date1=$(date +"%s")

fullfile1=$1
fullfile2=$2
QRNA=$3
sample=${fullfile1%_*_*.fastq.gz}
fname=$(basename $fullfile1)
#lib=${fname%.fastq.gz}
lib=${fullfile1%.fastq.gz}
rep=${lib##*_R*_}


if [ "${QRNA}" = "QRNA" ];
	then
		bowExtra=" -5 9 "
		starExtra=" --clip5pNbases 9 "
	else
		bowExtra=""
		starExtra=""
fi

# quantify rRNA processing
if [ ! -d rRNAproc/  ];
	then
		mkdir rRNAproc
fi

if [ ! -s rRNAproc/${sample}_${rep}_coverage.bed ] ;
	then
		echo "pipeStage: quantify rRNA processing, START @ $(date)"
		cd rRNAproc
		mkfifo ${sample}_${rep}_1
		mkfifo ${sample}_${rep}_2
		gunzip -c ../fastq/${fullfile1} > ${sample}_${rep}_1 &
		gunzip -c ../fastq/${fullfile2} > ${sample}_${rep}_2 &
		${SCRIPTS}bowtie/bowtie ${bowExtra} -p 4 -q -X 1000 --fr --best /data/ohler/Neel/Genomes/human/hg19/rRNA/bowtie/tx_rRNA  -1 ${sample}_${rep}_1 -2 ${sample}_${rep}_2 --sam ${sample}_${rep}_rRNA.sam
		mkdir ${sample}_${rep}_tmpdir
		${SCRIPTS}sambamba view -f bam -S ${sample}_${rep}_rRNA.sam > ${sample}_${rep}_rRNA_pre.bam 
		${SCRIPTS}sambamba sort --tmpdir=${sample}_${rep}_tmpdir -o ${sample}_${rep}_rRNA.bam  ${sample}_${rep}_rRNA_pre.bam 
		${BEDTOOLS}coverageBed -s -abam ${sample}_${rep}_rRNA.bam -b /data/ohler/Neel/Genomes/human/hg19/rRNA/tx_rRNA/tx_rRNA.bed > ${sample}_${rep}_coverage.bed
		rm ${sample}_${rep}_rRNA_pre.bam ${sample}_${rep}_rRNA.sam
		rm ${sample}_${rep}_1 ${sample}_${rep}_2
		rm -r ${sample}_${rep}_tmpdir
		cd ../
		
                echo "pipeStage: quantify rRNA processing, END @ $(date)"
		
	else
		wait
                echo "pipeStage: quantify rRNA processing, Already Complete @ $(date)"
fi






# Remove rRNA reads with bowtie
if [ ! -d filtFastq/  ];
	then
		mkdir filtFastq
fi

if [ ! -s filtFastq/filt_${sample}_${rep}_2.fastq ] ;
	then
		echo "pipeStage: Remove rRNA reads with bowtie, START @ $(date)"
		cd filtFastq
		mkfifo ${sample}_${rep}_1
		mkfifo ${sample}_${rep}_2
		gunzip -c ../fastq/${fullfile1} > ${sample}_${rep}_1 &
		gunzip -c ../fastq/${fullfile2} > ${sample}_${rep}_2 &
		${SCRIPTS}bowtie/bowtie ${bowExtra} -p 4 -q -X 1000 --fr -k 1 --best --un filt_${sample}_${rep}.fastq ${rRNA_bwt} -1 ${sample}_${rep}_1 -2 ${sample}_${rep}_2 --sam ${sample}_${rep}_rRNA.sam
		wait
                
		
		if [ "${QRNA}" = "QRNA" ];
			then
				cat filt_${sample}_${rep}_1.fastq | /data/ohler/Neel/bioTools/bin/fastx_trimmer -Q33 -f 10 -o trim_filt_${sample}_${rep}_1.fq
				cat filt_${sample}_${rep}_2.fastq | /data/ohler/Neel/bioTools/bin/fastx_trimmer -Q33 -f 10 -o trim_filt_${sample}_${rep}_2.fq
		
                fi
		rm ${sample}_${rep}_1 ${sample}_${rep}_2
                cd ../
		
                echo "pipeStage: Remove rRNA reads with bowtie, END @ $(date)"
		
	else
		wait
                echo "pipeStage: Remove rRNA reads with bowtie, Already Complete @ $(date)"
fi


# Perform alignment with STAR
if [ ! -d bam/  ];
	then
		mkdir bam
fi

if [ ! -s bam/${sample}_${rep}_Aligned.sortedByCoord.out.bam ] ;
	then
		echo "pipeStage: Perform Alignment with STAR, START @ $(date)"
		cd bam
	#    ${STAR} --genomeDir ${INDEX} --readFilesCommand zcat --readFilesIn ../$fullfile1 ../$fullfile2 ${strictPARS} ${starExtra}  --outFileNamePrefix ${sample}_${rep}_ 
	#    ${STAR} --genomeDir ${INDEX} --readFilesCommand zcat --readFilesIn ../fastq/$fullfile1 ../fastq/$fullfile2 ${strictPARS} ${starExtra}  --outFileNamePrefix ${sample}_${rep}_
		${STAR} --genomeDir ${INDEX} --readFilesIn ../filtFastq/filt_${sample}_${rep}_1.fastq ../filtFastq/filt_${sample}_${rep}_2.fastq ${strictPARS} ${starExtra} --outFileNamePrefix ${sample}_${rep}_ 
	#    ${QC}split_bam.py -i ${sample}_${rep}_Aligned.sortedByCoord.out.bam -r ${rRNA_bed}  -o ${sample}_${rep}_rRNA
		${SCRIPTS}sambamba index -t 4 ${sample}_${rep}_Aligned.sortedByCoord.out.bam
		cd ../
		echo "pipeStage: Perform Alignment with STAR, END @ $(date)"

	else
		wait
                echo "pipeStage: Perform Alignment with STAR, Already Complete @ $(date)"
fi

wait

# Determine the strandedness of the library: Important for numerous downstream analysis
${QC}infer_experiment.py -i bam/${sample}_${rep}_Aligned.sortedByCoord.out.bam -r ${QChouse} -s 500000 > ${sample}_${rep}.strand
if [ "$(grep 1++ ${sample}_${rep}.strand | cut -d " " -f 7)" \> "${strandSPthresh}" ];
    then
        ssjStrand=" -read1 0 -read2 1 "
        rsemStrand=" --forward-prob 1 "
        htseqStrand=" --stranded=yes "
        rseqStrand=$(grep 1++ ${sample}_${rep}.strand | cut -d " " -f 6 | sed 's/://g' | sed "s/\"//g")

elif [ "$(grep 2++ ${sample}_${rep}.strand | cut -d " " -f 7)" \> "${strandSPthresh}" ];
    then
        ssjStrand=" -read1 1 -read2 0 "
        rsemStrand=" --forward-prob 0 "
        htseqStrand=" --stranded=reverse "
        rseqStrand=$(grep 2++ ${sample}_${rep}.strand | cut -d " " -f 6 | sed 's/://g' | sed "s/\"//g")
wait
fi


# Create bigwig from bedgraph
if [ ! -d bigwig/  ];
	then
		mkdir bigwig
fi

if [ ! -s bigwig/${sample}_${rep}_fwd.bw ] ;
	then
		echo "pipeStage: Create bigwig from bedgraph, START @ $(date)"
		 if [ "$(grep 1++ ${sample}_${rep}.strand | cut -d " " -f 7)" \> "${strandSPthresh}" ];
			then    
				awk '{$4=$4*-1;print}' bam/${sample}_${rep}_Signal.Unique.str2.out.bg > bam/${sample}_${rep}_rev.bg      
				${SCRIPTS}bedGraphToBigWig bam/${sample}_${rep}_rev.bg ${CHRSIZE} bigwig/${sample}_${rep}_rev.bw
				${SCRIPTS}bedGraphToBigWig bam/${sample}_${rep}_Signal.Unique.str1.out.bg ${CHRSIZE} bigwig/${sample}_${rep}_fwd.bw

		elif [ "$(grep 2++ ${sample}_${rep}.strand | cut -d " " -f 7)" \> "${strandSPthresh}" ];
			then
				awk '{$4=$4*-1;print}' bam/${sample}_${rep}_Signal.Unique.str1.out.bg > bam/${sample}_${rep}_rev.bg      
				${SCRIPTS}bedGraphToBigWig bam/${sample}_${rep}_rev.bg ${CHRSIZE} bigwig/${sample}_${rep}_rev.bw
				${SCRIPTS}bedGraphToBigWig bam/${sample}_${rep}_Signal.Unique.str2.out.bg ${CHRSIZE} bigwig/${sample}_${rep}_fwd.bw
				echo "pipeStage: Create bigwig from bedgraph, END @ $(date)"
		fi

	else
		
                echo "pipeStage: Create bigwig from bedgraph, Already Complete @ $(date)"
fi



# Calculate primary transcript expression levels with RSEM
if [ ! -d pri_rsem/ ];
    then
        mkdir pri_rsem
fi

if [ ! -d pri_rsem/pri_${sample}_${rep} ];
    then
        mkdir pri_rsem/pri_${sample}_${rep}
fi
	
if [ ! -s pri_rsem/pri_${sample}_${rep}/pri_${sample}_${rep}.genes.results ] ;
    then
	echo "pipeStage: Calculate expression levels with pri_RSEM, START @ $(date)"

	if [ "${QRNA}" = "QRNA" ];
            then
                ${RSEM}rsem-calculate-expression --bowtie-path ${bowPATH} ${rsemStrand} ${rsemPARS} \
                filtFastq/trim_filt_${sample}_${rep}_1.fq filtFastq/trim_filt_${sample}_${rep}_2.fq \
                ${pri_rsemINDEX} pri_rsem/pri_${sample}_${rep}/pri_${sample}_${rep}
                pid2=$!
                wait
                ${RSEM}rsem-plot-model pri_rsem/pri_${sample}_${rep}/pri_${sample}_${rep} pri_rsem/pri_${sample}_${rep}/pri_${sample}_${rep}_model.pdf
            else
                ${RSEM}rsem-calculate-expression --bowtie-path ${bowPATH} ${rsemStrand} ${rsemPARS} \
                filtFastq/filt_${sample}_${rep}_1.fastq filtFastq/filt_${sample}_${rep}_2.fastq \
                ${pri_rsemINDEX} pri_rsem/pri_${sample}_${rep}/pri_${sample}_${rep}
                pid2=$!
                wait
                ${RSEM}rsem-plot-model pri_rsem/pri_${sample}_${rep}/pri_${sample}_${rep} pri_rsem/pri_${sample}_${rep}/pri_${sample}_${rep}_model.pdf
	fi
	echo "pipeStage:Calculate expression levels with pri_RSEM, END @ $(date)"

	else
    
	wait ${pid2}    
	echo "pipeStage: Calculate expression levels with pri_RSEM, Already Complete @ $(date)"
fi

wait

# QC analysis and metrics
if [ ! -d qc/  ];
	then
		mkdir qc
fi

if [ ! -s qc/${sample}_${rep}_cat.txt ] ;
	then
		echo "pipeStage: QC analysis, START @ $(date)"
		${QC}junction_annotation.py -i bam/${sample}_${rep}_Aligned.sortedByCoord.out.bam -r ${QCgene} -o qc/${sample}_${rep}
		${QC}geneBody_coverage.py -i bam/${sample}_${rep}_Aligned.sortedByCoord.out.bam -r ${QChouse} -o qc/${sample}_${rep}
		${QC}junction_saturation.py -i bam/${sample}_${rep}_Aligned.sortedByCoord.out.bam -r ${QCgene} -o qc/${sample}_${rep}
		${QC}read_distribution.py -i bam/${sample}_${rep}_Aligned.sortedByCoord.out.bam -r ${QCgene} > qc/${sample}_${rep}_cat.txt
		${QC}RPKM_count.py -i bam/${sample}_${rep}_Aligned.sortedByCoord.out.bam -r ${QCgene} -o qc/${sample}_${rep} -u -d "${rseqStrand}"
		echo "pipeStage: QC analysis, END @ $(date)"

	 else
		
		echo "pipeStage: QC analysis, Already Complete @ $(date)"
fi

# Calculate splice junction count and intron splicing info
if [ ! -d ssj/ ];
then
    mkdir ssj
fi

if [ ! -s ssj/${sample}_${rep}_BP.txt ];
	then
		echo "pipeStage: Calculate splice junction count junction info, START @ $(date)"
		${SCRIPTS}bam2ssj -cps ${CPS} -bam bam/${sample}_${rep}_Aligned.sortedByCoord.out.bam -out ssj/${sample}_${rep}.ssj -log ssj/${sample}_${rep}_ssj.log -v ${ssjStrand}

		cd ssj/
		DIR=$(pwd)
		perl ${SCRIPTS}parseSSJ.pl ${sample}_${rep}.ssj ${THRESHOLD}
		wait
		${BEDTOOLS}fastaFromBed -s -name -fi ${FASTA} -bed ${sample}_${rep}_5ss.bed -fo ${sample}_${rep}_5ss.fa
		wait
		${BEDTOOLS}fastaFromBed -s -name -fi ${FASTA} -bed ${sample}_${rep}_3ss.bed -fo ${sample}_${rep}_3ss.fa
		wait
		${BEDTOOLS}fastaFromBed -s -name -fi ${FASTA} -bed ${sample}_${rep}_BP.bed -fo ${sample}_${rep}_BP.fa
		wait
		${BEDTOOLS}nucBed -s -fi ${FASTA} -bed ${sample}_${rep}.bed  | cut -f 8 | sed '1s/8_pct_gc/GC/' >  ${sample}_${rep}_GC.txt 
		
		grep exon ../qc/${sample}_${rep}_read_count.xls |  awk '{if (  $6=="+" && $7 > 1 ) { print $1 "\t" $2 "\t" $3 "\t"$4 "\t" $7/sqrt($3-$2) "\t" $6}}' > ${sample}_${rep}_exons.bed
		grep exon ../qc/${sample}_${rep}_read_count.xls |  awk '{if (  $6=="-" && $8 > 1 ) { print $1 "\t" $2 "\t" $3 "\t"$4 "\t" $8/sqrt($3-$2) "\t" $6}}' >> ${sample}_${rep}_exons.bed
		closestBed -s -d -D a -id -io  -a ${sample}_${rep}.bed -b ${sample}_${rep}_exons.bed | awk '{if (  $13=="-1" ) { print $0}}' > ${sample}_${rep}_upstreamExontmp.bed
		closestBed -s -d -D a -iu -io  -a ${sample}_${rep}.bed -b ${sample}_${rep}_exons.bed | awk '{if (  $13=="1" ) { print $0}}' > ${sample}_${rep}_downstreamExontmp.bed
		Rscript ${SCRIPTS}selectAdjacentExon.R ${sample}_${rep}_upstreamExontmp.bed ${sample}_${rep}_downstreamExontmp.bed ${sample}_${rep}
		nucBed -s -fi /data/ohler/Neel/Genomes/human/hg19/genomes/v19/GRCh37.p13.genome.fa  -bed ${sample}_${rep}_upEx.bed   | awk ' { print $4 "\t" ($3-$2) "\t" $8}'  | sed '1s/8_pct_gc/GC/' | sed '1s/1/exonLength/' | sed '1s/4_usercol/seq_id/'  > ${sample}_${rep}_upEx_lGC.txt 
		nucBed -s -fi /data/ohler/Neel/Genomes/human/hg19/genomes/v19/GRCh37.p13.genome.fa  -bed ${sample}_${rep}_downEx.bed   | awk ' { print $4 "\t" ($3-$2) "\t" $8}'  | sed '1s/8_pct_gc/GC/' | sed '1s/1/exonLength/' | sed '1s/4_usercol/seq_id/' > ${sample}_${rep}_downEx_lGC.txt 
		

		fastaFromBed -s -name -tab -fi /data/ohler/Neel/Genomes/human/hg19/genomes/v19/GRCh37.p13.genome.fa -bed /data/ohler/Neel/metLab/data/labeling/ssj/Lab7_A_downEx.bed -fo Lab7_A_downEx.txt
		
		# Use MaxEnt to calculate 5' and 3' splice site scores
		cd ${SCRIPTS}fordownload/
		perl score5.pl ${DIR}/${sample}_${rep}_5ss.fa > ${DIR}/${sample}_${rep}_5ss.score
		wait
		perl score3.pl ${DIR}/${sample}_${rep}_3ss.fa > ${DIR}/${sample}_${rep}_3ss.score
		wait
		
		# Use svm-bpfinder to calculate polypyrimidine tract branchpoint scores
		cd ${SCRIPTS}svm-bpfinder/
		./svm_bpfinder.py ${DIR}/${sample}_${rep}_BP.fa Hsap 75 > ${DIR}/${sample}_${rep}_BP.txt
		
		
		
		cd ${DIR}
		sed -i '1s/^/5ssSeq\t5ssScore\n/' ${sample}_${rep}_5ss.score
		sed -i '1s/^/3ssSeq\t3ssScore\n/' ${sample}_${rep}_3ss.score
		paste ${sample}_${rep}.out ${sample}_${rep}_5ss.score ${sample}_${rep}_3ss.score  ${sample}_${rep}_GC.txt > ${sample}_${rep}.txt
		# rm ${sample}_${rep}*.bed rm ${sample}_${rep}*.ssj  
		#rm ${sample}_${rep}*streamExontmp.bed
		rm ${sample}_${rep}*.fa ${sample}_${rep}_GC.txt
		rm ${sample}_${rep}*.score ${sample}_${rep}.out
		cd ../
		echo "pipeStage:Calculate splice junction count junction info, END @ $(date)"

	else
		
                echo "pipeStage: Calculate splice junction count junction info, Already Complete @ $(date)"
fi

if [ ! -s ssj/${sample}_${rep}_downEx_ESR.txt ];
	then
	export PERLLIB=$PERLLIB:/data/ohler/Neel/bioTools/perl5:/data/ohler/Neel/bioTools/SROOGLE/web_server
	export PERL5LIB=$PERL5LIB:/data/ohler/Neel/bioTools/perl5:/data/ohler/Neel/bioTools/SROOGLE/web_server
	source ~/.bashrc
	echo "pipeStage: Retrieve ESR info, START @ $(date)"
		cd ssj/
		DIR=$(pwd)
		
		${BEDTOOLS}fastaFromBed -s -name -tab -fi ${FASTA} -bed ${sample}_${rep}_upEx.bed -fo ${sample}_${rep}_upEx.tmp
		${BEDTOOLS}fastaFromBed -s -name -tab -fi ${FASTA} -bed ${sample}_${rep}_downEx.bed -fo ${sample}_${rep}_downEx.tmp
		
		cd /data/ohler/Neel/bioTools/SROOGLE
		perl ${SCRIPTS}runSROOGLE.pl ${DIR}/${sample}_${rep}_upEx.tmp  ${sample}_${rep}_up > ${DIR}/${sample}_${rep}_upEx_ESR.txt
		perl ${SCRIPTS}runSROOGLE.pl ${DIR}/${sample}_${rep}_downEx.tmp  ${sample}_${rep}_down > ${DIR}/${sample}_${rep}_downEx_ESR.txt
		rm ${sample}_${rep}_up_tmp.html ${sample}_${rep}_down_tmp.html 
		cd ${DIR}
		rm ${sample}_${rep}_upEx.tmp rm ${sample}_${rep}_downEx.tmp 
		
		cd ../

		echo "pipeStage: Retrieve ESR info, END @ $(date)"

	else
		echo "pipeStage: Retrieve ESR info, Already Complete @ $(date)"
fi

wait

# Clean up big files
if [ -s filtFastq/${sample}_${rep}_rRNA.sam ];
	then
		echo "pipeStage: Clean up directories, START @ $(date)"
		cd bam

		echo "pipeStage: Compressing unmapped reads and begraph files, START @ $(date)"
		tar cvfz ${sample}_${rep}_Unmapped.tar.gz --wildcards ${sample}_${rep}_Unmapped.out.mate*fastq
		rm ${sample}_${rep}_rev.bg
		tar cvfz ${sample}_${rep}_bg.tar.gz --wildcards ${sample}_${rep}*.bg
		rm ${sample}_${rep}*.bg
		rm ${sample}_${rep}_Unmapped.out.mate*
		echo "pipeStage: Compressing unmapped reads and begraph files, END @ $(date)"

		echo "pipeStage: Compressing rRNA reads, START @ $(date)"
		cd ../filtFastq
		${SCRIPTS}sambamba view -f bam -h -t 4 -S ${sample}_${rep}_rRNA.sam -o ${sample}_${rep}_rRNA.bam
		${SCRIPTS}sambamba sort --tmpdir bam/${sample}_${rep}__STARtmp -t 4 -m 8G -o ${sample}_${rep}_rRNA_sorted.bam ${sample}_${rep}_rRNA.bam
		${SCRIPTS}sambamba index -t 4 -o ${sample}_${rep}_rRNA_sorted.bam
		rm ${sample}_${rep}_rRNA.sam
		echo "pipeStage: Compressing rRNA reads, END @ $(date)"
		
		cd ../
		rm -r bam/${sample}_${rep}__STARtmp
		echo "pipeStage: Clean up directories, END @ $(date)"
fi

rm ${sample}_${rep}.strand

date2=$(date +"%s")
diff=$(($date2-$date1))



echo "pairedEnd RNA-Seq pipeline for ${sample}_${rep} is finsihed"
echo "$(($diff / 3600)) hours and $((($diff % 3600)/60)) minutes elapsed."

exit 0
