Skip to content

Commit afbe502

Browse files
ahmigahmig
committed
Decouple diversity plot from data
1 parent a5612e7 commit afbe502

3 files changed

Lines changed: 146 additions & 120 deletions

File tree

‎workflow/rules/report.smk‎

Lines changed: 29 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -14,12 +14,11 @@ rule demix_plot_data:
1414
rule demix_plot:
1515
conda: "../envs/renv.yaml"
1616
params:
17-
design = config["PLOTS"]
18-
input:
19-
data = REPORT_DIR_TABLES/"demix.csv"
20-
params:
17+
design = config["PLOTS"],
2118
plot_width_mm = 159.2,
2219
plot_height_mm = 119.4,
20+
input:
21+
data = REPORT_DIR_TABLES/"demix.csv"
2322
output:
2423
plot = report(REPORT_DIR_PLOTS/"demix.png")
2524
log:
@@ -58,23 +57,38 @@ rule window:
5857
"../scripts/window.py"
5958

6059

61-
rule diversity:
60+
rule diversity_data:
6261
threads: 4
6362
conda: "../envs/renv.yaml"
6463
params:
65-
design = config["PLOTS"],
6664
bootstrap_reps = config["DIVERSITY_REPS"],
67-
plot_width = 159.2,
68-
plot_height = 119.4
65+
aln_reference = config["ALIGNMENT_REFERENCE"],
6966
input:
7067
study_fasta = OUTDIR/"nextalign"/f"{OUTPUT_NAME}.aligned.masked.fasta",
71-
context_fasta = OUTDIR/"context"/"nextalign"/"context_sequences.aligned.masked.fasta"
68+
context_fasta = OUTDIR/"context"/"nextalign"/"context_sequences.aligned.masked.fasta",
69+
output:
70+
divs = report(REPORT_DIR_TABLES/"diversity.txt"),
71+
json = temp(REPORT_DIR_TABLES/"diversity.json"),
72+
log:
73+
LOGDIR / "diversity_data" / "log.txt"
74+
script:
75+
"../scripts/report/diversity_data.R"
76+
77+
78+
rule diversity_plot:
79+
threads: 1
80+
conda: "../envs/renv.yaml"
81+
params:
82+
design = config["PLOTS"],
83+
plot_width_mm = 159.2,
84+
plot_height_mm = 119.4,
85+
input:
86+
divs = report(REPORT_DIR_TABLES/"diversity.txt"),
87+
json = REPORT_DIR_TABLES/"diversity.json",
7288
output:
73-
fig = report(REPORT_DIR_PLOTS/"figure_3.png"),
74-
json = temp(OUTDIR/"diversity.json"),
75-
table = REPORT_DIR_TABLES/"figure_3.csv"
89+
plot = report(REPORT_DIR_PLOTS/"diversity.png"),
7690
log:
77-
LOGDIR / "diversity" / "log.txt"
91+
LOGDIR / "diversity_plot" / "log.txt"
7892
script:
7993
"../scripts/report/diversity_plot.R"
8094

@@ -191,7 +205,7 @@ rule report:
191205
qmd = Path(config["REPORT_QMD"]).resolve(),
192206
demix = report(REPORT_DIR_PLOTS/"demix.png"),
193207
tree_ml = report(REPORT_DIR_PLOTS/"figure_2.png"),
194-
diversity = report(REPORT_DIR_PLOTS/"figure_3.png"),
208+
diversity = report(REPORT_DIR_PLOTS/"diversity.png"),
195209
fig_cor = report(REPORT_DIR_PLOTS/"figure_4.png"),
196210
SNV = report(REPORT_DIR_PLOTS/"figure_5a.png"),
197211
SNV_spike = report(REPORT_DIR_PLOTS/"figure_5b.png"),
@@ -203,7 +217,7 @@ rule report:
203217
evo = report(REPORT_DIR_PLOTS/"figure_11.png"),
204218
omega_plot = report(REPORT_DIR_PLOTS/"figure_12.png"),
205219
freyja_ts = OUTDIR/"demixing"/"freyja_data"/"last_barcode_update.txt",
206-
value = OUTDIR/"diversity.json",
220+
value = REPORT_DIR_TABLES/"diversity.json",
207221
stats_lm = OUTDIR/"stats.lm.json",
208222
table = OUTDIR/"summary_table.csv",
209223
sum_nv = OUTDIR/"summary_nv.json",

‎workflow/scripts/report/diversity_data.R‎

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
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(ape)
9+
library(pegas)
10+
library(future.apply)
11+
library(tidyverse)
12+
library(jsonlite)
13+
library(logger)
14+
log_threshold(INFO)
15+
16+
# Pi calculation
17+
nucleotide.diversity <- function(dna_object, record.names, sample.size) {
18+
sample <- sample(record.names, sample.size, replace = FALSE)
19+
dna_subset <- dna_object[record.names %in% sample]
20+
nuc.div(dna_subset)
21+
}
22+
23+
# Parallel bootstrapping
24+
boot.nd.parallel <- function(aln, sample.size, reps) {
25+
record.names <- names(aln)
26+
future_sapply(
27+
1:reps,
28+
function(x) {
29+
nucleotide.diversity(aln, record.names, sample.size)
30+
},
31+
future.seed = TRUE
32+
)
33+
}
34+
35+
# Read outgroup/context alignment
36+
log_info("Reading context")
37+
gene_ex <- read.dna(
38+
snakemake@input[["context_fasta"]],
39+
format = "fasta",
40+
as.matrix = FALSE
41+
)
42+
gene_ex <- gene_ex[
43+
!startsWith(names(gene_ex), snakemake@params$aln_reference)
44+
]
45+
46+
# Read target (study) alignment
47+
log_info("Reading target alignment")
48+
study_aln <- read.dna(
49+
snakemake@input[["study_fasta"]],
50+
format = "fasta",
51+
as.matrix = FALSE
52+
)
53+
study_aln <- study_aln[
54+
!startsWith(names(study_aln), snakemake@params$aln_reference)
55+
]
56+
57+
# Diversity value for our samples
58+
log_info("Calculating diversity value for studied samples")
59+
diversity <- nuc.div(study_aln)
60+
61+
# Perform bootstrap
62+
log_info("Performing calculation for nucleotide diversity in context samples")
63+
plan(multisession, workers = snakemake@threads)
64+
divs <- boot.nd.parallel(
65+
gene_ex,
66+
length(study_aln),
67+
snakemake@params[["bootstrap_reps"]]
68+
)
69+
plan(sequential)
70+
71+
# Test normality
72+
log_info("Normality test for nucleotide diversity values")
73+
st <- shapiro.test(divs)
74+
75+
# Calculate p-value (assuming normal distribution)
76+
log_info("Calculating p-value (assuming normal distribution)")
77+
test <- t.test(
78+
divs,
79+
alternative = "greater",
80+
mu = diversity,
81+
conf.level = 0.95
82+
)
83+
pvalue.norm <- test$p.value
84+
85+
# Estimate p-value empirically
86+
log_info("Estimating p-value empirically")
87+
empirical.probs <- ecdf(divs)
88+
pvalue.emp <- empirical.probs(diversity)
89+
90+
# Data for JSON file
91+
log_info("Building JSON data")
92+
p.value <- ifelse(st$p.value >= 0.05, pvalue.norm, pvalue.emp)
93+
list.div <- list(
94+
"diversity" = format(diversity, scientific = TRUE),
95+
"p.value" = ifelse(p.value >= 0.001, p.value, "< 0.001"),
96+
"normal.pvalue" = ifelse(st$p.value >= 0.001, p.value, "< 0.001"),
97+
"norm.text" = ifelse(st$p.value >= 0.05, "", "not"),
98+
"type.test" = ifelse(st$p.value >= 0.05, "", "empirical"),
99+
"boot.reps" = snakemake@params[["bootstrap_reps"]],
100+
"sample.size" = length(study_aln)
101+
)
102+
103+
log_info("Writing diversity distribution")
104+
write_lines(divs, snakemake@output[["divs"]])
105+
106+
log_info("Writing results JSON")
107+
write_json(list.div, snakemake@output[["json"]])

‎workflow/scripts/report/diversity_plot.R‎

Lines changed: 10 additions & 105 deletions
Original file line numberDiff line numberDiff line change
@@ -5,94 +5,22 @@ log <- file(snakemake@log[[1]], open = "wt")
55
sink(log, type = "message")
66
sink(log, type = "output")
77

8-
library(ape)
9-
library(pegas)
10-
library(future.apply)
118
library(tidyverse)
129
library(jsonlite)
1310
library(logger)
1411
log_threshold(INFO)
1512

16-
# Pi calculation
17-
nucleotide.diversity <- function(dna_object, record.names, sample.size) {
18-
sample <- sample(record.names, sample.size, replace = FALSE)
19-
dna_subset <- dna_object[record.names %in% sample]
20-
nuc.div(dna_subset)
21-
}
22-
23-
# Parallel bootstrapping
24-
boot.nd.parallel <- function(aln, sample.size = 12, reps = 100) {
25-
record.names <- names(aln)
26-
future_sapply(
27-
1:reps,
28-
function(x) {
29-
nucleotide.diversity(aln, record.names, sample.size)
30-
},
31-
future.seed = TRUE
32-
)
33-
}
34-
3513
# Import file with plots style
3614
source(snakemake@params[["design"]])
3715

38-
# Outgroup/context alignment
39-
gene_ex <- read.dna(
40-
snakemake@input[["context_fasta"]],
41-
format = "fasta",
42-
as.matrix = FALSE
43-
)
44-
gene_ex <- gene_ex[
45-
!startsWith(names(gene_ex), snakemake@config[["ALIGNMENT_REFERENCE"]])
46-
]
47-
48-
# Study alignment
49-
study_aln <- read.dna(
50-
snakemake@input[["study_fasta"]],
51-
format = "fasta",
52-
as.matrix = FALSE
53-
)
54-
study_aln <- study_aln[
55-
!startsWith(names(study_aln), snakemake@config[["ALIGNMENT_REFERENCE"]])
56-
]
57-
58-
# Diversity value for our samples
59-
log_info("Calculating diversity value for studied samples")
60-
diversity <- nuc.div(study_aln)
61-
62-
63-
# Perform bootstrap
64-
log_info("Performing calculation for nucleotide diversity in context samples")
65-
plan(multisession, workers = snakemake@threads)
66-
divs <- boot.nd.parallel(
67-
gene_ex,
68-
length(study_aln),
69-
snakemake@params[["bootstrap_reps"]]
70-
)
71-
plan(sequential)
72-
73-
# Test normality
74-
log_info("Normality test for nucleotide diversity values")
75-
st <- shapiro.test(divs)
76-
77-
# Calculate p-value (assuming normal distribution)
78-
79-
log_info("Calculating p-value (assuming normal distribution)")
80-
81-
test <- t.test(
82-
divs,
83-
alternative = "greater",
84-
mu = diversity,
85-
conf.level = 0.95
86-
)
87-
pvalue.norm <- test$p.value
88-
89-
# Estimate p-value empirically
90-
log_info("Estimating p-value empirically")
91-
empirical.probs <- ecdf(divs)
92-
pvalue.emp <- empirical.probs(diversity)
16+
# Read data
17+
log_info("Reading diversity results")
18+
divs <- read_lines(snakemake@input$divs) %>%
19+
as.numeric()
20+
json <- read_json(snakemake@input$json)
9321

9422
# Plot and save
95-
log_info("Plotting diversity plot")
23+
log_info("Plotting")
9624
p <- data.frame(pi = divs) %>%
9725
ggplot() +
9826
geom_density(
@@ -103,7 +31,7 @@ p <- data.frame(pi = divs) %>%
10331
color = "#eae2b7"
10432
) +
10533
geom_vline(
106-
aes(xintercept = diversity),
34+
xintercept = json$diversity,
10735
color = "#d62828"
10836
) +
10937
stat_function(
@@ -117,33 +45,10 @@ p <- data.frame(pi = divs) %>%
11745
)
11846

11947
ggsave(
120-
filename = snakemake@output[["fig"]],
48+
filename = snakemake@output[["plot"]],
12149
plot = p,
122-
width = snakemake@params[["plot_width"]],
123-
height = snakemake@params[["plot_height"]],
50+
width = snakemake@params[["plot_width_mm"]],
51+
height = snakemake@params[["plot_height_mm"]],
12452
units = "mm",
12553
dpi = 250
12654
)
127-
128-
# DATA JSON #####
129-
p.value <- ifelse(st$p.value >= 0.05, pvalue.norm, pvalue.emp)
130-
131-
list.div <- list(
132-
"diversity" = format(diversity, scientific = TRUE),
133-
"p.value" = ifelse(p.value >= 0.001, p.value, "< 0.001"),
134-
"normal.pvalue" = ifelse(st$p.value >= 0.001, p.value, "< 0.001"),
135-
"norm.text" = ifelse(st$p.value >= 0.05, "", "not"),
136-
"type.test" = ifelse(st$p.value >= 0.05, "", "empirical"),
137-
"boot.reps" = snakemake@params[["bootstrap_reps"]],
138-
"sample.size" = length(study_aln)
139-
)
140-
141-
json <- toJSON(list.div)
142-
write(json, snakemake@output[["json"]])
143-
144-
# PLOT TABLES
145-
data.frame(
146-
pi = divs,
147-
prop.value = diversity
148-
) %>%
149-
write.csv(snakemake@output[["table"]], row.names = FALSE)

0 commit comments

Comments
 (0)