Этот пример показывает ряд способов поиска шаблонов в профилях экспрессии генов.
Этот пример использует данные микромассивы экспрессии генов в дрожжах, опубликованные DeRisi и др. 1997 [1]. Авторы использовали ДНК микромассивов для изучения экспрессии височных генов почти всех генов Saccharomyces cerevisiae во время метаболического сдвига от ферментации к дыханию. Уровни экспрессии измеряли в семи временных точках во время диаксического сдвига. Полный набор данных можно загрузить с сайта Gene Expression Omnibus, http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE28.
Файл MAT yeastdata.mat содержит значения выражения (лог2 отношения CH2DN_MEAN и CH1DN_MEAN) из семи временных шагов эксперимента, имен генов и массива раз, в которые были измерены уровни экспрессии.
load yeastdata.mat
Чтобы получить представление о размере данных, которые вы можете использовать numel(genes) чтобы показать, сколько генов включено в набор данных.
numel(genes)
ans =
6400
Вы можете получить доступ к именам генов, сопоставленным с экспериментом, путем индексации переменной genes, массив ячеек, представляющая имена генов. Для примера - 15-й элемент в genes является YAL054C. Это указывает, что 15-я строка переменной yeastvalues содержит уровни выражений для YAL054C.
genes{15}
ans =
'YAL054C'
Простой график может использоваться, чтобы показать профиль выражения для этого ORF.
plot(times, yeastvalues(15,:)) xlabel('Time (Hours)'); ylabel('Log2 Relative Expression Level');
Можно также построить график фактических коэффициентов выражений, а не значений, преобразованных log2.
plot(times, 2.^yeastvalues(15,:)) xlabel('Time (Hours)'); ylabel('Relative Expression Level');
Ген, связанный с этим ORF, ACS1, по-видимому, сильно повышается во время диаксического сдвига. Можно сравнить экспрессию этого гена с экспрессией других генов путем построения нескольких линий на одном и том же рисунке.
hold on plot(times, 2.^yeastvalues(16:26,:)') xlabel('Time (Hours)'); ylabel('Relative Expression Level'); title('Profile Expression Levels');
Обычно набор данных экспрессии генов включает информацию, соответствующую генам, которые не показывают каких-либо интересных изменений во время эксперимента. Чтобы облегчить поиск интересных генов, можно уменьшить размер набора данных до некоторого подмножества, которое содержит только наиболее значимые гены.
Если вы просмотрите список генов, то увидите несколько пятен, отмеченных как 'EMPTY'. Это пустые пятна на массиве, и хотя с ними могут быть связаны данные, для целей этого примера можно считать эти точки шумом. Эти точки можно найти с помощью strcmp и удалите из набора данных с помощью команд индексации.
emptySpots = strcmp('EMPTY',genes);
yeastvalues(emptySpots,:) = [];
genes(emptySpots) = [];
numel(genes)
ans =
6314
Существует также несколько мест в наборе данных, где уровень выражения отмечен как NaN. Это указывает, что данные для этого места не были собраны на конкретном временном шаге. Одним из подходов к решению этих отсутствующих значений было бы вменение их с помощью среднего или медианы данных для конкретного гена с течением времени. Этот пример использует менее строгий подход простого выбрасывания данных для любых генов, где не был измерен один или несколько уровней экспрессии. Функция isnan используется для идентификации генов с отсутствующими данными, а команды индексации используются для удаления генов с отсутствующими данными.
nanIndices = any(isnan(yeastvalues),2); yeastvalues(nanIndices,:) = []; genes(nanIndices) = []; numel(genes)
ans =
6276
Если бы вы построили графики выражений всех остальных профилей, то увидели бы, что большинство профилей плоские и не значительно отличаются от других. Эти плоские данные, очевидно, полезны, поскольку это указывает на то, что гены, сопоставленные с этими профилями, не подвергаются значительному влиянию диаксического сдвига; однако в этом примере вас интересуют гены с большими изменениями в экспрессии, сопровождающими диаксический сдвиг. Можно использовать фильтрующие функции в Bioinformatics Toolbox™ для удаления генов с различными типами профилей, которые не дают полезной информации о генах, пораженных метаболическим изменением.
Можно использовать genevarfilter функция для фильтрации генов с небольшим отклонением с течением времени. Функция возвращает логический массив (т.е. маску) того же размера, что и переменная genes с соответствующими строкам yeastvalues с отклонением, большей, чем 10-й процентиль, и нулями, соответствующими таковым ниже порога. Можно использовать маску для индекса значений и удаления отфильтрованных генов.
mask = genevarfilter(yeastvalues); 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
Наконец, можно использовать функцию geneentropyfilter удалить гены, профили которых имеют низкую энтропию, для примера уровней энтропии в 15-м процентиле данных.
[mask,yeastvalues,genes] = geneentropyfilter(yeastvalues,genes,'prctile',15);
numel(genes)
ans = 614
Теперь, когда у вас есть управляемый список генов, можно искать отношения между профилями с помощью некоторых различных методов кластеризации из Statistics and Machine Learning Toolbox™. Для иерархической кластеризации функция pdist вычисляет парные расстояния между профилями и linkage создает иерархическое дерево кластеров.
corrDist = pdist(yeastvalues,'corr'); clusterTree = linkage(corrDist,'average');
The cluster функция вычисляет кластеры на основе расстояния отключения или максимального количества кластеров. В этом случае maxclust опция используется для идентификации 16 различных кластеров.
clusters = cluster(clusterTree,'maxclust',16);
Профили генов в этих кластерах могут быть построены вместе с помощью простого цикла и subplot команда.
figure for c = 1:16 subplot(4,4,c); plot(times,yeastvalues((clusters == c),:)'); axis tight end sgtitle('Hierarchical Clustering of Profiles');
Statistics and Machine Learning Toolbox также имеет функцию K-средних значений кластеризации. Снова найдено шестнадцать кластеров, но поскольку алгоритм отличается, они не обязательно будут такими же кластерами, как и те, которые найдены при иерархической кластеризации.
Инициализируйте состояние генератора случайных чисел, чтобы убедиться, что рисунки, сгенерированная этими командами, совпадают с рисунками в HTML версии этого примера.
rng('default'); [cidx, ctrs] = kmeans(yeastvalues,16,'dist','corr','rep',5,'disp','final'); figure for c = 1:16 subplot(4,4,c); plot(times,yeastvalues((cidx == c),:)'); axis tight end sgtitle('K-Means Clustering of Profiles');
Replicate 1, 21 iterations, total sum of distances = 23.4699. Replicate 2, 22 iterations, total sum of distances = 23.5615. Replicate 3, 10 iterations, total sum of distances = 24.823. Replicate 4, 28 iterations, total sum of distances = 23.4501. Replicate 5, 19 iterations, total sum of distances = 23.5109. Best total sum of distances = 23.4501
Вместо графического изображения всех профилей можно построить график как раз центроидов.
figure for c = 1:16 subplot(4,4,c); plot(times,ctrs(c,:)'); axis tight axis off end sgtitle('K-Means Clustering of Profiles');
Можно использовать clustergram функция для создания тепловой карты уровней выражения и дендрограммы из выхода иерархической кластеризации.
cgObj = clustergram(yeastvalues(:,2:end),'RowLabels',genes,'ColumnLabels',times(2:end));
Анализ основных компонентов (PCA) является полезным методом, который может использоваться, чтобы уменьшить размерность больших наборов данных, таких как наборы из микромассивов. PCA также может использоваться, чтобы найти сигналы в зашумленных данных. Функция mapcaplot вычисляет основные компоненты набора данных и создает графики поля точек результатов. Можно в интерактивном режиме выбрать точки данных на одном из графиков, и эти точки автоматически подсвечиваются на другом графике. Это позволяет вам визуализировать несколько размерностей одновременно.
h = mapcaplot(yeastvalues,genes);
Заметьте, что графики поля точек счетов первых двух главных компонентов показывают, что существует две отдельные области. Это не является неожиданным, поскольку процесс фильтрации удалил многие гены с низким отклонением или низкой информацией. Эти гены появились бы в середине графика поля точек.
Если необходимо просмотреть значения основных компонентов, pca функция в наборе инструментов Statistics and Machine Learning Toolbox используется для вычисления основных компонентов набора данных.
[pc, zscores, pcvars] = pca(yeastvalues);
Первый выход, pc, является матрицей основных компонентов yeastvalues данные. Первый столбец матрицы является первым основным компонентом, второй - вторым основным компонентом и так далее. Второй выход, zscores, состоит из основных счетов компонентов, то есть представления йаствалов в основном пространстве компонентов. Третий выход, pcvars, содержит отклонения основного компонента, которые дают меру того, какую часть дисперсии данных учитывает каждый из основных компонентов.
Понятно, что первый основной компонент составляет большую часть отклонения в модели. Можно вычислить точный процент отклонений, учитываемых каждым компонентом, как показано ниже.
pcvars./sum(pcvars) * 100
ans =
79.8316
9.5858
4.0781
2.6486
2.1723
0.9747
0.7089
Это означает, что почти 90% отклонения учитывается первыми двумя основными компонентами. Можно использовать cumsum команда для просмотра совокупной суммы отклонений.
cumsum(pcvars./sum(pcvars) * 100)
ans = 79.8316 89.4174 93.4955 96.1441 98.3164 99.2911 100.0000
Если вы хотите иметь больше контроля над графическим изображением основных компонентов, можно использовать scatter функция.
figure scatter(zscores(:,1),zscores(:,2)); xlabel('First Principal Component'); ylabel('Second Principal Component'); title('Principal Component Scatter Plot');
Альтернативный способ создать график поля точек - с функцией gscatter из набора инструментов Statistics and Machine Learning Toolbox. gscatter создает сгруппированный график поля точек, в которой точки из каждой группы имеют разный цвет или маркер. Можно использовать clusterdataили любая другая функция кластеризации, чтобы сгруппировать точки.
figure pcclusters = clusterdata(zscores(:,1:2),'maxclust',8,'linkage','av'); gscatter(zscores(:,1),zscores(:,2),pcclusters) xlabel('First Principal Component'); ylabel('Second Principal Component'); title('Principal Component Scatter Plot with Colored Clusters');
Если у вас есть Deep Learning Toolbox™, можно использовать самоорганизующуюся карту (SOM) для кластеризации данных.
% Check to see if the Deep Learning Toolbox is installed if ~exist(which('selforgmap'),'file') disp('The Self-Organizing Maps section of this example requires the Deep Learning Toolbox.') return end
The selforgmap функция создает новый сетевой объект SOM. В этом примере сгенерируется SOM с использованием первых двух основных компонентов.
P = zscores(:,1:2)'; net = selforgmap([4 4]);
Обучите сеть с помощью параметров по умолчанию.
net = train(net,P);
Использование plotsom отображение сети на графике поля точек данных. Обратите внимание, что алгоритм SOM использует случайные начальные точки, поэтому результаты будут варьироваться от запуска до запуска.
figure plot(P(1,:),P(2,:),'.g','markersize',20) hold on plotsom(net.iw{1,1},net.layers{1}.distances) hold off
Можно назначить кластеры с помощью SOM, найдя ближайший узел для каждой точки в наборе данных.
distances = dist(P',net.IW{1}');
[d,cndx] = min(distances,[],2); % cndx contains the cluster index
figure
gscatter(P(1,:),P(2,:),cndx); legend off;
hold on
plotsom(net.iw{1,1},net.layers{1}.distances);
hold off
Закройте все рисунки и приложения.
close('all'); delete(cgObj); delete(h); nntraintool('close');
[1] DeRisi, J.L., Iyer, V.R. and Brown, P.O., «Исследование метаболического и генетического контроля экспрессии генов на геномной шкале», Science, 278 (5338): 680-6, 1997.