#NEWLINE# def embed( args ) : #NEWLINE# #INDENT# '\n    %prog embed evidencefile scaffolds.fasta contigs.fasta\n\n    Use SSPACE evidencefile to scaffold contigs into existing scaffold\n    structure, as in `scaffolds.fasta`. Contigs.fasta were used by SSPACE\n    directly to scaffold.\n\n    Rules:\n    1. Only update existing structure by embedding contigs small enough to fit.\n    2. Promote singleton contigs only if they are big (>= min_length).\n    ' #NEWLINE# p = OptionParser ( embed . __doc__ ) #NEWLINE# p . set_mingap ( default = 10 ) #NEWLINE# p . add_option ( '--min_length' , default = 200 , type = 'int' , help = 'Minimum length to consider [default: %default]' ) #NEWLINE# ( opts , args ) = p . parse_args ( args ) #NEWLINE# if ( len ( args ) != 3 ) : #NEWLINE# #INDENT# sys . exit ( ( not p . print_help ( ) ) ) #NEWLINE# #UNINDENT# ( evidencefile , scaffolds , contigs ) = args #NEWLINE# min_length = opts . min_length #NEWLINE# ( splitfasta , oagp , cagp ) = gaps ( [ scaffolds , '--split' , '--mingap={0}' . format ( opts . mingap ) ] ) #NEWLINE# agp = AGP ( cagp ) #NEWLINE# p = agp . graph #NEWLINE# ef = EvidenceFile ( evidencefile , contigs ) #NEWLINE# sizes = ef . sz #NEWLINE# q = ef . graph #NEWLINE# logging . debug ( 'Reference graph: {0}' . format ( p ) ) #NEWLINE# logging . debug ( 'Patch graph: {0}' . format ( q ) ) #NEWLINE# newagp = deepcopy ( agp ) #NEWLINE# seen = set ( ) #NEWLINE# deleted = set ( ) #NEWLINE# for a in agp : #NEWLINE# #INDENT# if a . is_gap : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# name = a . component_id #NEWLINE# object = a . object #NEWLINE# if ( name in deleted ) : #NEWLINE# #INDENT# ( ( print >> sys . stderr ) , '* Skip {0}, already embedded' . format ( name ) ) #NEWLINE# continue #NEWLINE# #UNINDENT# seen . add ( name ) #NEWLINE# ( target_name , tag ) = get_target ( p , name ) #NEWLINE# path = q . get_path ( name , target_name , tag = tag ) #NEWLINE# path_size = ( sum ( [ sizes [ x . v ] for ( x , t ) in path ] ) if path else None ) #NEWLINE# status = NO_UPDATE #NEWLINE# if ( path and ( path_size > min_length ) and ( len ( path ) > 3 ) ) : #NEWLINE# #INDENT# path = None #NEWLINE# #UNINDENT# if ( not path ) : #NEWLINE# #INDENT# ( ( print >> sys . stderr ) , name , target_name , path , path_size , status ) #NEWLINE# continue #NEWLINE# #UNINDENT# backward = False #NEWLINE# for ( x , t ) in path : #NEWLINE# #INDENT# if ( x . v in seen ) : #NEWLINE# #INDENT# ( ( print >> sys . stderr ) , '* Does not allow backward patch on {0}' . format ( x . v ) ) #NEWLINE# backward = True #NEWLINE# break #NEWLINE# #UNINDENT# #UNINDENT# if backward : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# vv = q . get_node ( name ) #NEWLINE# path . appendleft ( ( vv , tag ) ) #NEWLINE# if ( tag == '>' ) : #NEWLINE# #INDENT# path . reverse ( ) #NEWLINE# status = INSERT_BEFORE #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( target_name is None ) : #NEWLINE# #INDENT# status = INSERT_AFTER #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# target = q . get_node ( target_name ) #NEWLINE# path . append ( ( target , tag ) ) #NEWLINE# status = INSERT_BETWEEN #NEWLINE# #UNINDENT# #UNINDENT# ( ( print >> sys . stderr ) , name , target_name , path , path_size , status ) #NEWLINE# lines = path_to_agp ( q , path , object , sizes , status ) #NEWLINE# if ( status == INSERT_BEFORE ) : #NEWLINE# #INDENT# lines = lines [ : ( - 1 ) ] #NEWLINE# td = newagp . insert_lines ( name , lines , delete = True , verbose = True ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( status == INSERT_AFTER ) : #NEWLINE# #INDENT# lines = lines [ 1 : ] #NEWLINE# td = newagp . insert_lines ( name , lines , after = True , delete = True , verbose = True ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# lines = lines [ 1 : ( - 1 ) ] #NEWLINE# td = newagp . update_between ( name , target_name , lines , delete = True , verbose = True ) #NEWLINE# #UNINDENT# #UNINDENT# deleted |= td #NEWLINE# seen |= td #NEWLINE# #UNINDENT# CUTOFF = opts . min_length #NEWLINE# for ( ctg , size ) in sizes . items ( ) : #NEWLINE# #INDENT# if ( ctg in seen ) : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# if ( size < CUTOFF ) : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# newagp . append ( AGPLine . cline ( ctg , ctg , sizes , '?' ) ) #NEWLINE# #UNINDENT# newagpfile = 'embedded.agp' #NEWLINE# newagp . print_to_file ( newagpfile , index = True ) #NEWLINE# tidy ( [ newagpfile , contigs ] )