Введение

Микромассивы содержат олигонуклеотид или зонды комплементарной ДНК для сравнения профиля выражения генов по геномной шкале. Определяя, являются ли изменения в экспрессии гена статистически значительными между различными условиями, e.g. два различных типа опухоли и определение биологической функции дифференцированно описанных генов, являются важными целями в эксперименте микромассивов.

Общедоступный набор данных, содержащий данные об экспрессии гена 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.

Следующее является несколькими картами трассы с генами в отрегулированном списке генов, подсветил:

Клеточный цикл

Еж Сигнальная трасса

mTor Сигнальная трасса

Ссылки

[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.