Картирование парных триммированных чтений на геном человека было выполнено с помощью команды bwa mem
bwa mem -t 10 ../../bwa/chr10.fna ../trim_pair/paired_out_1.fq.gz ../trim_pair/paired_out_2.fq.gz > map_16_ch10.sam 2> mem_log
a. Файлы формата .sam изначально были придуманы для хранения биологических последовательностей, выровненных с каким-либо референсом. Устройство формата довольно понятное: есть заголовки, начинающиеся с "@", и строки, содержащие информацию о ридах.
b. Информация о ридах, в свою очередь, разбивается по логическим колонкам: 11 обязательных, но иногда добавляется 12-я, содержащая дополнительную информацию. Подробнее о том, в каком столбце что содержится, написано в Таблице 1.
| № колонки | Поле | Тип данных | Краткое описание |
|---|---|---|---|
| 1 | QNAME | String | ID чтения |
| 2 | FLAG | Int | Комбинация значиений особых FLAGs (флагов), которые в битах дают характеристику каждому отдельному риду |
| 3 | RNAME | String | Название референса, на который картируется конкретный рид |
| 4 | POS | Int | Координата первого закартированного на референс нуклеотида. Равняется 0, если рид не отображен и не имеет координат |
| 5 | MAPQ | Int | Качество картирования, вычисляемое по формуле Q = -10lgP, где P - вероятность ошибки. Значение 255 свидетельствует о недоступности рида |
| 6 | CIGAR | String | Сжато кодируются события, замеченные в выравнивании (делеции/инсерции/замены) |
| 7 | RNEXT | String | Референсное название следующего закартированного рида. "=" означает, что RNEXT совпадает с RNAME |
| 8 | PNEXT | Int | Позиция картирования следующего рида |
| 9 | TLEN | Int | Расстояние между самой левой и самой правой позициями выравнивания |
| 10 | SEQ | String | Последовательность рида |
| 11 | QUAL | String | Качество прочтения каждого конкретного нкулеодтида (phred33) |
| 12 | Additional field | any |
NM:i - расстояние до генома (количество несовпадающих позиций) NH:i - количество картирований для каждого конкретного рида Общаий формат для дополнительного поля - TAG:TYPE:VALUE, где TAG задается двумя символами [A-Za-z][A-Za-z0-9], TYPE - одна буква, определяющая формат VALUE |
c. Полученный .sam файл весит 16 Гб, что слишком много, поэтому следующим шагом будет его конвертация в бинарный файл
a. Файл .sam был конвертирован в .bam командой
samtools sort -@ 10 -o map_16_chr10.bam map_16_ch10.sam 2> sort_log
b. Полученный .bam файл весит 4.6 Гб, файл .sam больше не нужен, а потому удален
c. Далее .bam файл был проиндексирован с помощью команды
samtools index -@ 10 map_16_chr10.bam 2> index_bam_log
a. Поскольку просто так взять и прочитать бинарный файл несколько затруднительно, можно использовать команду flagstat, которая позволяет собрать статистику по input-файлу и выдать результат
samtools flagstat map_16_chr10.bam > flagstat.txt
b. Теперь в файле flagstat.txt лежит вывод команды.
Картированных чтений обнаружено 14.56%, а в корректно картированных пара чтений картировано 10.75%. Полученные значения можно объяснить тем, что, во-первых, полный экзом картируется лишь на одну хромосому. Во-вторых, ожидается, что чтения из одной пары должны картироваться недалеко друг от друга и быть направлены навстречу друг другу, однако может быть такое, что какое-то чтение из пары могло закартироваться далеко от того места, где следовало бы, поскольку чтения - небольшие сиквенсы, а в геноме короткие участки последовательности могут повторяться. К тому же, чтения могут картироваться на геном не один раз.
Как уже было сказано в предыдущем пункте, если полный экзом картируется лишь на одну хромосому, ожидается множество некартированных чтений, что влечет за собой возникновение желания избавиться от них, оставив только правильно картированные риды.
a. Первый шаг - оставить только картированные на хромосому чтения. Команда ниже позволяет это сделать.
samtools view -h map_16_chr10.bam NC_000010.11 > reads.chr10.sam
NC_000010.11 - название хромосомы, найденное в .fna файле. Далее файл .sam был конвертирован в .bam.
samtools view -bS reads.chr10.sam > reads.chr10.bam
b. Затем с помощью команды
samtools flagstat reads.chr10.bam > 5b_flagstat.txt
был получен файл 5b_flagstat.txt, из которого было получено, что картированных чтений найдено 87.64%, а картированных чтений в корректно картированных парах 64.97%.
c. По сравнению с тем, что было получено в пункте 4b, заметно, что количество закартированных на конкретную (10-ю) хромосому ридов увеличилось, как и ожидалось, однако картированных чтений в корректно картированных парах всё ещё меньше.
d. Несмотря на то, что, казалось бы, раз программа отсекла все некартированные риды, то просто картированных должно быть 100%, однако такого почему-то не наблюдается. Возможно, это связано с тем, что samtools view не идеально распознаёт некартированные чтения из-за повторов, которые могут возникнуть при картировании, а поскольку не задается, какие именно выравнивания надо брать, в выход явно попадают и ошибки. Вероятно, эта программа нужна для того, чтобы отфильтровать точно некартированные риды, потому что после этого можно применять более жесткий фильтр по тому, какие выравнивания брать (в команде из 5е опция -f 0x2 будет выбирать те выравнивания, у которых в поле FLAG стоит соответствующий бит, а 0x2 ставится, если рид сопоставлен правильной паре).
e. Оставить только правильно картированные пары чтений можно с помощью команды:
samtools view -@ 10 -f 0x2 -bS reads.chr10.bam > correct_reads_chr10.bam
f. И снова надо проверить статистику в полученном файле:
samtools flagstat correct_reads_chr10.bam > 5f_flagstat.txt
Из полученного файла 5f_flagstat.txt было найдено, что картированных чтений найдено 100.00%, а картированных чтений в корректно картированных парах 100.00%.
g. При сравнении с результатами в пунке 5b возникает желание радоваться жизни, потому что наконец удалось получить файл, содержащий только картированные чтения корректно картированных пар.
h. Последний полученный .bam файл содержит только правильно спаренные картированыые чтения 10-й хромосомы человека. Данный файл необходимо проиндексировать с помощью команды:
samtools index -b -@ 10 correct_reads_chr10.bam correct_reads_chr10.bai
a. Дублированные чтения были промаркерованы с помощью команды:
picard MarkDuplicates -I correct_reads_chr10.bam -M metrix.file.txt -O marked_dupl.bam
b. Далее к полученному в 6а файлу была применена следующая команда:
samtools flagstat marked_dupl.bam > 6b_flagstat.txt
В результате был получен файл 6b_flagstat.txt, в котором указано, что нашлось 463748 дублированных чтений.