Введение

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

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

В этом примере используются 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 выполняет выборку данных с заменой. Эта функция также может случайным образом выборить данные без замены для выполнения анализа jackknife. Для этого анализа создается 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 ® уже запущен с дополнительными параллельными ресурсами. Сведения о настройке и выборе параллельных конфигураций см. в разделе «Программирование с пользовательскими конфигурациями» в документации по Toolbox™ параллельных вычислений. Если 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 методы из 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, и просмотреть эти соответствующие узлы в приложении Филогенетическое дерево. Можно также выбрать эти узлы ветвей в интерактивном режиме в приложении.

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

Ссылки

[1] Felsenstein, J., «Вывод филогении», Sinaur Associates, Inc., 2004.

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