Введение

Циркадные часы являются 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.