Введение

Разработка праймеров для ПЦР может быть сложной задачей. Хороший праймер обычно должен отвечать следующим критериям:

В этом примере используются MATLAB ® и Bioinformatics Toolbox для поиска праймеров ПЦР для гена гексосаминидазы человека.

Сначала загружают нуклеотидную последовательность гексосаминидазы из предоставленного MAT-файла hexosaminidase.mat. Последовательность ДНК, для которой вы хотите найти праймеры, находится в Sequence поле загруженной структуры.

load('hexosaminidase.mat','humanHEXA')
sequence = humanHEXA.Sequence;

Вы также можете использовать getgenbank функция для извлечения информации о последовательности из хранилища данных NCBI и ее загрузки в MATLAB. Ссылочная последовательность NCBI для HEXA имеет регистрационный номер NM_000520.

humanHEXA = getgenbank('NM_000520');

Расчет свойств олигонуклеотида

oligoprop функция является полезным инструментом для получения свойств последовательностей олигонуклеотидной ДНК. Эта функция вычисляет содержание ГХ, молекулярную массу, температуру плавления и информацию о вторичной структуре. oligoprop возвращает структуру, содержащую поля со связанными сведениями. Используйте help команда для просмотра других свойств oligoprop возвращает.

nt = oligoprop('AAGCTCAAAAACGCGCGGTATTCGACTGGCGTGATCTATTTTATGCT')

nt = 

  struct with fields:

                GC: 44.6809
           GCdelta: 0
          Hairpins: [3x47 char]
            Dimers: [9x47 char]
         MolWeight: 1.4468e+04
    MolWeightdelta: 0
                Tm: [68.9128 79.7752 85.3393 69.6497 68.2474 75.8931]
           Tmdelta: [0 0 0 0 0 0]
            Thermo: [4x3 double]
       Thermodelta: [4x3 double]

Поиск всех потенциальных праймеров вперед

Целью является создание списка всех потенциальных прямых праймеров длиной 20. Это можно сделать, выполнив итерацию по всей последовательности и заняв подпоследовательности в любой возможной позиции, или используя матрицу индексов. В следующем примере показано, как можно установить матрицу индексов и затем использовать ее для создания всех возможных прямых подпоследовательностей длиной 20, в данном случае подпоследовательности N-19 где N - длина целевой последовательности гексосаминидазы. Тогда вы можете использовать oligoprop чтобы получить свойства для каждого из потенциальных праймеров в списке.

N = length(sequence) % length of the target sequence
M = 20  % desired primer length
index = repmat((0:N-M)',1,M)+repmat(1:M,N-M+1,1);
fwdprimerlist = sequence(index);

for i = N-19:-1:1 % reverse order to pre-allocate structure
    fwdprimerprops(i) = oligoprop(fwdprimerlist(i,:));
end

N =

        2437


M =

    20

Поиск всех потенциальных обратных праймеров

После того, как у вас есть все потенциальные праймеры вперед, вы можете искать обратные праймеры аналогичным образом. Обратные праймеры находятся на комплементарной цепи. Для получения комплементарной нити используйте seqcomplement функция.

comp_sequence = seqcomplement(sequence);
revprimerlist = seqreverse(comp_sequence(index));

for i = N-19:-1:1 % reverse order to preallocate structure
    revprimerprops(i) = oligoprop(revprimerlist(i,:));
end

Фильтрация праймеров на основе содержимого GC

Информация о содержании GC для праймеров находится в структуре с полем GC. Чтобы исключить все потенциальные праймеры, которые не соответствуют критериям, изложенным выше (содержание GC от 45% до 55%), можно сделать вектор логического индексирования, который указывает, какие праймеры имеют содержание GC вне допустимого диапазона. Извлеките поле GC из структуры и преобразуйте его в числовой вектор.

fwdgc = [fwdprimerprops.GC]';
revgc = [revprimerprops.GC]';

bad_fwdprimers_gc = fwdgc < 45 | fwdgc > 55;
bad_revprimers_gc = revgc < 45 | revgc > 55;

Фильтрующие грунтовки на основе их температуры плавления

Температура плавления является значительной при разработке протоколов PCR. Создайте еще один вектор логической индексации для отслеживания праймеров с плохими температурами плавления. Температуры плавления от oligoprop оцениваются различными способами (базовый, с поправкой на соль, ближайший сосед). В следующем примере используются оценки ближайшего соседа для температур плавления с параметрами, установленными SantaLucia, Jr. [1]. Они хранятся в пятом элементе поля Tm возвращено oligoprop. Другие элементы этого поля представляют другие способы оценки температуры плавления. Вы также можете использовать mean функция для вычисления среднего значения по всем оценкам.

fwdtm = cell2mat({fwdprimerprops.Tm}');
revtm = cell2mat({revprimerprops.Tm}');
bad_fwdprimers_tm = fwdtm(:,5) < 50 | fwdtm(:,5) > 60;
bad_revprimers_tm = revtm(:,5) < 50 | revtm(:,5) > 60;

Поиск праймеров с самодимеризацией и образованием шпильки

Самосимеризация и образование шпильки могут предотвратить связывание праймера с последовательностью-мишенью. Как и выше, можно создать логические индексирующие векторы, чтобы указать, образуют ли потенциальные праймеры самоликвидаторы или шпильки.

bad_fwdprimers_dimers  = ~cellfun('isempty',{fwdprimerprops.Dimers}');
bad_fwdprimers_hairpin = ~cellfun('isempty',{fwdprimerprops.Hairpins}');
bad_revprimers_dimers  = ~cellfun('isempty',{revprimerprops.Dimers}');
bad_revprimers_hairpin = ~cellfun('isempty',{revprimerprops.Hairpins}');

Поиск праймеров без зажима GC

Прочное спаривание оснований на 3 '-конце праймера полезно. Найдите все праймеры, которые не заканчиваются на G или C. Помните, что все последовательности в списках 5 '- > 3'.

bad_fwdprimers_clamp = lower(fwdprimerlist(:,end)) == 'a' | lower(fwdprimerlist(:,end)) == 't';
bad_revprimers_clamp = lower(revprimerlist(:,end)) == 'a' | lower(revprimerlist(:,end)) == 't';

Поиск праймеров с нуклеотидными повторами

Праймеры, которые имеют растяжение повторных нуклеотидов, могут давать плохие результаты ПЦР. Это последовательности с низкой сложностью. Для устранения праймеров с растяжениями четырех и более повторяющихся оснований используйте функцию regexp.

fwdrepeats = regexpi(cellstr(fwdprimerlist),'a{4,}|c{4,}|g{4,}|t{4,}','ONCE');
revrepeats = regexpi(cellstr(revprimerlist),'a{4,}|c{4,}|g{4,}|t{4,}','ONCE');
bad_fwdprimers_repeats = ~cellfun('isempty',fwdrepeats);
bad_revprimers_repeats = ~cellfun('isempty',revrepeats);

Найти праймеры, удовлетворяющие всем критериям

Строки исходного списка подпоследовательностей соответствуют базовому номеру, с которого начинается каждая подпоследовательность. Можно использовать собранные до сих пор логические векторы индексирования и создать новый список праймеров, удовлетворяющих всем описанным выше критериям. На рисунке показано, как фильтруются прямые праймеры, где значения, равные 1, указывают на плохие праймеры, а значения, равные 0, указывают на хорошие праймеры.

bad_fwdprimers = [bad_fwdprimers_gc, bad_fwdprimers_tm,...
                  bad_fwdprimers_dimers, bad_fwdprimers_hairpin,...
                  bad_fwdprimers_clamp, bad_fwdprimers_repeats];
bad_revprimers = [bad_revprimers_gc, bad_revprimers_tm,...
                  bad_revprimers_dimers, bad_revprimers_hairpin,...
                  bad_revprimers_clamp, bad_revprimers_repeats];

good_fwdpos = find(all(~bad_fwdprimers,2));
good_fwdprimers = fwdprimerlist(good_fwdpos,:);
good_fwdprop = fwdprimerprops(good_fwdpos);
N_good_fwdprimers = numel(good_fwdprop)

good_revpos = find(all(~bad_revprimers,2));
good_revprimers = revprimerlist(good_revpos,:);
good_revprop = revprimerprops(good_revpos);
N_good_revprimers = numel(good_revprop)

figure
imagesc([bad_fwdprimers any(bad_fwdprimers,2)]);
title('Filtering candidate forward primers');
ylabel('Primer location');
xlabel('Criteria');
ax = gca;
ax.XTickLabel = char({'%GC','Tm','Dimers','Hairpin','GC clamp','Repeats','All'});
ax.XTickLabelRotation = 45;
colorbar

N_good_fwdprimers =

   140


N_good_revprimers =

   147

Проверка кросс-димеризации

Перекрестная димеризация может происходить между прямой и обратной праймерами, если они имеют значительную степень комплементарности. Праймеры не будут функционировать должным образом, если они димеризуются друг с другом. Чтобы проверить димеризацию, совместите каждый прямой праймер с каждым обратным праймером, используя swalign и сохранить пары праймеров с низким баллом. Эта информация может храниться в матрице со строками, представляющими прямые праймеры, и столбцами, представляющими обратные праймеры. Этот исчерпывающий расчет может занять довольно много времени. Однако нет смысла выполнять этот расчет на парах праймеров, где обратный праймер находится выше переднего праймера. Поэтому эти пары праймеров можно игнорировать. Изображение на рисунке показывает парные оценки перед порогом, низкие оценки (темно-синие) представляют пары праймеров, которые не димеризуются.

scr_mat = [-1,-1,-1,1;-1,-1,1,-1;-1,1,-1,-1;1,-1,-1,-1;];
scr = zeros(N_good_fwdprimers,N_good_revprimers);
for i = 1:N_good_fwdprimers
    for j = 1:N_good_revprimers
        if good_fwdpos(i) < good_revpos(j)
            scr(i,j) = swalign(good_fwdprimers(i,:), good_revprimers(j,:), ...
                              'SCORINGMATRIX',scr_mat,'GAPOPEN',5,'ALPHA','NT');
        else
            scr(i,j) = 13; % give a high score to ignore forward primers
                           % that are after reverse primers
        end
    end
end

figure
imagesc(scr)
title('Cross dimerization scores')
xlabel('Candidate reverse primers')
ylabel('Candidate forward primers')
colorbar

Пары праймеров с низкой оценкой идентифицируют как логическую в индикаторной матрице.

pairedprimers = scr<=3;

Визуализация потенциальных пар праймеров в области последовательности

Альтернативным способом представления этой информации является просмотр всех потенциальных комбинаций праймеров в области последовательности. Каждая точка на графике представляет возможную комбинацию между прямым и обратным праймерами после фильтрации всех этих случаев с потенциальной перекрестной димеризацией.

[f,r] = find(pairedprimers);
figure
plot(good_revpos(r),good_fwdpos(f),'r.','markersize',10)
axis([1 N 1 N])
title('Primer selection graph')
xlabel('Reverse primer positions')
ylabel('Forward primer positions')

Выбор пары праймеров для усиления определенной области

Можно использовать информацию, рассчитанную до сих пор, чтобы найти лучшие пары праймеров, которые позволяют усилить область 220 bp из положения 880 в положение 1100. Сначала вы найдете все пары, которые могут покрыть требуемую область с учетом длины грунтовки. Затем вычисляется евклидово расстояние между фактическими позициями и требуемыми и переупорядочивается список, начиная с ближайшего расстояния.

pairs = find(good_fwdpos(f)<(880-M) & good_revpos(r)>1100);
dist = (good_fwdpos(f(pairs))-(880-M)).^2 + (good_revpos(r(pairs))-(1100)).^2;
[dist,h] = sort(dist);
pairs = pairs(h);

hold on
plot(good_revpos(r(pairs)),good_fwdpos(f(pairs)),'b.','markersize',10)
plot([1100 1100],[1 880-M],'g')
plot([1100 N],[880-M 880-M],'g')

Извлечение пар праймеров

Используйте sprintf для создания отчета с десятью лучшими парами и связанной информацией. Эти пары праймеров затем могут быть проверены экспериментально. Эти праймеры также могут быть «BLASTED» с использованием blastncbi функция проверки специфичности.

Primers = sprintf('Fwd/Rev Primers      Start End   %%GC   mT   Length\n\n');
for i = 1:10
    fwd = f(pairs(i));
    rev = r(pairs(i));
    Primers = sprintf('%s%-21s%-6d%-6d%-4.4g%-4.4g\n%-21s%-6d%-6d%-4.4g%-7.4g%-6d\n\n', ...
    Primers, good_fwdprimers(fwd,:),good_fwdpos(fwd),good_fwdpos(fwd)+M-1,good_fwdprop(fwd).GC,good_fwdprop(fwd).Tm(5), ...
             good_revprimers(rev,:),good_revpos(rev)+M-1,good_revpos(rev),good_revprop(rev).GC,good_revprop(rev).Tm(5), ...
             good_revpos(rev) - good_fwdpos(fwd) );
end
disp(Primers)

Fwd/Rev Primers      Start End   %GC   mT   Length

tacatctcgccattacctgc 732   751   50  55.61
tcaacctcatctcctccaag 1181  1162  50  54.8   430   

atacatctcgccattacctg 731   750   45  52.87
tcaacctcatctcctccaag 1181  1162  50  54.8   431   

tacatctcgccattacctgc 732   751   50  55.61
aaatcaacctcatctcctcc 1184  1165  45  52.9   433   

tacatctcgccattacctgc 732   751   50  55.61
gaaatcaacctcatctcctc 1185  1166  45  51.08  434   

atacatctcgccattacctg 731   750   45  52.87
aaatcaacctcatctcctcc 1184  1165  45  52.9   434   

atacatctcgccattacctg 731   750   45  52.87
gaaatcaacctcatctcctc 1185  1166  45  51.08  435   

ggatacatctcgccattacc 729   748   50  53.45
tcaacctcatctcctccaag 1181  1162  50  54.8   433   

tacatctcgccattacctgc 732   751   50  55.61
gtgaaatcaacctcatctcc 1187  1168  45  51.63  436   

tacatctcgccattacctgc 732   751   50  55.61
ggtgaaatcaacctcatctc 1188  1169  45  51.63  437   

atacatctcgccattacctg 731   750   45  52.87
gtgaaatcaacctcatctcc 1187  1168  45  51.63  437   

Найти ферменты рестрикции, которые режут внутри праймера

Используйте rebasecuts функция для перечисления всех рестрикционных ферментов из базы данных REBASE ® [2], которые будут резать праймер. Эти рестрикционные ферменты могут быть использованы при разработке экспериментов клонирования. Например, это можно использовать для первой пары праймеров из только что рассчитанного списка возможных праймеров.

fwdprimer = good_fwdprimers(f(pairs(1)),:)
fwdcutter = unique(rebasecuts(fwdprimer))

revprimer = good_revprimers(r(pairs(1)),:)
revcutter = unique(rebasecuts(revprimer))

fwdprimer =

    'tacatctcgccattacctgc'


fwdcutter =

  14x1 cell array

    {'AbaSI' }
    {'Acc36I'}
    {'BfuAI' }
    {'BmeDI' }
    {'BspMI' }
    {'BveI'  }
    {'FspEI' }
    {'LpnPI' }
    {'MspJI' }
    {'RlaI'  }
    {'SetI'  }
    {'SgeI'  }
    {'SgrTI' }
    {'YkrI'  }


revprimer =

    'tcaacctcatctcctccaag'


revcutter =

  12x1 cell array

    {'AbaSI' }
    {'AspBHI'}
    {'BmeDI' }
    {'BsaXI' }
    {'FspEI' }
    {'MnlI'  }
    {'MspJI' }
    {'RlaI'  }
    {'SetI'  }
    {'SgeI'  }
    {'SgrTI' }
    {'YkrI'  }

Ссылки

[1] SantaLucia, J. Jr., «Унифицированный взгляд на термодинамику полимера, гантели и олигонуклеотидной ДНК ближайшего соседа», PNAS, 95 (4): 1460-5, 1998.

Roberts, R.J., et al., «REBASE - рестрикционные ферменты и ДНК-метилтрансферазы», Nucleic Acids Research, 33:D230-2, 2005.