Начальная загрузка филогенетических деревьев

То В этом примере показано, как сгенерировать начальную загрузку, реплицирует последовательностей ДНК. Данные, сгенерированные начальной загрузкой, используются, чтобы оценить доверие ветвей в филогенетическом дереве.

Введение

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

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

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

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

Начальная загрузка реплицирует, переставленное представление данных о последовательности ДНК. Чтобы заставить начальную загрузку реплицировать данных о приматах, основы произведены случайным образом от последовательностей с заменой и конкатенированы, чтобы сделать новые последовательности. То же количество основ как оригинал несколько, выравнивание используется в этом анализе, и затем разрывает, удалены, чтобы обеспечить новое попарное выравнивание. Функциональный 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

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

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

fun = @(x) seqlinkage(x,'average',{primates.Header});
boot_trees = cell(num_boots,1);
parpool('local');
Starting parallel pool (parpool) using the 'local' profile ...
parfor (n = 1:num_boots)
    dist_tmp = seqpdist(boots{n});
    boot_trees{n} = fun(dist_tmp);
end
delete(gcp('nocreate'));

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

Топология каждого дерева начальной загрузки по сравнению с тем из исходного дерева. Любая внутренняя ветвь, которая дает тот же раздел разновидностей, считается. Поскольку ветви могут быть упорядочены по-другому среди различных деревьев, но все еще представлять тот же раздел разновидностей, необходимо получить каноническую форму для каждого поддерева перед сравнением. Первый шаг должен получить канонические поддеревья исходного дерева с помощью 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., "Молекулярная эволюция и филогенетика", издательство Оксфордского университета. Глава 4, 2000.