Получение вариантов
Для получение вариантов воспользуемся командой:
$ bcftools mpileup -f ../bwa/chr10.fna mark_chr10.bam | bcftools call -mv -o variants.vcf
bcftools mpileup – программа, которая генерирует vcf файл с вероятностями разных вариантов на основании выравнивания.
-f ../bwa/chr10.fna – fasta-файл с референсным геномом.
mark_chr10.bam – файл с картированными ридами.
| – пайп)
bcftools call – команда, которая достаёт из вывода предыдущей программы нужные строки.
-mv – модель, которая ищет только многовариантные замены.
-o variants.vcf – выходной файл.
Полученный vcf-файл содержит заголовок, содержащий информацию о том, как файл был получен. После заголовка идут строки, которые описывают различные варианты. Строки вариантов содержат следующие поля: CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO.
Чтобы узнать, сколько всего каких замен, применим следующую команду:
$ bcftools stats variants.vcf > stats.txt
Всего было найдено 367042, из них 360206 однонуклеотидных полиморфизмов и 6836 инделей.
Фильтрация вариантов
Чтобы отфильтровать варианты, имеющие качество больше 30 и покрытие больше 50 применим следующую команду:
$ bcftools filter -i'%QUAL>30 && DP>50' variants.vcf -o filter_variants.vcfbcftools filter -i'%QUAL>30 && DP>50' variants.vcf -o filter_variants.vcf
Посмотрим на выдачу stats для этого файла:
$ bcftools stats filter_variants.vcf > stats2.txt
В новом файле 22854 варианта, из них 22568 однонуклеотидных полиморфизмов и 286 инделей. Числа сильно уменьшились по сравнению с предыдущим файлом.
Изучение покрытия
Получим участки генома, покрытые чтениями хотя бы один раз. Для этого применим следующую команду:
$ bedtools genomecov -bg -ibam mark_chr10.bam > coverage.bed
В полученном файле каждая строка содержит следующую информацию: имя последовательности; координаты в этой последовательности; сколько раз участок по этим координатам был покрыт. Выглядит это примерно вот так:
$ head coverage.bed
NC_000010.11 8317 8336 2
NC_000010.11 8787 8806 2
NC_000010.11 9999 10000 1
NC_000010.11 10000 10001 2
NC_000010.11 10001 10002 3
NC_000010.11 10002 10003 4
NC_000010.11 10003 10004 5
NC_000010.11 10004 10008 6
NC_000010.11 10008 10009 7
NC_000010.11 10009 10010 9
Теперь из полученного файла, отберём только такие участки, которые покрыты более 50 раз, и посчитаем сколько их:
$ awk '$4 > 50' coverage.bed > 50coverage.bed
$ wc -l 50coverage.bed
1313052 50coverage.bed
Всего получилось 1313052 участка, покрытых более 50 раз.