Исследование данных экспрессии генов микромассивов
Этот пример показывает, как идентифицировать дифференциально экспрессированные гены по данным микромассивов и использует Gene Ontology для определения значительных биологических функций, которые связаны с понижающими и повышающими регуляцию генами.
Введение
Микромассивы содержат олигонуклеотидные или кДНК-зонды для сравнения экспрессионного профиля генов в геномной шкале. Определение того, являются ли изменения экспрессии генов статистически значимыми между различными условиями, например, двумя различными типами опухолей, и определение биологической функции дифференциально экспрессируемых генов являются важными целями в эксперименте с микромассивами.
В данном примере используется общедоступный набор данных, содержащий данные экспрессии генов 42 опухолевых тканей эмбриональной центральной нервной системы (ЦНС) [1]. Файлы CEL можно загрузить здесь. Выборки гибридизовали на массивах Affymetrix ® HuGeneFL GeneChip ®. Необработанный набор данных был предварительно обработан процедурами Robust Multi-array Average (RMA) и GC Robust Multi-array Average (GCRMA). Для получения дополнительной информации о предварительной обработке микромассивов Affymetrix oligonucleotide, см. «Предварительная обработка данных микромассивов Affymetrix ® на уровне зонда».
Вы будете использовать t-тест и частоту ложных открытий, чтобы обнаружить дифференциально экспрессированные гены между двумя типами опухолей. Кроме того, вы рассмотрим условия Gene Ontology, связанные со значительно повышенными генами.
Загрузка данных выражения
Загрузите файл MAT cnsexpressiondata содержащий три объекта DataMatrix, сопоставленные со значениями экспрессии генов, предварительно обработанными с использованием RMA (expr_cns_rma), GCRMA с максимальной оценкой правдоподобия (expr_cns_gcrma_mle) и GCRMA с эмпирически-байесовской оценкой (expr_cns_gcrma_eb).
load cnsexpressiondata
В каждом объекте DataMatrix каждая строка соответствует набору зондов в массиве, и каждый столбец соответствует выборке. Объект DataMatrix expr_cns_gcrma_eb будет использоваться в этом примере, но также могут использоваться данные любой из двух других переменных выражения.
Извлеките свойства объекта DataMatrix expr_cns_gcrma_eb использование get команда.
get(expr_cns_gcrma_eb)
Name: ''
RowNames: {7129×1 cell}
ColNames: {1×42 cell}
NRows: 7129
NCols: 42
NDims: 2
ElementClass: 'single'
Определите количество генов и количество выборок путем доступа к количеству строк и количеству столбцов объекта DataMatrix соответственно.
nGenes = expr_cns_gcrma_eb.NRows nSamples = expr_cns_gcrma_eb.NCols
nGenes =
7129
nSamples =
42
Отображение между идентификатором набора зондов и соответствующим символом гена обеспечивается как объект Map в файле MAT HuGeneFL_GeneSymbol_Map.
load HuGeneFL_GeneSymbol_Map
Аннотируйте значения выражения в expr_cns_gcrma_eb с соответствующими символами генов путем создания массива ячеек генных символов из объекта Map и установки имен строк объекта Данных Матрицы.
huGenes = values(hu6800GeneSymbolMap, expr_cns_gcrma_eb.RowNames);
expr_cns_gcrma_eb = rownames(expr_cns_gcrma_eb, ':', huGenes);
Фильтрация данных выражения
Многие наборы зондов в этом примере не аннотируются, не выражаются или имеют небольшую изменчивость между выборками. Используйте следующие методы для фильтрации этих генов.
Удалите данные экспрессии генов с помощью пустых символов генов (в этом примере пустые символы помечены как '---').
expr_cns_gcrma_eb('---', :) = [];
Использование genelowvalfilter для фильтрации генов с очень низкими абсолютными значениями экспрессии.
[~, expr_cns_gcrma_eb] = genelowvalfilter(expr_cns_gcrma_eb);
Использование genevarfilter для фильтрации генов с небольшим отклонением между выборками.
[~, expr_cns_gcrma_eb] = genevarfilter(expr_cns_gcrma_eb);
Определите количество генов после фильтрации.
nGenes = expr_cns_gcrma_eb.NRows
nGenes =
5669
Идентификация дифференциальной экспрессии генов
Теперь можно сравнить значения экспрессии генов между двумя группами данных: опухолью злокачественных глиом (Mglio) ЦНС (MD) и ненейронным источником.
Из данных экспрессии всех 42 выборок в наборе данных извлеките данные 10 выборок MD и 10 выборок Mglio.
MDs = strncmp(expr_cns_gcrma_eb.ColNames,'Brain_MD', 8); Mglios = strncmp(expr_cns_gcrma_eb.ColNames,'Brain_MGlio', 11); MDData = expr_cns_gcrma_eb(:, MDs); get(MDData) MglioData = expr_cns_gcrma_eb(:, Mglios); get(MglioData)
Name: ''
RowNames: {5669×1 cell}
ColNames: {1×10 cell}
NRows: 5669
NCols: 10
NDims: 2
ElementClass: 'single'
Name: ''
RowNames: {5669×1 cell}
ColNames: {1×10 cell}
NRows: 5669
NCols: 10
NDims: 2
ElementClass: 'single'
Проведите t-тест для каждого гена, чтобы идентифицировать значительные изменения в значениях экспрессии между выборками MD и выборками Mglio. Можно просмотреть результаты теста из нормального квантильного графика t-значений и гистограмм t-значений и p-значений t-тестов.
[pvalues, tscores] = mattest(MDData, MglioData,... 'Showhist', true', 'Showplot', true);
В любой тестовой ситуации могут возникнуть два типа ошибок, ложноположительный путем объявления, что ген дифференциально экспрессируется, когда это не так, и ложноотрицательный, когда тест не идентифицирует действительно дифференциально экспрессируемый ген. В нескольких проверки гипотезы, которые одновременно тестируют нулевую гипотезу тысяч генов, каждый тест имеет специфическую ложноположительную частоту или частоту ложных открытий (FDR). Скорость ложного обнаружения определяется как ожидаемое отношение количества ложных срабатываний к общему количеству положительных вызовов в дифференциальном анализе выражения между двумя группами выборок [2].
В этом примере FDR вычисляется с помощью процедуры Storey-Tibshirani [2]. Процедура также вычисляет q-значение теста, который измеряет минимальный FDR, который происходит при вызове значимого теста. Оценка FDR зависит от действительно нулевого распределения нескольких тестов, что неизвестно. Сочетания методы могут использоваться, чтобы оценить истинное нулевое распределение тестовой статистики путем перестановки столбцов матрицы данных экспрессии генов [2] [3]. В зависимости от размера выборки может оказаться невозможным рассмотреть все возможные сочетания. Обычно случайное подмножество сочетаний рассматривается в случае большого размера выборки. Используйте nchoosek функция в Statistics and Machine Learning Toolbox™, чтобы узнать количество всех возможных сочетаний выборок в этом примере.
all_possible_perms = nchoosek(1:MDData.NCols+MglioData.NCols, MDData.NCols); size(all_possible_perms, 1)
ans =
184756
Выполните сочетание t-тест с помощью mattest и PERMUTE опция вычисления значений p 10 000 сочетаний путем перестановки столбцов матрицы данных экспрессии генов MDData и MglioData [3].
pvaluesCorr = mattest(MDData, MglioData, 'Permute', 10000);
Определите количество генов, которые считаются имеющими статистическую значимость, при срезе p-значения 0,05. Примечание: Вы можете получить разное количество генов из-за результата теста сочетания.
cutoff = 0.05; sum(pvaluesCorr < cutoff)
ans =
2121
Оцените FDR и q-значения для каждого теста с помощью mafdr. Количество pi0 является общей долей истинных нулевых гипотез в исследовании. Он оценивается из моделируемого распределения null через bootstrap или кубическую аппроксимацию полиномом. Примечание: Вы также можете вручную задать значение лямбды для оценки pi0.
figure;
[pFDR, qvalues] = mafdr(pvaluesCorr, 'showplot', true);
Определите количество генов, которые имеют q-значения меньше, чем значение среза. Примечание: Вы можете получить разное количество генов из-за теста сочетания и результатов bootstrap.
sum(qvalues < cutoff)
ans =
2173
Многие гены с низким FDR подразумевают, что две группы, MD и Mglio, биологически различны.
Можно также эмпирически оценить скорректированные по FDR значения p с помощью процедуры Бенджамини-Хохберга (BH) [4] путем установки mafdr входной параметр BHFDR к true.
pvaluesBH = mafdr(pvaluesCorr, 'BHFDR', true);
sum(pvaluesBH < cutoff)
ans =
1374
Можно хранить значения t, p-значения, pFDR, q-значения и скорректированные p-значения BH FDR вместе в качестве объекта DataMatrix.
testResults = [tscores pvaluesCorr pFDR qvalues pvaluesBH];
Обновите имя столбца для исправленных значений p BH FDR с помощью colnames метод объекта DataMatrix.
testResults = colnames(testResults, 5, {'FDR_BH'});
Можно отсортировать по p-значениям pvaluesCorr использование sortrows mathod.
testResults = sortrows(testResults, 2);
Отобразите первые 20 генов в testResults. Примечание. Ваши результаты могут отличаться от показанных ниже из-за теста сочетания и результатов bootstrap.
testResults(1:20, :)
ans =
t-scores p-values FDR q-values FDR_BH
PLEC1 -9.6223 6.7194e-09 1.3675e-05 7.171e-06 1.9974e-05
HNRPA1 9.359 1.382e-08 1.4063e-05 7.171e-06 1.9974e-05
FCGR2A -9.3548 1.394e-08 9.457e-06 7.171e-06 1.9974e-05
PLEC1 -9.3495 1.4094e-08 7.171e-06 7.171e-06 1.9974e-05
FBL 9.1518 1.9875e-08 8.0899e-06 7.1728e-06 1.998e-05
KIAA0367 -8.996 2.4324e-08 8.2509e-06 7.1728e-06 1.998e-05
ID2B -8.9285 2.6667e-08 7.7533e-06 7.1728e-06 1.998e-05
RBMX 8.8905 2.8195e-08 7.1728e-06 7.1728e-06 1.998e-05
PAFAH1B3 8.7561 3.5317e-08 7.9864e-06 7.9864e-06 2.2246e-05
H3F3A 8.6512 4.5191e-08 9.1973e-06 8.5559e-06 2.3832e-05
LRP1 -8.6465 4.6243e-08 8.5559e-06 8.5559e-06 2.3832e-05
PEA15 -8.3256 1.1419e-07 1.9367e-05 1.9367e-05 5.3947e-05
ID2B -8.1183 1.7041e-07 2.6679e-05 2.4793e-05 6.9059e-05
SFRS3 8.1166 1.7055e-07 2.4793e-05 2.4793e-05 6.9059e-05
HLA-DPA1 -7.8546 2.4004e-07 3.2569e-05 3.2569e-05 9.072e-05
C5orf13 7.7195 2.9229e-07 3.7179e-05 3.3475e-05 9.3243e-05
PTMA 7.7013 2.9658e-07 3.5506e-05 3.3475e-05 9.3243e-05
NAP1L1 7.674 3.0489e-07 3.4474e-05 3.3475e-05 9.3243e-05
HMGB2 7.6532 3.1251e-07 3.3475e-05 3.3475e-05 9.3243e-05
RAB31 -13.664 3.308e-07 3.3662e-05 3.3662e-05 9.3766e-05
Считается, что ген дифференциально экспрессируется между двумя группами выборок, если он показывает как статистическую, так и биологическую значимость. Этот пример сравнивает отношение экспрессии генов MD по сравнению с выборками опухоли Mglio. Поэтому повышенная регуляция гена в этом примере имеет более высокую экспрессию в MD, и пониженная регуляция гена имеет более высокую экспрессию в Mglio.
Постройте график -log10 значений p относительно биологического эффекта на графике вулкана. Примечание: Из графика вулкана UI, вы можете в интерактивном режиме изменить предел отключения и сворачивания значения p и экспортировать дифференциально экспрессированные гены.
diffStruct = mavolcanoplot(MDData, MglioData, pvaluesCorr)
diffStruct =
struct with fields:
Name: 'Differentially Expressed'
PVCutoff: 0.0500
FCThreshold: 2
GeneLabels: {327×1 cell}
PValues: [327×1 bioma.data.DataMatrix]
FoldChanges: [327×1 bioma.data.DataMatrix]
Щелкните по генам Ctrl в списках генов, чтобы пометить гены на графике. Как видно на графике вулкана, гены, специфичные для клеток гранул церебеллы на основе нейронов, таких как ZIC и NEUROD, находятся в списке регулируемых выше генов, в то время как гены, типичные для астроцитарной и олигодендроцитарной линии и дифференцировки клеток, такие как SOX2, PEA15 и ID2B, находятся в списке регулируемых вниз.
Определите количество дифференциально экспрессированных генов.
nDiffGenes = diffStruct.PValues.NRows
nDiffGenes = 327
В частности, определить список повышенных генов и список пониженных генов для MD по сравнению с Mglio.
up_geneidx = find(diffStruct.FoldChanges > 0); nUpGenes = length(up_geneidx) down_geneidx = find(diffStruct.FoldChanges < 0); nDownGenes = length(down_geneidx)
nUpGenes = 225 nDownGenes = 102
Аннотирование регулируемых генов с использованием генной онтологии
Можно использовать информацию Gene Ontology (GO) для аннотирования дифференциально экспрессированных генов, идентифицированных выше. Файл аннотации GO (GAF) для человека может быть загружен из Gene Ontology Current Annotations. Для удобства, карта между символами гена и связанными идентификаторами GO относительно поля аспекта Function включен в файл MAT goa_human.
load goa_human
Также можно запустить код ниже, чтобы загрузить базу данных Gene Ontology с последними аннотациями, прочитать загруженный файл аннотации Homo sapiens (назовите файл как goa_human.gaf) и создайте отображение между символами гена и связанными терминами GO.
% GO = geneont('live',true); % HGann = goannotread('goa_human.gaf',... % 'Aspect','F','Fields',{'DB_Object_Symbol','GOid'}); % HGmap = containers.Map(); % for i = 1:numel(HGann) % key = HGann(i).DB_Object_Symbol; % if isKey(HGmap,key) % HGmap(key) = [HGmap(key) HGann(i).GOid]; % else % HGmap(key) = HGann(i).GOid; % end % end
Найдите индексы повышенных генов для анализа генной онтологии.
up_genes = rownames(diffStruct.FoldChanges, up_geneidx); huGenes = rownames(expr_cns_gcrma_eb); for i = 1:nUpGenes up_geneidx(i) = find(strncmpi(huGenes, up_genes{i}, length(up_genes{i})), 1); end
Не все гены на чипе HuGeneFL аннотированы. Для каждого гена на чипе смотрите, аннотирован ли он путем сравнения его генного символа со списком генных символов от GO. Отслеживайте количество аннотированных генов и количество повышенных генов, сопоставленных с каждым термином GO. Обратите внимание, что данные в общедоступных хранилищах часто кураторятся и обновляются; поэтому результаты этого примера могут несколько отличаться при использовании современных наборов данных. Также возможно, что вы получаете предупреждения о недопустимых или устаревших идентификаторах из-за обновленного файла аннотации гена Homo sapiens.
m = GO.Terms(end).id; % gets the last term id chipgenesCount = zeros(m,1); % a vector of GO term counts for the entire chip. upgenesCount = zeros(m,1); % a vector of GO term counts for up-regulated genes. for i = 1:length(huGenes) if isKey(HGmap,huGenes{i}) goid = getrelatives(GO,HGmap(huGenes{i})); chipgenesCount(goid) = chipgenesCount(goid) + 1; if (any(i == up_geneidx)) upgenesCount(goid) = upgenesCount(goid) + 1; end end end
Определите статистически значимые члены GO с помощью гипергеометрического распределения вероятностей. Для каждого члена GO вычисляется p-значение, представляющее вероятность того, что количество аннотированных генов, связанных с ним, могло быть найдено случайно.
gopvalues = hygepdf(upgenesCount,max(chipgenesCount),...
max(upgenesCount),chipgenesCount);
[dummy, idx] = sort(gopvalues);
Сообщите о десяти наиболее значимых условиях GO следующим образом.
report = sprintf('GO Term p-value counts definition\n'); for i = 1:10 term = idx(i); report = sprintf('%s%s\t%-1.5f\t%3d / %3d\t%s...\n',... report, char(num2goid(term)), gopvalues(term),... upgenesCount(term), chipgenesCount(term),... GO(term).Term.definition(2:min(50,end))); end disp(report);
GO Term p-value counts definition GO:0005515 0.00000 131 / 3459 Interacting selectively and non-covalently with a... GO:0044822 0.00000 94 / 514 Interacting non-covalently with a poly(A) RNA, a ... GO:0003723 0.00000 95 / 611 Interacting selectively and non-covalently with a... GO:0003729 0.00000 82 / 460 Interacting selectively and non-covalently with m... GO:0003735 0.00000 54 / 159 The action of a molecule that contributes to the ... GO:0019843 0.00000 48 / 186 Interacting selectively and non-covalently with r... GO:0008135 0.00000 50 / 208 Functions during translation by interacting selec... GO:0000049 0.00000 47 / 188 Interacting selectively and non-covalently with t... GO:0000498 0.00000 46 / 179 Interacting selectively and non-covalently with r... GO:0001069 0.00000 46 / 179 Interacting selectively and non-covalently with a...
Выберите условия GO, связанные с конкретными функциями молекул, и создайте субонтологию, которая включает предков членов. Визуализируйте эту онтологию с помощью biograph функция. Можно окрасить графики узлы в соответствии с их значимостью. В этом примере красные узлы являются наиболее значимыми, в то время как синие узлы являются наименее значимыми терминами онтологии генов. Примечание: Возвращенные условия GO могут отличаться от показанных из-за частого обновления файла аннотации гена Homo sapiens.
fcnAncestors = GO(getancestors(GO,idx(1:5))); [cm,acc,rels] = getmatrix(fcnAncestors); BG = biograph(cm,get(fcnAncestors.Terms,'name')); for i = 1:numel(acc) pval = gopvalues(acc(i)); color = [(1-pval).^(1) pval.^(1/8) pval.^(1/8)]; BG.Nodes(i).Color = color; end view(BG)
Нахождение дифференциально экспрессированных генов в путях
Можно запросить информацию о путях дифференциально экспрессированных генов из базы данных путей KEGG через веб-сервис KEGG.
Ниже приведены несколько карт путей с генами в списке регулируемых генов, выделенных:
Ссылки
[1] Pomeroy, S.L., et al., «Предсказание исхода эмбриональной опухоли центральной нервной системы на основе экспрессии генов». Природа, 415 (6870): 436-42, 2001.
[2] Storey, J.D., and Tibshirani, R., «Statistical valuance for genomewide studies», PNAS, 100 (16): 9440-5, 2003.
[3] Dudoit, S., Shaffer, J.P., and Boldrick, J.C., «Множественная гипотеза, проверка в микромассив эксперименте», Statistical Science, 18 (1): 71-103, 2003.
[4] Benjamini, Y., and Hochberg, Y., «Controlling the false discovery rate: a practical and мощный подход к нескольким проверки», Journal of the Royal Statistical Society, Series B, 57 (1): 289-300, 1995.