feat: use all-sites quality filter for the evolutionary rate calculation

Author: ahmig

workflow/rules/distances.smk

@@ -1,3 +1,50 @@
+rule compile_fail_sites_vcf:
+    params:
+        header = ("#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"),
+        filter_text = "mask",
+        sub_text = "NA",
+        exc_text = "site_qual"
+    input:
+        sites = OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv",
+    output:
+        sites = temp(OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf"),
+    run:
+        import pandas as pd
+        sites = (
+            pd.read_table(input.sites, sep="\t")
+                .drop_duplicates(subset=("CHROM", "POS", "REF"))
+                .rename(columns={"CHROM": "#CHROM"})
+        )
+        sites["ID"] = "."
+        sites["ALT"] = "."
+        sites["QUAL"] = "."
+        sites["FILTER"] = params.filter_text
+        sites["INFO"] = f"SUB={params.sub_text};EXC={params.exc_text}"
+        sites[list(params.header)].to_csv(output.sites, sep="\t", index=False)
+
+
+rule merge_sites:
+    params:
+        header = ("#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO")
+    input:
+        lambda wildcards: select_problematic_vcf(),
+        OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf",
+    output:
+        sites = temp(OUTDIR / "all_mask_sites.vcf"),
+    run:
+        import pandas as pd
+        (
+            pd.concat(
+                [pd.read_table(path, sep="\t", comment="#", names=params.header) for path in input],
+                axis="rows",
+                ignore_index=True
+            )
+                .drop_duplicates(subset=("#CHROM", "POS", "FILTER"), keep="first")
+                .sort_values(list(params.header))
+                .to_csv(output.sites, sep="\t", index=False)
+        )
+
+
 rule extract_afwdist_variants:
     conda: "../envs/biopython.yaml"
     params:
@@ -8,7 +55,7 @@ rule extract_afwdist_variants:
         mask_class = ["mask"],
     input:
         variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
-        mask_vcf = lambda wildcards: select_problematic_vcf(),
+        mask_vcf = OUTDIR / "all_mask_sites.vcf",
         ancestor = OUTDIR/f"{OUTPUT_NAME}.ancestor.fasta",
         reference = OUTDIR/"reference.fasta",
     output:

workflow/rules/vaf.smk

@@ -237,7 +237,8 @@ rule filter_mpileup_all_sites:
     input:
         OUTDIR / "all_sites" / "{sample}.query.tsv",
     output:
-        temp(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv"),
+        sites_pass = temp(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv"),
+        sites_fail = temp(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv"),
     log:
         LOGDIR / "filter_mpileup_all_sites" / "{sample}.txt"
     run:
@@ -248,11 +249,13 @@ rule filter_mpileup_all_sites:
         df["TOTAL_AD"] = df.AD.str.split(",").apply(lambda values: sum(int(n) for n in values))
         df["TOTAL_ADF"] = df.ADF.str.split(",").apply(lambda values: sum(int(n) for n in values))
         df["TOTAL_ADR"] = df.ADR.str.split(",").apply(lambda values: sum(int(n) for n in values))
-        df[
+        mask = (
             (df.TOTAL_AD >= params.min_total_AD) &
             (df.TOTAL_ADF >= params.min_total_ADF) &
             (df.TOTAL_ADR >= params.min_total_ADR)
-        ].to_csv(output[0], sep="\t", index=False)
+        )
+        df[mask].to_csv(output.sites_pass, sep="\t", index=False)
+        df[~mask].to_csv(output.sites_fail, sep="\t", index=False)
 
 
 use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with:
@@ -262,6 +265,13 @@ use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with:
         OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv",
 
 
+use rule concat_vcf_fields as merge_fail_mpileup_all_sites with:
+    input:
+        expand(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv", sample=iter_samples()),
+    output:
+        OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv",
+
+
 rule fill_all_sites:
     conda: "../envs/renv.yaml"
     input: