Merge testing into master

Various changes required for the new harmony pipeline. The most important changes are listed below:

• Add testing, including various tests to make sure the program aborts on unsupported input
• Change the way the vcf files are compared to allow comparisons between different REF/ALT values
• Add support for comparing partial calls, e.g. A/A vs A/. is counted as 1 concordant allele, and 1 'no call'
• Add support for filtering on RGQ instead of GQ when a site is called as homref
• Write discordant variants from the 'call' VCF file to a valid VCF file

vtools/cli.py

@@ -9,6 +9,7 @@ vtools.cli
 import json
 import click
 from cyvcf2 import VCF, Writer
+import gzip
 
 from .evaluate import site_concordancy
 from .filter import FilterParams, FilterClass, Filterer
@@ -32,25 +33,40 @@ from .gcoverage import RefRecord, region_coverages
                    "May be called multiple times",
               required=True)
 @click.option("-s", "--stats", type=click.Path(writable=True),
-              help="Path to output stats json file", required=True)
-@click.option("-dc", "--discordant", type=click.Path(writable=True),
-              help="Path to output discordant VCF file",
+              help="Path to output stats json file")
+@click.option("-dvcf", "--discordant-vcf", type=click.Path(writable=True),
+              help="Path to output gzipped discordant vcf file",
               required=False)
-def evaluate_cli(call_vcf, positive_vcf, call_samples, positive_samples, stats,
-                 discordant):
+@click.option("-mq", "--min-qual", type=float,
+              help="Minimum quality of variants to consider", default=30)
+@click.option("-md", "--min-depth", type=int,
+              help="Minimum depth of variants to consider", default=0)
+def evaluate_cli(call_vcf, positive_vcf, call_samples, positive_samples,
+                 min_qual, min_depth, stats, discordant_vcf):
     c_vcf = VCF(call_vcf, gts012=True)
     p_vcf = VCF(positive_vcf, gts012=True)
     st, disc = site_concordancy(c_vcf, p_vcf, call_samples,
-                                positive_samples)
+                                positive_samples, min_qual, min_depth)
     # Write the stats json file
-    with open(stats, 'w') as fout:
-        print(json.dumps(st), file=fout)
-
-    # If specified, write the discordant variants
-    if discordant:
-        with open(discordant, 'w') as fout:
+    if stats is None:
+        print(json.dumps(st))
+    else:
+        with click.open_file(stats, 'w') as fout:
+            fout.write(json.dumps(st))
+
+    # If there were discordand records, and a discordant-vcf should be written
+    if len(disc) > 0 and discordant_vcf:
+        with click.open_file(discordant_vcf, 'w') as fout:
+            # First, we write the vcf header
+            with gzip.open(call_vcf, 'rt') as fin:
+                for line in fin:
+                    if line.startswith('#'):
+                        fout.write(line)
+                    else:
+                        break
+            # Then we write the vcf records that were discordant
             for record in disc:
-                print(record, file=fout, end='')
+                fout.write(str(record))
 
 
 @click.command()