#!/bin/bash

#gunzip *.gz

for i in {1,2,4,5,6,7,8}
do

#extract UMI
umi_tools extract --stdin=S${i}_S${i}_L001_R1_001.fastq --bc-pattern=NNNNNNNN --log=processed.log --stdout S${i}_R1_rmUMI.fastq

#put read 2 next to read 1 for separating them according to tags
paste S${i}_R1_rmUMI.fastq S${i}_S${i}_L001_R2_001.fastq > S${i}.fastq 

#separate reads according to tags and split them to read 1 and read 2
#H-tag
cat S${i}.fastq | grep -A 2 -B 1 "^CTAGGACGTACTC" | sed '/^--$/d' | cut -f1 > Htag_S${i}_R1.fastq    
cat S${i}.fastq | grep -A 2 -B 1 "^CTAGGACGTACTC" | sed '/^--$/d' | cut -f2 > Htag_S${i}_R2.fastq  

#U-tag
cat S${i}.fastq | grep -A 2 -B 1 "^TTAGCGACACGA" | sed '/^--$/d' | cut -f2 > Utag_S${i}_R2.fastq
cat S${i}.fastq | grep -A 2 -B 1 "^TTAGCGACACGA" | sed '/^--$/d' | cut -f1 > Utag_S${i}_R1.fastq

#M-tag
cat S${i}.fastq | grep -A 2 -B 1 "^AGCAGATGACGT" | sed '/^--$/d' | cut -f1 > Mtag_S${i}_R1.fastq
cat S${i}.fastq | grep -A 2 -B 1 "^AGCAGATGACGT" | sed '/^--$/d' | cut -f2 > Mtag_S${i}_R2.fastq

#Cutadapt to remove tags and adapter. [-g is for adapter on 5' end; -b is anywhere along the read; -n 2 to indicate must look for both adapters in each read not just stop after one] 
#remove the original 5' tag as well as 3' tag from NlaIII digestion from read 1
cutadapt -g "^TTAGCGACACGA;min_overlap=6" -b "CATGAGATCGGAAGAGC;min_overlap=6" -n 2 -o Utag_S${i}_R1_adaptcut.fastq Utag_S${i}_R1.fastq 
cutadapt -g "^AGCAGATGACGT;min_overlap=6" -b "CATGAGATCGGAAGAGC;min_overlap=6" -n 2 -o Mtag_S${i}_R1_adaptcut.fastq Mtag_S${i}_R1.fastq 
cutadapt -g "^CTAGGACGTACTC;min_overlap=6" -b "CATGAGATCGGAAGAGC;min_overlap=6" -n 2 -o Htag_S${i}_R1_adaptcut.fastq Htag_S${i}_R1.fastq

#remove reverse complement of Htag, Mtag and Utag that may enter read 2 for short reads
cutadapt -b "TCGTGTCGCTAA;min_overlap=6" -b "ACGTCATCTGCT;min_overlap=6" -b "GAGTACGTCCTAG;min_overlap=6" -o Utag_S${i}_R2_adaptcut.fastq Utag_S${i}_R2.fastq 
cutadapt -b "TCGTGTCGCTAA;min_overlap=6" -b "ACGTCATCTGCT;min_overlap=6" -b "GAGTACGTCCTAG;min_overlap=6" -o Mtag_S${i}_R2_adaptcut.fastq Mtag_S${i}_R2.fastq
cutadapt -b "TCGTGTCGCTAA;min_overlap=6" -b "ACGTCATCTGCT;min_overlap=6" -b "GAGTACGTCCTAG;min_overlap=6" -o Htag_S${i}_R2_adaptcut.fastq Htag_S${i}_R2.fastq

#perform alignment with bowtie2 [-x refer to reference; -q in samtolls refers to remove reads below quality score of two- multi mappers(aligned to more than one place); -b is to convert the default SAM file to BAM file for next step like dedup; finally the file needs to be sorted in order to index]
bowtie2 -p 13 -x ~/ramya/Ref/Mus_musculus_mm10/NCBI/GRCm38/Sequence/Bowtie2Index/genome -1 Utag_S${i}_R1_adaptcut.fastq -2 Utag_S${i}_R2_adaptcut.fastq  |samtools view -q 2 -b -| samtools sort -o  UtagS${i}.bam	 
bowtie2 -p 13 -x ~/ramya/Ref/Mus_musculus_mm10/NCBI/GRCm38/Sequence/Bowtie2Index/genome -1 Mtag_S${i}_R1_adaptcut.fastq -2 Mtag_S${i}_R2_adaptcut.fastq  |samtools view -q 2 -b -| samtools sort -o  MtagS${i}.bam		
bowtie2 -p 13 -x ~/ramya/Ref/Mus_musculus_mm10/NCBI/GRCm38/Sequence/Bowtie2Index/genome -1 Htag_S${i}_R1_adaptcut.fastq -2 Htag_S${i}_R2_adaptcut.fastq  |samtools view -q 2 -b -| samtools sort -o  HtagS${i}.bam		

#index alignments separately
samtools index UtagS${i}.bam 
samtools index MtagS${i}.bam
samtools index HtagS${i}.bam 

#perform deduplication using umi_tools
umi_tools dedup --stdin=UtagS${i}.bam --log=UtagS${i}logfile > UtagS${i}_dedup.bam
umi_tools dedup --stdin=MtagS${i}.bam --log=MtagS${i}logfile > MtagS${i}_dedup.bam
umi_tools dedup --stdin=HtagS${i}.bam --log=HtagS${i}logfile > HtagS${i}_dedup.bam

#convert bam file to bedgraph to measure coverage
bedtools genomecov -ibam UtagS${i}_dedup.bam -bg > UtagS${i}_dedup.bedGraph
bedtools genomecov -ibam MtagS${i}_dedup.bam -bg > MtagS${i}_dedup.bedGraph
bedtools genomecov -ibam HtagS${i}_dedup.bam -bg > HtagS${i}_dedup.bedGraph

#update format of bedgraph to include "chr" in front of chromosome numbers
cat UtagS${i}_dedup.bedGraph | awk '{print "chr"$1 "\t" $2 "\t"  $3 "\t"  $4}' > UtagS${i}_dedup2.bedGraph 
cat MtagS${i}_dedup.bedGraph | awk '{print "chr"$1 "\t" $2 "\t"  $3 "\t"  $4}' > MtagS${i}_dedup2.bedGraph 
cat HtagS${i}_dedup.bedGraph | awk '{print "chr"$1 "\t" $2 "\t"  $3 "\t"  $4}' > HtagS${i}_dedup2.bedGraph

#Only considers "assayable" CCGG sites, which are sequences in the reference genome that are flanked by CCGG and CATG
bedtools intersect -wb -a ~/ramya/Regions/mm10_CCGG_Assyable_Gcord -b UtagS${i}_dedup2.bedGraph -loj| cut -f1,2,3,7 >S${i}_Assayable_Utag
bedtools intersect -wb -a ~/ramya/Regions/mm10_CCGG_Assyable_Gcord -b MtagS${i}_dedup2.bedGraph -loj| cut -f7 > S${i}_Assayable_Mtag
bedtools intersect -wb -a ~/ramya/Regions/mm10_CCGG_Assyable_Gcord -b HtagS${i}_dedup2.bedGraph -loj| cut -f7 > S${i}_Assayable_Htag

#compile the numbers of U-tag, M-tag and H-tags at each assayable site
paste S${i}_Assayable_Utag S${i}_Assayable_Mtag S${i}_Assayable_Htag | sed 's/\./0/g' >S${i}_AllAssaybletags

#calculate the unmethylated, methylated and hydroxymethylated ratio at each site, and total coverage ($7)
cat S${i}_AllAssaybletags | awk '{if ($4+$5+$6>0) {print $1 "\t" $2 "\t" $3 "\t" $4/($4+$5+$6) "\t" $5/($4+$5+$6) "\t" $6/($4+$5+$6) "\t" $4+$5+$6} else {print $1 "\t" $2 "\t" $3 "\t" "NA" "\t" "NA" "\t" "NA" "\t" $4+$5+$6 }}' > S${i}_ratiotemp.txt  

#calculate 5hmC/5mC ratio ($7)
cat S${i}_ratiotemp.txt | awk '{if ($5!="NA" && $5!=0) {print $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6 "\t" $6/$5 "\t" $7} else {print $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6 "\t" "NA" "\t" $7}}' > S${i}_ratiotemp2.txt  

#calculate 5mC/(5C+5mC) ($7)
cat S${i}_ratiotemp2.txt | awk '{if ($4!="NA" && $4+$5!=0) {print $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6 "\t" $7 "\t" $5/($4+$5) "\t" $8} else {print $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6 "\t" $7" \t" "NA" "\t" $8}}' > S${i}_ratio.txt

#cleanup
rm Htag_S${i}_R1.fastq Htag_S${i}_R2.fastq Utag_S${i}_R1.fastq Utag_S${i}_R2.fastq  Mtag_S${i}_R1.fastq  Mtag_S${i}_R2.fastq
echo ${i}

done
