#$ -N Arima-AJ
#$ -cwd
#$ -o ../logs/$JOB_NAME.$JOB_ID.log
#$ -e ../logs/$JOB_NAME.$JOB_ID.err
#$ -pe smp 56
#$ -l h_vmem=4.4g

###################################################3
DIR=/mnt/data/arslan/Teas/01_Assembly/Arima-Yahs/Output-AJ
READ1=/mnt/data/arslan/Teas/00_Reads_analysis/AJ-HiC_f1.fq.gz
READ2=/mnt/data/arslan/Teas/00_Reads_analysis/AJ-HiC_r2.fq.gz
PICARD=/mnt/bin/picard/v3.0.0/picard/build/libs/picard.jar
CPU=$NSLOTS
MAPQ_FILTER=10
LABEL='AJ_genome'

export PATH=/mnt/bin/bwa/bwa-0.7.17/:$PATH
export PATH=/mnt/bin/samtools/samtools-1.16.1/bin/:$PATH
export PATH=/mnt/data/arslan/tool/mapping_pipeline/:$PATH
conda activate /mnt/bin/picard/v3.0.0/env

mkdir -p $DIR
cd $DIR

ln -s /mnt/data/arslan/Teas/01_Assembly/Purge/AJ_purge/AJ.curated.fasta AJ.curated.fasta

echo "##### Making reference genome Index #####"

bwa index -a bwtsw -p AJ AJ.curated.fasta

samtools faidx AJ.curated.fasta

echo "### Step 1.A: FASTQ to BAM (1st)"
bwa mem -t ${CPU} AJ ${READ1} | samtools view -@ $CPU -Sb - > AJ_f1.bam

echo "### Step 2.A: Filter 5-prime end (1st)"
samtools view -h AJ_f1.bam | filter_five_end.pl  | samtools view -Sb - > AJ_f1.filter.bam

echo "### Step 1.B: FASTQ to BAM (2nd)"
bwa mem -t ${CPU} AJ ${READ2} | samtools view -@ $CPU -Sb - > AJ_r2.bam

echo "### Step 2.B: Filter 5-prime end (2nd)"
samtools view -@ ${NSLOTS} -h AJ_r2.bam | filter_five_end.pl | samtools view -@ {NSLOTS} -Sb - > AJ_r2.filter.bam

echo "### Step 3A: Pair reads & mapping quality filter"
two_read_bam_combiner.pl AJ_f1.filter.bam AJ_r2.filter.bam samtools $MAPQ_FILTER | samtools view -bS -t AJ.curated.fasta.fai - | samtools sort -@ $CPU -o AJ.combine.bam -

echo "### Step 3.B: Add read group"
java -Xmx50G -Djava.io.tmpdir=temp/ -jar $PICARD AddOrReplaceReadGroups INPUT=AJ.combine.bam OUTPUT=AJ.combine.RG.bam ID=AJ LB=AJ SM=${LABEL} PL=ILLUMINA PU=none

echo "### Step 4: Mark duplicates"
java -Xmx50G -XX:-UseGCOverheadLimit -Djava.io.tmpdir=temp/ -jar $PICARD MarkDuplicates INPUT=AJ.combine.RG.bam OUTPUT=AJ.combine.markdup.bam METRICS_FILE=AJ.combine.markdup.txt TMP_DIR=tmp ASSUME_SORTED=TRUE VALIDATION_STRINGENCY=LENIENT REMOVE_DUPLICATES=TRUE

samtools index AJ.combine.RG.bam
samtools index AJ.combine.markdup.bam

get_stats.pl Tea2.combine.RG.bam > Tea2.combine.RG.bam.stats

echo "Finished Mapping Pipeline through Duplicate Removal"
