Для получения вариантов, воспользуемся следующей командой: bcftools mpileup -f ../bwa/chr8.fna ../bwa_mem/marked.bam | bcftools call -mv -o file.vcf.
Посмотрим как устроен полученный vcf-файл. Первые строчки начинаются с символов "##". В них содержится общая информация о теле файла, которое в свою очередь состоит из 9 обязательных столбцов. Описания колонок находится в таблице 1.
| Название столбца | Описание |
|---|---|
| CHROM | Имя последовательности, в которой есть вариант |
| POS | Позиция полиморфизма или инделя |
| ID | Индивидуальный номер варианта |
| REF | Референсный нуклеотид |
| ALT | Список аллельных вариантов |
| QUAL | Качество данного варианта |
| FILTER | Информация о наборе фильтров, который был применен к данному варианту |
| INFO | Основная информация о данном варианте |
| FORMAT | Описание образца |
Далее проанализируем наш файл с помощью команды bcftools stats file.vcf > file_stats.txt. Полученный файл находится тут. Из анализирующего файла узнаем: всего вариантов получилось 350182, из них 344080 SNPs и 6102 инделей.
Критериев фильтрации существует достаточно много. Применим к полученному vsf-файлу следующие команды:
bcftools filter -i'%QUAL>30 && DP>50' file.vcf -o filtered_file.vcf
bcftools stats filtered_file.vcf > filtered_stats.txt - сразу посмотрим на статистические значения получившегося файла.
В отфильтрованном файле оказалось 20285 записей. Из них 20030 о SNPs и 255 о инделях. Статистический файл можно найти тут.
Так как мы не знаем таргета, поступим следующим образом: найдем покрытые участки, отберем участки с хорошим покрытием, такие участки будем считать таргетом, проанализируем их.
С помощью команды bedtools genomecov -bg -ibam ../bwa_mem/marked.bam > genomecov_chr8.bed получим участки хромосомы, покрытые нашими чтениями хотя бы раз. Получившийся файл genomecov_chr8.bed устроен следующим образом: в первой колонке название референса (восьмой хромосомы в моем случае), во второй колонке позиция начала варианта, в третьей колонке позиция конца и в последней количество покрытий.
Теперь отберем те участки генома, которые покрыты более 50 раз: awk '{if($4 > 50) {print $1"\t"$2"\t"$3"\t"$4}}' genomecov_chr8.bed > 50_chr8.bed
В итоговом файле посчитаем количество оставшихся строк (то есть это количество вариантов, которые хорошо покрылись): wc -l 50_chr8.txt. Получилось 1116548 строк.
История команд лежит тут.