Skip to content

Commit c4bc02a

Browse files
ahmigahmig
committed
Decouple context phylogeny plot from data
1 parent 3e18cd6 commit c4bc02a

7 files changed

Lines changed: 242 additions & 102 deletions

File tree

‎config/design_plots.R‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -58,6 +58,17 @@ TREE_NODE_SIZE <- c(
5858
boot_alrt_pass = 0.8
5959
)
6060

61+
TREE_LEGEND_NAMES <- c(
62+
tip_label = "Target samples",
63+
boot_alrt_pass = sprintf(
64+
"UFBoot ≥ %s%s & SH-aLRT ≥ %s%s ",
65+
snakemake@params[["boot_th"]],
66+
"%",
67+
snakemake@params[["alrt_th"]],
68+
"%"
69+
)
70+
)
71+
6172
# Nucleotide variants classification colors and labels
6273
NV_TYPE_PALETTE <- c(
6374
Frameshift = "#568D63",

‎template.qmd‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ params:
3737
tree_ml: ""
3838
fig_cor_snp: ""
3939
stats_lm: ""
40+
stats_ml: ""
4041
table: ""
4142
sum_nv: ""
4243
heat_tab: ""
@@ -63,8 +64,11 @@ library(readr)
6364
# Diversity
6465
div_values <- fromJSON(params$div_value)
6566
66-
# Temporal signal
67-
stats <- fromJSON(params$stats_lm)
67+
# Temporal signal and context tree
68+
stats <- append(
69+
fromJSON(params$stats_lm),
70+
fromJSON(params$stats_ml)
71+
)
6872
correlation <- stats[["r2"]]
6973
sub_rate <- stats[["sub_rate"]]
7074
sub_rate <- round(sub_rate, digits = 2)

‎workflow/rules/distances.smk‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ rule weighted_distances:
1010
ancestor = OUTDIR/f"{OUTPUT_NAME}.ancestor.fasta",
1111
reference = OUTDIR/"reference.fasta"
1212
output:
13-
distances = REPORT_DIR_TABLES/f"figure_8.csv"
13+
distances = REPORT_DIR_TABLES/f"distances.csv"
1414
log:
1515
LOGDIR / "weighted_distances" / "log.txt"
1616
script:

‎workflow/rules/report.smk‎

Lines changed: 57 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -120,6 +120,59 @@ rule general_NV_description:
120120
"../scripts/report/NV_description.R"
121121

122122

123+
rule context_phylogeny_data:
124+
conda: "../envs/renv.yaml"
125+
params:
126+
design = config["PLOTS"],
127+
ref_name = config["ALIGNMENT_REFERENCE"],
128+
boot_th = 95.0,
129+
alrt_th = 80.0,
130+
input:
131+
target_fasta = OUTDIR/f"{OUTPUT_NAME}.fasta",
132+
tree = OUTDIR/f"tree_context/{OUTPUT_NAME}.treefile",
133+
output:
134+
json = REPORT_DIR_TABLES/"context_phylogeny.json",
135+
annotation = REPORT_DIR_TABLES/"context_phylogeny.csv",
136+
log:
137+
LOGDIR / "context_phylogeny_data" / "log.txt"
138+
script:
139+
"../scripts/report/context_phylogeny_data.R"
140+
141+
142+
rule context_phylogeny_plot:
143+
conda: "../envs/renv.yaml"
144+
params:
145+
design = config["PLOTS"],
146+
plot_height_mm = 119.4,
147+
plot_width_mm = 159.2,
148+
input:
149+
tree = OUTDIR/f"tree_context/{OUTPUT_NAME}.treefile",
150+
json = REPORT_DIR_TABLES/"context_phylogeny.json",
151+
annotation = REPORT_DIR_TABLES/"context_phylogeny.csv"
152+
output:
153+
plot = report(REPORT_DIR_PLOTS/"context_phylogeny.png"),
154+
log:
155+
LOGDIR / "context_phylogeny_plot" / "log.txt"
156+
script:
157+
"../scripts/report/context_phylogeny_plot.R"
158+
159+
160+
# rule allele_freq_phylogeny_data: # nj
161+
# ...
162+
163+
164+
# rule allele_freq_phylogeny_plot:
165+
# ...
166+
167+
168+
# rule time_signal_data: # tempest
169+
# ...
170+
171+
172+
# rule time_signal_plot: # tempest
173+
# ...
174+
175+
123176
rule phylo_plots:
124177
conda: "../envs/renv.yaml"
125178
params:
@@ -131,14 +184,13 @@ rule phylo_plots:
131184
plot_width_mm = 159.2,
132185
use_bionj = config["USE_BIONJ"]
133186
input:
134-
dist = REPORT_DIR_TABLES/f"figure_8.csv",
187+
dist = REPORT_DIR_TABLES/f"distances.csv",
135188
study_fasta = OUTDIR/f"{OUTPUT_NAME}.fasta",
136189
ml = OUTDIR/f"tree_context/{OUTPUT_NAME}.treefile",
137190
metadata = config["METADATA"]
138191
output:
139192
temest = report(REPORT_DIR_PLOTS/"figure_9.png"),
140193
tree = report(REPORT_DIR_PLOTS/"figure_8.png"),
141-
tree_ml = report(REPORT_DIR_PLOTS/"figure_2.png"),
142194
table = report(REPORT_DIR_TABLES/"figure_9.csv"),
143195
json = temp(OUTDIR/"stats.lm.json")
144196
log:
@@ -204,7 +256,7 @@ rule report:
204256
input:
205257
qmd = Path(config["REPORT_QMD"]).resolve(),
206258
demix = report(REPORT_DIR_PLOTS/"demix.png"),
207-
tree_ml = report(REPORT_DIR_PLOTS/"figure_2.png"),
259+
tree_ml = report(REPORT_DIR_PLOTS/"context_phylogeny.png"),
208260
diversity = report(REPORT_DIR_PLOTS/"diversity.png"),
209261
fig_cor = report(REPORT_DIR_PLOTS/"figure_4.png"),
210262
SNV = report(REPORT_DIR_PLOTS/"figure_5a.png"),
@@ -219,6 +271,7 @@ rule report:
219271
freyja_ts = OUTDIR/"demixing"/"freyja_data"/"last_barcode_update.txt",
220272
value = REPORT_DIR_TABLES/"diversity.json",
221273
stats_lm = OUTDIR/"stats.lm.json",
274+
stats_ml = REPORT_DIR_TABLES/"context_phylogeny.json",
222275
table = OUTDIR/"summary_table.csv",
223276
sum_nv = OUTDIR/"summary_nv.json",
224277
params:
@@ -258,6 +311,7 @@ rule report:
258311
"tree_ml = '{input.tree_ml}', "
259312
"fig_cor_snp = '{input.fig_cor}', "
260313
"stats_lm = '{input.stats_lm}', "
314+
"stats_ml = '{input.stats_ml}', "
261315
"table = '{input.table}', "
262316
"sum_nv = '{input.sum_nv}', "
263317
"heat_tab = '{input.heat_table}', "

‎workflow/scripts/report/context_phylogeny_data.R‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
#!/usr/bin/env Rscript
2+
3+
# Write stdout and stderr to log file
4+
log <- file(snakemake@log[[1]], open = "wt")
5+
sink(log, type = "message")
6+
sink(log, type = "output")
7+
8+
library(tidyverse)
9+
library(jsonlite)
10+
library(ape)
11+
library(ggtree)
12+
library(logger)
13+
14+
log_threshold(INFO)
15+
16+
# Import file with plots style
17+
source(snakemake@params[["design"]])
18+
19+
tree_ml <- read.tree(snakemake@input[["tree"]]) %>%
20+
root(
21+
snakemake@params[["ref_name"]],
22+
resolve.root = TRUE
23+
)
24+
25+
target_names <- read.dna(
26+
snakemake@input[["target_fasta"]],
27+
format = "fasta",
28+
as.matrix = FALSE,
29+
)
30+
31+
target_names <- target_names[
32+
!startsWith(names(target_names), snakemake@config[["ALIGNMENT_REFERENCE"]])
33+
]
34+
target_names <- names(target_names)
35+
36+
# ML tree with context data
37+
# Internal nodes color
38+
# Node labels contain SH-aLRT/UFboot values
39+
aLRT.values <- sapply(
40+
strsplit(tree_ml$node.label, "/"),
41+
function(x) {
42+
as.numeric(x[1])
43+
}
44+
)
45+
46+
bootstrap.values <- sapply(
47+
strsplit(tree_ml$node.label, "/"),
48+
function(x) {
49+
as.numeric(x[2])
50+
}
51+
)
52+
53+
aLRT.mask <- aLRT.values >= snakemake@params[["alrt_th"]]
54+
boot.mask <- bootstrap.values >= snakemake@params[["boot_th"]]
55+
56+
# MRCA for target samples
57+
log_info("Calculating MRCA of target samples")
58+
target.mrca <- getMRCA(tree_ml, target_names)
59+
target.mrca.clade <- extract.clade(tree_ml, target.mrca)
60+
target.mrca.clade.ntips <- Ntip(target.mrca.clade)
61+
62+
log_info("Building M-L tree annotation")
63+
tree_ml.ntips <- length(tree_ml$tip.label)
64+
tree_ml.nodes <- tree_ml$Nnode + tree_ml.ntips
65+
tree_ml.labels <- tree_ml$tip.label[1:tree_ml.nodes]
66+
tree_ml.node.pass <- c(rep(FALSE, tree_ml.ntips), aLRT.mask & boot.mask)
67+
68+
ml.tree.annot <- tibble(
69+
node = 1:tree_ml.nodes,
70+
) %>%
71+
mutate(
72+
Class = case_when(
73+
tree_ml.labels %in% target_names ~ TREE_LEGEND_NAMES["tip_label"],
74+
tree_ml.node.pass ~ TREE_LEGEND_NAMES["boot_alrt_pass"],
75+
TRUE ~ NA
76+
)
77+
)
78+
79+
# Write output files
80+
log_info("Writing tree annotation")
81+
write_csv(ml.tree.annot, snakemake@output$annotation)
82+
83+
log_info("Writing JSON data")
84+
target.node <- tree_ml$node.label[target.mrca - length(tree_ml$tip.label)]
85+
list(
86+
"boot" = strsplit(target.node, "/")[[1]][2] %>% as.numeric(),
87+
"alrt" = strsplit(target.node, "/")[[1]][1] %>% as.numeric(),
88+
"monophyly" = ifelse(
89+
is.monophyletic(tree_ml, target_names),
90+
"are",
91+
"are not"
92+
),
93+
"target_mrca" = target.mrca,
94+
"clade_tips" = target.mrca.clade.ntips,
95+
"max_tip_length" = max(node.depth.edgelength(tree_ml)[
96+
1:length(tree_ml$tip.label)
97+
]),
98+
"root" = snakemake@params[["ref_name"]]
99+
) %>%
100+
write_json(
101+
snakemake@output$json,
102+
auto_unbox = TRUE,
103+
digits = NA
104+
)

‎workflow/scripts/report/context_phylogeny_plot.R‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
#!/usr/bin/env Rscript
2+
3+
# Write stdout and stderr to log file
4+
log <- file(snakemake@log[[1]], open = "wt")
5+
sink(log, type = "message")
6+
sink(log, type = "output")
7+
8+
library(tidyverse)
9+
library(jsonlite)
10+
library(ape)
11+
library(ggtree)
12+
library(logger)
13+
14+
log_threshold(INFO)
15+
16+
# Import file with plots style
17+
source(snakemake@params[["design"]])
18+
19+
log_info("Reading JSON data")
20+
json <- read_json(snakemake@input$json)
21+
22+
log_info("Reading tree")
23+
tree_ml <- read.tree(snakemake@input[["tree"]]) %>%
24+
root(
25+
json$root,
26+
resolve.root = TRUE
27+
)
28+
29+
log_info("Reading tree annotation")
30+
annotation <- read_csv(snakemake@input$annotation)
31+
32+
log_info("Plotting M-L tree with context samples")
33+
p <- ggtree(tree_ml, layout = "circular") %<+% annotation +
34+
geom_highlight(node = json$target_mrca, colour = "red", alpha = 0) +
35+
geom_point(aes(color = Class, size = Class)) +
36+
geom_treescale(1.05 * json$max_tip_length) +
37+
geom_rootedge(0.05 * json$max_tip_length) +
38+
xlim(-json$max_tip_length / 3, NA) +
39+
scale_color_manual(
40+
values = setNames(
41+
TREE_PALETTE[names(TREE_LEGEND_NAMES)],
42+
TREE_LEGEND_NAMES
43+
),
44+
na.translate = FALSE
45+
) +
46+
scale_size_manual(
47+
values = setNames(
48+
TREE_NODE_SIZE[names(TREE_LEGEND_NAMES)],
49+
TREE_LEGEND_NAMES
50+
),
51+
na.translate = FALSE
52+
)
53+
54+
ggsave(
55+
filename = snakemake@output[["plot"]],
56+
plot = p,
57+
width = snakemake@params[["plot_width_mm"]],
58+
height = snakemake@params[["plot_height_mm"]],
59+
units = "mm",
60+
dpi = 250
61+
)

0 commit comments

Comments
 (0)