Последние недели 2-3 я довольно плотно занимался изучением возможностей определения STR (коротких тандемных повторов) на основании данных новых технологий высокопроизводительного секвенирования (NGS).
Напомню, что основной способ определения гаплотипов (набора локусов STR) подразумевает использования более традиционных технологий вроде капиллярного электрофореза, ПЦР или пирофореза. Именно так до сих пор типируются гаплотипы Y в научных и коммерческих лабораториях (например, в FTDNA).
Технологии NGS (next generation sequencing), особенно полногеномного сиквенса, были придуманы для других целей, поэтому технически определение STR на уровне условного железа (т.е. с помощью секвенатора) пока не представляется возможным. Поэтому единственное возможное решение — использование особых алгоритмов поиска коротких тандемных повторов в сиквенсе, причем как известных, так и неизвестных. Я не считаю себя дилетантом в области работы с сиквенсами (и их элайнментами), но по мере углубления в материал, я быстро понял всю сложность задачи. Основная сложность — выявление правильной периодичности повторов, т.е. вычисление числа самих повторов. Даже в природе, во время репликации ДНК, полимераза часто произвольно пробуксовывает и дает сбои именно на коротких тандемных повторах, и за счет этого типа мутаций аккумулируется изменчивость (вариативность) этого типа маркеров. То же самое касается и используемых алгоритмов, которые часто ошибаются не в мотиве тандемного повтора, а в числе повторов. Т.е. предположим что мотив повтора состоит из нуклеотидов AGAA. Допустим у человека этот мотив повторяется 12 раз подряд, но программа определяет вместо 12 повторов 11 или, наоборот 13.
Я изучил три программы, созданных для определения STR из данных NGS. Нужно отдать должное чувству юмора их создаталей, ибо названия программы образованы от аббревиатуры STR путем добавления какого-то смыслообразующего корня. Поэтому названия выглядят комично:
lobSTR (http://lobstr.teamerlich.org/)
HipSTR (https://github.com/tfwillems/HipSTR)
GangSTR (https://github.com/gymreklab/GangSTR)
Последную программу я пока так и не смог заставить работать, возможно в ее коде содержится некий баг. Большего успеха я добился с самой известной в списке программой lobSTR и похожей на нее HipSTR. Обе программы показали хорошие тестовые результаты на BAM файлах с парными ридами (paired reads) и высокую корреляцию с данными FTDNA.
Теперь о эксперимента. Для определения аккуратности определяемых этими программами локусов — STR — я взял тестовый BAM файл с сиквенсом Y хромосомы одного из клиентов FTDNA. Поскольку у этого клиента был сделан обычный STR-тест, можно было легко определить аккуратность алгоритма программа путем элементарного сравнения определенных lobSTR/HipSTR локусных значений STR со значениями соответствующих локусов STR, полученных в лаборатории традиционным способом — т.е. PCR и электрофорезом.
К сожалению, выдаваемый клиентам FTDNA bam файл с сиквенсом Y-хромосомы малопригоден в своем изначальном виде для определения STR. Я не знаю в чем дело, но эксеприменты с исходным BAM не дали достоверных результатов. Скорее всего, BAM содержит гибридные риды (парные и одиночные) сиквенса, а также непонятные HipSTR флаги ридов. Видимо, BAM собирался из FASTQ файлов, полученных разными сиквенаторами.
Кроме того, FTDNA или ее партнерская лаборатория, скорее всего использует какой-то кастомный или самописный ассемблер генома — и как следствие, вышеназванные программы очень плохо считывают входящие данные (ибо заточены на работу с BAM файлом сгенерированным классическими ассемблерами вроде BWA, и в меньшей степени, bowtie).
Поэтому пришлось заняться обратной разработкой BAM файла. Сначала я выделил из BAM файла парные риды и экспортировал их в формат FASTQ, а непарные удалил.
Далее я уже следовал рекомендуемой ведущими биоинформатиками процедуре из 12 промежуточных этапов(я не буду описывать все детали, скажу лишь что этот процесс великого делания включает в себя многочисленные фильтровки и рекалибровки нуклеотидных баз собираемого генома).
Пересобранный таким образом геном стал более доступным для нисходящей обработки в lobSTR/HipSTR, и после нескольких неудачных попыток я смог определить значения STR, которые оказались либо идентичными, либо близкими (с разницей в 1-2 повтора) типированным значениям STR.
Вот результы сравнения полученных в HipSTR/lobSTR значений DYS локусов с теми, что содержатся в отчете FTDNA
DYS marker lobSTR HipSTR FTDNA report
DYS389I 13 13 13
DYS389I 13 13 13
DYS389I 13 13 13
DYS389I 13 13 13
DYS390 24 24 24
DYS391 10 10 10
DYS392 11 11 11
DYS393 13 13 13
DYS426 11 11 11
DYS434 9 9 9
DYS435 11 11 11
DYS436 12 12 12
DYS437 15 15 15
DYS438 10 10 10
DYS439 11 12 12
DYS442 17 17 12
DYS444 10 10 10
DYS445 10 12 10
DYS446 13 13 13
DYS454 11 11 11
DYS458 17 17 17
DYS460 10 10 10
DYS461 12 12 12
DYS462 12 12 12
DYS472 8 8 8
DYS485 15 15 14
DYS492 12 12 12
DYS494 9 9 9
DYS511 9 9 9
DYS520 23 23 22
DYS522 12 11 11
DYS531 12 11 11
DYS533 13 13 13
DYS534 12 12 12
DYS537 11 11 11
DYS549 11 11 11
DYS556 11 11 11
DYS565 9 9 9
DYS570 18 18 18
DYS576 16 16 17
DYS578 8 8 8
DYS590 7 7 7
DYS594 10 10 10
DYS607 16 16 12
DYS635 23 23 23
DYS638 11 11 11
DYS641 10 10 10
DYS643 10 10 10
Видно что корреляция между результатами HipSTR и lobSTR выше (0.99) чем попарная корреляция между ними и результатами коммерческого тестирования в FTDNA (0.955 и 0.954). То есть результаты программ чаще согласуются друг с другом, чем с результатами FTDNA.
Обращает внимание то обстоятельства что полученные значения маркеров DYS607 и DYS442 в моем эксперименте существенно отличаются по числу повторов от референсных. Различие 4- 5 повтора. Но тут дело не в ошибке программе, а в разнице использзуемых номенклатур.
DYS442 has had changes in its nomenclature (http://www.hprg.com/hapest5/page2.html). FamilyTreeDNA reports a value 5 units shorter than NIST.
Заметки о генетике 