Введение

Циркадные часы - это 24-часовой цикл физиологических процессов, которые синхронизируются с внешним циклом день-ночь. Большая часть работ по циркадному осциллятору в растениях выполнена с использованием модельного растения Arabidopsis thaliana. В этом организме регуляция ряда генов, которые должны быть включены или выключены в определенное время дня и ночи, осуществляется небольшими регуляторными последовательностями, обнаруженными выше по потоку рассматриваемых генов. Один такой регулирующий мотив, AAAATATCT, также известный как Evening Element (EE), был определен в областях покровителя циркадных отрегулированных часами генов, которые показывают пиковое выражение вечером [2].

Загрузка восходящих областей генов, регулируемых часами

Рассмотрим три набора регулируемых часами генов, кластеризованных по времени суток, когда они максимально экспрессированы: множество 1 соответствует 1 КБ-длинным восходящим областям генов, пик экспрессии которых приходится на утро (8 часов утра -16 часов); набор 2 соответствует 1 KB-длинным восходящим областям генов, пик экспрессии которых приходится на вечер (4 п.м.-12); набор 3 соответствует 1 КБ-длинным восходящим областям генов, пик экспрессии которых приходится на ночь (12 часов ночи). Поскольку мы заинтересованы в регуляторном мотиве в вечерних генах, набор 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-меры (слова длины 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', также известный как мотив Evening Element (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

Фильтрация повторов

Помимо мотива ЭЭ, другие слова длины W = 9, по-видимому, чрезмерно представлены в целевом наборе. В частности, мы отмечаем наличие повторов, то есть слов, состоящих из одного нуклеотида или димера, повторяющихся для всей длины слова, таких как «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. Геномные области 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 из 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 Motif был идентифицирован и экспериментально подтвержден как регуляторный мотив для генов, экспрессия которых достигает пика в вечерние часы, он не разделяется всеми вечерними генами и не исключается из этих генов. Мы подсчитываем появление мотива 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 имеют мотив ЭЭ. Таким образом, не все гены в наборе 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] Криштианини, Н. и Хан, М.В., «Введение в вычислительную геномику: подход к тематическим исследованиям», Cambridge University Press, 2007.

[2] Harmer, S.L., et al., «Orchestrated Transcription of Key Pathways in Arabidopsis by the Circadian Clock», Science, 290 (5499): 2110-3, 2000.