Skip to content

Commit c87c87f

Browse files
ahmigahmig
committed
Decouple SNP AF trajectory panel plot from data
Updates color selection so discard low contrast colors
1 parent dd1e1ad commit c87c87f

5 files changed

Lines changed: 125 additions & 56 deletions

File tree

‎config/design_plots.R‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -101,3 +101,9 @@ DNDS_SHAPES <- c(
101101
dN = 2,
102102
dS = 4
103103
)
104+
105+
# Allele frequency trajectories panel color
106+
ALL.COLORS <- grDevices::colors()
107+
TRAJECTORY.PANEL.COLORS <- ALL.COLORS[
108+
!grepl("(gray|grey|white|snow|azure|beige)", ALL.COLORS)
109+
]

‎workflow/rules/report.smk‎

Lines changed: 17 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -240,7 +240,7 @@ rule dnds_plots:
240240
"../scripts/report/dnds_plots.R"
241241

242242

243-
rule af_time_correlation_per_snp_data:
243+
rule af_time_correlation_data:
244244
conda: "../envs/renv.yaml"
245245
params:
246246
cor_method = config["COR"]["METHOD"],
@@ -250,76 +250,44 @@ rule af_time_correlation_per_snp_data:
250250
metadata = config["METADATA"],
251251
output:
252252
fmt_variants = temp(OUTDIR/f"{OUTPUT_NAME}.variants.filled.dated.tsv"),
253-
correlations = report(REPORT_DIR_TABLES/"af_time_correlation_per_snp.csv"), # fig 6
253+
correlations = report(REPORT_DIR_TABLES/"af_time_correlation.csv"),
254+
subset = REPORT_DIR_TABLES/"af_time_correlation.subset.txt",
254255
log:
255-
LOGDIR / "af_time_correlation_per_snp_data" / "log.txt"
256+
LOGDIR / "af_time_correlation_data" / "log.txt"
256257
script:
257-
"../scripts/report/af_time_correlation_per_snp_data.R"
258+
"../scripts/report/af_time_correlation_data.R"
258259

259260

260-
rule af_time_correlation_per_snp_plot:
261+
rule af_time_correlation_plot:
261262
conda: "../envs/renv.yaml"
262263
params:
263264
design = config["PLOTS"],
264265
input:
265-
correlations = REPORT_DIR_TABLES/"af_time_correlation_per_snp.csv",
266+
correlations = REPORT_DIR_TABLES/"af_time_correlation.csv",
266267
output:
267-
plot = report(REPORT_DIR_PLOTS/"af_time_correlation_per_snp.png"), # fig 6
268+
plot = report(REPORT_DIR_PLOTS/"af_time_correlation.png"),
268269
log:
269-
LOGDIR / "af_time_correlation_per_snp_plot" / "log.txt"
270+
LOGDIR / "af_time_correlation_plot" / "log.txt"
270271
script:
271-
"../scripts/report/af_time_correlation_per_snp_plot.R"
272+
"../scripts/report/af_time_correlation_plot.R"
272273

273274

274-
rule af_trajectory_panel_data: # panel with AF trajectories in time
275-
# TODO
276-
conda: "../envs/renv.yaml"
277-
input:
278-
fmt_variants = OUTDIR/f"{OUTPUT_NAME}.variants.filled.dated.tsv",
279-
correlations = REPORT_DIR_TABLES/"af_time_correlation_per_snp.csv",
280-
output:
281-
table = report(REPORT_DIR_TABLES/"af_trajectory_panel.csv"), # fig 7
282-
log:
283-
LOGDIR / "af_trajectory_panel_data" / "log.txt"
284-
script:
285-
"../scripts/report/af_trajectory_panel_data.R"
286-
287-
288-
rule af_trajectory_panel_plot: # panel with AF trajectories in time
289-
# TODO
275+
rule af_trajectory_panel_plot:
290276
conda: "../envs/renv.yaml"
291277
params:
292278
design = config["PLOTS"],
279+
random_color_seed = 7291,
293280
input:
294-
table = REPORT_DIR_TABLES/"af_trajectory_panel.csv",
281+
fmt_variants = OUTDIR/f"{OUTPUT_NAME}.variants.filled.dated.tsv",
282+
subset = REPORT_DIR_TABLES/"af_time_correlation.subset.txt"
295283
output:
296-
plot = report(REPORT_DIR_PLOTS/"af_trajectory_panel.png"), # fig 7
284+
plot = report(REPORT_DIR_PLOTS/"af_trajectory_panel.png"),
297285
log:
298286
LOGDIR / "af_trajectory_panel_plot" / "log.txt"
299287
script:
300288
"../scripts/report/af_trajectory_panel_plot.R"
301289

302290

303-
# rule snp_plots:
304-
# conda: "../envs/renv.yaml"
305-
# params:
306-
# design = config["PLOTS"],
307-
# cor_method = config["COR"]["METHOD"],
308-
# cor_exact = config["COR"]["EXACT"]
309-
# input:
310-
# vcf = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
311-
# metadata = config["METADATA"]
312-
# output:
313-
# pseudovolcano = report(REPORT_DIR_PLOTS/"figure_6.png"),
314-
# snp_panel = report(REPORT_DIR_PLOTS/"figure_7.png"),
315-
# table_1 = report(REPORT_DIR_TABLES/"figure_6.csv"),
316-
# table_2 = report(REPORT_DIR_TABLES/"figure_7.csv")
317-
# log:
318-
# LOGDIR / "snp_plots" / "log.txt"
319-
# script:
320-
# "../scripts/report/snp_plots.R"
321-
322-
323291
rule summary_table:
324292
conda: "../envs/renv.yaml"
325293
input:
@@ -344,8 +312,8 @@ rule report:
344312
fig_cor = report(REPORT_DIR_PLOTS/"figure_4.png"),
345313
SNV = report(REPORT_DIR_PLOTS/"figure_5a.png"),
346314
SNV_spike = report(REPORT_DIR_PLOTS/"figure_5b.png"),
347-
volcano = report(REPORT_DIR_TABLES/"af_time_correlation_per_snp.png"),
348-
panel = report(REPORT_DIR_TABLES/"af_trajectory_panel.png"),
315+
volcano = report(REPORT_DIR_PLOTS/"af_time_correlation.png"),
316+
panel = report(REPORT_DIR_PLOTS/"af_trajectory_panel.png"),
349317
tree = report(REPORT_DIR_PLOTS/"allele_freq_tree.png"),
350318
temest = report(REPORT_DIR_PLOTS/"time_signal.png"),
351319
heat_table = report(REPORT_DIR_TABLES/"heatmap.csv"),

workflow/scripts/report/af_time_correlation_per_snp_data.R renamed to workflow/scripts/report/af_time_correlation_data.R

Lines changed: 33 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,7 @@ variants <- left_join(variants, metadata, by = c("SAMPLE" = "ID")) %>%
4545
CollectionDate
4646
) %>%
4747
mutate(
48-
interval = as.numeric(CollectionDate - min(CollectionDate))
48+
interval = as.numeric(difftime(CollectionDate, min(CollectionDate), units = "days"))
4949
)
5050

5151
# Save processed input
@@ -55,11 +55,7 @@ write_csv(variants, snakemake@output$fmt_variants)
5555
log_info("Calculating correlations")
5656
log_debug("Calculating unique SNPs")
5757
# Get list with all different polymorphisms
58-
unique.snps <- pull(
59-
variants,
60-
VARIANT_NAME
61-
) %>%
62-
unique()
58+
unique.snps <- unique(variants$VARIANT_NAME)
6359

6460
# Create an empty dataframe to be filled
6561
cor.df <- data.frame(
@@ -109,3 +105,34 @@ write_csv(
109105
correlations,
110106
snakemake@output[["correlations"]]
111107
)
108+
109+
log_info("Selecting variants whose allele frequency is significantly correlated with time")
110+
significant.variants <- correlations %>%
111+
filter(p.value.adj < 0.05) %>%
112+
pull(variant) %>%
113+
unique()
114+
115+
log_info("Significant: {significant.variants}")
116+
117+
log_info("Selecting variants in positions with more than one alternative allele")
118+
mult.alt.variants <- variants %>%
119+
select(
120+
VARIANT_NAME,
121+
POS
122+
) %>%
123+
distinct() %>%
124+
group_by(POS) %>%
125+
filter(n() > 1) %>%
126+
ungroup() %>%
127+
pull(VARIANT_NAME) %>%
128+
unique()
129+
130+
log_info("Mult all: {mult.alt.variants}")
131+
132+
# Build selected subset to represent
133+
variant.selection <- unique(c(significant.variants, mult.alt.variants))
134+
135+
log_info("Selection: {variant.selection}")
136+
137+
log_info("Writing selected variants subset")
138+
write_lines(variant.selection, snakemake@output$subset)

workflow/scripts/report/af_time_correlation_per_snp_plot.R renamed to workflow/scripts/report/af_time_correlation_plot.R

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,7 @@ p <- correlations %>%
3333
geom_hline(
3434
aes(yintercept = -log10(0.05)),
3535
linetype = 2,
36-
color = "orange"
36+
color = "red"
3737
) +
3838
labs(
3939
x = "Correlation coefficient",

‎workflow/scripts/report/af_trajectory_panel_plot.R‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
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+
set.seed(snakemake@params$random_color_seed) # seed for sampling colors
9+
10+
library(tidyverse)
11+
library(logger)
12+
log_threshold(INFO)
13+
14+
# Import file with plots style
15+
source(snakemake@params[["design"]])
16+
17+
log_info("Reading formatted variants table")
18+
variants <- read_csv(snakemake@input$fmt_variants)
19+
20+
log_info("Reading subset of variants to represent")
21+
selected.variants <- read_lines(snakemake@input$subset)
22+
23+
# Set plot height depending on the number of SNPs assuming 4 columns in the plot
24+
log_debug("Calculating plot height")
25+
plot.height <- ceiling(length(selected.variants) / 4) * 42
26+
27+
log_info("Plotting {length(selected.variants)} SNPs allele frequency trajectories in time")
28+
selected.colors <- sample(TRAJECTORY.PANEL.COLORS, length(selected.variants))
29+
log_debug("Selected color: {selected.colors}")
30+
p <- variants %>%
31+
filter(VARIANT_NAME %in% selected.variants) %>%
32+
ggplot() +
33+
aes(
34+
x = interval,
35+
y = ALT_FREQ,
36+
color = VARIANT_NAME
37+
) +
38+
scale_color_manual(values = selected.colors) +
39+
geom_point() +
40+
geom_line() +
41+
theme(
42+
legend.position = "bottom",
43+
legend.text = element_text(size = 9)
44+
) +
45+
labs(
46+
x = "Days since first sample",
47+
y = "Frequency",
48+
color = "NV"
49+
) +
50+
guides(color = guide_legend(ncol = 3))
51+
52+
if (length(selected.variants) > 1) {
53+
p <- p +
54+
facet_wrap(
55+
vars(POS),
56+
nrow = ceiling(length(selected.variants) / 4),
57+
ncol = 4
58+
)
59+
}
60+
61+
ggsave(
62+
filename = snakemake@output[["plot"]],
63+
plot = p,
64+
width = 159.2,
65+
height = max(100, plot.height),
66+
units = "mm",
67+
dpi = 250
68+
)

0 commit comments

Comments
 (0)