config_pr11-13.txt

#!/bin/sh

# Референсная последовательность хромосомы
inp_genome='../DATA/hg38/Homo_sapiens.GRCh38.dna.chromosome.6.fa'

# Прямые чтения
inp_reads_for='../DATA/dna_reads/SRR10720404_1.fastq.gz'

# Обратные чтения
inp_reads_rev='../DATA/dna_reads/SRR10720404_2.fastq.gz'

# Префикс FASTQ файла
srr='SRR10720404'

# Номер хромосомы
chr='6'

# Количество потоков
threads='8'





script_pr11-13.sh

#!/bin/sh

#########################################
#                SOFT                   #
#########################################
# trimmomatic 0.39                      #
# hisat2 2.2.1                          #
# samtools 1.17                         #
# bcftools 1.11                         #
#########################################

config_file=$1

. ${config_file}

# Создание необходимых директорий
mkdir genome
mkdir trimmed
mkdir mapped
mkdir variants

# Индексирование хромосомы для hisat и samtools
hisat2-build ${inp_genome} ./genome/index
samtools faidx ${inp_genome}

# Триммирование исходных чтений
TrimmomaticPE -threads ${threads} -phred33 ${inp_reads_for} ${inp_reads_rev} trimmed/${srr}_trimmomatic_forward_paired.fq.gz trimmed/${srr}_trimmomatic_forward_unpaired.fq.gz trimmed/${srr}_trimmomatic_reverse_paired.fq.gz trimmed/${srr}_trimmomatic_reverse_unpaired.fq.gz TRAILING:20 MINLEN:50

# Картирование триммированных чтений на референснвй геном
hisat2 -p ${threads} -x genome/index -1 trimmed/${srr}_trimmonatic_forward_paired.fq.gz -2 trimmed/${srr}_trimmonatic_reverse_paired.fq.gz --no-spliced-alignment -S mapped/${srr}.sam

# Конвертация sam в bam и сортировка картированных чтений
samtools sort -@ ${threads} -o mapped/${srr}.bam mapped/${srr}.sam
samtools index -@ ${threads} ${srr}.bam
samtools flagstat -@ ${threads} mapped/${srr}.bam > mapped/${srr}_flagstat.txt

# Получение чтений, картированных на нашу хромосому
samtools view -@ ${threads} -h -bS mapped/${srr}.bam ${chr} > mapped/${srr}_${chr}_Chr.bam
samtools view -@ ${threads} -f 0x2 -bS mapped/${srr}_${chr}_Chr.bam > mapped/${srr}_${chr}_Chr_sorted.bam
samtools flagstat -@ ${threads} mapped/${srr}_${chr}_Chr_sorted.bam > mapped/${srr}_${chr}_Chr_sorted_flagstat.txt
samtools index -@ ${threads} mapped/${srr}_${chr}_Chr_sorted.bam

# Получение и фильтрация вариантов
bcftools mpileup -f ${inp_genome} mapped/${srr}_${chr}_Chr_sorted.bam |  bcftools call -mv variants/${srr}_${chr}_Chr.vcf
bcftools stats variants/${srr}_${chr}_Chr.vcf > variants/${srr}_${chr}_Chr_stats.txt
bcftools filter -i '%QUAL>30 && DP>50' variants/${srr}_${chr}_Chr.vcf > variants/${srr}_${chr}_Chr_filtered.vcf
bcftools stats variants/${srr}_${chr}_Chr_filtered.vcf > variants/${srr}_${chr}_Chr_filtered_stats.txt