Update vcf filtering and clair3 command

Author: BioWilko

artic/minion.py

@@ -158,7 +158,7 @@ def run(parser, args):
             cmds.append(f"samtools index {args.sample}.{p}.trimmed.rg.sorted.bam")
 
             cmds.append(
-                f"run_clair3.sh --enable_long_indel --chunk_size=10000 --no_phasing_for_fa --bam_fn='{args.sample}.{p}.trimmed.rg.sorted.bam' --ref_fn='{ref}' --output='{args.sample}_rg_{p}' --threads='{args.threads}' --platform='ont' --model_path='{full_model_path}' --include_all_ctgs --min_coverage='{args.min_depth}'"
+                f"run_clair3.sh --enable_long_indel --chunk_size=10000 --haploid_precise --no_phasing_for_fa --bam_fn='{args.sample}.{p}.trimmed.rg.sorted.bam' --ref_fn='{ref}' --output='{args.sample}_rg_{p}' --threads='{args.threads}' --platform='ont' --model_path='{full_model_path}' --include_all_ctgs --min_coverage='{args.min_depth}'"
             )
 
             cmds.append(
@@ -183,7 +183,7 @@ def run(parser, args):
     fs_str = "--no-frameshifts" if args.no_frameshifts else ""
     indel_str = "--no-indels" if args.no_indels else ""
     cmds.append(
-        f"artic_vcf_filter {fs_str} {indel_str} --min-variant-quality {args.min_variant_quality} --min-depth {args.min_depth} {pre_filter_vcf}.gz {args.sample}.pass.vcf {args.sample}.fail.vcf"
+        f"artic_vcf_filter {fs_str} {indel_str} --min-depth {args.min_depth} {pre_filter_vcf}.gz {args.sample}.pass.vcf {args.sample}.fail.vcf"
     )
 
     # 9) get the depth of coverage for each readgroup, create a coverage mask and plots, and add failed variants to the coverage mask (artic_mask must be run before bcftools consensus)

artic/pipeline.py

@@ -104,12 +104,6 @@ def init_pipeline_parser():
         default=20,
         help="Minimum depth required to call a variant (default: %(default)d)",
     )
-    parser_minion.add_argument(
-        "--min-variant-quality",
-        type=int,
-        default=15,
-        help="Minimum quality required to call a variant (default: %(default)d)",
-    )
     parser_minion.add_argument(
         "--threads",
         type=int,

artic/vcf_filter.py

@@ -16,62 +16,59 @@ def in_frame(v):
     return False
 
 
-class NanoporeFilter:
-    def __init__(self, no_frameshifts):
+class Clair3Filter:
+    def __init__(self, no_frameshifts, min_depth):
         self.no_frameshifts = no_frameshifts
-        pass
+        self.min_depth = min_depth
+        self.min_variant_quality = 10
 
     def check_filter(self, v):
-        total_reads = float(v.INFO["TotalReads"])
         qual = v.QUAL
-        # strandbias = float(v.INFO["StrandFisherTest"])
 
-        if qual / total_reads < 3:
+        if qual < self.min_variant_quality:
             return False
 
         if self.no_frameshifts and not in_frame(v):
             return False
 
-        if v.is_indel:
-            strand_fraction_by_strand = v.INFO["SupportFractionByStrand"]
-            if float(strand_fraction_by_strand[0]) < 0.5:
-                return False
-
-            if float(strand_fraction_by_strand[1]) < 0.5:
-                return False
+        try:
+            # We don't really care about the depth here, just skip it if it isn't there
+            depth = v.INFO["DP"]
+        except KeyError:
+            depth = v.format("DP")[0][0]
 
-        if total_reads < 20:
+        if depth < self.min_depth:
             return False
 
         return True
 
 
-class MedakaFilter:
-    def __init__(self, no_frameshifts, min_depth, min_variant_quality):
-        self.no_frameshifts = no_frameshifts
-        self.min_depth = min_depth
-        self.min_variant_quality = min_variant_quality
+# class MedakaFilter:
+#     def __init__(self, no_frameshifts, min_depth):
+#         self.no_frameshifts = no_frameshifts
+#         self.min_depth = min_depth
+#         self.min_variant_quality = 20
 
-    def check_filter(self, v, min_depth):
-        try:
-            # We don't really care about the depth here, just skip it if it isn't there
-            depth = v.INFO["DP"]
-        except KeyError:
-            depth = v.format("DP")[0][0]
+#     def check_filter(self, v, min_depth):
+#         try:
+#             # We don't really care about the depth here, just skip it if it isn't there
+#             depth = v.INFO["DP"]
+#         except KeyError:
+#             depth = v.format("DP")[0][0]
 
-        if depth < min_depth:
-            return False
+#         if depth < min_depth:
+#             return False
 
-        if self.no_frameshifts and not in_frame(v):
-            return False
+#         if self.no_frameshifts and not in_frame(v):
+#             return False
 
-        if v.num_het:
-            return False
+#         if v.num_het:
+#             return False
 
-        if v.QUAL < self.min_variant_quality:
-            return False
+#         if v.QUAL < self.min_variant_quality:
+#             return False
 
-        return True
+#         return True
 
 
 def go(args):
@@ -80,7 +77,7 @@ def go(args):
     vcf_writer.write_header()
     vcf_writer_filtered = Writer(args.output_fail_vcf, vcf_reader, "w")
     vcf_writer_filtered.write_header()
-    filter = MedakaFilter(args.no_frameshifts)
+    filter = Clair3Filter(args.no_frameshifts, args.min_depth)
 
     variants = [v for v in vcf_reader]
 
@@ -101,15 +98,15 @@ def go(args):
             pass
 
         # now apply the filter to send variants to PASS or FAIL file
-        if filter.check_filter(v, args.min_depth):
+        if filter.check_filter(v):
             vcf_writer.write_record(v)
         else:
             variant_passes = False
 
             indx = "%s-%s" % (v.CHROM, v.POS)
             if len(group_variants[indx]) > 1:
                 for check_variant in group_variants[indx]:
-                    if filter.check_filter(check_variant, args.min_depth):
+                    if filter.check_filter(check_variant):
                         variant_passes = True
 
             if not variant_passes:
@@ -124,7 +121,6 @@ def main():
 
     parser = argparse.ArgumentParser()
     parser.add_argument("--no-frameshifts", action="store_true")
-    parser.add_argument("--min-variant-quality", type=int)
     parser.add_argument("--min-depth", type=int)
     parser.add_argument("inputvcf")
     parser.add_argument("output_pass_vcf")