Merge pull request #1 from LUMC/devel

Devel

README.md

@@ -154,6 +154,8 @@ Options:
                                 called multiple times  [required]
   -ps, --positive-samples TEXT  Sample(s) in positive-vcf to consider. May be
                                 called multiple times  [required]
+  -s, --stats PATH              Path to output stats json file
+  -dc, --discordant PATH        Path to output discordant VCF file
   --help                        Show this message and exit.
 ```
 

vtools/cli.py

@@ -31,11 +31,26 @@ from .gcoverage import RefRecord, region_coverages
               help="Sample(s) in positive-vcf to consider. "
                    "May be called multiple times",
               required=True)
-def evaluate_cli(call_vcf, positive_vcf, call_samples, positive_samples):
+@click.option("-s", "--stats", type=click.Path(writable=True),
+              help="Path to output stats json file", default='-')
+@click.option("-dc", "--discordant", type=click.Path(writable=True),
+              help="Path to output discordant VCF file",
+              required=False)
+def evaluate_cli(call_vcf, positive_vcf, call_samples, positive_samples, stats,
+                 discordant):
     c_vcf = VCF(call_vcf, gts012=True)
     p_vcf = VCF(positive_vcf, gts012=True)
-    evaluated = site_concordancy(c_vcf, p_vcf, call_samples, positive_samples)
-    print(json.dumps(evaluated))
+    st, disc = site_concordancy(c_vcf, p_vcf, call_samples,
+                                positive_samples)
+    # Write the stats json file
+    with click.open_file(stats, 'w') as fout:
+        print(json.dumps(st), file=fout)
+
+    # If specified, write the discordant variants
+    if discordant:
+        with click.open_file(discordant, 'w') as fout:
+            for record in disc:
+                print(record, file=fout, end='')
 
 
 @click.command()

vtools/evaluate.py

@@ -60,6 +60,8 @@ def site_concordancy(call_vcf: VCF,
         "alleles_no_call": 0,
         "alleles_low_qual": 0
     }
+    discordant_count = 0
+    discordant_records = list()
     for pos_record in positive_vcf:
         d['total_sites'] += 1
         query_str = "{0}:{1}-{2}".format(
@@ -130,4 +132,9 @@ def site_concordancy(call_vcf: VCF,
                 # the same has no call in one or more of the VCF files
                 d['alleles_no_call'] += 2
 
-    return d
+            # The current variant is discordant
+            if d['alleles_discordant'] > discordant_count:
+                discordant_count = d['alleles_discordant']
+                discordant_records.append(call_record)
+
+    return d, discordant_records