#!/bin/bash

# =============================================================================================================================================================
# ==                                                          SLURM parameters                                                                               ==
# =============================================================================================================================================================

# Decide how much computer power you need
# The programs used in this script can only be executed on 1 NODE
# You can use a max of 24 task per node. But the lower the number of tasks requested the shorter the queue when submitting a job
#SBATCH --nodes=1			
#SBATCH --ntasks-per-node=8
#SBATCH --mem=8GB
#SBATCH --time=24:00:00
#SBATCH --job-name=read_start_counts
#SBATCH --output=slurm-read_start_counts-%j.log



# =============================================================================================================================================================
# ==                                                          Load modules                                                                                   ==
# =============================================================================================================================================================

module load SAMtools/1.19.2-GCC-13.2.0

source /userapps/hb-molgen/conda_envs/venvs/PythonTensor/bin/activate

cpu=8
AnneTools=/userapps/hb-molgen/AnneTools



# =============================================================================================================================================================
# ==                                                          Functions                                                                                      ==
# =============================================================================================================================================================

function read_start_counts_forward {
	name=$1
	samtools view $datadir/$name.sorted.bam \
	  | awk '{if (!and($2,16)) print $3, $4, "+"}' \
	  | sort | uniq -c \
	  | awk '{print $2"\t"$3"\t"$4"\t"$1}' \
	  | sort -k1,1 -k2,2n > $resultsdir/$name.forward.read_start_counts.tsv
}

function read_start_counts_reverse {
	# reverse strand is more complex because we need the end of the alignment
	name=$1
	samtools view $datadir/$name.sorted.bam \
		| awk '{
			cigar = $6;
			ref_consumed = 0;
			while (match(cigar, /[0-9]+[MDN=X]/)) {
				ref_consumed += substr(cigar, RSTART, RLENGTH-1);
				cigar = substr(cigar, RSTART+RLENGTH);
			}
			if (and($2,16)) print $3, $4 + ref_consumed - 1, "-";
		}' \
		| sort | uniq -c \
		| awk '{print $2"\t"$3"\t"$4"\t"$1}' \
		| sort -k1,1 -k2,2n > $resultsdir/$name.reverse.read_start_counts.tsv
}

function map5prime {
	name=$1
	
	# map the read start points
	read_start_counts_forward $name
	read_start_counts_reverse $name
	
	# make bedgraph
	awk '{print $1"\t"$2-1"\t"$2"\t"$4}' $resultsdir/$name.forward.read_start_counts.tsv  > $resultsdir/$name.forward.read_start_counts.bedgraph
	awk '{print $1"\t"$2-1"\t"$2"\t"$4}' $resultsdir/$name.reverse.read_start_counts.tsv > $resultsdir/$name.reverse.read_start_counts.bedgraph

	# peak calling for exact 5'end
	python3 $AnneTools/PeakCallingEndSeq.py --strand forward  --basename $resultsdir/$name  
	python3 $AnneTools/PeakCallingEndSeq.py --strand reverse  --basename $resultsdir/$name

	intergenic=/scratch/p127804/PROJECTS/Wilfried/genomes/Bacillus_subtilis_subsp_subtilis_str_168_ASM904v1_pLS20.intergenic.table
	python3 $AnneTools/PeakFilterEndSeq.py --strand forward -a 10 -i $intergenic -p $resultsdir/$name.forward.peaks.tsv -o $resultsdir/$name.filtered.forward.peaks.tsv  
	python3 $AnneTools/PeakFilterEndSeq.py --strand reverse -a 10 -i $intergenic -p $resultsdir/$name.reverse.peaks.tsv -o $resultsdir/$name.filtered.reverse.peaks.tsv  

	cat $resultsdir/$name.filtered.forward.peaks.tsv $resultsdir/$name.filtered.reverse.peaks.tsv > $resultsdir/$name.filtered.peaks.tsv

	fna=/scratch/p127804/PROJECTS/Wilfried/genomes/Bacillus_subtilis_subsp_subtilis_str_168_ASM904v1_pLS20.g2d.fna
	python3 $AnneTools/GetRegionsFromFNA.py --fna $fna --table $resultsdir/$name.filtered.peaks.tsv -o $resultsdir/$name.regions  

	
}	

function FACoPv3 {
	bash /scratch/p127804/NGS/data/FACoPv3/30.FACoPv3_AnnotateGenomes.sh $1
}

# =============================================================================================================================================================
# ==                                                          Main script                                                                                    ==
# =============================================================================================================================================================

# Annotate genome using FACoPv3,, only need once if not already done
# FACoPv3 /scratch/p127804/PROJECTS/Wilfried/genomes/Bacillus_subtilis_subsp_subtilis_str_168_ASM904v1_pLS20.gff

datadir=/scratch/p127804/PROJECTS/Wilfried/202505_023/results/
resultsdir=/scratch/p127804/PROJECTS/Wilfried/202505_023/read_start_counts

mkdir -p $resultsdir

map5prime AJ023_03
#map5prime AJ023_04