#!/bin/bash
#qsub -l walltime=24:0:0 -q default -l select=1:ncpus=12:mem=40gb:scratch_local=200gb -I
#
# Alignment - PE RNA-Seq
# To increase sensitivity you might try to add --seedSearchStartLmax 30
# If you get many "Too short" unmapped reads you might adjust --outFilterScoreMinOverLread and --outFilterMatchNminOverLread and go lower (0.66->0.33 or even to 0)
#
# Requires STAR, pigz, multiQC
#
############################################################################################

### Variables

INPUT_DIR=/storage/brno11-elixir/home/vaclavhejret/Bi5444/RNA-Seq/data/preprocessed

mkdir -p $SCRATCH/mapping_result/results
mkdir -p $SCRATCH/mapping_result/results

USERNAME="bartonv"
OUTPUT_DIR=/storage/brno2/home/$USERNAME/RNAseq/mapping_result/results
OUTPUT_DIR_QC=/storage/brno2/home/$USERNAME/RNAseq/qc_result


APPENDIX1="1_trim.fastq.gz"
APPENDIX2="2_trim.fastq.gz"
APPENDIX="_trim.fastq.gz"

# Read length for sjdbOverhang; --sjdbOverhang 100 should work fine for most of the data, but more specific setting based on the real read length should be more sensitive https://groups.google.com/forum/#!msg/rna-star/h9oh10UlvhI/BfSPGivUHmsJ
RD_LENGTH=100 # 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

# Genome and annotation
GENOME=/storage/brno11-elixir/home/vaclavhejret/NGS_analysis/annotation/ensembl94/human/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
GENOME_INDEX=/storage/brno11-elixir/home/vaclavhejret/NGS_analysis/annotation/ensembl94/human/STAR_index/Homo_sapiens.GRCh38.dna.primary_assembly # Default genome index with --sjdbOverhang 100 - STAR default
GTF=/storage/brno11-elixir/home/vaclavhejret/NGS_analysis/annotation/ensembl94/human/Homo_sapiens.GRCh38.94.sorted.gtf.gz

THREADS=$PBS_NUM_PPN

# Binaries
module add star-2.5.2b
STAR=$(which STAR)
module add samtools-1.4
SAMTOOLS=$(which samtools)
module add python27-modules-gcc
#system variables needed to properly run MultiQC tool
export LC_ALL=C.UTF-8
export LANG=C.UTF-8
MULTIQC=$(which multiqc)

# Check the tools versions
which $STAR
which $SAMTOOLS
which $MULTIQC

####################################################################################################
### Copy inputs
cp $GTF $SCRATCH/
cp $GENOME $SCRATCH/

cp -r $GENOME_INDEX $SCRATCH/

cp $INPUT_DIR/*$APPENDIX $SCRATCH/

cd $SCRATCH/

####################################################################################################
### Genome and annotation preparation
GENOME=$(basename $GENOME)
GTF=$(basename $GTF)

unpigz -p $THREADS $GTF
unpigz -p $THREADS $GENOME

GENOME=${GENOME%.gz*}
GTF=${GTF%.gz*}
GEN_DIR=$(basename $GENOME_INDEX)

####################################################################################################
### Alignment
for i in *$APPENDIX1
do
	READ_FOR=$i
	READ_REV=${READ_FOR%$APPENDIX1*}$APPENDIX2
	OUTNAME=${READ_FOR%$APPENDIX1*}
	echo "Now I am processing PE reads $READ_FOR and $READ_REV - alignment"

	$STAR --runMode alignReads --runThreadN $THREADS --genomeDir $GEN_DIR \
	--readFilesIn $READ_FOR $READ_REV \
	--readFilesCommand zcat --sjdbOverhang $RD_LENGTH --sjdbGTFfile $GTF \
	--outSAMattributes All \
	--outFileNamePrefix ${OUTNAME}.${GTF%.*} \
	--outSAMattrRGline ID:$OUTNAME PL:Illumina PU:$OUTNAME SM:$OUTNAME \
	--outSAMheaderHD @HD VN:1.4 SO:coordinate \
	--outMultimapperOrder Random --outSAMtype BAM SortedByCoordinate

	echo "Done processing PE reads $READ_FOR and $READ_REV - alignment"
done

#creating of index of each aligned bam file, necessary for many tools and visualization
for i in *.bam
do
	${SAMTOOLS} index -@ ${THREADS} ${i}
done

# Move STAR logs
mkdir $SCRATCH/star_log
mv $SCRATCH/*.out $SCRATCH/star_log/

$MULTIQC -o $SCRATCH/star_log $SCRATCH/star_log/ 


####################################################################################################
### Finalize and copy results

mkdir -p $OUTPUT_DIR
mkdir -p $OUTPUT_DIR_QC

mkdir -p $SCRATCH/alignment/genome
mv $SCRATCH/*.bam $SCRATCH/alignment/genome/
mv $SCRATCH/*.bai $SCRATCH/alignment/genome/

rm $SCRATCH/*$APPENDIX
rm $SCRATCH/*.fa.gz
rm $SCRATCH/$GENOME*
rm $SCRATCH/$GTF*
rm -r $SCRATCH/$GEN_DIR
rm -r $SCRATCH/*_STARpass1
rm -r $SCRATCH/*_STARgenome
rm -r $SCRATCH/*_STARtmp

cp -r $SCRATCH/star_log $OUTPUT_DIR_QC/
rm -r $SCRATCH/star_log 
cp -r $SCRATCH/* $OUTPUT_DIR/

rm -r $SCRATCH/*
