#!/bin/bash
#
# Prepare STAR alignment index
#
# Requires STAR, pigz
#
############################################################################################
qsub -l walltime=6:0:0 -q default -l select=1:ncpus=24:mem=100gb:scratch_local=200gb -I
### Variables
GENOME=/storage/brno11-elixir/home/vaclavhejret/NGS_analysis/annotation/ensembl94/human/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
GTF=/storage/brno11-elixir/home/vaclavhejret/NGS_analysis/annotation/ensembl94/human/Homo_sapiens.GRCh38.94.sorted.gtf.gz

USERNAME="bartonv"
INDEX_OUTPUT_DIR=/storage/brno2/home/$USERNAME/RNAseq/results/STAR_index # Where you want to copy created index?

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

THREADS=$PBS_NUM_PPN # Number of threads to use

# Binaries
module add star-2.5.2b
STAR=$(which STAR)

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

cd $SCRATCH/

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

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

GENOME=${GENOME%.gz*}
GTF=${GTF%.gz*}
GEN_DIR=${GENOME%.fa*}

############################################################################################
### Indexing - RUN ONLY ONCE FOR ONE GENOME AND STAR VERSION AND READ LENGTH 
mkdir $SCRATCH/$GEN_DIR

$STAR --runMode genomeGenerate --runThreadN $THREADS --genomeDir $SCRATCH/$GEN_DIR --genomeFastaFiles $GENOME \
--sjdbGTFfile $GTF --sjdbOverhang $RD_LENGTH # --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

############################################################################################
### Copying results
mkdir -p $INDEX_OUTPUT_DIR
cp -r $SCRATCH/$GEN_DIR $INDEX_OUTPUT_DIR/

rm -r $SCRATCH/*
