1. Получение вариантов.
Получим варианты командой:
bcftools mpileup -f ./bwa/chr7.fna marked_chr7.bam | bcftools call -mv -o chr7.vcf
Где:
bcftools mpileup - алгоритм, создающий bcf/vcf файлы.
-f - референсный файл уже проиндексирован samtools faidx
marked_chr7.bam - дублированные маркированные чтения
-mv - опция вызова параметров для многоаллельного или редковариантного коллинга и для записи в файл только вариантов
-o - задание имени выходного файла
При открытии файла chr7.vcf сразу заметна шапка. Она отличается от тела файла тем, что строки начинаются с "##".
Ниже шапки(собственно в теле) расположена таблица. Некоторые из них:
CHROM - название полиморфизма с полиморшизмов или инделью
POS - позиция этого полиморфизма или инделя
ID - идентификатор полиморфизма или инделя
REF - нуклеотиды, находящиеся на этой позиции в референсе
ALT - найденные варианты
QUAL - качество вариантов
FILTER - какие фильтры смог пройти этот вариант
INFO - поле с дополнительной информацией
Проанализируем vcf файл с помощью команды:
bcftools stats chr7.vcf > stats.txt
Найдём общее количество вариантов - 388803, из них индели - 8183, а SNP - 380620.
2. Фильтрация вариантов.
Отфильтруем варианты. По глубине покрытия больше чем 50 и качеству большему 30. Команда:
bcftools filter -i'%QUAL>30 && DP>50' chr7.vcf -o chr7_filter.vcf
Применим уже знакому нам команду для того, что проанализировать файл:
bcftools stats chr7_filter.vcf > stats_filter.txt
В полученном файлк содержалось 22990 вариантов, причём инделей - 340, а SNP - 22650.
Очевидно, что в процесс фильтрации количество вариантов значительно сократилось (примерно в 17-18 раз).
3. Изучение покрытия.
Получим участки генома, покрытие нашими чтениями хотя бы 1 раз с помощью команды:
bedtools genomecov -bg -ibam marked_chr7.bam > chr7_cover.bed
В полученном файле было 4 колонки. Они представляют собой хромосомы (имя референса референса), начало и коней полиморфизма, покрытие.
Для фильтрации полиморфизмов с покрытием больше чем 50 была использована команда:
awk '{if($4 > 50) {print $1"\t"$2"\t"$3"\t"$4}}' chr7_cover.bed > chr7_filter.bed
Посчитаем количество строк командой wc -l chr7_filter.bed. Их оказалось 1597099.