Ресеквенирование. Поиск полиморфизмов у человека, продолжение

Практикум 13


1. Получение вариантов

a. BCFtools — набор утилит, позволяющий работать с файлами формата VCF (Variant Call Format) и его бинарным аналогом BCF.

Получим из dubble_mark.bam (см. 12 практикум) файл с отличиями от референсного генома в формате 13.vcf: bcftools mpileup --threads 10 -f ../chr8.fna dubble_mark.bam 2> 13_1.log | bcftools call -mv -o 13.vcf 2> 13_2.log.

  • bcftools mpileup генерирует VCF-файл;
  • -f ../chr8.fna задаёт имя faidx-проиндексированного референса в формате fasta;
  • bcftools call указывает/называет вариации (SNP/индели);
  • -mv: -m запускает альтернативную модель консенсусного указания (большие проблемы с переводом слова calling, не уверен, что смысл удаётся сохранить) для мультиаллельных и редких указаний. Это необходимо во избежание ограничений стандартного метода. -v позволяет вполучить на выходе только сайты с вариациями.
  • -o 13.vcf задаёт имя выходного файла;
  • --threads 10 задаёт количество задейтсвованных в процессе ядер;
  • 2> 13_?.log перенаправляет сообщения о работе программы в log-файл.

b. Файл 13.vcf состоит из шапки с общим описанием тела файла. Само тело файла из 8 обязательных столбцов и неограниченного числа дополнительных. Описание столбцов приведено в Таблице 1.


Таблица 1. Столбцы файла 13.vcf
№ Название столбца Описание
1 CHROM Имя референсной последовательности
2 POS Номер позиции вариации на референсной последовательности
3 ID Идентификатор вариации
4 REF Номер позиции первого закартированного на референс нуклеотида. Если картировать чтение не получилось, значение параметра равно нулю.
5 ALT Список альтернативных аллелей для этой позиции
6 QUAL Качество, соответствующее вариантам для этих аллелей
7 FILTER Показывает, через какой набор фильтров прошла вариация
8 INFO Список с элементами ключ=значение, описывающий вариацию
9 доп. FORMAT Список кодов, описывающих образцы (расшифровка в шапке файла)
10 доп. SAMPLEs Список значений для полей из FORMAT

c. Командой bcftools stats 13.vcf > 13_stats.txt получили файл со статистикой (выходной файл тоже формата .vcf, но я сглупил и написал .txt).

d. Всего получилось 419378 вариантов, 412349 однонуклеотидных замен, 7029 инделей.

2. Фильтрация вариантов

a. Командой bcftools filter -i'%QUAL>30 && DP>50' 13.vcf -o 13_filt.vcf отфильтруем варианты с качеством больше 30 и покрытием больше 50.

b. Воспользуемся командой bcftools stats 13_filt.vcf > 13_stats_filt.vcf. Из полученного файла видим, что теперь всего 20478 вариантов, 20272 SNP, 206 инделей, то есть все значения упали примерно в 20 раз по сравнению с имеющимся до фильтрации (п. 1d).

3. Изучение покрытия

a. Командой bedtools genomecov -bg -ibam dubble_mark.bam > gen_reg_1.bed 2> gen_reg_1.log получили участки генома, покрытые нашими чтениями хотя бы 1 раз.

b. Файл gen_reg_1.bed состоит из четырёх столбцов:

  1. имя референса;
  2. начало вариации;
  3. конец вариации;
  4. число покрытий.

c. Участки генома с покрытием больше 50 были отобраны командой awk '$4 > 50 {print}' gen_reg_1.bed > final.bed.

d. Команда wc -l final.bed помогла узнать, что всего таких участков 1008269.


Ссылки для себя на будущее