Этот пример показывает, как использовать Bioinformatics Toolbox™, чтобы найти потенциальные капсюли, которые могут использоваться для автоматизированного упорядочивания DNA.
Проект краткой информации для PCR может быть грандиозной задачей. Хорошая краткая информация должна обычно соответствовать следующим критериям:
Длина является 18-30 основами.
Температура плавления составляет 50-60 градусов Цельсия.
Содержимое GC между 45% и 55%.
Не формирует стабильные шпильки.
Не делает сам димеризация.
Не пересекает димеризацию с другими капсюлями в реакции.
Имеет зажим GC в 3' концах краткой информации.
Этот пример использует MATLAB® и Bioinformatics Toolbox, чтобы найти капсюли PCR для человеческого hexosaminidase гена.
Сначала загрузите hexosaminidase последовательность нуклеотида из обеспеченного MAT-файла hexosaminidase.mat. Последовательность DNA, для которой вы хотите найти капсюли, находится в поле Sequence загруженной структуры.
load('hexosaminidase.mat','humanHEXA') sequence = humanHEXA.Sequence;
Можно также использовать функцию getgenbank, чтобы получить информацию о последовательности из репозитория данных NCBI и загрузить его в MATLAB. Ссылочная последовательность NCBI для HEXA имеет инвентарный номер NM_000520.
humanHEXA = getgenbank('NM_000520');
Функция oligoprop является полезным инструментом, чтобы получить свойства последовательностей DNA олигонуклеотида. Эта функция вычисляет содержимое 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 является длиной цели hexosaminidase последовательность. Затем можно использовать функцию 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;
Температура плавления является значительной, когда вы разрабатываете протоколы PCR. Создайте другой логический вектор индексации, чтобы отслеживать капсюли с плохими температурами плавления. Температуры плавления от oligoprop оцениваются во множестве путей (основной, настроенный солью, ближайшего соседа). Следующий пример использует оценки ближайшего соседа для температур плавления с параметрами, установленными SantaLucia, младшим [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';
Капсюли, которые имеют фрагменты повторных нуклеотидов, могут дать плохие результаты PCR. Это последовательности с низкой сложностью. Чтобы устранить капсюли с фрагментами четырех или больше повторных основ, используйте функциональный 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 - 1 100. Во-первых, вы находите все пары, которые могут покрыть необходимую область, учтя длину краткой информации. Затем вы вычисляете Евклидово расстояние фактических положений к желаемым единицам и переупорядочиваете список начиная с самого близкого расстояния.
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, чтобы сгенерировать отчет с десятью лучшими парами и сопоставленной информацией. Эти пары краткой информации могут затем быть проверены экспериментально. Эти капсюли могут также быть 'УНИЧТОЖЕНЫ' с помощью функции 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. Младший, "Объединенное представление полимера, гантели и DNA олигонуклеотида термодинамика ближайшего соседа", PNAS, 95 (4):1460-5, 1998.
[2] Робертс, R.J., и др., "REBASE - ферменты ограничения и метилтрансферазы DNA", Исследование Нуклеиновых кислот, 33:D230-2, 2005.