plot PLINK2.0 output with ggtree

I've been working on an IBS matrix using PLINK2.0 and having the necessity to plot it using R. I've seen some libraries out there e.g. ape and ggtree and wanted to use the latter for the greater number of features.

So, this is how the matrix looks like:

1.53578 -0.105267       -0.0945413      -0.083604       -0.209443       -0.139285       -0.208677       -0.199941       -0.209932       -0.130887       -0.082297       -0.0719018
-0.105267       1.60124 -0.142996       -0.023115       -0.243876       -0.0435319      -0.243711       -0.256835       -0.255688       -0.132416       -0.102566       -0.0512363
-0.0945413      -0.142996       1.70865 -0.00168153     -0.229154       -0.0920099      -0.201994       -0.278662       -0.210596       -0.190707       -0.164444       -0.101866
-0.083604       -0.023115       -0.00168153     1.36444 -0.176534       -0.0313074      -0.178905       -0.256911       -0.220043       -0.199516       -0.16628        -0.0265423
-0.209443       -0.243876       -0.229154       -0.176534       1.66881 -0.223521       0.672806        -0.346079       -0.25505        -0.243328       -0.219959       -0.194674
-0.139285       -0.0435319      -0.0920099      -0.0313074      -0.223521       1.6054  -0.209163       -0.304027       -0.216554       -0.225491       -0.186448       0.0659375
-0.208677       -0.243711       -0.201994       -0.178905       0.672806        -0.209163       1.63191 -0.329851       -0.252808       -0.252563       -0.233733       -0.193315
-0.199941       -0.256835       -0.278662       -0.256911       -0.346079       -0.304027       -0.329851       2.72659 -0.253858       -0.0673777      -0.158547       -0.274507
-0.209932       -0.255688       -0.210596       -0.220043       -0.25505        -0.216554       -0.252808       -0.253858       2.45426 -0.233602       -0.159556       -0.186578
-0.130887       -0.132416       -0.190707       -0.199516       -0.243328       -0.225491       -0.252563       -0.0673777      -0.233602       1.91195 -0.0165982      -0.219464
-0.082297       -0.102566       -0.164444       -0.16628        -0.219959       -0.186448       -0.233733       -0.158547       -0.159556       -0.0165982      1.64659 -0.156163
-0.0719018      -0.0512363      -0.101866       -0.0265423      -0.194674       0.0659375       -0.193315       -0.274507       -0.186578       -0.219464       -0.156163       1.41031

Now, thanks to a post on SO I managed to get it into a format suitable for plotting using ape (see pic. below) — although it is specified the phylo format should be compatible also with ggtree. However, I can't manage to plot it as such...

Also, it appears the post says this is an unconventional way to go from matrix to tree and might not be reflective of the actual phylogenesis behind. So, is there another way around this, always using an IBS matrix generated by PLINK2.0? I'm happy to share the command I used if needed. Thanks in advance!

tree

ibs matrix plink ggtree

• 972 views • link

written 1.4 years ago by Matteo Ungaro •

I actually solved it simply loading the phylo object with ggtree. Below, the full code with some plot features:

###LOAD DATA AND WRANGLING
ibs_matrix = read.delim("/path/to/phylo_tree_header.phy", sep="\t", header=TRUE)
colnames(ibs_matrix)[1] <- ""
ibs_matrix[1] <- NULL
#ibs_matrix

ibs_matrix_t <- t(ibs_matrix)
#ibs_matrix_t


###ADD META INFO AND DF FORMATTING
variety <- c("wt", "lr", "wt", "wt", "lr", "wt", "wt", "wt", "lr", "lr", "lr", "lr")
location <- c("ESP", "ESP", "ITA", "ITA", "JOR", "TUR", "PAL", "GRC", "GRC", "ESP", "ESP", "ITA")

#meta <- list(variety="wt", "lr", "wt", "wt", "lr", "wt", "wt", "wt", "lr", "lr", "lr", "lr", location="ESP", "ESP", "ITA", "ITA", "JOR", "TUR", "PAL", "GRC", "GRC", "ESP", "ESP", "ITA")

meta_df <- data.frame(ibs_matrix_t[, 1], variety, location); meta_df <- meta_df[ -c(1) ]
meta_df$id <- rownames(meta_df); meta_df <- meta_df[,c(3,1,2)]
rownames(meta_df) <- NULL


###COMPUTE AND CALCULATE PIRWISE DISTANCES, PLUS CONVERSION TO PHYLO
pdist <- dist(ibs_matrix_t, method="euclidean", diag=TRUE, upper=TRUE, p=2)
distances <- distanceHadamard(as.matrix(pdist))
#lento(distances)

phylo <- as.phylo(distances)


###BASIC APE PLOT
#ape::plot.phylo(phylo)


###GGTREE PERSONALIZATION
t1 <- ggtree(phylo) %<+% meta_df +
  geom_tippoint(aes(color=variety)) + geom_tiplab(size=3)
t1

ggtree

• 0 views • link

written 1.4 years ago by Matteo Ungaro •

0 answers

No answers yet.