enable specification of java; utilize default java when available

Author: yunfeiguo

compare_vcf.py

@@ -43,7 +43,7 @@ def merge_results(outdir, varsim_tp, varsim_fn, vcfeval_tp,
 
 
 class VCFComparator(object):
-    def __init__(self, prefix, true_vcf, reference, regions, sample, vcfs, exclude_filtered, match_geno, log_to_file, opts):
+    def __init__(self, prefix, true_vcf, reference, regions, sample, vcfs, exclude_filtered, match_geno, log_to_file, opts, java = "java"):
         self.prefix = prefix
         self.true_vcf = true_vcf
         self.reference = reference
@@ -55,6 +55,7 @@ def __init__(self, prefix, true_vcf, reference, regions, sample, vcfs, exclude_f
         self.regions = regions
         self.opts = opts #additional options
         self.tp,self.tp_predict,self.fp,self.fn = None, None, None, None
+        self.java = java
 
     def run(self):
         '''
@@ -96,8 +97,8 @@ def get_fn(self):
         return self.fn
 
 class VarSimVCFComparator(VCFComparator):
-    def __init__(self, prefix, true_vcf, reference, regions, sample, vcfs, exclude_filtered, disallow_partial_fp, match_geno, log_to_file, opts):
-        VCFComparator.__init__(self, prefix, true_vcf, reference, regions, sample, vcfs, exclude_filtered, match_geno, log_to_file, opts)
+    def __init__(self, prefix, true_vcf, reference, regions, sample, vcfs, exclude_filtered, disallow_partial_fp, match_geno, log_to_file, opts, java = 'java'):
+        VCFComparator.__init__(self, prefix, true_vcf, reference, regions, sample, vcfs, exclude_filtered, match_geno, log_to_file, opts, java)
         self.disallow_partial_fp = disallow_partial_fp
     def get_tp_predict(self):
         '''
@@ -111,7 +112,7 @@ def run(self):
 
         :return:
         '''
-        cmd = ['java', utils.JAVA_XMX, '-jar', utils.VARSIMJAR, 'vcfcompare',
+        cmd = [self.java, utils.JAVA_XMX, '-jar', utils.VARSIMJAR, 'vcfcompare',
            '-prefix', self.prefix, '-true_vcf',
            self.true_vcf,
            '-reference', self.reference,
@@ -154,7 +155,7 @@ def run(self):
         #command example
         #rtg-tools-3.8.4-bdba5ea_install/rtg vcfeval --baseline truth.vcf.gz \
         #--calls compare1.vcf.gz -o vcfeval_split_snp -t ref.sdf --output-mode=annotate --sample xx --squash-ploidy --regions ?? \
-        cmd = ['java', utils.JAVA_XMX, '-jar', utils.RTGJAR, 'vcfeval',
+        cmd = [self.java, utils.JAVA_XMX, '-jar', utils.RTGJAR, 'vcfeval',
                '-o', self.prefix, '--baseline',
                self.true_vcf,
                '-t', self.reference,
@@ -211,7 +212,7 @@ def run(self):
                 raise Exception('{0} was not generated by vcfeval. Please check and rerun.'.format(i))
         self.tp, self.tp_predict, self.fn, self.fp = tp, tp_predict, fn, fp
 
-def generate_sdf(reference, log):
+def generate_sdf(reference, log, java = 'java'):
     '''
     take reference and generate SDF
     :param reference:
@@ -222,7 +223,7 @@ def generate_sdf(reference, log):
         LOGGER.info('{0} exists, doing nothing'.format(sdf))
         LOGGER.info('to rerun SDF generation, please remove or rename {0}'.format(sdf))
         return sdf
-    cmd = ['java', utils.JAVA_XMX, '-jar',utils.RTGJAR,'format',
+    cmd = [java, utils.JAVA_XMX, '-jar',utils.RTGJAR,'format',
            '-o', sdf, reference]
     if log:
         with utils.versatile_open(log, 'a') as logout:
@@ -237,6 +238,8 @@ def process(args):
     :param args:
     :return:
     '''
+    args.java = utils.get_java(args.java)
+    utils.check_java(args.java)
 
     # Setup logging
     FORMAT = '%(levelname)s %(asctime)-15s %(name)-20s %(message)s'
@@ -264,7 +267,7 @@ def process(args):
                sample = args.sample, vcfs = args.vcfs,
                exclude_filtered = args.exclude_filtered,
                disallow_partial_fp = args.disallow_partial_fp,
-               match_geno = args.match_geno, log_to_file= args.log_to_file, opts = args.vcfcompare_options)
+               match_geno = args.match_geno, log_to_file= args.log_to_file, opts = args.vcfcompare_options, java = args.java)
     varsim_tp, varsim_fn, varsim_fp = varsim_comparator.get_tp(), varsim_comparator.get_fn(), varsim_comparator.get_fp()
     varsim_tp = utils.sort_and_compress(varsim_tp)
     varsim_fn = utils.sort_and_compress(varsim_fn)
@@ -273,7 +276,7 @@ def process(args):
     sdf = args.sdf
     if not sdf:
         LOGGER.info("user did not supply SDF-formatted reference, trying to generate one...")
-        sdf = generate_sdf(args.reference, args.log_to_file)
+        sdf = generate_sdf(args.reference, args.log_to_file, java = args.java)
 
     '''for vcfeval
     sample column must be present, and not empty
@@ -290,25 +293,25 @@ def process(args):
                                             sample = args.sample, vcfs = [varsim_fp],
                                             exclude_filtered = args.exclude_filtered,
                                             match_geno = args.match_geno, log_to_file= args.log_to_file,
-                                          opts = args.vcfeval_options)
+                                          opts = args.vcfeval_options, java = args.java)
     vcfeval_tp, vcfeval_tp_predict = vcfeval_comparator.get_tp(), vcfeval_comparator.get_tp_predict()
     augmented_tp, augmented_fn, augmented_fp, augmented_t = merge_results(
                       outdir = args.out_dir,
                       varsim_tp = varsim_tp, varsim_fn = varsim_fn,
                       vcfeval_tp = vcfeval_tp, varsim_fp = varsim_fp, vcfeval_tp_predict = vcfeval_tp_predict)
     augmented_tp, augmented_fn, augmented_fp, augmented_t = summarize_results(os.path.join(args.out_dir,"augmented"), augmented_tp, augmented_fn, augmented_fp, augmented_t,
-                      var_types= args.var_types, sv_length= args.sv_length, regions = args.regions, bed_either = args.bed_either)
+                      var_types= args.var_types, sv_length= args.sv_length, regions = args.regions, bed_either = args.bed_either, java = args.java)
 
 
     if args.master_vcf and args.call_vcf:
-        match_false(augmented_fp, [args.call_vcf, args.master_vcf, augmented_fn], args.out_dir, args.sample, args.log_to_file, args.vcfeval_options, sdf)
-        match_false(augmented_fn, [args.call_vcf], args.out_dir, args.sample, args.log_to_file, args.vcfeval_options, sdf)
+        match_false(augmented_fp, [args.call_vcf, args.master_vcf, augmented_fn], args.out_dir, args.sample, args.log_to_file, args.vcfeval_options, sdf, args.java)
+        match_false(augmented_fn, [args.call_vcf], args.out_dir, args.sample, args.log_to_file, args.vcfeval_options, sdf, args.java)
 
     LOGGER.info("Variant comparison done.\nTrue positive: {0}\nFalse negative: {1}\nFalse positive: {2}\n".
                 format(augmented_tp, augmented_fn, augmented_fp))
 
 
-def match_false(augmented_file, files_to_pair_with, out_dir, sample, log_to_file, vcfeval_options, sdf):
+def match_false(augmented_file, files_to_pair_with, out_dir, sample, log_to_file, vcfeval_options, sdf, java = "java"):
     """Try to pair up each false call in a file (augmented_file) with a variant in the other files provided in a list (files_to_pair_with) to create an annotated version of the first file.
     By default the the first variant in the list is provided to get an AF, the 2nd to determine the simulated variant (for false positives) and the 3rd to determine if a false positive is
     a pure false positive (not simulated) or not (wrong genotype)"""
@@ -350,7 +353,7 @@ def match_false(augmented_file, files_to_pair_with, out_dir, sample, log_to_file
                                          exclude_filtered = False,
                                          match_geno = False,
                                          log_to_file= log_to_file,
-                                         opts = vcfeval_options)
+                                         opts = vcfeval_options, java = java)
 
                         equivalent_variant = utils.get_equivalent_variant(line_split, vcfeval_comparator.get_tp())
 
@@ -443,7 +446,7 @@ def parse_jsons(jsonfile, stats, count_sv = False, count_all = False):
                         print ("error in {}. No {} field".format(jsonfile, err))
                         stats[vt][mt] += 0
 
-def summarize_results(prefix, tp, fn, fp, t, var_types, sv_length = 100, regions = None, bed_either = False):
+def summarize_results(prefix, tp, fn, fp, t, var_types, sv_length = 100, regions = None, bed_either = False, java = 'java'):
     '''
     count variants by type and tabulate
     :param augmented_tp:
@@ -452,7 +455,7 @@ def summarize_results(prefix, tp, fn, fp, t, var_types, sv_length = 100, regions
     :param augmented_t:
     :return:
     '''
-    cmd = ['java', utils.JAVA_XMX, '-jar', utils.VARSIMJAR, 'vcfcompareresultsparser',
+    cmd = [java, utils.JAVA_XMX, '-jar', utils.VARSIMJAR, 'vcfcompareresultsparser',
            '-prefix', prefix, '-tp',tp,
            '-fn', fn, '-fp', fp,
            '-t', t, 
@@ -492,8 +495,6 @@ def summarize_results(prefix, tp, fn, fp, t, var_types, sv_length = 100, regions
 
 
 if __name__ == "__main__":
-    utils.check_java()
-
     main_parser = argparse.ArgumentParser(description="VarSim: A high-fidelity simulation validation framework",
                                           formatter_class=argparse.ArgumentDefaultsHelpFormatter)
     main_parser.add_argument("--reference", metavar="FASTA", help="reference filename", required=True, type=str)
@@ -520,6 +521,7 @@ def summarize_results(prefix, tp, fn, fp, t, var_types, sv_length = 100, regions
     main_parser.add_argument("--vcfeval_options", metavar="OPT", help="additional options for RTG vcfeval", default="", type = str)
     main_parser.add_argument("--bed_either", action = 'store_true', help="Use either break-end of the variant for filtering instead of both")
     main_parser.add_argument("--java_max_mem", metavar="XMX", help="max java memory", default="10g", type = str)
+    main_parser.add_argument("--java", metavar="PATH", help="path to java", default="java", type = str)
 
     args = main_parser.parse_args()
     process(args)

utils.py

@@ -11,11 +11,20 @@
 SORT_VCF = os.path.realpath(os.path.join(MY_DIR, "src","sort_vcf.sh"))
 BGZIP = os.path.realpath(os.path.join(MY_DIR, "opt","htslib-1.9_install/bin/bgzip"))
 JAVA_XMX = "-Xmx"
-
+DEFAULT_JAVA = os.path.realpath(os.path.join(MY_DIR, "opt",
+                                "jdk1.8.0_131", "bin", "java"))
 COMBINE_KEEP_ALL_DUPLICATE = 1
 COMBINE_KEEP_FIRST_DUPLICATE = 2
 COMBINE_KEEP_NO_DUPLICATE = 3
 
+def get_java(java = "java"):
+    '''
+    return default java if it exists, otherwise use user-specificed version
+    '''
+    if os.path.isfile(DEFAULT_JAVA):
+        return DEFAULT_JAVA
+    return java
+
 def count_variants(vcf):
     '''
     count number of variants
@@ -29,10 +38,10 @@ def count_variants(vcf):
                 count += 1
     return count
 
-def check_java():
+def check_java(java="java"):
     logger = logging.getLogger(check_java.__name__)
     try:
-        jv = subprocess.check_output("java -Xmx100m -version", stderr=subprocess.STDOUT, shell=True)
+        jv = subprocess.check_output("{} -Xmx100m -version".format(java), stderr=subprocess.STDOUT, shell=True)
         if "openjdk" in jv or "OpenJDK" in jv:
             raise EnvironmentError("Please replace OpenJDK with Oracle JDK")
         jv = filter(lambda x: x.startswith("java version"), jv.split("\n"))[0].split()[2].replace("\"", "")
@@ -42,8 +51,9 @@ def check_java():
     except subprocess.CalledProcessError:
         raise EnvironmentError("No java (>=1.8) found")
 
-def get_version():
-    return subprocess.check_output("java -jar {} -version".format(VARSIMJAR), shell=True).strip()
+def get_version(java="java"):
+    java = get_java(java)
+    return subprocess.check_output("{} -jar {} -version".format(java, VARSIMJAR), shell=True).strip()
 
 def run_shell_command(cmd, cmd_stdout, cmd_stderr, cmd_dir="."):
     '''