#! /usr/bin/bash

source confile.txt

#фильтрация чтений

java -jar /usr/share/java/trimmomatic.jar PE -threads 10 -phred33 ${SRR}_1.fastq.gz ${SRR}_2.fastq.gz ${SRR}_1_paired.fq.gz ${SRR}_1_unpaired.fq.gz ${SRR}_2_paired.fq.gz ${SRR}_2_unpaired.fq.gz TRAILING:20 MINLEN:50 2> trimming.log
mkdir paired
mkdir unpaired


#проверка качества в fastqc

fastqc ${SRR}_1_*.fq.gz  2>> log_fastqc.txt
fastqc ${SRR}_2_*.fq.gz  2>> log_fastqc.txt

mkdir fastqc12
mv *fastqc* fastqc12

mv ${SRR}_?_paired.fq.gz paired
mv ${SRR}_?_unpaired.fq.gz unpaired


#картирование чтений

hisat2 -p 1 --no-spliced-alignment -x reference/${chr}  -1 paired/${SRR}_1_paired.fq.gz -2 paired/${SRR}_2_paired.fq.gz -S ${SRR}.sam 2> log_hisat2.txt

#конвертация sam в bam

mkdir mapping
ls -lh ${SRR}.sam > log_samtools.txt
samtools sort -o ${SRR}.bam ${SRR}.sam 2>> log_samtools.txt
rm ${SRR}.sam
ls -lh ${SRR}.bam >> log_samtools.txt
samtools index ${SRR}.bam

#анализ bam-файла

samtools flagstat ${SRR}.bam > ${SRR}_bam_all_info.txt
mv ${SRR}_bam_all_info.txt mapping

#получение чтений, картированных на хромосому

samtools view -h -bS ${SRR}.bam ${chr} > ${SRR}_sec.bam
mv ${SRR}.bam mapping

#получение правильно картированных чтений

samtools view -f 0x2 -bS ${SRR}_sec.bam > ${SRR}_thi.bam
samtools flagstat ${SRR}_thi.bam > ${SRR}_bam_info.txt
samtools index ${SRR}_thi.bam
mv ${SRR}_bam_info.txt mapping
mv ${SRR}.bam.bai mapping

#получение вариантов

bcftools mpileup -f reference/${chr}.fa mapping/${SRR}.bam |  bcftools call -mv -o ${SRR}.vcf  2> log_bcftools.txt
bcftools stats ${SRR}.vcf > ${SRR}_vcf_info.txt

#фильтрация вариантов

bcftools filter -i'%QUAL>30 && DP>50' ${SRR}.vcf > ${SRR}_filtered.vcf
bcftools stats ${SRR}_filtered.vcf > ${SRR}_vcf_filtered_stats.txt

mkdir vcf
mv *vcf* ./vcf

#переместим все сообщения об ошибках в одну удобную папку

mkdir logs
mv *log* ./logs