#!/bin/bash
#PBS -N 04_alignment_STAR_PE
#PBS -l select=1:ncpus=12:mem=150gb:scratch_local=40gb
#PBS -l walltime=12:00:00 
#PBS -m ae
# The 4 lines above are options for scheduling system: job will run 1 hour at maximum, 1 machine with 4 processors + 4gb RAM memory + 10gb scratch memory are requested, email notification will be sent when the job aborts (a) or ends (e)

#
############################################################################################


### Variables

USERNAME="bartonv"
cd /storage/brno2/home/${USERNAME}/
THREADS=$PBS_NUM_PPN
INPUT_DIR=/storage/brno2/home/${USERNAME}/RNASeq/data/preprocessed
OUTPUT_DIR=/storage/brno2/home/${USERNAME}/RNASeq/data/mapped
OUTPUT_DIR_QC=/storage/brno2/home/${USERNAME}/RNASeq/results/qc/mapped


RD_LENGTH=75 # Read length from the sequencing. Illumina sometimes reports N+1 http://seqanswers.com/forums/archive/index.php/t-31154.html; in case you change this value uncomment next line as well


# Binaries
module add star

module add python27-modules-gcc # multiqc

module add samtools


####################################################################################################
### Genome and annotation preparation
GENOME=/storage/brno2/home/$USERNAME/RNASeq/results/STAR_index
GTF=/storage/brno2/home/${USERNAME}/RNASeq/data/genomic.gtf


####################################################################################################
### Alignment
for file in $INPUT_DIR/*_1.fastq.gz; do
	file_F=$(basename $file)
	file_R=${file_F%"_1.fastq.gz"}_2.fastq.gz
	
	  STAR --runThreadN $THREADS \
	    --genomeDir $GENOME \
	    --readFilesIn $INPUT_DIR/$file_F $INPUT_DIR/$file_R \
	    --readFilesCommand zcat \
	    --outFileNamePrefix $OUTPUT_DIR/${file_F%"_1.fastq.gz"}. \
	    --twopassMode Basic \
	    --outFilterType BySJout \
	    --outFilterMultimapNmax 20 \
	    --alignSJoverhangMin 8 \
	    --alignSJDBoverhangMin 1 \
	    --outFilterMismatchNmax 999 \
	    --outFilterMismatchNoverLmax 0.1 \
	    --alignIntronMin 20 \
	    --alignIntronMax 1000000 \
	    --alignMatesGapMax 1000000 \
	    --quantMode TranscriptomeSAM \
	    --outSAMattributes All \
	    --outSAMunmapped Within \
	    --outSAMtype BAM SortedByCoordinate

	echo "Done processing PE reads $file_F and $file_R - alignment"
done

####################################################################################################
### Indexing
for file in $OUTPUT_DIR/*.bam; do
	samtools index $file

done

####################################################################################################
### QC

multiqc -o $OUTPUT_DIR_QC $OUTPUT_DIR

