Этот пример показывает, как идентифицировать дифференциально экспрессируемые гены по данным микрочипов, и использует Gene Ontology для определения значимых биологических функций, которые связаны с пониженными и повышенными генами.
Микрочипы содержат олигонуклеотидные или кДНК зонды для сравнения профиля экспрессии генов в геномной шкале. Определение того, являются ли изменения в экспрессии генов статистически значимыми между различными состояниями, например двумя различными типами опухолей, и определение биологической функции дифференциально экспрессируемых генов, являются важными целями в эксперименте с микрочипами.
Для этого примера используют общедоступный набор данных, содержащий данные экспрессии генов 42 опухолевых тканей эмбриональной центральной нервной системы (ЦНС) [1]. Файлы CEL можно загрузить здесь. Образцы гибридизовали на массивах Affymetrix ® HuGeneFL GeneChip ®. Необработанный набор данных был предварительно обработан с помощью процедур Rability Multi-array Average (RMA) и GC Rustive Multi-array Average (GCRMA). Для получения дополнительной информации о предварительной обработке микрочипов Аффиметрикс олигонуклеотида см. Данные микрочипов аффиметрикс ® предварительной обработки на уровне зонда.
Вы будете использовать t-тест и частоту ложных открытий, чтобы обнаружить дифференциально экспрессируемые гены между двумя типами опухолей. Кроме того, вы посмотрите на термины генной онтологии, связанные со значительно повышенными генами.
Загрузка файла 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
Отображение между ID набора зондов и соответствующим генным символом предоставляется как объект Map в файле MAT. HuGeneFL_GeneSymbol_Map.
load HuGeneFL_GeneSymbol_Map
Аннотирование значений выражения в expr_cns_gcrma_eb с соответствующими символами генов путем создания клеточного массива символов генов из объекта Map и установки имен строк объекта Data Matrix.
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
Теперь можно сравнить значения экспрессии генов между двумя группами данных: медуллобластомы ЦНС (MD) и злокачественные глиомы ненейронного происхождения (Mglio) опухоли.
Из данных экспрессии всех 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, возникающий при вызове теста. Оценка БУР зависит от истинно нулевого распределения множественных тестов, которое неизвестно. Методы перестановки могут использоваться для оценки истинно нулевого распределения тестовой статистики путем перестановки столбцов матрицы данных экспрессии генов [2] [3]. В зависимости от размера выборки может оказаться невозможным рассмотреть все возможные перестановки. Обычно случайное подмножество перестановок рассматривается в случае большого размера выборки. Используйте nchoosek в Toolbox™ Статистика и машинное обучение, чтобы узнать количество всех возможных перестановок выборок в этом примере.
all_possible_perms = nchoosek(1:MDData.NCols+MglioData.NCols, MDData.NCols); size(all_possible_perms, 1)
ans =
184756
Выполнение t-теста перестановки с использованием mattest и PERMUTE возможность вычисления p-значений 10000 перестановок путем перестановки столбцов матрицы данных экспрессии генов MDData и MglioData [3].
pvaluesCorr = mattest(MDData, MglioData, 'Permute', 10000);
Определите количество генов, которые считаются имеющими статистическую значимость при отсечении p-значения 0,05. Примечание: Вы можете получить другое количество генов из-за результата теста перестановки.
cutoff = 0.05; sum(pvaluesCorr < cutoff)
ans =
2121
Оценка значений FDR и q для каждого теста с использованием mafdr. Величина pi0 является общей долей истинных нулевых гипотез в исследовании. Он оценивается по моделируемому нулевому распределению с помощью начальной загрузки или кубического полинома. Примечание: Вы также можете вручную установить значение лямбда для оценки pi0.
figure;
[pFDR, qvalues] = mafdr(pvaluesCorr, 'showplot', true);
Определите количество генов, которые имеют q-значения, меньшие, чем пороговое значение. Примечание: Вы можете получить другое количество генов из-за теста перестановки и результатов начальной загрузки.
sum(qvalues < cutoff)
ans =
2173
Многие гены с низким FDR предполагают, что две группы, MD и Mglio, биологически различны.
Вы также можете эмпирически оценить скорректированные значения p FDR с помощью процедуры Бенджамини-Хохберга (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 с помощью colnames метод объекта DataMatrix.
testResults = colnames(testResults, 5, {'FDR_BH'});
Можно сортировать по значениям p pvaluesCorr с использованием sortrows матод.
testResults = sortrows(testResults, 2);
Отображение первых 20 генов в testResults. Примечание.Ваши результаты могут отличаться от приведенных ниже из-за теста перестановки и результатов начальной загрузки.
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
Вы можете использовать информацию о генной онтологии (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] Стори, Дж. Д. и Тибширани, Р., «Статистическая значимость для исследований по всему геному», PNAS, 100 (16): 9440-5, 2003.
[3] Dudoit, S., Shaffer, J.P. и Boldrick, J.C., «Множественное тестирование гипотез в эксперименте с микрочипами», Statistical Science, 18 (1): 71-103, 2003.
[4] Benjamini, Y. и Hochberg, Y., «Контроль частоты ложных открытий: практический и мощный подход к множественному тестированию», Journal of the Royal Statistical Society, Series B, 57 (1): 289-300, 1995.