Add min-depth and min-variant-quality arguments

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"
+                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}'"
             )
 
             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} {pre_filter_vcf}.gz {args.sample}.pass.vcf {args.sample}.fail.vcf"
+        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"
     )
 
     # 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

@@ -98,6 +98,18 @@ def init_pipeline_parser():
         default=35,
         help="Allow fuzzy primer matching within this threshold (default: %(default)d)",
     )
+    parser_minion.add_argument(
+        "--min-depth",
+        type=int,
+        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/utils.py

@@ -822,7 +822,7 @@ def choose_model(read_file: str) -> dict:
                 file=sys.stderr,
             )
         )
-        sys.exit(1)
+        sys.exit(6)
 
     if len(tags["basecall_model_version_id"].split("_")) == 4:
         molecule, pore_type, kit_id, model = tags["basecall_model_version_id"].split(
@@ -845,7 +845,7 @@ def choose_model(read_file: str) -> dict:
                     file=sys.stderr,
                 )
             )
-            sys.exit(1)
+            sys.exit(6)
 
         if len(possible_models) == 1:
             print(
@@ -868,7 +868,7 @@ def choose_model(read_file: str) -> dict:
                     file=sys.stderr,
                 )
             )
-            sys.exit(1)
+            sys.exit(6)
 
         if len(possible_models) == 1:
             print(
@@ -902,7 +902,7 @@ def choose_model(read_file: str) -> dict:
                 file=sys.stderr,
             )
         )
-        sys.exit(1)
+        sys.exit(6)
 
     if len(possible_models) == 1:
         print(
@@ -924,7 +924,7 @@ def choose_model(read_file: str) -> dict:
                 file=sys.stderr,
             )
         )
-        sys.exit(1)
+        sys.exit(6)
 
     if len(possible_models) == 1:
         print(
@@ -955,4 +955,4 @@ def choose_model(read_file: str) -> dict:
         file=sys.stderr,
     )
 
-    sys.exit(1)
+    sys.exit(6)

artic/vcf_filter.py

@@ -50,14 +50,14 @@ class MedakaFilter:
     def __init__(self, no_frameshifts):
         self.no_frameshifts = no_frameshifts
 
-    def check_filter(self, v):
+    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 < 20:
+        
+        if depth < min_depth:
             return False
 
         if self.no_frameshifts and not in_frame(v):
@@ -94,12 +94,12 @@ def go(args):
         except KeyError:
             pass
 
-        if v.QUAL < 15:
+        if v.QUAL < args.min_variant_quality:
             print(f"Suppress variant {v.POS} due to low quality")
             continue
 
         # now apply the filter to send variants to PASS or FAIL file
-        if filter.check_filter(v):
+        if filter.check_filter(v, args.min_depth):
             vcf_writer.write_record(v)
         else:
             variant_passes = False
@@ -122,6 +122,8 @@ 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")
     parser.add_argument("output_fail_vcf")

environment.yml

@@ -18,7 +18,7 @@ dependencies:
   - pysam
   - pytest 
   - cyvcf2
-  - pyfaidx
+  - pyfaidx=0.6.0
   - requests
   - samtools
   - tqdm