Как сгенерировать вывод дерева Ньюика из матрицы попарных расстояний

Question

Я хотел бы получить филогенетические деревья из генетических данных.Я нашел несколько пакетов для рисования деревьев в R и Python, которые выглядят великолепно, например, ggtree в R. Но они требуют ввода данных, которые уже в древовидном формате, например, Newick.

Я думаю, что большинство людей начинают с vcfфайлы и производят файлы FASTA, но моя отправная точка - таблица генотипов - я работаю с гаплоидным организмом, поэтому каждая позиция - 0 (ref) или 1 (non-ref).Из этого я вычисляю попарно генетическое расстояние, используя dist () в R. Пример данных для 5 выборок, AE, с попарным расстоянием по десяти вариантным позициям:

# Generate dataframe with example genotypes
Variant <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
A <- c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0)
B <- c(1, 1, 0, 0, 1, 1, 0, 0, 1, 1)
C <- c(0, 0, 1, 1, 0, 0, 1, 1, 0, 1)
D <- c(1, 0, 1, 1, 0, 0, 1, 1, 0, 1)
E <- c(1, 0, 0, 0, 1, 1, 0, 0, 1, 1)
df = data.frame(Variant, A, B, C, D, E)
df
#  Remove first column with variant names
df$Variant <- NULL
# Transpose the columns and rows
df_t = t(df)
# Compute pairwise distance (Euclidean)
pdist = dist(df_t, method = "euclidean", diag = TRUE, upper = TRUE, p = 2)
pdist

Я хотел бы создать выходной файл иерархического дерева изpdist, например, в формате Newick, так что я могу подключить его к пакетам, таким как ggtree, для рисования красивых деревьев, например, круглых филограмм, раскрашенных ко-вариациями и т. д. Я пробовал искать, но не уверен, с чего начать.

РЕДАКТИРОВАТЬ/ ОБНОВЛЕНИЕ Этот сайт был полезен http://www.phytools.org/Cordoba2017/ex/2/Intro-to-phylogenies.html Я использовал пакеты: ape, phangorn, phytools, geiger

Этот код работает -

# Produce dendrogram
hclust = hclust(pdist)
# Check dendrogram looks sensible
plot(hclust)
class(hclust) # check that class is hclust
# Save to Newick file
my_tree <- as.phylo(hclust) 
write.tree(phy=my_tree, file="ExampleTree.newick") # Writes a Newick file
# Produce tree
plot(unroot(my_tree),type="unrooted",cex=1.5,
     use.edge.length=TRUE,lab4ut="axial",
     edge.width=2,
     no.margin=TRUE)

Дерево вывода:

Answer 1

Двоичные данные могут быть эффективно использованы для построения филогенетического дерева с помощью MrBayes, как упоминалось @ thomas-guillerme.Входной файл должен включать в себя двоичный data блок и mrbayes команды.

#nexus
begin data;
dimensions ntax = 5 nchar = 10;
format datatype = restriction;
matrix
A 0011001100
B 1100110011
C 0011001101
D 1011001101
E 1000110011;
end;

begin mrbayes;
lset coding = variable;
mcmc ngen = 1000000 samplefreq = 1000;
sump burnin = 200;
sumt burnin = 200;
end;

Длина прогона mcmc должна быть скорректирована с учетом сходимости цепи.Для начала код должен дать хорошее представление об отношениях, которые могут вывести данные.

Answer 2

Это нетривиальное задание.Чтобы построить дерево (как в раздвоенном) из матрицы расстояний, вам нужно будет использовать филогенетические алгоритмы и, возможно, лучше не делать это из матрицы расстояний (обратите внимание, что могут быть недостатки использования евклидова расстояния для двоичной матрицы).).

Однако, тем не менее, задача все еще может быть выполнена с использованием пакета phangorn .Например, вы можете создать спектры расщеплений из матрицы расстояний (то есть вероятных расщеплений, присутствующих в матрице ( подробнее здесь - paywalled ).

require(phangorn)
# Generate dataframe with example genotypes
Variant <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
A <- c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0)
B <- c(1, 1, 0, 0, 1, 1, 0, 0, 1, 1)
C <- c(0, 0, 1, 1, 0, 0, 1, 1, 0, 1)
D <- c(1, 0, 1, 1, 0, 0, 1, 1, 0, 1)
E <- c(1, 0, 0, 0, 1, 1, 0, 0, 1, 1)
df = data.frame(Variant, A, B, C, D, E)
df
#  Remove first column with variant names
df$Variant <- NULL
# Transpose the columns and rows
df_t = t(df)
# Compute pairwise distance (Euclidean)
pdist = dist(df_t, method = "euclidean", diag = TRUE, upper = TRUE, p = 2)

# calculate the Hadamard distance spectrum
distances <- distanceHadamard(as.matrix(pdist))
# representing the distances
lento(distances)
# Plotting the distances as a tree (a network actually)
plot(as.networx(distances), "2D")

Обратите внимание, что втакой же пакет neighborNet также доступен, но в руководстве подчеркивается, что эта функция является экспериментальной. Я предлагаю связаться с автором пакета для получения дополнительной информации.

Затем вы можете преобразовать свою сеть в "phylo", который может использоватьсяape и, вероятно, ggtree, приведя его в действие:

# Converting into a phylo object
phylo <- as.phylo(distances)

Но, опять же, обратите внимание, что это результирующее дерево, вероятно, неверно в филогенетическом смысле (т.е. предполагает спуск с модификацией), и я настоятельно рекомендую просто оценитьдерево, использующее модельный подход (например, с MrBayes или BEAST2 ).