Документация

Введение

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

Этот пример использует 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. Чтобы устранить все потенциальные капсюли, которые не соответствуют вышеизложенным критериям (содержимое 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}');

Нахождение капсюлей без зажима 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';

Нахождение капсюлей с повторениями нуклеотида

Капсюли, которые имеют фрагменты повторных нуклеотидов, могут дать плохие результаты 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.