Markov Ivan, Введение в анализ данных NGS (часть 3)

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

Осуществлено с помощью команды "bcftools mpileup" по шаблону, за исключением подключения 12 ядер (--threads 12), опция "-f" используется для обозначения референсного fasta-файла (chr7.fna), команда была применена к файлу маркированных дублированных чтений (marked.bam). Результат работы этой команды был перенаправлен черех пайп команде "bcftools call", у нее были подключены три опции: -m (альтернативная модель для работы с мультиаллельными и редкоими вариантами), -v (в выводе только сайты вариантов) и -o (название выходного файла). Итоговая команда: 

bcftools mpileup --threads 12 -f ../chr7.fna marked.bam |  bcftools call -mv -o 1.vcf

Верхние 35 строк, отделенные командой "head -n 35 1.vcf > vcfhead.txt", можно просмотреть по ссылке.

VCF - фйал тестового формата. Он содержит мета-информационные строки, строку заголовка и после строки данных, каждая из которых содержит информацию о позиции в геноме. Далее представлены поля строк данных и их короткое описание.

  1. CHROM: хромосома; идентификатор из референсного генома или ID в угловых скобках строки, указывающей на контиг в файле собрания.
  2. POS: позиция; референсная позиция, в котором первое основания имеет позицию 1.
  3. ID: идентификатор; отделенный полуколонкой список уникальных идентификаторов.
  4. REF: основание референса (или N, если неизвестно).
  5. ALT: альтернативное основание; разделенный запятыми список нереференсных аллелей.
  6. QUAL: качество; очки качества надежности по шкале Phred.
  7. FILTER: статус фильтра.
  8. INFO: дополнительные поля, построеные по конструкции =[,data].
  9. FORMAT: дополнительное поле, отражающее ,например, отсутствие GQ.

Данные взяты из мануала.

Просмотр статистики выполнен командой "bcftools stats":

bcftools stats 1.vcf > vcfstats.txt

Результат можно просмотреть по ссылке. В этом файле указано общее число вариантов - 402655, из них однонуклеотидных замен - 394681 и инделей - 7974.

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

Два фильтра, используемые в команде ниже - качество более 30 и покрытие более 50:

bcftools filter --threads 12 -i'%QUAL>30 && DP>50' 1.vcf -o 2.vcf

Для построения статистики вновь использавана "bcftools stats": 

bcftools filter --threads 12 -i'%QUAL>30 && DP>50' 1.vcf -o 2.vcf

Выходной файл находится по ссылке. Количество варинтов уменьшилось с 402655 до 23432.

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

Для того чтобы получить участки генома, покрытые хотя бы 1 раз была выполнена команда команда "bedtools genomecov": 

bedtools genomecov -ibam marked.bam -bg > genomecov.bed

По-умолчанию .bed файл содержит 54 полей: хромосому (или весь геном), начало варианта, конец варианта, и длину покрытия.

bg .bed формат
Рисунок 1. Поля bg .bed-файла

С помощью программы awk были выбраны участки геном, покрытые более 50 раз:

awk '$4 > 50' genomecov.bed > mt50.bed

В итоге образовался файл.

Для подсчета строк в файле была использована стандартная команда "wc":

wc -l mt50.bed
1492263 mt50.bed

Итого в нем 1492263 строки.