Skip to content

Commit 016d370

Browse files
ahmigahmig
committed
Decouple distance tree and time signal analysis plots from data
1 parent c4bc02a commit 016d370

6 files changed

Lines changed: 305 additions & 225 deletions

File tree

‎workflow/rules/report.smk‎

Lines changed: 52 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -157,46 +157,70 @@ rule context_phylogeny_plot:
157157
"../scripts/report/context_phylogeny_plot.R"
158158

159159

160-
# rule allele_freq_phylogeny_data: # nj
161-
# ...
162-
163-
164-
# rule allele_freq_phylogeny_plot:
165-
# ...
160+
rule allele_freq_tree_data:
161+
conda: "../envs/renv.yaml"
162+
params:
163+
use_bionj = config["USE_BIONJ"],
164+
ref_name = config["ALIGNMENT_REFERENCE"],
165+
input:
166+
dist = REPORT_DIR_TABLES/"distances.csv",
167+
output:
168+
tree = report(REPORT_DIR_TABLES/"allele_freq_tree.nwk"),
169+
log:
170+
LOGDIR / "allele_freq_tree_data" / "log.txt"
171+
script:
172+
"../scripts/report/allele_freq_tree_data.R"
166173

167174

168-
# rule time_signal_data: # tempest
169-
# ...
175+
rule allele_freq_tree_plot:
176+
conda: "../envs/renv.yaml"
177+
params:
178+
design = config["PLOTS"],
179+
ref_name = config["ALIGNMENT_REFERENCE"],
180+
plot_height_mm = 119.4,
181+
plot_width_mm = 159.2,
182+
input:
183+
tree = report(REPORT_DIR_TABLES/"allele_freq_tree.nwk"),
184+
study_fasta = OUTDIR/f"{OUTPUT_NAME}.fasta",
185+
metadata = config["METADATA"],
186+
output:
187+
plot = report(REPORT_DIR_PLOTS/"allele_freq_tree.png"),
188+
log:
189+
LOGDIR / "allele_freq_tree_plot" / "log.txt"
190+
script:
191+
"../scripts/report/allele_freq_tree_plot.R"
170192

171193

172-
# rule time_signal_plot: # tempest
173-
# ...
194+
rule time_signal_data:
195+
conda: "../envs/renv.yaml"
196+
params:
197+
ref_name = config["ALIGNMENT_REFERENCE"],
198+
input:
199+
tree = report(REPORT_DIR_TABLES/"allele_freq_tree.nwk"),
200+
metadata = config["METADATA"],
201+
output:
202+
table = report(REPORT_DIR_TABLES/"time_signal.csv"),
203+
json = REPORT_DIR_TABLES/"time_signal.json"
204+
log:
205+
LOGDIR / "time_signal_data" / "log.txt"
206+
script:
207+
"../scripts/report/time_signal_data.R"
174208

175209

176-
rule phylo_plots:
210+
rule time_signal_plot:
177211
conda: "../envs/renv.yaml"
178212
params:
179213
design = config["PLOTS"],
180-
ref_name = config["ALIGNMENT_REFERENCE"],
181-
boot_th = 95,
182-
alrt_th = 80,
183214
plot_height_mm = 119.4,
184215
plot_width_mm = 159.2,
185-
use_bionj = config["USE_BIONJ"]
186216
input:
187-
dist = REPORT_DIR_TABLES/f"distances.csv",
188-
study_fasta = OUTDIR/f"{OUTPUT_NAME}.fasta",
189-
ml = OUTDIR/f"tree_context/{OUTPUT_NAME}.treefile",
190-
metadata = config["METADATA"]
217+
table = report(REPORT_DIR_TABLES/"time_signal.csv"),
191218
output:
192-
temest = report(REPORT_DIR_PLOTS/"figure_9.png"),
193-
tree = report(REPORT_DIR_PLOTS/"figure_8.png"),
194-
table = report(REPORT_DIR_TABLES/"figure_9.csv"),
195-
json = temp(OUTDIR/"stats.lm.json")
219+
plot = report(REPORT_DIR_PLOTS/"time_signal.png"),
196220
log:
197-
LOGDIR / "phylo_plots" / "log.txt"
221+
LOGDIR / "time_signal_plot" / "log.txt"
198222
script:
199-
"../scripts/report/phylo_plots.R"
223+
"../scripts/report/time_signal_plot.R"
200224

201225

202226
rule evo_plots:
@@ -263,14 +287,14 @@ rule report:
263287
SNV_spike = report(REPORT_DIR_PLOTS/"figure_5b.png"),
264288
volcano = report(REPORT_DIR_PLOTS/"figure_6.png"),
265289
panel = report(REPORT_DIR_PLOTS/"figure_7.png"),
266-
tree = report(REPORT_DIR_PLOTS/"figure_8.png"),
267-
temest = report(REPORT_DIR_PLOTS/"figure_9.png"),
290+
tree = report(REPORT_DIR_PLOTS/"allele_freq_tree.png"),
291+
temest = report(REPORT_DIR_PLOTS/"time_signal.png"),
268292
heat_table = report(REPORT_DIR_TABLES/"heatmap.csv"),
269293
evo = report(REPORT_DIR_PLOTS/"figure_11.png"),
270294
omega_plot = report(REPORT_DIR_PLOTS/"figure_12.png"),
271295
freyja_ts = OUTDIR/"demixing"/"freyja_data"/"last_barcode_update.txt",
272296
value = REPORT_DIR_TABLES/"diversity.json",
273-
stats_lm = OUTDIR/"stats.lm.json",
297+
stats_lm = REPORT_DIR_TABLES/"time_signal.json",
274298
stats_ml = REPORT_DIR_TABLES/"context_phylogeny.json",
275299
table = OUTDIR/"summary_table.csv",
276300
sum_nv = OUTDIR/"summary_nv.json",

‎workflow/scripts/report/allele_freq_tree_data.R‎

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
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(data.table)
10+
library(ape)
11+
library(logger)
12+
13+
fix_negative_edge_length <- function(nj.tree) {
14+
edge_infos <- cbind(
15+
nj.tree$edge,
16+
nj.tree$edge.length
17+
) %>%
18+
as.data.table()
19+
colnames(edge_infos) <- c("from", "to", "length")
20+
nega_froms <- edge_infos[length < 0, sort(unique(from))]
21+
for (nega_from in nega_froms) {
22+
minus_length <- edge_infos[from == nega_from, ][order(length)][1, length]
23+
edge_infos[from == nega_from, length := length - minus_length]
24+
edge_infos[to == nega_from, length := length + minus_length]
25+
}
26+
nj.tree$edge.length <- edge_infos$length
27+
nj.tree
28+
}
29+
30+
log_threshold(INFO)
31+
32+
matrix <- read_csv(snakemake@input$dist)
33+
34+
# Tree construction
35+
if (snakemake@params$use_bionj) {
36+
log_info("Constructing BIONJ tree based on distances")
37+
tree_method <- bionj
38+
} else {
39+
log_info("Constructing NJ tree based on distances")
40+
tree_method <- nj
41+
}
42+
tree <- matrix %>%
43+
column_to_rownames(var = "...1") %>%
44+
as.matrix() %>%
45+
as.dist() %>%
46+
tree_method() %>%
47+
root(snakemake@params$ref_name, resolve.root = TRUE)
48+
49+
# Resolve possible negative edge lengths in tree
50+
log_info("Resolving negative edge lengths")
51+
tree <- fix_negative_edge_length(tree)
52+
53+
log_info("Saving tree")
54+
write.tree(tree, snakemake@output$tree)

‎workflow/scripts/report/allele_freq_tree_plot.R‎

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
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(ape)
10+
library(ggtree)
11+
library(logger)
12+
13+
log_threshold(INFO)
14+
15+
# Import file with plots style
16+
source(snakemake@params[["design"]])
17+
18+
log_info("Building tree annotation")
19+
20+
# Read metadata
21+
log_debug("Reading metadata")
22+
metadata <- read_csv(snakemake@input$metadata)
23+
24+
# Extract names from records
25+
log_debug("Reading FASTA")
26+
study_names <- read.dna(
27+
snakemake@input$study_fasta,
28+
format = "fasta",
29+
as.matrix = FALSE,
30+
)
31+
log_debug("Extracting names")
32+
study_names <- study_names[
33+
!startsWith(names(study_names), snakemake@params$ref_name)
34+
]
35+
study_names <- names(study_names)
36+
37+
# Obtain sample names ordered by CollectionDate
38+
log_debug("Sorting names by collection date")
39+
date_order <- metadata %>%
40+
arrange(CollectionDate) %>%
41+
filter(ID %in% study_names) %>%
42+
pull(ID) %>%
43+
unique()
44+
45+
log_debug("Building annotation dataframe")
46+
tree_tiplab <- data.frame(
47+
ID = date_order,
48+
order = seq(1, length(date_order), 1)
49+
) %>%
50+
rowwise() %>%
51+
mutate(
52+
tip_label = sprintf(
53+
"(%s)-%s",
54+
order,
55+
ID
56+
)
57+
) %>%
58+
ungroup() %>%
59+
add_row(
60+
ID = snakemake@params$ref_name,
61+
order = 0,
62+
tip_label = snakemake@params$ref_name
63+
)
64+
65+
# Read distance tree and root
66+
log_info("Reading tree")
67+
tree <- read.tree(snakemake@input$tree) %>%
68+
root(snakemake@params$ref_name, resolve.root = TRUE)
69+
70+
# Plot
71+
log_info("Plotting distance tree")
72+
max.tip.length <- max(node.depth.edgelength(tree)[1:length(tree$tip.label)])
73+
p <- ggtree(tree) %<+% tree_tiplab +
74+
geom_tiplab(aes(label = tip_label)) +
75+
geom_treescale(1.1 * max.tip.length) +
76+
geom_rootedge(0.05 * max.tip.length)
77+
78+
ggsave(
79+
filename = snakemake@output$plot,
80+
plot = p,
81+
width = snakemake@params$plot_width_mm,
82+
height = snakemake@params$plot_height_mm,
83+
units = "mm",
84+
dpi = 250
85+
)

0 commit comments

Comments
 (0)