#PBS -q dl560g10q@vm-pbs2
#PBS -l select=1:ncpus=5:mem=10g
#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=~/loxP_project/2022-06-14_rings/
cd $WORK_DIR; mkdir align; 
cd align;
# Declare a string array with type

declare -a files=("HAP8VIR_sub11" "HAP8VIR_sub17" "HAP8VIR_sub26" "HAP8VIR_sub30" "HAP8VIR_sub9" "HAP8VIR_sub14" "HAP8VIR_sub19" "HAP8VIR_sub27" "HAP8VIR_sub4" "HAP8VIR_sub15" "HAP8VIR_sub20" "HAP8VIR_sub28" "HAP8VIR_sub5" "HAP8VIR_sub16" "HAP8VIR_sub21" "HAP8VIR_sub29" "HAP8VIR_sub7")

for file in "${files[@]}"
do
echo ${file}
mkdir ${file}
cat $WORK_DIR/cutadapt/${file}_1.fq | grep ada1 -A 4 >> $WORK_DIR/cutadapt/${file}_ada1.fq
cat $WORK_DIR/cutadapt/${file}_2.fq | grep ada1 -A 4 >> $WORK_DIR/cutadapt/${file}_ada1.fq
cat $WORK_DIR/cutadapt/${file}_1.fq | grep ada2 -A 4 >> $WORK_DIR/cutadapt/${file}_ada2.fq
cat $WORK_DIR/cutadapt/${file}_2.fq | grep ada2 -A 4 >> $WORK_DIR/cutadapt/${file}_ada2.fq
## align ada1
bwa mem -t 5 -M -R '@RG\tID:${file}_ada1' /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}_ada1.fq | samtools view -bS -@ 5 | samtools sort > $WORK_DIR/align/${file}/${file}_ada1.sorted.bam;
samtools index $WORK_DIR/align/${file}/${file}_ada1.sorted.bam > $WORK_DIR/align/${file}/${file}_ada1.sorted.bam.bai;
## align ada2
bwa mem -t 5 -M -R '@RG\tID:${file}_ada2' /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}_ada2.fq | samtools view -bS -@ 5 | samtools sort > $WORK_DIR/align/${file}/${file}_ada2.sorted.bam;
samtools index $WORK_DIR/align/${file}/${file}_ada2.sorted.bam > $WORK_DIR/align/${file}/${file}_ada2.sorted.bam.bai;
## merge files
samtools merge $WORK_DIR/align/${file}/${file}_ada1_ada2.bam $WORK_DIR/align/${file}/${file}_ada1.sorted.bam $WORK_DIR/align/${file}/${file}_ada2.sorted.bam;
samtools index $WORK_DIR/align/${file}/${file}_ada1_ada2.bam > $WORK_DIR/align/${file}/${file}_ada1_ada2.bam.bai;
## 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
samtools index $WORK_DIR/align/${file}/${file}.sorted.bam > $WORK_DIR/align/${file}/${file}.sorted.bam.bai
## create coverage file
~/programms/bedtools2/bin/bedtools genomecov -ibam $WORK_DIR/align/${file}/${file}.sorted.bam -bg > $WORK_DIR/align/${file}/${file}.bedgraph
## create sam
samtools view $WORK_DIR/align/${file}/${file}.sorted.bam > $WORK_DIR/align/${file}/${file}.sorted.sam
python ~/loxP_project/plot_depth_hist.py $WORK_DIR/align/${file}/${file}.bedgraph $WORK_DIR/align/${file}/${file}.modifyed.bedgraph
done
