forked from YuLab-SMU/treedata-book
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path11_ggtree_gallery.Rmd
217 lines (170 loc) · 10.2 KB
/
11_ggtree_gallery.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
# Gallery of reproducible examples {#chapter11}
## Visualizing pairwise nucleotide sequence distance with phylogenetic tree
This example reproduces Fig. 1 of [@chen_ancient_2017]. It extracts accession numbers from tip labels of the HPV58 tree and calculates pairwise nucleotide sequence distances. The distance matrix is visualized as dot and line plots. This example demonstrates the abilities of adding multiple layers to a specific panel. As illustrated on Figure \@ref(fig:jv2017), the `facet_plot` function displays sequence distances as a dot plot and then adds a layer of line plot to the same panel, *i.e.* Sequence Distance. In addition, the tree in `facet_plot` can be fully annotated with multiple layers (clade labels, bootstrap support values, *etc*). The source code of this example was firstly published in Supplemental File of [@yu_two_2018].
```{r message=FALSE}
library(tibble)
library(tidyr)
library(Biostrings)
library(treeio)
library(ggplot2)
library(ggtree)
hpvtree <- paste0("https://raw.githubusercontent.com/GuangchuangYu/",
"plotting_tree_with_data/master/HPV58.tree")
tree <- read.tree(hpvtree)
clade <- c(A3 = 92, A1 = 94, A2 = 108, B1 = 156, B2 = 159, C = 163, D1 = 173, D2 = 176)
tree <- groupClade(tree, clade)
cols <- c(A1 = "#EC762F", A2 = "#CA6629", A3 = "#894418", B1 = "#0923FA",
B2 = "#020D87", C = "#000000", D1 = "#9ACD32",D2 = "#08630A")
## visualize the tree with tip labels and tree scale
p <- ggtree(tree, aes(color = group), ladderize = FALSE) %>% rotate(rootnode(tree)) +
geom_tiplab(aes(label = paste0("italic('", label, "')")), parse = TRUE, size = 2.5) +
geom_treescale(x = 0, y = 1, width = 0.002) +
scale_color_manual(values = c(cols, "black"), na.value = "black", name = "Lineage",
breaks = c("A1", "A2", "A3", "B1", "B2", "C", "D1", "D2")) +
guides(color = guide_legend(override.aes = list(size = 5, shape = 15))) +
theme_tree2(legend.position = c(.1, .88))
## Optional
## add labels for monophyletic (A, C and D) and paraphyletic (B) groups
p <- p + geom_cladelabel(94, "italic(A1)", color = cols[["A1"]], offset = .003, align = TRUE,
offset.text = -.001, barsize = 1.2, extend = c(0, 0.5), parse = TRUE) +
geom_cladelabel(108, "italic(A2)", color = cols[["A2"]], offset = .003, align = TRUE,
offset.text = -.001, barsize = 1.2, extend = 0.5, parse = TRUE) +
geom_cladelabel(131, "italic(A3)", color = cols[["A3"]], offset = .003, align = TRUE,
offset.text = -.001, barsize = 1.2, extend = c(0.5, 0), parse = TRUE) +
geom_cladelabel(92, "italic(A)", color = "darkgrey", offset = .00315, align = TRUE,
offset.text = 0.0002, barsize = 2, fontsize = 5, parse = TRUE) +
geom_cladelabel(156, "italic(B1)", color = cols[["B1"]], offset = .003, align = TRUE,
offset.text = -.001, barsize = 1.2, extend = c(0, 0.5), parse = TRUE) +
geom_cladelabel(159, "italic(B2)", color = cols[["B2"]], offset = .003, align = TRUE,
offset.text = -.001, barsize = 1.2, extend = c(0.5, 0), parse = TRUE) +
geom_strip(65, 71, "italic(B)", color = "darkgrey", offset = 0.00315, align = TRUE,
offset.text = 0.0002, barsize = 2, fontsize = 5, parse = TRUE) +
geom_cladelabel(163, "italic(C)", color = "darkgrey", offset = .0031, align = TRUE,
offset.text = 0.0002, barsize = 3.2, fontsize = 5, parse = TRUE) +
geom_cladelabel(173, "italic(D1)", color = cols[["D1"]], offset = .003, align = TRUE,
offset.text = -.001, barsize = 1.2, extend = c(0, 0.5), parse = TRUE) +
geom_cladelabel(176, "italic(D2)", color = cols[["D2"]], offset = .003, align = TRUE,
offset.text = -.001, barsize = 1.2, extend = c(0.5, 0), parse = TRUE) +
geom_cladelabel(172, "italic(D)", color = "darkgrey", offset = .00315, align = TRUE,
offset.text = 0.0002, barsize = 2, fontsize = 5, parse = TRUE)
## Optional
## display support values
p <- p + geom_nodelab(aes(subset = (node == 92), label = "*"),
color = "black", nudge_x = -.001, nudge_y = 1) +
geom_nodelab(aes(subset = (node == 155), label = "*"),
color = "black", nudge_x = -.0003, nudge_y = -1) +
geom_nodelab(aes(subset = (node == 158), label = "95/92/1.00"),
color = "black", nudge_x = -0.0001, nudge_y = -1, hjust = 1) +
geom_nodelab(aes(subset = (node == 162), label = "98/97/1.00"),
color = "black", nudge_x = -0.0001, nudge_y = -1, hjust = 1) +
geom_nodelab(aes(subset = (node == 172), label = "*"),
color = "black", nudge_x = -.0003, nudge_y = -1)
```
```{r eval=F}
## extract accession numbers from tip labels
tl <- tree$tip.label
acc <- sub("\\w+\\|", "", tl)
names(tl) <- acc
## read sequences from GenBank directly into R
## and convert the object to DNAStringSet
tipseq <- ape::read.GenBank(acc) %>% as.character %>%
lapply(., paste0, collapse = "") %>% unlist %>%
DNAStringSet
## align the sequences using muscle
tipseq_aln <- muscle::muscle(tipseq)
tipseq_aln <- DNAStringSet(tipseq_aln)
```
```{r echo=F}
## extract accession numbers from tip labels
tl <- tree$tip.label
acc <- sub("\\w+\\|", "", tl)
names(tl) <- acc
## writeXStringSet(tipseq_aln, file = "data/HPV58_aln.fas")
tipseq_aln <- readDNAStringSet("data/HPV58_aln.fas")
```
(ref:jv2017scap) Phylogeny of HPV58 complete genomes with dot and line plots of pairwise nucleotide sequence distances.
(ref:jv2017cap) **Phylogeny of HPV58 complete genomes with dot and line plots of pairwise nucleotide sequence distances**.
```{r jv2017, fig.width=12, fig.height=12, , fig.cap="(ref:jv2017cap)", fig.scap="(ref:jv2017scap)", warning=FALSE}
## calculate pairwise hamming distances among sequences
tipseq_dist <- stringDist(tipseq_aln, method = "hamming")
## calculate percentage of differences
tipseq_d <- as.matrix(tipseq_dist) / width(tipseq_aln[1]) * 100
## convert the matrix to tidy data frame for facet_plot
dd <- as_tibble(tipseq_d)
dd$seq1 <- rownames(tipseq_d)
td <- gather(dd,seq2, dist, -seq1)
td$seq1 <- tl[td$seq1]
td$seq2 <- tl[td$seq2]
g <- p$data$group
names(g) <- p$data$label
td$clade <- g[td$seq2]
## visualize the sequence differences using dot plot and line plot
## and align the sequence difference plot to the tree using facet_plot
p2 <- facet_plot(p, panel = "Sequence Distance", data = td, geom_point,
mapping = aes(x = dist, color = clade, shape = clade), alpha = .6) %>%
facet_plot(panel = "Sequence Distance", data = td, geom = geom_path,
mapping=aes(x = dist, group = seq2, color = clade), alpha = .6) +
scale_shape_manual(values = 1:8, guide = FALSE)
print(p2)
```
## Displaying different symbolic points for bootstrap values.
We can cut the bootstrap values into several intervals, *e.g.* to indicate whether the clade is high, moderate or low support. Then we can use these intervals as categorical variable to set different color or shape of symbolic points to indicate the bootstrap values belong to which category.
(ref:bpintervalscap) Partitioning bootstrap values.
(ref:bpintervalcap) **Partitioning bootstrap values**. Bootstrap values were divided into three categories and this information was used to color circle points.
```{r include = FALSE}
## phytools also have a read.newick function
read.newick <- treeio::read.newick
```
```{r bpinterval, fig.width=7.5, fig.height=8.6, fig.cap="(ref:bpintervalcap)", fig.scap="(ref:bpintervalscap)", dev="jpeg"}
library(treeio)
library(ggplot2)
library(ggtree)
tree <- read.newick("data/RMI.phy_phyml_tree_rooted_labeled.txt", node.label='support')
root <- rootnode(tree)
ggtree(tree, color="black", size=1.5, linetype=1, right=TRUE) +
geom_tiplab(size=4.5, hjust = -0.060, fontface="bold") + xlim(0, 0.09) +
geom_point2(aes(subset=!isTip & node != root,
fill=cut(support, c(0, 700, 900, 1000))),
shape=21, size=4) +
theme_tree(legend.position=c(0.2, 0.2)) +
scale_fill_manual(values=c("white", "grey", "black"), guide='legend',
name='Bootstrap Percentage(BP)',
breaks=c('(900,1e+03]', '(700,900]', '(0,700]'),
labels=expression(BP>=90,70 <= BP * " < 90", BP < 70))
```
## Phylogenetic tree with genome locus structure {#genome-locus}
The `geom_motif` is defined in `ggtree` and it is a wrapper layer of `gggenes::geom_gene_arrow`. The `geom_motif` can automatically adjust genomic alignment by selective gene (via the `on` parameter) and can label genes via the `label` parameter. In the following example, we use `example_genes` dataset provided by `r CRANpkg("gggenes")`. As the dataset only provide genomic coordinations of a set of genes, a phylogeny for the genomes need to be firstly constructed. We calculate jaccard similarity based on the ratio of overlapping genes among genomes and correspondingly determine genome distance. BioNJ algorithm was applied to construct the tree (Figure \@ref(fig:gggenes)).
(ref:gggenesscap) Genomic features with phylogenetic tree.
(ref:gggenescap) **Genomic features with phylogenetic tree.**
```{r gggenes, fig.width=9, fig.height=4, fig.cap="(ref:gggenescap)", fig.scap="(ref:gggenesscap)"}
library(dplyr)
library(ggplot2)
library(gggenes)
library(ggtree)
get_genes <- function(data, genome) {
filter(data, molecule == genome) %>% pull(gene)
}
g <- unique(example_genes[,1])
n <- length(g)
d <- matrix(nrow = n, ncol = n)
rownames(d) <- colnames(d) <- g
genes <- lapply(g, get_genes, data = example_genes)
for (i in 1:n) {
for (j in 1:i) {
jaccard_sim <- length(intersect(genes[[i]], genes[[j]])) /
length(union(genes[[i]], genes[[j]]))
d[j, i] <- d[i, j] <- 1 - jaccard_sim
}
}
tree <- ape::bionj(d)
p <- ggtree(tree, branch.length='none') +
geom_tiplab() + xlim_tree(5.5) +
geom_facet(mapping = aes(xmin = start, xmax = end, fill = gene),
data = example_genes, geom = geom_motif, panel = 'Alignment',
on = 'genE', label = 'gene', align = 'left') +
scale_fill_brewer(palette = "Set3") +
scale_x_continuous(expand=c(0,0)) +
theme(strip.text=element_blank(),
panel.spacing=unit(0, 'cm'))
facet_widths(p, widths=c(1,2))
```