bioinformatics.ru

Статистический анализ данных для биомедицины

R + DAVID + R + REViGO = функциональная классификация генов

Обновлено 01.05.2012 Автор: Сергей Мастицкий

Пример функциональной аннотации списков генов с помощью Gene Onthology, DAVID, REViGO и R.

Перепост из блога Сергея Мастицкого R: Анализ и визуализация данных (рекомендую).

Оригинал здесь - http://r-analytics.blogspot.com/2012/03/r-david-r-revigo.html

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

Входные данные такого типа обычно представляют собой огромную таблицу, строки которой соответствуют генам (несколько десятков тысяч), а столбцы - исследуемым образцам. На пересечениях строк и столбцов находятся числа, соответствующие уровням экспрессии генов. Типичная задача, которая обычно ставится перед аналитиком - выявить гены, уровни экспрессии которых различаются в экспериментальных группах. Уже несколько лет как de facto стандартом для выполнения подобного анализа является программное обеспечение, написанное на языке R. Речь, в частности, идет о большом количестве пакетов для R, созданных в рамках проекта Bioconductor. Отсюда первая буква R в названии этого сообщения - R используется для анализа изменения экспрессии генов. Примеры того, как имено это делается, я планирую привести в будущих сообщениях. Здесь же я хочу рассмотреть небольшие примеры того, что можно сделать со списками дифференциально экспрессированных генов уже после их обнаружения.

Предположим, в исследовании были задействованы три группы (контроль + две экспериментальные). Для каждой из этих групп при помощи R-пакетов из Bioconductor групп выявлены гены, уровни экспрессии которых статистически значимо отличались от контроля как в положительную (т.е. гены с повышенной экспрессией, "upregulated genes"), так и в отрицательную сторону (гены с пониженной экспрессией, "downregulated genes"). Следовательно, имеем четыре списка генов. Вопрос: что делать с этими списками дальше? Можно, конечно, попробовать проанализировать их "вручную" с целью обнаружить особенно интересные гены ("интересные" с точки зрения того, что о них уже известно в литературе) и затем сделать соответствующие ситуации выводы. Но если списки генов большие, сделать это будет весьма затруднительно. Возникает необходимость каким-то образом автоматизировать процесс извлечения полезной информации из таких списков. Для этого молекулярными биологами и биоинформатиками разработаны и постоянно обновляются специализированные базы данных, запросы из которых можно делать в режиме "онлайн".

Одна из наиболее известных таких баз данных - DAVID (Database for Annotation, Visualization and Integrated Discovery; Huang et al. 2009). DAVID предоставляет широкие возможности для анализа интересующих исследователя списков генов. Предположим, стоит задача выяснить, являются ли полученные нами в ходе предыдущего анализа списки генов лишь случайными наборами с точки зрения биологических процессов, в которых задействованы этих гены. Или так: содержат ли наши списки генов неслучайно (в статистическом смысле слова) большое число генов, связанных с выполнением определенных биологических функций?

Исторически сложилось несколько систем, задающих генам уникальные идентификаторы (например, системы GeneBank, Entrez, Ensembl и др.). База данных DAVID может принимать практически все существующие в настоящее время идентификаторы. Представим, что в ходе нашего анализа было получено четыре списка, содержащих идентификаторы генов из GeneBank. Пусть один из этих списков содержит 109 элементов (если вы намерены повторить описываемые мной примеры, вы можете скачать этот список в виде текстового файла здесь). Идентификаторы выглядят как NM_001081005, NM_138752, NM_019572 и т.д.

Заходим на сайт DAVID и выбираем опцию Start Analysis (Начать анализ) из главного меню:


В левой части страницы сайта появится поле Upload Gene List (Вставить список генов), в которое следует вставить скопированные из нашего текстового файла идентификаторы генов:


В качестве следующего шага следует "объяснить" программе, какую именно систему идентификаторов генов мы используем. Поэтому в выпадающем меню, расположенном ниже нашего списка (Step 2. Select Identifier), выбираем GENEBANK_ACCESSION:


Наконец, следует указать тип нашего списка (Step 3. List Type) - является ли он простым списком генов, выявленных в ходе определенного анализа (Gene List), или, например, списком всех генов, известных для изучаемого организма (Background). Выбираем первую опцию и нажимаем кнопку "Submit list":


В результате откроется новая страница, предлагающая несколько видов анализа. Нас интересует функциональный анализ генов, поэтому выбираем опцию Functional Annotation Tool (Инструмент для функциональной аннотации):


Функциональная аннотация генов - понятие широкое. Нас, в частности, интересует словарь биологических процессов из системы Gene Ontology (Генная Онтология - далее по тексту "ГО"). Поэтому на открывшейся странице выбираем Gene_Ontology, а затем жмем на кнопку "Chart" (Таблица) напротив опции GEOTERM_BP_5 (я опускаю здесь подробности по поводу выбора этой опции, оставляя Читателю возможность самому разобраться в тонкостях ГО):



Открывшееся новое окно браузера будет содержать таблицу, включащую список функционаьных категорий из ГО (столбец Term) и много другой информации. Для нас больше всего на данный момент интересны столбцы Term, а также P-value и Benjamini. Как следует из названия, в столбце P-value представлены P-значения для модифицированного точного теста Фишера (Fisher exact test), который используется для проверки нулевой гипотезы о том, что наш список генов содержит лишь случайное количество элементов, ассоциированных с конкретной функциональной категорией. Поскольку соответствующих категорий в ГО было обнаружено много (здесь 81), было выполнено много таких тестов. Как результат, возникла необходимость скорректировать полученные исходные P-значения с учетом количества тестов (т.н. проблема множественных сравнений). Столбец Benjamini содержит эти новые Р-значения, скорректированные по методу Benjamini-Hochberg (Benjamini & Hochberg 1995). Они-то нам и потребуются для дальнейшего анализа. Но как извлечь необходимые нам данные из окна браузера? К счастью, таблицу с результатами функционального анализа можно сохранить локально в виде текстового файла (правая кнопка мыши > Сохранить как...):



(Примечание: Если вы повторяете описываемые примеры, сохраните файл с именем, которое будет содержать фразу "_GO_terms", например, "List_1_GO_terms.txt" - пригодится в дальнейшем).

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

Полученные текстовые файлы с результатами ГО-анализа  удобнее всего открыть при помощи Excel. Сделав это, и обратив внимание на столбец Term, мы увидим, что многие перечисленные в нем функциональные категории синонимичны. Так, в примере ниже категории "GO:0050863~regulation of T cell activation", "GO:0042110~T cell activation" и "GO:0030217~T cell differentiation", по сути, говорят об одном и том же биологическом процессе - активации Т-лимфоцитов.



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

До недавнего времени с указанной избыточностью функциональных категорий бороться было достаточно сложно, и приходилось "вручную перелопачивать" результаты ГО-анализа с целью выбрать наиболее представительные категории. Однако ситуация изменилась с появлением нового инструмента - REViGO, что значит "Reduce and Visualize Gene Ontology" (Supek et al. 2011) (дословно "уменьши и визуализируй генную онтологию"). Алгоритм REViGO работает так, что избыточные ГО-термины из списка удаляются, а оставшиеся сходные термины объединяются в еще более крупные функционально-смысловые кластеры, что значительно облегчает их интерпретацию.

Входным типом данных для REViGO служит список ГО-идентификаторов (см. рисунок выше - содержат буквы GO, после которых следует уникальный код из 7 чисел) и соответствующих им Р-значений. Вся эта информация у нас уже имеется - мы получили ее при помощи базы данных DAVID. Но есть одна загвоздка: система DAVID автоматически "приклеила" названия ГО-категорий к их идентификаторам. К сожалению, в таком виде ГО-идентификаторы для REViGO непригодны. Конечно, можно потратить время и вручную очистить идентификаторы от их названий, но это не наш путь! Мы воспользуемся мощью R и напишем небольшой скрипт, который за пару секунд автоматически подготовит пригодные для REViGO файлы из тех DAVID-файлов, которые у нас уже есть. Отсюда вторая буква R в названии данного сообщения...

Итак, еще раз: у нас есть четыре (конечно, их может быть и другое количество - один и более) файла с таблицами, содержащими функциональные ГО-категории и еще много всего. Задача: одновременно загрузить все эти таблицы в R, выбрать из них только нужные нам столбцы (Term и Benjamini - см. выше) и сохранить результат в виде четырех отдельных текстовых файлов. Ниже приведен и прокомментирован код написанных мной для этих целей R-функций.

Файлы, созданные при помощи DAVID, должны находиться в рабочей папке R и содержать элемент "_GO_terms" в имени. Безусловно, это может быть и любой другой элемент - главное, чтобы он повторялся в именах файлов и помогал вам эти файлы идентифицировать. Следовательно, приведенный ниже код можно модифицировать в соответствии с вашими нуждами.

Для начала создадим текстовый вектор, который будет содержать имена загружаемых файлов (если какие-то из встречающихся ниже функций R вам неизвестны - щелкните по их названию, чтобы открыть справку): 
files.to.read = grep("GO_terms", dir(), value=TRUE)
Создадим список, в котором будут храниться загружаемые нами таблицы (сначала он будет пустым - мы заполним его позднее при помощи специально написанной для этого функции):

LIST = list()

Функция readGObatch(), код которой приведен ниже, будет работать следующим образом. У нее будет два аргумента - текстовый вектор (f), содержащий перечень считываемых файлов, и список (L), в который мы загрузим считываемые файлы. Из считанных файлов-таблиц функция readGObatch() извлечет только столбцы с заголовками Term и Benjamini (см. вектор cols.to.keep в коде), а также "очистит" ГО-идентификаторы от ненужного текста (см. команду c использованием функции substring() ). Функция будет возвращать список с уже готовыми "очищенными" таблицами, пригодными для анализа в REViGO:

readGObatch <- function(f, L, ...){
    for(i in 1:length(f)){
        L[[i]] = read.table(f[i], sep="\t", header=T)
        cols.to.keep = c(which(names(L[[i]]) == "Term"),
                         which(names(L[[i]]) == "Benjamini"))
    L[[i]] = L[[i]][cols.to.keep]
    L[[i]]$Term = substring(as.character(L[[i]]$Term), 1, 10)
    }
    names(L) = f
    L}

Заполнить созданный ранее список LIST при помощи функции readGObatch() теперь очень просто:

LIST = readGObatch(files.to.read, LIST)

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

filterByPvalue <- function(L, p){
    for(i in 1:length(L)){
        dat = as.data.frame(L[[i]])
        dat = dat[which(dat$Benjamini < p), ]
        L[[i]] = dat
    }
    rm(i)
    L
}

Применим эту новую функцию к списку LIST для отфильтровывания всех ГО-категорий, у которых P-значения, скорректированные по методу Benjamini-Hochberg, превышают, например, 0.10:

LIST = filterByPvalue(L = LIST, p = 0.10)

Теперь нам нужна функция, которая запишет все очищенные таблицы в виде отдельных текстовых файлов, и при этом добавит слово "Clean" ("чистый") в название этих файлов для облегчения их узнавания:

writeGObatch <- function(L){
    for(i in 1:length(L)){
        write.table(L[[i]], file =
            paste("Clean", names(L)[i], sep="_"),
                    sep = "\t", quote = F, row.names = F)
    }
}

Наконец, применим, эту последнюю функцию к нашему списку LIST:

writeGObatch(LIST)

В текущую рабочую папку R при этом будут записаны четыре (их может быть и другое количество) текстовых файла, содержащих таблицы вроде этой:




Приведенные на рисунке два столбца можно теперь скопировать (без заголовков!) и вставить в соответствующее поле на сайте REViGO:



Все настройки можно оставить без изменений и просто нажать кнопку "Start Revigo":



В результате откроется страница с графиком в духе Multidimensional Scaling (MDS) и расположенной ниже таблицей, в которой перечислены анализируемые ГО-идентификаторы и некоторые статистические параметры, рассчитанные алгоритмом REViGO (подробнее см. в оригинальной статье Supek et al. 2011). Хотя графика MDS уже вполне достаточно, чтобы сделать необходимые выводы, я лично предпочитаю другой тип графического представления результатов этого анализа - т.н. TreeMap:


Результат приведен ниже:



Из этого рисунка мгновенно видно, что анализируемый список "насыщен" (англ. enriched) генами, связанными с выполнением двух крупных биологических функций - активацией T-клеток и регуляцией передачи клеточных сигналов. Размер прямоугольников пропорционален статистической занчимости той или иной категории (если быть точным, то -log10(P-value)), так что активация Т-лимфоцитов, как видим, является более значимой, чем регуляция передачи сигналов.

Обратите внимание на то, что из большого количества ГО-категорий, поданных для обработки изначально (39 в рассматриваемом примере - после фильтрации по Р-значениям), алгоритм REViGO отфильтровал большинство из них, оставив только 4 неизбыточных категории,  которые были объединены в два функционально-смысловых кластера. Далее исследователь может вернуться к таблицам, выданным системой DAVID, и изучить внимательнее гены из "подсказанных" алгоритмом REViGO ГО-категорий. Проведя подобный анализ для всех списков дифференциально экспрессируемых генов, исследователь в итоге получит достаточно полное представление о биологических процессах, за которые эти гены отвечают.

< Предыдущая   Следующая >

Из последнего

Опрос

Какие материалы вам бы хотелось видеть на сайте?

Итоги