Skip to content

Commit 7978a48

Browse files
ahmigahmig
committed
Decouple NV panel plot from data
1 parent 8daeea0 commit 7978a48

10 files changed

Lines changed: 455 additions & 529 deletions

File tree

‎config/README.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -180,7 +180,7 @@ All of the following variables are pre-defined in [config.yaml](/config/config.y
180180
- `EXACT`: boolean flag indicating whether to compute an exact p-value when possible. This option applies only to certain methods and may be set to `null` (default) to let R decide automatically.
181181
- `LOG_PY_FMT`: logging format string for Python scripts.
182182
- `PLOTS`: path of the R script that sets the design and style of data visualizations.
183-
- `PLOT_GENOME_REGIONS`: path of a CSV file containing genome regions, e.g. SARS-CoV-2 non-structural protein (NSP) coordinates, for data visualization.
183+
- `PLOT_GENOME_REGIONS`: path of a CSV file containing genome regions, e.g. SARS-CoV-2 non-structural protein (NSP) coordinates, for data visualization (columns: `region`, `start`, `end`).
184184
- `REPORT_QMD`: path of the report template in Quarto markdown (QMD) format.
185185

186186
TODO: describe new entries in config.yaml

‎config/nsp_annotation.csv‎

Lines changed: 6 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,7 @@
11
region,start,end
2-
nsp1,266,805
3-
nsp2,806,2719
4-
nsp3,2720,8554
5-
nsp4,8555,10054
6-
nsp5,10055,10972
7-
nsp6,10973,11842
8-
nsp7,11843,12091
9-
nsp8,12092,12685
10-
nsp9,12686,13024
11-
nsp10,13025,13441
12-
nsp11,13442,13468
13-
nsp12,13468,16236
14-
nsp13,16237,18039
15-
nsp14,18040,19620
16-
nsp15,19621,20658
17-
nsp16,20659,21552
2+
NSP2,806,2719
3+
NSP3,2720,8554
4+
NSP4-11,8555,13468
5+
RdRP,13468,16236
6+
NSP13,16237,18039
7+
NSP14-16,18040,21552

‎workflow/rules/fetch.smk‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,3 +66,8 @@ rule fetch_problematic_vcf:
6666
select_problematic_vcf()
6767
shell:
6868
"curl {params.url} -o {output} -s 2> {log}"
69+
70+
71+
rule create_empty_file:
72+
output:
73+
temp(touch(OUTDIR/"empty.txt"))

‎workflow/rules/report.smk‎

Lines changed: 115 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,8 @@ rule window:
5151
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
5252
gb = OUTDIR/"reference.gb",
5353
output:
54-
window_df = temp(OUTDIR/f"{OUTPUT_NAME}.window.csv"),
54+
window_df = temp(REPORT_DIR_TABLES/"window.csv"),
55+
json = temp(REPORT_DIR_TABLES/"window.json"),
5556
log:
5657
LOGDIR / "window" / "log.txt"
5758
script:
@@ -109,32 +110,122 @@ rule extract_genbank_regions:
109110
"../scripts/report/extract_genbank_regions.py"
110111

111112

112-
rule general_NV_description:
113+
rule polymorphic_sites_over_time_data:
114+
conda: "../envs/renv.yaml"
115+
params:
116+
max_alt_freq = 1.0 - config["VC"]["IVAR_FREQ"],
117+
input:
118+
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
119+
metadata = config["METADATA"],
120+
output:
121+
table = REPORT_DIR_PLOTS/"polymorphic_sites_over_time.csv",
122+
json = temp(REPORT_DIR_TABLES/"polymorphic_sites_over_time.json"),
123+
log:
124+
LOGDIR / "polymorphic_sites_over_time_data" / "log.txt"
125+
script:
126+
"../scripts/report/polymorphic_sites_over_time_data.R"
127+
128+
129+
rule polymorphic_sites_over_time_plot:
113130
conda: "../envs/renv.yaml"
114131
params:
115-
samples = expand("{sample}", sample = iter_samples()),
116132
design = config["PLOTS"],
117-
window = config["WINDOW"]["WIDTH"],
118-
step = config["WINDOW"]["STEP"],
119-
max_alt_freq = 1.0 - config["VC"]["IVAR_FREQ"]
120133
input:
121-
coordinates = REPORT_DIR_TABLES/"genbank_regions.json",
122-
regions = config["PLOT_GENOME_REGIONS"],
123-
window = OUTDIR/f"{OUTPUT_NAME}.window.csv",
124-
vcf = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
125-
metadata = config["METADATA"]
134+
table = REPORT_DIR_PLOTS/"polymorphic_sites_over_time.csv",
126135
output:
127-
fig = report(REPORT_DIR_PLOTS/"figure_5a.png"), # panel
128-
fig_s = report(REPORT_DIR_PLOTS/"figure_5b.png"), # panel spike
129-
fig_cor = report(REPORT_DIR_PLOTS/"figure_4.png"), # SNPs vs time
130-
json = temp(OUTDIR/"summary_nv.json"),
131-
table_1 = report(REPORT_DIR_TABLES/"figure_5a.csv"),
132-
table_2 = report(REPORT_DIR_TABLES/"figure_5b.csv"),
133-
table_3 = report(REPORT_DIR_TABLES/"figure_4.csv")
136+
plot = report(REPORT_DIR_PLOTS/"polymorphic_sites_over_time.png"),
134137
log:
135-
LOGDIR / "general_NV_description" / "log.txt"
138+
LOGDIR / "polymorphic_sites_over_time_plot" / "log.txt"
136139
script:
137-
"../scripts/report/NV_description.R"
140+
"../scripts/report/polymorphic_sites_over_time_plot.R"
141+
142+
143+
rule nv_panel_data:
144+
conda: "../envs/renv.yaml"
145+
input:
146+
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
147+
metadata = config["METADATA"],
148+
output:
149+
table = REPORT_DIR_TABLES/"nv_panel.csv",
150+
json = temp(REPORT_DIR_TABLES/"nv_panel.json"),
151+
log:
152+
LOGDIR / "nv_panel_data" / "log.txt"
153+
script:
154+
"../scripts/report/nv_panel_data.R"
155+
156+
157+
rule nv_panel_zoom_on_feature_data:
158+
input:
159+
table = REPORT_DIR_TABLES/"nv_panel.csv",
160+
regions = REPORT_DIR_TABLES/"genbank_regions.json",
161+
output:
162+
table = temp(REPORT_DIR_TABLES/"nv_panel.{region_name}.csv"),
163+
log:
164+
LOGDIR / "nv_panel_zoom_on_feature_data" / "{region_name}.log.txt"
165+
script:
166+
"../scripts/report/nv_panel_zoom_on_feature_data.py"
167+
168+
169+
rule window_zoom_on_feature_data:
170+
input:
171+
table = REPORT_DIR_TABLES/"window.csv",
172+
regions = REPORT_DIR_TABLES/"genbank_regions.json",
173+
output:
174+
table = temp(REPORT_DIR_TABLES/"window.{region_name}.csv"),
175+
log:
176+
LOGDIR / "window_zoom_on_feature_data" / "{region_name}.log.txt"
177+
script:
178+
"../scripts/report/window_zoom_on_feature_data.py"
179+
180+
181+
rule nv_panel_plot:
182+
conda: "../envs/renv.yaml"
183+
params:
184+
design = config["PLOTS"],
185+
window_step = config["WINDOW"]["STEP"],
186+
plot_height_mm = 250,
187+
plot_width_mm = 240,
188+
input:
189+
panel = REPORT_DIR_TABLES/"nv_panel.csv",
190+
window = REPORT_DIR_TABLES/"window.csv",
191+
regions = REPORT_DIR_TABLES/"genbank_regions.json",
192+
highlight_window_regions = config["PLOT_GENOME_REGIONS"],
193+
output:
194+
plot = report(REPORT_DIR_PLOTS/"nv_panel.png"),
195+
log:
196+
LOGDIR / "nv_panel_plot" / "log.txt"
197+
script:
198+
"../scripts/report/nv_panel_plot.R"
199+
200+
201+
use rule nv_panel_plot as nv_panel_plot_S with:
202+
input:
203+
panel = REPORT_DIR_TABLES/"nv_panel.S.csv",
204+
window = REPORT_DIR_TABLES/"window.S.csv",
205+
regions = REPORT_DIR_TABLES/"genbank_regions.json",
206+
highlight_window_regions = OUTDIR/"empty.txt",
207+
output:
208+
plot = report(REPORT_DIR_PLOTS/"nv_panel.S.png"),
209+
log:
210+
LOGDIR / "nv_panel_plot_S" / "log.txt"
211+
212+
213+
rule merge_json_files:
214+
input:
215+
REPORT_DIR_TABLES/"nv_panel.json",
216+
REPORT_DIR_TABLES/"polymorphic_sites_over_time.json",
217+
REPORT_DIR_TABLES/"window.json",
218+
output:
219+
json = REPORT_DIR_TABLES/"nv_panel_summary.json",
220+
run:
221+
import json
222+
result = {}
223+
for path in input:
224+
with open(path) as f:
225+
d = json.load(f)
226+
result |= d # will replace existing keys
227+
with open(output.json, "w") as fw:
228+
json.dump(result, fw, indent=2)
138229

139230

140231
rule context_phylogeny_data:
@@ -329,9 +420,9 @@ rule report:
329420
demix = report(REPORT_DIR_PLOTS/"demix.png"),
330421
tree_ml = report(REPORT_DIR_PLOTS/"context_phylogeny.png"),
331422
diversity = report(REPORT_DIR_PLOTS/"diversity.png"),
332-
fig_cor = report(REPORT_DIR_PLOTS/"figure_4.png"),
333-
SNV = report(REPORT_DIR_PLOTS/"figure_5a.png"),
334-
SNV_spike = report(REPORT_DIR_PLOTS/"figure_5b.png"),
423+
fig_cor = report(REPORT_DIR_PLOTS/"polymorphic_sites_over_time.png"),
424+
SNV = report(REPORT_DIR_PLOTS/"nv_panel.png"),
425+
SNV_spike = report(REPORT_DIR_PLOTS/"nv_panel.S.png"),
335426
volcano = report(REPORT_DIR_PLOTS/"af_time_correlation.png"),
336427
panel = report(REPORT_DIR_PLOTS/"af_trajectory_panel.png"),
337428
tree = report(REPORT_DIR_PLOTS/"allele_freq_tree.png"),
@@ -344,7 +435,7 @@ rule report:
344435
stats_lm = REPORT_DIR_TABLES/"time_signal.json",
345436
stats_ml = REPORT_DIR_TABLES/"context_phylogeny.json",
346437
table = OUTDIR/"summary_table.csv",
347-
sum_nv = OUTDIR/"summary_nv.json",
438+
sum_nv = REPORT_DIR_TABLES/"nv_panel_summary.json",
348439
params:
349440
workflow_version = get_repo_version(BASE_PATH.as_posix(), __version__),
350441
min_ivar_freq = config["VC"]["IVAR_FREQ"],

0 commit comments

Comments
 (0)