Этот пример показывает ряд способов поиска паттернов в профилях экспрессии генов.
В этом примере используются данные исследования микрочипов экспрессии генов в дрожжах, опубликованные DeRisi et al. 1997 [1]. Авторы использовали ДНК-микрочипы для изучения временной экспрессии генов почти всех генов в Saccharomyces cerevisiae во время метаболического сдвига от ферментации к дыханию. Уровни экспрессии измеряли в семь временных точек во время диауксического сдвига. Полный набор данных можно загрузить с веб-сайта Gene Expression Omnibus, http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE28.
Файл MAT yeastdata.mat содержит значения выражения (log2 отношения 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
При построении графика профилей выражений для всех оставшихся профилей будет видно, что большинство профилей являются плоскими и существенно не отличаются от остальных. Эти плоские данные, очевидно, полезны, поскольку они указывают на то, что гены, связанные с этими профилями, не подвергаются существенному влиянию диауксического сдвига; однако в этом примере вас интересуют гены с большими изменениями экспрессии, сопровождающими диауксический сдвиг. Вы можете использовать функции фильтрации в Toolbox™ биоинформатики, чтобы удалить гены с различными типами профилей, которые не предоставляют полезную информацию о генах, затронутых метаболическим изменением.
Вы можете использовать genevarfilter функция фильтрации генов с небольшой дисперсией во времени. Функция возвращает логический массив (т.е. маску) того же размера, что и переменная genes с соответствующими строкам yeastvalues с дисперсией больше 10-го процентиля и нулями, соответствующими нулям ниже порогового значения. Маску можно использовать для индексирования значений и удаления отфильтрованных генов.
mask = genevarfilter(yeastvalues); 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 удалить гены, профили которых имеют низкую энтропию, например уровни энтропии в 15-м процентиле данных.
[mask,yeastvalues,genes] = geneentropyfilter(yeastvalues,genes,'prctile',15);
numel(genes)
ans = 614
Теперь, когда у вас есть управляемый список генов, вы можете искать связи между профилями, используя некоторые другие методы кластеризации из Toolbox™ статистики и машинного обучения. Для иерархической кластеризации функция pdist вычисляет попарные расстояния между профилями и linkage создает иерархическое дерево кластера.
corrDist = pdist(yeastvalues,'corr'); clusterTree = linkage(corrDist,'average');
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');
Инструментарий статистики и машинного обучения также имеет функцию кластеризации K-means. Опять же, шестнадцать кластеров найдены, но поскольку алгоритм отличается, они не обязательно будут теми же кластерами, которые найдены иерархической кластеризацией.
Инициализируйте состояние генератора случайных чисел, чтобы цифры, сгенерированные этой командой, соответствовали цифрам в 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 функция на панели инструментов статистики и машинного обучения используется для вычисления основных компонентов набора данных.
[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 на панели инструментов статистики и машинного обучения. 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');
При наличии Toolbox™ Deep Learning можно использовать самоорганизующуюся карту (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
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] ДеРиси, Дж. Л., Айер, В.Р. и Браун, П.О., «Исследование метаболического и генетического контроля экспрессии генов в геномном масштабе», Science, 278 (5338): 680-6, 1997.