Этот пример иллюстрирует простой подход к поиску потенциальных регулирующих мотивов в наборе коэкспрессированных геномных последовательностей путем идентификации значительно превалирующих несодержащих разрывы слов фиксированной длины. Обсуждение основано на тематическом исследовании, представленном в Главе 10 "Введения в Вычислительную Геномику. Подход Тематических исследований" [1].
Циркадные часы являются 24-часовым циклом физиологических процессов, которые синхронизируются с внешним круглосуточным циклом. Большая часть работы над циркадным генератором на объектах была выполнена с помощью объекта модели Arabidopsis thaliana. В этом организме регулирование серии генов, которые должны быть включены или выключены в определенное время дня и ночи, выполняется маленькими регулирующими последовательностями, найденными в восходящем направлении рассматриваемые гены. Один такой регулирующий мотив, AAAATATCT, также известный как Вечерний элемент (EE), был идентифицирован в областях промоутера циркадных отрегулированных часами генов, которые показывают пиковое выражение вечером [2].
Мы рассматриваем три набора отрегулированных часами генов, кластеризируемых согласно времени дня, когда они максимально описываются: установите 1, соответствует восходящим областям 1 Кбайт длиной генов, выражение которых достигают максимума утром (8:00 - 16:00); установите 2, соответствует восходящим областям 1 Кбайт длиной генов, выражение которых достигают максимума вечером (16:00 - 12:00); установите 3, соответствует восходящим областям 1 Кбайт длиной генов, выражение которых достигают максимума ночью (12:00 - 8:00). Поскольку мы интересуемся регулирующим мотивом в вечерних генах, установите 2, представляет наш целевой набор, в то время как установлено 1 и устанавливает 3, используются в качестве фона. В каждом наборе последовательности и их соответствующие противоположные дополнения конкатенированы друг другу с отдельными последовательностями, разделенными символом разрыва (-).
load evemotifdemodata.mat; % === concatenate both strands s1 = [[set1.Sequence] seqrcomplement([set1.Sequence])]; s2 = [[set2.Sequence] seqrcomplement([set2.Sequence])]; s3 = [[set3.Sequence] seqrcomplement([set3.Sequence])]; % === compute length and number of sequences in each set L1 = length(set1(1).Sequence); L2 = length(set2(1).Sequence); L3 = length(set3(1).Sequence); N1 = numel(set1) * 2; N2 = numel(set2) * 2; N3 = numel(set3) * 2; % === add separator between sequences seq1 = seqinsertgaps(s1, 1:L1:(L1*N1)+N1, 1); seq2 = seqinsertgaps(s2, 1:L2:(L2*N2)+N2, 1); seq3 = seqinsertgaps(s3, 1:L3:(L3*N3)+N3, 1);
Чтобы определить, который мотив кандидата превалирует в данном целевом наборе относительно фонового набора, мы идентифицируем весь возможный W-mers (слова длины W) в обоих наборах и вычисляем их частоту. Слово рассматривается превалирующим, если его частота в целевом наборе значительно выше, чем частота в фоновом режиме набор. Это различие также называется "полем".
type findOverrepresentedWords
function [nmersSorted, freqDiffSorted] = findOverrepresentedWords(seq, seq0, W) % FINDOVERREPRESENTEDWORDS helper for evemotifdemo % Copyright 2007 The MathWorks, Inc. %=== find and count words of length W nmers0 = nmercount(seq0, W); nmers = nmercount(seq, W); %=== compute frequency of words f = [nmers{:,2}]/(length(seq) - W + 1); f0 = [nmers0{:,2}]/(length(seq0) - W + 1); %=== determine words common to both set [nmersInt, i1, i2] = intersect(nmers(:,1),nmers0(:,1)); freqDiffInt = (f(i1) - f0(i2))'; %=== determine words specific to one set only [nmersXOr, i3, i4] = setxor(nmers(:,1),nmers0(:,1)); c0 = nmers(i3,1); d0 = nmers0(i4,1); nmersXOr = [c0; d0]; freqDiffXOr = [f(i3) -f0(i4)]'; %=== define all words and their difference in frequency (margin) nmersAll = [nmersInt; nmersXOr]; freqDiff = [freqDiffInt; freqDiffXOr]; %=== sort according to descending difference in frequency [freqDiffSorted, freqDiffSortedIndex] = sort(freqDiff, 'descend'); nmersSorted = nmersAll(freqDiffSortedIndex);
Если мы рассматриваем все слова длины W = 9, которые появляются более часто в целевом наборе (восходящая область генов, высоко описанных вечером) относительно фонового набора (восходящая область генов, высоко описанных утром и ночь), мы замечаем, что самое превалирующее слово является 'AAAATATCT', также известным как мотив Вечернего элемента (EE).
W = 9; [words, freqDiff] = findOverrepresentedWords(seq2, [seq1 seq3],W); words(1:10) freqDiff(1:10)
ans = 10x1 cell array {'AAAATATCT'} {'AGATATTTT'} {'CTCTCTCTC'} {'GAGAGAGAG'} {'AGAGAGAGA'} {'TCTCTCTCT'} {'AAATATCTT'} {'AAGATATTT'} {'AAAAATATC'} {'GATATTTTT'} ans = 1.0e-03 * 0.1439 0.1439 0.1140 0.1140 0.1074 0.1074 0.0713 0.0713 0.0695 0.0695
Помимо мотива EE, другие слова длины W = 9, кажется, превалируют в целевом наборе. В частности, мы замечаем присутствие повторений, i.e., слова, состоящие из одного нуклеотида или димера, повторенного для целого размера слова, такие как 'CTCTCTCTC'. Это явление довольно распространено в геномных последовательностях и обычно сопоставляется с нефункциональными компонентами. Поскольку в этом контексте повторения вряд ли будут биологически значительными, мы фильтруем их.
% === determine repeats wordsN = numel(words); r = zeros(wordsN,1); for i = 1:wordsN if (all(words{i}(1:2:end) == words{i}(1)) && ... % odd positions are the same all(words{i}(2:2:end) == words{i}(2))) % even positions are the same r(i) = 1; end end r = logical(r); % === filter out repeats words = words(~r); freqDiff = freqDiff(~r); % === consider the top 10 motifs motif = words(1:10) margin = freqDiff(1:10) EEMotif = motif{1} EEMargin = margin(1)
motif = 10x1 cell array {'AAAATATCT'} {'AGATATTTT'} {'AAATATCTT'} {'AAGATATTT'} {'AAAAATATC'} {'GATATTTTT'} {'AAATAAAAT'} {'ATTTTATTT'} {'TAAATAAAA'} {'TTTTATTTA'} margin = 1.0e-03 * 0.1439 0.1439 0.0713 0.0713 0.0695 0.0695 0.0656 0.0656 0.0600 0.0600 EEMotif = 'AAAATATCT' EEMargin = 1.4393e-04
После удаления повторений мы замечаем, что мотив EE ('AAAATATCT') и его противоположное дополнение ('AGATATTTT') наверху списка. Другие превалирующие слова являются или простыми вариантами мотива EE, такими как 'AAATATCTT', 'AAAAATATC', 'AAATATCTC', или их противоположными дополнениями, такими как 'AAGATATTT', 'GATATTTTT', 'GAGATATTT'.
Различные методы могут использоваться, чтобы оценить статистическое значение поля, вычисленного для мотива EE. Например, мы можем повторить анализ с помощью некоторых управляющих последовательностей и оценить получившиеся поля относительно поля EE. Геномные области Arabidopsis thaliana, которые являются еще дальше от записи, создают сайт, хорошие кандидаты с этой целью. В качестве альтернативы мы могли случайным образом разделить и переставить последовательности на рассмотрении и использовать их в качестве средств управления. Другое простое решение должно сгенерировать случайные последовательности согласно составу нуклеотида трех исходных наборов последовательностей, как показано ниже.
% === find base composition of each set bases1 = basecount(s1); bases2 = basecount(s2); bases3 = basecount(s3); % === generate random sequences according to base composition rs1 = randseq(length(s1),'fromstructure', bases1); rs2 = randseq(length(s2),'fromstructure', bases2); rs3 = randseq(length(s3),'fromstructure', bases3); % === add separator between sequences rseq1 = seqinsertgaps(rs1, 1:L1:(L1*N1)+N1, 1); rseq2 = seqinsertgaps(rs2, 1:L2:(L2*N2)+N2, 1); rseq3 = seqinsertgaps(rs3, 1:L3:(L3*N3)+N3, 1); % === compute margins for control set [words, freqDiff] = findOverrepresentedWords(rseq2, [rseq1 rseq3],W);
Переменная ctrlMargin
содержит предполагаемые поля главных мотивов для каждой из этих 100 управляющих последовательностей, сгенерированных аналогичный описанному выше. Распределение этих полей может быть аппроксимировано распределением экстремума. Мы используем функциональный gevfit
от Statistics and Machine Learning Toolbox™, чтобы оценить параметры (форма, шкала и местоположение) распределения экстремума и мы накладываем масштабированную версию его функции плотности вероятности, вычисленное использование gevpdf
, с гистограммой полей управляющих последовательностей.
% === estimate parameters of distribution nCtrl = length(ctrlMargin); buckets = ceil(nCtrl/10); parmhat = gevfit(ctrlMargin); k = parmhat(1); % shape parameter sigma = parmhat(2); % scale parameter mu = parmhat(3); % location parameter % === compute probability density function x = linspace(min(ctrlMargin), max([ctrlMargin EEMargin])); y = gevpdf(x, k, sigma, mu); % === scale probability density function [v, c] = hist(ctrlMargin,buckets); binWidth = c(2) - c(1); scaleFactor = nCtrl * binWidth; % === overlay figure() hold on; hist(ctrlMargin, buckets); h = findobj(gca,'Type','patch'); h.FaceColor = [.9 .9 .9]; plot(x, scaleFactor * y, 'r'); stem(EEMargin, 1, 'b'); xlabel('Margin'); ylabel('Number of sequences'); legend('Ctrl Margins', 'EVD pdf', 'EE Margin');
Поля управления являются различиями в частоте, которую мы ожидали бы находить, когда слово превалирует случайно один. Поле относительно мотива EE ясно значительно больше, чем поля управления и не соответствует в кривой плотности вероятности случайных средств управления. Поскольку поле EE больше, чем все 100 полей управления, мы можем прийти к заключению, что перепредставленность мотива EE в целевом наборе является статистически значительной, и оценка p-значения меньше 0.01.
Если мы повторяем поиск превалирующих слов длины W = 6... 11, мы замечаем, что все главные мотивы являются любой подстроками (если W <9) или суперстроки (если W> 9) мотива EE. Таким образом, как мы решаем то, что правильная длина этого мотива? Мы можем ожидать, что оптимальная длина максимизирует различие в частоте между мотивом в целевом наборе и тем же мотивом в фоновом режиме набор. Однако для того, чтобы сравнить поле через различные длины, поле должно быть нормировано с учетом естественного стремления более коротких слов, чтобы происходить более часто. Мы выполняем эту нормализацию путем деления каждого поля на поле, соответствующее самому превалирующему слову идентичной длины в случайном наборе последовательностей с составом нуклеотида, похожим на целевой набор. Для удобства, главных превалирующих слов для длины W = 6... 11 и их поля хранятся в переменных topMotif
и topMargin
. Точно так же главные превалирующие слова для длины W = 6... 11 и их поля в случайном наборе хранятся в переменных rTopMotif
и rTopMargin
.
% === top over-represented words, W = 6...11 in set 2 (evening) topMotif topMargin % === top over-represented words, W = 6...11 in random set rTopMotif rTopMargin % === compute score score = topMargin ./ rTopMargin; [bestScore, bestLength] = max(score); % === plot figure() plot(6:11, score(6:11)); xlabel('Motif length'); ylabel('Normalized margin'); title('Optimal motif length'); hold on line([bestLength bestLength], [0 bestScore], 'LineStyle', '-.')
topMotif = 11x1 cell array {0x0 double } {0x0 double } {0x0 double } {0x0 double } {0x0 double } {'AATATC' } {'AATATCT' } {'AAATATCT' } {'AAAATATCT' } {'AAAAATATCT' } {'AAAAAATATCT'} topMargin = 1.0e-03 * NaN NaN NaN NaN NaN 0.3007 0.2607 0.2074 0.1439 0.0648 0.0424 rTopMotif = 11x1 cell array {0x0 double } {0x0 double } {0x0 double } {0x0 double } {0x0 double } {'ATTATA' } {'TATAATA' } {'TTATTAAA' } {'GTTATTAAA' } {'ATTATATATC' } {'ATGTTATTATT'} rTopMargin = 1.0e-03 * NaN NaN NaN NaN NaN 0.5650 0.2374 0.0972 0.0495 0.0279 0.0183
Путем графического вывода нормированного поля по сравнению с длиной мотива мы находим, что длина W = 9 является самой информативной в различении превалирующих мотивов в целевой последовательности (вечерний набор) против фонового набора (утренние и ночные наборы).
Несмотря на то, что Мотив EE был идентифицирован и экспериментально подтвержден как регулирующий мотив для генов, выражение которых достигает максимума в вечерние часы, это не совместно используется всеми вечерними генами, и при этом это не исключительно из этих генов. Мы считаем случаи мотива EE в трех наборах последовательности и определяем, какая пропорция генов в каждом наборе содержат мотив.
EECount = zeros(3,1); % === determine positions where EE motif occurs loc1 = strfind(seq1, EEMotif); loc2 = strfind(seq2, EEMotif); loc3 = strfind(seq3, EEMotif); % === count occurrences EECount(1) = length(loc1); EECount(2) = length(loc2); EECount(3) = length(loc3); % === find proportions of genes with EE Motif NumGenes = [N1; N2; N3] / 2; EEProp = EECount ./ NumGenes; % === plot figure() bar(EEProp, 0.5); ylabel('Proportion of genes containing EE Motif'); xlabel('Gene set'); title('Presence of EE Motif'); ylim([0 1]) ax = gca; ax.XTickLabel = {'morning', 'evening', 'night'};
Появляется, как будто приблизительно 9% генов в наборе 1, 40% генов в наборе 2, и 13% генов в наборе 3 имеют мотив EE. Таким образом не все гены в наборе 2 имеют мотив, но это ясно обогащено в этой группе.
В отличие от многих других функциональных мотивов, мотив EE, кажется, не накапливается в определенных генных местоположениях в наборе анализируемых последовательностей. После определения местоположения каждого вхождения относительно записи создает сайт (TSS) мы наблюдаем относительно равномерное распределение случаев через восходящую область рассмотренных генов за возможным исключением средней области (между 400 и 500 базисами в восходящем направлении TSS).
offset = rem(loc2, 1001); figure(); hist(offset, 100); xlabel('Offset in upstream region (TSS = 0)'); ylabel('Number of sequences');
[1] Cristianini, N. и Хан, M.W., "Введение в вычислительную геномику: подход тематических исследований", издательство Кембриджского университета, 2007.
[2] Harmer, S.L., и др., "Организованная Запись Ключевых Трасс в Arabidopsis Циркадными Часами", Наука, 290 (5499):2110-3, 2000.