Этот пример демонстрирует поиск паттернов в профилях экспрессии генов в пекарских дрожжах с использованием нейронных сетей.
Цель состоит в том, чтобы получить некоторое понимание экспрессии генов в Saccharomyces cerevisiae, который обычно известен как пекарские дрожжи или пивные дрожжи. Именно гриб используется для выпечки хлеба и закваски вина из винограда.
Saccharomyces cerevisiae при введении в среду, богатую глюкозой, могут превращать глюкозу в этанол. Первоначально дрожжи превращают глюкозу в этанол метаболическим процессом, называемым «ферментацией». Однако после истощения запасов глюкозы дрожжи переходят от анаэробной ферментации глюкозы к аэробному дыханию этанолом. Этот процесс называется диауксическим сдвигом. Этот процесс представляет значительный интерес, поскольку он сопровождается серьезными изменениями в экспрессии генов.
В примере используются данные микрочипов ДНК для изучения временной экспрессии генов почти всех генов в Saccharomyces cerevisiae во время диауксического сдвига.
Для выполнения этого примера требуется 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, или log2 отношения 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
При построении графика профилей выражений для всех оставшихся профилей будет видно, что большинство профилей являются плоскими и существенно не отличаются от остальных. Эти плоские данные, очевидно, полезны, поскольку они указывают на то, что гены, связанные с этими профилями, не подвергаются существенному влиянию диауксического сдвига; однако в этом примере вас интересуют гены с большими изменениями экспрессии, сопровождающими диауксический сдвиг. Вы можете использовать функции фильтрации в 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 удаляет гены, которые имеют очень низкие абсолютные значения экспрессии. Обратите внимание, что функции генного фильтра также могут автоматически вычислять отфильтрованные данные и имена.
[mask, yeastvalues, genes] = ... genelowvalfilter(yeastvalues,genes,'absval',log2(3)); numel(genes)
ans = 822
Используйте geneentropyfilter для удаления генов, профили которых имеют низкую энтропию:
[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 создает сеть самоорганизующихся карт, которая затем может быть обучена с помощью функции поезда.
Входной размер равен 0, поскольку сеть еще не настроена на соответствие нашим входным данным. Это произойдет, когда сеть будет обучена.
net = selforgmap([5 3]); view(net)
Сейчас сеть готова к обучению.
Обучающий инструмент нейронной сети показывает обучаемую сеть и алгоритмы, используемые для ее обучения. Он также отображает состояние обучения во время обучения и критерии, которые остановили обучение, будут выделены зеленым цветом.
Кнопки внизу открывают полезные графики, которые можно открывать во время и после обучения. Ссылки рядом с именами алгоритмов и кнопками графика открывают документацию по этим темам.
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-средства, доступные в Toolbox™ Статистика и машинное обучение для кластерного анализа.
ORF - Открытая рамка считывания (ORF) - это часть последовательности гена, которая содержит последовательность оснований, непрерывную стоп-последовательностями, которая потенциально может кодировать белок.