diff draw_variants.py @ 31:a631c2f6d913

Update to Miller Lab devshed revision 3c4110ffacc3

--- a/draw_variants.py	Fri Jul 26 12:51:13 2013 -0400
+++ b/draw_variants.py	Fri Sep 20 13:25:27 2013 -0400
@@ -6,35 +6,81 @@
 
 ################################################################################
 
-if len(sys.argv) != 10:
-    gd_util.die('Usage')
-
-snp_input, indel_input, coverage_input, annotation_input, indiv_input, ref_name, min_coverage, output, ind_arg = sys.argv[1:]
+def load_pop(file, wrapped_dict):
+    if file == '/dev/null':
+        pop = None
+    else:
+        pop = Population()
+        pop.from_wrapped_dict(wrapped_dict)
+    return pop
 
-p_total = Population()
-p_total.from_wrapped_dict(ind_arg)
-
-p1 = Population()
-p1.from_population_file(indiv_input)
-if not p_total.is_superset(p1):
-    gd_util.die('There is an individual in the population individuals that is not in the SNP table')
+def append_tags(the_list, p, p_type, val):
+    if p is None:
+        return
+    for tag in p.tag_list():
+        column, name = tag.split(':')
+        if p_type == 'gd_genotype':
+            column = int(column) - 2
+        the_list.append('{0}:{1}:{2}'.format(val, column, name))
 
 ################################################################################
 
-prog = 'mk_Ji'
+if len(sys.argv) != 11:
+    gd_util.die('Usage')
+
+
+snp_file, snp_ext, snp_arg, indiv_input, annotation_input, cov_file, cov_ext, cov_arg, min_coverage, output = sys.argv[1:]
+
+p_snp = load_pop(snp_file, snp_arg)
+p_cov = load_pop(cov_file, cov_arg)
+
+if indiv_input == '/dev/null':
+    if p_snp is not None:
+        p_ind = p_snp
+    elif p_cov is not None:
+        p_ind = p_cov
+    else:
+        p_ind = None
+    order_p_ind = True
+else:
+    p_ind = Population()
+    p_ind.from_population_file(indiv_input)
+    order_p_ind = False
+
+## p ind must be from either p_snp or p_cov
+if p_snp is not None and p_cov is not None:
+    if not (p_snp.is_superset(p_ind) or p_cov.is_superset(p_ind)):
+        gd_util.die('There is an individual in the population individuals that is not in the SNP/Genotype or Coverage table')
+elif p_snp is not None:
+    if not p_snp.is_superset(p_ind):
+        gd_util.die('There is an individual in the population individuals that is not in the SNP/Genotype table')
+elif p_cov is not None:
+    if not p_cov.is_superset(p_ind):
+        gd_util.die('There is an individual in the population individuals that is not in the Coverage table')
+
+
+################################################################################
+
+prog = 'mito_draw'
 
 args = [ prog ]
-args.append(snp_input)
-args.append(indel_input)
-args.append(coverage_input)
+args.append(snp_file)
+args.append(cov_file)
 args.append(annotation_input)
 args.append(min_coverage)
-args.append(ref_name)
 
-for tag in p1.tag_list():
-    args.append(tag)
+if order_p_ind:
+    for column in sorted(p_ind.column_list()):
+        individual = p_ind.individual_with_column(column)
+        name = individual.name.split()[0]
+        args.append('{0}:{1}:{2}'.format(0, column, name))
+else:
+    append_tags(args, p_ind, 'gd_indivs', 0)
 
-with open('mk_Ji.out', 'w') as fh:
+append_tags(args, p_snp, snp_ext, 1)
+append_tags(args, p_cov, cov_ext, 2)
+
+with open('Ji.spec', 'w') as fh:
     gd_util.run_program(prog, args, stdout=fh)
 
 ################################################################################
@@ -48,9 +94,20 @@
 args.append(0.3)
 args.append('-g')
 args.append(0.2)
-args.append('mk_Ji.out')
+args.append('Ji.spec')
 
-with open(output, 'w') as fh:
+with open('Ji.svg', 'w') as fh:
     gd_util.run_program(prog, args, stdout=fh)
 
+################################################################################
+
+prog = 'convert'
+
+args = [ prog ]
+args.append('-density')
+args.append(100)
+args.append('Ji.svg')
+args.append('tiff:{0}'.format(output))
+
+gd_util.run_program(prog, args)
 sys.exit(0)