#NEWLINE# def processMali( mali , options ) : #NEWLINE# #INDENT# ncols = mali . getNumColumns ( ) #NEWLINE# if ( ncols == 0 ) : #NEWLINE# #INDENT# raise 'refusing to process empty alignment.' #NEWLINE# #UNINDENT# if ( options . block_size != None ) : #NEWLINE# #INDENT# if ( options . block_size < 1 ) : #NEWLINE# #INDENT# size = ( int ( ( ( float ( ncols ) / 3.0 ) * options . block_size ) ) * 3 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# size = ( int ( options . block_size ) * 3 ) #NEWLINE# #UNINDENT# size = min ( size , ncols ) #NEWLINE# mali . addAnnotation ( 'STATE' , ( ( 'N' * size ) + ( 'C' * ( ncols - size ) ) ) ) #NEWLINE# #UNINDENT# for id in mali . getIdentifiers ( ) : #NEWLINE# #INDENT# if ( options . separator in id ) : #NEWLINE# #INDENT# species = id . split ( options . separator ) [ 0 ] #NEWLINE# mali . rename ( id , species ) #NEWLINE# #UNINDENT# #UNINDENT# map_new2old = mali . mapIdentifiers ( ) #NEWLINE# map_old2new = IOTools . getInvertedDictionary ( map_new2old , make_unique = True ) #NEWLINE# ids = mali . getIdentifiers ( ) #NEWLINE# xgram = XGram . XGram ( ) #NEWLINE# if options . xrate_min_increment : #NEWLINE# #INDENT# xgram . setMinIncrement ( options . xrate_min_increment ) #NEWLINE# #UNINDENT# ( ninput , noutput , nskipped ) = ( 0 , 0 , 0 ) #NEWLINE# if options . clean_mali : #NEWLINE# #INDENT# mali . mGapChars = ( mali . mGapChars + ( 'n' , 'N' ) ) #NEWLINE# mali . removeGaps ( minimum_gaps = 1 , frame = 3 ) #NEWLINE# #UNINDENT# if options . input_filename_tree : #NEWLINE# #INDENT# nexus = TreeTools . Newick2Nexus ( open ( options . input_filename_tree , 'r' ) ) #NEWLINE# tree = nexus . trees [ 0 ] #NEWLINE# tree . relabel ( map_old2new ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# tree = None #NEWLINE# #UNINDENT# annotation = mali . getAnnotation ( 'STATE' ) #NEWLINE# chars = set ( list ( annotation ) ) #NEWLINE# for c in chars : #NEWLINE# #INDENT# assert ( c in ( 'N' , 'C' ) ) , "unknown annotation %s: only 'N' and 'C' are recognized" #NEWLINE# #UNINDENT# if ( len ( chars ) == 1 ) : #NEWLINE# #INDENT# if ( options . loglevel >= 1 ) : #NEWLINE# #INDENT# options . stdlog . write ( '# WARNING: only a single block' ) #NEWLINE# #UNINDENT# blocks = ( ( 'B0_' , chars [ 0 ] ) , ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# blocks = ( ( 'B0_' , 'N' ) , ( 'B1_' , 'C' ) ) #NEWLINE# #UNINDENT# ( result , mali , ids ) = prepareGrammar ( xgram , mali , tree , map_old2new , blocks , options ) #NEWLINE# trained_model = result . getModel ( ) #NEWLINE# ( pis , matrices ) = RateEstimation . getRateMatrix ( trained_model ) #NEWLINE# annotation = mali . getAnnotation ( 'STATE' ) #NEWLINE# for ( block , code ) in blocks : #NEWLINE# #INDENT# terminals = ( ( '%sCOD0' % block ) , ( '%sCOD1' % block ) , ( '%sCOD2' % block ) ) #NEWLINE# pi = pis [ terminals ] #NEWLINE# if ( options . shared_rates == 'all' ) : #NEWLINE# #INDENT# rate_prefix_rs = '' #NEWLINE# rate_prefix_rn = '' #NEWLINE# rate_prefix_ri = '' #NEWLINE# rate_prefix_rv = '' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . shared_rates == 'kappa' ) : #NEWLINE# #INDENT# rate_prefix_rs = block #NEWLINE# rate_prefix_rn = block #NEWLINE# rate_prefix_ri = '' #NEWLINE# rate_prefix_rv = '' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . shared_rates == 'kappa-ds' ) : #NEWLINE# #INDENT# rate_prefix_rs = '' #NEWLINE# rate_prefix_rn = block #NEWLINE# rate_prefix_ri = '' #NEWLINE# rate_prefix_rv = '' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . shared_rates == 'omega' ) : #NEWLINE# #INDENT# rate_prefix_rs = '' #NEWLINE# rate_prefix_rn = '' #NEWLINE# rate_prefix_ri = block #NEWLINE# rate_prefix_rv = block #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . shared_rates == 'omega-ds' ) : #NEWLINE# #INDENT# rate_prefix_rs = '' #NEWLINE# rate_prefix_rn = '' #NEWLINE# rate_prefix_ri = block #NEWLINE# rate_prefix_rv = '' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . shared_rates == 'ds' ) : #NEWLINE# #INDENT# rate_prefix_rs = '' #NEWLINE# rate_prefix_rn = block #NEWLINE# rate_prefix_ri = block #NEWLINE# rate_prefix_rv = block #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# rate_prefix_rs = block #NEWLINE# rate_prefix_rn = block #NEWLINE# rate_prefix_ri = block #NEWLINE# rate_prefix_rv = block #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# if options . shared_frequencies : #NEWLINE# #INDENT# frequency_prefix = '' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# frequency_prefix = block #NEWLINE# #UNINDENT# rs = trained_model . mGrammar . getParameter ( ( '%sRs' % rate_prefix_rs ) ) #NEWLINE# rn = trained_model . mGrammar . getParameter ( ( '%sRn' % rate_prefix_rn ) ) #NEWLINE# ri = trained_model . mGrammar . getParameter ( ( '%sRi' % rate_prefix_ri ) ) #NEWLINE# rv = trained_model . mGrammar . getParameter ( ( '%sRv' % rate_prefix_rv ) ) #NEWLINE# nchars = annotation . count ( code ) #NEWLINE# msg = ( 'iter=%i Rs=%6.4f Rn=%6.4f Ri=%6.4f Rv=%6.4f' % ( result . getNumIterations ( ) , rs , rn , ri , rv ) ) #NEWLINE# try : #NEWLINE# #INDENT# ( Q , t ) = RateEstimation . getQMatrix ( pi , Rsi = ( rs * ri ) , Rsv = ( rs * rv ) , Rni = ( rn * ri ) , Rnv = ( rn * rv ) ) #NEWLINE# avg_omega = ( ( rs + rn ) / 2.0 ) #NEWLINE# ( Q0 , t0 ) = RateEstimation . getQMatrix ( pi , Rsi = ( ri * avg_omega ) , Rsv = ( rv * avg_omega ) , Rni = ( ri * avg_omega ) , Rnv = ( rv * avg_omega ) ) #NEWLINE# avg_kappa = ( ( ri + rv ) / 2.0 ) #NEWLINE# ( Q1 , t1 ) = RateEstimation . getQMatrix ( pi , Rsi = ( rs * avg_kappa ) , Rsv = ( rs * avg_kappa ) , Rni = ( rn * avg_kappa ) , Rnv = ( rn * avg_kappa ) ) #NEWLINE# ( rI , rV , rS , rN ) = RateEstimation . countSubstitutions ( pi , Q ) #NEWLINE# ( rI0 , rV0 , rS0 , rN0 ) = RateEstimation . countSubstitutions ( pi , Q0 ) #NEWLINE# ( rI1 , rV1 , rS1 , rN1 ) = RateEstimation . countSubstitutions ( pi , Q1 ) #NEWLINE# dS = ( ( rS / ( 3 * rS0 ) ) * t ) #NEWLINE# dN = ( ( rN / ( 3 * rN0 ) ) * t ) #NEWLINE# o_kappa = ( options . value_format % ( ( ( rI / rI0 ) * rV0 ) / rV ) ) #NEWLINE# o_omega = ( options . value_format % ( dN / dS ) ) #NEWLINE# o_dn = ( options . value_format % dN ) #NEWLINE# o_ds = ( options . value_format % dS ) #NEWLINE# o_rn = ( options . value_format % rN ) #NEWLINE# o_rs = ( options . value_format % rS ) #NEWLINE# o_rn0 = ( options . value_format % rN0 ) #NEWLINE# o_rs0 = ( options . value_format % rS0 ) #NEWLINE# o_t = ( options . value_format % t ) #NEWLINE# o_t0 = ( options . value_format % t0 ) #NEWLINE# #UNINDENT# except ZeroDivisionError : #NEWLINE# #INDENT# o_kappa = 'na' #NEWLINE# o_omega = 'na' #NEWLINE# o_dn = 'na' #NEWLINE# o_ds = 'na' #NEWLINE# o_rn = 'na' #NEWLINE# o_rs = 'na' #NEWLINE# o_rn0 = 'na' #NEWLINE# o_rs0 = 'na' #NEWLINE# o_t = 'na' #NEWLINE# o_t0 = 'na' #NEWLINE# Q = None #NEWLINE# msg = 'insufficient data to estimate rate matrix.' #NEWLINE# #UNINDENT# options . stdout . write ( '\t' . join ( map ( str , ( code , block , o_dn , o_ds , o_omega , 'na' , 'na' , 'na' , 'na' , o_kappa , result . getLogLikelihood ( ) , 'na' , nchars ) ) ) ) #NEWLINE# if options . with_rho : #NEWLINE# #INDENT# options . stdout . write ( ( '\t' + '\t' . join ( map ( str , ( o_rn , o_rs , o_t , o_rn0 , o_rs0 , o_t0 ) ) ) ) ) #NEWLINE# #UNINDENT# options . stdout . write ( ( '\t%s\n' % msg ) )