#!/bin/bash
# RNA-Seq Pipeline
# 2023 Februari

# =============================================================================================================================================================
# ==                                                          Settings for SLURM                                                                             ==
# =============================================================================================================================================================

#SBATCH --nodes=1			
#SBATCH --ntasks-per-node=16
#SBATCH --mem=32GB
#SBATCH --job-name=RNAseq
#SBATCH --time=24:00:00

# Personal stuff; Change this to your personal settings 
#SBATCH --job-name=Bowtie2_Diego
#SBATCH --mail-type=END
#SBATCH --mail-user=anne.de.jong@.rug.nl
#SBATCH --output=slurm-Bowtie2-%j.log


# =============================================================================================================================================================
# ==                                                          Load modules                                                                                   ==
# =============================================================================================================================================================
module load Biopython/1.81-foss-2022b
module load SAMtools/1.16.1-GCC-11.3.0
module load BEDTools/2.30.0-GCC-11.3.0
module load Bowtie2/2.4.5-GCC-11.3.0

# setting global parameters
cpu=$SLURM_JOB_CPUS_PER_NODE
cpu=16
SOFTWAREDIR=/home1/p127804


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

function my_bowtie_SE {  # Map the reads (from the fastq file) to the reference genome ( .fna ) this will generate a SAM file (Sequence Alignment Map)
	name=$1
	forward=$2
	printf "\n\n========= Bowtie2 SE: "$name" =================================\n" >> $resultsdir/02.$name.bowtie2.log
	(bowtie2 --local -U $datadir/$forward -p $cpu -S $resultsdir/$name.sam -x $genome) 2>>$resultsdir/02.$name.bowtie2.log 
	SAM_2_BAM $name
}

function SAM_2_BAM {  # convert SAM to sorted BAM (binairy SAM) and remove unmapped reads with option -F 4
	name=$1
	printf "\n========== Convert "$name" ======================================\n"	
	samtools view -S -b -F 4 -o $resultsdir/$name.bam $resultsdir/$name.sam
	samtools sort -@ $cpu $resultsdir/$name.bam > $resultsdir/$name.sorted.bam
	samtools index $resultsdir/$name.sorted.bam
	rm $resultsdir/$name.sam  # save space by removing the SAM file 
	rm $resultsdir/$name.bam  # save space by removing the non sorted BAM file                                                                          # 3) save space by removing the SAM file 
}

function my_mkdir {
	if [ ! -d $1 ]; then
		mkdir $1 ;
	fi
}


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


# 1. FOLDERS
	# NOTE change the folder names to your own SCRATCHDIR and resultsdir
	# datadir=/data/pg-molgen/rnaseq/20230211_008/fastq
	resultsdir=/scratch/p127804/PROJECTS/Wilfried/202310_012/results
	datadir=/scratch/hb-molgen/NGSdata/20232410_012/WilfriedMeijer
	my_mkdir $resultsdir


# 2. 
	# NOTE to get all filenames use :  ls -aR >TOC.txt
	# The main function needs 4 parameters; Sample_name   Foward_Reads   Reverse_Reads   Genome_base_Name
		# 1) Sample Name: Choose a short name. This prefix will be used to name the result files
		# 2) The genome filename is normally stored in the $genome variable
		# 3) Foward_Reads:  Name of the fastq file containing the forward reads
		# 4) Reverse_Reads: Name of the fastq file containing the reverse reads, but only when Paired End (PE) sequencing is done 


# 3. For Transcriptomics usually Single End (SE) sequencing is done
	
genome=/scratch/hb-molgen/genomes/Bacillus_subtilis_subsp_subtilis_str_168_ASM904v1_pLS20

my_bowtie_SE WM012_32 WM012_32_S50_R1_001.fastq.gz
my_bowtie_SE WM012_38 WM012_38_S51_R1_001.fastq.gz
my_bowtie_SE WM012_46 WM012_46_S52_R1_001.fastq.gz
my_bowtie_SE WM012_48 WM012_48_S53_R1_001.fastq.gz
my_bowtie_SE WM012_51 WM012_51_S54_R1_001.fastq.gz
my_bowtie_SE WM012_53 WM012_53_S55_R1_001.fastq.gz
my_bowtie_SE WM012_63 WM012_63_S56_R1_001.fastq.gz
my_bowtie_SE WM012_64 WM012_64_S57_R1_001.fastq.gz
my_bowtie_SE WM012_65 WM012_65_S58_R1_001.fastq.gz
my_bowtie_SE WM012_66 WM012_66_S59_R1_001.fastq.gz
my_bowtie_SE WM012_67 WM012_67_S60_R1_001.fastq.gz
my_bowtie_SE WM012_68 WM012_68_S61_R1_001.fastq.gz
my_bowtie_SE WM012_69 WM012_69_S62_R1_001.fastq.gz
my_bowtie_SE WM012_70 WM012_70_S63_R1_001.fastq.gz
my_bowtie_SE WM012_71 WM012_71_S64_R1_001.fastq.gz
my_bowtie_SE WM012_72 WM012_72_S65_R1_001.fastq.gz



# Make a report summary
cat $resultsdir/*.bowtie2.log > $resultsdir/03.bowtie2_all.log
python3 $SOFTWAREDIR/AnneTools/Bowtie2LogTable.py -query $resultsdir/03.bowtie2_all.log > $resultsdir/03.bowtie2_table.txt