Skip to content

Commit 392bbd6

Browse files
ahmigahmig
committed
Merge branch 'main' into refactor/snakefiles
2 parents f89a3bc + 29d75e5 commit 392bbd6

11 files changed

Lines changed: 142 additions & 80 deletions

‎template.qmd‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,6 @@ stats <- append(
7777
)
7878
correlation <- stats[["r2"]]
7979
sub_rate <- stats[["sub_rate"]]
80-
sub_rate <- display.num(sub_rate, 2)
8180
p_value_lm <- stats[["pvalue"]]
8281
8382
# NV counts
@@ -262,7 +261,7 @@ frequency-weighted distances.](`r params$tree`){#fig-tree}
262261
To estimate the evolutionary rate, root-to-tip distances measured on the previous
263262
tree (@fig-tree) have been correlated with time, obtaining a $R^2$ of
264263
$`r display.num(correlation, 4)`$ and a p-value of $`r p_value_lm`$. The estimated evolutionary
265-
rate is $`r sub_rate`$ number of changes per year (@fig-tempest).
264+
rate is $`r display.num(sub_rate, 2)`$ changes per year (@fig-tempest).
266265

267266
![Scatterplot depicting the relationship between root-to-tip
268267
distances and the number of days passed since the first sample. The solid

‎workflow/core.smk‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,5 +17,6 @@ include: "rules/fetch.smk"
1717
include: "rules/fasta.smk"
1818
include: "rules/asr.smk"
1919
include: "rules/vaf.smk"
20+
include: "rules/sites.smk"
2021
include: "rules/distances.smk"
2122
include: "rules/evolution.smk"

‎workflow/rules/common.smk‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ def get_repo_version(base_dir: str, default: str, warn=False) -> str:
2929
"""
3030
try:
3131
last_tag_description = subprocess.check_output(
32-
f"git --git-dir={base_dir}/.git describe --always",
32+
f"git --git-dir={base_dir}/.git describe --tags --always",
3333
shell=True,
3434
stderr=subprocess.DEVNULL
3535
).strip().decode("utf-8")

‎workflow/rules/context.smk‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@ rule download_context:
2828
duplicate_accids = OUTDIR/"context"/"duplicate_accession_ids.txt"
2929
log:
3030
LOGDIR / "download_context" / "log.txt"
31+
retries: 2
3132
script:
3233
"../scripts/download_context.R"
3334

‎workflow/rules/distances.smk‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ rule extract_afwdist_variants:
88
mask_class = ["mask"],
99
input:
1010
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
11-
mask_vcf = lambda wildcards: select_problematic_vcf(),
11+
mask_vcf = OUTDIR / "all_mask_sites.vcf",
1212
ancestor = OUTDIR/f"{OUTPUT_NAME}.ancestor.fasta",
1313
reference = OUTDIR/"reference.fasta",
1414
output:

‎workflow/rules/sites.smk‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,123 @@
1+
rule bcftools_mpileup_all_sites:
2+
threads: 1
3+
conda: "../envs/var_calling.yaml"
4+
params:
5+
min_mq = 0,
6+
min_bq = config["VC"]["MIN_QUALITY"],
7+
mpileup_extra = "--no-BAQ"
8+
input:
9+
bam = get_input_bam,
10+
reference = OUTDIR/"vaf"/"{sample}.reference.fasta",
11+
output:
12+
mpileup = temp(OUTDIR / "all_sites" / "{sample}.mpileup.vcf"),
13+
query = temp(OUTDIR / "all_sites" / "{sample}.query.tsv"),
14+
log:
15+
mpileup = LOGDIR / "bcftools_mpileup_all_sites" / "{sample}.mpileup.txt",
16+
query = LOGDIR / "bcftools_mpileup_all_sites" / "{sample}.query.txt",
17+
shell:
18+
"bcftools mpileup {params.mpileup_extra} -a AD,ADF,ADR --fasta-ref {input.reference:q} --threads {threads} -Q {params.min_bq} -q {params.min_mq} -Ov -o {output.mpileup:q} {input.bam:q} >{log.mpileup:q} 2>&1 && "
19+
"echo 'CHROM\tPOS\tREF\tALT\tDP\tAD\tADF\tADR' >{output.query:q} && "
20+
"bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%DP\t[ %AD]\t[ %ADF]\t[ %ADR]\n' {output.mpileup:q} >>{output.query:q} 2>{log.query:q}"
21+
22+
23+
rule filter_mpileup_all_sites:
24+
threads: 1
25+
params:
26+
min_total_AD = config["VC"]["MIN_DEPTH"],
27+
min_total_ADF = 0,
28+
min_total_ADR = 0,
29+
input:
30+
OUTDIR / "all_sites" / "{sample}.query.tsv",
31+
output:
32+
sites_pass = temp(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv"),
33+
sites_fail = temp(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv"),
34+
log:
35+
LOGDIR / "filter_mpileup_all_sites" / "{sample}.txt"
36+
run:
37+
import pandas as pd
38+
df = pd.read_csv(input[0], sep="\t")
39+
df["SAMPLE"] = wildcards.sample
40+
df["REF_AD"] = df.AD.str.split(",").apply(lambda values: int(values[0]))
41+
df["TOTAL_AD"] = df.AD.str.split(",").apply(lambda values: sum(int(n) for n in values))
42+
df["TOTAL_ADF"] = df.ADF.str.split(",").apply(lambda values: sum(int(n) for n in values))
43+
df["TOTAL_ADR"] = df.ADR.str.split(",").apply(lambda values: sum(int(n) for n in values))
44+
mask = (
45+
(df.TOTAL_AD >= params.min_total_AD) &
46+
(df.TOTAL_ADF >= params.min_total_ADF) &
47+
(df.TOTAL_ADR >= params.min_total_ADR)
48+
)
49+
df[mask].to_csv(output.sites_pass, sep="\t", index=False)
50+
df[~mask].to_csv(output.sites_fail, sep="\t", index=False)
51+
52+
53+
use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with:
54+
input:
55+
expand(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv", sample=iter_samples()),
56+
output:
57+
OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
58+
59+
60+
use rule concat_vcf_fields as merge_fail_mpileup_all_sites with:
61+
input:
62+
expand(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv", sample=iter_samples()),
63+
output:
64+
OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv",
65+
66+
67+
rule fill_all_sites:
68+
conda: "../envs/renv.yaml"
69+
input:
70+
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
71+
sites = OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
72+
output:
73+
variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv",
74+
log:
75+
LOGDIR / "fill_all_sites" / "log.txt"
76+
script:
77+
"../scripts/fill_all_sites.R"
78+
79+
80+
rule compile_fail_sites_vcf:
81+
params:
82+
filter_text = "mask",
83+
sub_text = "NA",
84+
exc_text = "site_qual",
85+
input:
86+
sites = OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv",
87+
output:
88+
sites = temp(OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf"),
89+
run:
90+
import pandas as pd
91+
HEADER = ["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"]
92+
sites = (
93+
pd.read_table(input.sites, sep="\t")
94+
.drop_duplicates(subset=("CHROM", "POS", "REF"))
95+
.rename(columns={"CHROM": "#CHROM"})
96+
)
97+
sites["ID"] = "."
98+
sites["ALT"] = "."
99+
sites["QUAL"] = "."
100+
sites["FILTER"] = params.filter_text
101+
sites["INFO"] = f"SUB={params.sub_text};EXC={params.exc_text}"
102+
sites[HEADER].to_csv(output.sites, sep="\t", index=False)
103+
104+
105+
rule merge_mask_sites_vcf:
106+
input:
107+
lambda wildcards: select_problematic_vcf(),
108+
OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf",
109+
output:
110+
sites = temp(OUTDIR / "all_mask_sites.vcf"),
111+
run:
112+
import pandas as pd
113+
HEADER = ["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"]
114+
(
115+
pd.concat(
116+
[pd.read_table(path, sep="\t", comment="#", names=HEADER, dtype={"POS": "int64"}) for path in input],
117+
axis="rows",
118+
ignore_index=True
119+
)
120+
.drop_duplicates(subset=("#CHROM", "POS", "FILTER"), keep="first")
121+
.sort_values(by=["#CHROM", "POS"])
122+
.to_csv(output.sites, sep="\t", index=False)
123+
)

‎workflow/rules/vaf.smk‎

Lines changed: 0 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -206,75 +206,6 @@ use rule concat_vcf_fields as concat_variants with:
206206
OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
207207

208208

209-
rule bcftools_mpileup_all_sites:
210-
threads: 1
211-
conda: "../envs/var_calling.yaml"
212-
params:
213-
min_mq = 0,
214-
min_bq = config["VC"]["MIN_QUALITY"],
215-
mpileup_extra = "--no-BAQ"
216-
input:
217-
bam = get_input_bam,
218-
reference = OUTDIR/"vaf"/"{sample}.reference.fasta",
219-
output:
220-
mpileup = temp(OUTDIR / "all_sites" / "{sample}.mpileup.vcf"),
221-
query = temp(OUTDIR / "all_sites" / "{sample}.query.tsv"),
222-
log:
223-
mpileup = LOGDIR / "bcftools_mpileup_all_sites" / "{sample}.mpileup.txt",
224-
query = LOGDIR / "bcftools_mpileup_all_sites" / "{sample}.query.txt",
225-
shell:
226-
"bcftools mpileup {params.mpileup_extra} -a AD,ADF,ADR --fasta-ref {input.reference:q} --threads {threads} -Q {params.min_bq} -q {params.min_mq} -Ov -o {output.mpileup:q} {input.bam:q} >{log.mpileup:q} 2>&1 && "
227-
"echo 'CHROM\tPOS\tREF\tALT\tDP\tAD\tADF\tADR' >{output.query:q} && "
228-
"bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%DP\t[ %AD]\t[ %ADF]\t[ %ADR]\n' {output.mpileup:q} >>{output.query:q} 2>{log.query:q}"
229-
230-
231-
rule filter_mpileup_all_sites:
232-
threads: 1
233-
params:
234-
min_total_AD = config["VC"]["MIN_DEPTH"],
235-
min_total_ADF = 0,
236-
min_total_ADR = 0,
237-
input:
238-
OUTDIR / "all_sites" / "{sample}.query.tsv",
239-
output:
240-
temp(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv"),
241-
log:
242-
LOGDIR / "filter_mpileup_all_sites" / "{sample}.txt"
243-
run:
244-
import pandas as pd
245-
df = pd.read_csv(input[0], sep="\t")
246-
df["SAMPLE"] = wildcards.sample
247-
df["REF_AD"] = df.AD.str.split(",").apply(lambda values: int(values[0]))
248-
df["TOTAL_AD"] = df.AD.str.split(",").apply(lambda values: sum(int(n) for n in values))
249-
df["TOTAL_ADF"] = df.ADF.str.split(",").apply(lambda values: sum(int(n) for n in values))
250-
df["TOTAL_ADR"] = df.ADR.str.split(",").apply(lambda values: sum(int(n) for n in values))
251-
df[
252-
(df.TOTAL_AD >= params.min_total_AD) &
253-
(df.TOTAL_ADF >= params.min_total_ADF) &
254-
(df.TOTAL_ADR >= params.min_total_ADR)
255-
].to_csv(output[0], sep="\t", index=False)
256-
257-
258-
use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with:
259-
input:
260-
expand(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv", sample=iter_samples()),
261-
output:
262-
OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
263-
264-
265-
rule fill_all_sites:
266-
conda: "../envs/renv.yaml"
267-
input:
268-
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
269-
sites = OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
270-
output:
271-
variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv",
272-
log:
273-
LOGDIR / "fill_all_sites" / "log.txt"
274-
script:
275-
"../scripts/fill_all_sites.R"
276-
277-
278209
rule window_data:
279210
conda: "../envs/biopython.yaml"
280211
params:

‎workflow/scripts/calculate_dnds.R‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -28,9 +28,11 @@ variants <- read_delim(
2828
log_info("Reading metadata table")
2929
metadata <- read_delim(snakemake@input[["metadata"]]) %>%
3030
mutate(
31-
interval = as.numeric(
32-
as.Date(CollectionDate) - min(as.Date(CollectionDate))
33-
)
31+
interval = difftime(
32+
as.Date(CollectionDate),
33+
min(as.Date(CollectionDate)),
34+
units = "days"
35+
) |> as.numeric()
3436
) %>%
3537
select(ID, interval) %>%
3638
rename(SAMPLE = ID)

‎workflow/scripts/format_vcf_fields_longer.R‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,9 @@ read_tsv(
4545
# Rename "...[*]..." columns using the provided lookup via Snakemake config
4646
rename(all_of(unlist(snakemake@params$colnames_mapping))) %>%
4747

48+
# Ensure missing values are properly encoded
49+
mutate(across(where(is.character), ~ na_if(.x, "NA"))) %>%
50+
4851
# Separate &-delimited error column (more than one error/warning/info message per row is possible)
4952
mutate(split_errors = strsplit(ERRORS, "&")) %>%
5053
# Keep rows with none of the excluded ERRORS terms, if any

‎workflow/scripts/report/pairwise_trajectory_correlation_data.R‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ date_order <- read_csv(snakemake@input[["metadata"]]) %>%
2525

2626
log_info("Formatting variants")
2727
all_variants_wider <- variants %>%
28-
select(SAMPLE, VARIANT_NAME, ALT_FREQ) %>%
28+
distinct(SAMPLE, VARIANT_NAME, ALT_FREQ) %>%
2929
pivot_wider(
3030
names_from = VARIANT_NAME,
3131
values_from = ALT_FREQ

0 commit comments

Comments
 (0)