bcftools mpileup --threads 10 -f ../../bwa/chr6.fna ../samtools/marked_SRR10720416_chr6_f.bam
2> bcftools_mpileup.log | bcftools call --threads 10 -mv -o chr6.vcf 2> bcftools_call.log
Версия 1.9, htslib 1.9 (и далее та же версия для всех bcftools), log-файл mpileup, log-файл call.
Использование параметра --threads должно было ускорить выполнение пограммы, но, как выянилось, это работает только если выходной файл представляет собой архив, поэтому тут их использование было бессмысленным.
Параметр -f задаёт файл с референсом, с которым сравниваем, -v - оставить в выдаче только сайты с вариантами, -m - выбор другой модели для описания редких вариантов и вариантов в нескольких аллелях
В начале файла много общей информации о запущенной команде и отчасти о процессе выполнения (STDERR отделила, но эти на STDOUT). Дальше тело файла: таблица с информацией об отличиях: название референса; позиция варианта; ID варианта (у меня везде, где посмотрела они неизвестны и стоят просто точки); что на этом месте в референсе; что вместо этого в образце; качество; флаг, описывающий фильтры, которые мы наложили на варианты (если они были); различная информация о варианте в формате ключ=значение через точку с запятой; формат дополнительной информации (у меня GT:PL, т.е. генотип:оценка сходства); сама эта дополнительная информация для каждого образца).
bcftools stats chr6.vcf > chr6.stats 2> bcftools_stats.log
log-файл (пуст, поэтому в дальнейшем использовалось общее перенаправление), результат.
| Вариантов | 445 947 |
| Однонуклеотидных | 436 927 (98%) |
| Инделей | 9 020 (2%) |
bcftools filter -i'%QUAL>30 && DP>50' chr6.vcf -o chr6_filtered.vcf 2> bcftools_filter.log
log-файл(пуст)
bcftools stats chr6_filtered.vcf &> chr6_filtered.stats
| Вариантов | 23 797 |
| Однонуклеотидных | 23 493 (99%) |
| Инделей | 304 (1%) |
Количества уменьшились примерно в 20 раз, но соотношения однонуклеотидных и инделей изменилось незначительно (то есть варианты плохого качества были распределены между ними приблизительно равномерно).
3.1 Поиск хорошо покрытых участков
bedtools genomecov -bg -ibam ../samtools/marked_SRR10720416_chr6_f.bam > chr6.bed
2> bedtools_genomecov.log
log-файл(пуст)
Параметры: -bg - задаёт формат выходного файла как BedGraph, -ibam - даёт понять, что входным файлом будет bam
Файл bed содержит к каждой строке название хромосомы, позицию начала блока с одинаковой покрытостью, позицию конца, покрытость этого блока. Поскольку у нас не RNA-seq, проблем сплайсирования нет и -split использовать не нужно.
mawk '$4>50' chr6.bed > chr6_covered50.bed 2> mawk.log wc -l chr6_covered50.bed
Версия mawk 1.3.3, log-файл(пуст)
Выделили только строки с покрытием больше 50. В файле 1 352 481 строка, т.е. 1 352 481 участков, покрытых больше 50 раз.
3.2 Варианты в хорошо покрытых участках
bedtools intersect -a chr6_filtered.vcf -b chr6_covered50.bed -wa -header
> chr6_covered50.vcf 2> bedtools_intersect.log
log-файл(пуст)
bcftools stats chr6_covered50.vcf > chr6_covered50.stats
| Вариантов | 24 609 |
| Однонуклеотидных | 23 493 |
| Инделей | 1 116 |
То есть все однонуклеотидные варианты оказались именно на хорошо покрытых участках, число инделей даже увеличилось.
Чтобы понять, причину этого, я построила это всё для одних только инделей:
bcftools filter -i'TYPE="indel"' chr6_filtered.vcf -o chr6_indels.vcf
bedtools intersect -a chr6_indels.vcf -b chr6_covered50.bed -wo -header > chr6_indels_covered50.wo
В получившимся файле, содержащем для каждого пересечения записи как из vcf файла, так и из bed, можно легко видеть причину. Один индель при пересечении с bed файлом разрезается на несколько, и не только потому, что внутри него находится участок с покрытием меньше 50, но и потому, что последовательные участки с покрытием, скажем, 59 и 60 считаются за разные участки. Пример для одного исходного инделя (для удобства восприятия оставила только часть столбцов):
CHROM POS REF ALT bed: chrom start end cov overlap
NC_000006.12 4031699 CAGAAGAAGA CAGAAGA NC_000006.12 4031694 4031699 63 1
NC_000006.12 4031699 CAGAAGAAGA CAGAAGA NC_000006.12 4031702 4031704 60 2
NC_000006.12 4031699 CAGAAGAAGA CAGAAGA NC_000006.12 4031704 4031705 59 1
NC_000006.12 4031699 CAGAAGAAGA CAGAAGA NC_000006.12 4031705 4031706 58 1
NC_000006.12 4031699 CAGAAGAAGA CAGAAGA NC_000006.12 4031706 4031708 61 2
Для исправления этого я применила команду
bedtools merge -i chr6_covered50.bed > chr6_cov50_merged.bed
Теперь граничащие концами фрагменты объединились в один (потерялась информация о покрытии, но, поскольку эта операция производилась с файлом, в котором только участки с покрытием больше 50, это нестрашно). Видно, что всё в порядке и индель режентся только если в середине кусок меньше 50:
CHROM POS REF ALT bed: chrom start end overlap
NC_000006.12 4031699 CAGAAGAAGA CAGAAGA NC_000006.12 4031582 4031699 1
NC_000006.12 4031699 CAGAAGAAGA CAGAAGA NC_000006.12 4031702 4032768 6
bcftools stats chr6_indels_cov50_merged.wo > chr6_indels_cov50_merged.stats
Теперь инделей (точнее, их фрагментов) всего 359, что уже больше похоже на правду.
И можно произвести более верную оценку числа фрагментов, покрытых больше 50 раз (с учётом слияния соседей):
wc -l chr6_cov50_merged.bed
Результат 21 793. Если сравнить с изначальным (до слияния) значением 1 352 481, можно видеть, насколько сильно были "расщеплены" эти фрагменты.