#!/bin/bash
#qsub -l walltime=2:0:0 -q default -l select=1:ncpus=6:mem=40gb:scratch_local=200gb -I
#
# Gene/isoform counting - RNA-Seq
#
# Can handle both SE and PE but you should know the strandness of the experiment 
# Input is from STAR mapping to genome
#
# Requires MultiQC, featureCounts
#
############################################################################################
### Variables
INPUT_DIR=/storage/brno11-elixir/home/vaclavhejret/Bi5444/RNA-Seq/results/alignment/genome # Pointing to general mapping directory containing genome/ mapping

OUTPUT_DIR_QC=~/feature_count_qc
OUTPUT_DIR=~/gene_counts

COUNT_OVER="exon" # [exon, three_prime_utr] what feature to use for the summary? For QuantSeq REV it might be 3 UTR ("three_prime_utr" as for Ensembl annotation); for Lexogen REV we should add --readExtension5 100 (or similar) to help in case the 3' UTR annotation is incorrect and not complete ("For some gene annotations it might be an advantage to extend the 3’ UTR annotation further downstream in order to assign the mapped read correctly")
GTF=/storage/brno11-elixir/home/vaclavhejret/NGS_analysis/annotation/ensembl94/human/Homo_sapiens.GRCh38.94.sorted.gtf.gz # Reference GTF (.gz)

THREADS=$PBS_NUM_PPN


# Binaries
module add subread-1.5.2
FEATURE_COUNTS=$(which featureCounts)
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 $FEATURE_COUNTS
which $MULTIQC

############################################################################################
### Copy inputs
cp $GTF $SCRATCH/
cp -r $INPUT_DIR/* $SCRATCH/ 
cd $SCRATCH/

####################################################################################################
### annotation preparation

GTF=$(basename $GTF)
gunzip $GTF
GTF=${GTF%.gz*}


############################################################################################
### Gene counting
# Count reads using featureCounts
mkdir -p $SCRATCH/gene_counts/featureCounts
# -t which features in gtf should be counted
# -p if reads are pair end
# -P heck the fragment length if reads are pair end
# -C not counting chimeric fragments (each part of pair read became from different chromosome)
# -d minimum fragment length
# -D maximum fragment length
# -s 0 unstranded experiment
# -Q minimum mapping quality to count read 

$FEATURE_COUNTS -t $COUNT_OVER -g gene_id -p -P -C -s 0 -T $THREADS \
-Q 0 -d 1 -D 25000 -a $SCRATCH/$GTF -o $SCRATCH/gene_counts/featureCounts/Unstranded.featureCounts *.bam 
mkdir $SCRATCH/featureCounts_gc

#multiqc on logs from feature counts
$MULTIQC -o $SCRATCH/featureCounts_gc $SCRATCH/gene_counts/featureCounts/ 


############################################################################################
### Copying results
mkdir -p $OUTPUT_DIR
mkdir -p $OUTPUT_DIR_QC

cp -r $SCRATCH/gene_counts/* $OUTPUT_DIR/
cp -r $SCRATCH/featureCounts_gc $OUTPUT_DIR_QC/

rm -r $SCRATCH/*
