<s> from __future__ import division , print_function <EOL> from collections import defaultdict <EOL> from vgraph . bed import BedFile <EOL> from vgraph . norm import NormalizedLocus <EOL> from vgraph . intervals import union <EOL> from vgraph . iterstuff import sort_almost_sorted , is_empty_iter <EOL> from vgraph . linearmatch import generate_graph , generate_paths , generate_genotypes , intersect_paths <EOL> def valid_alleles ( alleles ) : <EOL> return not any ( '<STR_LIT:<>' in a or '<STR_LIT:[>' in a or '<STR_LIT:]>' in a for a in alleles ) <EOL> def is_alt_genotype ( record , name ) : <EOL> sample = record . samples [ name ] <EOL> indices = sample . allele_indices <EOL> return bool ( indices and - <NUM_LIT:1> not in indices and indices . count ( <NUM_LIT:0> ) != len ( indices ) ) <EOL> def records_to_loci ( ref , records , name ) : <EOL> last_left_stop = <NUM_LIT:0> <EOL> for recnum , record in enumerate ( records ) : <EOL> if valid_alleles ( record . alleles ) and is_alt_genotype ( record , name ) : <EOL> locus = NormalizedLocus ( recnum , record , ref , name , last_left_stop ) <EOL> last_left_stop = locus . left . stop <EOL> yield locus <EOL> def informative_chromosomes ( vars ) : <EOL> if not vars . index : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> return ( chrom for chrom in vars . index if not is_empty_iter ( vars . fetch ( chrom ) ) ) <EOL> def make_bedmap ( bedfile ) : <EOL> bedmap = defaultdict ( list ) <EOL> for record in bedfile : <EOL> bedmap [ record . contig ] . append ( record ) <EOL> return bedmap <EOL> def region_filter_include ( records , include ) : <EOL> for start , stop , ( rec , inc ) in union ( [ records , include ] ) : <EOL> if inc : <EOL> for record in rec : <EOL> yield record <EOL> def region_filter_exclude ( records , exclude ) : <EOL> for start , stop , ( rec , exc ) in union ( [ records , exclude ] ) : <EOL> if not exc : <EOL> for record in rec : <EOL> yield record <EOL> def filter_gq ( records , name , min_gq ) : <EOL> for record in records : <EOL> if '<STR_LIT>' in record . format : <EOL> gq = record . samples [ name ] [ '<STR_LIT>' ] <EOL> if gq is not None and gq < min_gq : <EOL> continue <EOL> yield record <EOL> def filter_records ( records , name , args ) : <EOL> if args . min_gq is not None : <EOL> records = filter_gq ( records , name , args . min_gq ) <EOL> if args . include_filter : <EOL> include = set ( f . strip ( ) for fs in args . include_filter for f in fs . split ( '<STR_LIT:U+002C>' ) ) <EOL> records = ( record for record in records if not include . isdisjoint ( record . filter ) ) <EOL> if args . exclude_filter : <EOL> exclude = set ( f . strip ( ) for fs in args . exclude_filter for f in fs . split ( '<STR_LIT:U+002C>' ) ) <EOL> records = ( record for record in records if exclude . isdisjoint ( record . filter ) ) <EOL> return records <EOL> def variants_by_chromosome ( refs , vars , names , args , get_all = False ) : <EOL> for var in vars : <EOL> if not var . index : <EOL> raise ValueError ( '<STR_LIT>' . format ( var . filename ) ) <EOL> chroms = [ set ( informative_chromosomes ( var ) ) for var in vars ] <EOL> chroms = set . union ( * chroms ) <EOL> if args . include_regions is not None : <EOL> include = make_bedmap ( BedFile ( args . include_regions ) ) <EOL> chroms &= set ( include ) <EOL> if args . exclude_regions is not None : <EOL> exclude = make_bedmap ( BedFile ( args . exclude_regions ) ) <EOL> if args . include_file_regions : <EOL> assert len ( args . include_file_regions ) == len ( vars ) <EOL> include_files = [ make_bedmap ( BedFile ( fn ) ) for fn in args . include_file_regions ] <EOL> if args . exclude_file_regions : <EOL> assert len ( args . exclude_file_regions ) == len ( vars ) <EOL> exclude_files = [ make_bedmap ( BedFile ( fn ) ) for fn in args . exclude_file_regions ] <EOL> for chrom in chroms : <EOL> ref = refs . fetch ( chrom ) . upper ( ) <EOL> loci = [ var . fetch ( chrom ) for var in vars ] <EOL> if get_all : <EOL> all_loci = loci = [ list ( l ) for l in loci ] <EOL> loci = [ filter_records ( l , name , args ) for l , name in zip ( loci , names ) ] <EOL> if args . include_file_regions : <EOL> loci = [ region_filter_include ( l , inc [ chrom ] ) for l , inc in zip ( loci , include_files ) ] <EOL> if args . exclude_file_regions : <EOL> loci = [ region_filter_exclude ( l , exl [ chrom ] ) for l , exl in zip ( loci , exclude_files ) ] <EOL> loci = [ records_to_loci ( ref , l , name ) for name , l in zip ( names , loci ) ] <EOL> loci = [ sort_almost_sorted ( l , key = NormalizedLocus . extreme_order_key ) for l in loci ] <EOL> if args . include_regions is not None : <EOL> loci = [ region_filter_include ( l , include [ chrom ] ) for l in loci ] <EOL> if args . exclude_regions is not None : <EOL> loci = [ region_filter_exclude ( l , exclude [ chrom ] ) for l in loci ] <EOL> if get_all : <EOL> yield chrom , ref , loci , all_loci <EOL> else : <EOL> yield chrom , ref , loci <EOL> def get_superlocus_bounds ( superloci ) : <EOL> start = min ( super [ <NUM_LIT:0> ] . left . start for super in superloci if super ) <EOL> stop = max ( super [ - <NUM_LIT:1> ] . left . stop for super in superloci if super ) <EOL> return start , stop <EOL> def locus_equal_trivial ( locus1 , locus2 ) : <EOL> left1 , left2 = locus1 . left , locus2 . left <EOL> if left1 . start != left2 . start : <EOL> return False <EOL> if left1 . stop != left2 . stop : <EOL> return False <EOL> alleles1 , alleles2 = left1 . alleles , left2 . alleles <EOL> g1 = tuple ( alleles1 [ i ] for i in locus1 . allele_indices ) <EOL> g2 = tuple ( alleles2 [ i ] for i in locus2 . allele_indices ) <EOL> if not locus1 . phased or not locus2 . phased : <EOL> g1 , g2 = tuple ( sorted ( g1 ) ) , tuple ( sorted ( g2 ) ) <EOL> return g1 == g2 <EOL> def superlocus_equal_trivial ( super1 , super2 ) : <EOL> if len ( super1 ) != len ( super2 ) : <EOL> return False <EOL> for locus1 , locus2 in zip ( super1 , super2 ) : <EOL> if not locus_equal_trivial ( locus1 , locus2 ) : <EOL> return False <EOL> return True <EOL> def superlocus_equal ( ref , start , stop , super1 , super2 , debug = False ) : <EOL> if superlocus_equal_trivial ( super1 , super2 ) : <EOL> return True , '<STR_LIT:T>' <EOL> start , stop = get_superlocus_bounds ( [ super1 , super2 ] ) <EOL> try : <EOL> graph1 , constraints1 = generate_graph ( ref , start , stop , super1 , debug ) <EOL> graph2 , constraints2 = generate_graph ( ref , start , stop , super2 , debug ) <EOL> paths1 = generate_paths ( graph1 , debug = debug ) <EOL> paths2 = generate_paths ( graph2 , feasible_paths = paths1 , debug = debug ) <EOL> paths1 , paths2 = intersect_paths ( paths1 , paths2 ) <EOL> genos1 = set ( generate_genotypes ( paths1 , constraints1 , debug ) ) <EOL> genos2 = set ( generate_genotypes ( paths2 , constraints2 , debug ) ) <EOL> if genos1 . isdisjoint ( genos2 ) : <EOL> status = False , '<STR_LIT:H>' <EOL> else : <EOL> status = True , '<STR_LIT:H>' <EOL> except ValueError : <EOL> status = None , '<STR_LIT:N>' <EOL> return status </s>
