#export PATH=$PATH:~/programms/bwa-0.7.17;
WORK_DIR=~/4c_plasmid/2021-10-20/
cd $WORK_DIR; #mkdir align; 
cd align;
# Declare a string array with type
declare -a files=("103_S103" "104_S104" "105_S105")

for file in "${files[@]}"
do
echo ${file}
mkdir ${file}
cat $WORK_DIR/cutadapt/${file}_1.fq | grep pr_plasmid1 -A 4 >> $WORK_DIR/cutadapt/${file}_pr_plasmid1.fq
cat $WORK_DIR/cutadapt/${file}_2.fq | grep pr_plasmid1 -A 4 >> $WORK_DIR/cutadapt/${file}_pr_plasmid1.fq
cat $WORK_DIR/cutadapt/${file}_1.fq | grep pr_plasmid2 -A 4 >> $WORK_DIR/cutadapt/${file}_pr_plasmid2.fq
cat $WORK_DIR/cutadapt/${file}_2.fq | grep pr_plasmid2 -A 4 >> $WORK_DIR/cutadapt/${file}_pr_plasmid2.fq
## align pr_plasmid1
bwa mem -t 5 -M -R '@RG\tID:${file}_pr_plasmid1' /home/polina/4c_plasmid/genome/hg38_plasmid $WORK_DIR/cutadapt/${file}_pr_plasmid1.fq | samtools view -bS -@ 5 | samtools sort > $WORK_DIR/align/${file}/${file}_pr_plasmid1.sorted.bam;
samtools index $WORK_DIR/align/${file}/${file}_pr_plasmid1.sorted.bam > $WORK_DIR/align/${file}/${file}_pr_plasmid1.sorted.bam.bai;
## align pr_plasmid2
bwa mem -t 5 -M -R '@RG\tID:${file}_pr_plasmid2' /home/polina/4c_plasmid/genome/hg38_plasmid $WORK_DIR/cutadapt/${file}_pr_plasmid2.fq | samtools view -bS -@ 5 | samtools sort > $WORK_DIR/align/${file}/${file}_pr_plasmid2.sorted.bam;
samtools index $WORK_DIR/align/${file}/${file}_pr_plasmid2.sorted.bam > $WORK_DIR/align/${file}/${file}_pr_plasmid2.sorted.bam.bai;
## merge files
samtools merge $WORK_DIR/align/${file}/${file}_pr_plasmid1_pr_plasmid2.bam $WORK_DIR/align/${file}/${file}_pr_plasmid1.sorted.bam $WORK_DIR/align/${file}/${file}_pr_plasmid2.sorted.bam;
samtools index $WORK_DIR/align/${file}/${file}_pr_plasmid1_pr_plasmid2.bam > $WORK_DIR/align/${file}/${file}_pr_plasmid1_pr_plasmid2.bam.bai;
## align paired reads
bwa mem -t 5 -M /home/polina/4c_plasmid/genome/hg38_plasmid $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
bedtools genomecov -ibam $WORK_DIR/align/${file}/${file}.sorted.bam -bg > $WORK_DIR/align/${file}/${file}.begraph
## create sam
samtools view $WORK_DIR/align/${file}/${file}.sorted.bam > $WORK_DIR/align/${file}/${file}.sorted.sam
done
