Задача: Анализ выражений генов в пекарских дрожжах (Saccharomyces Cerevisiae)

Цель состоит в том, чтобы получить некоторое понимание выражений генов в Saccharomyces cerevisiae, который обычно известен как пекарские дрожжи или пивные дрожжи. Именно грибок используется для выпечки хлеба и брожения вина из винограда.

Saccharomyces cerevisiae, при введении в среду, богатую глюкозой, может превращать глюкозу в этанол. Первоначально дрожжи превращают глюкозу в этанол метаболическим процессом, называемым «ферментацией». Однако после истощения глюкозы дрожжи переходят от анаэробной ферментации глюкозы к аэробному дыханию этанола. Этот процесс называется диаксическим сдвигом. Этот процесс представляет значительный интерес, поскольку он сопровождается серьезными изменениями в экспрессии генов.

Пример использует данные микромассивов ДНК для изучения экспрессии височных генов почти всех генов в Saccharomyces cerevisiae во время диаксического сдвига.

Вам нужно Bioinformatics Toolbox™, чтобы запустить этот пример.

if ~nnDependency.bioInfoAvailable
    errordlg('This example requires Bioinformatics Toolbox.');
    return;
end

Данные

Этот пример использует данные из DeRisi, JL, Iyer, VR, Brown, PO. «Исследование метаболического и генетического контроля экспрессии генов в геномной шкале». Наука. 1997 Окт 24; 278 (5338): 680-6. PMID: 9381177

Полный набор данных можно скачать с сайта Gene Expression Omnibus: https://www.yeastgenome.org

Начните, загрузив данные в MATLAB ®.

load yeastdata.mat

Уровни экспрессии генов измеряли в семи временных точках во время диаксического сдвига. Переменная times содержит время, в которое уровни экспрессии были измерены в эксперименте. Переменная genes содержит имена генов, уровни экспрессии которых были измерены. Переменная yeastvalues содержит данные «VALUE» или LOG_RAT2N_MEAN, или лог2 отношения CH2DN_MEAN и CH1DN_MEAN от семи временных шагов в эксперименте.

Чтобы получить представление о размере данных, вы можете использовать numel (гены), чтобы показать, сколько генов в наборе данных.

numel(genes)

ans =

        6400

гены являются массивом ячеек генных имен. Вы можете получить доступ к записям с помощью индексации массива ячеек MATLAB:

genes{15}

ans =

    'YAL054C'

Это указывает, что 15-я строка переменных yeastvalues содержит уровни экспрессии для ORF YAL054C.

Фильтрация генов

Набор данных довольно велик и много информации соответствует генам, которые не показывают каких-либо интересных изменений в ходе эксперимента. Чтобы облегчить поиск интересных генов, первое, что нужно сделать, это уменьшить размер набора данных путем удаления генов с профилями экспрессии, которые не показывают ничего интересного. Существует 6400 профилей выражений. Можно использовать ряд методов, чтобы свести это к некоторому подмножеству, которое содержит наиболее значимые гены.

Если вы просмотрите список генов, то увидите несколько пятен, отмеченных как 'EMPTY'. Это пустые пятна на массиве, и хотя с ними могут быть связаны данные, для целей этого примера можно считать эти точки шумом. Эти точки можно найти с помощью функции strcmp и удалить из набора данных с помощью команд индексации.

emptySpots = strcmp('EMPTY',genes);
yeastvalues(emptySpots,:) = [];
genes(emptySpots) = [];
numel(genes)

ans =

        6314

В данных yeastvalues вы также увидите несколько мест, где уровень выражения отмечен как NaN. Это указывает, что данные для этого места не были собраны на конкретном временном шаге. Одним из подходов к решению этих отсутствующих значений было бы вменение их с помощью среднего или медианы данных для конкретного гена с течением времени. Этот пример использует менее строгий подход простого выбрасывания данных для любых генов, где не был измерен один или несколько уровней экспрессии.

Функция isnan используется для идентификации генов с отсутствующими данными, а команды индексации используются для удаления генов с отсутствующими данными.

nanIndices = any(isnan(yeastvalues),2);
yeastvalues(nanIndices,:) = [];
genes(nanIndices) = [];
numel(genes)

ans =

        6276

Если бы вы построили графики выражений всех остальных профилей, то увидели бы, что большинство профилей плоские и не значительно отличаются от других. Эти плоские данные, очевидно, полезны, поскольку это указывает на то, что гены, сопоставленные с этими профилями, не подвергаются значительному влиянию диаксического сдвига; однако в этом примере вас интересуют гены с большими изменениями в экспрессии, сопровождающими диаксический сдвиг. Можно использовать фильтрующие функции в Bioinformatics Toolbox™ для удаления генов с различными типами профилей, которые не дают полезной информации о генах, пораженных метаболическим изменением.

Можно использовать функцию genevarfilter, чтобы отфильтровать гены с небольшим отклонением с течением времени. Функция возвращает логический массив того же размера, что и гены переменных с теми, которые соответствуют строкам дрожжевых значений с отклонением, большей, чем 10-й процентиль, и нулями, соответствующими тем, которые ниже порога.

mask = genevarfilter(yeastvalues);
% Use the mask as an index into the values to remove the filtered genes.
yeastvalues = yeastvalues(mask,:);
genes = genes(mask);
numel(genes)

ans =

        5648

Функция genelowvalfilter удаляет гены, которые имеют очень низкие абсолютные значения экспрессии. Обратите внимание, что функции filter генов также могут автоматически вычислять отфильтрованные данные и имена.

[mask, yeastvalues, genes] = ...
    genelowvalfilter(yeastvalues,genes,'absval',log2(3));
numel(genes)

ans =

   822

Используйте генеентропифильтер для удаления генов, профили которых имеют низкую энтропию:

[mask, yeastvalues, genes] = ...
    geneentropyfilter(yeastvalues,genes,'prctile',15);
numel(genes)

ans =

   614

Анализ основных компонентов

Теперь, когда у вас есть управляемый список генов, вы можете искать отношения между профилями.

Нормализация стандартного отклонения и среднего значения данных позволяет сети рассматривать каждый вход как одинаково важный для своей области значений значений.

Анализ основных компонентов (PCA) является полезным методом, который может использоваться, чтобы уменьшить размерность больших наборов данных, таких как наборы данных из микромассива анализа. Этот метод изолирует основные компоненты набора данных, устраняя те компоненты, которые вносят наименьший вклад в изменение набора данных.

Две переменные настройки могут использоваться, чтобы применить mapstd и processpca к новым данным, которые будут согласованными.

[x,std_settings] = mapstd(yeastvalues');  % Normalize data
[x,pca_settings] = processpca(x,0.15);    % PCA

Входные векторы сначала нормированы, используя mapstd, так что они имеют нулевое среднее и отклонение единства. processpca - функция, которая реализует алгоритм PCA. Второй аргумент передан processpca составляет 0,15. Это означает, что processpca устраняет основные компоненты, которые вносят менее 15% в общее изменение набора данных. Переменная pc теперь содержит основные компоненты данных о дрожжевых значениях.

Основные компоненты могут быть визуализированы с помощью функции рассеяния.

figure
scatter(x(1,:),x(2,:));
xlabel('First Principal Component');
ylabel('Second Principal Component');
title('Principal Component Scatter Plot');

Кластерный анализ: самоорганизующиеся карты

Теперь основные компоненты можно кластеризировать с помощью алгоритма кластеризации SOM.

Функция selforgmap создает сеть самоорганизующихся карт, которая затем может быть обучена с функцией train.

Размер входа 0, потому что сеть еще не сконфигурирована, чтобы соответствовать нашим входным данным. Это произойдет, когда сеть будет обучена.

net = selforgmap([5 3]);
view(net)

Сейчас сеть готова к обучению.

Neural Network Training Tool показывает обучаемую сеть и алгоритмы, используемые для ее обучения. В нем также отображается состояние обучения во время обучения, и критерии, которые остановили обучение, будут выделены зеленым цветом.

Кнопки внизу открывают полезные графики, которые можно открыть во время и после обучения. Ссылки рядом с именами алгоритмов и кнопками графика открывают документацию по этим темам.

net = train(net,x);
nntraintool
nntraintool('close')

Используйте plotsompos, чтобы отобразить сеть на графике поля точек первых двух размерностей данных.

figure
plotsompos(net,x);

Можно назначить кластеры с помощью SOM, найдя ближайший узел для каждой точки в наборе данных.

y = net(x);
cluster_indices = vec2ind(y);

Используйте plotsomhits, чтобы увидеть, сколько векторов назначено каждому из нейронов на карте.

figure
plotsomhits(net,x);

Можно также использовать другие алгоритмы кластеризации, такие как Иерархическая кластеризация и K-средних значений, доступные в Statistics and Machine Learning Toolbox™ для кластерного анализа.

Глоссарий

ORF - Открытая система координат считывания (ORF) является фрагментом последовательности гена, которая содержит последовательность основ, непрерывно прерываемых стоповыми последовательностями, которые потенциально могут кодировать белок.