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

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

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

Тем не менее, тем не менее, задача может быть выполнена с помощью 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).

Answer 2

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