<s> """<STR_LIT>""" <EOL> import sys <EOL> import os <EOL> import os . path as op <EOL> import csv <EOL> import logging <EOL> import numpy as np <EOL> from math import log , sqrt , pi , exp <EOL> from itertools import product , combinations <EOL> from functools import partial <EOL> from Bio import SeqIO <EOL> from Bio import AlignIO <EOL> from Bio . Align . Applications import ClustalwCommandline , MuscleCommandline <EOL> from jcvi . formats . base import LineFile , must_open <EOL> from jcvi . graphics . base import plt , savefig , AbstractLayout , markup <EOL> from jcvi . utils . table import write_csv <EOL> from jcvi . utils . cbook import gene_name <EOL> from jcvi . apps . base import OptionParser , ActionDispatcher , mkdir , sh , Popen , getpath , iglob <EOL> CLUSTALW_BIN = partial ( getpath , name = "<STR_LIT>" , warn = "<STR_LIT>" ) <EOL> MUSCLE_BIN = partial ( getpath , name = "<STR_LIT>" , warn = "<STR_LIT>" ) <EOL> PAL2NAL_BIN = partial ( getpath , name = "<STR_LIT>" , warn = "<STR_LIT>" ) <EOL> PAML_BIN = partial ( getpath , name = "<STR_LIT>" , warn = "<STR_LIT>" ) <EOL> class AbstractCommandline : <EOL> def run ( self ) : <EOL> r = Popen ( str ( self ) ) <EOL> return r . communicate ( ) <EOL> class YnCommandline ( AbstractCommandline ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , ctl_file , command = PAML_BIN ( "<STR_LIT>" ) ) : <EOL> self . ctl_file = ctl_file <EOL> self . parameters = [ ] <EOL> self . command = command <EOL> def __str__ ( self ) : <EOL> return self . command + "<STR_LIT>" % self . ctl_file <EOL> class MrTransCommandline ( AbstractCommandline ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , prot_align_file , nuc_file , output_file , outfmt = "<STR_LIT>" , <EOL> command = PAL2NAL_BIN ( "<STR_LIT>" ) ) : <EOL> self . prot_align_file = prot_align_file <EOL> self . nuc_file = nuc_file <EOL> self . output_file = output_file <EOL> self . outfmt = outfmt <EOL> self . command = command <EOL> self . parameters = [ ] <EOL> def __str__ ( self ) : <EOL> return self . command + "<STR_LIT>" % ( self . prot_align_file , self . nuc_file , self . outfmt , self . output_file ) <EOL> def main ( ) : <EOL> actions = ( <EOL> ( '<STR_LIT>' , '<STR_LIT>' ) , <EOL> ( '<STR_LIT>' , '<STR_LIT>' ) , <EOL> ( '<STR_LIT>' , '<STR_LIT>' ) , <EOL> ( '<STR_LIT>' , '<STR_LIT>' ) , <EOL> ( '<STR_LIT>' , '<STR_LIT>' ) , <EOL> ( '<STR_LIT>' , '<STR_LIT>' ) , <EOL> ( '<STR_LIT>' , '<STR_LIT>' ) , <EOL> ( '<STR_LIT>' , '<STR_LIT>' ) , <EOL> ) <EOL> p = ActionDispatcher ( actions ) <EOL> p . dispatch ( globals ( ) ) <EOL> def batch ( args ) : <EOL> """<STR_LIT>""" <EOL> from jcvi . apps . grid import MakeManager <EOL> p = OptionParser ( batch . __doc__ ) <EOL> opts , args = p . parse_args ( args ) <EOL> if len ( args ) < <NUM_LIT:2> : <EOL> sys . exit ( not p . print_help ( ) ) <EOL> cdsfile = args [ <NUM_LIT:0> ] <EOL> anchors = args [ <NUM_LIT:1> : ] <EOL> workdirs = [ "<STR_LIT:.>" . join ( op . basename ( x ) . split ( "<STR_LIT:.>" ) [ : <NUM_LIT:2> ] ) for x in anchors ] <EOL> for wd in workdirs : <EOL> mkdir ( wd ) <EOL> mm = MakeManager ( ) <EOL> for wd , ac in zip ( workdirs , anchors ) : <EOL> pairscdsfile = wd + "<STR_LIT>" <EOL> cmd = "<STR_LIT>" . format ( ac , cdsfile , pairscdsfile ) <EOL> mm . add ( ( ac , cdsfile ) , pairscdsfile , cmd ) <EOL> ksfile = wd + "<STR_LIT>" <EOL> cmd = "<STR_LIT>" . format ( pairscdsfile , ksfile , wd ) <EOL> mm . add ( pairscdsfile , ksfile , cmd ) <EOL> mm . write ( ) <EOL> class LayoutLine ( object ) : <EOL> def __init__ ( self , row , delimiter = '<STR_LIT:U+002C>' ) : <EOL> args = row . rstrip ( ) . split ( delimiter ) <EOL> args = [ x . strip ( ) for x in args ] <EOL> self . ksfile = args [ <NUM_LIT:0> ] <EOL> self . components = int ( args [ <NUM_LIT:1> ] ) <EOL> self . label = args [ <NUM_LIT:2> ] <EOL> self . color = args [ <NUM_LIT:3> ] <EOL> self . marker = args [ <NUM_LIT:4> ] <EOL> def __str__ ( self ) : <EOL> return "<STR_LIT:U+002CU+0020>" . join ( str ( x ) for x in ( self . ksfile , self . components , <EOL> self . label , self . color , self . marker ) ) <EOL> class Layout ( AbstractLayout ) : <EOL> def __init__ ( self , filename , delimiter = '<STR_LIT:U+002C>' ) : <EOL> super ( Layout , self ) . __init__ ( filename ) <EOL> if not op . exists ( filename ) : <EOL> ksfiles = iglob ( "<STR_LIT:.>" , "<STR_LIT>" ) <EOL> header = "<STR_LIT>" . split ( "<STR_LIT:|>" ) <EOL> contents = [ ] <EOL> for ksfile in ksfiles : <EOL> leg = op . basename ( ksfile ) . rsplit ( "<STR_LIT:.>" , <NUM_LIT:1> ) [ <NUM_LIT:0> ] <EOL> if leg . count ( "<STR_LIT:.>" ) == <NUM_LIT:1> : <EOL> leg = leg . replace ( "<STR_LIT:.>" , "<STR_LIT>" ) <EOL> contents . append ( ( ksfile , "<STR_LIT:1>" , leg , "<STR_LIT>" , "<STR_LIT>" ) ) <EOL> write_csv ( header , contents , comment = True , filename = filename ) <EOL> fp = open ( filename ) <EOL> for row in fp : <EOL> if row [ <NUM_LIT:0> ] == '<STR_LIT:#>' : <EOL> continue <EOL> self . append ( LayoutLine ( row , delimiter = delimiter ) ) <EOL> self . assign_colors ( ) <EOL> self . assign_markers ( ) <EOL> class KsPlot ( object ) : <EOL> def __init__ ( self , ax , ks_max , bins , legendp = '<STR_LIT>' ) : <EOL> self . ax = ax <EOL> self . ks_max = ks_max <EOL> self . interval = ks_max / bins <EOL> self . legendp = legendp <EOL> self . lines = [ ] <EOL> self . labels = [ ] <EOL> def add_data ( self , data , components = <NUM_LIT:1> , label = "<STR_LIT>" , <EOL> color = '<STR_LIT:r>' , marker = '<STR_LIT:.>' , fill = False , fitted = True ) : <EOL> ax = self . ax <EOL> ks_max = self . ks_max <EOL> interval = self . interval <EOL> line , line_mixture = plot_ks_dist ( ax , data , interval , components , ks_max , <EOL> color = color , marker = marker , <EOL> fill = fill , fitted = fitted ) <EOL> self . lines . append ( line ) <EOL> self . labels . append ( label ) <EOL> if fitted : <EOL> self . lines . append ( line_mixture ) <EOL> self . labels . append ( label + "<STR_LIT>" ) <EOL> self . labels = [ markup ( x ) for x in self . labels ] <EOL> def draw ( self , title = "<STR_LIT>" , filename = "<STR_LIT>" ) : <EOL> ax = self . ax <EOL> ks_max = self . ks_max <EOL> lines = self . lines <EOL> labels = self . labels <EOL> legendp = self . legendp <EOL> if len ( lines ) > <NUM_LIT:1> : <EOL> leg = ax . legend ( lines , labels , loc = legendp , <EOL> shadow = True , fancybox = True , prop = { "<STR_LIT:size>" : <NUM_LIT:10> } ) <EOL> leg . get_frame ( ) . set_alpha ( <NUM_LIT> ) <EOL> ax . set_xlim ( ( <NUM_LIT:0> , ks_max - self . interval ) ) <EOL> ax . set_title ( markup ( title ) , fontweight = "<STR_LIT>" ) <EOL> ax . set_xlabel ( markup ( '<STR_LIT>' ) ) <EOL> ax . set_ylabel ( '<STR_LIT>' ) <EOL> ax . set_xticklabels ( ax . get_xticks ( ) , family = '<STR_LIT>' ) <EOL> ax . set_yticklabels ( ax . get_yticks ( ) , family = '<STR_LIT>' ) <EOL> savefig ( filename , dpi = <NUM_LIT> ) <EOL> def multireport ( args ) : <EOL> """<STR_LIT>""" <EOL> p = OptionParser ( multireport . __doc__ ) <EOL> p . set_outfile ( outfile = "<STR_LIT>" ) <EOL> add_plot_options ( p ) <EOL> opts , args , iopts = p . set_image_options ( args , figsize = "<STR_LIT>" ) <EOL> if len ( args ) != <NUM_LIT:1> : <EOL> sys . exit ( not p . print_help ( ) ) <EOL> layoutfile , = args <EOL> ks_min = opts . vmin <EOL> ks_max = opts . vmax <EOL> bins = opts . bins <EOL> fill = opts . fill <EOL> layout = Layout ( layoutfile ) <EOL> print >> sys . stderr , layout <EOL> fig = plt . figure ( <NUM_LIT:1> , ( iopts . w , iopts . h ) ) <EOL> ax = fig . add_axes ( [ <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ] ) <EOL> kp = KsPlot ( ax , ks_max , bins , legendp = opts . legendp ) <EOL> for lo in layout : <EOL> data = KsFile ( lo . ksfile ) <EOL> data = [ x . ng_ks for x in data ] <EOL> data = [ x for x in data if ks_min <= x <= ks_max ] <EOL> kp . add_data ( data , lo . components , label = lo . label , color = lo . color , marker = lo . marker , <EOL> fill = fill , fitted = opts . fit ) <EOL> kp . draw ( title = opts . title , filename = opts . outfile ) <EOL> def get_GC3 ( cdsfile ) : <EOL> from jcvi . formats . fasta import Fasta <EOL> f = Fasta ( cdsfile , lazy = True ) <EOL> GC3 = { } <EOL> for name , rec in f . iteritems_ordered ( ) : <EOL> positions = rec . seq [ <NUM_LIT:2> : : <NUM_LIT:3> ] . upper ( ) <EOL> gc_counts = sum ( <NUM_LIT:1> for x in positions if x in "<STR_LIT>" ) <EOL> gc_ratio = gc_counts * <NUM_LIT:1.> / len ( positions ) <EOL> GC3 [ name ] = gc_ratio <EOL> return GC3 <EOL> def plot_GC3 ( GC3 , cdsfile , fill = "<STR_LIT>" ) : <EOL> from jcvi . graphics . histogram import histogram <EOL> numberfile = "<STR_LIT>" . format ( cdsfile ) <EOL> fw = must_open ( numberfile , "<STR_LIT:w>" ) <EOL> fw . write ( "<STR_LIT:\n>" . join ( map ( str , GC3 . values ( ) ) ) ) <EOL> fw . close ( ) <EOL> histogram ( numberfile , vmin = <NUM_LIT:0> , vmax = <NUM_LIT:1> , xlabel = "<STR_LIT>" , title = cdsfile , <EOL> bins = <NUM_LIT:50> , skip = <NUM_LIT:0> , ascii = False , log = <NUM_LIT:0> , fill = fill ) <EOL> logging . debug ( "<STR_LIT>" . format ( len ( GC3 ) , numberfile ) ) <EOL> def gc3 ( args ) : <EOL> """<STR_LIT>""" <EOL> p = OptionParser ( gc3 . __doc__ ) <EOL> p . add_option ( "<STR_LIT>" , default = False , action = "<STR_LIT:store_true>" , <EOL> help = "<STR_LIT>" ) <EOL> p . set_outfile ( ) <EOL> opts , args = p . parse_args ( args ) <EOL> outfile = opts . outfile <EOL> plot = opts . plot <EOL> if not <NUM_LIT:1> < len ( args ) < <NUM_LIT:4> : <EOL> sys . exit ( not p . print_help ( ) ) <EOL> ks_file , cdsfile = args [ : <NUM_LIT:2> ] <EOL> GC3 = get_GC3 ( cdsfile ) <EOL> if plot : <EOL> plot_GC3 ( GC3 , cdsfile , fill = "<STR_LIT>" ) <EOL> if len ( args ) == <NUM_LIT:3> : <EOL> cdsfile2 = args [ <NUM_LIT:2> ] <EOL> GC3_2 = get_GC3 ( cdsfile2 ) <EOL> GC3 . update ( GC3_2 ) <EOL> if plot : <EOL> plot_GC3 ( GC3_2 , cdsfile2 , fill = "<STR_LIT>" ) <EOL> data = KsFile ( ks_file ) <EOL> noriginals = len ( data ) <EOL> fw = must_open ( outfile , "<STR_LIT:w>" ) <EOL> writer = csv . writer ( fw ) <EOL> writer . writerow ( fields . split ( "<STR_LIT:U+002C>" ) ) <EOL> nlines = <NUM_LIT:0> <EOL> cutoff = <NUM_LIT> <EOL> for d in data : <EOL> a , b = d . name . split ( "<STR_LIT:;>" ) <EOL> aratio , bratio = GC3 [ a ] , GC3 [ b ] <EOL> if ( aratio + bratio ) / <NUM_LIT:2> > cutoff : <EOL> continue <EOL> writer . writerow ( d ) <EOL> nlines += <NUM_LIT:1> <EOL> logging . debug ( "<STR_LIT>" . format ( nlines , noriginals ) ) <EOL> def extract_pairs ( abed , bbed , groups ) : <EOL> """<STR_LIT>""" <EOL> agenome = op . basename ( abed . filename ) . split ( "<STR_LIT:.>" ) [ <NUM_LIT:0> ] <EOL> bgenome = op . basename ( bbed . filename ) . split ( "<STR_LIT:.>" ) [ <NUM_LIT:0> ] <EOL> aorder = abed . order <EOL> border = bbed . order <EOL> pairsfile = "<STR_LIT>" . format ( agenome , bgenome ) <EOL> fw = open ( pairsfile , "<STR_LIT:w>" ) <EOL> is_self = abed . filename == bbed . filename <EOL> npairs = <NUM_LIT:0> <EOL> for group in groups : <EOL> iter = combinations ( group , <NUM_LIT:2> ) if is_self else product ( group , repeat = <NUM_LIT:2> ) <EOL> for a , b in iter : <EOL> if a not in aorder or b not in border : <EOL> continue <EOL> print >> fw , "<STR_LIT:\t>" . join ( ( a , b ) ) <EOL> npairs += <NUM_LIT:1> <EOL> logging . debug ( "<STR_LIT>" . format ( pairsfile , npairs ) ) <EOL> def fromgroups ( args ) : <EOL> """<STR_LIT>""" <EOL> from jcvi . formats . bed import Bed <EOL> p = OptionParser ( fromgroups . __doc__ ) <EOL> opts , args = p . parse_args ( args ) <EOL> if len ( args ) < <NUM_LIT:2> : <EOL> sys . exit ( not p . print_help ( ) ) <EOL> groupsfile = args [ <NUM_LIT:0> ] <EOL> bedfiles = args [ <NUM_LIT:1> : ] <EOL> beds = [ Bed ( x ) for x in bedfiles ] <EOL> fp = open ( groupsfile ) <EOL> groups = [ row . strip ( ) . split ( "<STR_LIT:U+002C>" ) for row in fp ] <EOL> for b1 , b2 in product ( beds , repeat = <NUM_LIT:2> ) : <EOL> extract_pairs ( b1 , b2 , groups ) <EOL> def find_first_isoform ( a , f ) : <EOL> if a in f : <EOL> return a <EOL> for i in xrange ( <NUM_LIT:100> ) : <EOL> ia = "<STR_LIT:.>" . join ( ( a , str ( i ) ) ) <EOL> if ia in f : <EOL> return ia <EOL> return a <EOL> def prepare ( args ) : <EOL> """<STR_LIT>""" <EOL> from jcvi . formats . fasta import Fasta <EOL> p = OptionParser ( prepare . __doc__ ) <EOL> p . set_outfile ( ) <EOL> opts , args = p . parse_args ( args ) <EOL> outfile = opts . outfile <EOL> if len ( args ) == <NUM_LIT:2> : <EOL> pairsfile , cdsfile = args <EOL> pepfile = None <EOL> elif len ( args ) == <NUM_LIT:3> : <EOL> pairsfile , cdsfile , pepfile = args <EOL> else : <EOL> sys . exit ( not p . print_help ( ) ) <EOL> f = Fasta ( cdsfile ) <EOL> fp = open ( pairsfile ) <EOL> fw = must_open ( outfile , "<STR_LIT:w>" ) <EOL> if pepfile : <EOL> assert outfile != "<STR_LIT>" , "<STR_LIT>" <EOL> f2 = Fasta ( pepfile ) <EOL> fw2 = must_open ( outfile + "<STR_LIT>" , "<STR_LIT:w>" ) <EOL> for row in fp : <EOL> if row [ <NUM_LIT:0> ] == '<STR_LIT:#>' : <EOL> continue <EOL> a , b = row . split ( ) [ : <NUM_LIT:2> ] <EOL> if a == b : <EOL> logging . debug ( "<STR_LIT>" . format ( a , b ) ) <EOL> continue <EOL> if a not in f : <EOL> a = find_first_isoform ( a , f ) <EOL> assert a , a <EOL> if b not in f : <EOL> b = find_first_isoform ( b , f ) <EOL> assert b , b <EOL> acds = f [ a ] <EOL> bcds = f [ b ] <EOL> SeqIO . write ( ( acds , bcds ) , fw , "<STR_LIT>" ) <EOL> if pepfile : <EOL> apep = f2 [ a ] <EOL> bpep = f2 [ b ] <EOL> SeqIO . write ( ( apep , bpep ) , fw2 , "<STR_LIT>" ) <EOL> fw . close ( ) <EOL> if pepfile : <EOL> fw2 . close ( ) <EOL> def calc ( args ) : <EOL> """<STR_LIT>""" <EOL> from jcvi . formats . fasta import translate <EOL> p = OptionParser ( calc . __doc__ ) <EOL> p . add_option ( "<STR_LIT>" , action = "<STR_LIT:store_true>" , <EOL> help = "<STR_LIT>" "<STR_LIT>" ) <EOL> p . add_option ( "<STR_LIT>" , default = "<STR_LIT>" , choices = ( "<STR_LIT>" , "<STR_LIT>" ) , <EOL> help = "<STR_LIT>" ) <EOL> p . add_option ( "<STR_LIT>" , default = os . getcwd ( ) , help = "<STR_LIT>" ) <EOL> p . set_outfile ( ) <EOL> opts , args = p . parse_args ( args ) <EOL> if len ( args ) == <NUM_LIT:1> : <EOL> protein_file , dna_file = None , args [ <NUM_LIT:0> ] <EOL> elif len ( args ) == <NUM_LIT:2> : <EOL> protein_file , dna_file = args <EOL> else : <EOL> print >> sys . stderr , "<STR_LIT>" <EOL> sys . exit ( not p . print_help ( ) ) <EOL> output_h = must_open ( opts . outfile , "<STR_LIT:w>" ) <EOL> print >> output_h , fields <EOL> work_dir = op . join ( opts . workdir , "<STR_LIT>" ) <EOL> mkdir ( work_dir ) <EOL> if not protein_file : <EOL> protein_file = dna_file + "<STR_LIT>" <EOL> translate_args = [ dna_file , "<STR_LIT>" + protein_file ] <EOL> if opts . longest : <EOL> translate_args += [ "<STR_LIT>" ] <EOL> dna_file , protein_file = translate ( translate_args ) <EOL> prot_iterator = SeqIO . parse ( open ( protein_file ) , "<STR_LIT>" ) <EOL> dna_iterator = SeqIO . parse ( open ( dna_file ) , "<STR_LIT>" ) <EOL> for p_rec_1 , p_rec_2 , n_rec_1 , n_rec_2 in zip ( prot_iterator , prot_iterator , dna_iterator , dna_iterator ) : <EOL> print >> sys . stderr , "<STR_LIT>" , p_rec_1 . name , p_rec_2 . name <EOL> if opts . msa == "<STR_LIT>" : <EOL> align_fasta = clustal_align_protein ( ( p_rec_1 , p_rec_2 ) , work_dir ) <EOL> elif opts . msa == "<STR_LIT>" : <EOL> align_fasta = muscle_align_protein ( ( p_rec_1 , p_rec_2 ) , work_dir ) <EOL> mrtrans_fasta = run_mrtrans ( align_fasta , ( n_rec_1 , n_rec_2 ) , work_dir ) <EOL> if mrtrans_fasta : <EOL> ds_subs_yn , dn_subs_yn , ds_subs_ng , dn_subs_ng = find_synonymous ( mrtrans_fasta , work_dir ) <EOL> if ds_subs_yn is not None : <EOL> pair_name = "<STR_LIT>" % ( p_rec_1 . name , p_rec_2 . name ) <EOL> output_h . write ( "<STR_LIT>" % ( "<STR_LIT:U+002C>" . join ( str ( x ) for x in ( pair_name , <EOL> ds_subs_yn , dn_subs_yn , ds_subs_ng , dn_subs_ng ) ) ) ) <EOL> output_h . flush ( ) <EOL> sh ( "<STR_LIT>" ) <EOL> def find_synonymous ( input_file , work_dir ) : <EOL> """<STR_LIT>""" <EOL> ctl_file = op . join ( work_dir , "<STR_LIT>" ) <EOL> output_file = op . join ( work_dir , "<STR_LIT>" ) <EOL> ctl_h = open ( ctl_file , "<STR_LIT:w>" ) <EOL> ctl_h . write ( "<STR_LIT>" % <EOL> ( input_file , output_file ) ) <EOL> ctl_h . write ( "<STR_LIT>" ) <EOL> ctl_h . close ( ) <EOL> cl = YnCommandline ( ctl_file ) <EOL> print >> sys . stderr , "<STR_LIT>" , cl <EOL> r , e = cl . run ( ) <EOL> ds_value_yn = None <EOL> ds_value_ng = None <EOL> dn_value_yn = None <EOL> dn_value_ng = None <EOL> output_h = open ( output_file ) <EOL> row = output_h . readline ( ) <EOL> while row : <EOL> if row . find ( "<STR_LIT>" ) >= <NUM_LIT:0> : <EOL> for x in xrange ( <NUM_LIT:5> ) : <EOL> row = output_h . next ( ) <EOL> dn_value_ng , ds_value_ng = row . split ( '<STR_LIT:(>' ) [ <NUM_LIT:1> ] . split ( '<STR_LIT:)>' ) [ <NUM_LIT:0> ] . split ( ) <EOL> break <EOL> row = output_h . readline ( ) <EOL> output_h . close ( ) <EOL> output_h = open ( output_file ) <EOL> for line in output_h : <EOL> if line . find ( "<STR_LIT>" ) >= <NUM_LIT:0> and line . find ( "<STR_LIT>" ) == - <NUM_LIT:1> : <EOL> parts = line . split ( "<STR_LIT>" ) <EOL> ds_value_yn = extract_subs_value ( parts [ <NUM_LIT:1> ] ) <EOL> dn_value_yn = extract_subs_value ( parts [ <NUM_LIT:0> ] ) <EOL> if ds_value_yn is None or ds_value_ng is None : <EOL> h = open ( output_file ) <EOL> print >> sys . stderr , "<STR_LIT>" % h . read ( ) <EOL> return ds_value_yn , dn_value_yn , ds_value_ng , dn_value_ng <EOL> def extract_subs_value ( text ) : <EOL> """<STR_LIT>""" <EOL> parts = text . split ( ) <EOL> value = float ( parts [ - <NUM_LIT:1> ] ) <EOL> return value <EOL> def run_mrtrans ( align_fasta , recs , work_dir , outfmt = "<STR_LIT>" ) : <EOL> """<STR_LIT>""" <EOL> align_file = op . join ( work_dir , "<STR_LIT>" ) <EOL> nuc_file = op . join ( work_dir , "<STR_LIT>" ) <EOL> output_file = op . join ( work_dir , "<STR_LIT>" ) <EOL> align_h0 = open ( align_file + "<STR_LIT:0>" , "<STR_LIT:w>" ) <EOL> align_h0 . write ( str ( align_fasta ) ) <EOL> align_h0 . close ( ) <EOL> prot_seqs = { } <EOL> i = <NUM_LIT:0> <EOL> for rec in SeqIO . parse ( align_h0 . name , "<STR_LIT>" ) : <EOL> prot_seqs [ i ] = rec . seq <EOL> i += <NUM_LIT:1> <EOL> align_h = open ( align_file , "<STR_LIT:w>" ) <EOL> for i , rec in enumerate ( recs ) : <EOL> if len ( rec . id ) > <NUM_LIT:30> : <EOL> rec . id = rec . id [ : <NUM_LIT> ] + "<STR_LIT:_>" + str ( i ) <EOL> rec . description = "<STR_LIT>" <EOL> print >> align_h , "<STR_LIT>" . format ( rec . id , prot_seqs [ i ] ) <EOL> align_h . close ( ) <EOL> SeqIO . write ( recs , file ( nuc_file , "<STR_LIT:w>" ) , "<STR_LIT>" ) <EOL> cl = MrTransCommandline ( align_file , nuc_file , output_file , outfmt = outfmt ) <EOL> r , e = cl . run ( ) <EOL> if e is None : <EOL> print >> sys . stderr , "<STR_LIT>" , cl <EOL> return output_file <EOL> elif e . read ( ) . find ( "<STR_LIT>" ) >= <NUM_LIT:0> : <EOL> print >> sys . stderr , "<STR_LIT>" <EOL> return None <EOL> def clustal_align_protein ( recs , work_dir , outfmt = "<STR_LIT>" ) : <EOL> """<STR_LIT>""" <EOL> fasta_file = op . join ( work_dir , "<STR_LIT>" ) <EOL> align_file = op . join ( work_dir , "<STR_LIT>" ) <EOL> SeqIO . write ( recs , file ( fasta_file , "<STR_LIT:w>" ) , "<STR_LIT>" ) <EOL> clustal_cl = ClustalwCommandline ( cmd = CLUSTALW_BIN ( "<STR_LIT>" ) , <EOL> infile = fasta_file , outfile = align_file , outorder = "<STR_LIT>" , <EOL> type = "<STR_LIT>" ) <EOL> stdout , stderr = clustal_cl ( ) <EOL> aln_file = file ( clustal_cl . outfile ) <EOL> alignment = AlignIO . read ( aln_file , "<STR_LIT>" ) <EOL> print >> sys . stderr , "<STR_LIT>" % clustal_cl <EOL> if outfmt == "<STR_LIT>" : <EOL> return alignment . format ( "<STR_LIT>" ) <EOL> if outfmt == "<STR_LIT>" : <EOL> return alignment <EOL> def muscle_align_protein ( recs , work_dir , outfmt = "<STR_LIT>" , inputorder = True ) : <EOL> """<STR_LIT>""" <EOL> fasta_file = op . join ( work_dir , "<STR_LIT>" ) <EOL> align_file = op . join ( work_dir , "<STR_LIT>" ) <EOL> SeqIO . write ( recs , file ( fasta_file , "<STR_LIT:w>" ) , "<STR_LIT>" ) <EOL> muscle_cl = MuscleCommandline ( cmd = MUSCLE_BIN ( "<STR_LIT>" ) , <EOL> input = fasta_file , out = align_file , seqtype = "<STR_LIT>" , <EOL> clwstrict = True ) <EOL> stdout , stderr = muscle_cl ( ) <EOL> alignment = AlignIO . read ( muscle_cl . out , "<STR_LIT>" ) <EOL> if inputorder : <EOL> try : <EOL> muscle_inputorder ( muscle_cl . input , muscle_cl . out ) <EOL> except ValueError : <EOL> return "<STR_LIT>" <EOL> alignment = AlignIO . read ( muscle_cl . out , "<STR_LIT>" ) <EOL> print >> sys . stderr , "<STR_LIT>" % muscle_cl <EOL> if outfmt == "<STR_LIT>" : <EOL> return alignment . format ( "<STR_LIT>" ) <EOL> if outfmt == "<STR_LIT>" : <EOL> return alignment . format ( "<STR_LIT>" ) <EOL> def muscle_inputorder ( inputfastafile , alnfile , trunc_name = True ) : <EOL> """<STR_LIT>""" <EOL> sh ( "<STR_LIT>" . format ( alnfile ) , log = False ) <EOL> maxi = <NUM_LIT:30> if trunc_name else <NUM_LIT:1000> <EOL> aa = AlignIO . read ( alnfile , "<STR_LIT>" ) <EOL> alignment = dict ( ( a . id [ : maxi ] , a ) for a in aa ) <EOL> if trunc_name and len ( alignment ) < len ( aa ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> fw = must_open ( alnfile , "<STR_LIT:w>" ) <EOL> for rec in SeqIO . parse ( inputfastafile , "<STR_LIT>" ) : <EOL> a = alignment [ rec . id [ : maxi ] ] <EOL> fw . write ( "<STR_LIT>" . format ( a . id [ : maxi ] , a . seq ) ) <EOL> fw . close ( ) <EOL> sh ( "<STR_LIT>" . format ( alnfile ) , log = False ) <EOL> def subset ( args ) : <EOL> """<STR_LIT>""" <EOL> p = OptionParser ( subset . __doc__ ) <EOL> p . add_option ( "<STR_LIT>" , action = "<STR_LIT:store_true>" , <EOL> help = "<STR_LIT>" ) <EOL> p . add_option ( "<STR_LIT>" , action = "<STR_LIT:store_true>" , <EOL> help = "<STR_LIT>" ) <EOL> p . set_stripnames ( ) <EOL> p . set_outfile ( ) <EOL> opts , args = p . parse_args ( args ) <EOL> if len ( args ) != <NUM_LIT:2> : <EOL> sys . exit ( not p . print_help ( ) ) <EOL> pairsfile , ksfiles = args [ <NUM_LIT:0> ] , args [ <NUM_LIT:1> : ] <EOL> noheader = opts . noheader <EOL> block = opts . block <EOL> if block : <EOL> noheader = True <EOL> outfile = opts . outfile <EOL> ksvals = { } <EOL> for ksfile in ksfiles : <EOL> ksvals . update ( dict ( ( line . name , line ) for line in KsFile ( ksfile , strip_names = opts . strip_names ) ) ) <EOL> fp = open ( pairsfile ) <EOL> fw = must_open ( outfile , "<STR_LIT:w>" ) <EOL> if not noheader : <EOL> print >> fw , fields <EOL> i = j = <NUM_LIT:0> <EOL> for row in fp : <EOL> if row [ <NUM_LIT:0> ] == '<STR_LIT:#>' : <EOL> if block : <EOL> print >> fw , row . strip ( ) <EOL> continue <EOL> a , b = row . split ( ) [ : <NUM_LIT:2> ] <EOL> name = "<STR_LIT:;>" . join ( ( a , b ) ) <EOL> if name not in ksvals : <EOL> name = "<STR_LIT:;>" . join ( ( b , a ) ) <EOL> if name not in ksvals : <EOL> j += <NUM_LIT:1> <EOL> print >> fw , "<STR_LIT:\t>" . join ( ( a , b , "<STR_LIT:.>" , "<STR_LIT:.>" ) ) <EOL> continue <EOL> ksline = ksvals [ name ] <EOL> if block : <EOL> print >> fw , "<STR_LIT:\t>" . join ( str ( x ) for x in ( a , b , ksline . ks ) ) <EOL> else : <EOL> ksline . name = "<STR_LIT:;>" . join ( ( a , b ) ) <EOL> print >> fw , ksline <EOL> i += <NUM_LIT:1> <EOL> fw . close ( ) <EOL> logging . debug ( "<STR_LIT>" . format ( j ) ) <EOL> logging . debug ( "<STR_LIT>" . format ( i , outfile ) ) <EOL> return outfile <EOL> fields = "<STR_LIT>" <EOL> descriptions = { <EOL> '<STR_LIT:name>' : '<STR_LIT>' , <EOL> '<STR_LIT>' : '<STR_LIT>' , <EOL> '<STR_LIT>' : '<STR_LIT>' , <EOL> '<STR_LIT>' : '<STR_LIT>' , <EOL> '<STR_LIT>' : '<STR_LIT>' } <EOL> class KsLine : <EOL> def __init__ ( self , row , strip_names = False ) : <EOL> args = row . strip ( ) . split ( "<STR_LIT:U+002C>" ) <EOL> self . name = args [ <NUM_LIT:0> ] <EOL> self . yn_ks = self . get_float ( args [ <NUM_LIT:1> ] ) <EOL> self . yn_ka = self . get_float ( args [ <NUM_LIT:2> ] ) <EOL> self . ng_ks = self . get_float ( args [ <NUM_LIT:3> ] ) <EOL> self . ng_ka = self . get_float ( args [ <NUM_LIT:4> ] ) <EOL> self . ks = self . ng_ks <EOL> if "<STR_LIT:;>" in self . name : <EOL> self . gene_a , self . gene_b = self . name . split ( "<STR_LIT:;>" ) <EOL> if strip_names : <EOL> self . gene_a = gene_name ( self . gene_a ) <EOL> self . gene_b = gene_name ( self . gene_b ) <EOL> def get_float ( self , x ) : <EOL> try : <EOL> x = float ( x ) <EOL> except : <EOL> x = - <NUM_LIT:1> <EOL> return x <EOL> def __str__ ( self ) : <EOL> return "<STR_LIT:U+002C>" . join ( str ( x ) for x in ( self . name , self . yn_ks , self . yn_ka , <EOL> self . ng_ks , self . ng_ka ) ) <EOL> @ property <EOL> def anchorline ( self ) : <EOL> return "<STR_LIT:\t>" . join ( ( gene_name ( self . gene_a ) , gene_name ( self . gene_b ) , <EOL> "<STR_LIT>" . format ( self . ks ) ) ) <EOL> class KsFile ( LineFile ) : <EOL> def __init__ ( self , filename , strip_names = False ) : <EOL> super ( KsFile , self ) . __init__ ( filename ) <EOL> fp = open ( filename ) <EOL> for row in fp : <EOL> ksline = KsLine ( row , strip_names = strip_names ) <EOL> if ksline . name == "<STR_LIT:name>" : <EOL> continue <EOL> self . append ( ksline ) <EOL> logging . debug ( '<STR_LIT>' . format ( filename , len ( self ) ) ) <EOL> def print_to_anchors ( self , outfile ) : <EOL> fw = must_open ( outfile , "<STR_LIT:w>" ) <EOL> for row in self : <EOL> print >> fw , row . anchorline <EOL> fw . close ( ) <EOL> def my_hist ( ax , l , interval , max_r , color = '<STR_LIT:g>' , marker = '<STR_LIT:.>' , fill = False ) : <EOL> if not l : <EOL> return <EOL> n , p = [ ] , [ ] <EOL> total_len = len ( l ) <EOL> for i in np . arange ( <NUM_LIT:0> , max_r , interval ) : <EOL> xmin , xmax = i - <NUM_LIT> * interval , i + <NUM_LIT> * interval <EOL> nx = [ x for x in l if xmin <= x < xmax ] <EOL> n . append ( i ) <EOL> p . append ( len ( nx ) * <NUM_LIT> / total_len ) <EOL> if fill : <EOL> from pylab import poly_between <EOL> xs , ys = poly_between ( n , <NUM_LIT:0> , p ) <EOL> line = ax . fill ( xs , ys , fc = color , alpha = <NUM_LIT> ) <EOL> else : <EOL> line = ax . plot ( n , p , color = color , lw = <NUM_LIT:2> , ms = <NUM_LIT:3> , <EOL> marker = marker , mfc = "<STR_LIT:w>" , mec = color , mew = <NUM_LIT:2> ) <EOL> return line <EOL> def lognormpdf ( bins , mu , sigma ) : <EOL> return np . exp ( - ( np . log ( bins ) - mu ) ** <NUM_LIT:2> / ( <NUM_LIT:2> * sigma ** <NUM_LIT:2> ) ) / ( bins * sigma * sqrt ( <NUM_LIT:2> * pi ) ) <EOL> def lognormpdf_mix ( bins , probs , mus , sigmas , interval = <NUM_LIT> ) : <EOL> y = <NUM_LIT:0> <EOL> for prob , mu , sigma in zip ( probs , mus , sigmas ) : <EOL> y += prob * lognormpdf ( bins , mu , sigma ) <EOL> y *= <NUM_LIT:100> * interval <EOL> return y <EOL> def get_mixture ( data , components ) : <EOL> """<STR_LIT>""" <EOL> from jcvi . apps . base import popen <EOL> probs , mus , sigmas = [ ] , [ ] , [ ] <EOL> fw = must_open ( "<STR_LIT>" , "<STR_LIT:w>" ) <EOL> log_data = [ log ( x ) for x in data if x > <NUM_LIT> ] <EOL> data = "<STR_LIT:\n>" . join ( [ "<STR_LIT>" % x for x in log_data ] ) . replace ( "<STR_LIT>" , "<STR_LIT>" ) <EOL> fw . write ( data ) <EOL> fw . close ( ) <EOL> cmd = "<STR_LIT>" . format ( components , len ( log_data ) , fw . name ) <EOL> pipe = popen ( cmd ) <EOL> for row in pipe : <EOL> if row [ <NUM_LIT:0> ] != '<STR_LIT:#>' : <EOL> continue <EOL> atoms = row . split ( "<STR_LIT:U+002C>" ) <EOL> a , b , c = atoms [ <NUM_LIT:1> : <NUM_LIT:4> ] <EOL> a = float ( a ) <EOL> b = float ( b ) <EOL> c = float ( c ) <EOL> mus . append ( a ) <EOL> sigmas . append ( b ) <EOL> probs . append ( c ) <EOL> os . remove ( fw . name ) <EOL> return probs , mus , sigmas <EOL> def plot_ks_dist ( ax , data , interval , components , ks_max , <EOL> color = '<STR_LIT:r>' , marker = '<STR_LIT:.>' , fill = False , fitted = True ) : <EOL> line , = my_hist ( ax , data , interval , ks_max , <EOL> color = color , marker = marker , fill = fill ) <EOL> logging . debug ( "<STR_LIT>" . format ( len ( data ) ) ) <EOL> line_mixture = None <EOL> if fitted : <EOL> probs , mus , variances = get_mixture ( data , components ) <EOL> iv = <NUM_LIT> <EOL> bins = np . arange ( iv , ks_max , iv ) <EOL> y = lognormpdf_mix ( bins , probs , mus , variances , interval ) <EOL> line_mixture , = ax . plot ( bins , y , '<STR_LIT::>' , color = color , lw = <NUM_LIT:3> ) <EOL> for i in xrange ( components ) : <EOL> peak_val = exp ( mus [ i ] ) <EOL> mixline = lognormpdf_mix ( peak_val , probs , mus , variances , interval ) <EOL> ax . text ( peak_val , mixline , "<STR_LIT>" % peak_val , color = "<STR_LIT:w>" , size = <NUM_LIT:10> , bbox = dict ( ec = '<STR_LIT:w>' , fc = color , alpha = <NUM_LIT> , boxstyle = '<STR_LIT>' ) ) <EOL> return line , line_mixture <EOL> def add_plot_options ( p ) : <EOL> p . add_option ( "<STR_LIT>" , default = False , action = "<STR_LIT:store_true>" , <EOL> help = "<STR_LIT>" ) <EOL> p . add_option ( "<STR_LIT>" , default = <NUM_LIT:0.> , type = "<STR_LIT:float>" , <EOL> help = "<STR_LIT>" ) <EOL> p . add_option ( "<STR_LIT>" , default = <NUM_LIT> , type = "<STR_LIT:float>" , <EOL> help = "<STR_LIT>" ) <EOL> p . add_option ( "<STR_LIT>" , default = <NUM_LIT> , type = "<STR_LIT:int>" , <EOL> help = "<STR_LIT>" ) <EOL> p . add_option ( "<STR_LIT>" , default = "<STR_LIT>" , <EOL> help = "<STR_LIT>" ) <EOL> p . add_option ( "<STR_LIT>" , dest = "<STR_LIT>" , default = True , action = "<STR_LIT>" , <EOL> help = "<STR_LIT>" ) <EOL> p . add_option ( "<STR_LIT>" , default = "<STR_LIT>" , <EOL> help = "<STR_LIT>" ) <EOL> def report ( args ) : <EOL> '''<STR_LIT>''' <EOL> from jcvi . utils . cbook import SummaryStats <EOL> from jcvi . graphics . histogram import stem_leaf_plot <EOL> p = OptionParser ( report . __doc__ ) <EOL> p . add_option ( "<STR_LIT>" , default = False , action = "<STR_LIT:store_true>" , <EOL> help = "<STR_LIT>" ) <EOL> p . add_option ( "<STR_LIT>" , default = <NUM_LIT:1> , type = "<STR_LIT:int>" , <EOL> help = "<STR_LIT>" ) <EOL> add_plot_options ( p ) <EOL> opts , args , iopts = p . set_image_options ( args , figsize = "<STR_LIT>" ) <EOL> if len ( args ) != <NUM_LIT:1> : <EOL> sys . exit ( not p . print_help ( ) ) <EOL> ks_file , = args <EOL> data = KsFile ( ks_file ) <EOL> ks_min = opts . vmin <EOL> ks_max = opts . vmax <EOL> bins = opts . bins <EOL> for f in fields . split ( "<STR_LIT:U+002C>" ) [ <NUM_LIT:1> : ] : <EOL> columndata = [ getattr ( x , f ) for x in data ] <EOL> ks = ( "<STR_LIT>" in f ) <EOL> if not ks : <EOL> continue <EOL> columndata = [ x for x in columndata if ks_min <= x <= ks_max ] <EOL> st = SummaryStats ( columndata ) <EOL> title = "<STR_LIT>" . format ( descriptions [ f ] , ks_file ) <EOL> title += "<STR_LIT>" . format ( st . median , st . firstq , st . thirdq ) <EOL> title += "<STR_LIT>" . format ( st . mean , st . sd , st . size ) <EOL> tbins = ( <NUM_LIT:0> , ks_max , bins ) if ks else ( <NUM_LIT:0> , <NUM_LIT> , <NUM_LIT:10> ) <EOL> digit = <NUM_LIT:2> if ( ks_max * <NUM_LIT:1.> / bins ) < <NUM_LIT> else <NUM_LIT:1> <EOL> stem_leaf_plot ( columndata , * tbins , digit = digit , title = title ) <EOL> if not opts . pdf : <EOL> return <EOL> components = opts . components <EOL> data = [ x . ng_ks for x in data ] <EOL> data = [ x for x in data if ks_min <= x <= ks_max ] <EOL> fig = plt . figure ( <NUM_LIT:1> , ( iopts . w , iopts . h ) ) <EOL> ax = fig . add_axes ( [ <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ] ) <EOL> kp = KsPlot ( ax , ks_max , opts . bins , legendp = opts . legendp ) <EOL> kp . add_data ( data , components , fill = opts . fill , fitted = opts . fit ) <EOL> kp . draw ( title = opts . title ) <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> main ( ) </s>
