#PBS -q xl230g9q
#PBS -l select=1:ncpus=5:mem=20g
#PBS -l walltime=3:00:00
#PBS -l cput=30:00:00

source activate nn_Polina_tf2;
export PATH=$PATH:~/programms/bwa-0.7.17;
WORK_DIR=/mnt/storage/home/psbelokopytova/scratch/202312241147Polina/loxP_project/2023-05-18
cd $WORK_DIR; mkdir align; 
cd align;
# Declare a string array with type

#declare -a files=("sum_wt1" "sum_wt2" "sum_wt3")
declare -a files=("sum_rear1")

for file in "${files[@]}"
do
echo ${file}
mkdir ${file}


## align paired reads
bwa mem -t 5 -M /mnt/storage/home/psbelokopytova/insulatory_index/bds_pipeline_genome_data/hg38/bwa_index/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta $WORK_DIR/cutadapt/${file}_1.fq $WORK_DIR/cutadapt/${file}_2.fq | samtools view -bS -@ 5 | samtools sort > $WORK_DIR/align/${file}/${file}.sorted.bam
## get index for bam file
samtools index $WORK_DIR/align/${file}/${file}.sorted.bam > $WORK_DIR/align/${file}/${file}.sorted.bam.bai
## get reads avove mapq30  
samtools view -q 30 -bS $WORK_DIR/align/${file}/${file}.sorted.bam > $WORK_DIR/align/${file}/${file}_aboveMapq30.sorted.bam
samtools view $WORK_DIR/align/${file}/${file}_aboveMapq30.sorted.bam > $WORK_DIR/align/${file}/${file}_aboveMapq30.sorted.sam
## drop duplicates
python ~/scratch/202312241147Polina/loxP_project/drop_duplicates_by_read_segments.py $WORK_DIR/align/${file}/${file}_aboveMapq30.sorted.sam $WORK_DIR/align/${file}/${file}.nodups_reads2.txt
java -jar ~/programms/picard.jar FilterSamReads I=$WORK_DIR/align/${file}/${file}_aboveMapq30.sorted.bam O=$WORK_DIR/align/${file}/${file}_aboveMapq30_nodup2.sorted.bam READ_LIST_FILE=$WORK_DIR/align/${file}/${file}.nodups_reads2.txt FILTER=includeReadList
## get primary reads
samtools view -F 256 -bS $WORK_DIR/align/${file}/${file}_aboveMapq30_nodup2.sorted.bam > $WORK_DIR/align/${file}/${file}_aboveMapq30_nodup2_prim.sorted.bam 
samtools index $WORK_DIR/align/${file}/${file}_aboveMapq30_nodup2_prim.sorted.bam > $WORK_DIR/align/${file}/${file}_aboveMapq30_nodup2_prim.sorted.bam.bai

## create coverage file
~/programms/bedtools2/bin/bedtools genomecov -ibam $WORK_DIR/align/${file}/${file}_aboveMapq30_nodup2_prim.sorted.bam -bg > $WORK_DIR/align/${file}/${file}_aboveMapq30_nodup2_prim.bedgraph
python ~/scratch/202312241147Polina/loxP_project/plot_depth_hist.py $WORK_DIR/align/${file}/${file}_aboveMapq30_nodup2_prim.bedgraph $WORK_DIR/align/${file}/${file}_aboveMapq30_nodup2_prim.modified.bedgraph

done
