Исследование данных об экспрессии гена микромассивов
В этом примере показано, как идентифицировать дифференцированно выраженные гены из данных микромассивов и Генной Онтологии использования, чтобы определить значительные биологические функции, которые сопоставлены к вниз - и отрегулированные гены.
Введение
Микромассивы содержат олигонуклеотид или зонды комплементарной ДНК для сравнения профиля выражения генов по геномной шкале. Определение, если изменения в экспрессии гена являются статистически значительными между различными условиями, например, двумя различными типами опухоли и определением биологической функции дифференцированно выраженных генов, является важными целями в эксперименте микромассивов.
Общедоступный набор данных, содержащий данные об экспрессии гена 42 тканей опухоли эмбриональной центральной нервной системы (CNS) [1], используется в этом примере. Файлы CEL могут быть загружены здесь. Выборки были гибридизированы на массивах Affymetrix® HuGeneFL GeneChip®. Необработанный набор данных был предварительно обработан с Устойчивым средним значением мультимассивов (RMA) и GC Устойчивое Среднее значение Мультимассивов (GCRMA) процедуры. Для получения дополнительной информации о предварительной обработке олигонуклеотида микромассивов Affymetrix смотрите Предварительную обработку Affymetrix® Microarray Data на Тестовом Уровне.
Вы будете использовать 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
Идентификация дифференциальной экспрессии гена
Можно теперь сравнить значения экспрессии гена между двумя группами данных: CNS medulloblastomas (MD) и ненейронный источник злокачественные глиомы (Mglio) опухоль.
Из данных о выражении всех 42 выборок в наборе данных извлеките данные 10 выборок миллидня и 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].
В этом примере вы вычислите ФРГ с помощью процедуры [2] Яруса-Tibshirani. Процедура также вычисляет q-значение теста, который измеряет минимальный ФРГ, который происходит при вызове значительного теста. Оценка ФРГ зависит от действительно пустого распределения нескольких тестов, которое неизвестно. Методы сочетания могут использоваться, чтобы оценить действительно пустое распределение тестовой статистики путем перестановки столбцов матрицы данных экспрессии гена [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
Оцените ФРГ и q-значения для каждого теста с помощью mafdr. Количество pi0 является полной пропорцией истинных нулевых гипотез в исследовании. Это оценивается от симулированного пустого распределения через начальную загрузку или подгонку кубического полинома. Примечание: можно также вручную установить значение lambda для оценки pi0.
figure;
[pFDR, qvalues] = mafdr(pvaluesCorr, 'showplot', true);
Определите количество генов, которые имеют q-значения меньше, чем значение сокращения. Примечание: можно получить различное количество генов из-за теста сочетания и результатов начальной загрузки.
sum(qvalues < cutoff)
ans =
2173
Много генов с низким ФРГ подразумевают, что эти две группы, MD и Mglio, биологически отличны.
Можно также опытным путем оценить, что ФРГ настроил p-значения с помощью процедуры [4] Benjamini-Hochberg (BH) путем установки mafdr введите параметр BHFDR к истине.
pvaluesBH = mafdr(pvaluesCorr, 'BHFDR', true);
sum(pvaluesBH < cutoff)
ans =
1374
Можно сохранить t-баллы, p-значения, pFDRs, q-значения и ФРГ BH откорректировали p-значения вместе как объект DataMatrix.
testResults = [tscores pvaluesCorr pFDR qvalues pvaluesBH];
Обновитесь имя столбца для ФРГ BH откорректировало p-значения с помощью 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-значений против биологического эффекта в графике вулкана. Примечание: Из графика вулкана пользовательский интерфейс, можно в интерактивном режиме изменить сокращение 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, найдены в отрегулированном списке генов, в то время как гены, типичные для астроцитарного и oligodendrocytic происхождения и клеточной дифференцировки, такие как 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), чтобы аннотировать дифференцированно выраженные гены, идентифицированные выше. Пойти файл аннотации (GAF) для человека может быть загружен с Генной Онтологии Текущие Аннотации. Для удобства карта между названиями генов и сопоставленный ПЕРЕХОДИТ идентификаторы относительно к полю Function аспекта включен в файл MAT goa_human.
load goa_human
В качестве альтернативы можно запустить код ниже, чтобы загрузить базу данных Gene Ontology последними аннотациями, читать, загруженный файл аннотации Человека разумного (назовите файл как goa_human.gaf), и создайте отображение между названиями генов, и связанные ИДУТ условия.
% 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 аннотируются. Для каждого гена на чипе смотрите, аннотируется ли это путем сравнения его названия гена со списком названий генов от, ИДУТ. Отследите количество аннотируемых генов, и количество отрегулированных генов, сопоставленных с каждым, ИДУТ термин. Обратите внимание на то, что данные в общедоступных репозиториях часто курируются и обновляются; поэтому результаты этого примера могут немного отличаться, когда вы используете актуальные наборы данных. Также возможно, что вы получаете предупреждения о недопустимых или устаревших идентификаторах из-за обновленного генного файла аннотации Человека разумного.
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
Определите статистически значительный, ИДУТ условия с помощью гипергеометрического вероятностного распределения. Поскольку каждый ИДЕТ термин, p-значение вычисляется, представляя вероятность, что количество аннотируемых генов, сопоставленных с ним, возможно, было найдено случайно.
gopvalues = hygepdf(upgenesCount,max(chipgenesCount),...
max(upgenesCount),chipgenesCount);
[dummy, idx] = sort(gopvalues);
Сообщите, что лучшие старшие значащие десять ИДУТ условия можно следующим образом.
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 функция. Можно окрасить вершины графиков согласно их значению. В этом примере красные узлы старше значащие, в то время как синие узлы являются младшими значащими генными условиями онтологии. Примечание: ПОЙТИ возвращенные условия могут отличаться от показанных из-за частого обновления генного файла аннотации Человека разумного.
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., и др., "Прогноз центральной нервной системы эмбриональный результат опухоли на основе экспрессии гена". Природа, 415 (6870):436-42, 2001.
[2] Ярус, степень доктора юридических наук и Tibshirani, R., "Статистическое значение для genomewide исследований", PNAS, 100 (16):9440-5, 2003.
[3] Dudoit, S., Шаффер, J.P., и Boldrick, J.C., "Несколько тестирование гипотезы в эксперименте микромассивов", Статистическая Наука, 18 (1):71-103, 2003.
[4] Benjamini, Y. и Hochberg, Y., "Управляя ложным уровнем открытия: практический и мощный подход к нескольким тестирование", Журнал Королевского Статистического Общества, Серий B, 57 (1):289-300, 1995.