Skip to content

Commit dd1e1ad

Browse files
ahmigahmig
committed
Decouple SNP correlation (pseudovolcano) plot from data
1 parent 2b420c9 commit dd1e1ad

3 files changed

Lines changed: 234 additions & 13 deletions

File tree

‎workflow/rules/report.smk‎

Lines changed: 73 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -240,24 +240,84 @@ rule dnds_plots:
240240
"../scripts/report/dnds_plots.R"
241241

242242

243-
rule snp_plots:
243+
rule af_time_correlation_per_snp_data:
244244
conda: "../envs/renv.yaml"
245245
params:
246-
design = config["PLOTS"],
247246
cor_method = config["COR"]["METHOD"],
248-
cor_exact = config["COR"]["EXACT"]
247+
cor_exact = config["COR"]["EXACT"],
249248
input:
250-
vcf = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
251-
metadata = config["METADATA"]
249+
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
250+
metadata = config["METADATA"],
251+
output:
252+
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
254+
log:
255+
LOGDIR / "af_time_correlation_per_snp_data" / "log.txt"
256+
script:
257+
"../scripts/report/af_time_correlation_per_snp_data.R"
258+
259+
260+
rule af_time_correlation_per_snp_plot:
261+
conda: "../envs/renv.yaml"
262+
params:
263+
design = config["PLOTS"],
264+
input:
265+
correlations = REPORT_DIR_TABLES/"af_time_correlation_per_snp.csv",
266+
output:
267+
plot = report(REPORT_DIR_PLOTS/"af_time_correlation_per_snp.png"), # fig 6
268+
log:
269+
LOGDIR / "af_time_correlation_per_snp_plot" / "log.txt"
270+
script:
271+
"../scripts/report/af_time_correlation_per_snp_plot.R"
272+
273+
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
290+
conda: "../envs/renv.yaml"
291+
params:
292+
design = config["PLOTS"],
293+
input:
294+
table = REPORT_DIR_TABLES/"af_trajectory_panel.csv",
252295
output:
253-
pseudovolcano = report(REPORT_DIR_PLOTS/"figure_6.png"),
254-
snp_panel = report(REPORT_DIR_PLOTS/"figure_7.png"),
255-
table_1 = report(REPORT_DIR_TABLES/"figure_6.csv"),
256-
table_2 = report(REPORT_DIR_TABLES/"figure_7.csv")
296+
plot = report(REPORT_DIR_PLOTS/"af_trajectory_panel.png"), # fig 7
257297
log:
258-
LOGDIR / "snp_plots" / "log.txt"
298+
LOGDIR / "af_trajectory_panel_plot" / "log.txt"
259299
script:
260-
"../scripts/report/snp_plots.R"
300+
"../scripts/report/af_trajectory_panel_plot.R"
301+
302+
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"
261321

262322

263323
rule summary_table:
@@ -284,8 +344,8 @@ rule report:
284344
fig_cor = report(REPORT_DIR_PLOTS/"figure_4.png"),
285345
SNV = report(REPORT_DIR_PLOTS/"figure_5a.png"),
286346
SNV_spike = report(REPORT_DIR_PLOTS/"figure_5b.png"),
287-
volcano = report(REPORT_DIR_PLOTS/"figure_6.png"),
288-
panel = report(REPORT_DIR_PLOTS/"figure_7.png"),
347+
volcano = report(REPORT_DIR_TABLES/"af_time_correlation_per_snp.png"),
348+
panel = report(REPORT_DIR_TABLES/"af_trajectory_panel.png"),
289349
tree = report(REPORT_DIR_PLOTS/"allele_freq_tree.png"),
290350
temest = report(REPORT_DIR_PLOTS/"time_signal.png"),
291351
heat_table = report(REPORT_DIR_TABLES/"heatmap.csv"),

‎workflow/scripts/report/af_time_correlation_per_snp_data.R‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
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(logger)
10+
log_threshold(INFO)
11+
12+
# Read data
13+
log_info("Reading variants table")
14+
variants <- read_delim(
15+
snakemake@input[["variants"]],
16+
col_select = c(
17+
"VARIANT_NAME",
18+
"SAMPLE",
19+
"REGION",
20+
"ALT_FREQ",
21+
"POS"
22+
)
23+
) %>%
24+
# Fill positions without alt frequency with 0
25+
complete(
26+
nesting(REGION, VARIANT_NAME, POS),
27+
SAMPLE,
28+
fill = list(ALT_FREQ = 0)
29+
)
30+
31+
log_info("Reading metadata")
32+
metadata <- read_csv(snakemake@input[["metadata"]]) %>%
33+
filter(
34+
ID %in% variants$SAMPLE
35+
) %>%
36+
select(
37+
ID,
38+
CollectionDate
39+
)
40+
41+
log_info("Dating variants")
42+
variants <- left_join(variants, metadata, by = c("SAMPLE" = "ID")) %>%
43+
# Calculate days since first sample
44+
arrange(
45+
CollectionDate
46+
) %>%
47+
mutate(
48+
interval = as.numeric(CollectionDate - min(CollectionDate))
49+
)
50+
51+
# Save processed input
52+
log_info("Writing dated and frequency-filled variants")
53+
write_csv(variants, snakemake@output$fmt_variants)
54+
55+
log_info("Calculating correlations")
56+
log_debug("Calculating unique SNPs")
57+
# Get list with all different polymorphisms
58+
unique.snps <- pull(
59+
variants,
60+
VARIANT_NAME
61+
) %>%
62+
unique()
63+
64+
# Create an empty dataframe to be filled
65+
cor.df <- data.frame(
66+
variant = "",
67+
coefficient = 0,
68+
p.value = 0,
69+
p.value.adj = 0
70+
) %>%
71+
filter(p.value != 0)
72+
73+
log_debug("Calculating correlation using method = {snakemake@params$cor_method} and exact p-value = {snakemake@params$cor_exact}")
74+
correlations <- lapply(
75+
unique.snps,
76+
function(snp) {
77+
# Select SNP
78+
df <- filter(
79+
variants,
80+
VARIANT_NAME == snp
81+
)
82+
# Perform calculation
83+
test <- cor.test(
84+
df$ALT_FREQ,
85+
df$interval,
86+
method = snakemake@params$cor_method,
87+
exact = snakemake@params$cor_exact
88+
)
89+
# Adjust p-value
90+
p.value.adj <- p.adjust(
91+
test$p.value,
92+
method = "BH",
93+
n = length(unique.snps)
94+
)
95+
# Add row to dataframe
96+
add_row(
97+
cor.df,
98+
variant = snp,
99+
coefficient = test$estimate,
100+
p.value = test$p.value,
101+
p.value.adj = p.value.adj
102+
)
103+
}
104+
) %>%
105+
bind_rows()
106+
107+
log_info("Writing correlations table")
108+
write_csv(
109+
correlations,
110+
snakemake@output[["correlations"]]
111+
)

‎workflow/scripts/report/af_time_correlation_per_snp_plot.R‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
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(ggrepel)
10+
library(logger)
11+
log_threshold(INFO)
12+
13+
# Import file with plots style
14+
source(snakemake@params[["design"]])
15+
16+
log_info("Reading correlation data")
17+
correlations <- read_csv(snakemake@input$correlations)
18+
19+
log_info("Plotting coorrelation coefficients and p-values of each variants")
20+
p <- correlations %>%
21+
mutate(
22+
trans.p = -log10(p.value.adj),
23+
label = ifelse(p.value.adj < 0.05, variant, NA)
24+
) %>%
25+
ggplot() +
26+
aes(
27+
x = coefficient,
28+
y = trans.p
29+
) +
30+
geom_point() +
31+
geom_text_repel(aes(label = label), max.overlaps = 1000, direction = "x") +
32+
xlim(c(-1, 1)) +
33+
geom_hline(
34+
aes(yintercept = -log10(0.05)),
35+
linetype = 2,
36+
color = "orange"
37+
) +
38+
labs(
39+
x = "Correlation coefficient",
40+
y = "-log10(p-value)"
41+
)
42+
43+
ggsave(
44+
filename = snakemake@output$plot,
45+
plot = p,
46+
width = 159.2,
47+
height = 119.4,
48+
units = "mm",
49+
dpi = 250
50+
)

0 commit comments

Comments
 (0)