Введение

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

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

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

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

Помимо мотива EE, другие слова длины W = 9, по-видимому, чрезмерно представлены в целевом наборе. В частности, мы замечаем наличие повторов, то есть слов, состоящих из одного нуклеотида или димера, повторяемых на весь размер слова, таких как 'CTCTCTC'. Это явление довольно распространено в геномных последовательностях и в целом связано с нефункциональными компонентами. Поскольку в этом контексте повторы вряд ли будут биологически значимыми, мы их отфильтровываем.

% === 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 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 имеют мотив 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. and Hahn, M.W., «Introduction to Computational Genomics: A Case Studies Approach», 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.