Поиск пиков в данных ChiP-Seq

Поиск пиков (англ. Peak calling) — компьютерный метод поиска областей генома, обогащённых выровненными ридами из данных ChiP-Seq эксперимента. Такие области называют пиками, поскольку они являются локальными максимумами покрытия генома ридами (то есть количества прочтений, выровненных в данной позиции). Определённые таким образом регионы считаются местами связывания белка, с которым проводился опыт. Если речь идёт о транскрипционном факторе, то обогащённые области по сути являются сайтами связывания транскрипционного фактора. Определение областей пиков является важным шагом любого Chip-Seq эксперимента, поскольку именно на этом этапе происходит интерпретация данных — нужно отличить шумы от сигнала^[1].

Операция нахождения пиков содержит в себе несколько этапов, необходимых для получения итогового списка пиков^[2]:

1. Смещение ридов (важно в случае секвенирования с одного конца — single-end sequencing). Каждый рид, полученный после секвенирования, имеют длину 50-300 нуклеотидов и маркирует только один из концов фрагмента ДНК, очищенного иммунопреципитацией. Соответственно, он может не содержать в себе сайта связывания белка, если длина ДНК фрагментов больше длины рида, а белок был связан с центром фрагмента. На начальном этапе поиска идентифицируют области, обогащённые ридами. После этого находят отдельно максимумы обогащения для «+» и «-» цепей. Расстояние между максимумами считают средней длиной фрагментов ДНК. После этого все риды смещают на половину такой длины в сторону центра пика (то есть в сторону их 3’-конца), и определяют истинную область связывания. В случае секвенирования парных ридов (с обоих концов каждого фрагмента), эта операция не требуется, поскольку пара ридов точно маркирует каждый фрагмент, а значит сайт связывания заведомо находится между ними.
2. Оценка уровня шума. Чаще всего для оценки уровня шума проводят контрольные Chip-Seq эксперименты, в которых стадия иммунопреципитации пропущена или используются неспецифические антитела (IgG-контроли). Однако существуют программы, способные оценить шум по данным одного эксперимента.
3. Идентификация пиков. Чаще всего используется скользящее окно, которое сканирует геном и определяет регионы, в которых покрытие ридами превышает некий заданный порог, или отношение покрытия условного пика превышает окружающий его шум в заданное число раз.
4. Оценка качества определённых пиков. Многие программы подсчитывают P-value для каждого пика, его высоту и обогащения ридами относительно фона. Исходя из статистической значимости пиков можно сделать оценку ожидаемого уровня ложных предсказаний (FDR).
5. Удаление артефактов. В данных Chip-Seq встречаются артефакты ПЦР, возникающие при пробоподготовке библиотек для секвенирования. При фильтрации пиков на предмет артефактов, удаляются такие из них, которые содержат только дуплицированные (то есть идентичные) риды, поскольку считается, что такие пики возникли из-за ПЦР. Также критически относятся к ситуациям, когда риды картированы резко асимметрично (только на 3’ или 5’ цепь).

На качество определения пиков влияет множество параметров, начиная от способа картирования ридов на геном. Так, количество разрешённых при картировании несовпадений влияет на долю ридов, которые будут успешно картированы, а значит и на профиль покрытия генома. А операции по удалению ридов, способных выравниваться с несколькими местами в геноме, могут привести к потере пиков, расположенных в областях с повторами. Кроме того, в каждом из пяти перечисленных пунктов обобщённого алгоритма идентификации пиков существуют свои параметры, аккуратный подбор которых позволяет добиться максимальной селективности и специфичности^[3][4].

К настоящему времени существуют десятки программ для определения пиков, отличающиеся по своей применимости и используемым алгоритмам. В данной статье будет разобрана самая популярная программа — MACS, которая хорошо подходит для поиска относительно узких пиков, характерных для сайтов связывания транскрипционных факторов. Также будет рассмотрена программа SICER, целью которой является поиск слабых и диффузных сигналов (то есть широких и слабых пиков), характерных для распределения модификаций гистонов или отражающих работающие полимеразы^[5].

MACS (англ. Model-based Analysis of Chip-Seq) — программа анализа данных ChiP-Seq эксперимента. Используемый программой метод хорошо подходит для поиска относительно узких пиков, характерных для сайтов связывания транскрипционных факторов. MACS является наиболее часто используемым алгоритмом для поиска пиков.

1. Моделирование размера сдвига ChiP-Seq ридов

ChiP-Seq риды представляют собой концы фрагментов в ChiP-DNA библиотеке. По описанным выше причинам (см. Обобщённый алгоритм), они часто сдвинуты в сторону 5'-концов и размер этого сдвига неизвестен экспериментатору. Так как фрагменты ChiP-DNA равновероятно секвенированы с обоих концов, то плотность ридов вокруг истинного сайта связывания будет бимодальным: на «+» цепи плотность возрастает до сайта связывания, на «-» цепи — после (рис.1). MACS использует этот двойной паттерн для построения эмпирической модели размера сдвига ридов в сторону их 3’-конца для лучшего определения сайта связывания. Для первичного поиска регионов, богатых ридами, данный алгоритм использует 2 параметра: bandwidth (средняя длина фрагментов ДНК, получаемая при пробоподготовке, например, при облучении ультразвуком) и mfold (пороговое обогащение относительно уровня фона, используемого моделью, выше которого заданная область считается статистически обогащённой ридами). Скользя по геному «окном» размером 2*bandwidth алгоритм ищет регионы, mfold которых выше, чем mfold при случайном распределении ридов по геному. Затем MACS выбирает 1000 случайных, из определённых на предыдущем шаге областей, и отделяет в них риды «+» и «-» цепей. Если максимум плотности ридов «+» цепи лежит левее максимума для «-» цепи, то алгоритм выравнивает их относительно центров распределений. Расстояние между максимумами «+» и «-» цепей в выравнивании будет определяться как d (рис.1). MACS сдвигает все теги на расстояние d/2 в сторону 3'-конца для определения наиболее вероятного сайта взаимодействия ДНК с белком.

2. Поиск пиков

В экспериментах с контролем (см. Обобщённый алгоритм), MACS линейно увеличивает число контрольных ридов до общего количества ридов, полученных в ChiP-Seq эксперименте. В случае дуплицированных ридов, MACS способен эффективно удалять идентичные последовательности для снижения уровня шума.

Для статистической оценки значимости пиков программа использует предположение, что количество ридов, выровненных в данной позиции генома, распределено по Пуассону^[6]. Преимуществом данной модели является то, что один параметр — λ_BG, оценённый для всего генома по контрольным точкам, определяет среднее значение и дисперсию распределения. После того, как MACS сдвинул каждый относящийся к обогащённому региону рид на d/2, алгоритм ещё раз сканирует геном, теперь уже «окном» шириной 2d и находит статистически значимые пики (используя приближение о Пуассоновском распределении числа ридов с параметром λ_BG и уровнем значимости по умолчанию 10⁻⁵). Пересекающиеся пики сливаются, а участок их слияния называют summit. Он определяет локализацию сайта связывания.

В контрольных образцах часто наблюдаются флуктуации в распределении ридов. На это могут влиять многие факторы, такие как локальные структуры хроматина, шум, возникающий при амплификации и секвенировании и вариации в числе копий генома. Поэтому, MACS может использовать не только λ_BG, оценённый по всему геному, но и динамический параметр λ_local, который вычисляют для каждого пика как:

${\displaystyle \textstyle \lambda _{local}=\max {\bigl (}\lambda _{BG},\lambda _{1k},\lambda _{5k},\lambda _{10k}{\bigr )}}$

где λ_1k, λ_5k и λ_10k — это параметры, определённые для «окон» шириной в 1 кб, 5 кб и 10 кб соответственно, центрированных относительно рассматриваемого пика . Таким образом, λ_local отражает влияние местных флуктуаций.

Далее MACS использует λ_local для подсчёта p-value каждого потенциального пика и удаляет ложно положительные результаты, вызванные шумом в покрытии ридами. Потенциальные пики с p-value ниже заданного (по умолчанию 10⁻⁵) отмечаются, а отношение между количеством ридов ChiP-Seq, выровненных в данной позиции и λ_local называют fold enrichment. Для Chip-Seq эксперимента с контролем, MACS вычисляет ожидаемый уровень ложных предсказаний (FDR) для каждого определённого пика, используя процедуру, которая применялась и в ранних алгоритмах поиска, таких как MAT^[7] и MA2C^[8].

MACS сравнивали с другими алгоритмами по поиску пиков, например — ChIPSeq Peak Finder^[9], FindPeaks и QuEST. Для того чтобы сравнить специфичность их предсказаний, ChiP-Seq данные были поменяны местами с соответствующими им контролями, а по результатам программной обработки вычислили FDR каждого алгоритма. MACS каждый раз давал меньше ложноположительных результатов, нежели другие методы (рис.2).

SICER (англ. spatial identification of ChIP-enriched regions) — это компьютерный метод, разработанный для идентификации кластеров, обогащённых выровненными ридами на основании данных ChiP-Seq эксперимента с учётом оценки уровня шума. Хотя данный метод создан для анализа слабых и диффузных сигналов, характерных для модификаций гистонов, при правильном выборе параметров он может быть применён для поиска узких пиков, соответствующих сайтам связывания, например, транскрипционных факторов.

Программа представляет собой пакет, использующий в качестве интерфейса командную строку, подходит для операционных систем Unix/Linux, язык программирования: Python.

Главной особенностью метода является объединение сигналов от соседних нуклеосом, несущих одинаковую модификацию и расположенных близко друг к другу. Данный подход улучшает отношение сигнал-шум и особенно полезен при диффузных модификациях хроматина, которые могут иметь большую (килобазы или даже мегабазы) протяжённость в геноме. C помощью математического подхода идентифицируются статистически значимые области генома, обогащённые выровненными ридами. Далее проводится оценка уровня шума с использованием контрольной библиотеки (создаётся на основе контрольных ChIP-Seq экспериментов на тотальной клеточной ДНК без стадии иммунопреципитации)^[10].

На основании полученных данных идентифицируется набор локальных максимумов покрытия ридами (пиков), для которых отношения покрытия ридами превышает фоновый шум в определённое пороговое число раз. Геном делится на неперекрывающиеся участки фиксированной длины — «окна», статистически значимые области генома называются «острова». Статистически значимыми областями генома принято считать кластеры окон, для которых отношение их покрытия ридами к фоновому покрытию превышает некий заданный порог. Окна внутри островов, для которых велик фоновый шум или маленькое покрытие ридами — «гэпы» (рис. 3).

Острова представляют особый интерес. Гэпы разрешены в островах для учёта:

• Неточности ридов и неоднозначности данных секвенирования
• Повторяющихся геномных участков и дуплицированных ридов
• Немодифицированных нуклеосом

Размер гэпа может быть скорректирован в зависимости от природы модификации хроматина. Статистический вес острова высчитывается, основываясь на степени обогащения ридами на всем идентифицированном участке, а не на степени обогащения ридами каждого отдельного пика^[11].

На вход программе подаётся библиотека выровненных ридов из данных ChIP-Seq эксперимента, выбирается размер окна, эффективная длина генома (наибольшая длина генома, покрытая выровненными ридами) и размер гэпа. Размер гэпа должен быть кратным выбранному размеру окна. Как правильно, чем большую протяжённость имеют домены, тем большого размера следует выбирать гэп. Для точной оценки размера гэпа можно представить суммарный статистический вес всех значимых островов как функцию от размера гэпа. Таким образом, лучшим размером гэпа будет тот, который соответствует наибольшему весу островов. На рис. 4 рассмотрен пример такого графика, в данном случае размер гэпа следуют выбирать равным шести окнам.

При наличии контрольной библиотеки идентифицированный набор островов проверяется на способность кластеризоваться с мягким порогом Е-value (основан на предварительном анализе обогащения ридами относительно фона в контрольной библиотеке). Далее SICER вычисляет p-value и q-value для каждого идентифицированного острова с использованием распределения Пуассона и FDR уже с учётом удаления артефактов. При отсутствии контрольной библиотеки тоже можно получить статистически значимые острова, используя различные пороги E-value. На выходе программы получается набор островов, высоко обогащённых выровненными ридами^[12] (рис. 5).

Особенность данной программы состоит в том, что для её алгоритма локальная степень обогащения сильно зависима от контекста. Иными словами, если окно является обогащённым, этого ещё недостаточно, чтобы оно стало статистически значимым. Однако в присутствии соседних обогащённых окон оно (наряду с другими членами кластера) становится таковым. Это и отличает SICER от подходов, основанных на локальной статистике, которые бы идентифицировали бы его как не значимый в обоих случаях. Данный метод, благодаря нивелированию общего шума, подходит для таких задач, как статистический анализ данных ChIP-Seq экспериментов, нормализация данных и масштабный анализ больших последовательностей с диффузным обогащением^[10].

Математика

1. Схема подсчёта статистических весов

Геном эффективной длины ${\displaystyle L}$ делится на неперекрывающиеся окна размером ${\displaystyle w}$. Статистический вес ${\displaystyle s}$ для окна с ${\displaystyle l}$ ридами идентифицируется как ${\displaystyle s(l)=-\log {P(l,\lambda )}}$, где ${\displaystyle P(l,\lambda )}$ — распределение Пуассона с параметром ${\displaystyle \lambda =wN/L}$, учитывающим среднее число ридов в окне; ${\displaystyle N}$ - общее число ридов в библиотеке ChIP-Seq. Статистические веса окон представляют собой отрицательный логарифм вероятности нахождения ${\displaystyle l}$ ридов в окне при условии, что риды могут охватывать любой участок в геноме с равной вероятностью (модель случайного фона). Статистические веса кластеров окон являются аддитивными, представляющими собой отрицательный логарифм совместной вероятности нахождения наблюдаемой конфигурации в модели случайного фона. Чем выше статистический вес, тем меньше вероятность того, что наблюдаемый профиль является случайным.

2. Идентификация острова

Для числа ридов в каждом окне устанавливается порог ${\displaystyle l_{0}}$_, полученный из значения P-value, основанного на распределении Пуассона:

${\displaystyle \sum _{l=l_{0}}^{\infty }{P(l,\lambda )\leq {p_{0}}}}$

Более того, ${\displaystyle l_{0}}$ зависит от размера ChIP-Seq библиотеки.

3. Рекуррентное соотношение для вероятности острова с данным статистическим весом на случайном фоне

Сначала высчитывается вероятность распределения статистических весов отдельного окна:

${\displaystyle p(s)=\sum _{l\geq {l_{0}}}{\delta (s-s(l))P(b,\lambda )}}$

где ${\displaystyle \delta ()}$ — дельта-функция Дирака. Вероятность ${\displaystyle t}$ того, что окно является «неподходящим» равна:

${\displaystyle t=P(0,\lambda )+P(1,\lambda )+...+P(l_{0}-1,\lambda )}$

Далее вводится параметр ${\displaystyle g}$, отражающий число гэпов в окне. Острова идентифицируются как кластеры окон, разделённые гэпами, не большего размера, чем ${\displaystyle g}$. На основании этого высчитывается гэп-фактор ${\displaystyle G}$. Статистический вес острова является совокупным статистическим весом всех окон на данном острове.

${\displaystyle G=1+t+t^{2}+...+t^{g}}$

${\displaystyle M(s)}$ — вероятность нахождения острова со статистическим весом ${\displaystyle s}$, начиная от заданной позиции вдоль генома; зависит от ${\displaystyle \lambda }$, ${\displaystyle l_{0}}$ и ${\displaystyle g}$ через ${\displaystyle p(s)}$, ${\displaystyle t}$ и ${\displaystyle G}$. Так как остров должен быть окружён гэпами размеров не меньше, чем ${\displaystyle g+1}$, ${\displaystyle M(s)}$ можно записать как сумму:

${\displaystyle M(s)=t^{g+1}\mathrm {M} (s)t^{g+1}}$

Тогда рекуррентное соотношение:

${\displaystyle \mathrm {M} (s)=G(\lambda ,l_{0},g)\int \limits _{s_{0}}^{s}ds'\mathrm {M} (s-s')p(s')}$

Граничные условия: ${\displaystyle \mathrm {M} (0)=G(\lambda ,l_{0},g)^{-1}}$ и ${\displaystyle s_{0}=-ln{P(l_{0},\lambda )}}$.

Число островов со статистическим весом: ${\displaystyle s-LM(s)}$.

4. Асимптотика для распределения статистических весов островов на случайном фоне

Подставим анзац в рекуррентное соотношение:

${\displaystyle \mathrm {M} (s)=\alpha {exp(-\beta {s})}}$

Тогда

${\displaystyle G(\lambda ,l_{0},g)\sum _{l\geq {l_{0}}}P(l,\lambda )^{1-\beta }=1}$

5. Нахождение значимых островов без использования контрольной библиотеки

Распределение статистических весов островов на случайном фоне позволяет выявить пороговое значение статистического веса ${\displaystyle s_{T}}$_,. Данный параметр используется для поиска статистически значимых островов. ${\displaystyle s_{T}}$ определяется c учётом того, что ожидаемое число островов со статистическим весом выше, чем порог ${\displaystyle s_{T}}$_, должно быть не больше порогового значения E-value ${\displaystyle e}$:

${\displaystyle \sum _{s\geq {s_{T}}}LM(s)\leq {e}}$

6.Нахождение значимых островов c использованием контрольной библиотеки

Сначала используется мягкий порог E-value для идентификации набора островов-кандидатов, которые демонстрируют кластеризацию при модели случайного фона. Затем для каждого кандидата высчитывается число картированных ридов ${\displaystyle n_{s}}$ и ридов контрольной библиотеки ${\displaystyle n_{c}}$. P-value вычисляется как:

${\displaystyle \textstyle \sum _{n=n_{s}}^{\infty }\displaystyle P(n_{s},cn_{c})}$

где с равно соотношению размера библиотеки ChIP к размеру контрольной библиотеки (${\displaystyle c=N_{s}/N_{c}}$). Острова-кандидаты с ${\displaystyle n_{s}\leq {cn_{c}}}$ отбрасываются, потому что необходимо рассматривать только случаи обогащения. Значимые острова могут быть идентифицированы с помощью порога P-value с использованием поправки Бонферрони^[10].