Введение

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

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

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

Создание Bootstrap репликатов из данных

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

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

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

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

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] Felsenstein, J., «Infering Phylogenies», Sinaur Associates, Inc., 2004.

[2] Nei, M. and Kumar, S., «Molecular Evolution and Phylogenetics», Oxford University Press. Глава 4, 2000.