RNAcontacts, вычислительный конвейер для предсказания РНК-РНК-контактов по данным конформационного секвенирования РНК in situ

Полный текст

СПИСОК СОКРАЩЕНИЙ КС – конформационное секвенирование; RNA-seq – секвенирование РНК (RNA sequencing); ЛПБМ – лигирование пространственно близких молекул.

ВВЕДЕНИЕ

Быстрое развитие технологий высокопроизводительного секвенирования позволило выявлять пространственные контакты между нуклеиновыми кислотами, включая контакты ДНК в трехмерной структуре хроматина [1–3], функциональные взаимодействия энхансеров с промоторами [4, 5], а также хроматин-ассоциированные контакты РНК с ДНК [6, 7]. Эти методы основаны на принципе лигирования пространственно сближенных молекул (ЛПБМ, proximity ligation assays) – разрезании и последующем стохастическом лигировании нуклеиновых кислот, входящих в состав макромолекулярных комплексов, которое преимущественно происходит между пространственно близкими участками. Глубокое секвенирование полученных химерных фрагментов дает сотни миллионов чтений, содержащих последовательности взаимодействующих локусов.

В некоторых недавно разработанных методах для изучения взаимодействий РНК-РНК in vivo и in vitro применяется ЛПБМ [8, 9]. В некоторых из этих методов, таких, как PARIS [10], LIGR-seq [11], SPLASH [12] и COMRADES [13], используются производные псоралена для того, чтобы вызвать обратимое сшивание между дуплексами РНК для определения комплементарных взаимодействий. В протоколе конформационного секвенирования РНК in situ (RIC-seq) цепи РНК сшиваются через РНК-связывающие белки (РСБ) [14], что позволяет не только устанавливать вторичные и третичные структуры РНК, но и помогает воссоздать трехмерные карты взаимодействий РНК с РСБ. Во всех этих случаях взаимодействия закодированы в химерных последовательностях РНК, полученных путем разрезания и последующего лигирования.

В отличие от взаимодействий ДНК с ДНК, проявляющихся в данных КС в виде чтений с разрывами, которые картируются только на два пространственно близких геномных локуса, взаимодействия РНК с РНК дают чтения, которые могут картироваться более сложным образом, поскольку пре-мРНК подвергаются сплайсингу. В частности, химерные фрагменты могут содержать как экзон-экзонные соединения (ЭЭС), так и продукты сшивки пространственно близких РНК, в результате чего образуются как чтения с каноническими интронными GT/AG-разрывами, возникающими при сплайсинге, так и чтения с другими разрывами, возникающими при лигировании (рис. 1A). Точное картирование таких чтений представляет сложную задачу, поскольку большинство картировщиков может работать только с одним типом разрывов. Так, картировщик, имеющий только одну модель разрывов, был бы вынужден либо увеличивать штраф для разрывов без GT/AG-мотива для более точной идентификации координат канонических интронов, либо ослаблять требования к наличию GT/AG-мотива с целью более точного нахождения координат РНК-контактов. Поэтому разработка вычислительного метода, позволяющего картировать короткие чтения с двумя различными типами разрывов, является актуальной задачей биоинформатики. В данной работе представлен вычислительный конвейер, который позволяет достичь этой цели без разработки специализированного картировщика.

ЭКСПЕРИМЕНТАЛЬНАЯ ЧАСТЬ

Геномы и их аннотации

Сборка генома человека hg19 (февраль 2009) и аннотация транскриптома GENCODE v34lift37 были загружены с веб-сайта Genome Reference Consortium [15] и веб-сайта GENCODE [16] соответственно. Координаты интронов получены из выходных данных программы STAR (см. ниже).

Данные высокопроизводительного секвенирования

Две реплики данных RIC-seq с деплецией рРНК (GSM3629915 и GSM3629916) в клеточной линии HeLa [14] были загружены из Gene Expression Omnibus в формате FASTQ (GSE127188). Соответствующий набор контрольных данных секвенирования РНК (RNA-seq) в клеточной линии HeLa был загружен из консорциума ENCODE (ENCLB555ASI и ENCLB555ASJ). На первом проходе данные RNA-seq картировали на геном человека с помощью программы-картировщика STAR версии 2.7.3a в парноконцевом режиме со следующими дополнительными параметрами:

--runMode alignReads --outSAMtype BAM SortedByCoordinate --chimOutType Junctions.

На втором проходе данные RIC-seq картировали на геном человека с помощью той же версии STAR со следующими дополнительными параметрами:

--chimSegmentMin 15 --chimJunctionOverhangMin 15 --chimScoreJunctionNonGTAG -1 --scoreGapNoncan -1 --scoreGapATAC -1 --scoreGapGCAG -1 --chimSegmentReadGapMax 3 --outFilterMatchNminOverLread 0.5 --outFilterScoreMinOverLread 0.5

Параметр --chimSegmentReadGapMax 3 введен для пропуска при картировании дополнительного биотинилированного остатка цитозина в RIC-seq [14]. При втором проходе штраф снижается до -1 для всех типов неканонических разрывов.

Реализация метода

RIC-contacts реализован в популярной системе управления задачами Snakemake [17] и находится в свободном доступе на GitHub [18]. Пути к входным файлам указываются в конфигурационном файле в формате yaml, который также содержит настройки STAR и дополнительные параметры, контролирующие минимальное расстояние между двумя точками лигирования в кластере и порог разброса, для визуализации коллинеарных разрывов через UCSC Genome Browser [19]. Коллинеарные чтения извлечены из BAM-файлов с помощью отдельных программ (neo.pl в репозитории RNAcontacts) и пакета samtools v1.14 [20]. Пакет bedops v2.4.41 использовался для кластеризации точек разрыва [21]. Количество поддерживающих чтений подсчитано с помощью пакета bedtools v2.29.0 [22].

Визуализация результатов

Для визуализации карт контактов списки контактов были конвертированы в формат ‘cool‘ с помощью пакета ‘cooler‘ v0.8.11 с разрешением 100 п.н. Полученные карты визуализировали с помощью пакета pygenometracks v3.7. Визуализация чтений с разрывами выполнена с помощью IGV v2.11.2 и UCSC Genome Browser [19]. По умолчанию в UCSC Genome Browser визуализируются только коллинеарные контакты с разбросом не более 50000 нуклеотидов (см. также руководство [18]).

РЕЗУЛЬТАТЫ

Метод RNAcontacts

Для анализа данных КС РНК нами разработан вычислительный конвейер RNAcontacts, который обходит проблему нескольких различных типов разрывов путем выравнивания коротких чтений в двухпроходном режиме (рис. 1Б). Метод основан на использовании картировщика STAR [23], который на первом этапе выравнивает набор данных секвенирования РНК (RNA-seq) в парноконцевом режиме, чтобы определить интроны, которые экспрессируются в данном биологическом образце, используя строгий штраф для не-GT/AG-разрывов. Здесь RNA-seq представляет собой контрольный эксперимент, не содержащий химерных фрагментов, полученных в результате ЛПБМ. На втором проходе чтения сгенерированные в эксперименте КС РНК выравниваются с использованием ослабленного штрафа за не-GT/AG-разрывы, при этом на вход подается список разрывов, идентифицированных на первом проходе как интроны, так что картировщик будет преимущественно делать разрывы в чтениях по координатам из предоставленного списка. Поскольку данные КС РНК могут содержать химерные чтения на произвольном геномном расстоянии или in trans, выравнивание на втором проходе выполняется в одноконцевом режиме. Все выравнивания с разрывами, полученные на втором проходе, анализируются для извлечения контактов РНК–РНК и исключения интронов, полученных на первом проходе.

 

Рис. 1. A – в протоколе RIC-seq [14] несплайсированный участок РНК может быть лигирован через РНК-связывающий белок (РСБ) с другим участком, содержащим экзон-экзонные соединения (ЭЭС). Последовательность, образованная в результате лигирования, выравнивается на геном с двумя разрывами: не-GT/AG-разрывом, соответствующим точке лигирования, и каноническим GT/AG-разрывом, соответствующим сплайсингу. Б – схема RNAcontacts. На первом проходе чтения из контрольного эксперимента RNA-seq выравниваются на референсный геном для выявления экспрессируемых ЭЭС. Эти ЭЭС используются на втором проходе как известные интроны при выравнивании данных КС для обнаружения разрывов, которые кодируют контакты РНК–РНК

 

Программы выравнивания чтений с разрывами обычно генерируют два отдельных выходных файла, соответствующих коллинеарным и неколлинеарным разрывам. В частности, программа-картировщик STAR сообщает о коллинеарных разрывах (одна и та же хромосома, одна и та же цепь, согласованный порядок картирования фрагментов) в стандартном выходном файле SAM/BAM, в то время как неколлинеарные разрывы помещаются в отдельный выходной файл с химерными разрывами, так как BAM-формат не позволяет представить их с помощью одной CIGAR-строки [23]. RNAcontacts извлекает координаты неоразрывов, т.е. коллинеарных разрывов, найденных на втором проходе, из результатов SAM/BAM и объединяет их с химерными разрывами. Следует отметить, что не только trans-, но и cis-контакты могут быть закодированы как в нео-, так и в химерных разрывах. Объединенный результат второго прохода состоит из нео- и химерных разрывов в чтениях, которые соответствуют точкам лигирования РНК.

В ходе эксперимента RIC-seq в клеточной линии HeLa [14] были картированы 94.3% из ~224 млн чтений из двух реплик, где 72.0% картированных чтений картированы уникально (табл. S1). При этом 18.5% уникально картированных чтений содержали по крайней мере один разрыв, в то время как значения для LIGR-seq, PARIS и SPLASH, как сообщалось ранее, составили 3.5, 2 и 0.5% соответственно [24]. Отметим, что результаты работы программ картирования коротких чтений, учитывающих сплайсинг, могут различаться. При использовании протокола RIC-seq небольшие различия в координатах разрывов могут возникать даже при сопоставлении чтений, которые перекрывают одну и ту же точку лигирования, поскольку картирование одной и той же последовательности на одну или другую цепь может приводить к смещению координат из-за отсутствия консенсусных последовательностей у точек лигирования (рис. S1). Кроме того, различные копии одной и той же РНК разрезаются и повторно лигируются стохастически, что приводит к еще большей вариабельности. Учитывая эту техническую и биологическую вариативность, мы ожидаем выявление кластеров точек лигирования, а не четко определенных сайтов, как это имеет место в случае сайтов сплайсинга в GT/AG-интронах.

Распределение расстояний между двумя последовательными точками разрыва быстро убывает, причем приблизительно 50% расстояний составляет меньше 9 нуклеотидов и 90% расстояний – меньше 21 нуклеотида (рис. S2). Поэтому точки разрыва были кластеризованы с использованием односвязной (single-linkage) кластеризации с порогом по расстоянию (δ) 10 и 20 нуклеотидов (рис. S3). Контакты определены как пары кластеров, соединяемые хотя бы одним разрывом, а число чтений, поддерживающих контакт, определяли как суммарное число чтений, поддерживающих все разрывы из кластера.

Для каждого значения δ мы разделили контакты на три группы: внутригенные (оба конца контакта принадлежат аннотированному гену), контакты in cis (на одной хромосоме, но не в одном гене) и контакты in trans (на разных хромосомах). Число контактов (n), длина кластера (s), расстояние между контактирующими кластерами (d, которое определяется только для внутригенных и cis-контактов) и число поддерживающих чтений (r) незначительно различались для двух значений δ (табл. 1). В среднем мы обнаружили на 30% больше внутригенных контактов, чем контактов in cis, и более чем двукратное обогащение контактов in trans по отношению к двум другим группам. При δ = 10 большинство кластеров имели длину 10 нуклеотидов (рис. S4), что указывает на то, что они состоят только из одной индивидуальной точки лигирования, окруженной с обеих сторон участками из 5 нуклеотидов.

 

Таблица 1. Характеристика кластеров контактов РНК-РНК

δ	Метрика	Внутригенные кластерыконтактов	in cis	in trans
10	n	1369158	1061470	4881920
	s	10.1±2.2	10.1±2.1	10.2±2.2
	log2d	10.8±3.3	17.6±4.4	нет
	log2r	0.8±0.8	0.5±0.7	0.5±0.7
20	n	1313727	1035656	4851697
	s	20.2±2.4	20.2±2.3	20.3±2.4
	log2d	11.0±3.3	17.6±4.4	нет
	log2r	0.8±0.8	0.5±0.7	0.5±0.7

Примечание. Расстояние между кластерами, δ. Количество контактов, n. Длина кластера, s. Расстояние между контактирующими кластерами, d. Количество чтений, поддерживающих контакт, r. Показано среднее значение ± стандартное отклонение.

 

Расстояния между контактирующими кластерами по-разному распределены в нео- и химерных разрывах в чтениях как по количеству собственно контактов, так и при взвешивании по количеству поддерживающих чтений (рис. S5). Примечательно, что распределения имеют две моды, причем первая мода d-~1000 соответствует внутригенным контактам, задаваемым как нео-, так и химерными разрывами. Химерные разрывы могут задавать внутригенные контакты, если лигирование произошло в обратной ориентации, как например, в кольцевых РНК [25]. Вторая мода для неоразрывов обусловлена условием d≤250000, которое накладывает картировщик STAR на коллинеарные разрывы, однако более длинные контакты in cis захвачены химерным распределением. При этом большинство контактов in cis и in trans поддерживается только одним чтением, в то время как большинство внутригенных контактов поддерживается двумя чтениями (рис. S6). Таким образом, поддержка чтениями в отдельных экспериментах RIC-seq довольно скудна даже после объединения контактов в кластеры.

Чувствительность RNAcontacts

Чтобы сравнить производительность RNAcontacts с производительностью RICpipe (метод, изначально разработанный для анализа данных RIC-seq), мы проанализировали внутрихромосомные разрывы в чтениях с разбросом 50 нуклеотидов или более. Мы исключили разрывы в чтениях, картируемых на рРНК, из результатов RNAcontacts, поскольку RICpipe также удаляет чтения рРНК [14]. Только 40% (соответственно, 45%) разрывов, идентифицированных с помощью RNAcontacts (соответственно, RICpipe), имели точно такие же координаты, как и разрывы, идентифицированные другим методом, что указывает на различия в процедурах картирования (рис. 2A). Однако RNAcontacts удалось выровнять больше чтений, поддерживающих идентифицированные разрывы, по сравнению с RICpipe, что указывает на увеличение чувствительности примерно на 40% (рис. 2Б). При проведении сравнения с использованием окон размером 100 нуклеотидов, т.е. без точного совпадения координат, мы заметили, что результаты двух методов в значительной степени согласованы, о чем также свидетельствуют визуально похожие карты контактов с немного большим количеством контактов, полученных RNAcontacts, по сравнению с RICpipe (рис. 2В).

 

Рис. 2. Сравнение RNAcontacts и RICpipe. A – диаграмма Венна для координат разрывов, полученных по RNAcontacts и RICpipe. Б – то же, но координаты разрывов взвешены по их поддержке чтениями. В – карты контактов в гене PLEC-210, полученные из RNAcontacts (вверху) и RICpipe (внизу)

 

Затем мы проверили производительность RNAcontacts на данных RIC-seq в клеточной линии HeLa с первым проходом картирования и без него. Для этого мы отдельно выполнили второй проход RNAcontacts, предоставляя только аннотированные интроны из GENCODE [16], не добавляя специфичные для клеточной линии HeLa интроны, найденные на первом проходе. В результате получен примерно 1% ложных соединений лигирования, которые соответствуют неаннотированным интронам в HeLa. Мы также обнаружили, что 16809 из ~3.5 млн разрывов в чтениях, идентифицированных с помощью RICpipe, соответствуют экзон-экзонным соединениям. Хотя количество таких разрывов невелико, они поддерживаются значительной долей (>30%) чтений. Таким образом, двухпроходный метод обеспечивает более высокую специфичность (меньший процент ложноположительных результатов) при обнаружении контактов РНК, особенно в условиях, когда экспрессируемый транскриптом значительно отличается от аннотированного.

ОБСУЖДЕНИЕ И ВЫВОДЫ

В данной работе представлено концептуальное решение проблемы картирования чтений с двумя типами разрывов, характерными для экспериментов КС РНК. Несмотря на то, что представленное решение задействует картировщик STAR, предложенный подход не ограничивается использованием только этой программы, и вместо нее может быть использована любая другая программа для выравниваний чтений с разрывами [26]. Мы продемонстрировали, что экзон-экзонные соединения составляют большую часть чтений с разрывами в данных RIC-seq, и RNAcontacts позволяет обнаруживать разрывы, соответствующие точкам лигирования в экспериментах RIC-seq, с большей чувствительностью, чем RICpipe. Реализация RNAcontacts в воспроизводимой и масштабируемой системе управления рабочим процессом Snakemake позволяет быстро и единообразно обрабатывать множество наборов однотипных данных КС.

Данные экспериментов конформационного секвенирования РНК похожи по характеру на данные Hi-C, но имеют важные отличия, связанные с разрешением. Если для Hi-C обычной практикой является усреднение контактов хроматина на уровне тысяч или миллионов пар нуклеотидов, то обнаружение контактов РНК с помощью КС по своей сути нацелено на масштаб отдельных нуклеотидов. При этом поддержка известных контактов чтениями RIC-seq, например, в известной вторичной структуре в гене SF1, очень слаба (рис. 3). Большинство контактов RIC-seq in cis и in trans поддерживаются только одним чтением, в связи с чем возникает проблема оценки статистической значимости наблюдаемых контактов. Эта проблема, возможно, будет решена в исследованиях, анализирующих большие наборы экспериментов RIC-seq. Мы ожидаем, что в скором времени появится гораздо больше наборов данных экспериментов КС РНК, аналогичных RIC-seq, которые могут быть проанализированы представленным методом.

 

Рис. 3. Разрывные чтения, поддерживающие структуру РНК в гене SF1 человека [28]. Комплементарные цепи показаны оранжевым цветом. Ниже показаны разрывные чтения (чтения из двух реплик показаны синим и зеленым цветом)

 

Таким образом, RNAcontacts реализует вычислительный конвейер анализа контактов РНК–РНК, учитывающий два типа разрывов, специфичных для метода КС РНК. Несмотря на то, что программное обеспечение разработано для протокола RIC-seq, область его применения распространяется на любой аналогичный метод КС, в котором одним из взаимодействующих партнеров является РНК.

Авторы выражают благодарность Т.М. Иванову за дополнительное тестирование программы.

Работа поддержана грантами Министерства науки и образования Российской Федерации (075-10-2021-116) и исследовательского гранта Национальной ключевой программы исследований и разработок Китая (2021YFE0114900).

Авторы заявляют, что у них нет конкурирующих интересов.

Д.Д.П. разработал план и руководил исследованием; Г.Л. и Ч.К. предоставили данные и участвовали в их обработке; С.Д.М. провел анализ данных и разработал программное обеспечение. Все авторы участвовали в написании текста статьи, прочитали и одобрили его.

ДОПОЛНИТЕЛЬНЫЕ МАТЕРИАЛЫ

Дополнительные данные и материалы доступны в репозитории Zenodo (https://zenodo.org/record/7027475) [27].

Об авторах

Сергей Д. Маргасюк

Сколковский институт науки и технологий

Email: d.pervouchine@skoltech.ru
Россия, Москва, 121205

Мария А. Власенок

Сколковский институт науки и технологий

Email: d.pervouchine@skoltech.ru
Россия, Москва, 121205

Гуо Ли

Колледж наук о жизни, Чжэцзянский университет

Email: d.pervouchine@skoltech.ru
Китайская республика, Ханчжоу, Чжэцзян, ZJ310058

Чанчан Као

Ключевая лаборатория биологии РНК, Институт биофизики Китайской академии наук

Email: d.pervouchine@skoltech.ru
Китайская республика, Пекин, 100101

Дмитрий Д. Первушин

Сколковский институт науки и технологий

Автор, ответственный за переписку.
Email: d.pervouchine@skoltech.ru
Россия, Москва, 121205

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

Дополнительные файлы