Введение

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

Этот пример использует 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');

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

The oligoprop функция является полезным инструментом для получения свойств олигонуклеотидных последовательностей ДНК. Эта функция вычисляет содержание GC, молекулярный вес, температуру плавления и информацию о вторичной структуре. 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;

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

Температура плавления значительна, когда вы разрабатываете протоколы ПЦР. Создайте другой логический вектор индексации, чтобы отслеживать праймеры с плохой температурой плавления. Температура плавления от 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')

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

Можно использовать информацию, вычисленную до сих пор, чтобы найти лучшие пары праймеров, которые позволяют усилить область 220bp от положения 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] Санта-Люсия-младший, «Унифицированное представление полимера, гантели и олигонуклеотидной ДНК, ближайшей соседней термодинамики», PNAS, 95 (4): 1460-5, 1998.

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