Ресеквенирование. Поиск полиморфизмов у человека

Очистка чтений и анализ их качества

Описание команд и ссылки на результаты доступны в таблице (все необходимые изображения расположены на html странице, ссылка в таблице). В исходном файле до чистки было 10410 последовательностей длиной 40-100 нуклеотидов. В результате чистки с помощью прорамммы Trimmomatic были удалены нуклеотиды с плохим (ниже Q=28, Q=-10lg(вероятности ошибки)) качеством, а также риды короче 50 нуклеотидов. Поэтому на картинке качества сиквенсов "пропали" риды с длинными усами, а коэффициент эксцесса пика распределения качества ридов увеличится, он станет более острым. В результате файл стал содержать 10191 последовательностей. Однако в обоих ридах стоит флажок "failure" по критерию содержания GC на последовательность (ссылки на результаты в виде html страниц смотри в таблице). Этот флажок загорается, если имеются отклонения от нормального распределения содержания GC-пар в более чем 30% ридов. В распределениях, полученных на основе данных моих файлов (ссылка в таблице), имеется сдвиг влево, что является индикатором автоматической ошибки, не зависящей от позиции основания. Кроме того на левой половине пика имеется много остроконечных локальных максимумов, которые служат индикатором проблем с библиотекой. Обычно это обозначает наличие специфического загрязнения, например димеры адаптеров. Но в наших файлах адаптеры уже должны были быть удалены (это видно в поле Adapter Content), поэтому сложно говорить о причинах таких сдвигов. Другое поле, содержащее предупредительные сигналы - распределение длины сиквенсов. В моих файлах пости все последовательности превышают в длину 97 нуклеотидов, однако некоторые содержат меньшую длину. В описании программы указано, что индикаторы из этого раздела можно игнорировать, так как для некоторых секвенаторов наличие ридов неодинаковой длины является нормой. По остальным критериям и изначалный и конечный файл удовлетворяют параметрам качества.
Другими интересными полями выдачи являются проверка распределения нуклеотидов по сиквенсам. В норме все 4 линии должны быть параллельны друг другу и разница в содержании между А и Т и между G и С не должна превышать 10 процентов (правило Чаргоффа). Еще есть поле, описывающее содержание N-неопознанного нуклеотида в ридах. Оно показывает качество секвенирования, насколько прибор смог различить последовательность нуклеотидов. Другим немаловажным показателям качеста секвенирвоания является количество дуплицированных сиквенсов. Маленький уровень дупликация является индикатором высокого уровня покрытия таргетной последовательности, а высокий - ошибки в проведении ПЦР (что почему-то каких-то последовательностей сильно больше, чем других).

Картирование чтений и анализ выравнивания

Используемые команды и программы указаны в таблице 1. Очищенный файл содержал 10191 ридов, из них с помощью программы Hisat2 все были непарными, 10140 (99.50%) выровнены ровно 1 раз, 4 (0.04%) выровнены больше 1 раза и 47 (0.46%) вообще не удалось картировать на хромосому.
С помощью команды samtools depth chr2sorted.bam >depth.tsv я получила табличку с информацией о том, сколько раз был покрыт каждый нуклеотид. C помощью программы excel я отсортировала данные, выбрала нуклеотид, который покрыт много раз (234204021). Ссылка на табличку. Затем с помощью программы genome brouser я определила положение этого нуклеотида на хромосоме. Он попал на ген ATG16L1 . Файл, содержащий экзоны этого гена (в том числе 3' и 5 ' нетранслируемые экзоны доступен по ссылке в формтае .fasta. Для примера был взят экзон chr2:234202903-234202996. еще раз была запущена комнада samtools depth chr2sorted.bam -r chr2:234202903-234202996 >depth2.tsv. Среднее число покрытий 10,21, медиана -11,мода -9 (расчеты в том же файле excel). Также была построена гистограмма распределения глубины покрытия нуклеотида от координаты генома (см рисунок 1). Я думаю, можно утверждать, что покрытие данного экзона довольно равномерное.

Дальше мне стало интересен общий профиль покрытия всей хромосомы. Среднее покрытие составляло 26,7, мода -1, медиана -9. То есть половина нуклеотидов генома покрыта меньше 9 раз. Если взглянуть на гистограмму распределния покрытий по хромосоме видно, что на протяжении всего геноа чередуются участки с высоким и с низким покрытием (см. рисунок 2).

На рисунке 3ниже представлена визуализация участка гена ATG16L1 (chr2:234,180,994-234,183,542). Видно, что бОльшая часть ридов лежит в окрестностях экзонов.

Анализ SNP

Необходимый команды описаны в таблице 1. Описание полиморфизмов содержится в таблице 2.

Таблица 2. Описание полиморфизмов

Координата	Тип полиморфизма	Референс	Чтение	Глубина покрытия	Качество чтения
сhr2:55516588	Замена	G	C	23	184,999
сhr2:55523309	Замена	C	A, G	36	211,003
сhr2:238454154	Замена	C	A	101	225,009

На рисунке 3 представлена визуализация полиморфизмов с помощью IGV. На рисунке B) видно, что программа позволяет визуализировать инсерции в 1-2 нуклеотида (фиолетовые поля), разночтения между ридами, представляет, как легли покрытия на референс. Сначала файл VCF был переведен в формат, пригодны для работы Annovar., в нем не оказалось инделей, только 72 полиморфизма. Первым делом был проведен анализ наличия rs у данных полиморфизмов. Команду смотри в таблице. Необходимая информация содержится в файле rs.hg19_snp138_dropped, то есть те SNP, которые есть в базе dpSNP. Оказалось, что rs есть у 62 нуклеотидов.
Затем я проаннотировала эти SNP в refgene, базе данных, основанной на геномной разметке. Команда указана в таблице. В результате на выход мы получаем три файла: refgene.variant_function содержащий информацию о том, в какое место попал SNP (3' и 5' нетранслируемые области, интроны, экзоны обычных генов, экзоны некодирующих РНК, области, близкие к соединениям спалйсинга, upstream - в пределах 1кб до начала транскрипции, downstream-в пределах 1 кб после сайта окончания транскрипции, интергены), refgene.exonic_variant_function(содержит подробную информацию о том, синонимичная ли замена или нет, какая аминокислота на какую поменялась), rs.log(stdout). В экзоны попало 6, в интроны - 57, ncRNA-exonic -1, UTR3-6, UTR5-2. Гены- ATG16L1,MLPH, CCDC88A, MIR6811. Подробнее о SNP, попавших в экзоны смотри в таблице 3.

Таблица 3. SNP в экзонах

Ген	Координата	АК референс	АК чтение	Нуклеотид референс	Нуклеотид чтение
ATG16L1	сhr2:234183368	T	A	A	G
MLPH	сhr2:238427194	L	P	T	C
MLPH	сhr2:238427251	G	D	G	A
MLPH	сhr2:238443226	H	R	A	G
MLPH	сhr2:238449007	V	A	T	C
MLPH	сhr2:238449107	E	E синонимичная замена	A	G

Затем SNP были проаннотированы с помощью базы данных Gwas, основанная на полногеномном поиске ассоциаций между генотипическими и фенотипическими признаками. Всего нашлось 4 SNP: 55516323 связан с ростом человека, 234173503 и 234183368 с болезнью Крона, а 238443226 с раком простаты.
Чтобы определить, имеют ли SNP этого пациента какую-либо клиническую аннотацию, был произведен поиск по базе данных Clinvar, содержащую инфрмацию о генотипической вариабильности и ее связи с человеческим здоровьем. Оказалось, что нуклеотид с номером 234183368 (замена A на G) связан с заменой 300-ой АК треонина на аланин, идентификатор в этой базе данных:RCV000001189.2. , и это приводит к предрасположенности к Inflammatory bowel disease ( воспаление кишечника), к которым как раз и относится болезнь Крона. Чаще всего заболевание развивается до 30 лет [1]. Данная мутация случается в гене ATG16L1 - autophagy related 16 like 1( Gene ID: 55054), белок, закодированный в этом гене входит в состав большого комплекса, ответственного за аутофагию - процесс, направляющий компоненты внутриклеточного матрикса на деградацию в лизосомы [2]. На рисунке 5 представлена визуализация этого полиморфизма с помощью IGV. Ссылка на сводную табличку с полиморфизмами.

Список литературы

[1] https://en.wikipedia.org/wiki/Inflammatory_bowel_disease [2] https://www.ncbi.nlm.nih.gov/gene/55054