Skip to content

Commit a0b3c90

Browse files
ahmigahmig
authored
refactor: move data analysis rules to core snakefiles (#36)
* refactor: move window data rule to VAF snakefile * chore: remove unused column from selection * perf: use set for masked site lookup and validate sequence lengths * refactor: move NV panel data rules to VCF snakefile * refactor: move correlation rules to VCF snakefile * refactor: move rate analyses to core snakefiles * feat: calculate confidence interval of evolutionary rate
1 parent 7dd4216 commit a0b3c90

6 files changed

Lines changed: 157 additions & 148 deletions

File tree

‎workflow/rules/distances.smk‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,3 +50,35 @@ rule format_afwdist_results:
5050
LOGDIR/"format_afwdist_results"/"log.txt"
5151
script:
5252
"../scripts/format_afwdist_results.py"
53+
54+
55+
rule allele_freq_tree_data:
56+
conda: "../envs/renv.yaml"
57+
params:
58+
use_bionj = config["USE_BIONJ"],
59+
outgroup_id = config["ALIGNMENT_REFERENCE"],
60+
input:
61+
dist = OUTDIR/f"{OUTPUT_NAME}.distances.csv",
62+
output:
63+
tree = REPORT_DIR_TABLES/"allele_freq_tree.nwk",
64+
log:
65+
LOGDIR / "allele_freq_tree_data" / "log.txt"
66+
script:
67+
"../scripts/report/allele_freq_tree_data.R"
68+
69+
70+
rule time_signal_data:
71+
conda: "../envs/renv.yaml"
72+
params:
73+
outgroup_id = config["ALIGNMENT_REFERENCE"],
74+
confidence_interval = 0.95,
75+
input:
76+
tree = report(REPORT_DIR_TABLES/"allele_freq_tree.nwk"),
77+
metadata = config["METADATA"],
78+
output:
79+
table = report(REPORT_DIR_TABLES/"time_signal.csv"),
80+
json = REPORT_DIR_TABLES/"time_signal.json",
81+
log:
82+
LOGDIR / "time_signal_data" / "log.txt"
83+
script:
84+
"../scripts/report/time_signal_data.R"

‎workflow/rules/report.smk‎

Lines changed: 0 additions & 139 deletions
Original file line numberDiff line numberDiff line change
@@ -78,22 +78,6 @@ rule extract_genbank_regions:
7878
"../scripts/report/extract_genbank_regions.py"
7979

8080

81-
rule polymorphic_sites_over_time_data:
82-
conda: "../envs/renv.yaml"
83-
params:
84-
max_alt_freq = 1.0 - config["VC"]["MIN_FREQ"],
85-
input:
86-
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
87-
metadata = config["METADATA"],
88-
output:
89-
table = REPORT_DIR_PLOTS/"polymorphic_sites_over_time.csv",
90-
json = temp(REPORT_DIR_TABLES/"polymorphic_sites_over_time.json"),
91-
log:
92-
LOGDIR / "polymorphic_sites_over_time_data" / "log.txt"
93-
script:
94-
"../scripts/report/polymorphic_sites_over_time_data.R"
95-
96-
9781
rule polymorphic_sites_over_time_plot:
9882
conda: "../envs/renv.yaml"
9983
params:
@@ -110,63 +94,6 @@ rule polymorphic_sites_over_time_plot:
11094
"../scripts/report/polymorphic_sites_over_time_plot.R"
11195

11296

113-
rule window_data:
114-
conda: "../envs/biopython.yaml"
115-
params:
116-
window = config["WINDOW"]["WIDTH"],
117-
step = config["WINDOW"]["STEP"],
118-
features = config.get("GB_FEATURES", {}),
119-
gb_qualifier_display = "gene"
120-
input:
121-
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
122-
gb = OUTDIR/"reference.gb",
123-
output:
124-
window_df = REPORT_DIR_TABLES/"window.csv",
125-
json = temp(REPORT_DIR_TABLES/"window.json"),
126-
log:
127-
LOGDIR / "window_data" / "log.txt"
128-
script:
129-
"../scripts/report/window_data.py"
130-
131-
132-
rule nv_panel_data:
133-
conda: "../envs/renv.yaml"
134-
input:
135-
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
136-
metadata = config["METADATA"],
137-
output:
138-
table = REPORT_DIR_TABLES/"nv_panel.csv",
139-
json = temp(REPORT_DIR_TABLES/"nv_panel.json"),
140-
log:
141-
LOGDIR / "nv_panel_data" / "log.txt"
142-
script:
143-
"../scripts/report/nv_panel_data.R"
144-
145-
146-
rule nv_panel_zoom_on_feature_data:
147-
input:
148-
table = REPORT_DIR_TABLES/"nv_panel.csv",
149-
regions = REPORT_DIR_TABLES/"genbank_regions.json",
150-
output:
151-
table = temp(REPORT_DIR_TABLES/"nv_panel.{region_name}.csv"),
152-
log:
153-
LOGDIR / "nv_panel_zoom_on_feature_data" / "{region_name}.log.txt"
154-
script:
155-
"../scripts/report/nv_panel_zoom_on_feature_data.py"
156-
157-
158-
rule window_zoom_on_feature_data:
159-
input:
160-
table = REPORT_DIR_TABLES/"window.csv",
161-
regions = REPORT_DIR_TABLES/"genbank_regions.json",
162-
output:
163-
table = temp(REPORT_DIR_TABLES/"window.{region_name}.csv"),
164-
log:
165-
LOGDIR / "window_zoom_on_feature_data" / "{region_name}.log.txt"
166-
script:
167-
"../scripts/report/window_zoom_on_feature_data.py"
168-
169-
17097
rule nv_panel_plot:
17198
conda: "../envs/renv.yaml"
17299
params:
@@ -256,21 +183,6 @@ rule context_phylogeny_plot:
256183
"../scripts/report/context_phylogeny_plot.R"
257184

258185

259-
rule allele_freq_tree_data:
260-
conda: "../envs/renv.yaml"
261-
params:
262-
use_bionj = config["USE_BIONJ"],
263-
outgroup_id = config["ALIGNMENT_REFERENCE"],
264-
input:
265-
dist = OUTDIR/f"{OUTPUT_NAME}.distances.csv",
266-
output:
267-
tree = REPORT_DIR_TABLES/"allele_freq_tree.nwk",
268-
log:
269-
LOGDIR / "allele_freq_tree_data" / "log.txt"
270-
script:
271-
"../scripts/report/allele_freq_tree_data.R"
272-
273-
274186
rule allele_freq_tree_plot:
275187
conda: "../envs/renv.yaml"
276188
params:
@@ -290,22 +202,6 @@ rule allele_freq_tree_plot:
290202
"../scripts/report/allele_freq_tree_plot.R"
291203

292204

293-
rule time_signal_data:
294-
conda: "../envs/renv.yaml"
295-
params:
296-
outgroup_id = config["ALIGNMENT_REFERENCE"],
297-
input:
298-
tree = report(REPORT_DIR_TABLES/"allele_freq_tree.nwk"),
299-
metadata = config["METADATA"],
300-
output:
301-
table = report(REPORT_DIR_TABLES/"time_signal.csv"),
302-
json = REPORT_DIR_TABLES/"time_signal.json",
303-
log:
304-
LOGDIR / "time_signal_data" / "log.txt"
305-
script:
306-
"../scripts/report/time_signal_data.R"
307-
308-
309205
rule time_signal_plot:
310206
conda: "../envs/renv.yaml"
311207
params:
@@ -339,24 +235,6 @@ rule dnds_plots:
339235
"../scripts/report/dnds_plots.R"
340236

341237

342-
rule af_time_correlation_data:
343-
conda: "../envs/renv.yaml"
344-
params:
345-
cor_method = config["COR"]["METHOD"],
346-
cor_exact = config["COR"]["EXACT"],
347-
input:
348-
variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv",
349-
metadata = config["METADATA"],
350-
output:
351-
fmt_variants = temp(REPORT_DIR_TABLES/"variants.filled.dated.tsv"),
352-
correlations = report(REPORT_DIR_TABLES/"af_time_correlation.csv"),
353-
subset = REPORT_DIR_TABLES/"af_time_correlation.subset.txt",
354-
log:
355-
LOGDIR / "af_time_correlation_data" / "log.txt"
356-
script:
357-
"../scripts/report/af_time_correlation_data.R"
358-
359-
360238
rule af_time_correlation_plot:
361239
conda: "../envs/renv.yaml"
362240
params:
@@ -392,23 +270,6 @@ rule af_trajectory_panel_plot:
392270
"../scripts/report/af_trajectory_panel_plot.R"
393271

394272

395-
rule pairwise_trajectory_correlation_data:
396-
conda: "../envs/renv.yaml"
397-
params:
398-
cor_method = config["COR"]["METHOD"],
399-
cor_use = "pairwise.complete.obs",
400-
input:
401-
variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv",
402-
metadata = config["METADATA"],
403-
output:
404-
table = REPORT_DIR_TABLES/"pairwise_trajectory_frequency_data.csv",
405-
matrix = report(REPORT_DIR_TABLES/"pairwise_trajectory_correlation_matrix.csv"),
406-
log:
407-
LOGDIR / "pairwise_trajectory_correlation_data" / "log.txt"
408-
script:
409-
"../scripts/report/pairwise_trajectory_correlation_data.R"
410-
411-
412273
rule summary_table:
413274
conda: "../envs/renv.yaml"
414275
input:

‎workflow/rules/vaf.smk‎

Lines changed: 115 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -167,7 +167,8 @@ rule format_vcf_fields_longer:
167167
filter_exclude=config["ANNOTATION"]["FILTER_EXCLUDE"],
168168
variant_name_pattern=lambda wildcards: config["ANNOTATION"][
169169
"VARIANT_NAME_PATTERN"
170-
], # lambda to deactivate automatic wildcard expansion in pattern
170+
],
171+
# lambda to deactivate automatic wildcard expansion in pattern
171172
sep=",",
172173
input:
173174
tsv=OUTDIR / "vaf" / "fields" / "{sample}.tsv",
@@ -189,7 +190,6 @@ rule concat_vcf_fields:
189190
run:
190191
import pandas as pd
191192
from functools import reduce
192-
193193
reduce(
194194
lambda a, b: pd.concat((a, b), axis="rows", ignore_index=True),
195195
(pd.read_csv(path, sep=params.sep) for path in input),
@@ -219,3 +219,116 @@ use rule concat_vcf_fields as concat_variants with:
219219
expand(OUTDIR / "vaf" / "variants" / "{sample}.tsv", sample=iter_samples()),
220220
output:
221221
OUTDIR / f"{OUTPUT_NAME}.variants.tsv",
222+
223+
224+
rule window_data:
225+
conda:
226+
"../envs/biopython.yaml"
227+
params:
228+
window=config["WINDOW"]["WIDTH"],
229+
step=config["WINDOW"]["STEP"],
230+
features=config.get("GB_FEATURES", {}),
231+
gb_qualifier_display="gene",
232+
input:
233+
variants=OUTDIR / f"{OUTPUT_NAME}.variants.tsv",
234+
gb=OUTDIR / "reference.gb",
235+
output:
236+
window_df=REPORT_DIR_TABLES / "window.csv",
237+
json=temp(REPORT_DIR_TABLES / "window.json"),
238+
log:
239+
LOGDIR / "window_data" / "log.txt",
240+
script:
241+
"../scripts/report/window_data.py"
242+
243+
244+
rule nv_panel_data:
245+
conda:
246+
"../envs/renv.yaml"
247+
input:
248+
variants=OUTDIR / f"{OUTPUT_NAME}.variants.tsv",
249+
metadata=config["METADATA"],
250+
output:
251+
table=REPORT_DIR_TABLES / "nv_panel.csv",
252+
json=temp(REPORT_DIR_TABLES / "nv_panel.json"),
253+
log:
254+
LOGDIR / "nv_panel_data" / "log.txt",
255+
script:
256+
"../scripts/report/nv_panel_data.R"
257+
258+
259+
rule nv_panel_zoom_on_feature_data:
260+
input:
261+
table=REPORT_DIR_TABLES / "nv_panel.csv",
262+
regions=REPORT_DIR_TABLES / "genbank_regions.json",
263+
output:
264+
table=temp(REPORT_DIR_TABLES / "nv_panel.{region_name}.csv"),
265+
log:
266+
LOGDIR / "nv_panel_zoom_on_feature_data" / "{region_name}.log.txt",
267+
script:
268+
"../scripts/report/nv_panel_zoom_on_feature_data.py"
269+
270+
271+
rule window_zoom_on_feature_data:
272+
input:
273+
table=REPORT_DIR_TABLES / "window.csv",
274+
regions=REPORT_DIR_TABLES / "genbank_regions.json",
275+
output:
276+
table=temp(REPORT_DIR_TABLES / "window.{region_name}.csv"),
277+
log:
278+
LOGDIR / "window_zoom_on_feature_data" / "{region_name}.log.txt",
279+
script:
280+
"../scripts/report/window_zoom_on_feature_data.py"
281+
282+
283+
rule af_time_correlation_data:
284+
conda:
285+
"../envs/renv.yaml"
286+
params:
287+
cor_method=config["COR"]["METHOD"],
288+
cor_exact=config["COR"]["EXACT"],
289+
input:
290+
variants=OUTDIR / f"{OUTPUT_NAME}.variants.all_sites.tsv",
291+
metadata=config["METADATA"],
292+
output:
293+
fmt_variants=temp(REPORT_DIR_TABLES / "variants.filled.dated.tsv"),
294+
correlations=report(REPORT_DIR_TABLES / "af_time_correlation.csv"),
295+
subset=REPORT_DIR_TABLES / "af_time_correlation.subset.txt",
296+
log:
297+
LOGDIR / "af_time_correlation_data" / "log.txt",
298+
script:
299+
"../scripts/report/af_time_correlation_data.R"
300+
301+
302+
rule pairwise_trajectory_correlation_data:
303+
conda:
304+
"../envs/renv.yaml"
305+
params:
306+
cor_method=config["COR"]["METHOD"],
307+
cor_use="pairwise.complete.obs",
308+
input:
309+
variants=OUTDIR / f"{OUTPUT_NAME}.variants.all_sites.tsv",
310+
metadata=config["METADATA"],
311+
output:
312+
table=REPORT_DIR_TABLES / "pairwise_trajectory_frequency_data.csv",
313+
matrix=report(REPORT_DIR_TABLES / "pairwise_trajectory_correlation_matrix.csv"),
314+
log:
315+
LOGDIR / "pairwise_trajectory_correlation_data" / "log.txt",
316+
script:
317+
"../scripts/report/pairwise_trajectory_correlation_data.R"
318+
319+
320+
rule polymorphic_sites_over_time_data:
321+
conda:
322+
"../envs/renv.yaml"
323+
params:
324+
max_alt_freq=1.0 - config["VC"]["MIN_FREQ"],
325+
input:
326+
variants=OUTDIR / f"{OUTPUT_NAME}.variants.tsv",
327+
metadata=config["METADATA"],
328+
output:
329+
table=REPORT_DIR_PLOTS / "polymorphic_sites_over_time.csv",
330+
json=temp(REPORT_DIR_TABLES / "polymorphic_sites_over_time.json"),
331+
log:
332+
LOGDIR / "polymorphic_sites_over_time_data" / "log.txt",
333+
script:
334+
"../scripts/report/polymorphic_sites_over_time_data.R"

‎workflow/scripts/calculate_dnds.R‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,6 @@ variants <- read_delim(
1919
col_select = c(
2020
"SAMPLE",
2121
"VARIANT_NAME",
22-
"CHROM",
2322
"ALT_FREQ",
2423
"SYNONYMOUS",
2524
"POS"

0 commit comments

Comments
 (0)