Этот пример показывает, как использовать Bioinformatics Toolbox™ для поиска потенциальных праймеров, которые могут использоваться для автоматического секвенирования ДНК.
Проект праймера для ПЦР может быть сложной задачей. Хороший праймер обычно должен соответствовать следующим критериям:
Длина составляет 18 - 30 основы.
Температура плавления 50-60 степени Цельсия.
Содержимое GC составляет от 45% до 55%.
Не образует устойчивых шпильек.
Не димеризуется сам.
Не пересекает димеризацию с другими праймерами в реакции.
Имеет зажим GC на 3 '-конце грунтовки.
Этот пример использует 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 от 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}');
Полезно сильная основа, спаривающееся на 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.