#!/bin/bash
#qsub -l walltime=5:0:0 -q default -l select=1:ncpus=1:mem=30gb:scratch_local=200gb
# QC of RNA-Seq experiment (after alignment)
#
# Requires Picard (java), UCSC scripts, multiqc
#
############################################################################################
### Variables
#please change these paths to your directories
INPUT_DIR=/storage/brno11-elixir/home/vaclavhejret/Bi5444/RNA-Seq/results/alignment/genome # Pointing to genome/ mapping directory

mkdir ~/alignment_qc
OUTPUT_DIR=~/alignment_qc

GTF=/storage/brno11-elixir/home/vaclavhejret/NGS_analysis/annotation/ensembl94/human/Homo_sapiens.GRCh38.94.sorted.gtf.gz 
#can copy the folder with UCSC tools into your folder
UCSC_SCRIPTS=/storage/brno11-elixir/home/vaclavhejret/NGS_analysis/tools/ucsc # http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/
APPENDIX=".sortedByCoord.out.bam"
THREADS=$PBS_NUM_PPN

#initiate necessary tools
module add RSeQC-2.6.1
module add picard-2.8.1
PICARD_RUN="java -jar $PICARD281"
# Check the tools versions
which $PICARD_RUN

####################################################################################################
### Copying inputs and preparing

cd $SCRATCH/
cp $GTF $SCRATCH/
#copy just one sample bam file as an example
cp $INPUT_DIR/*$APPENDIX $SCRATCH/

GTF=$(basename $GTF)
unpigz -p $THREADS $GTF
GTF=${GTF%.gz}

### Prepare QC input bed12 from gtf for RSeQC; https://gist.github.com/gireeshkbogu/f478ad8495dca56545746cd391615b93
#no direct coversion beween gtf and bed12, need more steps
#convert gtf to genepred format 
$UCSC_SCRIPTS/gtfToGenePred $GTF ${GTF%.*}.genePred
#convert genepred to bed12
$UCSC_SCRIPTS/genePredToBed ${GTF%.*}.genePred ${GTF%.*}.bed12
#sorting of bed12 by chromosome name, and chromosome start position 
sort -k1,1 -k2,2n ${GTF%.*}.bed12 > tmp.bed12
mv tmp.bed12 Homo_sapiens.GRCh38.94.sorted.bed12

# Prepare refFlat from gtf for Picard CollectRnaSeqMetrics
#convert gtf to extended genePred 
$UCSC_SCRIPTS/gtfToGenePred -genePredExt $GTF ${GTF%.*}.refFlat.txt.tmp
#using 12th column of extended genePred (gene name) as a first column to create refFlat
paste <(cut -f 12 ${GTF%.*}.refFlat.txt.tmp) <(cut -f 1-10 ${GTF%.*}.refFlat.txt.tmp) > Homo_sapiens.GRCh38.94.sorted.refFlat

REF_FLAT=${SCRATCH}/Homo_sapiens.GRCh38.94.sorted.refFlat
REF_BED=${SCRATCH}/Homo_sapiens.GRCh38.94.sorted.bed12

####################################################################################################
### Alignment QC
### QC by Picard
mkdir $SCRATCH/picard

# MINIMUM_LENGTH When calculating coverage based values (e.g. CV of coverage) only use transcripts of this length or greater.
# VALIDATION_STRINGENCY if any problem within the anotation occurs, print just warning and continue
for i in *$APPENDIX
do
	$PICARD_RUN CollectRnaSeqMetrics \
	I=$i \
	O=$SCRATCH/picard/${i%.*}.output.RNA_Metrics.txt \
	REF_FLAT=$REF_FLAT \
	STRAND=NONE \
	MINIMUM_LENGTH=200 \
	CHART_OUTPUT=${i%.*}.npc.pdf \
	VALIDATION_STRINGENCY=LENIENT 
done

mv $SCRATCH/*.pdf $SCRATCH/picard

### QC by RSeQC
mkdir $SCRATCH/RSeQC

for i in *$APPENDIX
do
	echo "Processing file $i"

	# Distribution of reads along the genome
	echo "Read distribution"
	read_distribution.py -i $i -r $REF_BED  > $SCRATCH/RSeQC/${i%.bam*}.read_distribution.txt

	# Check strandeness of the experiment
	echo "Infer experiment"
	infer_experiment.py -r $REF_BED -i $i > $SCRATCH/RSeQC/${i%.bam*}.infer_experiment.txt
done

#there is conflict within python dependencies on metacentrum between RSeQC and multiqc, we need to first remove RSeQC module and then load module for multiqc
module remove RSeQC-2.6.1
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)

#run multiqc on output files from Picard CollectRnaSeqMetrics and RSeQC
$MULTIQC -o $SCRATCH/alignment_qc $SCRATCH/RSeQC $SCRATCH/picard/ 

############################################################################################
### Copying results
cp -r $SCRATCH/picard $OUTPUT_DIR/
cp -r $SCRATCH/RSeQC $OUTPUT_DIR/
cp -r $SCRATCH/alignment_qc $OUTPUT_DIR/

rm -r $SCRATCH/*
