signature:
def newconsensus ( data , sample , tmpchunk , optim ) :
    """new faster replacement to consensus"""

implementation:
# # do reference map funcs ?
isref = "reference" in data . paramsdict [ "assembly_method" ]
# # temporarily store the mean estimates to Assembly
data . _este = data . stats . error_est . mean ( )
data . _esth = data . stats . hetero_est . mean ( )
# # get number relative to tmp file
tmpnum = int ( tmpchunk . split ( "." ) [ - 1 ] )
# # prepare data for reading
clusters = open ( tmpchunk , 'rb' )
pairdealer = itertools . izip ( * [ iter ( clusters ) ] * 2 )
maxlen = data . _hackersonly [ "max_fragment_length" ]
# # write to tmp cons to file to be combined later
consenshandle = os . path . join ( data . dirs . consens , sample . name + "_tmpcons." + str ( tmpnum ) )
tmp5 = consenshandle . replace ( "_tmpcons." , "_tmpcats." )
with h5py . File ( tmp5 , 'w' ) as io5 :
    io5 . create_dataset ( "cats" , ( optim , maxlen , 4 ) , dtype = np . uint32 )
    io5 . create_dataset ( "alls" , ( optim , ) , dtype = np . uint8 )
    io5 . create_dataset ( "chroms" , ( optim , 3 ) , dtype = np . int64 )
    # # local copies to use to fill the arrays
    catarr = io5 [ "cats" ] [ : ]
    nallel = io5 [ "alls" ] [ : ]
    refarr = io5 [ "chroms" ] [ : ]
# # if reference - mapped then parse the fai to get index number of chroms
if isref :
    fai = pd . read_csv ( data . paramsdict [ "reference_sequence" ] + ".fai" , names = [ 'scaffold' , 'size' , 'sumsize' , 'a' , 'b' ] , sep = "\t" )
    faidict = { j : i for i , j in enumerate ( fai . scaffold ) }
# # store data for stats counters
counters = { "name" : tmpnum , "heteros" : 0 , "nsites" : 0 , "nconsens" : 0 }
# # store data for what got filtered
filters = { "depth" : 0 , "maxh" : 0 , "maxn" : 0 }
# # store data for writing
storeseq = { }
# # set max limits
if 'pair' in data . paramsdict [ "datatype" ] :
    maxhet = sum ( data . paramsdict [ "max_Hs_consens" ] )
    maxn = sum ( data . paramsdict [ "max_Ns_consens" ] )
else :
    maxhet = data . paramsdict [ "max_Hs_consens" ] [ 0 ]
    maxn = data . paramsdict [ "max_Ns_consens" ] [ 0 ]
# # load the refmap dictionary if refmapping
done = 0
while not done :
    try :
        done , chunk = clustdealer ( pairdealer , 1 )
    except IndexError :
        raise IPyradError ( "clustfile formatting error in %s" , chunk )
    if chunk : # # get names and seqs
        piece = chunk [ 0 ] . strip ( ) . split ( "\n" )
        names = piece [ 0 : : 2 ]
        seqs = piece [ 1 : : 2 ]
        # # pull replicate read info from seqs
        reps = [ int ( sname . split ( ";" ) [ - 2 ] [ 5 : ] ) for sname in names ]
        # # IF this is a reference mapped read store the chrom and pos info
        # # - 1 defaults to indicating an anonymous locus , since we are using
        # # the faidict as 0 indexed . If chrompos fails it defaults to - 1
        ref_position = ( - 1 , 0 , 0 )
        if isref :
            try : # # parse position from name string
                name , _ , _ = names [ 0 ] . rsplit ( ";" , 2 )
                chrom , pos0 , pos1 = name . rsplit ( ":" , 2 )
                # # pull idx from . fai reference dict
                chromint = faidict [ chrom ] + 1
                ref_position = ( int ( chromint ) , int ( pos0 ) , int ( pos1 ) )
            except Exception as inst :
                LOGGER . debug ( "Reference sequence chrom/pos failed for {}" . format ( names [ 0 ] ) )
                LOGGER . debug ( inst )
        # # apply read depth filter
        if nfilter1 ( data , reps ) : # # get stacks of base counts
            sseqs = [ list ( seq ) for seq in seqs ]
            arrayed = np . concatenate ( [ [ seq ] * rep for seq , rep in zip ( sseqs , reps ) ] )
            arrayed = arrayed [ : , : maxlen ]
            # # get consens call for each site , applies paralog - x - site filter
            # consens = np . apply _ along _ axis ( basecall , 0 , arrayed , data )
            consens = basecaller ( arrayed , data . paramsdict [ "mindepth_majrule" ] , data . paramsdict [ "mindepth_statistical" ] , data . _esth , data . _este , )
            # # apply a filter to remove low coverage sites / Ns that
            # # are likely sequence repeat errors . This is only applied to
            # # clusters that already passed the read - depth filter ( 1)
            if "N" in consens :
                try :
                    consens , arrayed = removerepeats ( consens , arrayed )
                except ValueError as _ :
                    LOGGER . info ( "Caught a bad chunk w/ all Ns. Skip it." )
                    continue
            # # get hetero sites
            hidx = [ i for ( i , j ) in enumerate ( consens ) if j in list ( "RKSYWM" ) ]
            nheteros = len ( hidx )
            # # filter for max number of hetero sites
            if nfilter2 ( nheteros , maxhet ) : # # filter for maxN , & minlen
                if nfilter3 ( consens , maxn ) : # # counter right now
                    current = counters [ "nconsens" ]
                    # # get N alleles and get lower case in consens
                    consens , nhaps = nfilter4 ( consens , hidx , arrayed )
                    # # store the number of alleles observed
                    nallel [ current ] = nhaps
                    # # store a reduced array with only CATG
                    catg = np . array ( [ np . sum ( arrayed == i , axis = 0 ) for i in list ( "CATG" ) ] , dtype = 'uint32' ) . T
                    catarr [ current , : catg . shape [ 0 ] , : ] = catg
                    refarr [ current ] = ref_position
                    # # store the seqdata for tmpchunk
                    storeseq [ counters [ "name" ] ] = "" . join ( list ( consens ) )
                    counters [ "name" ] += 1
                    counters [ "nconsens" ] += 1
                    counters [ "heteros" ] += nheteros
                else : # LOGGER . debug ( " @ haplo " )
                    filters [ 'maxn' ] += 1
            else : # LOGGER . debug ( " @ hetero " )
                filters [ 'maxh' ] += 1
        else : # LOGGER . debug ( " @ depth " )
            filters [ 'depth' ] += 1
# # close infile io
clusters . close ( )
# # write final consens string chunk
if storeseq :
    with open ( consenshandle , 'wb' ) as outfile :
        outfile . write ( "\n" . join ( [ ">" + sample . name + "_" + str ( key ) + "\n" + str ( storeseq [ key ] ) for key in storeseq ] ) )
# # write to h5 array , this can be a bit slow on big data sets and is not
# # currently convered by progressbar movement .
with h5py . File ( tmp5 , 'a' ) as io5 :
    io5 [ "cats" ] [ : ] = catarr
    io5 [ "alls" ] [ : ] = nallel
    io5 [ "chroms" ] [ : ] = refarr
del catarr
del nallel
del refarr
# # return stats
counters [ 'nsites' ] = sum ( [ len ( i ) for i in storeseq . itervalues ( ) ] )
return counters , filters


signature:
def basecaller ( arrayed , mindepth_majrule , mindepth_statistical , estH , estE ) :
    """call all sites in a locus array ."""

implementation:
# # an array to fill with consensus site calls
cons = np . zeros ( arrayed . shape [ 1 ] , dtype = np . uint8 )
cons . fill ( 78 )
arr = arrayed . view ( np . uint8 )
# # iterate over columns
for col in xrange ( arr . shape [ 1 ] ) : # # the site of focus
    carr = arr [ : , col ]
    # # make mask of N and - sites
    mask = carr == 45
    mask += carr == 78
    marr = carr [ ~ mask ]
    # # skip if only empties ( e . g . , N - )
    if not marr . shape [ 0 ] :
        cons [ col ] = 78
    # # skip if not variable
    elif np . all ( marr == marr [ 0 ] ) :
        cons [ col ] = marr [ 0 ]
    # # estimate variable site call
    else : # # get allele freqs ( first - most , second , third = p , q , r )
        counts = np . bincount ( marr )
        pbase = np . argmax ( counts )
        nump = counts [ pbase ]
        counts [ pbase ] = 0
        qbase = np . argmax ( counts )
        numq = counts [ qbase ]
        counts [ qbase ] = 0
        rbase = np . argmax ( counts )
        numr = counts [ rbase ]
        # # based on biallelic depth
        bidepth = nump + numq
        if bidepth < mindepth_majrule :
            cons [ col ] = 78
        else : # # if depth is too high , reduce to sampled int
            if bidepth > 500 :
                base1 = int ( 500 * ( nump / float ( bidepth ) ) )
                base2 = int ( 500 * ( numq / float ( bidepth ) ) )
            else :
                base1 = nump
                base2 = numq
            # # make statistical base call
            if bidepth >= mindepth_statistical :
                ishet , prob = get_binom ( base1 , base2 , estE , estH )
                # LOGGER . info ( " ishet , prob , b1 , b2 : % s % s % s % s " , ishet , prob , base1 , base2)
                if prob < 0.95 :
                    cons [ col ] = 78
                else :
                    if ishet :
                        cons [ col ] = TRANS [ ( pbase , qbase ) ]
                    else :
                        cons [ col ] = pbase
            # # make majrule base call
            else : # if bidepth > = mindepth _ majrule :
                if nump == numq :
                    cons [ col ] = TRANS [ ( pbase , qbase ) ]
                else :
                    cons [ col ] = pbase
return cons . view ( "S1" )


signature:
def nfilter1 ( data , reps ) :
    """applies read depths filter"""

implementation:
if sum ( reps ) >= data . paramsdict [ "mindepth_majrule" ] and sum ( reps ) <= data . paramsdict [ "maxdepth" ] :
    return 1
else :
    return 0


signature:
def nfilter4 ( consens , hidx , arrayed ) :
    """applies max haplotypes filter returns pass and consens"""

implementation:
# # if less than two Hs then there is only one allele
if len ( hidx ) < 2 :
    return consens , 1
# # store base calls for hetero sites
harray = arrayed [ : , hidx ]
# # remove any reads that have N or - base calls at hetero sites
# # these cannot be used when calling alleles currently .
harray = harray [ ~ np . any ( harray == "-" , axis = 1 ) ]
harray = harray [ ~ np . any ( harray == "N" , axis = 1 ) ]
# # get counts of each allele ( e . g . , AT : 2 , CG : 2)
ccx = Counter ( [ tuple ( i ) for i in harray ] )
# # Two possibilities we would like to distinguish , but we can ' t . Therefore ,
# # we just throw away low depth third alleles that are within seq . error .
# # 1 ) a third base came up as a sequencing error but is not a unique allele
# # 2 ) a third or more unique allele is there but at low frequency
# # remove low freq alleles if more than 2 , since they may reflect
# # sequencing errors at hetero sites , making a third allele , or a new
# # allelic combination that is not real .
if len ( ccx ) > 2 :
    totdepth = harray . shape [ 0 ]
    cutoff = max ( 1 , totdepth // 10 )
    alleles = [ i for i in ccx if ccx [ i ] > cutoff ]
else :
    alleles = ccx . keys ( )
# # how many high depth alleles ?
nalleles = len ( alleles )
# # if 2 alleles then save the phase using lowercase coding
if nalleles == 2 :
    try :
        consens = storealleles ( consens , hidx , alleles )
    except ( IndexError , KeyError ) : # # the H sites do not form good alleles
        LOGGER . info ( "failed at phasing loc, skipping" )
        LOGGER . info ( """
    consens %s
    hidx %s
    alleles %s
                """ , consens , hidx , alleles )
    return consens , nalleles
# # just return the info for later filtering
else :
    return consens , nalleles


signature:
def storealleles ( consens , hidx , alleles ) :
    """store phased allele data for diploids"""

implementation:
# # find the first hetero site and choose the priority base
# # example , if W : then priority base in A and not T . PRIORITY = ( order : CATG )
bigbase = PRIORITY [ consens [ hidx [ 0 ] ] ]
# # find which allele has priority based on bigbase
bigallele = [ i for i in alleles if i [ 0 ] == bigbase ] [ 0 ]
# # uplow other bases relative to this one and the priority list
# # e . g . , if there are two hetero sites ( WY ) and the two alleles are
# # AT and TC , then since bigbase of ( W ) is A second hetero site should
# # be stored as y , since the ordering is swapped in this case ; the priority
# # base ( C versus T ) is C , but C goes with the minor base at h site 1.
# consens = list ( consens )
for hsite , pbase in zip ( hidx [ 1 : ] , bigallele [ 1 : ] ) :
    if PRIORITY [ consens [ hsite ] ] != pbase :
        consens [ hsite ] = consens [ hsite ] . lower ( )
# # return consens
return consens


signature:
def cleanup ( data , sample , statsdicts ) :
    """cleaning up . optim is the size ( nloci ) of tmp arrays"""

implementation:
LOGGER . info ( "in cleanup for: %s" , sample . name )
isref = 'reference' in data . paramsdict [ "assembly_method" ]
# # collect consens chunk files
combs1 = glob . glob ( os . path . join ( data . dirs . consens , sample . name + "_tmpcons.*" ) )
combs1 . sort ( key = lambda x : int ( x . split ( "." ) [ - 1 ] ) )
# # collect tmpcat files
tmpcats = glob . glob ( os . path . join ( data . dirs . consens , sample . name + "_tmpcats.*" ) )
tmpcats . sort ( key = lambda x : int ( x . split ( "." ) [ - 1 ] ) )
# # get shape info from the first cat , ( optim , maxlen , 4)
with h5py . File ( tmpcats [ 0 ] , 'r' ) as io5 :
    optim , maxlen , _ = io5 [ 'cats' ] . shape
# # save as a chunked compressed hdf5 array
handle1 = os . path . join ( data . dirs . consens , sample . name + ".catg" )
with h5py . File ( handle1 , 'w' ) as ioh5 :
    nloci = len ( tmpcats ) * optim
    dcat = ioh5 . create_dataset ( "catg" , ( nloci , maxlen , 4 ) , dtype = np . uint32 , chunks = ( optim , maxlen , 4 ) , compression = "gzip" )
    dall = ioh5 . create_dataset ( "nalleles" , ( nloci , ) , dtype = np . uint8 , chunks = ( optim , ) , compression = "gzip" )
    # # only create chrom for reference - aligned data
    if isref :
        dchrom = ioh5 . create_dataset ( "chroms" , ( nloci , 3 ) , dtype = np . int64 , chunks = ( optim , 3 ) , compression = "gzip" )
    # # Combine all those tmp cats into the big cat
    start = 0
    for icat in tmpcats :
        io5 = h5py . File ( icat , 'r' )
        end = start + optim
        dcat [ start : end ] = io5 [ 'cats' ] [ : ]
        dall [ start : end ] = io5 [ 'alls' ] [ : ]
        if isref :
            dchrom [ start : end ] = io5 [ 'chroms' ] [ : ]
        start += optim
        io5 . close ( )
        os . remove ( icat )
# # store the handle to the Sample
sample . files . database = handle1
# # record results
xcounters = { "nconsens" : 0 , "heteros" : 0 , "nsites" : 0 }
xfilters = { "depth" : 0 , "maxh" : 0 , "maxn" : 0 }
# # merge finished consens stats
for counters , filters in statsdicts : # # sum individual counters
    for key in xcounters :
        xcounters [ key ] += counters [ key ]
    for key in xfilters :
        xfilters [ key ] += filters [ key ]
# # merge consens read files
handle1 = os . path . join ( data . dirs . consens , sample . name + ".consens.gz" )
with gzip . open ( handle1 , 'wb' ) as out :
    for fname in combs1 :
        with open ( fname ) as infile :
            out . write ( infile . read ( ) + "\n" )
        os . remove ( fname )
sample . files . consens = [ handle1 ]
# # set Sample stats _ dfs values
if int ( xcounters [ 'nsites' ] ) :
    prop = int ( xcounters [ "heteros" ] ) / float ( xcounters [ 'nsites' ] )
else :
    prop = 0
sample . stats_dfs . s5 . nsites = int ( xcounters [ "nsites" ] )
sample . stats_dfs . s5 . nhetero = int ( xcounters [ "heteros" ] )
sample . stats_dfs . s5 . filtered_by_depth = xfilters [ 'depth' ]
sample . stats_dfs . s5 . filtered_by_maxH = xfilters [ 'maxh' ]
sample . stats_dfs . s5 . filtered_by_maxN = xfilters [ 'maxn' ]
sample . stats_dfs . s5 . reads_consens = int ( xcounters [ "nconsens" ] )
sample . stats_dfs . s5 . clusters_total = sample . stats_dfs . s3 . clusters_total
sample . stats_dfs . s5 . heterozygosity = float ( prop )
# # set the Sample stats summary value
sample . stats . reads_consens = int ( xcounters [ "nconsens" ] )
# # save state to Sample if successful
if sample . stats . reads_consens :
    sample . stats . state = 5
else :
    print ( "No clusters passed filtering in Sample: {}" . format ( sample . name ) )
return sample


signature:
def chunk_clusters ( data , sample ) :
    """split job into bits and pass to the client"""

implementation:
# # counter for split job submission
num = 0
# # set optim size for chunks in N clusters . The first few chunks take longer
# # because they contain larger clusters , so we create 4X as many chunks as
# # processors so that they are split more evenly .
optim = int ( ( sample . stats . clusters_total // data . cpus ) + ( sample . stats . clusters_total % data . cpus ) )
# # break up the file into smaller tmp files for each engine
# # chunking by cluster is a bit trickier than chunking by N lines
chunkslist = [ ]
# # open to clusters
with gzip . open ( sample . files . clusters , 'rb' ) as clusters : # # create iterator to sample 2 lines at a time
    pairdealer = itertools . izip ( * [ iter ( clusters ) ] * 2 )
    # # Use iterator to sample til end of cluster
    done = 0
    while not done : # # grab optim clusters and write to file .
        done , chunk = clustdealer ( pairdealer , optim )
        chunkhandle = os . path . join ( data . dirs . clusts , "tmp_" + str ( sample . name ) + "." + str ( num * optim ) )
        if chunk :
            chunkslist . append ( ( optim , chunkhandle ) )
            with open ( chunkhandle , 'wb' ) as outchunk :
                outchunk . write ( "//\n//\n" . join ( chunk ) + "//\n//\n" )
            num += 1
return chunkslist


signature:
def get_subsamples ( data , samples , force ) :
    """Apply state , ncluster , and force filters to select samples to be run ."""

implementation:
subsamples = [ ]
for sample in samples :
    if not force :
        if sample . stats . state >= 5 :
            print ( """\
    Skipping Sample {}; Already has consens reads. Use force arg to overwrite.\
    """ . format ( sample . name ) )
        elif not sample . stats . clusters_hidepth :
            print ( """\
    Skipping Sample {}; No clusters found.""" . format ( sample . name , int ( sample . stats . clusters_hidepth ) ) )
        elif sample . stats . state < 4 :
            print ( """\
    Skipping Sample {}; not yet finished step4 """ . format ( sample . name ) )
        else :
            subsamples . append ( sample )
    else :
        if not sample . stats . clusters_hidepth :
            print ( """\
    Skipping Sample {}; No clusters found in {}.""" . format ( sample . name , sample . files . clusters ) )
        elif sample . stats . state < 4 :
            print ( """\
    Skipping Sample {}; not yet finished step4""" . format ( sample . name ) )
        else :
            subsamples . append ( sample )
if len ( subsamples ) == 0 :
    raise IPyradWarningExit ( """
    No samples to cluster, exiting.
    """ )
# # if sample is already done skip
if "hetero_est" not in data . stats :
    print ( "  No estimates of heterozygosity and error rate. Using default " "values" )
    for sample in subsamples :
        sample . stats . hetero_est = 0.001
        sample . stats . error_est = 0.0001
if data . _headers :
    print ( u"""\
  Mean error  [{:.5f} sd={:.5f}]
  Mean hetero [{:.5f} sd={:.5f}]""" . format ( data . stats . error_est . mean ( ) , data . stats . error_est . std ( ) , data . stats . hetero_est . mean ( ) , data . stats . hetero_est . std ( ) ) )
return subsamples


signature:
def run ( data , samples , force , ipyclient ) :
    """checks if the sample should be run and passes the args"""

implementation:
# # prepare dirs
data . dirs . consens = os . path . join ( data . dirs . project , data . name + "_consens" )
if not os . path . exists ( data . dirs . consens ) :
    os . mkdir ( data . dirs . consens )
# # zap any tmp files that might be leftover
tmpcons = glob . glob ( os . path . join ( data . dirs . consens , "*_tmpcons.*" ) )
tmpcats = glob . glob ( os . path . join ( data . dirs . consens , "*_tmpcats.*" ) )
for tmpfile in tmpcons + tmpcats :
    os . remove ( tmpfile )
# # filter through samples for those ready
samples = get_subsamples ( data , samples , force )
# # set up parallel client : how many cores ?
lbview = ipyclient . load_balanced_view ( )
data . cpus = data . _ipcluster [ "cores" ]
if not data . cpus :
    data . cpus = len ( ipyclient . ids )
# # wrap everything to ensure destruction of temp files
inst = ""
try : # # calculate depths , if they changed .
    samples = calculate_depths ( data , samples , lbview )
    # # chunk clusters into bits for parallel processing
    lasyncs = make_chunks ( data , samples , lbview )
    # # process chunks and cleanup
    process_chunks ( data , samples , lasyncs , lbview )
except KeyboardInterrupt as inst :
    raise inst
finally : # # if process failed at any point delete tmp files
    tmpcons = glob . glob ( os . path . join ( data . dirs . clusts , "tmp_*.[0-9]*" ) )
    tmpcons += glob . glob ( os . path . join ( data . dirs . consens , "*_tmpcons.*" ) )
    tmpcons += glob . glob ( os . path . join ( data . dirs . consens , "*_tmpcats.*" ) )
    for tmpchunk in tmpcons :
        os . remove ( tmpchunk )
    # # Finished step 5 . Set step 6 checkpoint to 0 to force
    # # re - running from scratch .
    data . _checkpoint = 0


signature:
def calculate_depths ( data , samples , lbview ) :
    """check whether mindepth has changed , and thus whether clusters _ hidepth
 needs to be recalculated , and get new maxlen for new highdepth clusts .
 if mindepth not changed then nothing changes ."""

implementation:
# # send jobs to be processed on engines
start = time . time ( )
printstr = " calculating depths    | {} | s5 |"
recaljobs = { }
maxlens = [ ]
for sample in samples :
    recaljobs [ sample . name ] = lbview . apply ( recal_hidepth , * ( data , sample ) )
# # block until finished
while 1 :
    ready = [ i . ready ( ) for i in recaljobs . values ( ) ]
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( len ( ready ) , sum ( ready ) , printstr . format ( elapsed ) , spacer = data . _spacer )
    time . sleep ( 0.1 )
    if len ( ready ) == sum ( ready ) :
        print ( "" )
        break
# # check for failures and collect results
modsamples = [ ]
for sample in samples :
    if not recaljobs [ sample . name ] . successful ( ) :
        LOGGER . error ( "  sample %s failed: %s" , sample . name , recaljobs [ sample . name ] . exception ( ) )
    else :
        modsample , _ , maxlen , _ , _ = recaljobs [ sample . name ] . result ( )
        modsamples . append ( modsample )
        maxlens . append ( maxlen )
# # reset global maxlen if something changed
data . _hackersonly [ "max_fragment_length" ] = int ( max ( maxlens ) ) + 4
return samples


signature:
def make_chunks ( data , samples , lbview ) :
    """calls chunk _ clusters and tracks progress ."""

implementation:
# # first progress bar
start = time . time ( )
printstr = " chunking clusters     | {} | s5 |"
elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
progressbar ( 10 , 0 , printstr . format ( elapsed ) , spacer = data . _spacer )
# # send off samples to be chunked
lasyncs = { }
for sample in samples :
    lasyncs [ sample . name ] = lbview . apply ( chunk_clusters , * ( data , sample ) )
# # block until finished
while 1 :
    ready = [ i . ready ( ) for i in lasyncs . values ( ) ]
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( len ( ready ) , sum ( ready ) , printstr . format ( elapsed ) , spacer = data . _spacer )
    time . sleep ( 0.1 )
    if len ( ready ) == sum ( ready ) :
        print ( "" )
        break
# # check for failures
for sample in samples :
    if not lasyncs [ sample . name ] . successful ( ) :
        LOGGER . error ( "  sample %s failed: %s" , sample . name , lasyncs [ sample . name ] . exception ( ) )
return lasyncs


signature:
def process_chunks ( data , samples , lasyncs , lbview ) :
    """submit chunks to consens func and . . ."""

implementation:
# # send chunks to be processed
start = time . time ( )
asyncs = { sample . name : [ ] for sample in samples }
printstr = " consens calling       | {} | s5 |"
# # get chunklist from results
for sample in samples :
    clist = lasyncs [ sample . name ] . result ( )
    for optim , chunkhandle in clist :
        args = ( data , sample , chunkhandle , optim )
        # asyncs [ sample . name ] . append ( lbview . apply _ async ( consensus , * args ) )
        asyncs [ sample . name ] . append ( lbview . apply_async ( newconsensus , * args ) )
        elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
        progressbar ( 10 , 0 , printstr . format ( elapsed ) , spacer = data . _spacer )
# # track progress
allsyncs = list ( itertools . chain ( * [ asyncs [ i . name ] for i in samples ] ) )
while 1 :
    ready = [ i . ready ( ) for i in allsyncs ]
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( len ( ready ) , sum ( ready ) , printstr . format ( elapsed ) , spacer = data . _spacer )
    time . sleep ( 0.1 )
    if len ( ready ) == sum ( ready ) :
        break
# # get clean samples
casyncs = { }
for sample in samples :
    rlist = asyncs [ sample . name ]
    statsdicts = [ i . result ( ) for i in rlist ]
    casyncs [ sample . name ] = lbview . apply ( cleanup , * ( data , sample , statsdicts ) )
while 1 :
    ready = [ i . ready ( ) for i in casyncs . values ( ) ]
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( 10 , 10 , printstr . format ( elapsed ) , spacer = data . _spacer )
    time . sleep ( 0.1 )
    if len ( ready ) == sum ( ready ) :
        print ( "" )
        break
# # check for failures :
for key in asyncs :
    asynclist = asyncs [ key ]
    for async in asynclist :
        if not async . successful ( ) :
            LOGGER . error ( "  async error: %s \n%s" , key , async . exception ( ) )
for key in casyncs :
    if not casyncs [ key ] . successful ( ) :
        LOGGER . error ( "  casync error: %s \n%s" , key , casyncs [ key ] . exception ( ) )
# # get samples back
subsamples = [ i . result ( ) for i in casyncs . values ( ) ]
for sample in subsamples :
    data . samples [ sample . name ] = sample
# # build Assembly stats
data . stats_dfs . s5 = data . _build_stat ( "s5" )
# # write stats file
data . stats_files . s5 = os . path . join ( data . dirs . consens , 's5_consens_stats.txt' )
with io . open ( data . stats_files . s5 , 'w' ) as out : # out . write ( data . stats _ dfs . s5 . to _ string ( ) )
    data . stats_dfs . s5 . to_string ( buf = out , formatters = { 'clusters_total' : '{:.0f}' . format , 'filtered_by_depth' : '{:.0f}' . format , 'filtered_by_maxH' : '{:.0f}' . format , 'filtered_by_maxN' : '{:.0f}' . format , 'reads_consens' : '{:.0f}' . format , 'nsites' : '{:.0f}' . format , 'nhetero' : '{:.0f}' . format , 'heterozygosity' : '{:.5f}' . format } )


signature:
def make ( data , samples ) :
    """reads in . loci and builds alleles from case characters"""

implementation:
# read in loci file
outfile = open ( os . path . join ( data . dirs . outfiles , data . name + ".alleles" ) , 'w' )
lines = open ( os . path . join ( data . dirs . outfiles , data . name + ".loci" ) , 'r' )
# # Get the longest sample name for pretty printing
longname = max ( len ( x ) for x in data . samples . keys ( ) )
# # Padding between name and sequence in output file . This should be the
# # same as write _ outfiles . write _ tmp _ loci . name _ padding
name_padding = 5
writing = [ ]
loc = 0
for line in lines :
    if ">" in line :
        name , seq = line . split ( " " ) [ 0 ] , line . split ( " " ) [ - 1 ]
        allele1 , allele2 = splitalleles ( seq . strip ( ) )
        # # Format the output string . the " - 2 " below accounts for the additional
        # # 2 characters added to the sample name that don ' t get added to the
        # # snpsites line , so you gotta bump this line back 2 to make it
        # # line up right .
        writing . append ( name + "_0" + " " * ( longname - len ( name ) - 2 + name_padding ) + allele1 )
        writing . append ( name + "_1" + " " * ( longname - len ( name ) - 2 + name_padding ) + allele2 )
    else :
        writing . append ( line . strip ( ) )
    loc += 1
    # # print every 10K loci "
    if not loc % 10000 :
        outfile . write ( "\n" . join ( writing ) + "\n" )
        writing = [ ]
outfile . write ( "\n" . join ( writing ) )
outfile . close ( )


signature:
def make ( data , samples ) :
    """builds snps output"""

implementation:
# # get attr
ploidy = data . paramsdict [ "max_alleles_consens" ]
names = [ i . name for i in samples ]
longname = max ( [ len ( i ) for i in names ] )
# # TODO : use iter cuz of super huge files
inloci = open ( os . path . join ( data . dirs . outfiles , data . name + ".loci" ) , 'r' ) . read ( )
# # Potential outfiles
snpsout = os . path . join ( data . dirs . outfiles , data . name + ".snps" )
usnpsout = os . path . join ( data . dirs . outfiles , data . name + ".usnps" )
structout = os . path . join ( data . dirs . outfiles , data . name + ".str" )
genoout = os . path . join ( data . dirs . outfiles , data . name + ".snps.geno" )
ugenoout = os . path . join ( data . dirs . outfiles , data . name + ".usnps.geno" )
# # Output file for writing some stats
statsfile = os . path . join ( data . dirs . outfiles , data . name + ".snps.stats" )
# # The output formats to write
formats = data . paramsdict [ "output_formats" ]
seed = data . _hackersonly [ "random_seed" ]
np . random . seed ( int ( seed ) )
# # output . snps and . unlinked _ snps "
S = { }
# # snp dict
Si = { }
# # unlinked snp dict
for name in list ( names ) :
    S [ name ] = [ ]
    Si [ name ] = [ ]
# # record bi - allelic snps "
bis = 0
# # for each locus select out the SNPs "
for loc in inloci . strip ( ) . split ( "|" ) [ : - 1 ] :
    pis = ""
    ns = [ ]
    ss = [ ]
    cov = { }
    # # record coverage for each SNP
    for line in loc . split ( "\n" ) :
        if ">" in line :
            ns . append ( line . split ( ) [ 0 ] . replace ( ">" , "" ) )
            ss . append ( line . split ( ) [ - 1 ] )
        else :
            pis = [ i [ 0 ] for i in enumerate ( line ) if i [ 1 ] in list ( '*-' ) ]
    # # assign snps to S , and record coverage for usnps "
    for tax in S :
        if tax in ns :
            if pis :
                for snpsite in pis :
                    snpsite -= ( longname + 5 )
                    S [ tax ] . append ( ss [ ns . index ( tax ) ] [ snpsite ] )
                    if snpsite not in cov :
                        cov [ snpsite ] = 1
                    else :
                        cov [ snpsite ] += 1
                    # # downweight selection of gap sites "
                    if ss [ ns . index ( tax ) ] [ snpsite ] != '-' :
                        cov [ snpsite ] += 1
        else :
            if pis :
                for snpsite in pis :
                    S [ tax ] . append ( "N" )
                Si [ tax ] . append ( "N" )
    # # randomly select among snps w / greatest coverage for unlinked snp "
    maxlist = [ ]
    for j , k in cov . items ( ) :
        if k == max ( cov . values ( ) ) :
            maxlist . append ( j )
    # # Is bi - allelic ? "
    bisnps = [ ]
    for i in maxlist :
        if len ( set ( [ ss [ ns . index ( tax ) ] [ i ] for tax in S if tax in ns ] ) ) < 3 :
            bisnps . append ( i )
    # rando = pis [ np . random . randint ( len ( pis ) ) ]
    # rando - = ( longname + 5)
    if bisnps :
        rando = bisnps [ np . random . randint ( len ( bisnps ) ) ]
    elif maxlist :
        rando = maxlist [ np . random . randint ( len ( maxlist ) ) ]
    tbi = 0
    for tax in S :
        if tax in ns :
            if pis : # # if none are bi - allelic "
                if not bisnps :
                    tbi += 1
                Si [ tax ] . append ( ss [ ns . index ( tax ) ] [ rando ] )
        if pis : # # add spacer between loci "
            S [ tax ] . append ( " " )
        else : # # invariable locus "
            S [ tax ] . append ( "_ " )
    bis += tbi
# # names
SF = list ( S . keys ( ) )
SF . sort ( )
# # Write linked snps format
if "snps" in formats :
    with open ( snpsout , 'w' ) as outfile :
        print >> outfile , "## %s taxa, %s loci, %s snps" % ( len ( S ) , len ( "" . join ( S . values ( ) [ 0 ] ) . split ( " " ) ) - 1 , len ( "" . join ( S [ SF [ 0 ] ] ) . replace ( " " , "" ) ) )
        for i in SF :
            print >> outfile , i + ( " " * ( longname - len ( i ) + 3 ) ) + "" . join ( S [ i ] )
# # Write unlinked snps format
if "usnps" in formats :
    with open ( usnpsout , 'w' ) as outfile :
        print >> outfile , len ( Si ) , len ( "" . join ( Si . values ( ) [ 0 ] ) )
        for i in SF :
            print >> outfile , i + ( " " * ( longname - len ( i ) + 3 ) ) + "" . join ( Si [ i ] )
    with open ( statsfile , 'a' ) as statsout :
        print >> statsout , "sampled unlinked SNPs=" , len ( Si . values ( ) [ 0 ] )
        print >> statsout , "sampled unlinked bi-allelic SNPs=" , len ( Si . values ( ) [ 0 ] ) - bis
# # Write STRUCTURE format
if "str" in formats :
    with open ( structout , 'w' ) as outfile :
        B = { 'A' : '0' , 'T' : '1' , 'G' : '2' , 'C' : '3' , 'N' : '-9' , '-' : '-9' }
        if ploidy > 1 :
            for line in SF :
                print >> outfile , line + ( " " * ( longname - len ( line ) + 3 ) ) + "\t" * 6 + "\t" . join ( [ B [ unstruct ( j ) [ 0 ] ] for j in Si [ line ] ] )
                print >> outfile , line + ( " " * ( longname - len ( line ) + 3 ) ) + "\t" * 6 + "\t" . join ( [ B [ unstruct ( j ) [ 1 ] ] for j in Si [ line ] ] )
        else :
            for line in SF :
                print >> outfile , line + ( " " * ( longname - len ( line ) + 3 ) ) + "\t" * 6 + "\t" . join ( [ B [ unstruct ( j ) [ 1 ] ] for j in Si [ line ] ] )
# # Do linked and unlinked snps in . geno format
if "geno" in formats :
    with open ( ugenoout , 'w' ) as outfile :
        for i in range ( len ( Si . values ( ) [ 0 ] ) ) :
            getref = 0
            ref = "N"
            while ref == "N" :
                ref = unstruct ( Si [ SF [ getref ] ] [ i ] ) [ 0 ]
                getref += 1
            SNProw = "" . join ( map ( str , [ unstruct ( Si [ j ] [ i ] ) . count ( ref ) if Si [ j ] [ i ] != "N" else "9" for j in SF ] ) )
            # # print ref , SNProw
            if len ( set ( SNProw ) ) > 1 :
                print >> outfile , SNProw
    with open ( genoout , 'w' ) as outfile :
        for i in range ( len ( S . values ( ) [ 0 ] ) ) :
            if S [ SF [ 0 ] ] [ i ] . strip ( "_" ) . strip ( ) :
                getref = 0
                ref = "N"
                while ref == "N" : # print i , S [ SF [ 0 ] ] [ i ]
                    ref = unstruct ( S [ SF [ getref ] ] [ i ] ) [ 0 ]
                    getref += 1
                    SNProw = "" . join ( map ( str , [ unstruct ( S [ j ] [ i ] ) . count ( ref ) if S [ j ] [ i ] != "N" else "9" for j in SF ] ) )
                # # print ref , SNProw
                if len ( set ( SNProw ) ) > 1 :
                    print >> outfile , SNProw


signature:
def cluster_info ( ipyclient , spacer = "" ) :
    """reports host and engine info for an ipyclient"""

implementation:
# # get engine data , skips busy engines .
hosts = [ ]
for eid in ipyclient . ids :
    engine = ipyclient [ eid ]
    if not engine . outstanding :
        hosts . append ( engine . apply ( _socket . gethostname ) )
# # report it
hosts = [ i . get ( ) for i in hosts ]
result = [ ]
for hostname in set ( hosts ) :
    result . append ( "{}host compute node: [{} cores] on {}" . format ( spacer , hosts . count ( hostname ) , hostname ) )
print "\n" . join ( result )


signature:
def _debug_on ( ) :
    """Turns on debugging by creating hidden tmp file
 This is only run by the _ _ main _ _ engine ."""

implementation:
# # make tmp file and set loglevel for top - level init
with open ( __debugflag__ , 'w' ) as dfile :
    dfile . write ( "wat" )
__loglevel__ = "DEBUG"
_LOGGER . info ( "debugging turned on and registered to be turned off at exit" )
_set_debug_dict ( __loglevel__ )


signature:
def _set_debug_dict ( __loglevel__ ) :
    """set the debug dict"""

implementation: _lconfig . dictConfig ( { 'version' : 1 , 'disable_existing_loggers' : False , 'formatters' : { 'standard' : { 'format' : "%(asctime)s \t" + "pid=%(process)d \t" + "[%(filename)s]\t" + "%(levelname)s \t" + "%(message)s" } , } , 'handlers' : { __name__ : { 'level' : __loglevel__ , 'class' : 'logging.FileHandler' , 'filename' : __debugfile__ , 'formatter' : "standard" , 'mode' : 'a+' } } , 'loggers' : { __name__ : { 'handlers' : [ __name__ ] , 'level' : __loglevel__ , 'propogate' : True } } } )


signature:
def _debug_off ( ) :
    """turns off debugging by removing hidden tmp file"""

implementation:
if _os . path . exists ( __debugflag__ ) :
    _os . remove ( __debugflag__ )
__loglevel__ = "ERROR"
_LOGGER . info ( "debugging turned off" )
_set_debug_dict ( __loglevel__ )


signature:
def _cmd_exists ( cmd ) :
    """check if dependency program is there"""

implementation: return _subprocess . call ( "type " + cmd , shell = True , stdout = _subprocess . PIPE , stderr = _subprocess . PIPE ) == 0


signature:
def _getbins ( ) :
    """gets the right version of vsearch , muscle , and smalt
 depending on linux vs osx"""

implementation:
# Return error if system is 32 - bit arch .
# This is straight from the python docs :
# https : / / docs . python . org / 2 / library / platform . html # cross - platform
if not _sys . maxsize > 2 ** 32 :
    _sys . exit ( "ipyrad requires 64bit architecture" )
# # get platform mac or linux
_platform = _sys . platform
# # get current location
if 'VIRTUAL_ENV' in _os . environ :
    ipyrad_path = _os . environ [ 'VIRTUAL_ENV' ]
else :
    path = _os . path . abspath ( _os . path . dirname ( __file__ ) )
    ipyrad_path = _os . path . dirname ( path )
# # find bin directory
ipyrad_path = _os . path . dirname ( path )
bin_path = _os . path . join ( ipyrad_path , "bin" )
# # get the correct binaries
if 'linux' in _platform :
    vsearch = _os . path . join ( _os . path . abspath ( bin_path ) , "vsearch-linux-x86_64" )
    muscle = _os . path . join ( _os . path . abspath ( bin_path ) , "muscle-linux-x86_64" )
    smalt = _os . path . join ( _os . path . abspath ( bin_path ) , "smalt-linux-x86_64" )
    bwa = _os . path . join ( _os . path . abspath ( bin_path ) , "bwa-linux-x86_64" )
    samtools = _os . path . join ( _os . path . abspath ( bin_path ) , "samtools-linux-x86_64" )
    bedtools = _os . path . join ( _os . path . abspath ( bin_path ) , "bedtools-linux-x86_64" )
    qmc = _os . path . join ( _os . path . abspath ( bin_path ) , "QMC-linux-x86_64" )
else :
    vsearch = _os . path . join ( _os . path . abspath ( bin_path ) , "vsearch-osx-x86_64" )
    muscle = _os . path . join ( _os . path . abspath ( bin_path ) , "muscle-osx-x86_64" )
    smalt = _os . path . join ( _os . path . abspath ( bin_path ) , "smalt-osx-x86_64" )
    bwa = _os . path . join ( _os . path . abspath ( bin_path ) , "bwa-osx-x86_64" )
    samtools = _os . path . join ( _os . path . abspath ( bin_path ) , "samtools-osx-x86_64" )
    bedtools = _os . path . join ( _os . path . abspath ( bin_path ) , "bedtools-osx-x86_64" )
    # # only one compiled version available , works for all ?
    qmc = _os . path . join ( _os . path . abspath ( bin_path ) , "QMC-osx-x86_64" )
# Test for existence of binaries
assert _cmd_exists ( muscle ) , "muscle not found here: " + muscle
assert _cmd_exists ( vsearch ) , "vsearch not found here: " + vsearch
assert _cmd_exists ( smalt ) , "smalt not found here: " + smalt
assert _cmd_exists ( bwa ) , "bwa not found here: " + bwa
assert _cmd_exists ( samtools ) , "samtools not found here: " + samtools
assert _cmd_exists ( bedtools ) , "bedtools not found here: " + bedtools
# assert _ cmd _ exists ( qmc ) , " wQMC not found here : " + qmc
return vsearch , muscle , smalt , bwa , samtools , bedtools , qmc


signature:
def nworker ( data , chunk ) :
    """Worker to distribute work to jit funcs . Wraps everything on an
 engine to run single - threaded to maximize efficiency for
 multi - processing ."""

implementation:
# # set the thread limit on the remote engine
oldlimit = set_mkl_thread_limit ( 1 )
# # open seqarray view , the modified arr is in bootstarr
with h5py . File ( data . database . input , 'r' ) as io5 :
    seqview = io5 [ "bootsarr" ] [ : ]
    maparr = io5 [ "bootsmap" ] [ : , 0 ]
    smps = io5 [ "quartets" ] [ chunk : chunk + data . _chunksize ]
    # # create an N - mask array of all seq cols
    nall_mask = seqview [ : ] == 78
# # init arrays to fill with results
rquartets = np . zeros ( ( smps . shape [ 0 ] , 4 ) , dtype = np . uint16 )
rinvariants = np . zeros ( ( smps . shape [ 0 ] , 16 , 16 ) , dtype = np . uint16 )
# # fill arrays with results as we compute them . This iterates
# # over all of the quartet sets in this sample chunk . It would
# # be nice to have this all numbified .
for idx in xrange ( smps . shape [ 0 ] ) :
    sidx = smps [ idx ]
    seqs = seqview [ sidx ]
    # # these axis calls cannot be numbafied , but I can ' t
    # # find a faster way that is JIT compiled , and I ' ve
    # # really , really , really tried . Tried again now that
    # # numba supports axis args for np . sum . Still can ' t
    # # get speed improvements by numbifying this loop .
    nmask = np . any ( nall_mask [ sidx ] , axis = 0 )
    nmask += np . all ( seqs == seqs [ 0 ] , axis = 0 )
    # # here are the jitted funcs
    bidx , invar = calculate ( seqs , maparr , nmask , TESTS )
    # # store results
    rquartets [ idx ] = smps [ idx ] [ bidx ]
    rinvariants [ idx ] = invar
# # reset thread limit
set_mkl_thread_limit ( oldlimit )
# # return results . . .
return rquartets , rinvariants


signature:
def store_all ( self ) :
    """Populate array with all possible quartets . This allows us to
 sample from the total , and also to continue from a checkpoint"""

implementation:
with h5py . File ( self . database . input , 'a' ) as io5 :
    fillsets = io5 [ "quartets" ]
    # # generator for all quartet sets
    qiter = itertools . combinations ( xrange ( len ( self . samples ) ) , 4 )
    i = 0
    while i < self . params . nquartets : # # sample a chunk of the next ordered N set of quartets
        dat = np . array ( list ( itertools . islice ( qiter , self . _chunksize ) ) )
        end = min ( self . params . nquartets , dat . shape [ 0 ] + i )
        fillsets [ i : end ] = dat [ : end - i ]
        i += self . _chunksize
        # # send progress update to stdout on engine
        print ( min ( i , self . params . nquartets ) )


signature:
def store_random ( self ) :
    """Populate array with random quartets sampled from a generator .
 Holding all sets in memory might take a lot , but holding a very
 large list of random numbers for which ones to sample will fit
 into memory for most reasonable sized sets . So we ' ll load a
 list of random numbers in the range of the length of total
 sets that can be generated , then only keep sets from the set
 generator if they are in the int list . I did several tests to
 check that random pairs are as likely as 0 & 1 to come up together
 in a random quartet set ."""

implementation:
with h5py . File ( self . database . input , 'a' ) as io5 :
    fillsets = io5 [ "quartets" ]
    # # set generators
    qiter = itertools . combinations ( xrange ( len ( self . samples ) ) , 4 )
    rand = np . arange ( 0 , n_choose_k ( len ( self . samples ) , 4 ) )
    np . random . shuffle ( rand )
    rslice = rand [ : self . params . nquartets ]
    rss = np . sort ( rslice )
    riter = iter ( rss )
    del rand , rslice
    # # print progress update 1 to the engine stdout
    print ( self . _chunksize )
    # # set to store
    rando = riter . next ( )
    tmpr = np . zeros ( ( self . params . nquartets , 4 ) , dtype = np . uint16 )
    tidx = 0
    while 1 :
        try :
            for i , j in enumerate ( qiter ) :
                if i == rando :
                    tmpr [ tidx ] = j
                    tidx += 1
                    rando = riter . next ( )
                # # print progress bar update to engine stdout
                if not i % self . _chunksize :
                    print ( min ( i , self . params . nquartets ) )
        except StopIteration :
            break
    # # store into database
    fillsets [ : ] = tmpr
    del tmpr


signature:
def store_equal ( self ) :
    """Takes a tetrad class object and populates array with random
 quartets sampled equally among splits of the tree so that
 deep splits are not overrepresented relative to rare splits ,
 like those near the tips ."""

implementation:
with h5py . File ( self . database . input , 'a' ) as io5 :
    fillsets = io5 [ "quartets" ]
    # # require guidetree
    if not os . path . exists ( self . files . tree ) :
        raise IPyradWarningExit ( "To use sampling method 'equal' requires a guidetree" )
    tre = ete3 . Tree ( self . files . tree )
    tre . unroot ( )
    tre . resolve_polytomy ( recursive = True )
    # # randomly sample internals splits
    splits = [ ( [ self . samples . index ( z . name ) for z in i ] , [ self . samples . index ( z . name ) for z in j ] ) for ( i , j ) in tre . get_edges ( ) ]
    # # only keep internal splits , not single tip edges
    splits = [ i for i in splits if all ( [ len ( j ) > 1 for j in i ] ) ]
    # # how many min quartets shoudl be equally sampled from each split
    squarts = self . params . nquartets // len ( splits )
    # # keep track of how many iterators are saturable .
    saturable = 0
    # # turn each into an iterable split sampler
    # # if the nquartets for that split is small , then sample all ,
    # # if it is big then make it a random sampler for that split .
    qiters = [ ]
    # # iterate over splits sampling quartets evenly
    for idx , split in enumerate ( splits ) : # # if small number at this split then sample all possible sets
    # # we will exhaust this quickly and then switch to random for
    # # the larger splits .
        total = n_choose_k ( len ( split [ 0 ] ) , 2 ) * n_choose_k ( len ( split [ 1 ] ) , 2 )
        if total < squarts * 2 :
            qiter = ( i + j for ( i , j ) in itertools . product ( itertools . combinations ( split [ 0 ] , 2 ) , itertools . combinations ( split [ 1 ] , 2 ) ) )
            saturable += 1
        # # else create random sampler across that split , this is slower
        # # because it can propose the same split repeatedly and so we
        # # have to check it against the ' sampled ' set .
        else :
            qiter = ( random_product ( split [ 0 ] , split [ 1 ] ) for _ in xrange ( self . params . nquartets ) )
        # # store all iterators into a list
        qiters . append ( ( idx , qiter ) )
    # # create infinite cycler of qiters
    qitercycle = itertools . cycle ( qiters )
    # # store visited quartets
    sampled = set ( )
    # # fill chunksize at a time
    i = 0
    empty = set ( )
    edge_targeted = 0
    random_targeted = 0
    # # keep filling quartets until nquartets are sampled .
    while i < self . params . nquartets : # # grab the next iterator
        cycle , qiter = qitercycle . next ( )
        # # sample from iterators , store sorted set .
        try :
            qrtsamp = tuple ( sorted ( qiter . next ( ) ) )
            if qrtsamp not in sampled :
                sampled . add ( qrtsamp )
                edge_targeted += 1
                i += 1
                # # print progress bar update to engine stdout
                if not i % self . _chunksize :
                    print ( min ( i , self . params . nquartets ) )
        except StopIteration :
            empty . add ( cycle )
            if len ( empty ) == saturable :
                break
    # # if array is not full then add random samples
    while i <= self . params . nquartets :
        newset = tuple ( sorted ( np . random . choice ( range ( len ( self . samples ) ) , 4 , replace = False ) ) )
        if newset not in sampled :
            sampled . add ( newset )
            random_targeted += 1
            i += 1
            # # print progress bar update to engine stdout
            if not i % self . _chunksize :
                print ( min ( i , self . params . nquartets ) )
    # # store into database
    print ( self . params . nquartets )
    fillsets [ : ] = np . array ( tuple ( sampled ) )
    del sampled


signature:
def random_combination ( nsets , n , k ) :
    """Returns nsets unique random quartet sets sampled from
 n - choose - k without replacement combinations ."""

implementation:
sets = set ( )
while len ( sets ) < nsets :
    newset = tuple ( sorted ( np . random . choice ( n , k , replace = False ) ) )
    sets . add ( newset )
return tuple ( sets )


signature:
def random_product ( iter1 , iter2 ) :
    """Random sampler for equal _ splits functions"""

implementation:
iter4 = np . concatenate ( [ np . random . choice ( iter1 , 2 , replace = False ) , np . random . choice ( iter2 , 2 , replace = False ) ] )
return iter4


signature:
def resolve_ambigs ( tmpseq ) :
    """Randomly resolve ambiguous bases . This is applied to each boot
 replicate so that over reps the random resolutions don ' t matter .
 Sites are randomly resolved , so best for unlinked SNPs since
 otherwise linked SNPs are losing their linkage information . . .
 though it ' s not like we ' re using it anyways ."""

implementation:
# # the order of rows in GETCONS
for aidx in xrange ( 6 ) : # np . uint ( [ 82 , 75 , 83 , 89 , 87 , 77 ] ) :
    ambig , res1 , res2 = GETCONS [ aidx ]
    # # get true wherever tmpseq is ambig
    idx , idy = np . where ( tmpseq == ambig )
    halfmask = np . random . choice ( np . array ( [ True , False ] ) , idx . shape [ 0 ] )
    for col in xrange ( idx . shape [ 0 ] ) :
        if halfmask [ col ] :
            tmpseq [ idx [ col ] , idy [ col ] ] = res1
        else :
            tmpseq [ idx [ col ] , idy [ col ] ] = res2
return tmpseq


signature:
def set_mkl_thread_limit ( cores ) :
    """set mkl thread limit and return old value so we can reset
 when finished ."""

implementation:
if "linux" in sys . platform :
    mkl_rt = ctypes . CDLL ( 'libmkl_rt.so' )
else :
    mkl_rt = ctypes . CDLL ( 'libmkl_rt.dylib' )
oldlimit = mkl_rt . mkl_get_max_threads ( )
mkl_rt . mkl_set_num_threads ( ctypes . byref ( ctypes . c_int ( cores ) ) )
return oldlimit


signature:
def get_total ( tots , node ) :
    """get total number of quartets possible for a split"""

implementation:
if ( node . is_leaf ( ) or node . is_root ( ) ) :
    return 0
else : # # Get counts on down edges .
# # How to treat polytomies here ?
    if len ( node . children ) > 2 :
        down_r = node . children [ 0 ]
        down_l = node . children [ 1 ]
        for child in node . children [ 2 : ] :
            down_l += child
    else :
        down_r , down_l = node . children
    lendr = sum ( 1 for i in down_r . iter_leaves ( ) )
    lendl = sum ( 1 for i in down_l . iter_leaves ( ) )
    # # get count on up edge sister
    up_r = node . get_sisters ( ) [ 0 ]
    lenur = sum ( 1 for i in up_r . iter_leaves ( ) )
    # # everyone else
    lenul = tots - ( lendr + lendl + lenur )
    # # return product
    return lendr * lendl * lenur * lenul


signature:
def get_sampled ( data , totn , node ) :
    """get total number of quartets sampled for a split"""

implementation:
# # convert tip names to ints
names = sorted ( totn )
cdict = { name : idx for idx , name in enumerate ( names ) }
# # skip some nodes
if ( node . is_leaf ( ) or node . is_root ( ) ) :
    return 0
else : # # get counts on down edges
    if len ( node . children ) > 2 :
        down_r = node . children [ 0 ]
        down_l = node . children [ 1 ]
        for child in node . children [ 2 : ] :
            down_l += child
    else :
        down_r , down_l = node . children
    lendr = set ( cdict [ i ] for i in down_r . get_leaf_names ( ) )
    lendl = set ( cdict [ i ] for i in down_l . get_leaf_names ( ) )
    # # get count on up edge sister
    up_r = node . get_sisters ( ) [ 0 ]
    lenur = set ( cdict [ i ] for i in up_r . get_leaf_names ( ) )
    # # everyone else
    lenul = set ( cdict [ i ] for i in totn ) - set . union ( lendr , lendl , lenur )
idx = 0
sampled = 0
with h5py . File ( data . database . output , 'r' ) as io5 :
    end = io5 [ "quartets" ] . shape [ 0 ]
    while 1 : # # break condition
        if idx >= end :
            break
        # # counts matches
        qrts = io5 [ "quartets" ] [ idx : idx + data . _chunksize ]
        for qrt in qrts :
            sqrt = set ( qrt )
            if all ( [ sqrt . intersection ( i ) for i in [ lendr , lendl , lenur , lenul ] ] ) :
                sampled += 1
        # # increase span
        idx += data . _chunksize
return sampled


signature:
def consensus_tree ( trees , names = None , cutoff = 0.0 ) :
    """An extended majority rule consensus function for ete3.
 Modelled on the similar function from scikit - bio tree module . If
 cutoff = 0.5 then it is a normal majority rule consensus , while if
 cutoff = 0.0 then subsequent non - conflicting clades are added to the tree ."""

implementation:
# # find which clades occured with freq > cutoff
namedict , clade_counts = find_clades ( trees , names = names )
# # filter out the < cutoff clades
fclade_counts = filter_clades ( clade_counts , cutoff )
# # build tree
consens_tree , _ = build_trees ( fclade_counts , namedict )
# # make sure no singleton nodes were left behind
return consens_tree , clade_counts


signature:
def find_clades ( trees , names ) :
    """A subfunc of consensus _ tree ( ) . Traverses trees to count clade occurrences .
 Names are ordered by names , else they are in the order of the first
 tree ."""

implementation:
# # index names from the first tree
if not names :
    names = trees [ 0 ] . get_leaf_names ( )
ndict = { j : i for i , j in enumerate ( names ) }
namedict = { i : j for i , j in enumerate ( names ) }
# # store counts
clade_counts = defaultdict ( int )
# # count as bitarray clades in each tree
for tree in trees :
    tree . unroot ( )
    for node in tree . traverse ( 'postorder' ) : # bits = bitarray ( ' 0 ' * len ( tree ) )
        bits = np . zeros ( len ( tree ) , dtype = np . bool_ )
        for child in node . iter_leaf_names ( ) :
            bits [ ndict [ child ] ] = True
        # # if parent is root then mirror flip one child ( where bit [ 0 ] = 0)
        # if not node . is _ root ( ) :
        # if node . up . is _ root ( ) :
        # if bits [ 0 ] :
        # bits . invert ( )
        bitstring = "" . join ( [ np . binary_repr ( i ) for i in bits ] )
        clade_counts [ bitstring ] += 1
# # convert to freq
for key , val in clade_counts . items ( ) :
    clade_counts [ key ] = val / float ( len ( trees ) )
# # return in sorted order
clade_counts = sorted ( clade_counts . items ( ) , key = lambda x : x [ 1 ] , reverse = True )
return namedict , clade_counts


signature:
def build_trees ( fclade_counts , namedict ) :
    """A subfunc of consensus _ tree ( ) . Build an unrooted consensus tree
 from filtered clade counts ."""

implementation:
# # storage
nodes = { }
idxarr = np . arange ( len ( fclade_counts [ 0 ] [ 0 ] ) )
queue = [ ]
# # create dict of clade counts and set keys
countdict = defaultdict ( int )
for clade , count in fclade_counts :
    mask = np . int_ ( list ( clade ) ) . astype ( np . bool )
    ccx = idxarr [ mask ]
    queue . append ( ( len ( ccx ) , frozenset ( ccx ) ) )
    countdict [ frozenset ( ccx ) ] = count
while queue :
    queue . sort ( )
    ( clade_size , clade ) = queue . pop ( 0 )
    new_queue = [ ]
    # search for ancestors of clade
    for ( _ , ancestor ) in queue :
        if clade . issubset ( ancestor ) : # update ancestor such that , in the following example :
        # ancestor = = { 1 , 2 , 3 , 4}
        # clade = = { 2 , 3}
        # new _ ancestor = = { 1 , { 2 , 3 } , 4}
            new_ancestor = ( ancestor - clade ) | frozenset ( [ clade ] )
            countdict [ new_ancestor ] = countdict . pop ( ancestor )
            ancestor = new_ancestor
        new_queue . append ( ( len ( ancestor ) , ancestor ) )
    # if the clade is a tip , then we have a name
    if clade_size == 1 :
        name = list ( clade ) [ 0 ]
        name = namedict [ name ]
    else :
        name = None
    # the clade will not be in nodes if it is a tip
    children = [ nodes . pop ( c ) for c in clade if c in nodes ]
    node = ete3 . Tree ( name = name )
    # node = toytree . tree ( name = name ) . tree
    for child in children :
        node . add_child ( child )
    if not node . is_leaf ( ) :
        node . dist = int ( round ( 100 * countdict [ clade ] ) )
        node . support = int ( round ( 100 * countdict [ clade ] ) )
    else :
        node . dist = int ( 100 )
        node . support = int ( 100 )
    nodes [ clade ] = node
    queue = new_queue
tre = nodes . values ( ) [ 0 ]
tre . unroot ( )
# # return the tree and other trees if present
return tre , list ( nodes . values ( ) )


signature:
def _refresh ( self ) :
    """Remove all existing results files and reinit the h5 arrays
 so that the tetrad object is just like fresh from a CLI start ."""

implementation:
# # clear any existing results files
oldfiles = [ self . files . qdump ] + self . database . __dict__ . values ( ) + self . trees . __dict__ . values ( )
for oldfile in oldfiles :
    if oldfile :
        if os . path . exists ( oldfile ) :
            os . remove ( oldfile )
# # store old ipcluster info
oldcluster = copy . deepcopy ( self . _ipcluster )
# # reinit the tetrad object data .
self . __init__ ( name = self . name , data = self . files . data , mapfile = self . files . mapfile , workdir = self . dirs , method = self . params . method , guidetree = self . files . tree , resolve_ambigs = self . params . resolve_ambigs , save_invariants = self . params . save_invariants , nboots = self . params . nboots , nquartets = self . params . nquartets , initarr = True , quiet = True , cli = self . kwargs . get ( "cli" ) )
# # retain the same ipcluster info
self . _ipcluster = oldcluster


signature:
def _store_N_samples ( self , start , ipyclient , quiet = False ) :
    """Find all quartets of samples and store in a large array
 A chunk size is assigned for sampling from the array of quartets
 based on the number of cpus available . This should be relatively
 large so that we don ' t spend a lot of time doing I / O , but small
 enough that jobs finish often for checkpointing ."""

implementation:
breaks = 2
if self . params . nquartets < 5000 :
    breaks = 1
if self . params . nquartets > 100000 :
    breaks = 8
if self . params . nquartets > 500000 :
    breaks = 16
if self . params . nquartets > 5000000 :
    breaks = 32
# # chunk up the data
ncpus = len ( ipyclient )
self . _chunksize = ( self . params . nquartets // ( breaks * ncpus ) + ( self . params . nquartets % ( breaks * ncpus ) ) )
# # create h5 OUT empty arrays
# # ' quartets ' stores the inferred quartet relationship ( 1 x 4)
# # This can get huge , so we need to choose the dtype wisely .
# # the values are simply the index of the taxa , so uint16 is good .
with h5py . File ( self . database . output , 'w' ) as io5 :
    io5 . create_dataset ( "quartets" , ( self . params . nquartets , 4 ) , dtype = np . uint16 , chunks = ( self . _chunksize , 4 ) )
    # # group for bootstrap invariant matrices ( ( 16 , 16 ) , uint32)
    # # these store the actual matrix counts . dtype uint32 can store
    # # up to 4294967295 . More than enough . uint16 max is 65535.
    # # the 0 boot is the original seqarray .
    io5 . create_group ( "invariants" )
# # append to h5 IN array ( which has the seqarray , bootsarr , maparr )
# # and fill it with all of the quartet sets we will ever sample .
# # the samplign method will vary depending on whether this is random ,
# # all , or equal splits ( in a separate but similar function ) .
with h5py . File ( self . database . input , 'a' ) as io5 :
    try :
        io5 . create_dataset ( "quartets" , ( self . params . nquartets , 4 ) , dtype = np . uint16 , chunks = ( self . _chunksize , 4 ) , compression = 'gzip' )
    except RuntimeError :
        raise IPyradWarningExit ( "database file already exists for this analysis, " + "you must run with the force flag to overwrite" )
# # submit store job to write into self . database . input
if self . params . method == "all" :
    async = ipyclient [ 0 ] . apply ( store_all , self )
elif self . params . method == "random" :
    async = ipyclient [ 0 ] . apply ( store_random , self )
elif self . params . method == "equal" :
    async = ipyclient [ 0 ] . apply ( store_equal , self )
# # progress bar
printstr = "generating q-sets | {} | "
prog = 0
while 1 :
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    if not quiet :
        if async . stdout :
            prog = int ( async . stdout . strip ( ) . split ( ) [ - 1 ] )
        progressbar ( self . params . nquartets , prog , printstr . format ( elapsed ) , spacer = "" )
    if not async . ready ( ) :
        time . sleep ( 0.1 )
    else :
        break
if not async . successful ( ) :
    raise IPyradWarningExit ( async . result ( ) )
if not quiet :
    print ( "" )


signature:
def _dump_qmc ( self ) :
    """Writes the inferred quartet sets from the database to a text
 file to be used as input for QMC . Quartets that had no information
 available ( i . e . , no SNPs ) were written to the database as 0,0,0,0
 and are excluded here from the output ."""

implementation:
# # open the h5 database
with h5py . File ( self . database . output , 'r' ) as io5 : # # create an output file for writing
    self . files . qdump = os . path . join ( self . dirs , self . name + ".quartets.txt" )
    with open ( self . files . qdump , 'w' ) as qdump : # # pull from db
        for idx in xrange ( 0 , self . params . nquartets , self . _chunksize ) :
            qchunk = io5 [ "quartets" ] [ idx : idx + self . _chunksize , : ]
            quarts = [ tuple ( j ) for j in qchunk if np . any ( j ) ]
            # # shuffle and format for qmc
            np . random . shuffle ( quarts )
            chunk = [ "{},{}|{},{}" . format ( * i ) for i in quarts ]
            qdump . write ( "\n" . join ( chunk ) + "\n" )


signature:
def _run_qmc ( self , boot ) :
    """Runs quartet max - cut QMC on the quartets qdump file ."""

implementation:
# # build command
self . _tmp = os . path . join ( self . dirs , ".tmptre" )
cmd = [ ip . bins . qmc , "qrtt=" + self . files . qdump , "otre=" + self . _tmp ]
# # run it
proc = subprocess . Popen ( cmd , stderr = subprocess . STDOUT , stdout = subprocess . PIPE )
res = proc . communicate ( )
if proc . returncode :
    raise IPyradWarningExit ( res [ 1 ] )
# # parse tmp file written by qmc into a tree and rename it
with open ( self . _tmp , 'r' ) as intree :
    tre = ete3 . Tree ( intree . read ( ) . strip ( ) )
    names = tre . get_leaves ( )
    for name in names :
        name . name = self . samples [ int ( name . name ) ]
    tmptre = tre . write ( format = 9 )
# # save the tree to file
if boot :
    self . trees . boots = os . path . join ( self . dirs , self . name + ".boots" )
    with open ( self . trees . boots , 'a' ) as outboot :
        outboot . write ( tmptre + "\n" )
else :
    self . trees . tree = os . path . join ( self . dirs , self . name + ".tree" )
    with open ( self . trees . tree , 'w' ) as outtree :
        outtree . write ( tmptre )
# # save the file
self . _save ( )


signature:
def _compute_stats ( self , start , ipyclient , quiet = False ) :
    """Compute sampling stats and consens trees ."""

implementation:
# # get name indices
names = self . samples
# # make a consensus from bootstrap reps .
if self . checkpoint . boots :
    tre = ete3 . Tree ( self . trees . tree , format = 0 )
    tre . unroot ( )
    with open ( self . trees . boots , 'r' ) as inboots :
        bb = [ ete3 . Tree ( i . strip ( ) , format = 0 ) for i in inboots . readlines ( ) ]
        bb = [ tre ] + bb
    # # calculate consensus supports
    ctre , counts = consensus_tree ( bb , names = names )
    self . trees . cons = os . path . join ( self . dirs , self . name + ".cons" )
    with open ( self . trees . cons , 'w' ) as ocons :
        ocons . write ( ctre . write ( format = 0 ) )
else :
    ctre = ete3 . Tree ( self . trees . tree , format = 0 )
    ctre . unroot ( )
# # build stats file and write trees
self . trees . nhx = os . path . join ( self . dirs , self . name + ".nhx" )
lbview = ipyclient . load_balanced_view ( )
qtots = { }
qsamp = { }
tots = sum ( 1 for i in ctre . iter_leaves ( ) )
totn = set ( ctre . get_leaf_names ( ) )
# # iterate over node traversal
for node in ctre . traverse ( ) : # # this is slow , needs to look at every sampled quartet
# # so we send it to be processed on engines
    qtots [ node ] = lbview . apply ( get_total , * ( tots , node ) )
    qsamp [ node ] = lbview . apply ( get_sampled , * ( self , totn , node ) )
# # wait for jobs to finish ( + 1 to lenjob is for final progress printer )
alljobs = qtots . values ( ) + qsamp . values ( )
lenjobs = len ( alljobs ) + 1
printstr = "calculating stats | {} | "
done = 0
while 1 :
    if not quiet :
        done = sum ( [ i . ready ( ) for i in alljobs ] )
        elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
        progressbar ( lenjobs , done , printstr . format ( elapsed ) , spacer = "" )
    if ( lenjobs - 1 ) == done :
        break
    else :
        time . sleep ( 0.1 )
# # store results in the tree object
for node in ctre . traverse ( ) :
    total = qtots [ node ] . result ( )
    sampled = qsamp [ node ] . result ( )
    node . add_feature ( "quartets_total" , total )
    node . add_feature ( "quartets_sampled" , sampled )
features = [ "quartets_total" , "quartets_sampled" ]
# # update final progress
elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
progressbar ( 1 , 1 , printstr . format ( elapsed ) , spacer = "" )
if not quiet :
    print ( "" )
# # write tree in NHX format
with open ( self . trees . nhx , 'w' ) as outtre :
    outtre . write ( ctre . write ( format = 0 , features = features ) )


signature:
def _load ( self , name , workdir , quiet = False ) :
    """Load a JSON serialized tetrad instance to continue from a checkpoint ."""

implementation:
# # load the JSON string and try with name + . json
path = os . path . join ( workdir , name )
if not path . endswith ( ".tet.json" ) :
    path += ".tet.json"
# # expand user
path = path . replace ( "~" , os . path . expanduser ( "~" ) )
# # load the json file as a dictionary
try :
    with open ( path , 'r' ) as infile :
        fullj = _byteify ( json . loads ( infile . read ( ) , object_hook = _byteify ) , ignore_dicts = True )
except IOError :
    raise IPyradWarningExit ( """\
        Cannot find checkpoint (.tet.json) file at: {}""" . format ( path ) )
# # set old attributes into new tetrad object
self . name = fullj [ "name" ]
self . files . data = fullj [ "files" ] [ "data" ]
self . files . mapfile = fullj [ "files" ] [ "mapfile" ]
self . dirs = fullj [ "dirs" ]
self . _init_seqarray ( quiet = quiet )
self . _parse_names ( )
# # fill in the same attributes
for key in fullj : # # fill Params a little different
    if key in [ "files" , "params" , "database" , "trees" , "stats" , "checkpoint" ] :
        filler = fullj [ key ]
        for ikey in filler :
            self . __dict__ [ key ] . __setattr__ ( ikey , fullj [ key ] [ ikey ] )
    else :
        self . __setattr__ ( key , fullj [ key ] )


signature:
def _inference ( self , start , ipyclient , quiet ) :
    """Sends slices of quartet sets to parallel engines for computing ,
 enters results into output database , sends finished quartet sets
 to QMC for tree inference , and prints progress bars ."""

implementation:
# # load - balancer for single - threaded execution jobs
lbview = ipyclient . load_balanced_view ( )
# # an iterator that grabs quartet chunk start positions
jobs = range ( self . checkpoint . arr , self . params . nquartets , self . _chunksize )
# # if this is a bootstrap then init a new boot array in the database
# # max val is 65535 in here if uint16
bootkey = "boot{}" . format ( self . checkpoint . boots )
with h5py . File ( self . database . output , 'r+' ) as io5 :
    if bootkey not in io5 [ "invariants" ] . keys ( ) :
        io5 [ "invariants" ] . create_dataset ( bootkey , ( self . params . nquartets , 16 , 16 ) , dtype = np . uint16 , chunks = ( self . _chunksize , 16 , 16 ) )
# # start progress bar if new or skip if bootstrapping
elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
if self . checkpoint . boots :
    printstr = "bootstrap trees   | {} | "
else :
    printstr = "initial tree      | {} | "
    if not quiet :
        progressbar ( 1 , 0 , printstr . format ( elapsed ) , spacer = "" )
# # submit jobs distriuted across the cluster .
asyncs = { }
for job in jobs :
    asyncs [ job ] = lbview . apply ( nworker , * ( self , job ) )
# # wait for jobs to finish , catch results as they return and
# # enter them into the HDF5 database to keep memory low .
done = 0
while 1 : # # gather finished jobs
    finished = [ i for i , j in asyncs . iteritems ( ) if j . ready ( ) ]
    # # iterate over finished list
    for key in finished :
        async = asyncs [ key ]
        if async . successful ( ) : # # store result
            done += 1
            results = async . result ( )
            self . _insert_to_array ( key , results )
            # # purge from memory
            del asyncs [ key ]
        else :
            raise IPyradWarningExit ( async . result ( ) )
    # # progress bar is different if first vs boot tree
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    if not self . checkpoint . boots :
        if not quiet :
            progressbar ( len ( jobs ) , done , printstr . format ( elapsed ) , spacer = "" )
    else :
        if not quiet :
            progressbar ( self . params . nboots , self . checkpoint . boots , printstr . format ( elapsed ) , spacer = "" )
    # # done is counted on finish , so this means we ' re done
    if len ( asyncs ) == 0 :
        break
    else :
        time . sleep ( 0.1 )
# # dump quartets into a text file for QMC
self . _dump_qmc ( )
# # send to QMC
if not self . checkpoint . boots :
    self . _run_qmc ( 0 )
else :
    self . _run_qmc ( 1 )
# # reset the checkpoint arr
self . checkpoint . arr = 0
# # print spacer if finished first tree or last boot .
if ( not self . checkpoint . boots ) and ( not quiet ) :
    print ( "" )
elif ( self . checkpoint . boots == self . params . nboots ) and ( not quiet ) :
    print ( "" )


signature:
def _insert_to_array ( self , chunk , results ) :
    """Enters results arrays into the HDF5 database ."""

implementation:
# # two result arrs
chunksize = self . _chunksize
qrts , invs = results
# # enter into db
with h5py . File ( self . database . output , 'r+' ) as io5 :
    io5 [ 'quartets' ] [ chunk : chunk + chunksize ] = qrts
    # # entered as 0 - indexed !
    if self . params . save_invariants :
        if self . checkpoint . boots :
            key = "invariants/boot{}" . format ( self . checkpoint . boots )
            io5 [ key ] [ chunk : chunk + chunksize ] = invs
        else :
            io5 [ "invariants/boot0" ] [ chunk : chunk + chunksize ] = invs


signature:
def run ( self , force = False , quiet = False , ipyclient = None ) :
    """Parameters
 force ( bool ) :
 Overwrite existing results for object with the same name
 and workdir as this one .
 verbose ( int ) :
 0 = primt nothing ; 1 = print progress bars ; 2 = print pringress
 bars and cluster information .
 ipyclient ( ipyparallel . Client object ) :
 A connected ipyclient object . If ipcluster instance is
 not running on the default profile then . . ."""

implementation:
# # force overwrite needs to clear out the HDF5 database
if force :
    self . _refresh ( )
# # print nquartet statement
if not quiet :
    print ( "inferring {} quartet tree sets" . format ( self . params . nquartets ) )
# # wrap the run in a try statement to ensure we properly shutdown
# # and cleanup on exit or interrupt .
inst = None
try : # # find and connect to an ipcluster instance given the information
# # in the _ ipcluster dictionary if a connected client was not given .
    if not ipyclient :
        args = self . _ipcluster . items ( ) + [ ( "spacer" , "" ) ]
        ipyclient = ip . core . parallel . get_client ( ** dict ( args ) )
    # # print the cluster connection information
    if not quiet :
        ip . cluster_info ( ipyclient )
    # # store ipyclient engine pids to the dict so we can
    # # hard - interrupt them later if assembly is interrupted .
    # # Only stores pids of engines that aren ' t busy at this moment ,
    # # otherwise it would block here while waiting to find their pids .
    self . _ipcluster [ "pids" ] = { }
    for eid in ipyclient . ids :
        engine = ipyclient [ eid ]
        if not engine . outstanding :
            pid = engine . apply ( os . getpid ) . get ( )
            self . _ipcluster [ "pids" ] [ eid ] = pid
    # # fill the input array with quartets to sample - - - - -
    start = time . time ( )
    if ( not self . checkpoint . boots ) and ( not self . trees . tree ) :
        self . _store_N_samples ( start , ipyclient , quiet = quiet )
    # # calculate invariants for the full seqarray - - - - -
    start = time . time ( )
    if not self . trees . tree :
        self . _inference ( start , ipyclient , quiet = quiet )
    else :
        if not quiet :
            print ( "initial tree already inferred" )
    # # calculate invariants for each bootstrap rep - - - - -
    start = time . time ( )
    if self . params . nboots :
        if self . checkpoint . boots : # < = self . params . nboots :
            if not quiet :
                print ( "{} bootstrap trees already inferred" . format ( self . checkpoint . boots ) )
        while self . checkpoint . boots < self . params . nboots : # # resample the bootseq array
            if self . files . mapfile :
                self . _sample_bootseq_array_map ( )
            else :
                self . _sample_bootseq_array ( )
            # # start boot inference
            self . checkpoint . boots += 1
            self . _inference ( start , ipyclient , quiet = quiet )
    # # write output stats - - - - -
    # self . files . stats = os . path . join ( self . dirs , self . name + " _ stats . txt " )
    start = time . time ( )
    self . _compute_stats ( start , ipyclient , quiet = quiet )
# # handle exceptions so they will be raised after we clean up below
except KeyboardInterrupt as inst :
    print ( "\nKeyboard Interrupt by user. Cleaning up..." )
except IPyradWarningExit as inst :
    print ( "\nError encountered: {}" . format ( inst ) )
except Exception as inst :
    print ( "\nUnknown exception encountered: {}" . format ( inst ) )
# # close client when done or interrupted
finally :
    try : # # save the Assembly
        self . _save ( )
        # # can ' t close client if it was never open
        if ipyclient : # # send SIGINT ( 2 ) to all engines
            ipyclient . abort ( )
            time . sleep ( 1 )
            for engine_id , pid in self . _ipcluster [ "pids" ] . items ( ) :
                if ipyclient . queue_status ( ) [ engine_id ] [ "tasks" ] :
                    os . kill ( pid , 2 )
                time . sleep ( 0.25 )
            # # if CLI , stop jobs and shutdown
            if 'ipyrad-cli' in self . _ipcluster [ "cluster_id" ] :
                ipyclient . shutdown ( hub = True , block = False )
                ipyclient . close ( )
            else :
                if not ipyclient . outstanding :
                    ipyclient . purge_everything ( )
                else : # # nanny : kill everything , something bad happened
                    ipyclient . shutdown ( hub = True , block = False )
                    ipyclient . close ( )
                    print ( "\nwarning: ipcluster shutdown and must be restarted" )
        # # reraise the error now that we ' re cleaned up
        # if inst :
        # raise inst
    # # if exception during shutdown then we really screwed up
    except Exception as inst2 :
        print ( "warning: error during shutdown:\n{}" . format ( inst2 ) )


signature:
def start_ipcluster ( data ) :
    """Start ipcluster"""

implementation:
# # if MPI argument then use - - ip arg to view all sockets
iparg = ""
if "MPI" in data . _ipcluster [ "engines" ] :
    iparg = "--ip=*"
# # make ipcluster arg call
standard = """
        ipcluster start
                  --daemonize
                  --cluster-id={}
                  --engines={}
                  --profile={}
                  --n={}
                  {}""" . format ( data . _ipcluster [ "cluster_id" ] , data . _ipcluster [ "engines" ] , data . _ipcluster [ "profile" ] , data . _ipcluster [ "cores" ] , iparg )
# # wrap ipcluster start
try :
    LOGGER . info ( shlex . split ( standard ) )
    subprocess . check_call ( shlex . split ( standard ) , stderr = subprocess . STDOUT , stdout = subprocess . PIPE )
except subprocess . CalledProcessError as inst :
    LOGGER . debug ( "  ipcontroller already running." )
    raise
except Exception as inst :
    sys . exit ( "  Error launching ipcluster for parallelization:\n({})\n" . format ( inst ) )


signature:
def register_ipcluster ( data ) :
    """The name is a unique id that keeps this _ _ init _ _ of ipyrad distinct
 from interfering with other ipcontrollers . Run statements are wrapped
 so that ipcluster will be killed on exit ."""

implementation:
# # check if this pid already has a running cluster
data . _ipcluster [ "cluster_id" ] = "ipyrad-cli-" + str ( os . getpid ( ) )
start_ipcluster ( data )
return data


signature:
def get_client ( cluster_id , profile , engines , timeout , cores , quiet , spacer , ** kwargs ) :
    """Creates a client to view ipcluster engines for a given profile and
 returns it with at least one engine spun up and ready to go . If no
 engines are found after nwait amount of time then an error is raised .
 If engines = = MPI it waits a bit longer to find engines . If the number
 of engines is set then it waits even longer to try to find that number
 of engines ."""

implementation:
# # save stds for later , we ' re gonna hide them to prevent external printing
save_stdout = sys . stdout
save_stderr = sys . stderr
sys . stdout = cStringIO . StringIO ( )
sys . stderr = cStringIO . StringIO ( )
# # get cluster _ info print string
connection_string = "{}establishing parallel connection:" . format ( spacer )
# # wrapped search for ipcluster
try : # # are we looking for a running ipcluster instance ?
    if profile not in [ None , "default" ] :
        args = { 'profile' : profile , "timeout" : timeout }
    else :
        clusterargs = [ cluster_id , profile , timeout ]
        argnames = [ "cluster_id" , "profile" , "timeout" ]
        args = { key : value for key , value in zip ( argnames , clusterargs ) }
    # # get connection within timeout window of wait time and hide messages
    ipyclient = ipp . Client ( ** args )
    sys . stdout = save_stdout
    sys . stderr = save_stderr
    # # check that all engines have connected
    if ( engines == "MPI" ) or ( "ipyrad-cli-" in cluster_id ) :
        if not quiet :
            print ( connection_string )
    for _ in range ( 6000 ) :
        initid = len ( ipyclient )
        time . sleep ( 0.01 )
        # # If MPI then wait for all engines to start so we can report
        # # how many cores are on each host . If Local then only wait for
        # # one engine to be ready and then just go .
        if ( engines == "MPI" ) or ( "ipyrad-cli-" in cluster_id ) : # # wait for cores to be connected
            if cores :
                time . sleep ( 0.1 )
                if initid == cores :
                    break
            if initid :
                time . sleep ( 3 )
                if len ( ipyclient ) == initid :
                    break
        else :
            if cores :
                if initid == cores :
                    break
            else :
                if initid :
                    break
except KeyboardInterrupt as inst : # # ensure stdout is reset even if Exception was raised
    sys . stdout = save_stdout
    sys . stderr = save_stderr
    raise inst
# # This is raised if ipcluster is not running - - - - -
except IOError as inst : # # ensure stdout is reset even if Exception was raised
    sys . stdout = save_stdout
    sys . stderr = save_stderr
    if "ipyrad-cli-" in cluster_id :
        raise IPyradWarningExit ( NO_IPCLUSTER_CLI )
    else :
        raise IPyradWarningExit ( NO_IPCLUSTER_API )
except ( ipp . TimeoutError , ipp . NoEnginesRegistered ) as inst : # # raised by ipp if no connection file is found for ' nwait ' seconds
    sys . stdout = save_stdout
    sys . stderr = save_stderr
    raise inst
except Exception as inst : # # if any other exceptions were missed . . .
    sys . stdout = save_stdout
    sys . stderr = save_stderr
    raise inst
finally : # # ensure that no matter what we reset the stds
    sys . stdout = save_stdout
    sys . stderr = save_stderr
return ipyclient


signature:
def ambigcutters ( seq ) :
    """Returns both resolutions of a cut site that has an ambiguous base in
 it , else the single cut site"""

implementation:
resos = [ ]
if any ( [ i in list ( "RKSYWM" ) for i in seq ] ) :
    for base in list ( "RKSYWM" ) :
        if base in seq :
            resos . append ( seq . replace ( base , AMBIGS [ base ] [ 0 ] ) )
            resos . append ( seq . replace ( base , AMBIGS [ base ] [ 1 ] ) )
    return resos
else :
    return [ seq , "" ]


signature:
def splitalleles ( consensus ) :
    """takes diploid consensus alleles with phase data stored as a mixture
 of upper and lower case characters and splits it into 2 alleles"""

implementation:
# # store two alleles , allele1 will start with bigbase
allele1 = list ( consensus )
allele2 = list ( consensus )
hidx = [ i for ( i , j ) in enumerate ( consensus ) if j in "RKSWYMrkswym" ]
# # do remaining h sites
for idx in hidx :
    hsite = consensus [ idx ]
    if hsite . isupper ( ) :
        allele1 [ idx ] = PRIORITY [ hsite ]
        allele2 [ idx ] = MINOR [ hsite ]
    else :
        allele1 [ idx ] = MINOR [ hsite . upper ( ) ]
        allele2 [ idx ] = PRIORITY [ hsite . upper ( ) ]
# # convert back to strings
allele1 = "" . join ( allele1 )
allele2 = "" . join ( allele2 )
return allele1 , allele2


signature:
def comp ( seq ) :
    """returns a seq with complement . Preserves little n ' s for splitters ."""

implementation:
# # makes base to its small complement then makes upper
return seq . replace ( "A" , 't' ) . replace ( 'T' , 'a' ) . replace ( 'C' , 'g' ) . replace ( 'G' , 'c' ) . replace ( 'n' , 'Z' ) . upper ( ) . replace ( "Z" , "n" )


signature:
def fullcomp ( seq ) :
    """returns complement of sequence including ambiguity characters ,
 and saves lower case info for multiple hetero sequences"""

implementation:
# # this is surely not the most efficient . . .
seq = seq . replace ( "A" , 'u' ) . replace ( 'T' , 'v' ) . replace ( 'C' , 'p' ) . replace ( 'G' , 'z' ) . replace ( 'u' , 'T' ) . replace ( 'v' , 'A' ) . replace ( 'p' , 'G' ) . replace ( 'z' , 'C' )
# # No complement for S & W b / c complements are S & W , respectively
seq = seq . replace ( 'R' , 'u' ) . replace ( 'K' , 'v' ) . replace ( 'Y' , 'b' ) . replace ( 'M' , 'o' ) . replace ( 'u' , 'Y' ) . replace ( 'v' , 'M' ) . replace ( 'b' , 'R' ) . replace ( 'o' , 'K' )
seq = seq . replace ( 'r' , 'u' ) . replace ( 'k' , 'v' ) . replace ( 'y' , 'b' ) . replace ( 'm' , 'o' ) . replace ( 'u' , 'y' ) . replace ( 'v' , 'm' ) . replace ( 'b' , 'r' ) . replace ( 'o' , 'k' )
return seq


signature:
def fastq_touchup_for_vsearch_merge ( read , outfile , reverse = False ) :
    """option to change orientation of reads and sets Qscore to B"""

implementation:
counts = 0
with open ( outfile , 'w' ) as out : # # read in paired end read files 4 lines at a time
    if read . endswith ( ".gz" ) :
        fr1 = gzip . open ( read , 'rb' )
    else :
        fr1 = open ( read , 'rb' )
    quarts = itertools . izip ( * [ iter ( fr1 ) ] * 4 )
    # # a list to store until writing
    writing = [ ]
    while 1 :
        try :
            lines = quarts . next ( )
        except StopIteration :
            break
        if reverse :
            seq = lines [ 1 ] . strip ( ) [ : : - 1 ]
        else :
            seq = lines [ 1 ] . strip ( )
        writing . append ( "" . join ( [ lines [ 0 ] , seq + "\n" , lines [ 2 ] , "B" * len ( seq ) ] ) )
        # # write to disk
        counts += 1
        if not counts % 1000 :
            out . write ( "\n" . join ( writing ) + "\n" )
            writing = [ ]
    if writing :
        out . write ( "\n" . join ( writing ) )
out . close ( )
fr1 . close ( )


signature:
def merge_pairs_after_refmapping ( data , two_files , merged_out ) :
    """A function to merge fastq files produced by bam2fq ."""

implementation:
# # create temp files
nonmerged1 = tempfile . NamedTemporaryFile ( mode = 'wb' , dir = data . dirs . edits , suffix = "_nonmerged_R1_.fastq" ) . name
nonmerged2 = tempfile . NamedTemporaryFile ( mode = 'wb' , dir = data . dirs . edits , suffix = "_nonmerged_R2_.fastq" ) . name
# # get the maxn and minlen values
minlen = str ( max ( 32 , data . paramsdict [ "filter_min_trim_len" ] ) )
try :
    maxn = sum ( data . paramsdict [ 'max_low_qual_bases' ] )
except TypeError :
    maxn = data . paramsdict [ 'max_low_qual_bases' ]
# # set the quality scores abritrarily high and orient R2 correctly
tmp1 = two_files [ 0 ] [ 0 ]
tmp2 = two_files [ 0 ] [ 1 ]
fastq_touchup_for_vsearch_merge ( tmp1 , tmp1 + ".tu" , False )
fastq_touchup_for_vsearch_merge ( tmp2 , tmp2 + ".tu" , True )
# # command string to call vsearch
cmd = [ ipyrad . bins . vsearch , "--fastq_mergepairs" , tmp1 + ".tu" , "--reverse" , tmp2 + ".tu" , "--fastqout" , merged_out , "--fastqout_notmerged_fwd" , nonmerged1 , "--fastqout_notmerged_rev" , nonmerged2 , "--fasta_width" , "0" , "--fastq_minmergelen" , minlen , "--fastq_maxns" , str ( maxn ) , "--fastq_minovlen" , "10" , "--fastq_maxdiffs" , "4" , "--label_suffix" , "_m1" , "--fastq_qmax" , "1000" , "--threads" , "2" , "--fastq_allowmergestagger" ]
# # run vsearch but allow kbd
proc = sps . Popen ( cmd , stderr = sps . STDOUT , stdout = sps . PIPE )
try :
    res = proc . communicate ( ) [ 0 ]
except KeyboardInterrupt :
    proc . kill ( )
# # cleanup tmp files if job failed or stopped
if proc . returncode :
    LOGGER . error ( "Error: %s %s" , cmd , res )
    raise IPyradWarningExit ( "Error merge pairs:\n %s\n%s" , cmd , res )
# # record how many read pairs were merged
with open ( merged_out , 'r' ) as tmpf :
    nmerged = sum ( 1 for i in tmpf . readlines ( ) ) // 4
# # Concat unmerged pairs with a ' nnnn ' separator
with open ( merged_out , 'ab' ) as combout : # # read in paired end read files 4 lines at a time
    fr1 = open ( nonmerged1 , 'rb' )
    quart1 = itertools . izip ( * [ iter ( fr1 ) ] * 4 )
    fr2 = open ( nonmerged2 , 'rb' )
    quart2 = itertools . izip ( * [ iter ( fr2 ) ] * 4 )
    quarts = itertools . izip ( quart1 , quart2 )
    # # a list to store until writing
    writing = [ ]
    counts = 0
    # # iterate until done
    while 1 :
        try :
            read1s , read2s = quarts . next ( )
        except StopIteration :
            break
        # # store the read
        writing . append ( "" . join ( [ read1s [ 0 ] , read1s [ 1 ] . strip ( ) + "nnnn" + read2s [ 1 ] , read1s [ 2 ] , read1s [ 3 ] . strip ( ) + "nnnn" + read2s [ 3 ] , ] ) )
        # # count up until time to write
        counts += 1
        if not counts % 10 :
            combout . write ( "" . join ( writing ) )
            writing = [ ]
    # # write the remaining
    if writing :
        combout . write ( "" . join ( writing ) )
    # # close handles
    fr1 . close ( )
    fr2 . close ( )
    combout . close ( )
# # remove temp files ( or do this later )
rmfiles = [ nonmerged1 , nonmerged2 , tmp1 , tmp2 , tmp1 + ".tu" , tmp2 + ".tu" ]
for rmfile in rmfiles :
    if os . path . exists ( rmfile ) :
        os . remove ( rmfile )
return nmerged


signature:
def merge_after_pysam ( data , clust ) :
    """This is for pysam post - flight merging . The input is a cluster
 for an individual locus . We have to split the clusters , write
 R1 and R2 to files then call merge _ pairs ( ) . This is not ideal ,
 it ' s slow , but it works . This is the absolute worst way to do this ,
 it bounces all the files for each locus off the disk . I / O _ hog _ ."""

implementation:
try :
    r1file = tempfile . NamedTemporaryFile ( mode = 'wb' , delete = False , dir = data . dirs . edits , suffix = "_R1_.fastq" )
    r2file = tempfile . NamedTemporaryFile ( mode = 'wb' , delete = False , dir = data . dirs . edits , suffix = "_R2_.fastq" )
    r1dat = [ ]
    r2dat = [ ]
    for locus in clust :
        sname , seq = locus . split ( "\n" )
        # # Have to munge the sname to make it look like fastq format
        sname = "@" + sname [ 1 : ]
        r1 , r2 = seq . split ( "nnnn" )
        r1dat . append ( "{}\n{}\n{}\n{}" . format ( sname , r1 , "+" , "B" * ( len ( r1 ) ) ) )
        r2dat . append ( "{}\n{}\n{}\n{}" . format ( sname , r2 , "+" , "B" * ( len ( r2 ) ) ) )
    r1file . write ( "\n" . join ( r1dat ) )
    r2file . write ( "\n" . join ( r2dat ) )
    r1file . close ( )
    r2file . close ( )
    # # Read in the merged data and format to return as a clust
    merged_file = tempfile . NamedTemporaryFile ( mode = 'wb' , dir = data . dirs . edits , suffix = "_merged.fastq" ) . name
    clust = [ ]
    merge_pairs ( data , [ ( r1file . name , r2file . name ) ] , merged_file , 0 , 1 )
    with open ( merged_file ) as infile :
        quarts = itertools . izip ( * [ iter ( infile ) ] * 4 )
        while 1 :
            try :
                sname , seq , _ , _ = quarts . next ( )
                # # Vsearch expects R2 oriented how it would be in a raw data file
                # # i . e . revcomp , and that ' s also how it returns it
                # # but we want to maintain the genomic orientation so R1 and
                # # R2 are both on the + strand and both in ascending positional order
                # # so here if the reads don ' t merge we have to revcomp R2
                if not "_m1" in sname . rsplit ( ";" , 1 ) [ 1 ] :
                    try :
                        R1 , R2 = seq . split ( "nnnn" )
                        seq = R1 + "nnnn" + revcomp ( R2 )
                    except ValueError as inst :
                        LOGGER . error ( "Failed merge_after_pysam: {} {}" . format ( sname , seq ) )
                        raise
            except StopIteration :
                break
            # # put sname back
            sname = ">" + sname [ 1 : ]
            clust . extend ( [ sname . strip ( ) , seq . strip ( ) ] )
except :
    LOGGER . info ( "Error in merge_pairs post-refmap." )
    raise
finally :
    for i in [ r1file . name , r2file . name , merged_file ] :
        if os . path . exists ( i ) :
            log_level = logging . getLevelName ( LOGGER . getEffectiveLevel ( ) )
            # # uncomment this to shit ALL over the filesystem .
            # # if not log _ level = = " DEBUG " :
            os . remove ( i )
return clust


signature:
def merge_pairs ( data , two_files , merged_out , revcomp , merge ) :
    """Merge PE reads . Takes in a list of unmerged files [ r1 , r2 ] and the
 filehandle to write merged data to , and it returns the number of reads
 that were merged ( overlapping ) . If merge = = 0 then only concat pairs ( nnnn ) ,
 no merging in vsearch .
 Parameters
 two _ files ( tuple ) :
 A list or tuple of the [ r1 , r2 ] files to be merged .
 merged _ out ( str ) :
 A string file handle for the merged data to be written to .
 revcomp ( bool ) :
 Whether or not to revcomp the R2s .
 merge ( bool ) :
 Whether or not to perform vsearch merging . If not then reads are simply
 concatenated with a ' nnnn ' separator .
 Returns
 If merge is on then the func will return the number of pairs
 successfully merged , else it returns - 1."""

implementation:
LOGGER . debug ( "Entering merge_pairs()" )
# # Return the number of merged pairs
nmerged = - 1
# # Check input files from inside list - tuple [ ( r1 , r2 ) ]
for fhandle in two_files [ 0 ] :
    if not os . path . exists ( fhandle ) :
        raise IPyradWarningExit ( """
    Attempting to merge a file that doesn't exist - {}""" . format ( fhandle ) )
# # If it already exists , clean up the old merged file
if os . path . exists ( merged_out ) :
    os . remove ( merged_out )
# # if merge then catch nonmerged in a separate file
if merge :
    nonmerged1 = tempfile . NamedTemporaryFile ( mode = 'wb' , dir = data . dirs . edits , suffix = "_nonmerged_R1_.fastq" ) . name
    nonmerged2 = tempfile . NamedTemporaryFile ( mode = 'wb' , dir = data . dirs . edits , suffix = "_nonmerged_R2_.fastq" ) . name
# # if not merging then the nonmerged reads will come from the normal edits
else :
    nonmerged1 = two_files [ 0 ] [ 0 ]
    nonmerged2 = two_files [ 0 ] [ 1 ]
# # get the maxn and minlen values
try :
    maxn = sum ( data . paramsdict [ 'max_low_qual_bases' ] )
except TypeError :
    maxn = data . paramsdict [ 'max_low_qual_bases' ]
minlen = str ( max ( 32 , data . paramsdict [ "filter_min_trim_len" ] ) )
# # we need to gunzip the files if they are zipped ( at least for now )
if merge and two_files [ 0 ] [ 0 ] . endswith ( ".gz" ) :
    LOGGER . info ( "gunzipping pairs" )
    tmp1 = os . path . splitext ( two_files [ 0 ] [ 0 ] ) [ 0 ] + ".tmp1"
    tmp2 = os . path . splitext ( two_files [ 0 ] [ 1 ] ) [ 0 ] + ".tmp2"
    out1 = open ( tmp1 , 'w' )
    out2 = open ( tmp2 , 'w' )
    gun1 = sps . Popen ( [ "gunzip" , "-c" , two_files [ 0 ] [ 0 ] ] , stderr = sps . STDOUT , stdout = out1 , close_fds = True )
    gun2 = sps . Popen ( [ "gunzip" , "-c" , two_files [ 0 ] [ 1 ] ] , stderr = sps . STDOUT , stdout = out2 , close_fds = True )
    _ = gun1 . communicate ( )
    _ = gun2 . communicate ( )
    out1 . close ( )
    out2 . close ( )
else :
    tmp1 = two_files [ 0 ] [ 0 ]
    tmp2 = two_files [ 0 ] [ 1 ]
try : # # If we are actually mergeing and not just joining then do vsearch
    if merge : # # create tmp files with high quality scores and with R2 oriented
        cmd = [ ipyrad . bins . vsearch , "--fastq_mergepairs" , tmp1 , "--reverse" , tmp2 , "--fastqout" , merged_out , "--fastqout_notmerged_fwd" , nonmerged1 , "--fastqout_notmerged_rev" , nonmerged2 , "--fasta_width" , "0" , "--fastq_minmergelen" , minlen , "--fastq_maxns" , str ( maxn ) , "--fastq_minovlen" , "20" , "--fastq_maxdiffs" , "4" , "--label_suffix" , "_m1" , "--fastq_qmax" , "1000" , "--threads" , "2" , "--fastq_allowmergestagger" ]
        LOGGER . debug ( "merge cmd: %s" , " " . join ( cmd ) )
        proc = sps . Popen ( cmd , stderr = sps . STDOUT , stdout = sps . PIPE )
        try :
            res = proc . communicate ( ) [ 0 ]
        except KeyboardInterrupt :
            proc . kill ( )
        if proc . returncode :
            LOGGER . error ( "Error: %s %s" , cmd , res )
            # # remove temp files
            rmfiles = [ os . path . splitext ( two_files [ 0 ] [ 0 ] ) [ 0 ] + ".tmp1" , os . path . splitext ( two_files [ 0 ] [ 1 ] ) [ 0 ] + ".tmp2" , nonmerged1 , nonmerged2 ]
            for rmfile in rmfiles :
                if os . path . exists ( rmfile ) :
                    os . remove ( rmfile )
            raise IPyradWarningExit ( "Error merge pairs:\n %s\n%s" , cmd , res )
        # # record how many read pairs were merged
        with open ( merged_out , 'r' ) as tmpf : # nmerged = len ( tmpf . readlines ( ) ) / / 4
            nmerged = sum ( 1 for i in tmpf . readlines ( ) ) // 4
    # # Combine the unmerged pairs and append to the merge file
    with open ( merged_out , 'ab' ) as combout : # # read in paired end read files 4 lines at a time
        if nonmerged1 . endswith ( ".gz" ) :
            fr1 = gzip . open ( nonmerged1 , 'rb' )
        else :
            fr1 = open ( nonmerged1 , 'rb' )
        quart1 = itertools . izip ( * [ iter ( fr1 ) ] * 4 )
        if nonmerged2 . endswith ( ".gz" ) :
            fr2 = gzip . open ( nonmerged2 , 'rb' )
        else :
            fr2 = open ( nonmerged2 , 'rb' )
        quart2 = itertools . izip ( * [ iter ( fr2 ) ] * 4 )
        quarts = itertools . izip ( quart1 , quart2 )
        # # a list to store until writing
        writing = [ ]
        counts = 0
        # # iterate until done
        while 1 :
            try :
                read1s , read2s = quarts . next ( )
            except StopIteration :
                break
            if revcomp :
                writing . append ( "" . join ( [ read1s [ 0 ] , read1s [ 1 ] . strip ( ) + "nnnn" + comp ( read2s [ 1 ] . strip ( ) [ : : - 1 ] ) + "\n" , read1s [ 2 ] , read1s [ 3 ] . strip ( ) + "nnnn" + read2s [ 3 ] . strip ( ) [ : : - 1 ] + "\n" , ] ) )
            else :
                writing . append ( "" . join ( [ read1s [ 0 ] , read1s [ 1 ] . strip ( ) + "nnnn" + read2s [ 1 ] , read1s [ 2 ] , read1s [ 3 ] . strip ( ) + "nnnn" + read2s [ 3 ] , ] ) )
            counts += 1
            if not counts % 10 :
                combout . write ( "" . join ( writing ) )
                writing = [ ]
        if writing :
            combout . write ( "" . join ( writing ) )
    # # close handles
    fr1 . close ( )
    fr2 . close ( )
    combout . close ( )
except Exception as inst :
    LOGGER . error ( "Exception in merge_pairs - {}" . format ( inst ) )
    raise
# # No matter what happens please clean up the temp files .
finally : # # if merged then delete the nonmerge tmp files
    if merge : # # remove temp files
        rmfiles = [ nonmerged1 , nonmerged2 , os . path . splitext ( two_files [ 0 ] [ 0 ] ) [ 0 ] + ".tmp1" , os . path . splitext ( two_files [ 0 ] [ 1 ] ) [ 0 ] + ".tmp2" ]
        for rmfile in rmfiles :
            if os . path . exists ( rmfile ) :
                os . remove ( rmfile )
return nmerged


signature:
def revcomp ( sequence ) :
    "returns reverse complement of a string"

implementation:
sequence = sequence [ : : - 1 ] . strip ( ) . replace ( "A" , "t" ) . replace ( "T" , "a" ) . replace ( "C" , "g" ) . replace ( "G" , "c" ) . upper ( )
return sequence