Mercurial > repos > miller-lab > genome_diversity
diff draw_variants.py @ 31:a631c2f6d913
Update to Miller Lab devshed revision 3c4110ffacc3
author | Richard Burhans <burhans@bx.psu.edu> |
---|---|
date | Fri, 20 Sep 2013 13:25:27 -0400 |
parents | 8997f2ca8c7a |
children | ea52b23f1141 |
line wrap: on
line diff
--- 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)