CC BY
45УДК 577.21:004.02:004.94
Ю.Л. Орлов12, Г.В. Васильев1, Е.В. Кулакова1-2, Н.А. Колчанов1-2
CHIP-SEQ ДАННЫЕ И ИХ АНАЛИЗ Обзорная статья
1Институт цитологии и генетики Сибирского отделения Российской академии наук Российская Федерация, 630090, г. Новосибирск, пр-т ак. Лаврентьева, 10 2Новосибирский государственный университет Российская Федерация, 630090, г. Новосибирск, ул. Пирогова, д. 2.
Секвенирование геномов и исследование регуляции транскрипции
С начала 2000-х годов, после секвенирова-ния первых полных геномов, в молекулярной генетике произошла технологическая революция, связанная с появлением экспрессионных микрочипов высокой плотности и технологий высокопроизводительного секвенирования ДНК. В настоящее время становится возможной детальная функциональная аннотация ре-гуляторных геномных последовательностей по экспериментальным данным, полученным в результате массового параллельного секвенирования ДНК [30]. Кроме определения положения и структуры белок-кодирующих генов, полногеномная аннотация включает описание некодирующих РНК, выделение регуляторных районов генов, исследование хромосомных аномалий и однонуклеотидных полиморфизмов, определение функции белков, предсказание их вторичной и пространственной структуры [9].
Исследование механизмов регуляции транскрипции генов - самостоятельная фундаментальная задача. Определяющую роль в контроле экспрессии генов на уровне транскрипции играют транскрипционные факторы - белки, специфически взаимодействующие с регуля-торными районами генов и белками транскрипционной машины. Экспрессия гена в значительной степени определяется и другими условиями, такими как состояние хроматина в данном районе генома, уровень метилирования ДНК, плотность нуклеосомной упаковки ДНК [14]. Компьютерные модели транскрипции необходимы как для успешного предсказания особенностей экспрессии генов, так и для выполнения прикладных исследований, например, реконструкции регуляторных сетей, создания генетических конструкций с задан-
ными свойствами, исследования механизмов заболеваний, канцерогенеза, поиска мишеней для лекарств, токсикологических исследований, выявления ключевых биомаркеров [30].
Реконструировать генные сети (регуля-торные контуры) взаимодействующих генов транскрипционных факторов можно через определение их сайтов связывания в промоторах генов - мишеней воздействия этих транскрипционных факторов. Регуляторные районы групп коэкспрессирующихся генов зачастую имеют общие черты организации, что выражается в наличии регуляторных паттернов (цис-регуляторных модулей), состоящих из устойчивых сочетаний сайтов связывания нескольких транскрипционных факторов различных типов и других нуклеотидных мотивов [4]. Выявление и анализ таких регуляторных паттернов является основой для построения обобщенных моделей регуляторных районов группы коэкспрессирующихся генов.
В настоящее время накоплен колоссальный объем данных в области регуляции экспрессии генов эукариот, наблюдается их непрерывный рост как первичных данных, доступных в частности в репозитории данных GEO NCBI (Gene Expression Omnibus) (www.ncbi. nlm.nih.gov/geo/). Все большую актуальность приобретает формализация описания механизмов регуляции транскрипции и разработка на этой основе методов интеграции гетерогенной информации об особенностях регуляции экспрессии генов, использующая как технологии определения сайтов связывания транскрипционных факторов в геноме, так и оценки уровней экспрессии генов на микрочипах.
Период после завершения проекта «Геном человека» принято называть «постгеномной эрой». Черновой вариант генома был опубликован в 2001 г., окончательным завершением
проекта можно считать опубликование первичной структуры последней из секвениро-ванных хромосом - хромосомы I [11]. Однако действительно фундаментальные изменения в молекулярной биологии были внесены не столько данными о нуклеотидных последовательностях первых секвенированных геномов, сколько появлением новых революционных методов секвенирования ДНК. Первые геномные проекты выполнялись с использованием капиллярных секвенаторов, т.е. технологии первого поколения. Секвенирование по Сэнгеру дает последовательности очень высокого качества и достаточно длинные (500-800 п.н.) для эффективной сборки в протяженные контиги, однако требует очень высоких затрат времени, денег и высококвалифицированного труда. Радикальное удешевление секвенирования ДНК явилось результатом появления приборов второго поколения, основанных на принципах массового параллельного секвенирования. Совершенствование этих систем привело к огромному повышению производительности приборов и одновременно фантастическому падению стоимости их работы. Так, расчетная стоимость ре-секвенирования человеческого генома, по данным National Human Genome Research Institute (http://genome.gov/sequencingcosts), за период с 2001 по 2011 год упала со 100 миллионов долларов до 10 тысяч, т.е. в среднем снижалась в 1000 раз за год. Заметим, что данная стоимость рассчитана для крупных геномных центров, фактически представляющих из себя хорошо организованное индустриальное производство. К 2014 году эта стоимость составила 3-5 тысяч долларов. Доступность секвенирования, помимо огромного вклада в фундаментальные научные исследования, позволяет использовать секвенирование индивидуального генома как стандартный диагностический тест, что исключительно важно для медицинских приложений.
Ключевой особенностью технологий второго поколения является получение молекулярных колоний на твердофазном носителе с помощью ПЦР с одной молекулы и дальнейшее параллельное одновременное секвенирование сотен миллионов и миллиардов этих кластеров в ходе синтеза новых цепей. В настоящий момент существуют четыре технологии второго поколения, каждая из них реализована в нескольких вариантах приборов. Технология 454 от фирмы
Roche основана на пиросеквенировании с оптической детекцией света, испускаемого люцифе-разой. SOLiD от Applied Biosystems выполняет секвенирование путем последовательного ли-гирования флюоресцентно меченых олигону-клеотидов. Технология Illumina основана на использовании обратимой терминации синтеза флюоресцентно мечеными нуклеотидами. Прибор IonTorrent от Applied Biosystems работает на принципе полупроводникового секвениро-вания с детекцией изменений рН при синтезе новых цепей ДНК внутри микролунок чипа. Эти технологии обеспечивают секвенирование достаточно коротких фрагментов ДНК - 75 п.н. у SOLiD, 100-300 у Illumina HiSeq, 200-400 для IonTorrent, 400-700 п.н. для Roche FLX Titanium [3]. Огромные объёмы данных, представленные короткими последовательностями, ставят технически более сложные задачи биоинформационного анализа, выравнивания, картирования и ассемблирования коротких последовательностей ДНК. Будет справедливо сказать, что на сегодняшний день корректный анализ данных секвенирования ДНК является большей проблемой, чем их получение.
Серьёзные надежды в течение десятилетия возлагались на разработку третьего поколения приборов - «одномолекулярных» секвенаторов. Их ключевой особенностью является отсутствие ПЦР в процессе пробоподготовки, при этом секвенируется каждая единичная молекула ДНК, а не кластеры идентичных копий. Это является достоинством метода, поскольку именно наличие ПЦР во время пробоподготовки вызывает большую часть системных ошибок, характерных для секвена-торов второго поколения. За последние годы было анонсировано несколько десятков фирм, готовых предложить работающую модель одномолекулярного секвенатора в самое ближайшее время. Однако эти обещания так и не претворились в реальность. Представителем реально работающих в настоящее время одно-молекулярных секвенаторов является прибор фирмы Pacific Biosciences (PacBio), основанный на детекции включения в синтезируемую цепь флуоресцентно меченых нуклеотидов. Его особенностью является возможность непрерывного чтения молекул ДНК длиной до нескольких тысяч нуклеотидов при крайне невысокой его точности - более 10% ошибок
[28]. В 2014 году продемонстрирована работоспособность в реальных условиях прибора Oxford Technologies Nanopore, основанного на прямой детекции оснований ДНК при проходе через нанопоры в электронном чипе и позволяющего сверхдлинные чтения (до 50 т.п.н.). Для использования в геномных проектах любая технология должна обладать тремя условиями: высокая производительность, низкая себестоимость работы и высокая точность прочтения. К сожалению, ни один одномоле-кулярный секвенатор до сих пор не обладает одновременно всеми тремя достоинствами и поэтому на настоящий момент безнадежно проигрывает приборам второго поколения. Однако большая длина чтения (десятки т.п.н.) и отсутствие ПЦР делают эти технологии крайне полезными для узкого круга задач.
Технология иммунопреципитации хроматина
Среди методов, позволяющих изучать связывание транскрипционных факторов (ТФ) с ДНК, наибольшими перспективами для полногеномных исследований обладает метод иммунопреципитации хроматина (Chromatin ImmunoPrecipitation - ChIP) с последующим секвенированием, ChIP-seq.
Ядро клетки
Рис.
Процедура ChIP-seq на первом этапе требует фиксации контактирующих молекул ДНК и белков в клетке с помощью формальдегида, что вызывает образование ковалентных сшивок между ДНК и белками. Далее хроматин дробится ультразвуком на фрагменты (существует термин «соникация» или, дословно «озвучивание» хроматина). Альтернативно разделение ДНК может выполняться с помощью разрезания ферментами рестрикции. Затем с помощью иммунопреципитации со специфическими антителами выделяется ДНК, с которыми физически связан интересующий исследователя белок. Белковая фракция отмывается, а ДНК (относительно короткие фрагменты не более нескольких сотен оснований) направляется на секвенирование с помощью имеющегося оборудования массового параллельного секвени-рования ДНК (Roche 454, Illumina или SOLiD).
Прочитывание ДНК выполняется не для всей последовательности экстрагированного фрагмента, составляющей до нескольких сотен нуклеотидов, а для крайнего участка -от 20 до 75 п.о. Если выполняется лигирова-ние концов фрагмента ДНК, они могут секве-нироваться попарно (так называемые парные концы, или PET - Paired End Tags). На рис. 1 представлена схема определения кластеров
Фиксация формальдегидом
Дефрагментация хроматина
ТГмкГТ
Иммунопреципитации связных фрагментов ДНК антителами
Секвенирование
y^V-iTÏTTTinr
Определение кластеров фрагментов ДНК на хромосоме
5' :
Расчет ChIP профиля, определение пиков
3'
Картирование на
референсный геном
Пики профиля размером 3 и 4
Си игл тон
yd
J
L
= 3'
"Y"
Общий размах пика СЫР-$еч
1. Схема иммунопреципитации хроматина, картирования фрагментов ДНК и пример получаемого ступенчатого профиля CЫP-seq
фрагментов ДНК на хромосоме для парных фрагментов (парных концов метода ChIP-PET) и для одиночных фрагментов. Задача картирования секвенированных фрагментов (прочтений, или «ридов» (от. англ. reads)) технически выполняется с помощью выравнивания и быстрого поиска совпадений с последовательностями ДНК хромосом генома.
Картирование секвенированных последовательностей технически требует до нескольких часов машинного времени на персональном компьютере для достаточно больших по размеру геномов эукариот. Для ее решения применяются компьютерные методы оптимизации выравнивания и быстрого поиска совпадений в нуклеотидных последовательностях. Существует набор программ картирования прочтений для форматов данных Illumina и форматов цветовой кодировки SOLiD (см. обзор программ на SEQanswers, http://seqanswers. com). Далее, используя координаты секвени-рованных фрагментов на хромосомах рефе-
ренсного генома, строится так называемый геномный профиль и выполняется его анализ (рис. 1 и рис. 2).
При построении профиля связывания координаты прочтений удлиняются до размера исходного фрагмента ДНК (150-200 нт) и получается «лестница» или «ступени» фрагментов, наложенных друг на друга в геномных координатах. При использовании парных концов длина фрагментов определяется точно, при одностороннем секвенировании длина фрагмента оценивается по среднему значению (например, 150 нт). Из наложенных друг на друга фрагментов строится ступенчатый профиль. Затем определяются пики такого геномного профиля. Пик - это наиболее высокая точка локального профиля. Пик более вероятно содержит сайт связывания транскрипционного фактора, чем окружающий геномный район, поскольку все фрагменты ДНК, составляющие «ступени» пика, должны быть связаны с белком в СЫР-эксперименте.
АСАСААССССССССААСОСС СААСАТСТСССААСАТСТТС
Рис. 2. Полногеномное представление сайтов связывания транскрипционного фактора р53 на хромосоме 6 человека (вверху) и выделенный стрелкой район генома, содержащий ген CDKN1A и его 5' район в большем масштабе. Два участка профиля, сформированного фрагментами ДНК, прошедшими иммунопреципитацию, образуют пики, содержащие узнаваемый консенсусный мотив связывания р53 AGACATGCCCAGACATGCCC. Внизу показана частотная матрица мотива связывания, восстановленная по полногеномным данным [31]
Исторически метод иммунопреципитации хроматина предназначался для анализа взаимодействий белок-ДНК на единичной или ограниченной выборке данных (несколько промоторных участков). Массовое использование олигонуклеотидных микрочипов позволило усовершенствовать технологию и получать гибридизационный сигнал связывания исследуемой последовательности с заранее подготовленными пробами. Технология получила название ChIP-on-chip, или ChIP-chip (то есть хроматин-иммунопреципитация на микрочипе) [7, 12]. Такой анализ возможен для достаточно больших наборов проб, включая, например, отдельные хромосомы или все про-моторные районы генов, известные в геноме. Существуют варианты микрочиповых технологий для определения связывания с белками в специализированных вариантах - метод DamID (использующий белок Dam - DNA adenine methyltransferase), метод DIP-seq (DNA ImmunoPrecipitation) для исследования связывания с ДНК без хроматина [26].
Общий набор методов выделения сайтов связывания, основанных на иммунопреци-питации, по-английски называют "ChIP" [7]. Технология ChIP-seq применялась для поиска сайтов связывания ТФ мыши [4], человека [5, 31, 32], анализа модификаций гистонов в различных тканях [14].
Исследование профиля ChIP-seq требует специализированных компьютерных программ и ориентированных на конкретную задачу методов, в зависимости от технологий секвенирования (коррекция на специфические ошибки), размера и особенностей генома (наличие повторенных последовательностей, детали аннотации - от компактного, хорошо изученного генома дрожжей до большого по размеру генома человека и, например, геномов растений, где нет референсной последовательности). Первичные данные секвениро-вания поступают в формате bed-файлов либо в FASTA формате для картирования на геном (тысячи и миллионы последовательностей).
В результате стандартного эксперимента ChIP-seq [4] по извлечению белка, специфически связанного с ДНК, получается набор последовательностей, упорядоченное расположение которых в геноме (в хромосомных координатах) дает ступенчатый профиль. Короткие
секвенированные фрагменты уже на геномной карте удлиняют до 150-300 нуклеотидов, что в среднем соответствует размеру экстрагированной в ChIP-seq ДНК (и соответствует размеру одной-двух нуклеосом). Заметим, что практически достаточно секвенировать первые 25-35 нуклеотидов, чтобы однозначно найти координаты фрагмента в геноме.
Образующийся численный профиль связывания вдоль хромосомы (одномерная координата) содержит интересующие исследователя места в геноме - участки (сайты) связывания с белком. Их видно (например, при увеличении масштаба в геномном браузере) как пики (горки) профиля (рис. 2). Пики получаются в результате перекрывания нескольких фрагментов (графически как наложенные друг на друга полоски в одномерных координатах хромосомы). Пик профиля как сигнал можно математически описать координатами начала и конца участка, указать высоту пика, положение вершины (саммита), оценить вероятность получения пика по случайным причинам. При близком расположении нескольких участков, связанных с исследуемым белком, пик может иметь несколько вершин (мультимодальный пик), что задает неоднозначность при выборе позиции специфичного сайта.
Секвенирование фрагментов ДНК может выполняться с двух сторон с использованием технологии секвенирования парных концов, что позволяет более точно картировать сайты [6]. Отметим, что относительно недавние работы 2005-2007 годов использовавшие секвенирование парных концов (PET - Paired End Tags), что вошло в название ChIP-PET, включали клонирование выделенных фрагментов ДНК как необходимый технологический шаг перед секвенированием [31, 32]. Клонирование увеличивало общее время эксперимента и вело к частичной потере данных (при неравномерном клонировании фрагментов ДНК по длине). В настоящее время метод ChIP-seq, основанный на прямом секвенировании полученных после иммунопреципитации фрагментов ДНК, включает возможность использования протокола секвенирования парных концов [6, 22].
В целом, секвенирование обладает рядом преимуществ по сравнению с микрочиповыми технологиями исследования связывания белков в регуляторных районах ChIP-on-chip [12].
Прежде всего, результатом ChIP-seq является картина полногеномного распределения сайтов связывания транскрипционных факторов (ССТФ). Кроме того, полученный результат свободен от предварительной селекции исследуемых районов, которые могут существенно искажать конечный результат. Результатом ChIP-seq являются не только относительные уровни сигнала связывания, а конкретные участки последовательностей ДНК, что дает больше возможностей для теоретического анализа и точного определения сайтов.
Метод иммунопреципитации хроматина с последующим секвенированием широко используется для получения картины полногеномного распределения сайтов связывания различных транскрипционных факторов [4, 5]. Более 160 профилей связывания для генома человека представлены в данных проекта ENCODE [9], десятки профилей ChIP-seq доступны для генома мыши, других модельных организмов.
Обработка найденного набора сайтов связывания исследуемого транскрипционного фактора, найденных в геноме при помощи ChIP-seq, требует разработки компьютерных программ сравнения положения таких сайтов и расположения генов в геноме, определения потенциальных генов мишеней действия этого транскрипционного фактора [6]. Еще более сложной становится задача описания регуля-торных районов генов и определения генов-мишеней при анализе данных трехмерных контактов. Данными являются уже не линейные координаты последовательностей, а двумерные - наборы пар контактов и матрицы контактов.
Математическая задача анализа профиля секвенирования
В типичной задаче анализа профиля ChIP-seq геном человека (размер около 3 Гб) разделен по хромосомам от 20 до 200 Мб; профили строятся для каждой хромосомы, стандартно каждой позиции нуклеотида ставится в соответствие высота профиля, линейно увеличивая размер файла в зависимости от размера генома. Однако при более компактной записи профиля, размер файла чуть меньше, порядка 100 Мб. При дополнительной аннотации размеры файлов значительно увеличиваются,
фактически на пределе возможностей работы стандартного персонального компьютера. Для других организмов геном может иметь как меньший (дрожжи, нематода), так и существенно больший размер (геномы растений, таких как пшеница); практически объем данных ChIP-seq составляет от 100 Мб до 1 Гб.
В результате секвенирования ChIP-seq имеем набор коротких фрагментов (прочтений, или «ридов»), размером от 25 до 100 нуклеотидов. Эти фрагменты распределены по геному неравномерно. Размер библиотек (наборов данных секвенирования, получаемых в одном эксперименте) колеблется от 2 до 40 миллионов таких коротких последовательностей. Для анализа мотива связывания транскрипционного фактора используют нуклеотидные последовательности пиков геномного профиля (см. рис. 3), сравнение с известными базами данных сайтов связывания [19].
Уникальность (однозначность) картирования короткой последовательности на геном представляет особую проблему анализа данных. Если в геноме был повтор (два участка ДНК в разных местах генома, достаточно длинных, скажем, 100 нт и более), и наш короткий фрагмент ДНК попал в повтор, то мы не можем его однозначно (уникально) картировать. Пример таких затруднений - картирование регуляторных сайтов для генов, имеющих псевдогены. Общее свойство протяженных геномных последовательностей -иметь участки, недоступные для однозначного картирования короткими фрагментами из-за длинных совершенных повторов. Существуют термины «картируемость» (тарраЬШу) и «уникома» (uniqueome) - как часть генома, которая однозначно (уникально) определяется короткими фрагментами заданной длины [21]. Для каждой длины секвенированных ридов существует своя «уникома» - например, для фрагментов размером 50 нт некартируемых участков гораздо меньше, чем для фрагментов 25 нт. Наборы таких карт можно рассчитать, существуют и готовые разметки «уникальности» для нескольких референсных геномов, в частности геномов человека и мыши [21]. Отметим, что программы разметки «уникальности» требуют высокопроизводительных компьютерных вычислений - фактически, поиска всех повторов в геноме.
Рис. 3. Расположение секвенированных парных концов СЫР-РЕТ и анализ нуклеотидной последовательности связывания Мус. Известный сайт связывания Мус в первом интроне гена NPM1 определен кластером РЕТ-4. Показаны перекрывающиеся РЕТ-фрагменты (вверху), образующийся ступенчатый профиль (пик) пересечения фрагментов, координаты в геноме человека. Внизу показана 86 нт последовательность, соответствующая пересечению всех фрагментов ДНК в кластере с выделенными элементами CACGТG, соответствующими каноническому Е-боксу [32]
Поскольку эксперимент ChIP-seq проводится на пуле (большом наборе) клеток, каждый акт связывания, представленный фрагментом ДНК, не зависит от других последовательностей. Перекрывание фрагментов в геноме в районе пика подтверждает специфичность связывания. Истинные (специфичные) пики профиля имеют большую высоту (сформированы большим числом фрагментов ДНК), что можно показать статистически. Для определения порогового значения пика, когда с уверенностью можно говорить про специфичность, удобно упорядочить все полученные пики по высоте и рассчитать число пиков в зависимости от высоты (рис. 4). Полученное распределение всегда очень неравномерно - большинство пиков имеет небольшую высоту, затем число высоких пиков быстро (почти экспоненциально) убывает.
Задача определения набора пиков в геноме решается статистически с помощью сравнения экспериментально полученного распределения набора пиков к ожидаемому по случайным причинам (рис. 4). Был предложен компьютерный алгоритм определения пиков профиля и их последующей фильтрации от «шума» и ошибок
секвенирования, основанный на статистике распределения числа сайтов в геноме в зависимости от высоты пика ChIP-seq и сравнении специфического профиля с контрольным. Контрольное секвенирование выполняется после разделения геномной ДНК ультразвуком без антител, или с иммунопреципитацией к неспецифическому белку - GFP или ^О. Такой подход в разных вариантах реализован в компьютерных программах определения пиков профиля ChIP-seq [33].
Распределение вероятности наблюдения РоЫ(Х=т) сайтов (пиков профиля) X фиксированной высоты т может быть представлено как взвешенная сумма специфического и неспецифического распределений:
РоЪ(Х=т) = а*Рр(Х=т) + (1 - а )*Рт(Х=т),
где РоЫ - функция вероятности наблюдаемого распределения встречаемости ЗДШ пиков; X - пики;
т=1,2,3,... - высота пика; Р - вероятность специфичных пиков; 0 < а < 1 - доля специфичных пиков в общем распределении пиков в геноме;
Рис. 4. Распределение числа сайтов в геноме в зависимости от высоты пика ChIP-seq (для связывания ТФ Nanog
мыши) [4]
Pm - вероятность неспецифичных пиков (шумовой сигнал) в общем распределении пиков по профилю в геноме [18].
На основе экспериментальных СЫР-данных мы можем построить эмпирическую функцию распределения пиков в геномном профиле. Распределение числа специфичных пиков Psp в зависимости от их высоты может быть промоделировано с помощью распределения Пуассона или распределения Парето, начиная с большой высоты пика, где шумового сигнала уже нет (например, с m=10 на рис. 4) [18]. Неспецифичное распределение Pns может быть оценено с помощью компьютерной симуляции по появлению кластеров (пиков) при случайном (равномерном) распределении прочтений (секвенированных фрагментов ДНК) вдоль хромосомы. Разработанная компьютерная модель симуляции случайных пиков в геноме требует достаточно интенсивных пересчетов (несколько часов работы на персональном компьютере). Позднее в качестве контроля использовалось число пиков неспецифичного секвенирования (без иммунопреципитации) всей доступной ДНК или с использованием неспецифичных белков (таких как GFP, IgG). Специфичность связывания P зависит от силы связывания ДНК с белком, что может быть подтверждено далее с помощью выборочного тестирования методом qPCR [14].
Для опубликованных ТФ в геноме получается от 1000 до 20 000 сайтов связывания, подтвержденных (выборочно) независимым тестированием [4, 14]. Стандартно получалось около 5000 сайтов (мест на хромосоме), каждый со своей высотой пика, характеризующей его «силу связывания». Сила связывания определяется по короткой последовательности ДНК - те самые 6-8 нт мотива связывания ТФ. Показано, что высота пика профиля ChIP-seq в геномном эксперименте должна коррелировать со связыванием в отдельных проверочных экспериментах qPCR.
Наиболее высокие пики геномного профиля, как правило, содержат мотив сайта связывания - последовательность, похожую на известную консенсусную последовательность ДНК (рис. 2). Мотив может допускать несовпадения в нуклеотидной последовательности, описываемые позиционной весовой матрицей для данного сайта. Такие весовые матрицы для многих транскрипционных факторов есть в базах данных, таких как TRANSFAC, TRRD (http://wwwmgs.bionet.nsc.ru/mgs/gnw/trrd/), JASPAR (http://jaspar.genereg.net/) (результаты «одиночных» разрозненных экспериментов или более старых технологий, собранные по литературе). Высокий процент содержания в пиках ChIP-seq исследуемого мотива, оцененного весовой матрицей, показывает, что положение пиков определено правильно.
Другая оценка качества сигнала связывания в профиле ChIP-seq - это относительная высота пика по отношению к контролю - отношение числа специфических фрагментов ДНК в исследуемой точке генома к числу неспецифических, полученном в контрольном эксперименте секвенирования без им-мунопреципитации [2]. Контрольный профиль получают в результате эксперимента по секвенированию в тех же условиях, но без специфических белков («пустой» прогон). Обычно требуют превышения высоты пика по отношению к контролю в 3-5 раз. Используя контрольный профиль секвенирова-ния можно получить оценку ошибки ложного предсказания для каждого пика, FDR (False Discovery Rate).
Есть ряд опубликованных программ анализа пиков ChIP-seq, например GLItR (GLobal Identifier of Target Regions), MACS [33], HPeak, PeakFinder, GLITR, QuEST, CisGenome, USeq и PICS (см. для обзора [2, 19]).
Для картирования длинных фрагментов ДНК хорошо подходят программы MUMmer и BLAT. Для картирования коротких фрагментов ДНК набор программ достаточно велик: MAQ (Mapping and Assembly with Quality), SOAP (Short Oligonucleotide Alignment) [24], использующая индексы для представления референсной последовательности; программа ELAND, ориентированная на стандарт данных Illumina (http://www.illumina.com/ systems.ilmn).
Распространены также программы Bowtie, использующая индекс Burrows-Wheeler индекс (BWT) [20], SeqMap, RMAP, ZOOM [2]. Сегодня необходимо развитие методов для более широкого их использования для геномов с недостаточной аннотацией.
В целом, задачи компьютерной обработки полногеномных данных секвенирования можно разделить на следующие направления:
• первичная фильтрация данных, картирование;
• анализ профилей секвенирования ДНК, сопряженного с иммунопреципитацией (ChIP-seq) с выделением как точечных участков связывания, так и протяженных участков (модификации гистонов);
• разметка сайтов связывания нуклеосом по нуклеотидной последовательности;
• интеграция данных секвенирования между собой (кластеры сайтов связывания транскрипционных факторов);
• интеграция данных с геномной аннотацией (определение генов-мишеней).
Программа HPeak (Hidden Markov model Peak) использует скрытые марковские модели для определения геномных участков, контактирующих с белками. Используется статистическая модель, учитывающая реалистичное распределение вероятностей для прочтений ДНК.
Программа MACS [33] исходно была ориентирована на технологию Illumina Solexa (данные Solexa Genome Analyzer). Программа использует фрагменты («риды») в противоположных ориентациях, чтобы определить так называемый размер сдвига - близость между «ридами», содержащими сайты связывания. Преимуществом MACS является локальное моделирование «шумового», или контрольного секвенирования с помощью распределения Пуассона по участкам хромосом (участки, размер которых задается пользователем, тысяча -десять тысяч пар нуклеотидов).
Рассмотрим пример серии экспериментов ChIP-seq для одного и то же типа клеток. Масштабное исследование профилей связывания 13 различных транскрипционных факторов (ТФ) в геноме мыши было выполнено в работе [4] на эмбриональных стволовых клетках. Исследовались ТФ, значимые для поддержания плюрипотентности и развития организма, такие как Oct4, Nanog, Sox2. Использовалась платформа секвенирования Illumina Solexa. Для каждого транскрипционного фактора было получено от 5 до 12 миллионов последовательностей, в среднем 65% из них картировались однозначно на референсный геном мыши, что соответствует стандартам экспериментов ChIP-seq. Для контроля использовалась иммунопреципитация к неспецифическому белку GFP (Green Fluorescent Protein). Определение значимых кластеров (пиков) профиля выполняется с помощью сравнения профилей с контрольным секве-нированием и последующей нормализацией. Обычно получается несколько тысяч сайтов на геном для одного транскрипционного фактора. Так, в работе Chen et al. было определено от одной до 40 тысяч сайтов связывания для этих 13 ТФ [4].
Рис. 5. Профили связывания 13 различных транскрипционных факторов в геноме мыши для гена Pou5f1. Внизу представлен профиль контрольного секвениро-вания (для неспецифического белка вБР) [4]
Для сильно «зашумленных» данных или недостаточно глубокого секвенирования могут быть специфические статистические проблемы. Например, при неглубоком секвенировании (0,5-1 миллион ридов вместо 5-10 миллионов),
геном человека не покрыт полностью секвени-рованными фрагментами и трудно гарантировать, что все сайты в эксперименте выявлены. Ранее считалось (эмпирически), что даже 5 миллионов фрагментов уже достаточно, чтобы выявить все сайты в геноме [4]. Сейчас стандарт эксперимента ChIP-seq - это 20-40 миллионов прочтений ДНК. Можно показать статистически, что при дальнейшем увеличении глубины секвенирования просто повышается порог высоты пика при выделении пиков из профиля, и уже не детектируются новые сайты в геноме (в данных условиях эксперимента, фиксированного антитела), что обосновывает достаточность использования меньшего числа прочтений [4].
Следует отметить, что описанные возможности применимы для геномов хорошо изученных организмов, геномы которых не только известны, но и хорошо аннотированы (человек, мышь, A. thaliana, нематода C. elegance и др.). Однако существует большое число важных для изучения организмов, геномы которых могут быть плохо аннотированы (рыбка D. rerio) или вообще отсутствовать (слабо изученные геномы, например, паразитические черви). Тем не менее, для модельных медицинских и биотехнологических задач изучение этих организмов очень актуально, и эксперименты ChIP-seq выполнимы.
Полногеномное секвенирование коротких фрагментов без специфического белка может использоваться и для определения нуклеосом-ной упаковки (позиционирования всех нуклео-сом) в полном геноме эукариот [15]. Исследовалась, в частности, нуклеосомная упаковка в геноме дрожжей и ее взаимодействие с сайтами связывания транскрипционных факторов. При этом анализ профиля секвенирования ДНК (прямое секвенирование нуклеосомных фрагментов без иммунопреципитации) интегрировался с данными о сайтах связывания, определенными первоначально посредством ChIP-chip технологии [15].
Секвенирование и трехмерная структура хромосом
Проблема анализа трехмерной структуры генома активно исследуется различными методами, использующими геномное секвениро-вание. Интересно отметить технологии Hi-C [8] и ChIA-PET [10], позволяющие получать
информацию о пространственной структуре хромосом через секвенирование. Метод Hi-C [8] клетки без специфических взаимодействий, не используя иммунопреципитацию. Метод показывает, как молекула ДНК, линейный размер которой многократно превышает размеры интерфазного ядра, компактизуется и помещается в малый объем ядра. Hi-C основан на технологии 3C (Chromosome Conformation Capture) и технологиях массового параллельного секве-нирования. 3C - метод, который позволяет исследовать пространственное взаимодействие двух специфичных районов генома. Методика измеряет вероятность взаимодействия двух районов генома в большой популяции клеток (105-106). На первом этапе формальдегидом фиксируют белок-белковые, белок-ДНК и белок-РНК-взаимодействия за счет образования ковалентных связей. Затем ДНК фрагмен-тируется с помощью рестриктаз. После этого проводят лигирование (соединения концов двухцепочечной ДНК) в условия разбавления. В таких условиях лигируются только концы молекул ДНК, сближенных в пространстве. Далее нужные молекулы выделяются (обычно с помощью мечения биотином). Затем такие лиги-рованные концы секвенируются, и полученные парные последовательности ДНК картируются на геном, аналогично описанным процедурам ChIP-seq. Создается библиотека (большой набор) попарно взаимодействующих молекул ДНК. Исследование таких парных контактов представляет новую задачу биоинформатики.
Исследование Hi-C показало, что укладка молекул ДНК в ядре достигается путем упаковки по модели «фрактальной глобулы». Такая модель позволяет избежать большого числа сложных узлов. Однако в значительной степени остаются неисследованными важнейшие вопросы о хромосомных контактах промоторов генов в трехмерном пространстве ядра клетки и расположения функциональных групп генов в хромосомных доменах, которые важны в контексте аннотации и исследования генных сетей.
Рассматриваются пространственные домены на хромосоме. Такие домены выделяются при обработке данных Hi-C, и уже представлены в ряде интернет-доступных баз данных (например 3DGD [29] и Mouse Encode Project at Ren Lab [27]).
Метод ChIA-PET (Chromatin Immunopre-cipitation Analysis - Paired End Tags) [10, 23], использующий иммунопреципитацию хроматина, позволяет определять контактирующие участки хромосом, контакты которых опосредованы белками или белковыми комплексами. Отличие от метода Hi-C состоит в том, что в процедуру обработки хроматиновых контактов добавлена стадия иммунопреципитации, так же как в методе ChIP-seq. Таким образом, определяется структура хромосомных контактов, опосредованных специфическим белком (например, транскрипционным фактором ERa или фактором CTCF), а не все контакты хромосом. Метод ChIA-PET предназначен для исследования хромосомных петель, прежде всего между регуляторными районами генов, контактов между удаленными энхансерами и промоторами. Контактирующие участки хромосом выделяются более специфично, в связи с исследуемым белком - транскрипционным фактором, поэтому этот метод можно считать расширением метода ChIP-seq.
Для белка ERa - рецептора эстрогенов -выполнялся анализ карт контактов на хромосомах по методу ChIA-PET [10]. Данные хромосомных контактов, полученные с помощью полногеномного секвенирования, независимо подтверждались с помощью экспериментов по технологии 3C (Chromosome conformation capture) и флуоресцентной гибридизации in situ (FISH).
Эксперименты ChIA-PET по определению хромосомных контактов в ядре клетки, опосредованных уже не отдельным транскрипционным фактором, а целым транскрипционным комплексом РНК полимеразы II представлены в работе [23]. Был проведен компьютерный анализ распределения хромосомных контактов относительно генов и относительно участков генома, ассоциированных с модификациями хроматина (модификации гистона H3, связанные с открытым состоянием хроматина) [17]. Предложена классификация моделей промоторных, энха-серных и мультигенных контактов, опосредованных комплексом РНК полимеразы II [23], включающая классы: базальный промотор, промотор-энхансер и мультигенный ло-кус. Модель базального промотора включает только локальные петли ДНК в промоторе,
без удаленных взаимодействий. Модель одиночного гена включает только петли в районе гена - между энхансером и промотором, возможно между 5' и 3' районами гена, но без других белок-кодирующих генов. И, наконец, мультигенная модель включает сразу несколько генов, расположенных рядом друг с другом и контактирующих промоторными районами. Введен термин «хромоперон» (chromoperon) или «хромосомный оперон». В этой модели также возможен контакт промоторов с удаленными энхансерами.
Для связывания транскрипционных факторов в геноме важно состояние хроматина, его открытость (отсутствие нуклеосомной упаковки), доступность ДНК. Маркеры модификаций гистонов, прежде всего гистона Н3, модификации лизина в позициях 4, 14, 36, включающие метилирование и ацетилиро-вания, связаны с доступностью ДНК и могут быть определены экспериментально также по технологии ChIP-seq [14]. Связь хромосомных контактов и модификаций хроматина (метилирование и ацетилирование) гистонов (гистона Н3) исследовалась статистически
в масштабе генома [23]. Исследование контактов хромосом, опосредованных комплексом РНК полимеразы RNA Pol II в культурах клеток человека, подтвердило ассоциацию с сайтами связывания факторов транскрипции и модификациями хроматина [23]. Было показано статистически, что хроматин в контактирующих участках открыт - гистоны имеют маркеры модификации активной транскрипции - H3K4me3, H3K9ac (соответствуют пикам профиля ChIP-seq), в то же время модификации репрессии транскрипции H3K27me3 не имеют пиков.
Экспериментально полученные профили модификаций хроматина позволяют предсказывать сайты связывания транскрипционных факторов в полногеномной шкале, что было показано детально для связывания рецептора эстрогенов ERa [14]. Более того, уже методом ChIA-PET показано, что ассоциации хромосомных контактов, опосредованных белком ERa, также связаны с открытым (активированным) состоянием хроматина, в частности модификациями H3K4me3, H3K4me1 [10] в линии клеток MCF-7.
Рис. 6. ChIA-PET. Модели промоторных, энхасерных и мультигенных контактов, опосредованных комплексом
РНК полимеразы II [23]
Заключение и обсуждение
Технологии секвенирования, сопряженные с иммунопреципитацией хроматина (ChIP), позволяют исследовать механизмы регуляции транскрипции в масштабе генома эукариот [7, 30]. В данном обзоре кратко представлены особенности моделирования такого рода полногеномных данных, связанные с анализом профилей ChIP-seq, выделением сайтов связывания и дальнейшим анализом их ну-клеотидных последовательностей [4, 14, 23]. Технология ChIP-seq позволяет исследовать профили связывания различных транскрипционных факторов и ко-факторов, сопоставлять карты их взаимодействий с профилями модификаций хроматина и нуклеосомной упаковки [14]. Расширение технологии секвенирования ДНК на пространственно-контактирующие участки хромосом (ChIA-PET) позволяет изучать уже пространственную организацию генома в ядре клетки, также в связи с регуляцией транскрипции генов [23]. Рост объемов таких гетерогенных экспериментальных данных требует разработки новых биоинформационных решений [2, 19].
Авторы благодарны Y.Ruan и G.Li за научное обсуждение. Работа выполнена при поддержке бюджетного проекта ИЦиГ СО РАН VI.61.1.2, Интеграционного проект СО РАН №130, РФФИ 14-04-01906.
Список использованных источников
1. 3C-методы в исследованиях пространственной организации генома / Н.Р. Баттулин [и др.] // Вавиловский журнал генетики и селекции. - 2012. -Т. 16. № 4/2. - С. 872-876.
2. Practical Guidelines for the Comprehensive Analysis of ChIP-seq Data / T. Bailey [et al.] // PLoS Comput Biol. - 2013. - 9(11):e1003326.
3. Barba, M. Historical Perspective, Development and Applications of Next-Generation Sequencing in Plant Virology / M. Barba, H. Czosnek, A. Hadidi // Viruses. - 2014. -Vol. 6. - P. 106-136.
4. Integration of external signaling pathways with the core transcriptional network in embryonic stem cells / X. Chen [et al.] // Cell. - 2008. -Vol. 133(6). P. 1106-1117.
5. A genome-wide RNAi screen identifies PRDM14 as a regulator of POU5F1 and human
embryonic stem cell identity / N.-Y. Chia [et al.] // Nature. - 2010. - Vol. 468(7321). - P. 316-320.
6. dPeak: high resolution identification of transcription factor binding sites from PET and SET ChlP-Seq data / D. Chung [et al.] // PLoS Comput Biol. - 2013. - 9(10):e1003246.
7. Collas, P. Chop it, ChIP it, check it: the current status of chromatin immunoprecipitation / P. Collas, J.A. Dahl // Front Biosci. -2008. -Vol. 13. - P. 929-943.
8. Dekker, J. Exploring the three-dimensional organization of genomes: interpreting chromatin interaction data / J. Dekker, M.A. Marti-Re-nom, L A. Mirny // Nat Rev Genet. - 2013. -Vol. 14(6). P. 390-403.
9. An integrated encyclopedia of DNA elements in the human genome / ENCODE Project Consortium B.E. Bernstein [et al.] // Nature. -2012. - Vol. 489(7414). - P. 57-74.
10. An oestrogen-receptor-alpha-bound human chromatin interactome / M.J. Fullwood [et al.] // Nature. - 2009. - Vol. 462(7269). -P. 58-64.
11. The DNA sequence and biological annotation of human chromosome 1 / S.G. Gregory [et al.] // Nature. - 2006. - Vol. 441. - P. 315-321.
12. ChlP-chip versus ChIP-seq: lessons for experimental design and data analysis / J.W. Ho [et al.] // BMC Genomics. - 2011. - 28;12:134.
13. Iterative correction of Hi-C data reveals hallmarks of chromosome organization / M. Imakaev [et al.] // Nat Methods. - 2012. -Vol. 9(10). - P. 999-1003.
14. Integrative model of genomic factors for determining binding site selection by estrogen receptor a / R. Joseph [et al.] // Mol Syst Biol. -2010. - Vol. 6. - P. 456.
15. The DNA-encoded nucleosome organization of a eukaryotic genome / N. Kaplan [et al.] // Nature. - 2009. - Vol. 458(7236). - P. 362-366.
16. Kedes, L. The new date, new format, new goals and new sponsor of the Archon Genomics X PRIZE Competition / L. Kedes, G. Campany // Nature Genetics. - 2011. - Vol. 43. - P. 1055-1058.
17. Interactome maps of mouse gene regulatory domains reveal basic principles of transcriptional regulation / K.R. Kieffer-Kwon [et al.] // Cell. - 2013. - Vol. 155(7). - P. 1507-1520.
18. Computational analysis and modeling of genome-scale avidity distribution of transcription factor binding sites in chip-pet experiments
/ V.A. Kuznetsov [et al.] // Genome Inform. -
2007. - Vol. 19. - P. 83-94.
19. A practical comparison of methods for detecting transcription factor binding sites in ChlP-seq experiments / T.D. Laajala [et al.] // BMC Genomics. -2009. - Vol. 10. - P. 618.
20. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome / B. Langmead [et al.] // Genome Biol. - 2009. -Vol. 10(3):R25.
21. Lee, H. Genomic dark matter: the reliability of short read mapping illustrated by the genome mappability score / H. Lee, M.C. Schatz // Bio-informatics. - 2012. - Vol. 28(16). P. 2097-3105.
22. Differential motif enrichment analysis of paired ChIP-seq experiments / T. Lesluyes [et al.] // BMC Genomics. - 2014. - Vol. 15(1). -P. 752.
23. Extensive promoter-centered chromatin interactions provide a topological basis for transcription regulation / G. Li [et al.] // Cell. -2012. - Vol. 148(1-2). - P. 84-98.
24. SOAP: short oligonucleotide alignment program / R. Li [et al.] // Bioinformatics. -
2008. - Vol. 24. - P. 713-714.
25. Comprehensive mapping of long-range interactions reveals folding principles of the human genome / E. Lieberman-Aiden [et al.] // Science. - 2009. - Vol. 326(5950). - P. 289-293
26. DIP-chip: rapid and accurate determination of DNA-binding specificity / X. Liu [et al.] // Genome Res. -2005. - Vol. 15(3). - P. 421-427.
27. Mouse Encode Project at Ren Lab [Электронный ресурс] / Yue Lab at Penn State. - Mode of access: http://yuelab.org/hi-c/.
28. A tale of three next generation sequencing platforms: comparison of Ion torrent, pacific biosciences and illumina MiSeq sequencers / M. Quail [et al.] // BMC Genomics. - 2012. -Vol. 13. - P. 341-354.
29. Three dimensional genome database (3DGD) [Электронный ресурс] - Mode of access: http://3dgd.biosino.org/protein/page/view-Pattern.jsp.
30. Tucker, T. Massively parallel sequencing: the next big thing in genetic medicine / T. Tucker, M. Marra, J.M. Friedman // Am J Hum Genet. - 2009. - Vol. 85(2). - P. 142-154.
31. A global map of p53 transcription-factor binding sites in the human genome / C.L. Wei [et al.] // Cell. - 2006. - P. 124(1). - P. 207-219.
32. Global mapping of c-Myc binding sites and target gene networks in human B cells / K.I. Zeller [et al.] // Proc Natl Acad Sci USA. -2006. - Vol. 103. - P. 17 834-17 839.
33. Model-based analysis of ChIP-Seq (MACS) / Y. Zhang [et al.] // Genome Biol. -2008. - Vol. 9(9):R137.
Дата поступления статьи 22 сентбря 2014 г.
