signature:
def loci2bpp ( name , locifile , imap , guidetree , minmap = None , maxloci = None , infer_sptree = 0 , infer_delimit = 0 , delimit_alg = ( 0 , 5 ) , seed = 12345 , burnin = 1000 , nsample = 10000 , sampfreq = 2 , thetaprior = ( 5 , 5 ) , tauprior = ( 4 , 2 , 1 ) , traits_df = None , nu = 0 , kappa = 0 , useseqdata = 1 , usetraitdata = 1 , cleandata = 0 , wdir = None , finetune = ( 0.01 , 0.01 , 0.01 , 0.01 , 0.01 , 0.01 , 0.01 , 0.01 ) , verbose = 0 ) :
    """Converts loci file format to bpp file format , i . e . , concatenated phylip - like
 format , and produces imap and ctl input files for bpp .
 Parameters :
 name :
 A prefix name for output files that will be produced
 locifile :
 A . loci file produced by ipyrad .
 imap :
 A Python dictionary with ' species ' names as keys , and lists of sample
 names for the values . Any sample that is not included in the imap
 dictionary will be filtered out of the data when converting the . loci
 file into the bpp formatted sequence file . Each species in the imap
 dictionary must also be present in the input ' guidetree ' .
 guidetree :
 A newick string species tree hypothesis [ e . g . , ( ( ( a , b ) , ( c , d ) ) , e ) ; ]
 All species in the imap dictionary must also be present in the guidetree
 Optional parameters :
 infer _ sptree :
 Default = 0 , only infer parameters on a fixed species tree . If 1 , then the
 input tree is treated as a guidetree and tree search is employed to find
 the best tree . The results will include support values for the inferred
 topology .
 infer _ delimit :
 Default = 0 , no delimitation . If 1 then splits in the tree that separate
 ' species ' will be collapsed to test whether fewer species are a better
 fit to the data than the number in the input guidetree .
 delimit _ alg :
 Species delimitation algorithm . This is a tuple . The first value
 is the algorithm ( 0 or 1 ) and the following values are arguments
 for the given algorithm . See other ctl files for examples of what the
 delimitation line looks like . This is where you can enter the params
 ( e . g . , alpha , migration ) for the two different algorithms .
 For example , the following args would produce the following ctl lines :
 alg = 0 , epsilon = 5
 > delimit _ alg = ( 0 , 5)
 speciesdelimitation = 1 0 5
 alg = 1 , alpha = 2 , migration = 1
 > delimit _ alg = ( 1 , 2 , 1)
 speciesdelimitation = 1 1 2 1
 alg = 1 , alpha = 2 , migration = 1 , diagnosis = 0 , ? = 1
 > delimit _ alg = ( 1 , 2 , 1 , 0 , 1)
 speciesdelimitation = 1 1 2 1 0 1
 seed :
 A random number seed at start of analysis .
 burnin :
 Number of burnin generations in mcmc
 nsample :
 Number of mcmc generations to run .
 sampfreq :
 How often to sample from the mcmc chain .
 thetaprior :
 Prior on theta ( 4Neu ) , gamma distributed . mean = a / b . e . g . , ( 5 , 5)
 tauprior
 Prior on root tau , gamma distributed mean = a / b . Last number is
 dirichlet prior for other taus . e . g . , ( 4 , 2 , 1)
 traits _ df :
 A pandas DataFrame with trait data properly formatted . This means only
 quantitative traits are included , and missing values are NaN .
 The first column contains sample names , with " Indiv " as the header .
 The following columns have a header row with trait names . This script
 will write a CSV trait file with trait values mean - standardized , with
 NaN replaced by " NA " , and with sample not present in IMAP removed .
 nu :
 A prior on phenotypic trait variance ( 0 ) for iBPP analysis .
 kappa :
 A prior on phenotypic trait mean ( 0 ) for iBPP analysis .
 useseqdata :
 If false inference proceeds without sequence data ( can be used to test
 the effect of priors on the tree distributions ) .
 usetraitdata :
 If false inference proceeds without trait data ( can be used to test
 the effect of priors on the trait distributions ) .
 cleandata :
 If 1 then sites with missing or hetero characters are removed .
 wdir :
 A working directory to write files to .
 finetune :
 See bpp documentation .
 verbose :
 If verbose = 1 the ctl file text will also be written to screen ( stderr ) ."""

implementation:
# # check args
if not imap :
    raise IPyradWarningExit ( IMAP_REQUIRED )
if minmap :
    if minmap . keys ( ) != imap . keys ( ) :
        raise IPyradWarningExit ( KEYS_DIFFER )
# # working directory , make sure it exists
if wdir :
    wdir = os . path . abspath ( wdir )
    if not os . path . exists ( wdir ) :
        raise IPyradWarningExit ( " working directory (wdir) does not exist" )
else :
    wdir = os . path . curdir
# # if traits _ df then we make ' . ibpp ' files
prog = 'bpp'
if isinstance ( traits_df , pd . DataFrame ) :
    prog = 'ibpp'
outfile = OPJ ( wdir , "{}.{}.seq.txt" . format ( name , prog ) )
mapfile = OPJ ( wdir , "{}.{}.imap.txt" . format ( name , prog ) )
# # open outhandles
fout = open ( outfile , 'w' )
fmap = open ( mapfile , 'w' )
# # parse the loci file
with open ( locifile , 'r' ) as infile : # # split on " / / " for legacy compatibility
    loci = infile . read ( ) . strip ( ) . split ( "|\n" )
    nloci = len ( loci )
# # all samples
samples = list ( itertools . chain ( * imap . values ( ) ) )
# # iterate over loci , printing to outfile
nkept = 0
for iloc in xrange ( nloci ) :
    lines = loci [ iloc ] . split ( "//" ) [ 0 ] . split ( )
    names = lines [ : : 2 ]
    names = [ "^" + i for i in names ]
    seqs = [ list ( i ) for i in lines [ 1 : : 2 ] ]
    seqlen = len ( seqs [ 0 ] )
    # # whether to skip this locus based on filters below
    skip = 0
    # # if minmap filter for sample coverage
    if minmap :
        covd = { }
        for group , vals in imap . items ( ) :
            covd [ group ] = sum ( [ "^" + i in names for i in vals ] )
        # # check that coverage is good enough
        if not all ( [ covd [ group ] >= minmap [ group ] for group in minmap ] ) :
            skip = 1
    # # too many loci ?
    if maxloci :
        if nkept >= maxloci :
            skip = 1
    # # build locus as a string
    if not skip : # # convert to phylip with caret starter and replace - with N .
        data = [ "{:<30} {}" . format ( i , "" . join ( k ) . replace ( "-" , "N" ) ) for ( i , k ) in zip ( names , seqs ) if i [ 1 : ] in samples ]
        # # if not empty , write to the file
        if data :
            fout . write ( "{} {}\n\n{}\n\n" . format ( len ( data ) , seqlen , "\n" . join ( data ) ) )
            nkept += 1
# # close up shop
fout . close ( )
# # write the imap file :
data = [ "{:<30} {}" . format ( val , key ) for key in sorted ( imap ) for val in imap [ key ] ]
fmap . write ( "\n" . join ( data ) )
fmap . close ( )
# # write ctl file
write_ctl ( name , imap , guidetree , nkept , infer_sptree , infer_delimit , delimit_alg , seed , burnin , nsample , sampfreq , thetaprior , tauprior , traits_df , nu , kappa , cleandata , useseqdata , usetraitdata , wdir , finetune , verbose )
# # print message ?
sys . stderr . write ( "new files created ({} loci, {} species, {} samples)\n" . format ( nkept , len ( imap . keys ( ) ) , sum ( [ len ( i ) for i in imap . values ( ) ] ) ) )
sys . stderr . write ( "  {}.{}.seq.txt\n" . format ( name , prog ) )
sys . stderr . write ( "  {}.{}.imap.txt\n" . format ( name , prog ) )
sys . stderr . write ( "  {}.{}.ctl.txt\n" . format ( name , prog ) )
if isinstance ( traits_df , pd . DataFrame ) :
    sys . stderr . write ( "  {}.{}.traits.txt\n" . format ( name , prog ) )
# # return the ctl file string
return os . path . abspath ( "{}.{}.ctl.txt" . format ( OPJ ( wdir , name ) , prog ) )


signature:
def write_ctl ( name , imap , guidetree , nloci , infer_sptree , infer_delimit , delimit_alg , seed , burnin , nsample , sampfreq , thetaprior , tauprior , traits_df , nu0 , kappa0 , cleandata , useseqdata , usetraitdata , wdir , finetune , verbose ) :
    """write outfile with any args in argdict"""

implementation:
# # A string to store ctl info
ctl = [ ]
# # check the tree ( can do this better once we install ete3 w / ipyrad )
if not guidetree . endswith ( ";" ) :
    guidetree += ";"
# # if traits _ df then we make ' . ibpp ' files
prog = 'bpp'
if isinstance ( traits_df , pd . DataFrame ) :
    prog = 'ibpp'
# # write the top header info
ctl . append ( "seed = {}" . format ( seed ) )
ctl . append ( "seqfile = {}.{}.seq.txt" . format ( OPJ ( wdir , name ) , prog ) )
ctl . append ( "Imapfile = {}.{}.imap.txt" . format ( OPJ ( wdir , name ) , prog ) )
ctl . append ( "mcmcfile = {}.{}.mcmc.txt" . format ( OPJ ( wdir , name ) , prog ) )
ctl . append ( "outfile = {}.{}.out.txt" . format ( OPJ ( wdir , name ) , prog ) )
if isinstance ( traits_df , pd . DataFrame ) :
    ctl . append ( "traitfile = {}.{}.traits.txt" . format ( OPJ ( wdir , name ) , prog ) )
# # number of loci ( checks that seq file exists and parses from there )
ctl . append ( "nloci = {}" . format ( nloci ) )
ctl . append ( "usedata = {}" . format ( useseqdata ) )
ctl . append ( "cleandata = {}" . format ( cleandata ) )
# # infer species tree
if infer_sptree :
    ctl . append ( "speciestree = 1 0.4 0.2 0.1" )
else :
    ctl . append ( "speciestree = 0" )
# # infer delimitation ( with algorithm 1 by default )
ctl . append ( "speciesdelimitation = {} {} {}" . format ( infer_delimit , delimit_alg [ 0 ] , " " . join ( [ str ( i ) for i in delimit_alg [ 1 : ] ] ) ) )
# # if using iBPP ( if not traits _ df , we assume you ' re using bpp ( v . 3.3 + )
if isinstance ( traits_df , pd . DataFrame ) : # # check that the data frame is properly formatted
    try :
        traits_df . values . astype ( float )
    except Exception :
        raise IPyradWarningExit ( PDREAD_ERROR )
    # # subsample to keep only samples that are in IMAP , we do not need to
    # # standarize traits b / c ibpp does that for us .
    samples = sorted ( list ( itertools . chain ( * imap . values ( ) ) ) )
    didx = [ list ( traits_df . index ) . index ( i ) for i in traits_df . index if i not in samples ]
    dtraits = traits_df . drop ( traits_df . index [ didx ] )
    # # mean standardize traits values after excluding samples
    straits = dtraits . apply ( lambda x : ( x - x . mean ( ) ) / ( x . std ( ) ) )
    # # convert NaN to " NA " cuz that ' s what ibpp likes , and write to file
    ftraits = straits . fillna ( "NA" )
    traitdict = ftraits . T . to_dict ( "list" )
    # # get reverse imap dict
    rev = { val : key for key in sorted ( imap ) for val in imap [ key ] }
    # # write trait file
    traitfile = "{}.{}.traits.txt" . format ( os . path . join ( wdir , name ) , prog )
    with open ( traitfile , 'w' ) as tout :
        tout . write ( "Indiv\n" )
        tout . write ( "\t" . join ( [ 'Species' ] + list ( ftraits . columns ) ) + "\n" )
        # for key in sorted ( traitdict ) :
        # tout . write ( " \ t " . join ( [ key , rev [ key ] ] + \
        # [ " ^ " + str ( i ) for i in traitdict [ key ] ] ) + " \ n "
        nindT = 0
        for ikey in sorted ( imap . keys ( ) ) :
            samps = imap [ ikey ]
            for samp in sorted ( samps ) :
                if samp in traitdict :
                    tout . write ( "\t" . join ( [ samp , rev [ samp ] ] + [ str ( i ) for i in traitdict [ samp ] ] ) + "\n" )
                    nindT += 1
    # tout . write ( " Indiv \ n " + " \ t " . join ( [ " Species " ] + \
    # [ " t _ { } " . format ( i ) for i in range ( len ( traitdict . values ( ) [ 0 ] ) ) ] ) + " \ n " )
    # for key in sorted ( traitdict ) :
    # print > > tout , " \ t " . join ( [ key , rev [ key ] ] + \
    # [ str ( i ) for i in traitdict [ key ] ] )
    # ftraits . to _ csv ( traitfile )
    # # write ntraits and nindT and traitfilename
    ctl . append ( "ntraits = {}" . format ( traits_df . shape [ 1 ] ) )
    ctl . append ( "nindT = {}" . format ( nindT ) )
    # traits _ df . shape [ 0 ] ) )
    ctl . append ( "usetraitdata = {}" . format ( usetraitdata ) )
    ctl . append ( "useseqdata = {}" . format ( useseqdata ) )
    # # trait priors
    ctl . append ( "nu0 = {}" . format ( nu0 ) )
    ctl . append ( "kappa0 = {}" . format ( kappa0 ) )
    # # remove ibpp incompatible options
    ctl . remove ( "usedata = {}" . format ( useseqdata ) )
    ctl . remove ( "speciestree = {}" . format ( infer_sptree ) )
# # get tree values
nspecies = str ( len ( imap ) )
species = " " . join ( sorted ( imap ) )
ninds = " " . join ( [ str ( len ( imap [ i ] ) ) for i in sorted ( imap ) ] )
# # write the tree
ctl . append ( """\
species&tree = {} {}
                 {}
                 {}""" . format ( nspecies , species , ninds , guidetree ) )
# # priors
ctl . append ( "thetaprior = {} {}" . format ( * thetaprior ) )
ctl . append ( "tauprior = {} {} {}" . format ( * tauprior ) )
# # other values , fixed for now
ctl . append ( "finetune = 1: {}" . format ( " " . join ( [ str ( i ) for i in finetune ] ) ) )
# CTL . append ( " finetune = 1 : 1 0.002 0.01 0.01 0.02 0.005 1.0 " )
ctl . append ( "print = 1 0 0 0" )
ctl . append ( "burnin = {}" . format ( burnin ) )
ctl . append ( "sampfreq = {}" . format ( sampfreq ) )
ctl . append ( "nsample = {}" . format ( nsample ) )
# # write out the ctl file
with open ( "{}.{}.ctl.txt" . format ( OPJ ( wdir , name ) , prog ) , 'w' ) as out :
    out . write ( "\n" . join ( ctl ) )
# # if verbose print ctl
if verbose :
    sys . stderr . write ( "ctl file\n--------\n" + "\n" . join ( ctl ) + "\n--------\n\n" )


signature:
def _collapse_outgroup ( tree , taxdicts ) :
    """collapse outgroup in ete Tree for easier viewing"""

implementation:
# # check that all tests have the same outgroup
outg = taxdicts [ 0 ] [ "p4" ]
if not all ( [ i [ "p4" ] == outg for i in taxdicts ] ) :
    raise Exception ( "no good" )
# # prune tree , keep only one sample from outgroup
tre = ete . Tree ( tree . write ( format = 1 ) )
# tree . copy ( method = " deepcopy " )
alltax = [ i for i in tre . get_leaf_names ( ) if i not in outg ]
alltax += [ outg [ 0 ] ]
tre . prune ( alltax )
tre . search_nodes ( name = outg [ 0 ] ) [ 0 ] . name = "outgroup"
tre . ladderize ( )
# # remove other ougroups from taxdicts
taxd = copy . deepcopy ( taxdicts )
newtaxdicts = [ ]
for test in taxd : # test [ " p4 " ] = [ outg [ 0 ] ]
    test [ "p4" ] = [ "outgroup" ]
    newtaxdicts . append ( test )
return tre , newtaxdicts


signature:
def _decompose_tree ( ttree , orient = 'right' , use_edge_lengths = True ) :
    """decomposes tree into component parts for plotting"""

implementation:
# # set attributes
ttree . _orient = orient
ttree . _use_edge_lengths = use_edge_lengths
ult = use_edge_lengths == False
# # map numeric values to internal nodes from root to tips
names = { }
idx = 0
for node in ttree . tree . traverse ( "preorder" ) :
    if not node . is_leaf ( ) :
        if node . name :
            names [ idx ] = node . name
        else :
            names [ idx ] = idx
            node . name = str ( idx )
        node . idx = idx
        idx += 1
# # map number to the tips , these will be the highest numbers
for node in ttree . tree . get_leaves ( ) :
    names [ idx ] = node . name
    node . idx = idx
    idx += 1
# # create empty edges and coords arrays
ttree . node_labels = names
ttree . tip_labels = ttree . tree . get_leaf_names ( )
# self . tip _ labels = self . tree . get _ leaf _ names ( ) [ : : - 1]
# self . node _ labels = self . names
ttree . edges = np . zeros ( ( idx - 1 , 2 ) , dtype = int )
ttree . verts = np . zeros ( ( idx , 2 ) , dtype = float )
ttree . _lines = [ ]
# np . zeros ( ( ntips - 1 ) , dtype = int )
ttree . _coords = [ ]
# np . zeros ( ( idx * 2 - ntips ) , dtype = float )
# # postorder : first children and then parents . This moves up the list .
nidx = 0
tip_num = len ( ttree . tree . get_leaves ( ) ) - 1
# # tips to root to fill in the verts and edges
for node in ttree . tree . traverse ( "postorder" ) :
    if node . is_leaf ( ) : # # set the xy - axis positions of the tips
        node . y = ttree . tree . get_distance ( node )
        if ult :
            node . y = 0.
        node . x = tip_num
        tip_num -= 1
        # # edges connect this vert to
        ttree . verts [ node . idx ] = [ node . x , node . y ]
        ttree . edges [ nidx ] = [ node . up . idx , node . idx ]
    elif node . is_root ( ) :
        node . y = ttree . tree . get_distance ( node )
        if ult :
            node . y = - 1 * node . get_farthest_leaf ( True ) [ 1 ] - 1
        node . x = sum ( i . x for i in node . children ) / float ( len ( node . children ) )
        ttree . verts [ node . idx ] = [ node . x , node . y ]
    else : # # create new nodes left and right
        node . y = ttree . tree . get_distance ( node )
        if ult :
            node . y = - 1 * node . get_farthest_leaf ( True ) [ 1 ] - 1
        node . x = sum ( i . x for i in node . children ) / float ( len ( node . children ) )
        ttree . edges [ nidx , : ] = [ node . up . idx , node . idx ]
        ttree . verts [ node . idx ] = [ node . x , node . y ]
    nidx += 1
# # root to tips to fill in the coords and lines
cidx = 0
for node in ttree . tree . traverse ( ) : # # add yourself
    if not node . is_leaf ( ) :
        ttree . _coords += [ [ node . x , node . y ] ]
        pidx = cidx
        cidx += 1
        for child in node . children : # # add children
            ttree . _coords += [ [ child . x , node . y ] , [ child . x , child . y ] ]
            ttree . _lines += [ [ pidx , cidx ] ]
            # # connect yourself to newx
            ttree . _lines += [ [ cidx , cidx + 1 ] ]
            # # connect newx to child
            cidx += 2
ttree . _coords = np . array ( ttree . _coords , dtype = float )
ttree . _lines = np . array ( ttree . _lines , dtype = int )
# # invert for sideways trees
if ttree . _orient in [ 'up' , 0 ] :
    pass
if ttree . _orient in [ 'left' , 1 ] :
    ttree . verts [ : , 1 ] = ttree . verts [ : , 1 ] * - 1
    ttree . verts = ttree . verts [ : , [ 1 , 0 ] ]
    ttree . _coords [ : , 1 ] = ttree . _coords [ : , 1 ] * - 1
    ttree . _coords = ttree . _coords [ : , [ 1 , 0 ] ]
if ttree . _orient in [ 'down' , 0 ] :
    ttree . verts [ : , 1 ] = ttree . verts [ : , 1 ] * - 1
    ttree . _coords [ : , 1 ] = ttree . _coords [ : , 1 ] * - 1
if ttree . _orient in [ 'right' , 3 ] :
    ttree . verts = ttree . verts [ : , [ 1 , 0 ] ]
    ttree . _coords = ttree . _coords [ : , [ 1 , 0 ] ]


signature:
def draw ( self , show_tip_labels = True , show_node_support = False , use_edge_lengths = False , orient = "right" , print_args = False , * args , ** kwargs ) :
    """plot the tree using toyplot . graph .
 Parameters :
 show _ tip _ labels : bool
 Show tip names from tree .
 use _ edge _ lengths : bool
 Use edge lengths from newick tree .
 show _ node _ support : bool
 Show support values at nodes using a set of default
 options ."""

implementation:
# # re - decompose tree for new orient and edges args
self . _decompose_tree ( orient = orient , use_edge_lengths = use_edge_lengths )
# # update kwargs with entered args and all other kwargs
dwargs = { }
dwargs [ "show_tip_labels" ] = show_tip_labels
dwargs [ "show_node_support" ] = show_node_support
dwargs . update ( kwargs )
# # pass to panel plotter
canvas , axes , panel = tree_panel_plot ( self , print_args , ** dwargs )
return canvas , axes , panel


signature:
def tree_panel_plot ( ttree , print_args = False , * args , ** kwargs ) :
    """signature . . ."""

implementation:
# # create Panel plot object and set height & width
panel = Panel ( ttree )
# tree , edges , verts , names )
if not kwargs . get ( "width" ) :
    panel . kwargs [ "width" ] = min ( 1000 , 25 * len ( panel . tree ) )
if not kwargs . get ( "height" ) :
    panel . kwargs [ "height" ] = panel . kwargs [ "width" ]
# # update defaults with kwargs & update size based on ntips & ntests
panel . kwargs . update ( kwargs )
# # magic node label arguments overrides others
if panel . kwargs [ "show_node_support" ] :
    nnodes = sum ( 1 for i in panel . tree . traverse ( ) ) - len ( panel . tree )
    # # set node values
    supps = [ int ( panel . tree . search_nodes ( idx = j ) [ 0 ] . support ) for j in range ( nnodes ) ]
    if not panel . kwargs [ "vsize" ] :
        panel . kwargs [ "vsize" ] = 20
    sizes = [ panel . kwargs [ "vsize" ] for j in range ( nnodes ) ]
    # # add leaf values
    supps += [ "" ] * len ( panel . tree )
    sizes += [ 0 ] * len ( panel . tree )
    # # override args
    panel . kwargs [ "vlabel" ] = supps
    panel . kwargs [ "vsize" ] = sizes
    panel . kwargs [ "vlshow" ] = True
    # panel . kwargs [ " vmarker " ] = ' s ' # # square
    # # if unrooted then hide root node scores
    if len ( panel . tree . children ) > 2 :
        supps [ 0 ] = ""
        sizes [ 0 ] = 0
    # print ( panel . kwargs [ " vlabels " ] )
    # print ( panel . kwargs [ " vsize " ] )
elif panel . kwargs . get ( "vlabel" ) :
    panel . kwargs [ "vlabel" ] = panel . kwargs [ "vlabel" ]
    panel . kwargs [ "vlshow" ] = True
else :
    panel . kwargs [ "vlabel" ] = panel . node_labels . keys ( )
# names . keys ( )
# # debugger / see all options
if print_args :
    print ( panel . kwargs )
# # maybe add panels for plotting tip traits in the future
# # create a canvas and a single cartesian coord system
canvas = toyplot . Canvas ( height = panel . kwargs [ 'height' ] , width = panel . kwargs [ 'width' ] )
axes = canvas . cartesian ( bounds = ( "10%" , "90%" , "10%" , "90%" ) )
axes . show = panel . kwargs [ "show_axes" ]
# # add panel plots to the axis
panel . _panel_tree ( axes )
if panel . kwargs [ "show_tip_labels" ] :
    panel . _panel_tip_labels ( axes )
return canvas , axes , panel


signature:
def get_quick_depths ( data , sample ) :
    """iterate over clustS files to get data"""

implementation:
# # use existing sample cluster path if it exists , since this
# # func can be used in step 4 and that can occur after merging
# # assemblies after step3 , and if we then referenced by data . dirs . clusts
# # the path would be broken .
# # If branching at step 3 to test different clust thresholds , the
# # branched samples will retain the samples . files . clusters of the
# # parent ( which have the clust _ threshold value of the parent ) , so
# # it will look like nothing has changed . If we call this func
# # from step 3 then it indicates we are in a branch and should
# # reset the sample . files . clusters handle to point to the correct
# # data . dirs . clusts directory . See issue # 229.
# # Easier to just always trust that samples . files . clusters is right ,
# # no matter what step ?
# if sample . files . clusters and not sample . stats . state = = 3:
# pass
# else :
# # # set cluster file handles
sample . files . clusters = os . path . join ( data . dirs . clusts , sample . name + ".clustS.gz" )
# # get new clustered loci
fclust = data . samples [ sample . name ] . files . clusters
clusters = gzip . open ( fclust , 'r' )
pairdealer = itertools . izip ( * [ iter ( clusters ) ] * 2 )
# # storage
depths = [ ]
maxlen = [ ]
# # start with cluster 0
tdepth = 0
tlen = 0
# # iterate until empty
while 1 : # # grab next
    try :
        name , seq = pairdealer . next ( )
    except StopIteration :
        break
    # # if not the end of a cluster
    # print name . strip ( ) , seq . strip ( )
    if name . strip ( ) == seq . strip ( ) :
        depths . append ( tdepth )
        maxlen . append ( tlen )
        tlen = 0
        tdepth = 0
    else :
        tdepth += int ( name . split ( ";" ) [ - 2 ] [ 5 : ] )
        tlen = len ( seq )
# # return
clusters . close ( )
return np . array ( maxlen ) , np . array ( depths )


signature:
def sample_cleanup ( data , sample ) :
    """stats , cleanup , and link to samples"""

implementation:
# # get maxlen and depths array from clusters
maxlens , depths = get_quick_depths ( data , sample )
try :
    depths . max ( )
except ValueError : # # If depths is an empty array max ( ) will raise
    print ( "    no clusters found for {}" . format ( sample . name ) )
    return
# # Test if depths is non - empty , but just full of zeros .
if depths . max ( ) : # # store which min was used to calculate hidepth here
    sample . stats_dfs . s3 [ "hidepth_min" ] = data . paramsdict [ "mindepth_majrule" ]
    # # If our longest sequence is longer than the current max _ fragment _ length
    # # then update max _ fragment _ length . For assurance we require that
    # # max len is 4 greater than maxlen , to allow for pair separators .
    hidepths = depths >= data . paramsdict [ "mindepth_majrule" ]
    maxlens = maxlens [ hidepths ]
    # # Handle the case where there are no hidepth clusters
    if maxlens . any ( ) :
        maxlen = int ( maxlens . mean ( ) + ( 2. * maxlens . std ( ) ) )
    else :
        maxlen = 0
    if maxlen > data . _hackersonly [ "max_fragment_length" ] :
        data . _hackersonly [ "max_fragment_length" ] = maxlen + 4
    # # make sense of stats
    keepmj = depths [ depths >= data . paramsdict [ "mindepth_majrule" ] ]
    keepstat = depths [ depths >= data . paramsdict [ "mindepth_statistical" ] ]
    # # sample summary stat assignments
    sample . stats [ "state" ] = 3
    sample . stats [ "clusters_total" ] = depths . shape [ 0 ]
    sample . stats [ "clusters_hidepth" ] = keepmj . shape [ 0 ]
    # # store depths histogram as a dict . Limit to first 25 bins
    bars , bins = np . histogram ( depths , bins = range ( 1 , 26 ) )
    sample . depths = { int ( i ) : v for i , v in zip ( bins , bars ) if v }
    # # sample stat assignments
    # # Trap numpy warnings ( " mean of empty slice " ) printed by samples
    # # with few reads .
    with warnings . catch_warnings ( ) :
        warnings . simplefilter ( "ignore" , category = RuntimeWarning )
        sample . stats_dfs . s3 [ "merged_pairs" ] = sample . stats . reads_merged
        sample . stats_dfs . s3 [ "clusters_total" ] = depths . shape [ 0 ]
        try :
            sample . stats_dfs . s3 [ "clusters_hidepth" ] = int ( sample . stats [ "clusters_hidepth" ] )
        except ValueError : # # Handle clusters _ hidepth = = NaN
            sample . stats_dfs . s3 [ "clusters_hidepth" ] = 0
        sample . stats_dfs . s3 [ "avg_depth_total" ] = depths . mean ( )
        LOGGER . debug ( "total depth {}" . format ( sample . stats_dfs . s3 [ "avg_depth_total" ] ) )
        sample . stats_dfs . s3 [ "avg_depth_mj" ] = keepmj . mean ( )
        LOGGER . debug ( "mj depth {}" . format ( sample . stats_dfs . s3 [ "avg_depth_mj" ] ) )
        sample . stats_dfs . s3 [ "avg_depth_stat" ] = keepstat . mean ( )
        sample . stats_dfs . s3 [ "sd_depth_total" ] = depths . std ( )
        sample . stats_dfs . s3 [ "sd_depth_mj" ] = keepmj . std ( )
        sample . stats_dfs . s3 [ "sd_depth_stat" ] = keepstat . std ( )
else :
    print ( "    no clusters found for {}" . format ( sample . name ) )
# # Get some stats from the bam files
# # This is moderately hackish . samtools flagstat returns
# # the number of reads in the bam file as the first element
# # of the first line , this call makes this assumption .
if not data . paramsdict [ "assembly_method" ] == "denovo" :
    refmap_stats ( data , sample )
log_level = logging . getLevelName ( LOGGER . getEffectiveLevel ( ) )
if not log_level == "DEBUG" : # # Clean up loose files only if not in DEBUG
# # - edits / * derep , utemp , * utemp . sort , * htemp , * clust . gz
    derepfile = os . path . join ( data . dirs . edits , sample . name + "_derep.fastq" )
    mergefile = os . path . join ( data . dirs . edits , sample . name + "_merged_.fastq" )
    uhandle = os . path . join ( data . dirs . clusts , sample . name + ".utemp" )
    usort = os . path . join ( data . dirs . clusts , sample . name + ".utemp.sort" )
    hhandle = os . path . join ( data . dirs . clusts , sample . name + ".htemp" )
    clusters = os . path . join ( data . dirs . clusts , sample . name + ".clust.gz" )
    for f in [ derepfile , mergefile , uhandle , usort , hhandle , clusters ] :
        try :
            os . remove ( f )
        except :
            pass


signature:
def persistent_popen_align3 ( clusts , maxseqs = 200 , is_gbs = False ) :
    """keeps a persistent bash shell open and feeds it muscle alignments"""

implementation:
# # create a separate shell for running muscle in , this is much faster
# # than spawning a separate subprocess for each muscle call
proc = sps . Popen ( [ "bash" ] , stdin = sps . PIPE , stdout = sps . PIPE , universal_newlines = True )
# # iterate over clusters in this file until finished
aligned = [ ]
for clust in clusts : # # new alignment string for read1s and read2s
    align1 = ""
    align2 = ""
    # # don ' t bother aligning if only one seq
    if clust . count ( ">" ) == 1 :
        aligned . append ( clust . replace ( ">" , "" ) . strip ( ) )
    else : # # do we need to split the alignment ? ( is there a PE insert ? )
        try : # # make into list ( only read maxseqs lines , 2X cuz names )
            lclust = clust . split ( ) [ : maxseqs * 2 ]
            # # try to split cluster list at nnnn separator for each read
            lclust1 = list ( itertools . chain ( * zip ( lclust [ : : 2 ] , [ i . split ( "nnnn" ) [ 0 ] for i in lclust [ 1 : : 2 ] ] ) ) )
            lclust2 = list ( itertools . chain ( * zip ( lclust [ : : 2 ] , [ i . split ( "nnnn" ) [ 1 ] for i in lclust [ 1 : : 2 ] ] ) ) )
            # # put back into strings
            clust1 = "\n" . join ( lclust1 )
            clust2 = "\n" . join ( lclust2 )
            # # Align the first reads .
            # # The muscle command with alignment as stdin and / / as splitter
            cmd1 = "echo -e '{}' | {} -quiet -in - ; echo {}" . format ( clust1 , ipyrad . bins . muscle , "//" )
            # # send cmd1 to the bash shell
            print ( cmd1 , file = proc . stdin )
            # # read the stdout by line until splitter is reached
            # # meaning that the alignment is finished .
            for line in iter ( proc . stdout . readline , '//\n' ) :
                align1 += line
            # # Align the second reads .
            # # The muscle command with alignment as stdin and / / as splitter
            cmd2 = "echo -e '{}' | {} -quiet -in - ; echo {}" . format ( clust2 , ipyrad . bins . muscle , "//" )
            # # send cmd2 to the bash shell
            print ( cmd2 , file = proc . stdin )
            # # read the stdout by line until splitter is reached
            # # meaning that the alignment is finished .
            for line in iter ( proc . stdout . readline , '//\n' ) :
                align2 += line
            # # join up aligned read1 and read2 and ensure names order matches
            la1 = align1 [ 1 : ] . split ( "\n>" )
            la2 = align2 [ 1 : ] . split ( "\n>" )
            dalign1 = dict ( [ i . split ( "\n" , 1 ) for i in la1 ] )
            dalign2 = dict ( [ i . split ( "\n" , 1 ) for i in la2 ] )
            align1 = [ ]
            try :
                keys = sorted ( dalign1 . keys ( ) , key = DEREP , reverse = True )
            except ValueError as inst : # # Lines is empty . This means the call to muscle alignment failed .
            # # Not sure how to handle this , but it happens only very rarely .
                LOGGER . error ( "Muscle alignment failed: Bad clust - {}\nBad lines - {}" . format ( clust , lines ) )
                continue
            # # put seed at top of alignment
            seed = [ i for i in keys if i . split ( ";" ) [ - 1 ] [ 0 ] == "*" ] [ 0 ]
            keys . pop ( keys . index ( seed ) )
            keys = [ seed ] + keys
            for key in keys :
                align1 . append ( "\n" . join ( [ key , dalign1 [ key ] . replace ( "\n" , "" ) + "nnnn" + dalign2 [ key ] . replace ( "\n" , "" ) ] ) )
            # # append aligned cluster string
            aligned . append ( "\n" . join ( align1 ) . strip ( ) )
        # # Malformed clust . Dictionary creation with only 1 element will raise .
        except ValueError as inst :
            LOGGER . debug ( "Bad PE cluster - {}\nla1 - {}\nla2 - {}" . format ( clust , la1 , la2 ) )
        # # Either reads are SE , or at least some pairs are merged .
        except IndexError : # # limit the number of input seqs
            lclust = "\n" . join ( clust . split ( ) [ : maxseqs * 2 ] )
            # # the muscle command with alignment as stdin and / / as splitter
            cmd = "echo -e '{}' | {} -quiet -in - ; echo {}" . format ( lclust , ipyrad . bins . muscle , "//" )
            # # send cmd to the bash shell ( TODO : PIPE could overflow here ! )
            print ( cmd , file = proc . stdin )
            # # read the stdout by line until / / is reached . This BLOCKS .
            for line in iter ( proc . stdout . readline , '//\n' ) :
                align1 += line
            # # remove ' > ' from names , and ' \ n ' from inside long seqs
            lines = align1 [ 1 : ] . split ( "\n>" )
            try : # # find seed of the cluster and put it on top .
                seed = [ i for i in lines if i . split ( ";" ) [ - 1 ] [ 0 ] == "*" ] [ 0 ]
                lines . pop ( lines . index ( seed ) )
                lines = [ seed ] + sorted ( lines , key = DEREP , reverse = True )
            except ValueError as inst : # # Lines is empty . This means the call to muscle alignment failed .
            # # Not sure how to handle this , but it happens only very rarely .
                LOGGER . error ( "Muscle alignment failed: Bad clust - {}\nBad lines - {}" . format ( clust , lines ) )
                continue
            # # format remove extra newlines from muscle
            aa = [ i . split ( "\n" , 1 ) for i in lines ]
            align1 = [ i [ 0 ] + '\n' + "" . join ( [ j . replace ( "\n" , "" ) for j in i [ 1 : ] ] ) for i in aa ]
            # # trim edges in sloppy gbs / ezrad data . Maybe relevant to other types too . . .
            if is_gbs :
                align1 = gbs_trim ( align1 )
            # # append to aligned
            aligned . append ( "\n" . join ( align1 ) . strip ( ) )
# cleanup
proc . stdout . close ( )
if proc . stderr :
    proc . stderr . close ( )
proc . stdin . close ( )
proc . wait ( )
# # return the aligned clusters
return aligned


signature:
def gbs_trim ( align1 ) :
    """No reads can go past the left of the seed , or right of the least extended
 reverse complement match . Example below . m is a match . u is an area where
 lots of mismatches typically occur . The cut sites are shown .
 Original locus *
 Seed TGCAG * * * * * - - - - -
 Forward - match TGCAGmmmmm - - - - -
 Forward - match TGCAGmmmmm - - - - -
 Forward - match TGCAGmmmmm - - - - -
 Revcomp - match - - - - - mmmmmCTGCAuuuuu
 Revcomp - match - - - - - mmmmmCTGCAuuuuu
 Revcomp - match - - - - - mmmmmCTGCA
 Revcomp - match - - - - - mmmmmCTGCAuuuuu
 Trimmed locus *
 Seed TGCAG * * * * * - - - - -
 Forward - match TGCAGmmmmm - - - - -
 Forward - match TGCAGmmmmm - - - - -
 Forward - match TGCAGmmmmm - - - - -
 Revcomp - match - - - - - mmmmm
 Revcomp - match - - - - - mmmmmCTGCA
 Revcomp - match - - - - - mmmmm
 Revcomp - match - - - - - mmmmm"""

implementation:
leftmost = rightmost = None
dd = { k : v for k , v in [ j . rsplit ( "\n" , 1 ) for j in align1 ] }
seed = [ i for i in dd . keys ( ) if i . rsplit ( ";" ) [ - 1 ] [ 0 ] == "*" ] [ 0 ]
leftmost = [ i != "-" for i in dd [ seed ] ] . index ( True )
revs = [ i for i in dd . keys ( ) if i . rsplit ( ";" ) [ - 1 ] [ 0 ] == "-" ]
if revs :
    subright = max ( [ [ i != "-" for i in seq [ : : - 1 ] ] . index ( True ) for seq in [ dd [ i ] for i in revs ] ] )
else :
    subright = 0
rightmost = len ( dd [ seed ] ) - subright
# # if locus got clobbered then print place - holder NNN
names , seqs = zip ( * [ i . rsplit ( "\n" , 1 ) for i in align1 ] )
if rightmost > leftmost :
    newalign1 = [ n + "\n" + i [ leftmost : rightmost ] for i , n in zip ( seqs , names ) ]
else :
    newalign1 = [ n + "\nNNN" for i , n in zip ( seqs , names ) ]
return newalign1


signature:
def align_and_parse ( handle , max_internal_indels = 5 , is_gbs = False ) :
    """much faster implementation for aligning chunks"""

implementation:
# # data are already chunked , read in the whole thing . bail if no data .
try :
    with open ( handle , 'rb' ) as infile :
        clusts = infile . read ( ) . split ( "//\n//\n" )
        # # remove any empty spots
        clusts = [ i for i in clusts if i ]
        # # Skip entirely empty chunks
        if not clusts :
            raise IPyradError
except ( IOError , IPyradError ) :
    LOGGER . debug ( "skipping empty chunk - {}" . format ( handle ) )
    return 0
# # count discarded clusters for printing to stats later
highindels = 0
# # iterate over clusters sending each to muscle , splits and aligns pairs
try :
    aligned = persistent_popen_align3 ( clusts , 200 , is_gbs )
except Exception as inst :
    LOGGER . debug ( "Error in handle - {} - {}" . format ( handle , inst ) )
    # raise IPyradWarningExit ( " error hrere { } " . format ( inst ) )
    aligned = [ ]
# # store good alignments to be written to file
refined = [ ]
# # filter and trim alignments
for clust in aligned : # # check for too many internal indels
    filtered = aligned_indel_filter ( clust , max_internal_indels )
    # # reverse complement matches . No longer implemented .
    # filtered = overshoot _ filter ( clust )
    # # finally , add to outstack if alignment is good
    if not filtered :
        refined . append ( clust )
        # " \ n " . join ( stack ) )
    else :
        highindels += 1
# # write to file after
if refined :
    outhandle = handle . rsplit ( "." , 1 ) [ 0 ] + ".aligned"
    with open ( outhandle , 'wb' ) as outfile :
        outfile . write ( "\n//\n//\n" . join ( refined ) + "\n" )
# # remove the old tmp file
log_level = logging . getLevelName ( LOGGER . getEffectiveLevel ( ) )
if not log_level == "DEBUG" :
    os . remove ( handle )
return highindels


signature:
def aligned_indel_filter ( clust , max_internal_indels ) :
    """checks for too many internal indels in muscle aligned clusters"""

implementation:
# # make into list
lclust = clust . split ( )
# # paired or not
try :
    seq1 = [ i . split ( "nnnn" ) [ 0 ] for i in lclust [ 1 : : 2 ] ]
    seq2 = [ i . split ( "nnnn" ) [ 1 ] for i in lclust [ 1 : : 2 ] ]
    intindels1 = [ i . rstrip ( "-" ) . lstrip ( "-" ) . count ( "-" ) for i in seq1 ]
    intindels2 = [ i . rstrip ( "-" ) . lstrip ( "-" ) . count ( "-" ) for i in seq2 ]
    intindels = intindels1 + intindels2
    if max ( intindels ) > max_internal_indels :
        return 1
except IndexError :
    seq1 = lclust [ 1 : : 2 ]
    intindels = [ i . rstrip ( "-" ) . lstrip ( "-" ) . count ( "-" ) for i in seq1 ]
    if max ( intindels ) > max_internal_indels :
        return 1
return 0


signature:
def build_clusters ( data , sample , maxindels ) :
    """Combines information from . utemp and . htemp files to create . clust files ,
 which contain un - aligned clusters . Hits to seeds are only kept in the
 cluster if the number of internal indels is less than ' maxindels ' .
 By default , we set maxindels = 6 for this step ( within - sample clustering ) ."""

implementation:
# # If reference assembly then here we ' re clustering the unmapped reads
if "reference" in data . paramsdict [ "assembly_method" ] :
    derepfile = os . path . join ( data . dirs . edits , sample . name + "-refmap_derep.fastq" )
else :
    derepfile = os . path . join ( data . dirs . edits , sample . name + "_derep.fastq" )
# # i / o vsearch files
uhandle = os . path . join ( data . dirs . clusts , sample . name + ".utemp" )
usort = os . path . join ( data . dirs . clusts , sample . name + ".utemp.sort" )
hhandle = os . path . join ( data . dirs . clusts , sample . name + ".htemp" )
# # create an output file to write clusters to
sample . files . clusters = os . path . join ( data . dirs . clusts , sample . name + ".clust.gz" )
clustsout = gzip . open ( sample . files . clusters , 'wb' )
# # Sort the uhandle file so we can read through matches efficiently
cmd = [ "sort" , "-k" , "2" , uhandle , "-o" , usort ]
proc = sps . Popen ( cmd , close_fds = True )
_ = proc . communicate ( ) [ 0 ]
# # load ALL derep reads into a dictionary ( this can be a few GB of RAM )
# # and is larger if names are larger . We are grabbing two lines at a time .
alldereps = { }
with open ( derepfile , 'rb' ) as ioderep :
    dereps = itertools . izip ( * [ iter ( ioderep ) ] * 2 )
    for namestr , seq in dereps :
        nnn , sss = [ i . strip ( ) for i in namestr , seq ]
        alldereps [ nnn [ 1 : ] ] = sss
# # store observed seeds ( this could count up to > million in bad data sets )
seedsseen = set ( )
# # Iterate through the usort file grabbing matches to build clusters
with open ( usort , 'rb' ) as insort : # # iterator , seed null , seqlist null
    isort = iter ( insort )
    lastseed = 0
    fseqs = [ ]
    seqlist = [ ]
    seqsize = 0
    while 1 : # # grab the next line
        try :
            hit , seed , _ , ind , ori , _ = isort . next ( ) . strip ( ) . split ( )
            LOGGER . debug ( ">{} {} {}" . format ( hit , seed , ori , seq ) )
        except StopIteration :
            break
        # # same seed , append match
        if seed != lastseed :
            seedsseen . add ( seed )
            # # store the last cluster ( fseq ) , count it , and clear fseq
            if fseqs : # # sort fseqs by derep after pulling out the seed
                fseqs = [ fseqs [ 0 ] ] + sorted ( fseqs [ 1 : ] , key = lambda x : int ( x . split ( ";size=" ) [ 1 ] . split ( ";" ) [ 0 ] ) , reverse = True )
                seqlist . append ( "\n" . join ( fseqs ) )
                seqsize += 1
                fseqs = [ ]
            # # occasionally write / dump stored clusters to file and clear mem
            if not seqsize % 10000 :
                if seqlist :
                    clustsout . write ( "\n//\n//\n" . join ( seqlist ) + "\n//\n//\n" )
                    # # reset list and counter
                    seqlist = [ ]
            # # store the new seed on top of fseq list
            fseqs . append ( ">{}*\n{}" . format ( seed , alldereps [ seed ] ) )
            lastseed = seed
        # # add match to the seed
        # # revcomp if orientation is reversed ( comp preserves nnnn )
        if ori == "-" :
            seq = comp ( alldereps [ hit ] ) [ : : - 1 ]
        else :
            seq = alldereps [ hit ]
        # # only save if not too many indels
        if int ( ind ) <= maxindels :
            fseqs . append ( ">{}{}\n{}" . format ( hit , ori , seq ) )
        else :
            LOGGER . info ( "filtered by maxindels: %s %s" , ind , seq )
# # write whatever is left over to the clusts file
if fseqs :
    seqlist . append ( "\n" . join ( fseqs ) )
if seqlist :
    clustsout . write ( "\n//\n//\n" . join ( seqlist ) + "\n//\n//\n" )
# # now write the seeds that had no hits . Make dict from htemp
with open ( hhandle , 'rb' ) as iotemp :
    nohits = itertools . izip ( * [ iter ( iotemp ) ] * 2 )
    seqlist = [ ]
    seqsize = 0
    while 1 :
        try :
            nnn , _ = [ i . strip ( ) for i in nohits . next ( ) ]
        except StopIteration :
            break
        # # occasionally write to file
        if not seqsize % 10000 :
            if seqlist :
                clustsout . write ( "\n//\n//\n" . join ( seqlist ) + "\n//\n//\n" )
                # # reset list and counter
                seqlist = [ ]
        # # append to list if new seed
        if nnn [ 1 : ] not in seedsseen :
            seqlist . append ( "{}*\n{}" . format ( nnn , alldereps [ nnn [ 1 : ] ] ) )
            seqsize += 1
# # write whatever is left over to the clusts file
if seqlist :
    clustsout . write ( "\n//\n//\n" . join ( seqlist ) )
# + " \ n / / \ n / / \ n " )
# # close the file handle
clustsout . close ( )
del alldereps


signature:
def setup_dirs ( data ) :
    """sets up directories for step3 data"""

implementation:
# # make output folder for clusters
pdir = os . path . realpath ( data . paramsdict [ "project_dir" ] )
data . dirs . clusts = os . path . join ( pdir , "{}_clust_{}" . format ( data . name , data . paramsdict [ "clust_threshold" ] ) )
if not os . path . exists ( data . dirs . clusts ) :
    os . mkdir ( data . dirs . clusts )
# # make a tmpdir for align files
data . tmpdir = os . path . abspath ( os . path . expanduser ( os . path . join ( pdir , data . name + '-tmpalign' ) ) )
if not os . path . exists ( data . tmpdir ) :
    os . mkdir ( data . tmpdir )
# # If ref mapping , init samples and make the refmapping output directory .
if not data . paramsdict [ "assembly_method" ] == "denovo" : # # make output directory for read mapping process
    data . dirs . refmapping = os . path . join ( pdir , "{}_refmapping" . format ( data . name ) )
    if not os . path . exists ( data . dirs . refmapping ) :
        os . mkdir ( data . dirs . refmapping )


signature:
def new_apply_jobs ( data , samples , ipyclient , nthreads , maxindels , force ) :
    """Create a DAG of prealign jobs to be run in order for each sample . Track
 Progress , report errors . Each assembly method has a slightly different
 DAG setup , calling different functions ."""

implementation:
# # is datatype gbs ? used in alignment - trimming by align _ and _ parse ( )
is_gbs = bool ( "gbs" in data . paramsdict [ "datatype" ] )
# # Two view objects , threaded and unthreaded
lbview = ipyclient . load_balanced_view ( )
start = time . time ( )
elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
firstfunc = "derep_concat_split"
printstr = " {}    | {} | s3 |" . format ( PRINTSTR [ firstfunc ] , elapsed )
# printstr = " { } | { } | s3 | " . format ( PRINTSTR [ ] , elapsed )
progressbar ( 10 , 0 , printstr , spacer = data . _spacer )
# # TODO : for HPC systems this should be done to make sure targets are spread
# # among different nodes .
if nthreads :
    if nthreads < len ( ipyclient . ids ) :
        thview = ipyclient . load_balanced_view ( targets = ipyclient . ids [ : : nthreads ] )
    elif nthreads == 1 :
        thview = ipyclient . load_balanced_view ( )
    else :
        if len ( ipyclient ) > 40 :
            thview = ipyclient . load_balanced_view ( targets = ipyclient . ids [ : : 4 ] )
        else :
            thview = ipyclient . load_balanced_view ( targets = ipyclient . ids [ : : 2 ] )
# # get list of jobs / dependencies as a DAG for all pre - align funcs .
dag , joborder = build_dag ( data , samples )
# # dicts for storing submitted jobs and results
results = { }
# # submit jobs to the engines in single or threaded views . The topological
# # sort makes sure jobs are input with all dependencies found .
for node in nx . topological_sort ( dag ) : # # get list of async results leading to this job
    deps = [ results . get ( n ) for n in dag . predecessors ( node ) ]
    deps = ipp . Dependency ( dependencies = deps , failure = True )
    # # get func , sample , and args for this func ( including [ data , sample ] )
    funcstr , chunk , sname = node . split ( "-" , 2 )
    func = FUNCDICT [ funcstr ]
    sample = data . samples [ sname ]
    # # args vary depending on the function
    if funcstr in [ "derep_concat_split" , "cluster" ] :
        args = [ data , sample , nthreads , force ]
    elif funcstr in [ "mapreads" ] :
        args = [ data , sample , nthreads , force ]
    elif funcstr in [ "build_clusters" ] :
        args = [ data , sample , maxindels ]
    elif funcstr in [ "muscle_align" ] :
        handle = os . path . join ( data . tmpdir , "{}_chunk_{}.ali" . format ( sample . name , chunk ) )
        args = [ handle , maxindels , is_gbs ]
    else :
        args = [ data , sample ]
    # submit and store AsyncResult object . Some jobs are threaded .
    if nthreads and ( funcstr in THREADED_FUNCS ) : # LOGGER . info ( ' submitting % s to % s - threaded view ' , funcstr , nthreads )
        with thview . temp_flags ( after = deps , block = False ) :
            results [ node ] = thview . apply ( func , * args )
    else : # LOGGER . info ( ' submitting % s to single - threaded view ' , funcstr )
        with lbview . temp_flags ( after = deps , block = False ) :
            results [ node ] = lbview . apply ( func , * args )
# # track jobs as they finish , abort if someone fails . This blocks here
# # until all jobs are done . Keep track of which samples have failed so
# # we only print the first error message .
sfailed = set ( )
for funcstr in joborder + [ "muscle_align" , "reconcat" ] :
    errfunc , sfails , msgs = trackjobs ( funcstr , results , spacer = data . _spacer )
    LOGGER . info ( "{}-{}-{}" . format ( errfunc , sfails , msgs ) )
    if errfunc :
        for sidx in xrange ( len ( sfails ) ) :
            sname = sfails [ sidx ]
            errmsg = msgs [ sidx ]
            if sname not in sfailed :
                print ( "  sample [{}] failed. See error in ./ipyrad_log.txt" . format ( sname ) )
                LOGGER . error ( "sample [%s] failed in step [%s]; error: %s" , sname , errfunc , errmsg )
                sfailed . add ( sname )
# # Cleanup of successful samples , skip over failed samples
badaligns = { }
for sample in samples : # # The muscle _ align step returns the number of excluded bad alignments
    for async in results :
        func , chunk , sname = async . split ( "-" , 2 )
        if ( func == "muscle_align" ) and ( sname == sample . name ) :
            if results [ async ] . successful ( ) :
                badaligns [ sample ] = int ( results [ async ] . get ( ) )
# # for the samples that were successful :
for sample in badaligns : # # store the result
    sample . stats_dfs . s3 . filtered_bad_align = badaligns [ sample ]
    # # store all results
    try :
        sample_cleanup ( data , sample )
    except Exception as inst :
        msg = "  Sample {} failed this step. See ipyrad_log.txt.\
                  " . format ( sample . name )
        print ( msg )
        LOGGER . error ( "%s - %s" , sample . name , inst )
# # store the results to data
data_cleanup ( data )


signature:
def build_dag ( data , samples ) :
    """build a directed acyclic graph describing jobs to be run in order ."""

implementation:
# # Create DAGs for the assembly method being used , store jobs in nodes
snames = [ i . name for i in samples ]
dag = nx . DiGraph ( )
# # get list of pre - align jobs from globals based on assembly method
joborder = JOBORDER [ data . paramsdict [ "assembly_method" ] ]
# # WHICH JOBS TO RUN : iterate over the sample names
for sname in snames : # # append pre - align job for each sample to nodes list
    for func in joborder :
        dag . add_node ( "{}-{}-{}" . format ( func , 0 , sname ) )
    # # append align func jobs , each will have max 10
    for chunk in xrange ( 10 ) :
        dag . add_node ( "{}-{}-{}" . format ( "muscle_align" , chunk , sname ) )
    # # append final reconcat jobs
    dag . add_node ( "{}-{}-{}" . format ( "reconcat" , 0 , sname ) )
# # ORDER OF JOBS : add edges / dependency between jobs : ( first - this , then - that )
for sname in snames :
    for sname2 in snames : # # enforce that clust / map cannot start until derep is done for ALL
    # # samples . This is b / c . . .
        dag . add_edge ( "{}-{}-{}" . format ( joborder [ 0 ] , 0 , sname2 ) , "{}-{}-{}" . format ( joborder [ 1 ] , 0 , sname ) )
    # # add remaining pre - align jobs
    for idx in xrange ( 2 , len ( joborder ) ) :
        dag . add_edge ( "{}-{}-{}" . format ( joborder [ idx - 1 ] , 0 , sname ) , "{}-{}-{}" . format ( joborder [ idx ] , 0 , sname ) )
    # # Add 10 align jobs , none of which can start until all chunker jobs
    # # are finished . Similarly , reconcat jobs cannot start until all align
    # # jobs are finished .
    for sname2 in snames :
        for chunk in range ( 10 ) :
            dag . add_edge ( "{}-{}-{}" . format ( "muscle_chunker" , 0 , sname2 ) , "{}-{}-{}" . format ( "muscle_align" , chunk , sname ) )
            # # add that the final reconcat job can ' t start until after
            # # each chunk of its own sample has finished aligning .
            dag . add_edge ( "{}-{}-{}" . format ( "muscle_align" , chunk , sname ) , "{}-{}-{}" . format ( "reconcat" , 0 , sname ) )
# # return the dag
return dag , joborder


signature:
def _plot_dag ( dag , results , snames ) :
    """makes plot to help visualize the DAG setup . For developers only ."""

implementation:
try :
    import matplotlib . pyplot as plt
    from matplotlib . dates import date2num
    from matplotlib . cm import gist_rainbow
    # # first figure is dag layout
    plt . figure ( "dag_layout" , figsize = ( 10 , 10 ) )
    nx . draw ( dag , pos = nx . spring_layout ( dag ) , node_color = 'pink' , with_labels = True )
    plt . savefig ( "./dag_layout.png" , bbox_inches = 'tight' , dpi = 200 )
    # # second figure is times for steps
    pos = { }
    colors = { }
    for node in dag : # jobkey = " { } - { } " . format ( node , sample )
        mtd = results [ node ] . metadata
        start = date2num ( mtd . started )
        # runtime = date2num ( md . completed ) # - start
        # # sample id to separate samples on x - axis
        _ , _ , sname = node . split ( "-" , 2 )
        sid = snames . index ( sname )
        # # 1e6 to separate on y - axis
        pos [ node ] = ( start + sid , start * 1e6 )
        colors [ node ] = mtd . engine_id
    # # x just spaces out samples ;
    # # y is start time of each job with edge leading to next job
    # # color is the engine that ran the job
    # # all jobs were submitted as 3 second wait times
    plt . figure ( "dag_starttimes" , figsize = ( 10 , 16 ) )
    nx . draw ( dag , pos , node_list = colors . keys ( ) , node_color = colors . values ( ) , cmap = gist_rainbow , with_labels = True )
    plt . savefig ( "./dag_starttimes.png" , bbox_inches = 'tight' , dpi = 200 )
except Exception as inst :
    LOGGER . warning ( inst )


signature:
def trackjobs ( func , results , spacer ) :
    """Blocks and prints progress for just the func being requested from a list
 of submitted engine jobs . Returns whether any of the jobs failed .
 func = str
 results = dict of asyncs"""

implementation:
# # TODO : try to insert a better way to break on KBD here .
LOGGER . info ( "inside trackjobs of %s" , func )
# # get just the jobs from results that are relevant to this func
asyncs = [ ( i , results [ i ] ) for i in results if i . split ( "-" , 2 ) [ 0 ] == func ]
# # progress bar
start = time . time ( )
while 1 : # # how many of this func have finished so far
    ready = [ i [ 1 ] . ready ( ) for i in asyncs ]
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    printstr = " {}    | {} | s3 |" . format ( PRINTSTR [ func ] , elapsed )
    progressbar ( len ( ready ) , sum ( ready ) , printstr , spacer = spacer )
    time . sleep ( 0.1 )
    if len ( ready ) == sum ( ready ) :
        print ( "" )
        break
sfails = [ ]
errmsgs = [ ]
for job in asyncs :
    if not job [ 1 ] . successful ( ) :
        sfails . append ( job [ 0 ] )
        errmsgs . append ( job [ 1 ] . result ( ) )
return func , sfails , errmsgs


signature:
def declone_3rad ( data , sample ) :
    """3rad uses random adapters to identify pcr duplicates . We will
 remove pcr dupes here . Basically append the radom adapter to
 each sequence , do a regular old vsearch derep , then trim
 off the adapter , and push it down the pipeline . This will
 remove all identical seqs with identical random i5 adapters ."""

implementation:
LOGGER . info ( "Entering declone_3rad - {}" . format ( sample . name ) )
# # Append i5 adapter to the head of each read . Merged file is input , and
# # still has fq qual score so also have to append several qscores for the
# # adapter bases . Open the merge file , get quarts , go through each read
# # and append the necessary stuff .
adapter_seqs_file = tempfile . NamedTemporaryFile ( mode = 'wb' , delete = False , dir = data . dirs . edits , suffix = "_append_adapters_.fastq" )
try :
    with open ( sample . files . edits [ 0 ] [ 0 ] ) as infile :
        quarts = itertools . izip ( * [ iter ( infile ) ] * 4 )
        # # a list to store until writing
        writing = [ ]
        counts = 0
        while 1 :
            try :
                read = quarts . next ( )
            except StopIteration :
                break
            # # Split on + , get [ 1 ] , split on " _ " ( can be either _ r1 or
            # # _ m1 if merged reads ) and get [ 0 ] for the i5
            # # prepend " EEEEE " as qscore for the adapters
            i5 = read [ 0 ] . split ( "+" ) [ 1 ] . split ( "_" ) [ 0 ]
            # # If any non ACGT in the i5 then drop this sequence
            if 'N' in i5 :
                continue
            writing . append ( "\n" . join ( [ read [ 0 ] . strip ( ) , i5 + read [ 1 ] . strip ( ) , read [ 2 ] . strip ( ) , "E" * 8 + read [ 3 ] . strip ( ) ] ) )
            # # Write the data in chunks
            counts += 1
            if not counts % 1000 :
                adapter_seqs_file . write ( "\n" . join ( writing ) + "\n" )
                writing = [ ]
        if writing :
            adapter_seqs_file . write ( "\n" . join ( writing ) )
            adapter_seqs_file . close ( )
    tmp_outfile = tempfile . NamedTemporaryFile ( mode = 'wb' , delete = False , dir = data . dirs . edits , suffix = "_decloned_w_adapters_.fastq" )
    # # Close the tmp file bcz vsearch will write to it by name , then
    # # we will want to reopen it to read from it .
    tmp_outfile . close ( )
    # # Derep the data ( adapters + seq )
    derep_and_sort ( data , adapter_seqs_file . name , os . path . join ( data . dirs . edits , tmp_outfile . name ) , 2 )
    # # Remove adapters from head of sequence and write out
    # # tmp _ outfile is now the input file for the next step
    # # first vsearch derep discards the qscore so we iterate
    # # by pairs
    with open ( tmp_outfile . name ) as infile :
        with open ( os . path . join ( data . dirs . edits , sample . name + "_declone.fastq" ) , 'wb' ) as outfile :
            duo = itertools . izip ( * [ iter ( infile ) ] * 2 )
            # # a list to store until writing
            writing = [ ]
            counts2 = 0
            while 1 :
                try :
                    read = duo . next ( )
                except StopIteration :
                    break
                # # Peel off the adapters . There ' s probably a faster
                # # way of doing this .
                writing . append ( "\n" . join ( [ read [ 0 ] . strip ( ) , read [ 1 ] . strip ( ) [ 8 : ] ] ) )
                # # Write the data in chunks
                counts2 += 1
                if not counts2 % 1000 :
                    outfile . write ( "\n" . join ( writing ) + "\n" )
                    writing = [ ]
            if writing :
                outfile . write ( "\n" . join ( writing ) )
                outfile . close ( )
    LOGGER . info ( "Removed pcr duplicates from {} - {}" . format ( sample . name , counts - counts2 ) )
except Exception as inst :
    raise IPyradError ( "    Caught error while decloning " + "3rad data - {}" . format ( inst ) )
finally : # # failed samples will cause tmp file removal to raise .
# # just ignore it .
    try : # # Clean up temp files
        if os . path . exists ( adapter_seqs_file . name ) :
            os . remove ( adapter_seqs_file . name )
        if os . path . exists ( tmp_outfile . name ) :
            os . remove ( tmp_outfile . name )
    except Exception as inst :
        pass


signature:
def derep_and_sort ( data , infile , outfile , nthreads ) :
    """Dereplicates reads and sorts so reads that were highly replicated are at
 the top , and singletons at bottom , writes output to derep file . Paired
 reads are dereplicated as one concatenated read and later split again .
 Updated this function to take infile and outfile to support the double
 dereplication that we need for 3rad ( 5/29/15 iao ) ."""

implementation:
# # datatypes options
strand = "plus"
if "gbs" in data . paramsdict [ "datatype" ] or "2brad" in data . paramsdict [ "datatype" ] :
    strand = "both"
# # pipe in a gzipped file
if infile . endswith ( ".gz" ) :
    catcmd = [ "gunzip" , "-c" , infile ]
else :
    catcmd = [ "cat" , infile ]
# # do dereplication with vsearch
cmd = [ ipyrad . bins . vsearch , "--derep_fulllength" , "-" , "--strand" , strand , "--output" , outfile , "--threads" , str ( nthreads ) , "--fasta_width" , str ( 0 ) , "--fastq_qmax" , "1000" , "--sizeout" , "--relabel_md5" , ]
LOGGER . info ( "derep cmd %s" , " " . join ( cmd ) )
# # run vsearch
proc1 = sps . Popen ( catcmd , stderr = sps . STDOUT , stdout = sps . PIPE , close_fds = True )
proc2 = sps . Popen ( cmd , stdin = proc1 . stdout , stderr = sps . STDOUT , stdout = sps . PIPE , close_fds = True )
try :
    errmsg = proc2 . communicate ( ) [ 0 ]
except KeyboardInterrupt :
    LOGGER . info ( "interrupted during dereplication" )
    raise KeyboardInterrupt ( )
if proc2 . returncode :
    LOGGER . error ( "error inside derep_and_sort %s" , errmsg )
    raise IPyradWarningExit ( errmsg )


signature:
def data_cleanup ( data ) :
    """cleanup / statswriting function for Assembly obj"""

implementation:
data . stats_dfs . s3 = data . _build_stat ( "s3" )
data . stats_files . s3 = os . path . join ( data . dirs . clusts , "s3_cluster_stats.txt" )
with io . open ( data . stats_files . s3 , 'w' ) as outfile :
    data . stats_dfs . s3 . to_string ( buf = outfile , formatters = { 'merged_pairs' : '{:.0f}' . format , 'clusters_total' : '{:.0f}' . format , 'clusters_hidepth' : '{:.0f}' . format , 'filtered_bad_align' : '{:.0f}' . format , 'avg_depth_stat' : '{:.2f}' . format , 'avg_depth_mj' : '{:.2f}' . format , 'avg_depth_total' : '{:.2f}' . format , 'sd_depth_stat' : '{:.2f}' . format , 'sd_depth_mj' : '{:.2f}' . format , 'sd_depth_total' : '{:.2f}' . format } )


signature:
def concat_multiple_edits ( data , sample ) :
    """if multiple fastq files were appended into the list of fastqs for samples
 then we merge them here before proceeding ."""

implementation:
# # if more than one tuple in fastq list
if len ( sample . files . edits ) > 1 : # # create a cat command to append them all ( doesn ' t matter if they
# # are gzipped , cat still works ) . Grab index 0 of tuples for R1s .
    cmd1 = [ "cat" ] + [ i [ 0 ] for i in sample . files . edits ]
    # # write to new concat handle
    conc1 = os . path . join ( data . dirs . edits , sample . name + "_R1_concatedit.fq.gz" )
    with open ( conc1 , 'w' ) as cout1 :
        proc1 = sps . Popen ( cmd1 , stderr = sps . STDOUT , stdout = cout1 , close_fds = True )
        res1 = proc1 . communicate ( ) [ 0 ]
    if proc1 . returncode :
        raise IPyradWarningExit ( "error in: %s, %s" , cmd1 , res1 )
    # # Only set conc2 if R2 actually exists
    conc2 = 0
    if os . path . exists ( str ( sample . files . edits [ 0 ] [ 1 ] ) ) :
        cmd2 = [ "cat" ] + [ i [ 1 ] for i in sample . files . edits ]
        conc2 = os . path . join ( data . dirs . edits , sample . name + "_R2_concatedit.fq.gz" )
        with gzip . open ( conc2 , 'w' ) as cout2 :
            proc2 = sps . Popen ( cmd2 , stderr = sps . STDOUT , stdout = cout2 , close_fds = True )
            res2 = proc2 . communicate ( ) [ 0 ]
        if proc2 . returncode :
            raise IPyradWarningExit ( "error in: %s, %s" , cmd2 , res2 )
    # # store new file handles
    sample . files . edits = [ ( conc1 , conc2 ) ]
return sample . files . edits


signature:
def cluster ( data , sample , nthreads , force ) :
    """Calls vsearch for clustering . cov varies by data type , values were chosen
 based on experience , but could be edited by users"""

implementation:
# # get the dereplicated reads
if "reference" in data . paramsdict [ "assembly_method" ] :
    derephandle = os . path . join ( data . dirs . edits , sample . name + "-refmap_derep.fastq" )
    # # In the event all reads for all samples map successfully then clustering
    # # the unmapped reads makes no sense , so just bail out .
    if not os . stat ( derephandle ) . st_size : # # In this case you do have to create empty , dummy vsearch output
    # # files so building _ clusters will not fail .
        uhandle = os . path . join ( data . dirs . clusts , sample . name + ".utemp" )
        usort = os . path . join ( data . dirs . clusts , sample . name + ".utemp.sort" )
        hhandle = os . path . join ( data . dirs . clusts , sample . name + ".htemp" )
        for f in [ uhandle , usort , hhandle ] :
            open ( f , 'a' ) . close ( )
        return
else :
    derephandle = os . path . join ( data . dirs . edits , sample . name + "_derep.fastq" )
# # create handles for the outfiles
uhandle = os . path . join ( data . dirs . clusts , sample . name + ".utemp" )
temphandle = os . path . join ( data . dirs . clusts , sample . name + ".htemp" )
# # If derep file doesn ' t exist then bail out
if not os . path . isfile ( derephandle ) :
    LOGGER . warn ( "Bad derephandle - {}" . format ( derephandle ) )
    raise IPyradError ( "Input file for clustering doesn't exist - {}" . format ( derephandle ) )
# # testing one sample fail
# if sample . name = = " 1C _ 0 " :
# # datatype specific optimization
# # minsl : the percentage of the seed that must be matched
# # smaller values for RAD / ddRAD where we might want to combine , say 50bp
# # reads and 100bp reads in the same analysis .
# # query _ cov : the percentage of the query sequence that must match seed
# # smaller values are needed for gbs where only the tips might overlap
# # larger values for pairgbs where they should overlap near completely
# # small minsl and high query cov allows trimmed reads to match to untrim
# # seed for rad / ddrad / pairddrad .
strand = "plus"
cov = 0.75
minsl = 0.5
if data . paramsdict [ "datatype" ] in [ "gbs" , "2brad" ] :
    strand = "both"
    cov = 0.5
    minsl = 0.5
elif data . paramsdict [ "datatype" ] == 'pairgbs' :
    strand = "both"
    cov = 0.75
    minsl = 0.75
# # If this value is not null ( which is the default ) then override query cov
if data . _hackersonly [ "query_cov" ] :
    cov = str ( data . _hackersonly [ "query_cov" ] )
    assert float ( cov ) <= 1 , "query_cov must be <= 1.0"
# # get call string
cmd = [ ipyrad . bins . vsearch , "-cluster_smallmem" , derephandle , "-strand" , strand , "-query_cov" , str ( cov ) , "-id" , str ( data . paramsdict [ "clust_threshold" ] ) , "-minsl" , str ( minsl ) , "-userout" , uhandle , "-userfields" , "query+target+id+gaps+qstrand+qcov" , "-maxaccepts" , "1" , "-maxrejects" , "0" , "-threads" , str ( nthreads ) , "-notmatched" , temphandle , "-fasta_width" , "0" , "-fastq_qmax" , "100" , "-fulldp" , "-usersort" ]
# # not sure what the benefit of this option is exactly , needs testing ,
# # might improve indel detection on left side , but we don ' t want to enforce
# # aligning on left side if not necessarily , since quality trimmed reads
# # might lose bases on left side in step2 and no longer align .
# if data . paramsdict [ " datatype " ] in [ " rad " , " ddrad " , " pairddrad " ] :
# cmd + = [ " - leftjust " ]
# # run vsearch
LOGGER . debug ( "%s" , cmd )
proc = sps . Popen ( cmd , stderr = sps . STDOUT , stdout = sps . PIPE , close_fds = True )
# # This is long running so we wrap it to make sure we can kill it
try :
    res = proc . communicate ( ) [ 0 ]
except KeyboardInterrupt :
    proc . kill ( )
    raise KeyboardInterrupt
# # check for errors
if proc . returncode :
    LOGGER . error ( "error %s: %s" , cmd , res )
    raise IPyradWarningExit ( "cmd {}: {}" . format ( cmd , res ) )


signature:
def muscle_chunker ( data , sample ) :
    """Splits the muscle alignment into chunks . Each chunk is run on a separate
 computing core . Because the largest clusters are at the beginning of the
 clusters file , assigning equal clusters to each file would put all of the
 large cluster , that take longer to align , near the top . So instead we
 randomly distribute the clusters among the files . If assembly method is
 reference then this step is just a placeholder and nothing happens ."""

implementation:
# # log our location for debugging
LOGGER . info ( "inside muscle_chunker" )
# # only chunk up denovo data , refdata has its own chunking method which
# # makes equal size chunks , instead of uneven chunks like in denovo
if data . paramsdict [ "assembly_method" ] != "reference" : # # get the number of clusters
    clustfile = os . path . join ( data . dirs . clusts , sample . name + ".clust.gz" )
    with iter ( gzip . open ( clustfile , 'rb' ) ) as clustio :
        nloci = sum ( 1 for i in clustio if "//" in i ) // 2
        # tclust = clustio . read ( ) . count ( " / / " ) / / 2
        optim = ( nloci // 20 ) + ( nloci % 20 )
        LOGGER . info ( "optim for align chunks: %s" , optim )
    # # write optim clusters to each tmp file
    clustio = gzip . open ( clustfile , 'rb' )
    inclusts = iter ( clustio . read ( ) . strip ( ) . split ( "//\n//\n" ) )
    # # splitting loci so first file is smaller and last file is bigger
    inc = optim // 10
    for idx in range ( 10 ) : # # how big is this chunk ?
        this = optim + ( idx * inc )
        left = nloci - this
        if idx == 9 : # # grab everything left
            grabchunk = list ( itertools . islice ( inclusts , int ( 1e9 ) ) )
        else : # # grab next chunks - worth of data
            grabchunk = list ( itertools . islice ( inclusts , this ) )
            nloci = left
        # # write the chunk to file
        tmpfile = os . path . join ( data . tmpdir , sample . name + "_chunk_{}.ali" . format ( idx ) )
        with open ( tmpfile , 'wb' ) as out :
            out . write ( "//\n//\n" . join ( grabchunk ) )
    # # write the chunk to file
    # grabchunk = list ( itertools . islice ( inclusts , left ) )
    # if grabchunk :
    # tmpfile = os . path . join ( data . tmpdir , sample . name + " _ chunk _ 9 . ali " )
    # with open ( tmpfile , ' a ' ) as out :
    # out . write ( " \ n / / \ n / / \ n " . join ( grabchunk ) )
    clustio . close ( )


signature:
def reconcat ( data , sample ) :
    """takes aligned chunks ( usually 10 ) and concatenates them"""

implementation:
try : # # get chunks
    chunks = glob . glob ( os . path . join ( data . tmpdir , sample . name + "_chunk_[0-9].aligned" ) )
    # # sort by chunk number , cuts off last 8 = ( aligned )
    chunks . sort ( key = lambda x : int ( x . rsplit ( "_" , 1 ) [ - 1 ] [ : - 8 ] ) )
    LOGGER . info ( "chunk %s" , chunks )
    # # concatenate finished reads
    sample . files . clusters = os . path . join ( data . dirs . clusts , sample . name + ".clustS.gz" )
    # # reconcats aligned clusters
    with gzip . open ( sample . files . clusters , 'wb' ) as out :
        for fname in chunks :
            with open ( fname ) as infile :
                dat = infile . read ( )
                # # avoids mess if last chunk was empty
                if dat . endswith ( "\n" ) :
                    out . write ( dat + "//\n//\n" )
                else :
                    out . write ( dat + "\n//\n//\n" )
            os . remove ( fname )
except Exception as inst :
    LOGGER . error ( "Error in reconcat {}" . format ( inst ) )
    raise


signature:
def derep_concat_split ( data , sample , nthreads , force ) :
    """Running on remote Engine . Refmaps , then merges , then dereplicates ,
 then denovo clusters reads ."""

implementation:
# # report location for debugging
LOGGER . info ( "INSIDE derep %s" , sample . name )
# # MERGED ASSEMBIES ONLY :
# # concatenate edits files within Samples . Returns a new sample . files . edits
# # with the concat file . No change if not merged Assembly .
mergefile = os . path . join ( data . dirs . edits , sample . name + "_merged_.fastq" )
if not force :
    if not os . path . exists ( mergefile ) :
        sample . files . edits = concat_multiple_edits ( data , sample )
    else :
        LOGGER . info ( "skipped concat_multiple_edits: {} exists" . format ( mergefile ) )
else :
    sample . files . edits = concat_multiple_edits ( data , sample )
# # PAIRED DATA ONLY :
# # Denovo : merge or concat fastq pairs [ sample . files . pairs ]
# # Reference : only concat fastq pairs [ ]
# # Denovo + Reference : . . .
if 'pair' in data . paramsdict [ 'datatype' ] : # # the output file handle for merged reads
# # modify behavior of merging vs concating if reference
    if "reference" in data . paramsdict [ "assembly_method" ] :
        nmerged = merge_pairs ( data , sample . files . edits , mergefile , 0 , 0 )
    else :
        nmerged = merge_pairs ( data , sample . files . edits , mergefile , 1 , 1 )
    # # store results
    sample . files . edits = [ ( mergefile , ) ]
    sample . stats . reads_merged = nmerged
# # 3rad uses random adapters to identify pcr duplicates . We will
# # remove pcr dupes here . Basically append the radom adapter to
# # each sequence , do a regular old vsearch derep , then trim
# # off the adapter , and push it down the pipeline . This will
# # remove all identical seqs with identical random i5 adapters .
if "3rad" in data . paramsdict [ "datatype" ] :
    declone_3rad ( data , sample )
    derep_and_sort ( data , os . path . join ( data . dirs . edits , sample . name + "_declone.fastq" ) , os . path . join ( data . dirs . edits , sample . name + "_derep.fastq" ) , nthreads )
else : # # convert fastq to fasta , then derep and sort reads by their size .
# # we pass in only one file b / c paired should be merged by now .
    derep_and_sort ( data , sample . files . edits [ 0 ] [ 0 ] , os . path . join ( data . dirs . edits , sample . name + "_derep.fastq" ) , nthreads )


signature:
def run ( data , samples , noreverse , maxindels , force , ipyclient ) :
    """run the major functions for clustering within samples"""

implementation:
# # list of samples to submit to queue
subsamples = [ ]
# # if sample is already done skip
for sample in samples : # # If sample not in state 2 don ' t try to cluster it .
    if sample . stats . state < 2 :
        print ( """\
    Sample not ready for clustering. First run step2 on sample: {}""" . format ( sample . name ) )
        continue
    if not force :
        if sample . stats . state >= 3 :
            print ( """\
    Skipping {}; aleady clustered. Use force to re-cluster""" . format ( sample . name ) )
        else :
            if sample . stats . reads_passed_filter :
                subsamples . append ( sample )
    else : # # force to overwrite
        if sample . stats . reads_passed_filter :
            subsamples . append ( sample )
# # run subsamples
if not subsamples :
    print ( "  No Samples ready to be clustered. First run step2()." )
else : # # arguments to apply _ jobs , inst catches exceptions
    try : # # make dirs that are needed including tmpdir
        setup_dirs ( data )
        # # if refmapping make filehandles that will be persistent
        if not data . paramsdict [ "assembly_method" ] == "denovo" :
            for sample in subsamples :
                refmap_init ( data , sample , force )
                # # set thread - count to 2 for paired - data
                nthreads = 2
        # # set thread - count to 1 for single - end data
        else :
            nthreads = 1
        # # overwrite nthreads if value in _ ipcluster dict
        if "threads" in data . _ipcluster . keys ( ) :
            nthreads = int ( data . _ipcluster [ "threads" ] )
            # # if more CPUs than there are samples then increase threads
            _ncpus = len ( ipyclient )
            if _ncpus > 2 * len ( data . samples ) :
                nthreads *= 2
        # # submit jobs to be run on cluster
        args = [ data , subsamples , ipyclient , nthreads , maxindels , force ]
        new_apply_jobs ( * args )
    finally : # # this can fail if jobs were not stopped properly and are still
    # # writing to tmpdir . don ' t cleanup if debug is on .
        try :
            log_level = logging . getLevelName ( LOGGER . getEffectiveLevel ( ) )
            if not log_level == "DEBUG" :
                if os . path . exists ( data . tmpdir ) :
                    shutil . rmtree ( data . tmpdir )
                # # get all refmap _ derep . fastqs
                rdereps = glob . glob ( os . path . join ( data . dirs . edits , "*-refmap_derep.fastq" ) )
                # # Remove the unmapped fastq files
                for rmfile in rdereps :
                    os . remove ( rmfile )
        except Exception as _ :
            LOGGER . warning ( "failed to cleanup files/dirs" )


signature:
def parse_params ( args ) :
    """Parse the params file args , create and return Assembly object ."""

implementation:
# # check that params . txt file is correctly formatted .
try :
    with open ( args . params ) as paramsin :
        plines = paramsin . readlines ( )
except IOError as _ :
    sys . exit ( "  No params file found" )
# # check header : big version changes can be distinguished by the header
legacy_version = 0
try : # # try to update the Assembly . . .
    legacy_version = 1
    if not len ( plines [ 0 ] . split ( ) [ 0 ] ) == 7 :
        raise IPyradWarningExit ( """
        Error: file '{}' is not compatible with ipyrad v.{}.
        Please create and update a new params file using the -n argument.
        For info on which parameters have changed see the changelog:
        (http://ipyrad.readthedocs.io/releasenotes.html)
        """ . format ( args . params , ip . __version__ ) )
except IndexError :
    raise IPyradWarningExit ( """
        Error: Params file should not have any empty lines at the top
        of the file. Verify there are no blank lines and rerun ipyrad.
        Offending file - {}
        """ . format ( args . params ) )
# # update and backup
if legacy_version : # which version . . .
# update _ to _ 6 ( )
    pass
# # make into a dict . Ignore blank lines at the end of file
# # Really this will ignore all blank lines
items = [ i . split ( "##" ) [ 0 ] . strip ( ) for i in plines [ 1 : ] if not i . strip ( ) == "" ]
# keys = [ i . split ( " ] " ) [ - 2 ] [ - 1 ] for i in plines [ 1 : ] ]
# keys = range ( len ( plines ) - 1)
keys = ip . Assembly ( 'null' , quiet = True ) . paramsdict . keys ( )
parsedict = { str ( i ) : j for i , j in zip ( keys , items ) }
return parsedict


signature:
def showstats ( parsedict ) :
    """loads assembly or dies , and print stats to screen"""

implementation:
# project _ dir = parsedict [ ' 1 ' ]
project_dir = parsedict [ "project_dir" ]
if not project_dir :
    project_dir = "./"
# # Be nice if somebody also puts in the file extension
# assembly _ name = parsedict [ ' 0 ' ]
assembly_name = parsedict [ "assembly_name" ]
my_assembly = os . path . join ( project_dir , assembly_name )
# # If the project _ dir doesn ' t exist don ' t even bother trying harder .
if not os . path . isdir ( project_dir ) :
    msg = """
    Trying to print stats for Assembly ({}) that doesn't exist. You must
    first run steps before you can show results.
    """ . format ( project_dir )
    sys . exit ( msg )
if not assembly_name :
    msg = """
    Assembly name is not set in params.txt, meaning it was either changed or
    erased since the Assembly was started. Please restore the original name.
    You can find the name of your Assembly in the "project dir": {}.
    """ . format ( project_dir )
    raise IPyradError ( msg )
data = ip . load_json ( my_assembly , quiet = True , cli = True )
print ( "\nSummary stats of Assembly {}" . format ( data . name ) + "\n------------------------------------------------" )
if not data . stats . empty :
    print ( data . stats )
    print ( "\n\nFull stats files" + "\n------------------------------------------------" )
    fullcurdir = os . path . realpath ( os . path . curdir )
    for i in range ( 1 , 8 ) : # enumerate ( sorted ( data . stats _ files ) ) :
        key = "s" + str ( i )
        try :
            val = data . stats_files [ key ]
            val = val . replace ( fullcurdir , "." )
            print ( "step {}: {}" . format ( i , val ) )
        except ( KeyError , AttributeError ) :
            print ( "step {}: None" . format ( i ) )
    print ( "\n" )
else :
    print ( "No stats to display" )


signature:
def branch_assembly ( args , parsedict ) :
    """Load the passed in assembly and create a branch . Copy it
 to a new assembly , and also write out the appropriate params . txt"""

implementation:
# # Get the current assembly
data = getassembly ( args , parsedict )
# # get arguments to branch command
bargs = args . branch
# # get new name , trim off . txt if it was accidentally added
newname = bargs [ 0 ]
if newname . endswith ( ".txt" ) :
    newname = newname [ : - 4 ]
# # look for subsamples
if len ( bargs ) > 1 : # # Branching and subsampling at step 6 is a bad idea , it messes up
# # indexing into the hdf5 cluster file . Warn against this .
    if any ( [ x . stats . state == 6 for x in data . samples . values ( ) ] ) :
        pass
        # # TODODODODODO
        # print ( " wat " )
    # # are we removing or keeping listed samples ?
    subsamples = bargs [ 1 : ]
    # # drop the matching samples
    if bargs [ 1 ] == "-" : # # check drop names
        fails = [ i for i in subsamples [ 1 : ] if i not in data . samples . keys ( ) ]
        if any ( fails ) :
            raise IPyradWarningExit ( "\
                    \n  Failed: unrecognized names requested, check spelling:\n  {}" . format ( "\n  " . join ( [ i for i in fails ] ) ) )
        print ( "  dropping {} samples" . format ( len ( subsamples ) - 1 ) )
        subsamples = list ( set ( data . samples . keys ( ) ) - set ( subsamples ) )
    # # If the arg after the new param name is a file that exists
    if os . path . exists ( bargs [ 1 ] ) :
        new_data = data . branch ( newname , infile = bargs [ 1 ] )
    else :
        new_data = data . branch ( newname , subsamples )
# # keeping all samples
else :
    new_data = data . branch ( newname , None )
print ( "  creating a new branch called '{}' with {} Samples" . format ( new_data . name , len ( new_data . samples ) ) )
print ( "  writing new params file to {}" . format ( "params-" + new_data . name + ".txt\n" ) )
new_data . write_params ( "params-" + new_data . name + ".txt" , force = args . force )


signature:
def merge_assemblies ( args ) :
    """merge all given assemblies into a new assembly . Copies the params
 from the first passed in extant assembly . this function is called
 with the ipyrad - m flag . You must pass it at least 3 values , the first
 is a new assembly name ( a new ` param - newname . txt ` will be created ) .
 The second and third args must be params files for currently existing
 assemblies . Any args beyond the third must also be params file for
 extant assemblies ."""

implementation:
print ( "\n  Merging assemblies: {}" . format ( args . merge [ 1 : ] ) )
# # Make sure there are the right number of args
if len ( args . merge ) < 3 :
    sys . exit ( _WRONG_NUM_CLI_MERGE )
# # Make sure the first arg isn ' t a params file , i could see someone doing it
newname = args . merge [ 0 ]
if os . path . exists ( newname ) and "params-" in newname :
    sys . exit ( _WRONG_ORDER_CLI_MERGE )
# # Make sure first arg will create a param file that doesn ' t already exist
if os . path . exists ( "params-" + newname + ".txt" ) and not args . force :
    sys . exit ( _NAME_EXISTS_MERGE . format ( "params-" + newname + ".txt" ) )
# # Make sure the rest of the args are params files that already exist
assemblies_to_merge = args . merge [ 1 : ]
for assembly in assemblies_to_merge :
    if not os . path . exists ( assembly ) :
        sys . exit ( _DOES_NOT_EXIST_MERGE . format ( assembly ) )
# # Get assemblies for each of the passed in params files .
# # We ' re recycling some of the machinery for loading assemblies here
assemblies = [ ]
for params_file in args . merge [ 1 : ] :
    args . params = params_file
    parsedict = parse_params ( args )
    assemblies . append ( getassembly ( args , parsedict ) )
# # Do the merge
merged_assembly = ip . merge ( newname , assemblies )
# # Write out the merged assembly params file and report success
merged_assembly . write_params ( "params-{}.txt" . format ( newname ) , force = args . force )
print ( "\n  Merging succeeded. New params file for merged assembly:" )
print ( "\n    params-{}.txt\n" . format ( newname ) )


signature:
def getassembly ( args , parsedict ) :
    """loads assembly or creates a new one and set its params from
 parsedict . Does not launch ipcluster ."""

implementation:
# # Creating an assembly with a full path in the name will " work "
# # but it is potentially dangerous , so here we have assembly _ name
# # and assembly _ file , name is used for creating new in cwd , file is
# # used for loading existing .
# # Be nice if the user includes the extension .
# project _ dir = ip . core . assembly . _ expander ( parsedict [ ' 1 ' ] )
# assembly _ name = parsedict [ ' 0 ' ]
project_dir = ip . core . assembly . _expander ( parsedict [ 'project_dir' ] )
assembly_name = parsedict [ 'assembly_name' ]
assembly_file = os . path . join ( project_dir , assembly_name )
# # Assembly creation will handle error checking on
# # the format of the assembly _ name
# # make sure the working directory exists .
if not os . path . exists ( project_dir ) :
    os . mkdir ( project_dir )
try : # # If 1 and force then go ahead and create a new assembly
    if ( '1' in args . steps ) and args . force :
        data = ip . Assembly ( assembly_name , cli = True )
    else :
        data = ip . load_json ( assembly_file , cli = True )
        data . _cli = True
except IPyradWarningExit as _ : # # if no assembly is found then go ahead and make one
    if '1' not in args . steps :
        raise IPyradWarningExit ( "  Error: You must first run step 1 on the assembly: {}" . format ( assembly_file ) )
    else : # # create a new assembly object
        data = ip . Assembly ( assembly_name , cli = True )
# # for entering some params . . .
for param in parsedict : # # trap assignment of assembly _ name since it is immutable .
    if param == "assembly_name" : # # Raise error if user tried to change assembly name
        if parsedict [ param ] != data . name :
            data . set_params ( param , parsedict [ param ] )
    else : # # all other params should be handled by set _ params
        try :
            data . set_params ( param , parsedict [ param ] )
        except IndexError as _ :
            print ( "  Malformed params file: {}" . format ( args . params ) )
            print ( "  Bad parameter {} - {}" . format ( param , parsedict [ param ] ) )
            sys . exit ( - 1 )
return data


signature:
def _check_version ( ) :
    """Test if there ' s a newer version and nag the user to upgrade ."""

implementation:
import urllib2
from distutils . version import LooseVersion
header = "\n -------------------------------------------------------------" + "\n  ipyrad [v.{}]" . format ( ip . __version__ ) + "\n  Interactive assembly and analysis of RAD-seq data" + "\n -------------------------------------------------------------"
try :
    htmldat = urllib2 . urlopen ( "https://anaconda.org/ipyrad/ipyrad" ) . readlines ( )
    curversion = next ( ( x for x in htmldat if "subheader" in x ) , None ) . split ( ">" ) [ 1 ] . split ( "<" ) [ 0 ]
    if LooseVersion ( ip . __version__ ) < LooseVersion ( curversion ) :
        msg = """
  A new version of ipyrad is available (v.{}). To upgrade run:

    conda install -c ipyrad ipyrad\n""" . format ( curversion )
        print ( header + "\n" + msg )
    else :
        pass
        # print ( " You are up to date " )
except Exception as inst : # # Always fail silently
    pass


signature:
def main ( ) :
    """main function"""

implementation:
# # turn off traceback for the CLI
ip . __interactive__ = 0
# # Check for a new version on anaconda
_check_version ( )
# # parse params file input ( returns to stdout if - - help or - - version )
args = parse_command_line ( )
# # Turn the debug output written to ipyrad _ log . txt up to 11!
# # Clean up the old one first , it ' s cleaner to do this here than
# # at the end ( exceptions , etc )
if os . path . exists ( ip . __debugflag__ ) :
    os . remove ( ip . __debugflag__ )
if args . debug :
    print ( "\n  ** Enabling debug mode ** " )
    ip . _debug_on ( )
    atexit . register ( ip . _debug_off )
# # create new paramsfile if - n
if args . new : # # Create a tmp assembly , call write _ params to make default params . txt
    try :
        tmpassembly = ip . Assembly ( args . new , quiet = True , cli = True )
        tmpassembly . write_params ( "params-{}.txt" . format ( args . new ) , force = args . force )
    except Exception as inst :
        print ( inst )
        sys . exit ( 2 )
    print ( "\n  New file 'params-{}.txt' created in {}\n" . format ( args . new , os . path . realpath ( os . path . curdir ) ) )
    sys . exit ( 2 )
# # if params then must provide action argument with it
if args . params :
    if not any ( [ args . branch , args . results , args . steps ] ) :
        print ( """
    Must provide action argument along with -p argument for params file.
    e.g., ipyrad -p params-test.txt -r              ## shows results
    e.g., ipyrad -p params-test.txt -s 12           ## runs steps 1 & 2
    e.g., ipyrad -p params-test.txt -b newbranch    ## branch this assembly
    """ )
        sys . exit ( 2 )
if not args . params :
    if any ( [ args . branch , args . results , args . steps ] ) :
        print ( """
    Must provide params file for branching, doing steps, or getting results.
    e.g., ipyrad -p params-test.txt -r              ## shows results
    e.g., ipyrad -p params-test.txt -s 12           ## runs steps 1 & 2
    e.g., ipyrad -p params-test.txt -b newbranch    ## branch this assembly
    """ )
# # if branching , or merging do not allow steps in same command
# # print spacer
if any ( [ args . branch , args . merge ] ) :
    args . steps = ""
    print ( "" )
# # always print the header when doing steps
header = "\n -------------------------------------------------------------" + "\n  ipyrad [v.{}]" . format ( ip . __version__ ) + "\n  Interactive assembly and analysis of RAD-seq data" + "\n -------------------------------------------------------------"
# # Log the current version . End run around the LOGGER
# # so it ' ll always print regardless of log level .
with open ( ip . __debugfile__ , 'a' ) as logfile :
    logfile . write ( header )
    logfile . write ( "\n  Begin run: {}" . format ( time . strftime ( "%Y-%m-%d %H:%M" ) ) )
    logfile . write ( "\n  Using args {}" . format ( vars ( args ) ) )
    logfile . write ( "\n  Platform info: {}" . format ( os . uname ( ) ) )
# # if merging just do the merge and exit
if args . merge :
    print ( header )
    merge_assemblies ( args )
    sys . exit ( 1 )
# # if download data do it and then exit . Runs single core in CLI .
if args . download :
    if len ( args . download ) == 1 :
        downloaddir = "sra-fastqs"
    else :
        downloaddir = args . download [ 1 ]
    sratools_download ( args . download [ 0 ] , workdir = downloaddir , force = args . force )
    sys . exit ( 1 )
# # create new Assembly or load existing Assembly , quit if args . results
elif args . params :
    parsedict = parse_params ( args )
    if args . branch :
        branch_assembly ( args , parsedict )
    elif args . steps : # # print header
        print ( header )
        # # Only blank the log file if we ' re actually going to run a new
        # # assembly . This used to be in _ _ init _ _ , but had the side effect
        # # of occasionally blanking the log file in an undesirable fashion
        # # for instance if you run a long assembly and it crashes and
        # # then you run ` - r ` and it blanks the log , it ' s crazymaking .
        if os . path . exists ( ip . __debugfile__ ) :
            if os . path . getsize ( ip . __debugfile__ ) > 50000000 :
                with open ( ip . __debugfile__ , 'w' ) as clear :
                    clear . write ( "file reset" )
        # # run Assembly steps
        # # launch or load assembly with custom profile / pid
        data = getassembly ( args , parsedict )
        # # set CLI ipcluster terms
        data . _ipcluster [ "threads" ] = args . threads
        # # if ipyclient is running ( and matched profile ) then use that one
        if args . ipcluster :
            ipyclient = ipp . Client ( profile = args . ipcluster )
            data . _ipcluster [ "cores" ] = len ( ipyclient )
        # # if not then we need to register and launch an ipcluster instance
        else : # # set CLI ipcluster terms
            ipyclient = None
            data . _ipcluster [ "cores" ] = args . cores if args . cores else detect_cpus ( )
            data . _ipcluster [ "engines" ] = "Local"
            if args . MPI :
                data . _ipcluster [ "engines" ] = "MPI"
                if not args . cores :
                    raise IPyradWarningExit ( "must provide -c argument with --MPI" )
            # # register to have a cluster - id with " ip - name "
            data = register_ipcluster ( data )
        # # set to print headers
        data . _headers = 1
        # # run assembly steps
        steps = list ( args . steps )
        data . run ( steps = steps , force = args . force , show_cluster = 1 , ipyclient = ipyclient )
    if args . results :
        showstats ( parsedict )


signature:
def get_binom ( base1 , base2 , estE , estH ) :
    """return probability of base call"""

implementation:
prior_homo = ( 1. - estH ) / 2.
prior_hete = estH
# # calculate probs
bsum = base1 + base2
hetprob = scipy . misc . comb ( bsum , base1 ) / ( 2. ** ( bsum ) )
homoa = scipy . stats . binom . pmf ( base2 , bsum , estE )
homob = scipy . stats . binom . pmf ( base1 , bsum , estE )
# # calculate probs
hetprob *= prior_hete
homoa *= prior_homo
homob *= prior_homo
# # final
probabilities = [ homoa , homob , hetprob ]
bestprob = max ( probabilities ) / float ( sum ( probabilities ) )
# # return
if hetprob > homoa :
    return True , bestprob
else :
    return False , bestprob


signature:
def removerepeats ( consens , arrayed ) :
    """Checks for interior Ns in consensus seqs and removes those that are at
 low depth , here defined as less than 1/3 of the average depth . The prop 1/3
 is chosen so that mindepth = 6 requires 2 base calls that are not in [ N , - ] ."""

implementation:
# # default trim no edges
consens = "" . join ( consens ) . replace ( "-" , "N" )
# # split for pairs
try :
    cons1 , cons2 = consens . split ( "nnnn" )
    split = consens . index ( "nnnn" )
    arr1 = arrayed [ : , : split ]
    arr2 = arrayed [ : , split + 4 : ]
except ValueError :
    cons1 = consens
    cons2 = ""
    arr1 = arrayed
# # trim from left and right of cons1
edges = [ None , None ]
lcons = len ( cons1 )
cons1 = cons1 . lstrip ( "N" )
edges [ 0 ] = lcons - len ( cons1 )
# # trim from right if nonzero
lcons = len ( cons1 )
cons1 = cons1 . rstrip ( "N" )
if lcons - len ( cons1 ) :
    edges [ 1 ] = - 1 * ( lcons - len ( cons1 ) )
# # trim same from arrayed
arr1 = arr1 [ : , edges [ 0 ] : edges [ 1 ] ]
# # trim from left and right of cons2 if present
if cons2 : # # trim from left and right of cons1
    edges = [ None , None ]
    lcons = len ( cons2 )
    cons2 = cons2 . lstrip ( "N" )
    edges [ 0 ] = lcons - len ( cons2 )
    # # trim from right if nonzero
    lcons = len ( cons2 )
    cons2 = cons2 . rstrip ( "N" )
    if lcons - len ( cons2 ) :
        edges [ 1 ] = - 1 * ( lcons - len ( cons2 ) )
    # # trim same from arrayed
    arr2 = arr2 [ : , edges [ 0 ] : edges [ 1 ] ]
    # # reconstitute pairs
    consens = cons1 + "nnnn" + cons2
    consens = np . array ( list ( consens ) )
    sep = np . array ( arr1 . shape [ 0 ] * [ list ( "nnnn" ) ] )
    arrayed = np . hstack ( [ arr1 , sep , arr2 ] )
# # if single - end . . .
else :
    consens = np . array ( list ( cons1 ) )
    arrayed = arr1
# # get column counts of Ns and - s
ndepths = np . sum ( arrayed == 'N' , axis = 0 )
idepths = np . sum ( arrayed == '-' , axis = 0 )
# # get proportion of bases that are N - at each site
nons = ( ( ndepths + idepths ) / float ( arrayed . shape [ 0 ] ) ) >= 0.75
# # boolean of whether base was called N
isn = consens == "N"
# # make ridx
ridx = nons * isn
# # apply filter
consens = consens [ ~ ridx ]
arrayed = arrayed [ : , ~ ridx ]
return consens , arrayed