Исследование данных об экспрессии гена микромассивов
Этот пример показывает, как идентифицировать дифференцированно выраженные гены от микроданных массива и использует Генную Онтологию, чтобы определить значительные биологические функции, которые сопоставлены к вниз - и отрегулированные гены.
Введение
Микромассивы содержат олигонуклеотид или зонды комплементарной ДНК для сравнения профиля выражения генов в геномной шкале. Определение, если изменения в экспрессии гена являются статистически значительными между различными условиями, например, двумя различными типами опухоли и определением биологической функции дифференцированно выраженных генов, является важными целями в эксперименте микромассивов.
Общедоступный набор данных, содержащий данные об экспрессии гена 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 каждая строка соответствует тестовому набору на массиве, и каждый столбец соответствует выборке. expr_cns_gcrma_eb объекта DataMatrix будет использоваться в этом примере, но данные из любой из других двух переменных выражения могут использоваться также.
Получите свойства объекта DataMatrix expr_cns_gcrma_eb с помощью команды get.
get(expr_cns_gcrma_eb)
Name: ''
RowNames: {7129x1 cell}
ColNames: {1x42 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: {5669x1 cell}
ColNames: {1x10 cell}
NRows: 5669
NCols: 10
NDims: 2
ElementClass: 'single'
Name: ''
RowNames: {5669x1 cell}
ColNames: {1x10 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 =
2167
Много генов с низким ФРГ подразумевают, что эти две группы, MD и Mglio, биологически отличны.
Можно также опытным путем оценить, что ФРГ настроил p-значения с помощью процедуры [4] Benjamini-Hochberg (BH) путем установки параметра входа mafdr BHFDR на истину.
pvaluesBH = mafdr(pvaluesCorr, 'BHFDR', true);
sum(pvaluesBH < cutoff)
ans =
1373
Можно сохранить t-очки, p-значения, pFDRs, q-значения и ФРГ BH исправили p-значения вместе как объект 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
RAB31 -13.664 5.1444e-11 1.0471e-07 1.0471e-07 2.9164e-07
PLEC1 -9.6223 1.0596e-07 0.00010784 8.8536e-05 0.00024658
HNRPA1 9.359 2.0854e-07 0.0001415 8.8536e-05 0.00024658
FCGR2A -9.3548 2.2318e-07 0.00011357 8.8536e-05 0.00024658
PLEC1 -9.3495 2.6476e-07 0.00010778 8.8536e-05 0.00024658
FBL 9.1518 3.1671e-07 0.00010744 8.8536e-05 0.00024658
KIAA0367 -8.996 3.8835e-07 0.00011293 8.8536e-05 0.00024658
ID2B -8.9285 4.0623e-07 0.00010336 8.8536e-05 0.00024658
RBMX 8.8905 4.2322e-07 9.5717e-05 8.8536e-05 0.00024658
PAFAH1B3 8.7561 4.5893e-07 9.3415e-05 8.8536e-05 0.00024658
H3F3A 8.6512 5.0812e-07 9.4024e-05 8.8536e-05 0.00024658
LRP1 -8.6465 5.2195e-07 8.8536e-05 8.8536e-05 0.00024658
PEA15 -8.3256 7.5427e-07 0.0001181 0.0001149 0.00032
ID2B -8.1183 8.4551e-07 0.00012293 0.0001149 0.00032
SFRS3 8.1166 8.4672e-07 0.0001149 0.0001149 0.00032
HLA-DPA1 -7.8546 1.0944e-06 0.00013923 0.00013717 0.00038204
C5orf13 7.7195 1.1781e-06 0.00014106 0.00013717 0.00038204
PTMA 7.7013 1.213e-06 0.00013717 0.00013717 0.00038204
NAP1L1 7.674 1.3307e-06 0.00014255 0.000138 0.00038434
HMGB2 7.6532 1.3559e-06 0.000138 0.000138 0.00038434
Ген считается дифференцированно выраженным между двумя группами выборок, если он показывает и статистическое и биологическое значение. Этот пример сравнивает отношение экспрессии гена MD по выборкам опухоли Mglio. Поэтому отрегулированный ген в этом примере имеет более высокое выражение в MD, и вниз отрегулированный ген имеет более высокое выражение в Mglio.
Постройте-log10 p-значений против биологического эффекта в графике вулкана. Примечание: Из графика вулкана пользовательский интерфейс, можно в интерактивном режиме изменить сокращение p-значения и свернуть предел изменения и экспортировать дифференцированно выраженные гены.
diffStruct = mavolcanoplot(MDData, MglioData, pvaluesCorr)
diffStruct =
Name: 'Differentially Expressed'
PVCutoff: 0.0500
FCThreshold: 2
GeneLabels: {327x1 cell}
PValues: [327x1 bioma.data.DataMatrix]
FoldChanges: [327x1 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), чтобы аннотировать дифференцированно выраженные гены, идентифицированные выше. Файл аннотации для Человека разумного (gene_association.goa_human.gz) может быть загружен с Генной Онтологии Текущие Аннотации. Для удобства ИДЕТ карта между названиями генов и сопоставленный, идентификатор относительно к полю Function аспекта включен в файл MAT goa_human.
load goa_human
Также можно запустить код ниже, чтобы загрузить базу данных Gene Ontology последними аннотациями, считать загруженный файл аннотации Человека разумного и создать отображение между названиями генов, и связанные ИДУТ условия.
% GO = geneont('live',true); % HGann = goannotread('gene_association.goa_human',... % '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.