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

Введение

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

Загрузка данных о последовательности и создание исходного дерева

Этот пример использует 12 предварительно выровненных последовательностей, изолированных от различных hominidae разновидностей и сохраненных в FASTA-отформатированном файле. Филогенетическое дерево создается при помощи метода UPGMA с попарными расстояниями. А именно, функция seqpdist вычисляет попарные расстояния среди продуманных последовательностей, и затем функциональный seqlinkage создает дерево и возвращает данные в объекте phytree. Можно использовать функцию phytreeviewer, чтобы визуализировать и исследовать дерево.

primates = fastaread('primatesaligned.fa');
num_seqs = length(primates)

num_seqs = 12

orig_primates_dist = seqpdist(primates);
orig_primates_tree =  seqlinkage(orig_primates_dist,'average',primates);
phytreeviewer(orig_primates_tree);

Создание начальной загрузки реплицирует от данных

Начальная загрузка реплицирует, переставленное представление данных о последовательности DNA. Чтобы заставить начальную загрузку реплицировать данных о приматах, основы выбраны случайным образом от последовательностей с заменой и конкатенированы, чтобы сделать новые последовательности. То же количество основ как оригинал несколько, выравнивание используется в этом анализе, и затем разрывает, удалены, чтобы обеспечить новое попарное выравнивание. Функциональный randsample выбирает данные с заменой. Эта функция может также выбрать данные случайным образом без замены, чтобы выполнить анализ складного ножа. Для этого анализа 100 начальных загрузок реплицируют для каждой последовательности, создаются.

num_boots = 100;
seq_len = length(primates(1).Sequence);

boots = cell(num_boots,1);
for n = 1:num_boots
    reorder_index = randsample(seq_len,seq_len,true);
    for i = num_seqs:-1:1 %reverse order to preallocate memory
        bootseq(i).Header = primates(i).Header;
        bootseq(i).Sequence = strrep(primates(i).Sequence(reorder_index),'-',''); 
    end
    boots{n} = bootseq; 
end

Вычисление расстояний между начальными загрузками и филогенетической реконструкцией

Определение расстояний между последовательностями DNA для большого набора данных и создание филогенетических деревьев могут быть длительными. Распределение этих вычислений по нескольким машинам/ядрам уменьшает время вычисления. Этот пример принимает, что вы уже запустили пул MATLAB® с дополнительных параллельных ресурсов. Для получения информации о подготовке и выборе параллельных настроек, см. "Программирование с Пользовательскими Настройками" в документации Parallel Computing Toolbox™. Если у вас нет Parallel Computing Toolbox™, следующий цикл PARFOR выполняется последовательно без дальнейшей модификации.

fun = @(x) seqlinkage(x,'average',{primates.Header});
boot_trees = cell(num_boots,1);
parfor (n = 1:num_boots)
    dist_tmp = seqpdist(boots{n});
    boot_trees{n} = fun(dist_tmp);
end

Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 12).

Подсчет ответвлений с подобной топологией

Топология каждого дерева начальной загрузки по сравнению с тем из исходного дерева. Любое внутреннее ответвление, которое дает тот же раздел разновидностей, считается. Поскольку ответвления могут быть упорядочены по-другому среди различных деревьев, но все еще представлять тот же раздел разновидностей, необходимо получить каноническую форму для каждого поддерева перед сравнением. Первый шаг должен получить канонические поддеревья исходного дерева с помощью subtree и методов getcanonical от Bioinformatics Toolbox™.

for i = num_seqs-1:-1:1  % for every branch, reverse order to preallocate
    branch_pointer = i + num_seqs;
    sub_tree = subtree(orig_primates_tree,branch_pointer);
    orig_pointers{i} = getcanonical(sub_tree);
    orig_species{i} = sort(get(sub_tree,'LeafNames'));
end

Теперь можно получить канонические поддеревья для всех ответвлений каждого дерева начальной загрузки.

for j = num_boots:-1:1
    for i = num_seqs-1:-1:1  % for every branch
        branch_ptr = i + num_seqs;
        sub_tree = subtree(boot_trees{j},branch_ptr);
        clusters_pointers{i,j} = getcanonical(sub_tree);
        clusters_species{i,j} = sort(get(sub_tree,'LeafNames'));
    end
end

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

count = zeros(num_seqs-1,1);
for i = 1 : num_seqs -1  % for every branch
    for j = 1 : num_boots * (num_seqs-1)
        if isequal(orig_pointers{i},clusters_pointers{j})
            if isequal(orig_species{i},clusters_species{j})
                count(i) = count(i) + 1;
            end
        end
    end
end

Pc = count ./ num_boots   % confidence probability (Pc)

Pc = 11×1

    1.0000
    1.0000
    0.9900
    0.9900
    0.5400
    0.5400
    1.0000
    0.4300
    0.3900
    0.3900
      ⋮

Визуализация значений уверенности в исходном дереве

Информация об уверенности, сопоставленная с каждым узлом ответвления, может храниться в дереве путем аннотирования имен узла. Таким образом можно создать новое дерево, эквивалентное исходному дереву приматов, и аннотировать имена ответвления, чтобы включать доверительные уровни, вычисленные выше. phytreeviewer отображает эти данные во всплывающих подсказках, когда мышь наведена по внутренним узлам дерева.

[ptrs,dist,names] = get(orig_primates_tree,'POINTERS','DISTANCES','NODENAMES');

for i = 1:num_seqs -1  % for every branch
    branch_ptr = i + num_seqs;
    names{branch_ptr} = [names{branch_ptr} ', confidence: ' num2str(100*Pc(i)) ' %'];
end

tr = phytree(ptrs,dist,names)

    Phylogenetic tree object with 12 leaves (11 branches)

Можно выбрать узлы ответвления с доверительным уровнем, больше, чем данный порог, например, 0.9, и просмотреть эти соответствующие узлы в приложении Phylogenetic Tree. Можно также выбрать эти узлы ответвления в интерактивном режиме в рамках приложения.

high_conf_branch_ptr = find(Pc > 0.9) + num_seqs;
view(tr, high_conf_branch_ptr)

Ссылки

[1] Фелзенштайн, J., "выводя филогении", Sinaur Associates, Inc., 2004.

[2] Nei, M. и Кумар, S., "Молекулярный Evolution и Phylogenetics", издательство Оксфордского университета. Глава 4, 2000.