#!/bin/bash

# Preprocessing - RNA-Seq PE
#
# Requires FastQC, Trimmomatic (java), multiQC, seqtk, fastp
#
qsub -l walltime=2:0:0 -q default -l select=1:ncpus=6:mem=30gb:scratch_local=150gb -I
#

############################################################################################
### Variables
USERNAME="bartonv"
INPUT_DIR=/storage/brno2/home/$USERNAME/RNAseq/data/raw

ADAPTERS=/storage/brno11-elixir/home/vaclavhejret/Bi5444/RNA-Seq/results/qc/minion/adapters_seq.txt

mkdir -p $SCRATCH/RNAseq/data/preprocessed
#mkdir -p $SCRATCH/RNAseq/results/qc


 # Adapters to trim. Fasta format where /1 or /2 at the end \
#of sequence name signalizes which read is scanned for the adapter. If nothing is specified both reads are checked. You can also use _rc to scan reverse-complement. Please see Trimmomatic manual http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf

APPENDIX=".fastq.gz"
APPENDIX1="_1.fastq.gz"
APPENDIX2="_2.fastq.gz"

PHRED_TRIM=5 # Trim the 3' end of read if four consequent bases have average base quality smaller than this value
LEN_FILTER=30 # Filter sequences shorter than this value after quality trimming
RD_LENGTH=70 # Read length from the sequencing OR trimming down to this length. Illumina sometimes reports N+1 http://seqanswers.com/forums/archive/index.php/t-31154.html

THREADS=$PBS_NUM_PPN # Number of threads to use

module add fastQC-0.11.5
FASTQC=$(which fastqc)

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)

module add jdk-8
module add trimmomatic-0.36
TRIMMOMATIC="java -jar /storage/brno2/home/ebudin/Trimmomatic-0.39/trimmomatic-0.39.jar"


############################################################################################
### Copy inputs
cd $SCRATCH/
cp $INPUT_DIR/*$APPENDIX $SCRATCH/
cp $ADAPTERS $SCRATCH/


ADAPTERS=$(basename $ADAPTERS)

############################################################################################
### Preprocess and qc
mkdir -p $SCRATCH/RNAseq/results/trimmomatic

for i in *$APPENDIX1
do
	READ_FOR=$i
	READ_REV=${i%$APPENDIX1*}$APPENDIX2

	echo "Now I am processing PE reads $READ_FOR and $READ_REV - Trimming"

	$TRIMMOMATIC PE -threads $THREADS -phred33 \
	$READ_FOR $READ_REV \
	$SCRATCH/${READ_FOR%$APPENDIX}_trim.fastq.gz $SCRATCH/${READ_FOR%$APPENDIX}_unpaired.fastq.gz \
	$SCRATCH/${READ_REV%$APPENDIX}_trim.fastq.gz $SCRATCH/${READ_REV%$APPENDIX}_unpaired.fastq.gz \
	ILLUMINACLIP:${ADAPTERS}:2:30:5:3:true CROP:$RD_LENGTH LEADING:3 TRAILING:3 SLIDINGWINDOW:4:$PHRED_TRIM MINLEN:$LEN_FILTER &> $SCRATCH/RNAseq/results/trimmomatic/${READ_FOR%$APPENDIX1}_trim.log

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

#fastqc on preprocessed reads
$FASTQC --outdir $SCRATCH/RNAseq/data/preprocessed --format fastq --threads $THREADS $SCRATCH/*_trim.fastq.gz
#multiqc on fastqc files
$MULTIQC -o $SCRATCH/RNAseq/data/preprocessed $SCRATCH/RNAseq/data/preprocessed/
#mutliqc on trimmomatic logs
$MULTIQC -o $SCRATCH/RNAseq/results/trimmomatic $SCRATCH/RNAseq/results/trimmomatic/

############################################################################################
### Copy results

OUTPUT_DIR=/storage/brno2/home/$USERNAME/RNAseq/data/preprocessed
OUTPUT_DIR_QC=/storage/brno2/home/$USERNAME/RNAseq/results/qc

mkdir -p $OUTPUT_DIR
mkdir -p $OUTPUT_DIR_QC/fastqc

cp -r $SCRATCH/RNAseq/data/preprocessed $OUTPUT_DIR_QC/fastqc/ 
cp -r $SCRATCH/RNAseq/results/trimmomatic $OUTPUT_DIR_QC/
cp $SCRATCH/*_trim.fastq.gz $OUTPUT_DIR/

rm -r $SCRATCH/*
