#$ -S /bin/bash
#$ -N merge_remove-primers_quality-filter_uniques_find-joint
#$ -o /workdir/users/pd378/oilPCR/bulk_process_all_OIL_for_paper/read_process_out.txt
#$ -e /workdir/users/pd378/oilPCR/bulk_process_all_OIL_for_paper/read_process_err.txt
#$ -wd /workdir/users/pd378/oilPCR/bulk_process_all_OIL_for_paper #Your working directory
#$ -pe parenv 1
#$ -l h_vmem=30G
#$ -t 1-46
#$ -q short.q@cbsubrito

#export usearch
export PATH=/programs/usearch11.0.667:$PATH
##Set dirs
WRK=/workdir/users/pd378/oilPCR/bulk_process_all_OIL_for_paper
##Path to the list of file names
LIST=$WRK//file_lists/gene_lists/CTXlist.txt

####MERGE VARIABLES###
MERGE_IN=$WRK/rawdata ## raw reads folder
MERGE_OUT=$WRK/s1_merge_reads ##where to put the merged reads
MIN_MERGE=337       #cutoff for minimum merge length
MAX_MERGE=367       #cutoff for maximum merge length 
MAX_DIFFS=13       #How many difference are allowed for overlap
PCTID=85 
####PRIMER STRIP VARIABLES####  set variables in call
PSTRIP_OUT=$WRK/s2_strip_primers #output for the stripped reads

####QUALITY FILTER VARIABLES####
QFILTER_OUT=$WRK/s3_quality_filter
MAX_ERROR=1.0
####FILTER UNIQUES VARIABLES####
UNIQUE_OUT=$WRK/s4_filter_uniques

####CUTADAPT
TARGET_OUT=$WRK/s5_split_reads/target_reads
RIBO_OUT=$WRK/s5_split_reads/ribo_reads
FUSION="TTCTGGTSACYTAYTTYACCCACAGCMGCCGCGGTAATWC"

##Make the design file
        DESIGN=$(sed -n "${SGE_TASK_ID}p" $LIST)
        NAME=`basename "$DESIGN"`

##############RUN USEARCH MERGE PAIRS ####################
mkdir $MERGE_OUT
mkdir $MERGE_OUT/stats
FORWARD=$MERGE_IN/${NAME}_R1.fastq
REVERSE=$MERGE_IN/${NAME}_R2.fastq
M_OUT=$MERGE_OUT/${NAME}_merged.fq
M_STATS=$MERGE_OUT/stats/${NAME}.txt

usearch -fastq_mergepairs $FORWARD -reverse $REVERSE -fastq_maxdiffs $MAX_DIFFS -fastq_pctid $PCTID -fastq_minmergelen $MIN_MERGE -fastq_maxmergelen $MAX_MERGE -fastqout $M_OUT 2> $M_STATS

echo $NAME merge complete

############ RUN PRIMER TRUNCATE #####################
mkdir $PSTRIP_OUT
mkdir $PSTRIP_OUT/stats
P_OUT=$PSTRIP_OUT/${NAME}_stripped.fq
P_STATS=$PSTRIP_OUT/stats/${NAME}.txt

usearch -fastx_truncate $M_OUT -stripleft 20 -stripright 20 -fastqout $P_OUT 2> $P_STATS

echo $NAME primer strip complete

########### RUN QUALITY FILTER ###################
mkdir $QFILTER_OUT
mkdir $QFILTER_OUT/stats
Q_OUT=$QFILTER_OUT/${NAME}_filtered.fa
Q_STATS=$QFILTER_OUT/stats/${NAME}.txt

usearch -fastq_filter $P_OUT -fastq_maxee $MAX_ERROR -fastaout $Q_OUT 2> $Q_STATS

echo $NAME quality filter compete

########### RUN FILTER UNIQUES ##################
mkdir $UNIQUE_OUT
mkdir $UNIQUE_OUT/stats
U_OUT=$UNIQUE_OUT/${NAME}_unique.fa
U_STATS=$UNIQUE_OUT/stats/${NAME}.txt

usearch -fastx_uniques $Q_OUT -fastaout $U_OUT -relabel Uniq -sizeout 2> $U_STATS

echo $NAME filter uniques complete

############# RUN CUTADAPT #######################
mkdir $WRK/s5_split_reads
mkdir $TARGET_OUT
mkdir $RIBO_OUT
mkdir $TARGET_OUT/stats
mkdir $RIBO_OUT/stats
T_OUT=$TARGET_OUT/${NAME}_targ.fa
R_OUT=$RIBO_OUT/${NAME}_ribo.fa


cutadapt -g $FUSION -o $R_OUT --discard-untrimmed   $U_OUT > $RIBO_OUT/stats/${NAME}_ribo.stats
cutadapt -a $FUSION -o $T_OUT --discard-untrimmed   $U_OUT > $TARGET_OUT/stats/${NAME}_targ.stats