Поиск вариантов

Задание 1: Получение вариантов

Получим, собственно варианты, воспользовавшись средствами bcftools:

$ bcftools mpileup -f ../bwa/chr3.fna ../cartir/aln-se.mark.dupl.bam | bcftools call -mv -o vars.vcf

• bcftools mpileup - обращение к программе mpileup из набора инструментов bcftools;
• -f ../bwa/chr3.fna - такой опцией мы поясняем программе, что референс у нас будет в fasta-формате + пишем путь к референсу;
• ../cartir/aln-se.mark.dupl.bam - наш файл с "корректно" картированными чтениями, в котором ещё и дубликаты промаркированны;
• Далее через пайп мы выход bcftools mpileup подаём программе bcftools call;
• bcftools call - обращение к программе call из набора инструментов bcftools;
• -mv - опция, которая говорит программе выдавать только отличия нашей сборки от референса, причём в некой альтернативной модели вызова мультиаллелей и редких вариантов замен (тут не совсем понятно, что имеется ввиду);
• -o vars.vcf - задаём выходной файл.

Внутри полученный vcf-файл выглядит таким образом:

Первые несколько строк начинаются с символа "#" и несут служебную информацию, объясняя версию программы, показывая использованную команду, приводя расшифорвки для используемых в теле файла сокращений. Само тело файла представляет из себя таблицу, в которой приведена различная информация о разного рода заменах в нашей сборке по сравнению с референсом.

Далее проанализируем полученный нами vcf-файл:

$ bcftools stats vars.vcf > vars.vcf.stats

• Получилось вариантов: 425305;
• Количество SNP: 416801 (98%);
• Количество инделей: 8504 (2%).

Задание 2: Фильтрация вариантов

Собственно, сейчас нам нужно отфильтровать наше огромное количество вариантов по необходимым нам критериям. Сделаем мы это также с помощью bcftools:

$ bcftools filter -i'%QUAL>30 && DP>50' vars.vcf > vars_filtered.vcf

В нашей ситуации мы отфильтровали наши варианты по качеству (QUAL) и по количеству чтений, которые наш вариант покрывают (DP).

Проанализируем получившийся файл с отфильтрованными вариантами:

$ bcftools stats vars_filtered.vcf > vars_filtered.vcf.stats

• Получилось вариантов: 25814;
• Количество SNP: 25460 (98.6%);
• Количество инделей: 354 (1.37%).

Как мы можем заметить, количество вариантов что в общем, что в частных случаях SNP и инделей, уменьшилось, причём почти в 20 раз, и при этом доли SNP и инделей в общем количестве вариантов также изменились: доля SNP стала больше, а доля инделей - меньше. Уменьшение количества вариантов довольно чётко говорит нам о том, что фильтрация прошла успешно, так как большая часть вариантов (93.93%) отсеялась по нашим довольно жёстким критериям.

Задание 3: Изучение покрытия

У нас есть небольшая проблема - мы не знаем нашего таргета. Исходя из этого, поступим следующим образом: найдем покрытые участки, отберем участки с хорошим покрытием и такие участки [с хорошим покрытием] и будем считать таргетом, а после проанализируем их.

Сначала найдём в принципе покрытые участки. Воспользуемся для этого инструментами bedtools:

$ bedtools genomecov -ibam ../cartir/aln-se.mark.dupl.bam -bg > aln-se_cov.bed

• bedtools genomecov - обращение к программе genomecov из набора bedtools;
• -ibam ../cartir/aln-se.mark.dupl.bam - подача входного bam-файла (опция поясняет программе, что входной файл у меня именно .bam);
• -bg - выход программы представить в удобном формате BedGraph;
• > aln-se_cov.bed - перенаправление выхода в файл.

Внутри полученный файл aln-se_cov.bed устроен таким образом:

То есть, у нас каждая строчка представлена так: имя хромосомы, координаты полученного после разбиения участка (через табулятор, 2 штуки: первая и последняя) и покрытие (кол-во ридов, которые покрывают этот участок) этого участка (как я понял, каждого нуклеотида в участке).

Далее я отберу такие участки генома, которые покрыты более 50 раз:

$ awk '{if($4 > 50) {print $1"\t"$2"\t"$3"\t"$4}}' aln-se_cov.bed > aln-se_cov50.bed

Посчитаю количество строк в получившемся файле:

$ wc -l aln-se_cov50.bed
1624352 aln-se_cov50.bed

И попробую найти варианты, которые попали в хорошо покрытые участки:

$ bedtools intersect -a vars_filtered.vcf -b aln-se_cov50.bed > vars_cov50.vcf

• bedtools intersect - обращение к программе intersect из набора bedtools;
• -a vars_filtered.vcf - подача файла с отфильтрованными вариантами;
• -b aln-se_cov50.bed - подача файла с участками с хорошим покрытием из нашей сборки;
• > vars_cov50.vcf - перенаправление выхода в файл.

Отлично, получилось! Варианты, находящиеся на участках с прекрасным покрытием найдены! =)