Анализ взвешенных сетей коэкспрессии генов

Анализ взвешенной сети корреляций (англ. weighted correlation network analysis), также известный как анализ взвешенной сети коэкспрессии генов (англ. weighted gene coexpression network analysis, WGCNA), это широко используемый, особенно при изучении биологических сетей, метод глубинного анализа данных, основанный на попарных корреляциях между переменными. Метод может быть использован на большинстве многомерных наборов данных, но наиболее широкое распространение он получил в геномике. Метод позволяет определять модули (кластеры коэкспрессирующихся генов), межмодульные хабы и узлы сети относительно принадлежности к модулю, изучать отношения между модулями коэкспрессии и сравнивать топологии различных сетей. WGCNA может быть использован как метод снижения размерности данных (связанный с непрямым факторным анализом), как метод кластеризации, как метод отбора признаков (например, для скрининга генов).

История

Метод WGCNA был разработан Стивом Хорватом, профессором Калифорнийского университета в Лос-Анджелесе и членами (в том числе бывшими) его лаборатории (в частности, Питером Лангфелдером, Бином Чжаноми и Джуном Донгом). Большая часть метода появилась в ходе прикладных исследований. В частности, взвешенные корреляционные сети были разработаны в рамках совместных дискуссий с исследователями рака Полом Мишелем, Стэнли Ф. Нельсоном, и нейробиологами Дэниэлем Гешвиндом и Майклом Олдхэмом.

Алгоритм работы

Построение первичной сети

Исходные данные представляют собой матрицу экспрессии генов в нескольких образцах $X_{n\times m}=[x_{il}]$, где $n$ — число генов, $m$ — число образцов. Для построения сети коэкспрессии для каждой пары генов $(i,j)$ вычисляется параметр сходства, $s_{ij}$, который должен лежать в $[0;1]$. $s_{ij}$ по умолчанию равен абсолютному значению коэффициента корреляции Пирсона: 

$s_{ij}=|cor(x_{i},x_{j})|$,

где профили коэкспрессии генов $x_{i}$ и $x_{j}$ состоят из экспрессии генов i и j среди множества образцов. Для достижения большей устойчивости к выбросам можно использовать “урезанный” коэффициент корреляции, а для того, чтобы сохранить знак коэффициента можно использовать простое преобразование корреляции: $s_{ij}={\frac {1+cor(x_{i},x_{j})}{2}}$, так как использование абсолютного значения корреляции может привести к потере биологически значимой информации, поскольку при этом нельзя различить репрессию и активацию генов. Получается матрица сходства $S_{n\times n}=[s_{ij}]$.

Фильтрация ребер по весу с применением мягкого безмасштабного критерия

Для вычисления матрицы смежности $A_{n\times n}=[a_{ij}]$ сети по матрице сходства необходима функция смежности, которая отображает интервал $[0;1]$ в интервал $[0;1]$. Традиционной функцией смежности является сигнум-функция с жестким порогом:

$a_{ij}=signum(a_{ij},\tau )\equiv {\begin{cases}1,&s_{ij}\geq \tau \\0,&s_{ij}<\tau \end{cases}}$

Такая жесткая фильтрация ребер применяется при построении невзвешенных сетей, а результат согласован с интуитивным пониманием концепции сети (связность узла совпадает с числом связанных с ним соседей). Однако такой подход часто приводит к потере информации: например, если установить $\tau =0.8$, то между узлами с параметром сходства 0.79 не будет никакой связи. Таким образом такие сети очень чувствительны к выбору гиперпараметра.

Взвешенные сети лишены этого недостатка. В алгоритме WGCNA в качестве функции смежности используется степенная функция:

$a_{ij}=power(s_{ij},\beta )\equiv |s_{ij}|^{\beta }$

Значение $\beta$ определяется с помощью топологического критерия безмасштабности сети для целочисленных значений $\beta$. Для разных целочисленных значений $\beta$ строится линейная модель зависимости логарифма доли вершин сети со степенью $k$ $\log(p(k))$ от логарифма степени $log(k)$. Выбирается наименьшее значение $\beta$, при котором коэффициент детерминации $R^{2}$ соответствующей линейной модели превосходит 0.8.

«Мягкая» матрица смежности позволяет лишь проранжировать узлы сети согласно силе их связи с рассматриваемым узлом. Если необходимо определить ограниченный список соседей, вводится порог по силе связи. Такой способ «мягкой» фильтрации ребер сети и называется взвешиванием сети.

Сглаживание взвешенной сети

Для определения модулей сети и удаления шума производится операция сглаживания сети. Сначала вычисляется матрица топологического сходства (topological overlap measure) $\Omega _{n\times n}=[\omega _{ij}]$:

$\omega _{ij}={\frac {l_{ij}+a_{ij}}{\min\{k_{i},k_{j}\}+1-a_{ij}}}$,

где $a_{ij}$ — элемент матрицы смежности, $l_{ij}=\sum _{u}a_{iu}a_{uj}$, $k_{i}=\sum _{u}a_{iu}$.

$\omega _{ij}$ равна $1$, если узел $i$ с меньшим числом соседей связан с узлом $j$ с бо́льшим числом соседей и все соседи узла $i$ являются соседями узла $j$. $\omega _{ij}$ равна $0$, если узлы $i$ и $j$ не связаны и не имеют общих соседей. Так как $0\leq a_{ij}\leq 1$, то и $0\leq \omega _{ij}\leq 1$. Топологическое сходство двух узлов отражает их относительную внутрисвязанность.

На основе матрицы топологического сходства строится матрица несходства $D_{n\times n}=[d_{ij}^{\omega }]$:

$d_{ij}^{\omega }=1-\omega _{ij}$.

Модули сети выделяются согласно матрице несходства.

Иерархическая кластеризация и выделение модулей

Для выделения модулей сети используется иерархическая кластеризация матрицы несходства. Модули генов представляют собой ветви полученной дендрограммы. Для выделения модулей производится обрезка ветвей дендрограммы. Статический способ обрезки, который выделяет разветвления ниже определенного порога как отдельный кластер, производит модули, которые легко выделяются визуально, но не соответствуют строгому определению понятия «модуль».

В WGCNA используется метод динамического обрезания дендрограммы (англ. Dynamic Tree Cut), который представляет собой адаптивный итеративный процесс разбиения и комбинации модулей, и останавливается, когда число модулей остается стабильным. Исходным набором модулей являются большие модули, определяемые статическим способом. Далее происходит рекурсивное разбиение на подмодули на основе наличия паттерна флуктуаций длины ветвей внутри одного исходного модуля. Маленькие модули объединяются с соседними, чтобы избежать чрезмерного дробления.

Более поздние версии алгоритма используют гибридный динамический (англ. Dynamic Hybrid) подход, который позволяет более успешно детектировать выбросы в каждом кластере. В качестве зачатков модулей используются ветви, удовлетворяющие следующим критериям:

• Содержат не менее установленного порога дочерних листьев (объектов);
• Объекты, расположенные слишком далеко от модуля исключаются из него, даже если они принадлежат одной дочерней ветви на дендрограмме;
• Каждый модуль должен быть различимым от его окружения;
• Ядро модуля должно быть сильно связано.

После определения таких зачаточных модулей все остальные объекты, не попавшие в модулей на первом шаге, по возможности включаются в состав образованных модулей методом ближайших медоидов. Объекты в составе одного итогового модуля могут не находиться рядом на дендрограмме, однако эта несогласованность является обманчивой и следует из ограничений отображения реальных данных в виде дендрограммы.

Валидация модулей

Для валидации модулей используется бутстрэп-анализ на неполных выборках образцов. Дальнейшему анализу подвергаются только те модули, которые были обнаружены в большом числе бутстрэп-реплик.

Интерпретация результатов

Обобщение профиля экспрессии генов модуля

Модуль коэкспрессии $q$ представляется собственным геном (англ. eigengene) $E^{(q)}$, который является правым сингулярным вектором, соответствующим наибольшему правому сингулярному значению при сингулярном разложении матрицы экспрессии генов этого модуля $X_{n^{(q)}\times m}^{(q)}=[x_{ij}^{(q)}]$.

Соотнесение модулей с внешними характеристиками

Для определения значимости влияния экспрессии одного гена $x_{i}$ на проявление черты $T=(T_{1},...,T_{m})$ можно использовать модуль коэффициента корреляции $T:GS_{i}=|cor(x_{i},T)|$ или p-value соответствующего корреляционного или регрессионного анализа $GS_{i}=-\log p_{i}$. Для определения значимости влияния экспрессии генов модуля $q$ используются разные величины:

• Среднее значение модулей коэффициентов корреляции между экспрессиями генов и проявления черты по данному модулю $T:GS^{(q)}={\frac {1}{n^{(q)}}}\sum _{i:x_{i}\in q}^{n^{(q)}}T:GS_{i}$;
• Модуль коэффициента корреляции собственного гена и проявления черты $T:GS^{(q)}=|cor(E^{(q)},T)|$ и соответствующий p-value $GS^{(q)}=-\log p^{(q)}$ регрессионного анализа.

Чем больше величина значимости, тем более значимым является данный ген $x_{i}$ или модуль генов $q$ для проявления данной черты $T$.

Взаимодействие модулей

Для оценки взаимодействия (коэкспрессии) модулей используются как попарные корреляции между собственными генами модулей, так и построение мета-сети на основе матрицы коэкспрессии собственных генов модулей и выделение мета-модулей в этой мета-сети, в том числе с использованием информации о проявлении черт для определения значимости модулей.

Выделение важных генов

Одним из ключевых моментов в анализе построенной сети является выделение центральных узлов (хабов) — генов, которые высоко скоррелированы со многими другими генами сети. Две метрики связности узла: стандартная $k_{i}=\sum _{j=1}^{n}a_{ij}$ и топологическая $\omega _{i}=\sum _{j=1}^{n}\omega _{ij}$ — могут быть вычислены как по всей сети (полносетевая связность, англ. whole-network connectivity), так и только по генам того же модуля (внутримодальная связность, англ. intramodular connectivity). Стандартная и топологическая полносетевая и внутримодальная связности обозначаются как $kTotal$, $\omega Total$, $kIM$ и $\omega IM$, соответственно.

Использование внутримодульной связности (нормированной на максимальное значение связности внутри модуля) является более предпочтительным по сравнению с полносетевой связностью при сравнении генов из различных модулей, так как сильно связный ген но из маленького по размеру, но важного модуля может иметь гораздо меньшее значение полносетевой связности по сравнению с геном, имеющим среднюю степень связности из большого по размеру, но не очень важного модуля.

В свою очередь, внутримодульная связность гена отражает то, насколько связан или коэкспрессирован рассматриваемый ген по отношению к генам этого модуля. Эта метрика может быть интерпретирована как мера принадлежности к модулю.

Также выделяют метрики $kOut=kTotal-kIM$ и $kDiff=kIn-kOut$. Высокое среднее значение $kIM$ для генов модуля (даже превышающее среднее значение $kTotal$) может означать то, что этот модуль является ядром сети и важным с биологической точки зрения в изучаемых условиях. Высокое значение $kOut$, напротив, означает, что общая связанность больше, чем связанность внутри модуля, то есть гены, входящие в состав этого модуля, стабильны и слабо меняют экспрессию в исследуемых условиях.

Нечёткая мера принадлежности к модулю

Бинарная мера принадлежности к модулю может быть не самой подходящей мерой для ряда применений. Непрерывная величина может служить нечёткой мерой принадлежности гена к модулю. Такой величиной может быть $kIM$ или мера связности, основанная на собственных векторах модулей, $kME$. Последняя определяется как корреляция гена $x_{i}$ и собственного гена модуля $q$ $E^{(q)}$: $kME\equiv K_{cor,i}^{(q)}:=cor(x_{i},E^{(q)})$. Высокое значение $kME$ (близкое к $1$ или $-1$) говорит о сильной связанности гена $x_{i}$ и модуля $q$, значение, близкое к $0$, говорит о том, что ген $x_{i}$ не принадлежит модулю $q$. Метрики $kIM$ и $kME$ связанны между собой: внутримодальные хабы обладают тенденцией иметь высокое значение $kME$ для соответствующего модуля.

Сравнение взвешенных и невзвешенных сетей корреляций

Взвешенная корреляционная сеть может рассматриваться как частный случай взвешенной сети, сети зависимостей или корреляционной сети. Анализ взвешенных корреляционных сетей имеет следующие достоинства:

• Построение сети (на основе мягкого порога коэффициента корреляции) сохраняет непрерывный характер исходной информации о корреляции. Например, взвешенные корреляционные сети, построенные на основе корреляций между числовыми переменными не требуют выбора жесткого порога. Дихотомическое деление информации и (жесткий) выбор порога может привести к потере информации;
• Взвешенные корреляционные сети облегчают геометрическую интерпретацию на основе угловой интерпретации корреляции;
• Полученная статистика сети может быть использована для улучшения результатов стандартных методов глубокого анализа данных, таких как кластерный анализ, так как меры несхожести зачастую могут быть преобразованы во взвешенные сети;
• WGCNA предоставляет обширные статистические данные о консервативности модулей, которые могут быть использованы для количественного представления возможности встретить их в другом состоянии. Также статистика консервативности модулей позволяет исследовать различия между модульной структурой сетей;
• Взвешенные и корреляционные сети часто можно аппроксимировать «факторизуемыми» сетями. Таких приближений бывает сложно добиться для разреженных, невзвешенных сетей. Таким образом, взвешенные (корреляционные) сети позволяют использовать экономную параметризацию (в терминах модулей и принадлежности к модулям) .

Применение

WGCNA широко применяется для анализа данных об экспрессии генов, например для поиска межмодульных хабов.

Эта методика часто используется в качестве шага снижения размерности данных в применении к системной генетике, где модули представлены собственными генами. Собственные гены модулей могут быть использованы для расчёта корреляции с экспериментальными данными.

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

WGCNA можно применять для метаанализа данных, полученных из разных источников.

Программная реализация алгоритма

Пакет программного обеспечения R

Пакет WGCNA программного обеспечения R предоставляет функции для осуществления всех элементов анализа взвешенных сетей (построение модулей, выбор генов хабов, статистика консервативности модулей, дифференциальный анализ сети, статистика сети). Пакет WGCNA доступен в CRAN, стандартном репозитории пакетов дополнений для R, однако для работы WGCNA необходимы пакеты, доступные в репозитории Bioconductor. Пакету посвящён отдельный сайт, где опубликованы краткие руководства с демонстрацией возможностей пакета, а также существует блог одного из авторов метода, в котором публикуются тематические статьи, руководства и новости пакета.

Экспорт результатов для анализа в дочерних программах

Пакет содержит функции для визуализации полученных сетей в VisANT и Cytoscape. Также разработан пакет anRichment для расширенной функциональной аннотации генов в модулях, выделенных WGCNA.

iterativeWGCNA

Существует расширение базового пакета R под названием iterativeWGCNA, реализованное на языке Python. Расширение позволяет усилить устойчивость детектируемых модулей и уменьшить потерю информации. Для этого гены, не попавшие в модули, а также гены со слабым значением связности к своим модулям снова подвергаются WGCNA-анализу и для них заново определяются модули.