Этот пример показывает много способов искать шаблоны в профилях экспрессии гена.
Этот пример использует данные из исследования микромассивов экспрессии гена в дрожжах, опубликованных DeRisi, и др. 1997 [1]. Авторы использовали микромассивы ДНК, чтобы изучить временную экспрессию гена почти всех генов в Saccharomyces cerevisiae во время метаболического сдвига от ферментации до дыхания. Уровни экспрессии были измерены в семи моментах времени во время сдвига diauxic. Полный набор данных может быть загружен с веб-сайта Gene Expression Omnibus, https://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, кажется, строго отрегулирован во время сдвига diauxic. Можно сравнить выражение этого гена к выражению других генов путем строения нескольких графиков на той же фигуре.
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
Если бы необходимо было построить профили выражения всех остающихся профилей, вы видели бы, что большинство профилей является плоским и не существенно отличается от других. Эти плоские данные, очевидно, имеют применение, когда это указывает, что гены, сопоставленные с этими профилями, не значительно затронуты сдвигом diauxic; однако, в этом примере, вы интересуетесь генами с большими изменениями в выражении, сопровождающем сдвиг diauxic. Можно использовать функции фильтрации в 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');
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, состоит из баллов основного компонента, т.е. представления yeastvalues на пробеле основного компонента. Третий выход, 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
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);
[1] DeRisi, J.L., Iyer, В.Р. и Браун, отделение связи, "Исследуя метаболическое и генетическое управление экспрессии гена по геномной шкале", Наука, 278 (5338):680-6, 1997.