Skip to content

Commit 28d55d2

Browse files
ahmigahmig
committed
refactor: organize site quality filter in a new snakefile
Reverts marking failed sites in filter with a different string (back to just "mask")
1 parent c14600f commit 28d55d2

4 files changed

Lines changed: 125 additions & 126 deletions

File tree

‎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/distances.smk‎

Lines changed: 1 addition & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -1,57 +1,11 @@
1-
rule compile_fail_sites_vcf:
2-
params:
3-
filter_text = "fail_site",
4-
sub_text = "NA",
5-
exc_text = "site_qual",
6-
input:
7-
sites = OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv",
8-
output:
9-
sites = temp(OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf"),
10-
run:
11-
import pandas as pd
12-
HEADER = ["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"]
13-
sites = (
14-
pd.read_table(input.sites, sep="\t")
15-
.drop_duplicates(subset=("CHROM", "POS", "REF"))
16-
.rename(columns={"CHROM": "#CHROM"})
17-
)
18-
sites["ID"] = "."
19-
sites["ALT"] = "."
20-
sites["QUAL"] = "."
21-
sites["FILTER"] = params.filter_text
22-
sites["INFO"] = f"SUB={params.sub_text};EXC={params.exc_text}"
23-
sites[HEADER].to_csv(output.sites, sep="\t", index=False)
24-
25-
26-
rule merge_mask_sites_vcf:
27-
input:
28-
lambda wildcards: select_problematic_vcf(),
29-
OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf",
30-
output:
31-
sites = temp(OUTDIR / "all_mask_sites.vcf"),
32-
run:
33-
import pandas as pd
34-
HEADER = ["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"]
35-
(
36-
pd.concat(
37-
[pd.read_table(path, sep="\t", comment="#", names=HEADER, dtype={"POS": "int64"}) for path in input],
38-
axis="rows",
39-
ignore_index=True
40-
)
41-
.drop_duplicates(subset=("#CHROM", "POS", "FILTER"), keep="first")
42-
.sort_values(by=["#CHROM", "POS"])
43-
.to_csv(output.sites, sep="\t", index=False)
44-
)
45-
46-
471
rule extract_afwdist_variants:
482
conda: "../envs/biopython.yaml"
493
params:
504
sample_col = "SAMPLE",
515
position_col = "POS",
526
sequence_col = "ALT",
537
frequency_col = "ALT_FREQ",
54-
mask_class = ["mask", "fail_site"],
8+
mask_class = ["mask"],
559
input:
5610
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
5711
mask_vcf = OUTDIR / "all_mask_sites.vcf",

‎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 & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -204,82 +204,3 @@ use rule concat_vcf_fields as concat_variants with:
204204
expand(OUTDIR/"vaf"/"{sample}.variants.tsv", sample=iter_samples()),
205205
output:
206206
OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
207-
208-
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-
sites_pass = temp(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv"),
241-
sites_fail = temp(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv"),
242-
log:
243-
LOGDIR / "filter_mpileup_all_sites" / "{sample}.txt"
244-
run:
245-
import pandas as pd
246-
df = pd.read_csv(input[0], sep="\t")
247-
df["SAMPLE"] = wildcards.sample
248-
df["REF_AD"] = df.AD.str.split(",").apply(lambda values: int(values[0]))
249-
df["TOTAL_AD"] = df.AD.str.split(",").apply(lambda values: sum(int(n) for n in values))
250-
df["TOTAL_ADF"] = df.ADF.str.split(",").apply(lambda values: sum(int(n) for n in values))
251-
df["TOTAL_ADR"] = df.ADR.str.split(",").apply(lambda values: sum(int(n) for n in values))
252-
mask = (
253-
(df.TOTAL_AD >= params.min_total_AD) &
254-
(df.TOTAL_ADF >= params.min_total_ADF) &
255-
(df.TOTAL_ADR >= params.min_total_ADR)
256-
)
257-
df[mask].to_csv(output.sites_pass, sep="\t", index=False)
258-
df[~mask].to_csv(output.sites_fail, sep="\t", index=False)
259-
260-
261-
use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with:
262-
input:
263-
expand(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv", sample=iter_samples()),
264-
output:
265-
OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
266-
267-
268-
use rule concat_vcf_fields as merge_fail_mpileup_all_sites with:
269-
input:
270-
expand(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv", sample=iter_samples()),
271-
output:
272-
OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv",
273-
274-
275-
rule fill_all_sites:
276-
conda: "../envs/renv.yaml"
277-
input:
278-
variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
279-
sites = OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
280-
output:
281-
variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv",
282-
log:
283-
LOGDIR / "fill_all_sites" / "log.txt"
284-
script:
285-
"../scripts/fill_all_sites.R"

0 commit comments

Comments
 (0)