signature:
def multicat ( data , samples , ipyclient ) :
    """Runs singlecat and cleanup jobs for each sample .
 For each sample this fills its own hdf5 array with catg data & indels .
 This is messy , could use simplifiying ."""

implementation:
# # progress ticker
start = time . time ( )
printstr = " indexing clusters     | {} | s6 |"
elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
progressbar ( 20 , 0 , printstr . format ( elapsed ) )
# # parallel client
lbview = ipyclient . load_balanced_view ( )
# # First submit a sleeper job as temp _ flag for cleanups
last_sample = 0
cleanups = { }
cleanups [ last_sample ] = lbview . apply ( time . sleep , 0.0 )
# # get samples and names , sorted
snames = [ i . name for i in samples ]
snames . sort ( )
# # Build an array for quickly indexing consens reads from catg files .
# # save as a npy int binary file .
uhandle = os . path . join ( data . dirs . across , data . name + ".utemp.sort" )
bseeds = os . path . join ( data . dirs . across , data . name + ".tmparrs.h5" )
# # send as first async1 job
async1 = lbview . apply ( get_seeds_and_hits , * ( uhandle , bseeds , snames ) )
async2 = lbview . apply ( fill_dups_arr , data )
# # progress bar for seed / hit sorting
while not ( async1 . ready ( ) and async2 . ready ( ) ) :
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( 20 , 0 , printstr . format ( elapsed ) )
    time . sleep ( 0.1 )
if not async1 . successful ( ) :
    raise IPyradWarningExit ( "error in get_seeds: %s" , async1 . exception ( ) )
if not async2 . successful ( ) :
    raise IPyradWarningExit ( "error in fill_dups: %s" , async2 . exception ( ) )
# # make a limited njobs view based on mem limits
# # is using smallview necessary ? ( yes , it is for bad libraries )
smallview = ipyclient . load_balanced_view ( targets = ipyclient . ids [ : : 2 ] )
# # make sure there are no old tmp . h5 files
smpios = [ os . path . join ( data . dirs . across , sample . name + '.tmp.h5' ) for sample in samples ]
for smpio in smpios :
    if os . path . exists ( smpio ) :
        os . remove ( smpio )
# # send ' singlecat ( ) ' jobs to engines
jobs = { }
for sample in samples :
    sidx = snames . index ( sample . name )
    jobs [ sample . name ] = smallview . apply ( singlecat , * ( data , sample , bseeds , sidx ) )
# # check for finished and submit disk - writing job when finished
alljobs = len ( jobs )
while 1 : # # check for finished jobs
    curkeys = jobs . keys ( )
    for key in curkeys :
        async = jobs [ key ]
        if async . ready ( ) :
            if async . successful ( ) : # # submit cleanup for finished job
                args = ( data , data . samples [ key ] , snames . index ( key ) )
                with lbview . temp_flags ( after = cleanups [ last_sample ] ) :
                    cleanups [ key ] = lbview . apply ( write_to_fullarr , * args )
                last_sample = key
                del jobs [ key ]
            else :
                err = jobs [ key ] . exception ( )
                errmsg = "singlecat error: {} {}" . format ( key , err )
                raise IPyradWarningExit ( errmsg )
    # # print progress or break
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( alljobs , alljobs - len ( jobs ) , printstr . format ( elapsed ) )
    time . sleep ( 0.1 )
    if not jobs :
        break
# # add the dask _ chroms func for reference data
if 'reference' in data . paramsdict [ "assembly_method" ] :
    with lbview . temp_flags ( after = cleanups . values ( ) ) :
        cleanups [ 'ref' ] = lbview . apply ( dask_chroms , * ( data , samples ) )
# # wait for " write _ to _ fullarr " jobs to finish
print ( "" )
start = time . time ( )
printstr = " building database     | {} | s6 |"
while 1 :
    finished = [ i for i in cleanups . values ( ) if i . ready ( ) ]
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( len ( cleanups ) , len ( finished ) , printstr . format ( elapsed ) )
    time . sleep ( 0.1 )
    # # break if one failed , or if finished
    if not all ( [ i . successful ( ) for i in finished ] ) :
        break
    if len ( cleanups ) == len ( finished ) :
        break
# # check for errors
for job in cleanups :
    if cleanups [ job ] . ready ( ) :
        if not cleanups [ job ] . successful ( ) :
            err = " error in write_to_fullarr ({}) {}" . format ( job , cleanups [ job ] . result ( ) )
            LOGGER . error ( err )
            raise IPyradWarningExit ( err )
# # remove large indels array file and singlecat tmparr file
ifile = os . path . join ( data . dirs . across , data . name + ".tmp.indels.hdf5" )
if os . path . exists ( ifile ) :
    os . remove ( ifile )
if os . path . exists ( bseeds ) :
    os . remove ( bseeds )
for sh5 in [ os . path . join ( data . dirs . across , i . name + ".tmp.h5" ) for i in samples ] :
    os . remove ( sh5 )
# # print final progress
elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
progressbar ( 10 , 10 , printstr . format ( elapsed ) )
print ( "" )


signature:
def singlecat ( data , sample , bseeds , sidx , nloci ) :
    """Orders catg data for each sample into the final locus order . This allows
 all of the individual catgs to simply be combined later . They are also in
 the same order as the indels array , so indels are inserted from the indel
 array that is passed in ."""

implementation:
LOGGER . info ( "in single cat here" )
# # enter ref data ?
isref = 'reference' in data . paramsdict [ "assembly_method" ]
# # grab seeds and hits info for this sample
with h5py . File ( bseeds , 'r' ) as io5 : # # get hits just for this sample and sort them by sample order index
    hits = io5 [ "uarr" ] [ : ]
    hits = hits [ hits [ : , 1 ] == sidx , : ]
    # hits = hits [ hits [ : , 2 ] . argsort ( ) ]
    # # get seeds just for this sample and sort them by sample order index
    seeds = io5 [ "seedsarr" ] [ : ]
    seeds = seeds [ seeds [ : , 1 ] == sidx , : ]
    # seeds = seeds [ seeds [ : , 2 ] . argsort ( ) ]
    full = np . concatenate ( ( seeds , hits ) )
    full = full [ full [ : , 0 ] . argsort ( ) ]
# # still using max + 20 len limit , rare longer merged reads get trimmed
# # we need to allow room for indels to be added too
maxlen = data . _hackersonly [ "max_fragment_length" ] + 20
# # we ' ll fill a new catg and alleles arr for this sample in locus order ,
# # which is known from seeds and hits
ocatg = np . zeros ( ( nloci , maxlen , 4 ) , dtype = np . uint32 )
onall = np . zeros ( nloci , dtype = np . uint8 )
ochrom = np . zeros ( ( nloci , 3 ) , dtype = np . int64 )
# # grab the sample ' s data and write to ocatg and onall
if not sample . files . database :
    raise IPyradWarningExit ( "missing catg file - {}" . format ( sample . name ) )
with h5py . File ( sample . files . database , 'r' ) as io5 : # # get it and delete it
    catarr = io5 [ "catg" ] [ : ]
    tmp = catarr [ full [ : , 2 ] , : maxlen , : ]
    del catarr
    ocatg [ full [ : , 0 ] , : tmp . shape [ 1 ] , : ] = tmp
    del tmp
    # # get it and delete it
    nall = io5 [ "nalleles" ] [ : ]
    onall [ full [ : , 0 ] ] = nall [ full [ : , 2 ] ]
    del nall
    # # fill the reference data
    if isref :
        chrom = io5 [ "chroms" ] [ : ]
        ochrom [ full [ : , 0 ] ] = chrom [ full [ : , 2 ] ]
        del chrom
# # get indel locations for this sample
ipath = os . path . join ( data . dirs . across , data . name + ".tmp.indels.hdf5" )
with h5py . File ( ipath , 'r' ) as ih5 :
    indels = ih5 [ "indels" ] [ sidx , : , : maxlen ]
# # insert indels into ocatg
newcatg = inserted_indels ( indels , ocatg )
del ocatg , indels
# # save individual tmp h5 data
smpio = os . path . join ( data . dirs . across , sample . name + '.tmp.h5' )
with h5py . File ( smpio , 'w' ) as oh5 :
    oh5 . create_dataset ( "icatg" , data = newcatg , dtype = np . uint32 )
    oh5 . create_dataset ( "inall" , data = onall , dtype = np . uint8 )
    if isref :
        oh5 . create_dataset ( "ichrom" , data = ochrom , dtype = np . int64 )


signature:
def write_to_fullarr ( data , sample , sidx ) :
    """writes arrays to h5 disk"""

implementation:
# # enter ref data ?
# isref = ' reference ' in data . paramsdict [ " assembly _ method " ]
LOGGER . info ( "writing fullarr %s %s" , sample . name , sidx )
# # save big arrays to disk temporarily
with h5py . File ( data . clust_database , 'r+' ) as io5 : # # open views into the arrays we plan to fill
    chunk = io5 [ "catgs" ] . attrs [ "chunksize" ] [ 0 ]
    catg = io5 [ "catgs" ]
    nall = io5 [ "nalleles" ]
    # # adding an axis to newcatg makes it write about 1000X faster .
    smpio = os . path . join ( data . dirs . across , sample . name + '.tmp.h5' )
    with h5py . File ( smpio ) as indat : # # grab all of the data from this sample ' s arrays
        newcatg = indat [ "icatg" ]
        onall = indat [ "inall" ]
        # # enter it into the full array one chunk at a time
        for cidx in xrange ( 0 , catg . shape [ 0 ] , chunk ) :
            end = cidx + chunk
            catg [ cidx : end , sidx : sidx + 1 , : ] = np . expand_dims ( newcatg [ cidx : end , : ] , axis = 1 )
            nall [ : , sidx : sidx + 1 ] = np . expand_dims ( onall , axis = 1 )


signature:
def dask_chroms ( data , samples ) :
    """A dask relay function to fill chroms for all samples"""

implementation:
# # example concatenating with dask
h5s = [ os . path . join ( data . dirs . across , s . name + ".tmp.h5" ) for s in samples ]
handles = [ h5py . File ( i ) for i in h5s ]
dsets = [ i [ '/ichrom' ] for i in handles ]
arrays = [ da . from_array ( dset , chunks = ( 10000 , 3 ) ) for dset in dsets ]
stack = da . stack ( arrays , axis = 2 )
# # max chrom ( should we check for variable hits ? if so , things can get wonk )
maxchrom = da . max ( stack , axis = 2 ) [ : , 0 ]
# # max pos
maxpos = da . max ( stack , axis = 2 ) [ : , 2 ]
# # min pos
mask = stack == 0
stack [ mask ] = 9223372036854775807
# # max int64 value
minpos = da . min ( stack , axis = 2 ) [ : , 1 ]
final = da . stack ( [ maxchrom , minpos , maxpos ] , axis = 1 )
final . to_hdf5 ( data . clust_database , "/chroms" )
# # close the h5 handles
_ = [ i . close ( ) for i in handles ]


signature:
def inserted_indels ( indels , ocatg ) :
    """inserts indels into the catg array"""

implementation:
# # return copy with indels inserted
newcatg = np . zeros ( ocatg . shape , dtype = np . uint32 )
# # iterate over loci and make extensions for indels
for iloc in xrange ( ocatg . shape [ 0 ] ) : # # get indels indices
    indidx = np . where ( indels [ iloc , : ] ) [ 0 ]
    if np . any ( indidx ) : # # which new ( empty ) rows will be added
        allrows = np . arange ( ocatg . shape [ 1 ] )
        mask = np . ones ( allrows . shape [ 0 ] , dtype = np . bool_ )
        for idx in indidx :
            mask [ idx ] = False
        not_idx = allrows [ mask == 1 ]
        # # fill in new data into all other spots
        newcatg [ iloc ] [ not_idx ] = ocatg [ iloc , : not_idx . shape [ 0 ] ]
    else :
        newcatg [ iloc ] = ocatg [ iloc ]
return newcatg


signature:
def fill_superseqs ( data , samples ) :
    """Fills the superseqs array with seq data from cat . clust
 and fill the edges array with information about paired split locations ."""

implementation:
# # load super to get edges
io5 = h5py . File ( data . clust_database , 'r+' )
superseqs = io5 [ "seqs" ]
splits = io5 [ "splits" ]
# # samples are already sorted
snames = [ i . name for i in samples ]
LOGGER . info ( "snames %s" , snames )
# # get maxlen again
maxlen = data . _hackersonly [ "max_fragment_length" ] + 20
LOGGER . info ( "maxlen inside fill_superseqs is %s" , maxlen )
# # data has to be entered in blocks
infile = os . path . join ( data . dirs . across , data . name + "_catclust.gz" )
clusters = gzip . open ( infile , 'r' )
pairdealer = itertools . izip ( * [ iter ( clusters ) ] * 2 )
# # iterate over clusters
chunks = superseqs . attrs [ "chunksize" ]
chunksize = chunks [ 0 ]
done = 0
iloc = 0
cloc = 0
chunkseqs = np . zeros ( chunks , dtype = "|S1" )
chunkedge = np . zeros ( chunksize , dtype = np . uint16 )
while 1 :
    try :
        done , chunk = clustdealer ( pairdealer , 1 )
    except IndexError :
        raise IPyradWarningExit ( "clustfile formatting error in %s" , chunk )
    # # if chunk is full put into superseqs and reset counter
    if cloc == chunksize :
        LOGGER . info ( "cloc chunk writing %s" , cloc )
        superseqs [ iloc - cloc : iloc ] = chunkseqs
        splits [ iloc - cloc : iloc ] = chunkedge
        # # reset chunkseqs , chunkedge , cloc
        cloc = 0
        chunkseqs = np . zeros ( ( chunksize , len ( samples ) , maxlen ) , dtype = "|S1" )
        chunkedge = np . zeros ( ( chunksize ) , dtype = np . uint16 )
    # # get seq and split it
    if chunk :
        try :
            fill = np . zeros ( ( len ( samples ) , maxlen ) , dtype = "|S1" )
            fill . fill ( "N" )
            piece = chunk [ 0 ] . strip ( ) . split ( "\n" )
            names = piece [ 0 : : 2 ]
            seqs = np . array ( [ list ( i ) for i in piece [ 1 : : 2 ] ] )
            # # fill in the separator if it exists
            separator = np . where ( np . all ( seqs == 'n' , axis = 0 ) ) [ 0 ]
            if np . any ( separator ) :
                chunkedge [ cloc ] = separator . min ( )
            # # fill in the hits
            # # seqs will be ( 5 , ) IF the seqs are variable lengths , which
            # # can happen if it had duplicaes AND there were indels , and
            # # so the indels did not get aligned
            try :
                shlen = seqs . shape [ 1 ]
            except IndexError as inst :
                shlen = min ( [ len ( x ) for x in seqs ] )
            for name , seq in zip ( names , seqs ) :
                sidx = snames . index ( name . rsplit ( "_" , 1 ) [ 0 ] )
                # fill [ sidx , : shlen ] = seq [ : maxlen ]
                fill [ sidx , : shlen ] = seq [ : shlen ]
            # # PUT seqs INTO local ARRAY
            chunkseqs [ cloc ] = fill
        except Exception as inst :
            LOGGER . info ( inst )
            LOGGER . info ( "\nfill: %s\nshlen %s\nmaxlen %s" , fill . shape , shlen , maxlen )
            LOGGER . info ( "dupe chunk \n{}" . format ( "\n" . join ( chunk ) ) )
        # # increase counters if there was a chunk
        cloc += 1
        iloc += 1
    if done :
        break
# # write final leftover chunk
superseqs [ iloc - cloc : , ] = chunkseqs [ : cloc ]
splits [ iloc - cloc : ] = chunkedge [ : cloc ]
# # close super
io5 . close ( )
clusters . close ( )
# # edges is filled with splits for paired data .
LOGGER . info ( "done filling superseqs" )


signature:
def count_seeds ( usort ) :
    """uses bash commands to quickly count N seeds from utemp file"""

implementation:
with open ( usort , 'r' ) as insort :
    cmd1 = [ "cut" , "-f" , "2" ]
    cmd2 = [ "uniq" ]
    cmd3 = [ "wc" ]
    proc1 = sps . Popen ( cmd1 , stdin = insort , stdout = sps . PIPE , close_fds = True )
    proc2 = sps . Popen ( cmd2 , stdin = proc1 . stdout , stdout = sps . PIPE , close_fds = True )
    proc3 = sps . Popen ( cmd3 , stdin = proc2 . stdout , stdout = sps . PIPE , close_fds = True )
    res = proc3 . communicate ( )
    nseeds = int ( res [ 0 ] . split ( ) [ 0 ] )
    proc1 . stdout . close ( )
    proc2 . stdout . close ( )
    proc3 . stdout . close ( )
return nseeds


signature:
def sort_seeds ( uhandle , usort ) :
    """sort seeds from cluster results"""

implementation:
cmd = [ "sort" , "-k" , "2" , uhandle , "-o" , usort ]
proc = sps . Popen ( cmd , close_fds = True )
proc . communicate ( )


signature:
def build_clustbits ( data , ipyclient , force ) :
    """Reconstitutes clusters from . utemp and htemp files and writes them
 to chunked files for aligning in muscle ."""

implementation:
# # If you run this step then we clear all tmp . fa and . indel . h5 files
if os . path . exists ( data . tmpdir ) :
    shutil . rmtree ( data . tmpdir )
    os . mkdir ( data . tmpdir )
# # parallel client
lbview = ipyclient . load_balanced_view ( )
start = time . time ( )
printstr = " building clusters     | {} | s6 |"
elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
progressbar ( 3 , 0 , printstr . format ( elapsed ) , spacer = data . _spacer )
uhandle = os . path . join ( data . dirs . across , data . name + ".utemp" )
usort = os . path . join ( data . dirs . across , data . name + ".utemp.sort" )
async1 = ""
# # skip usorting if not force and already exists
if not os . path . exists ( usort ) or force : # # send sort job to engines . Sorted seeds allows us to work through
# # the utemp file one locus at a time instead of reading all into mem .
    LOGGER . info ( "building reads file -- loading utemp file into mem" )
    async1 = lbview . apply ( sort_seeds , * ( uhandle , usort ) )
    while 1 :
        elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
        progressbar ( 3 , 0 , printstr . format ( elapsed ) , spacer = data . _spacer )
        if async1 . ready ( ) :
            break
        else :
            time . sleep ( 0.1 )
# # send count seeds job to engines .
async2 = lbview . apply ( count_seeds , usort )
while 1 :
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( 3 , 1 , printstr . format ( elapsed ) , spacer = data . _spacer )
    if async2 . ready ( ) :
        break
    else :
        time . sleep ( 0.1 )
# # wait for both to finish while printing progress timer
nseeds = async2 . result ( )
# # send the clust bit building job to work and track progress
async3 = lbview . apply ( sub_build_clustbits , * ( data , usort , nseeds ) )
while 1 :
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( 3 , 2 , printstr . format ( elapsed ) , spacer = data . _spacer )
    if async3 . ready ( ) :
        break
    else :
        time . sleep ( 0.1 )
elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
progressbar ( 3 , 3 , printstr . format ( elapsed ) , spacer = data . _spacer )
print ( "" )
# # check for errors
for job in [ async1 , async2 , async3 ] :
    try :
        if not job . successful ( ) :
            raise IPyradWarningExit ( job . result ( ) )
    except AttributeError : # # If we skip usorting then async1 = = " " so the call to
    # # successful ( ) raises , but we can ignore it .
        pass


signature:
def sub_build_clustbits ( data , usort , nseeds ) :
    """A subfunction of build _ clustbits to allow progress tracking . This func
 splits the unaligned clusters into bits for aligning on separate cores ."""

implementation:
# # load FULL concat fasta file into a dict . This could cause RAM issues .
# # this file has iupac codes in it , not ambigs resolved , and is gzipped .
LOGGER . info ( "loading full _catcons file into memory" )
allcons = { }
conshandle = os . path . join ( data . dirs . across , data . name + "_catcons.tmp" )
with gzip . open ( conshandle , 'rb' ) as iocons :
    cons = itertools . izip ( * [ iter ( iocons ) ] * 2 )
    for namestr , seq in cons :
        nnn , sss = [ i . strip ( ) for i in namestr , seq ]
        allcons [ nnn [ 1 : ] ] = sss
# # set optim to approximately 4 chunks per core . Smaller allows for a bit
# # cleaner looking progress bar . 40 cores will make 160 files .
optim = ( ( nseeds // ( data . cpus * 4 ) ) + ( nseeds % ( data . cpus * 4 ) ) )
LOGGER . info ( "building clustbits, optim=%s, nseeds=%s, cpus=%s" , optim , nseeds , data . cpus )
# # iterate through usort grabbing seeds and matches
with open ( usort , 'rb' ) as insort : # # iterator , seed null , and seqlist null
    isort = iter ( insort )
    loci = 0
    lastseed = 0
    fseqs = [ ]
    seqlist = [ ]
    seqsize = 0
    while 1 : # # grab the next line
        try :
            hit , seed , ori = isort . next ( ) . strip ( ) . split ( )
        except StopIteration :
            break
        try : # # if same seed , append match
            if seed != lastseed : # # store the last fseq , count it , and clear it
                if fseqs :
                    seqlist . append ( "\n" . join ( fseqs ) )
                    seqsize += 1
                    fseqs = [ ]
                # # occasionally write to file
                if seqsize >= optim :
                    if seqlist :
                        loci += seqsize
                        with open ( os . path . join ( data . tmpdir , data . name + ".chunk_{}" . format ( loci ) ) , 'w' ) as clustsout :
                            LOGGER . debug ( "writing chunk - seqsize {} loci {} {}" . format ( seqsize , loci , clustsout . name ) )
                            clustsout . write ( "\n//\n//\n" . join ( seqlist ) + "\n//\n//\n" )
                        # # reset list and counter
                        seqlist = [ ]
                        seqsize = 0
                # # store the new seed on top of fseq
                fseqs . append ( ">{}\n{}" . format ( seed , allcons [ seed ] ) )
                lastseed = seed
            # # add match to the seed
            seq = allcons [ hit ]
            # # revcomp if orientation is reversed
            if ori == "-" :
                seq = fullcomp ( seq ) [ : : - 1 ]
            fseqs . append ( ">{}\n{}" . format ( hit , seq ) )
        except KeyError as inst : # # Caught bad seed or hit ? Log and continue .
            LOGGER . error ( "Bad Seed/Hit: seqsize {}\tloci {}\tseed {}\thit {}" . format ( seqsize , loci , seed , hit ) )
# # write whatever is left over to the clusts file
if fseqs :
    seqlist . append ( "\n" . join ( fseqs ) )
    seqsize += 1
    loci += seqsize
if seqlist :
    with open ( os . path . join ( data . tmpdir , data . name + ".chunk_{}" . format ( loci ) ) , 'w' ) as clustsout :
        clustsout . write ( "\n//\n//\n" . join ( seqlist ) + "\n//\n//\n" )
# # final progress and cleanup
del allcons
clustbits = glob . glob ( os . path . join ( data . tmpdir , data . name + ".chunk_*" ) )
# # return stuff
return clustbits , loci


signature:
def build_input_file ( data , samples , randomseed ) :
    """[ This is run on an ipengine ]
 Make a concatenated consens file with sampled alleles ( no RSWYMK / rswymk ) .
 Orders reads by length and shuffles randomly within length classes"""

implementation:
# # get all of the consens handles for samples that have consens reads
# # this is better than using sample . files . consens for selecting files
# # b / c if they were moved we only have to edit data . dirs . consens
# # scratch the statement above , people shouldn ' t be moving files ,
# # they should be using merge / branch , and so sample . files . consens
# # is needed to keep track of samples from different dirs if they
# # are later merged into the same assembly .
# conshandles = [ os . path . join ( data . dirs . consens , sample . name + " . consens . gz " ) \
# for sample in samples if \
# sample . stats . reads _ consens ]
conshandles = [ sample . files . consens [ 0 ] for sample in samples if sample . stats . reads_consens ]
conshandles . sort ( )
assert conshandles , "no consensus files found"
# # concatenate all of the gzipped consens files
cmd = [ 'cat' ] + conshandles
# allcons = os . path . join ( data . dirs . consens , data . name + " _ catcons . tmp " )
allcons = os . path . join ( data . dirs . across , data . name + "_catcons.tmp" )
LOGGER . debug ( " " . join ( cmd ) )
with open ( allcons , 'w' ) as output :
    call = sps . Popen ( cmd , stdout = output , close_fds = True )
    call . communicate ( )
# # a string of sed substitutions for temporarily replacing hetero sites
# # skips lines with ' > ' , so it doesn ' t affect taxon names
subs = [ "/>/!s/W/A/g" , "/>/!s/w/A/g" , "/>/!s/R/A/g" , "/>/!s/r/A/g" , "/>/!s/M/A/g" , "/>/!s/m/A/g" , "/>/!s/K/T/g" , "/>/!s/k/T/g" , "/>/!s/S/C/g" , "/>/!s/s/C/g" , "/>/!s/Y/C/g" , "/>/!s/y/C/g" ]
subs = ";" . join ( subs )
# # impute pseudo - haplo information to avoid mismatch at hetero sites
# # the read data with hetero sites is put back into clustered data later .
# # pipe passed data from gunzip to sed .
cmd1 = [ "gunzip" , "-c" , allcons ]
cmd2 = [ "sed" , subs ]
LOGGER . debug ( " " . join ( cmd1 ) )
proc1 = sps . Popen ( cmd1 , stdout = sps . PIPE , close_fds = True )
allhaps = allcons . replace ( "_catcons.tmp" , "_cathaps.tmp" )
with open ( allhaps , 'w' ) as output :
    LOGGER . debug ( " " . join ( cmd2 ) )
    proc2 = sps . Popen ( cmd2 , stdin = proc1 . stdout , stdout = output , close_fds = True )
    proc2 . communicate ( )
proc1 . stdout . close ( )
# # now sort the file using vsearch
allsort = allcons . replace ( "_catcons.tmp" , "_catsort.tmp" )
cmd1 = [ ipyrad . bins . vsearch , "--sortbylength" , allhaps , "--fasta_width" , "0" , "--output" , allsort ]
LOGGER . debug ( " " . join ( cmd1 ) )
proc1 = sps . Popen ( cmd1 , close_fds = True )
proc1 . communicate ( )
# # shuffle sequences within size classes . Tested seed ( 8/31/2016)
# # shuffling works repeatably with seed .
random . seed ( randomseed )
# # open an iterator to lengthsorted file and grab two lines at at time
allshuf = allcons . replace ( "_catcons.tmp" , "_catshuf.tmp" )
outdat = open ( allshuf , 'w' )
indat = open ( allsort , 'r' )
idat = itertools . izip ( iter ( indat ) , iter ( indat ) )
done = 0
chunk = [ idat . next ( ) ]
while not done : # # grab 2 - lines until they become shorter ( unless there ' s only one )
    oldlen = len ( chunk [ - 1 ] [ - 1 ] )
    while 1 :
        try :
            dat = idat . next ( )
        except StopIteration :
            done = 1
            break
        if len ( dat [ - 1 ] ) == oldlen :
            chunk . append ( dat )
        else : # # send the last chunk off to be processed
            random . shuffle ( chunk )
            outdat . write ( "" . join ( itertools . chain ( * chunk ) ) )
            # # start new chunk
            chunk = [ dat ]
            break
# # do the last chunk
random . shuffle ( chunk )
outdat . write ( "" . join ( itertools . chain ( * chunk ) ) )
indat . close ( )
outdat . close ( )


signature:
def clean_and_build_concat ( data , samples , randomseed , ipyclient ) :
    """STEP 6-1:
 Clears dirs and databases and calls ' build _ input _ file ( ) '"""

implementation:
# # but check for new clust database name if this is a new branch
cleanup_tempfiles ( data )
catclust = os . path . join ( data . dirs . across , data . name + "_catclust.gz" )
if os . path . exists ( catclust ) :
    os . remove ( catclust )
if os . path . exists ( data . clust_database ) :
    os . remove ( data . clust_database )
# # get parallel view
start = time . time ( )
printstr = " concat/shuffle input  | {} | s6 |"
# # make a vsearch input fasta file with all samples reads concat
async = ipyclient [ 0 ] . apply ( build_input_file , * [ data , samples , randomseed ] )
while 1 :
    ready = int ( async . ready ( ) )
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( 1 , ready , printstr . format ( elapsed ) , spacer = data . _spacer )
    if ready :
        break
    else :
        time . sleep ( 0.1 )
print ( "" )
# # store that this step was successful
if not async . successful ( ) :
    raise IPyradWarningExit ( async . result ( ) )


signature:
def run ( data , samples , noreverse , force , randomseed , ipyclient , ** kwargs ) :
    """For step 6 the run function is sub divided a bit so that users with really
 difficult assemblies can possibly interrupt and restart the step from a
 checkpoint .
 Substeps that are run :
 1 . build concat consens file ,
 2 . cluster all consens ,
 3 . split clusters into bits ,
 4 . align bits ,
 5 . build indel array
 6 . build h5 array .
 7 . Enter seq data & cleanup"""

implementation:
# # if force then set checkpoint to zero and run all substeps for just
# # the user specified steps .
if force :
    data . _checkpoint = 0
    if kwargs . get ( 'substeps' ) :
        substeps = kwargs . get ( 'substeps' )
    else :
        substeps = range ( 1 , 8 )
# # if { data } . _ checkpoint attribute exists then find the checkpoint where
# # this assembly left off ( unless force ) and build step list from there .
else :
    if kwargs . get ( 'substeps' ) :
        substeps = kwargs . get ( 'substeps' )
    else :
        if hasattr ( data , '_checkpoint' ) :
            substeps = range ( max ( 1 , data . _checkpoint ) , 8 )
        else :
            data . _checkpoint = 0
            substeps = range ( 1 , 8 )
# # build substeps list to subset which funtions need to be run
if isinstance ( substeps , ( int , float , str ) ) :
    substeps = [ substeps ]
    substeps = [ int ( i ) for i in substeps ]
# # print continuation message
if substeps [ 0 ] != 1 :
    print ( "{}Continuing from checkpoint 6.{}" . format ( data . _spacer , substeps [ 0 ] ) )
LOGGER . info ( "checkpoint = %s" , data . _checkpoint )
LOGGER . info ( "substeps = %s" , substeps )
# # Set variables on data that are needed for all steps ;
data . dirs . across = os . path . realpath ( os . path . join ( data . paramsdict [ "project_dir" ] , data . name + "_across" ) )
data . tmpdir = os . path . join ( data . dirs . across , data . name + "-tmpalign" )
data . clust_database = os . path . join ( data . dirs . across , data . name + ".clust.hdf5" )
if not os . path . exists ( data . dirs . across ) :
    os . mkdir ( data . dirs . across )
if not os . path . exists ( data . tmpdir ) :
    os . mkdir ( data . tmpdir )
data . cpus = data . _ipcluster [ "cores" ]
if not data . cpus :
    data . cpus = len ( ipyclient )
# # STEP 6-1 : Clean database and build input concat file for clustering
if 1 in substeps :
    clean_and_build_concat ( data , samples , randomseed , ipyclient )
    data . _checkpoint = 1
# # STEP 6-2 : Cluster across w / vsearch ; uses all threads on largest host
if 2 in substeps :
    call_cluster ( data , noreverse , ipyclient )
    data . _checkpoint = 2
# # builds consens cluster bits and writes them to the tmp directory . These
# # will not be deleted until either step 6-6 is complete , or the force flag
# # is used . This will clear the tmpdir if it is run .
if 3 in substeps :
    build_clustbits ( data , ipyclient , force )
    data . _checkpoint = 3
# # muscle align the cluster bits and create tmp hdf5 indel arrays for the
# # next step . These will not be deleted until . . .
if 4 in substeps :
    multi_muscle_align ( data , samples , ipyclient )
    data . _checkpoint = 4
# # fill the indel array with the indel tmp arrays from aligning step .
if 5 in substeps :
    build_indels ( data , samples , ipyclient )
    data . _checkpoint = 5
if 6 in substeps : # # builds the final HDF5 array which includes three main keys
# # / catg - - contains all indiv catgs and has indels inserted
# # . attr [ ' samples ' ] = [ samples ]
# # / filters - - filled for dups , left empty for others until step 7.
# # . attr [ ' filters ' ] = [ f1 , f2 , f3 , f4 , f5]
# # / seqs - - contains the clustered sequence data as string arrays
# # . attr [ ' samples ' ] = [ samples ]
# # / edges - - gets the paired split locations for now .
# # / snps - - left empty for now
# # FILL SUPERCATG and fills dupfilter , indfilter , and nalleles
# # this function calls singlecat ( ) on each sample and enters their
# # resulting arrays into the superarray . If all singlecats are built
# # then it will continue to enter them into the database .
    LOGGER . info ( "multicat -- building full database" )
    new_multicat ( data , samples , ipyclient )
    data . _checkpoint = 6
if 7 in substeps : # # FILL SUPERSEQS and fills edges ( splits ) for paired - end data
    fill_superseqs ( data , samples )
    data . _checkpoint = 7
    # # remove files but not dir ( used in step 1 too )
    cleanup_tempfiles ( data )
    # # remove the tmpdir
    if os . path . exists ( data . tmpdir ) :
        shutil . rmtree ( data . tmpdir )
    # # set sample states
    for sample in samples :
        sample . stats . state = 6
    print ( "" )


signature:
def cleanup_tempfiles ( data ) :
    """Function to remove older files . This is called either in substep 1 or after
 the final substep so that tempfiles are retained for restarting interrupted
 jobs until we ' re sure they ' re no longer needed ."""

implementation:
# # remove align - related tmp files
tmps1 = glob . glob ( os . path . join ( data . tmpdir , "*.fa" ) )
tmps2 = glob . glob ( os . path . join ( data . tmpdir , "*.npy" ) )
for tmp in tmps1 + tmps2 :
    if os . path . exists ( tmp ) :
        os . remove ( tmp )
# # remove cluster related files
removal = [ os . path . join ( data . dirs . across , data . name + ".utemp" ) , os . path . join ( data . dirs . across , data . name + ".htemp" ) , os . path . join ( data . dirs . across , data . name + "_catcons.tmp" ) , os . path . join ( data . dirs . across , data . name + "_cathaps.tmp" ) , os . path . join ( data . dirs . across , data . name + "_catshuf.tmp" ) , os . path . join ( data . dirs . across , data . name + "_catsort.tmp" ) , os . path . join ( data . dirs . across , data . name + ".tmparrs.h5" ) , os . path . join ( data . dirs . across , data . name + ".tmp.indels.hdf5" ) , ]
for rfile in removal :
    if os . path . exists ( rfile ) :
        os . remove ( rfile )
# # remove singlecat related h5 files
smpios = glob . glob ( os . path . join ( data . dirs . across , '*.tmp.h5' ) )
for smpio in smpios :
    if os . path . exists ( smpio ) :
        os . remove ( smpio )


signature:
def assembly_cleanup ( data ) :
    """cleanup for assembly object"""

implementation:
# # build s2 results data frame
data . stats_dfs . s2 = data . _build_stat ( "s2" )
data . stats_files . s2 = os . path . join ( data . dirs . edits , 's2_rawedit_stats.txt' )
# # write stats for all samples
with io . open ( data . stats_files . s2 , 'w' , encoding = 'utf-8' ) as outfile :
    data . stats_dfs . s2 . fillna ( value = 0 ) . astype ( np . int ) . to_string ( outfile )


signature:
def parse_single_results ( data , sample , res1 ) :
    """parse results from cutadapt into sample data"""

implementation:
# # set default values
# sample . stats _ dfs . s2 [ " reads _ raw " ] = 0
sample . stats_dfs . s2 [ "trim_adapter_bp_read1" ] = 0
sample . stats_dfs . s2 [ "trim_quality_bp_read1" ] = 0
sample . stats_dfs . s2 [ "reads_filtered_by_Ns" ] = 0
sample . stats_dfs . s2 [ "reads_filtered_by_minlen" ] = 0
sample . stats_dfs . s2 [ "reads_passed_filter" ] = 0
# # parse new values from cutadapt results output
lines = res1 . strip ( ) . split ( "\n" )
for line in lines :
    if "Total reads processed:" in line :
        value = int ( line . split ( ) [ 3 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "reads_raw" ] = value
    if "Reads with adapters:" in line :
        value = int ( line . split ( ) [ 3 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "trim_adapter_bp_read1" ] = value
    if "Quality-trimmed" in line :
        value = int ( line . split ( ) [ 1 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "trim_quality_bp_read1" ] = value
    if "Reads that were too short" in line :
        value = int ( line . split ( ) [ 5 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "reads_filtered_by_minlen" ] = value
    if "Reads with too many N" in line :
        value = int ( line . split ( ) [ 5 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "reads_filtered_by_Ns" ] = value
    if "Reads written (passing filters):" in line :
        value = int ( line . split ( ) [ 4 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "reads_passed_filter" ] = value
# # save to stats summary
if sample . stats_dfs . s2 . reads_passed_filter :
    sample . stats . state = 2
    sample . stats . reads_passed_filter = sample . stats_dfs . s2 . reads_passed_filter
    sample . files . edits = [ ( OPJ ( data . dirs . edits , sample . name + ".trimmed_R1_.fastq.gz" ) , 0 ) ]
    # # write the long form output to the log file .
    LOGGER . info ( res1 )
else :
    print ( "{}No reads passed filtering in Sample: {}" . format ( data . _spacer , sample . name ) )


signature:
def parse_pair_results ( data , sample , res ) :
    """parse results from cutadapt for paired data"""

implementation:
LOGGER . info ( "in parse pair mod results\n%s" , res )
# # set default values
sample . stats_dfs . s2 [ "trim_adapter_bp_read1" ] = 0
sample . stats_dfs . s2 [ "trim_adapter_bp_read2" ] = 0
sample . stats_dfs . s2 [ "trim_quality_bp_read1" ] = 0
sample . stats_dfs . s2 [ "trim_quality_bp_read2" ] = 0
sample . stats_dfs . s2 [ "reads_filtered_by_Ns" ] = 0
sample . stats_dfs . s2 [ "reads_filtered_by_minlen" ] = 0
sample . stats_dfs . s2 [ "reads_passed_filter" ] = 0
lines = res . strip ( ) . split ( "\n" )
qprimed = 0
for line in lines : # # set primer to catch next line
    if "Quality-trimmed" in line :
        qprimed = 1
    # # grab read1 and read2 lines when qprimed
    if "Read 1:" in line :
        if qprimed :
            value = int ( line . split ( ) [ 2 ] . replace ( "," , "" ) )
            sample . stats_dfs . s2 [ "trim_quality_bp_read1" ] = value
    if "Read 2:" in line :
        if qprimed :
            value = int ( line . split ( ) [ 2 ] . replace ( "," , "" ) )
            sample . stats_dfs . s2 [ "trim_quality_bp_read2" ] = value
            qprimed = 0
    if "Read 1 with adapter:" in line :
        value = int ( line . split ( ) [ 4 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "trim_adapter_bp_read1" ] = value
    if "Read 2 with adapter:" in line :
        value = int ( line . split ( ) [ 4 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "trim_adapter_bp_read2" ] = value
    if "Total read pairs processed:" in line :
        value = int ( line . split ( ) [ 4 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "reads_raw" ] = value
    if "Pairs that were too short" in line :
        value = int ( line . split ( ) [ 5 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "reads_filtered_by_minlen" ] = value
    if "Pairs with too many N" in line :
        value = int ( line . split ( ) [ 5 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "reads_filtered_by_Ns" ] = value
    if "Pairs written (passing filters):" in line :
        value = int ( line . split ( ) [ 4 ] . replace ( "," , "" ) )
        sample . stats_dfs . s2 [ "reads_passed_filter" ] = value
# # save to stats summary
if sample . stats_dfs . s2 . reads_passed_filter :
    sample . stats . state = 2
    sample . stats . reads_passed_filter = sample . stats_dfs . s2 . reads_passed_filter
    sample . files . edits = [ ( OPJ ( data . dirs . edits , sample . name + ".trimmed_R1_.fastq.gz" ) , OPJ ( data . dirs . edits , sample . name + ".trimmed_R2_.fastq.gz" ) ) ]
else :
    print ( "No reads passed filtering in Sample: {}" . format ( sample . name ) )


signature:
def cutadaptit_single ( data , sample ) :
    """Applies quality and adapter filters to reads using cutadapt . If the ipyrad
 filter param is set to 0 then it only filters to hard trim edges and uses
 mintrimlen . If filter = 1 , we add quality filters . If filter = 2 we add
 adapter filters ."""

implementation:
sname = sample . name
# # if ( GBS , ddRAD ) we look for the second cut site + adapter . For single - end
# # data we don ' t bother trying to remove the second barcode since it ' s not
# # as critical as with PE data .
if data . paramsdict [ "datatype" ] == "rad" :
    adapter = data . _hackersonly [ "p3_adapter" ]
else : # # if GBS then the barcode can also be on the other side .
    if data . paramsdict [ "datatype" ] == "gbs" : # # make full adapter ( - revcompcut - revcompbarcode - adapter )
    # # and add adapter without revcompbarcode
        if data . barcodes :
            adapter = fullcomp ( data . paramsdict [ "restriction_overhang" ] [ 1 ] ) [ : : - 1 ] + fullcomp ( data . barcodes [ sample . name ] ) [ : : - 1 ] + data . _hackersonly [ "p3_adapter" ]
            # # add incomplete adapter to extras ( - recompcut - adapter )
            data . _hackersonly [ "p3_adapters_extra" ] . append ( fullcomp ( data . paramsdict [ "restriction_overhang" ] [ 1 ] ) [ : : - 1 ] + data . _hackersonly [ "p3_adapter" ] )
        else :
            LOGGER . warning ( "No barcode information present, and is therefore not " + "being used for adapter trimming of SE gbs data." )
            # # else no search for barcodes on 3'
            adapter = fullcomp ( data . paramsdict [ "restriction_overhang" ] [ 1 ] ) [ : : - 1 ] + data . _hackersonly [ "p3_adapter" ]
    else :
        adapter = fullcomp ( data . paramsdict [ "restriction_overhang" ] [ 1 ] ) [ : : - 1 ] + data . _hackersonly [ "p3_adapter" ]
# # get length trim parameter from new or older version of ipyrad params
trim5r1 = trim3r1 = [ ]
if data . paramsdict . get ( "trim_reads" ) :
    trimlen = data . paramsdict . get ( "trim_reads" )
    # # trim 5 ' end
    if trimlen [ 0 ] :
        trim5r1 = [ "-u" , str ( trimlen [ 0 ] ) ]
    if trimlen [ 1 ] < 0 :
        trim3r1 = [ "-u" , str ( trimlen [ 1 ] ) ]
    if trimlen [ 1 ] > 0 :
        trim3r1 = [ "--length" , str ( trimlen [ 1 ] ) ]
else :
    trimlen = data . paramsdict . get ( "edit_cutsites" )
    trim5r1 = [ "--cut" , str ( trimlen [ 0 ] ) ]
# # testing new ' trim _ reads ' setting
cmdf1 = [ "cutadapt" ]
if trim5r1 :
    cmdf1 += trim5r1
if trim3r1 :
    cmdf1 += trim3r1
cmdf1 += [ "--minimum-length" , str ( data . paramsdict [ "filter_min_trim_len" ] ) , "--max-n" , str ( data . paramsdict [ "max_low_qual_bases" ] ) , "--trim-n" , "--output" , OPJ ( data . dirs . edits , sname + ".trimmed_R1_.fastq.gz" ) , sample . files . concat [ 0 ] [ 0 ] ]
if int ( data . paramsdict [ "filter_adapters" ] ) : # # NEW : only quality trim the 3 ' end for SE data .
    cmdf1 . insert ( 1 , "20" )
    cmdf1 . insert ( 1 , "-q" )
    cmdf1 . insert ( 1 , str ( data . paramsdict [ "phred_Qscore_offset" ] ) )
    cmdf1 . insert ( 1 , "--quality-base" )
# # if filter _ adapters = = 3 then p3 _ adapters _ extra will already have extra
# # poly adapters added to its list .
if int ( data . paramsdict [ "filter_adapters" ] ) > 1 : # # first enter extra cuts ( order of input is reversed )
    for extracut in list ( set ( data . _hackersonly [ "p3_adapters_extra" ] ) ) [ : : - 1 ] :
        cmdf1 . insert ( 1 , extracut )
        cmdf1 . insert ( 1 , "-a" )
    # # then put the main cut so it appears first in command
    cmdf1 . insert ( 1 , adapter )
    cmdf1 . insert ( 1 , "-a" )
# # do modifications to read1 and write to tmp file
LOGGER . info ( cmdf1 )
proc1 = sps . Popen ( cmdf1 , stderr = sps . STDOUT , stdout = sps . PIPE , close_fds = True )
try :
    res1 = proc1 . communicate ( ) [ 0 ]
except KeyboardInterrupt :
    proc1 . kill ( )
    raise KeyboardInterrupt
# # raise errors if found
if proc1 . returncode :
    raise IPyradWarningExit ( " error in {}\n {}" . format ( " " . join ( cmdf1 ) , res1 ) )
# # return result string to be parsed outside of engine
return res1


signature:
def cutadaptit_pairs ( data , sample ) :
    """Applies trim & filters to pairs , including adapter detection . If we have
 barcode information then we use it to trim reversecut + bcode + adapter from
 reverse read , if not then we have to apply a more general cut to make sure
 we remove the barcode , this uses wildcards and so will have more false
 positives that trim a little extra from the ends of reads . Should we add
 a warning about this when filter _ adapters = 2 and no barcodes ?"""

implementation:
LOGGER . debug ( "Entering cutadaptit_pairs - {}" . format ( sample . name ) )
sname = sample . name
# # applied to read pairs
# trim _ r1 = str ( data . paramsdict [ " edit _ cutsites " ] [ 0 ] )
# trim _ r2 = str ( data . paramsdict [ " edit _ cutsites " ] [ 1 ] )
finput_r1 = sample . files . concat [ 0 ] [ 0 ]
finput_r2 = sample . files . concat [ 0 ] [ 1 ]
# # Get adapter sequences . This is very important . For the forward adapter
# # we don ' t care all that much about getting the sequence just before the
# # Illumina adapter , b / c it will either be random ( in RAD ) , or the reverse
# # cut site of cut1 or cut2 ( gbs or ddrad ) . Either way , we can still trim it
# # off later in step7 with trim overhang if we want . And it should be invar -
# # iable unless the cut site has an ambiguous char . The reverse adapter is
# # super important , however b / c it can contain the inline barcode and
# # revcomp cut site . We def want to trim out the barcode , and ideally the
# # cut site too to be safe . Problem is we don ' t always know the barcode if
# # users demultiplexed their data elsewhere . So , if barcode is missing we
# # do a very fuzzy match before the adapter and trim it out .
# # this just got more complicated now that we allow merging technical
# # replicates in step 1 since a single sample might have multiple barcodes
# # associated with it and so we need to search for multiple adapter + barcode
# # combinations .
# # We will assume that if they are ' linking _ barcodes ( ) ' here then there are
# # no technical replicates in the barcodes file . If there ARE technical
# # replicates , then they should run step1 so they are merged , in which case
# # the sample specific barcodes will be saved to each Sample under its
# # . barcode attribute as a list .
if not data . barcodes : # # try linking barcodes again in case user just added a barcodes path
# # after receiving the warning . We assume no technical replicates here .
    try :
        data . _link_barcodes ( )
    except Exception as inst :
        LOGGER . warning ( "  error adding barcodes info: %s" , inst )
# # barcodes are present meaning they were parsed to the samples in step 1.
if data . barcodes :
    try :
        adapter1 = fullcomp ( data . paramsdict [ "restriction_overhang" ] [ 1 ] ) [ : : - 1 ] + data . _hackersonly [ "p3_adapter" ]
        if isinstance ( sample . barcode , list ) :
            bcode = fullcomp ( sample . barcode [ 0 ] ) [ : : - 1 ]
        elif isinstance ( data . barcodes [ sample . name ] , list ) :
            bcode = fullcomp ( data . barcodes [ sample . name ] [ 0 ] [ : : - 1 ] )
        else :
            bcode = fullcomp ( data . barcodes [ sample . name ] ) [ : : - 1 ]
        # # add full adapter ( - revcompcut - revcompbcode - adapter )
        adapter2 = fullcomp ( data . paramsdict [ "restriction_overhang" ] [ 0 ] ) [ : : - 1 ] + bcode + data . _hackersonly [ "p5_adapter" ]
    except KeyError as inst :
        msg = """
    Sample name does not exist in the barcode file. The name in the barcode file
    for each sample must exactly equal the raw file name for the sample minus
    `_R1`. So for example a sample called WatDo_PipPrep_R1_100.fq.gz must
    be referenced in the barcode file as WatDo_PipPrep_100. The name in your
    barcode file for this sample must match: {}
    """ . format ( sample . name )
        LOGGER . error ( msg )
        raise IPyradWarningExit ( msg )
else :
    print ( NO_BARS_GBS_WARNING )
    # adapter1 = fullcomp ( data . paramsdict [ " restriction _ overhang " ] [ 1 ] ) [ : : - 1 ] + \
    # data . _ hackersonly [ " p3 _ adapter " ]
    # adapter2 = " XXX "
    adapter1 = data . _hackersonly [ "p3_adapter" ]
    adapter2 = fullcomp ( data . _hackersonly [ "p5_adapter" ] )
# # parse trim _ reads
trim5r1 = trim5r2 = trim3r1 = trim3r2 = [ ]
if data . paramsdict . get ( "trim_reads" ) :
    trimlen = data . paramsdict . get ( "trim_reads" )
    # # trim 5 ' end
    if trimlen [ 0 ] :
        trim5r1 = [ "-u" , str ( trimlen [ 0 ] ) ]
    if trimlen [ 1 ] < 0 :
        trim3r1 = [ "-u" , str ( trimlen [ 1 ] ) ]
    if trimlen [ 1 ] > 0 :
        trim3r1 = [ "--length" , str ( trimlen [ 1 ] ) ]
    # # legacy support for trimlen = 0,0 default
    if len ( trimlen ) > 2 :
        if trimlen [ 2 ] :
            trim5r2 = [ "-U" , str ( trimlen [ 2 ] ) ]
    if len ( trimlen ) > 3 :
        if trimlen [ 3 ] :
            if trimlen [ 3 ] < 0 :
                trim3r2 = [ "-U" , str ( trimlen [ 3 ] ) ]
            if trimlen [ 3 ] > 0 :
                trim3r2 = [ "--length" , str ( trimlen [ 3 ] ) ]
else : # # legacy support
    trimlen = data . paramsdict . get ( "edit_cutsites" )
    trim5r1 = [ "-u" , str ( trimlen [ 0 ] ) ]
    trim5r2 = [ "-U" , str ( trimlen [ 1 ] ) ]
# # testing new ' trim _ reads ' setting
cmdf1 = [ "cutadapt" ]
if trim5r1 :
    cmdf1 += trim5r1
if trim3r1 :
    cmdf1 += trim3r1
if trim5r2 :
    cmdf1 += trim5r2
if trim3r2 :
    cmdf1 += trim3r2
cmdf1 += [ "--trim-n" , "--max-n" , str ( data . paramsdict [ "max_low_qual_bases" ] ) , "--minimum-length" , str ( data . paramsdict [ "filter_min_trim_len" ] ) , "-o" , OPJ ( data . dirs . edits , sname + ".trimmed_R1_.fastq.gz" ) , "-p" , OPJ ( data . dirs . edits , sname + ".trimmed_R2_.fastq.gz" ) , finput_r1 , finput_r2 ]
# # additional args
if int ( data . paramsdict [ "filter_adapters" ] ) < 2 : # # add a dummy adapter to let cutadapt know whe are not using legacy - mode
    cmdf1 . insert ( 1 , "XXX" )
    cmdf1 . insert ( 1 , "-A" )
if int ( data . paramsdict [ "filter_adapters" ] ) :
    cmdf1 . insert ( 1 , "20,20" )
    cmdf1 . insert ( 1 , "-q" )
    cmdf1 . insert ( 1 , str ( data . paramsdict [ "phred_Qscore_offset" ] ) )
    cmdf1 . insert ( 1 , "--quality-base" )
if int ( data . paramsdict [ "filter_adapters" ] ) > 1 : # # if technical replicates then add other copies
    if isinstance ( sample . barcode , list ) :
        for extrabar in sample . barcode [ 1 : ] :
            data . _hackersonly [ "p5_adapters_extra" ] += fullcomp ( data . paramsdict [ "restriction_overhang" ] [ 0 ] ) [ : : - 1 ] + fullcomp ( extrabar ) [ : : - 1 ] + data . _hackersonly [ "p5_adapter" ]
            data . _hackersonly [ "p5_adapters_extra" ] += fullcomp ( data . paramsdict [ "restriction_overhang" ] [ 1 ] ) [ : : - 1 ] + data . _hackersonly [ "p3_adapter" ]
    # # first enter extra cuts
    zcut1 = list ( set ( data . _hackersonly [ "p3_adapters_extra" ] ) ) [ : : - 1 ]
    zcut2 = list ( set ( data . _hackersonly [ "p5_adapters_extra" ] ) ) [ : : - 1 ]
    for ecut1 , ecut2 in zip ( zcut1 , zcut2 ) :
        cmdf1 . insert ( 1 , ecut1 )
        cmdf1 . insert ( 1 , "-a" )
        cmdf1 . insert ( 1 , ecut2 )
        cmdf1 . insert ( 1 , "-A" )
    # # then put the main cut first
    cmdf1 . insert ( 1 , adapter1 )
    cmdf1 . insert ( 1 , '-a' )
    cmdf1 . insert ( 1 , adapter2 )
    cmdf1 . insert ( 1 , '-A' )
# # do modifications to read1 and write to tmp file
LOGGER . debug ( " " . join ( cmdf1 ) )
# sys . exit ( )
try :
    proc1 = sps . Popen ( cmdf1 , stderr = sps . STDOUT , stdout = sps . PIPE , close_fds = True )
    res1 = proc1 . communicate ( ) [ 0 ]
except KeyboardInterrupt :
    proc1 . kill ( )
    LOGGER . info ( "this is where I want it to interrupt" )
    raise KeyboardInterrupt ( )
# # raise errors if found
if proc1 . returncode :
    raise IPyradWarningExit ( " error [returncode={}]: {}\n{}" . format ( proc1 . returncode , " " . join ( cmdf1 ) , res1 ) )
LOGGER . debug ( "Exiting cutadaptit_pairs - {}" . format ( sname ) )
# # return results string to be parsed outside of engine
return res1


signature:
def run2 ( data , samples , force , ipyclient ) :
    """Filter for samples that are already finished with this step , allow others
 to run , pass them to parallel client function to filter with cutadapt ."""

implementation:
# # create output directories
data . dirs . edits = os . path . join ( os . path . realpath ( data . paramsdict [ "project_dir" ] ) , data . name + "_edits" )
if not os . path . exists ( data . dirs . edits ) :
    os . makedirs ( data . dirs . edits )
# # get samples
subsamples = choose_samples ( samples , force )
# # only allow extra adapters in filters = = 3,
# # and add poly repeats if not in list of adapters
if int ( data . paramsdict [ "filter_adapters" ] ) == 3 :
    if not data . _hackersonly [ "p3_adapters_extra" ] :
        for poly in [ "A" * 8 , "T" * 8 , "C" * 8 , "G" * 8 ] :
            data . _hackersonly [ "p3_adapters_extra" ] . append ( poly )
    if not data . _hackersonly [ "p5_adapters_extra" ] :
        for poly in [ "A" * 8 , "T" * 8 , "C" * 8 , "G" * 8 ] :
            data . _hackersonly [ "p5_adapters_extra" ] . append ( poly )
else :
    data . _hackersonly [ "p5_adapters_extra" ] = [ ]
    data . _hackersonly [ "p3_adapters_extra" ] = [ ]
# # concat is not parallelized ( since it ' s disk limited , generally )
subsamples = concat_reads ( data , subsamples , ipyclient )
# # cutadapt is parallelized by ncores / 2 because cutadapt spawns threads
lbview = ipyclient . load_balanced_view ( targets = ipyclient . ids [ : : 2 ] )
run_cutadapt ( data , subsamples , lbview )
# # cleanup is . . .
assembly_cleanup ( data )


signature:
def concat_reads ( data , subsamples , ipyclient ) :
    """concatenate if multiple input files for a single samples"""

implementation:
# # concatenate reads if they come from merged assemblies .
if any ( [ len ( i . files . fastqs ) > 1 for i in subsamples ] ) : # # run on single engine for now
    start = time . time ( )
    printstr = " concatenating inputs  | {} | s2 |"
    finished = 0
    catjobs = { }
    for sample in subsamples :
        if len ( sample . files . fastqs ) > 1 :
            catjobs [ sample . name ] = ipyclient [ 0 ] . apply ( concat_multiple_inputs , * ( data , sample ) )
        else :
            sample . files . concat = sample . files . fastqs
    # # wait for all to finish
    while 1 :
        finished = sum ( [ i . ready ( ) for i in catjobs . values ( ) ] )
        elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
        progressbar ( len ( catjobs ) , finished , printstr . format ( elapsed ) , spacer = data . _spacer )
        time . sleep ( 0.1 )
        if finished == len ( catjobs ) :
            print ( "" )
            break
    # # collect results , which are concat file handles .
    for async in catjobs :
        if catjobs [ async ] . successful ( ) :
            data . samples [ async ] . files . concat = catjobs [ async ] . result ( )
        else :
            error = catjobs [ async ] . result ( )
            # exception ( )
            LOGGER . error ( "error in step2 concat %s" , error )
            raise IPyradWarningExit ( "error in step2 concat: {}" . format ( error ) )
else :
    for sample in subsamples : # # just copy fastqs handles to concat attribute
        sample . files . concat = sample . files . fastqs
return subsamples


signature:
def run_cutadapt ( data , subsamples , lbview ) :
    """sends fastq files to cutadapt"""

implementation:
# # choose cutadapt function based on datatype
start = time . time ( )
printstr = " processing reads      | {} | s2 |"
finished = 0
rawedits = { }
# # sort subsamples so that the biggest files get submitted first
subsamples . sort ( key = lambda x : x . stats . reads_raw , reverse = True )
LOGGER . info ( [ i . stats . reads_raw for i in subsamples ] )
# # send samples to cutadapt filtering
if "pair" in data . paramsdict [ "datatype" ] :
    for sample in subsamples :
        rawedits [ sample . name ] = lbview . apply ( cutadaptit_pairs , * ( data , sample ) )
else :
    for sample in subsamples :
        rawedits [ sample . name ] = lbview . apply ( cutadaptit_single , * ( data , sample ) )
# # wait for all to finish
while 1 :
    finished = sum ( [ i . ready ( ) for i in rawedits . values ( ) ] )
    elapsed = datetime . timedelta ( seconds = int ( time . time ( ) - start ) )
    progressbar ( len ( rawedits ) , finished , printstr . format ( elapsed ) , spacer = data . _spacer )
    time . sleep ( 0.1 )
    if finished == len ( rawedits ) :
        print ( "" )
        break
# # collect results , report failures , and store stats . async = sample . name
for async in rawedits :
    if rawedits [ async ] . successful ( ) :
        res = rawedits [ async ] . result ( )
        # # if single cleanup is easy
        if "pair" not in data . paramsdict [ "datatype" ] :
            parse_single_results ( data , data . samples [ async ] , res )
        else :
            parse_pair_results ( data , data . samples [ async ] , res )
    else :
        print ( "  found an error in step2; see ipyrad_log.txt" )
        LOGGER . error ( "error in run_cutadapt(): %s" , rawedits [ async ] . exception ( ) )


signature:
def choose_samples ( samples , force ) :
    """filter out samples that are already done with this step , unless force"""

implementation:
# # hold samples that pass
subsamples = [ ]
# # filter the samples again
if not force :
    for sample in samples :
        if sample . stats . state >= 2 :
            print ( """\
    Skipping Sample {}; Already filtered. Use force argument to overwrite.\
    """ . format ( sample . name ) )
        elif not sample . stats . reads_raw :
            print ( """\
    Skipping Sample {}; No reads found in file {}\
    """ . format ( sample . name , sample . files . fastqs ) )
        else :
            subsamples . append ( sample )
else :
    for sample in samples :
        if not sample . stats . reads_raw :
            print ( """\
    Skipping Sample {}; No reads found in file {}\
    """ . format ( sample . name , sample . files . fastqs ) )
        else :
            subsamples . append ( sample )
return subsamples


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

implementation:
# # if more than one tuple in fastq list
if len ( sample . files . fastqs ) > 1 : # # create a cat command to append them all ( doesn ' t matter if they
# # are gzipped , cat still works ) . Grab index 0 of tuples for R1s .
    cmd1 = [ "cat" ] + [ i [ 0 ] for i in sample . files . fastqs ]
    isgzip = ".gz"
    if not sample . files . fastqs [ 0 ] [ 0 ] . endswith ( ".gz" ) :
        isgzip = ""
    # # write to new concat handle
    conc1 = os . path . join ( data . dirs . edits , sample . name + "_R1_concat.fq{}" . format ( isgzip ) )
    with open ( conc1 , 'w' ) as cout1 :
        proc1 = sps . Popen ( cmd1 , stderr = sps . STDOUT , stdout = cout1 , close_fds = True )
        res1 = proc1 . communicate ( ) [ 0 ]
    if proc1 . returncode :
        raise IPyradWarningExit ( "error in: {}, {}" . format ( cmd1 , res1 ) )
    # # Only set conc2 if R2 actually exists
    conc2 = 0
    if "pair" in data . paramsdict [ "datatype" ] :
        cmd2 = [ "cat" ] + [ i [ 1 ] for i in sample . files . fastqs ]
        conc2 = os . path . join ( data . dirs . edits , sample . name + "_R2_concat.fq{}" . format ( isgzip ) )
        with open ( conc2 , 'w' ) as cout2 :
            proc2 = sps . Popen ( cmd2 , stderr = sps . STDOUT , stdout = cout2 , close_fds = True )
            res2 = proc2 . communicate ( ) [ 0 ]
        if proc2 . returncode :
            raise IPyradWarningExit ( "Error concatenating fastq files. Make sure all " + "these files exist: {}\nError message: {}" . format ( cmd2 , proc2 . returncode ) )
    # # store new file handles
    sample . files . concat = [ ( conc1 , conc2 ) ]
return sample . files . concat


signature:
def make ( data , samples ) :
    """Convert vcf from step6 to . loci format to facilitate downstream format conversion"""

implementation:
invcffile = os . path . join ( data . dirs . consens , data . name + ".vcf" )
outlocifile = os . path . join ( data . dirs . outfiles , data . name + ".loci" )
importvcf ( invcffile , outlocifile )


signature:
def importvcf ( vcffile , locifile ) :
    """Function for importing a vcf file into loci format . Arguments
 are the input vcffile and the loci file to write out ."""

implementation:
try : # # Get names of all individuals in the vcf
    with open ( invcffile , 'r' ) as invcf :
        for line in invcf :
            if line . split ( ) [ 0 ] == "#CHROM" : # # This is maybe a little clever . The names in the vcf are everything after
            # # the " FORMAT " column , so find that index , then slice everything after it .
                names_col = line . split ( ) . index ( "FORMAT" ) + 1
                names = line . split ( ) [ names_col : ]
                LOGGER . debug ( "Got names - %s" , names )
                break
        print ( "wat" )
    # # Get the column to start reading at
except Exception :
    print ( "wat" )


signature:
def get_targets ( ipyclient ) :
    """A function to find 2 engines per hostname on the ipyclient .
 We ' ll assume that the CPUs are hyperthreaded , which is why
 we grab two . If they are not then no foul . Two multi - threaded
 jobs will be run on each of the 2 engines per host ."""

implementation:
# # fill hosts with async [ gethostname ]
hosts = [ ]
for eid in ipyclient . ids :
    engine = ipyclient [ eid ]
    if not engine . outstanding :
        hosts . append ( engine . apply ( socket . gethostname ) )
# # capture results of asyncs
hosts = [ i . get ( ) for i in hosts ]
hostset = set ( hosts )
hostzip = zip ( hosts , ipyclient . ids )
hostdict = { host : [ i [ 1 ] for i in hostzip if i [ 0 ] == host ] for host in hostset }
targets = list ( itertools . chain ( * [ hostdict [ i ] [ : 2 ] for i in hostdict ] ) )
# # return first two engines from each host
return targets


signature:
def compute_tree_stats ( self , ipyclient ) :
    """compute stats for stats file and NHX tree features"""

implementation:
# # get name indices
names = self . samples
# # get majority rule consensus tree of weighted Q bootstrap trees
if self . params . nboots : # # Tree object
    fulltre = ete3 . Tree ( self . trees . tree , format = 0 )
    fulltre . unroot ( )
    # # only grab as many boots as the last option said was max
    with open ( self . trees . boots , 'r' ) as inboots :
        bb = [ ete3 . Tree ( i . strip ( ) , format = 0 ) for i in inboots . readlines ( ) ]
        wboots = [ fulltre ] + bb [ - self . params . nboots : ]
    # # infer consensus tree and write to file
    wctre , wcounts = consensus_tree ( wboots , names = names )
    self . trees . cons = os . path . join ( self . dirs , self . name + ".cons" )
    with open ( self . trees . cons , 'w' ) as ocons :
        ocons . write ( wctre . write ( format = 0 ) )
else :
    wctre = ete3 . Tree ( self . trees . tree , format = 0 )
    wctre . unroot ( )
# # build stats file and write trees
self . trees . nhx = os . path . join ( self . dirs , self . name + ".nhx" )
with open ( self . files . stats , 'w' ) as ostats : # # print Tetrad info
# ostats . write ( STATS _ STRING . format ( * * self . stats ) )
# # print bootstrap splits
    if self . params . nboots :
        ostats . write ( "## splits observed in {} trees\n" . format ( len ( wboots ) ) )
        for i , j in enumerate ( self . samples ) :
            ostats . write ( "{:<3} {}\n" . format ( i , j ) )
        ostats . write ( "\n" )
        for split , freq in wcounts :
            if split . count ( '1' ) > 1 :
                ostats . write ( "{}   {:.2f}\n" . format ( split , round ( freq , 2 ) ) )
        ostats . write ( "\n" )
# # parallelized this function because it can be slogging
lbview = ipyclient . load_balanced_view ( )
# # store results in dicts
qtots = { }
qsamp = { }
tots = sum ( 1 for i in wctre . iter_leaves ( ) )
totn = set ( wctre . get_leaf_names ( ) )
# # iterate over node traversal .
for node in wctre . traverse ( ) : # # this is slow , needs to look at every sampled quartet
# # so we send it be processed on an engine
    qtots [ node ] = lbview . apply ( _get_total , * ( tots , node ) )
    qsamp [ node ] = lbview . apply ( _get_sampled , * ( self , totn , node ) )
# # wait for jobs to finish
ipyclient . wait ( )
# # put results into tree
for node in wctre . traverse ( ) : # # this is fast , just calcs n _ choose _ k
    total = qtots [ node ] . result ( )
    sampled = qsamp [ node ] . result ( )
    # # store the results to the tree
    node . add_feature ( "quartets_total" , total )
    node . add_feature ( "quartets_sampled" , sampled )
features = [ "quartets_total" , "quartets_sampled" ]
# # return as NHX format with extra info
with open ( self . trees . nhx , 'w' ) as outtre :
    outtre . write ( wctre . write ( format = 0 , features = features ) )


signature:
def random_combination ( iterable , nquartets ) :
    """Random selection from itertools . combinations ( iterable , r ) .
 Use this if not sampling all possible quartets ."""

implementation:
pool = tuple ( iterable )
size = len ( pool )
indices = random . sample ( xrange ( size ) , nquartets )
return tuple ( pool [ i ] for i in indices )


signature:
def random_product ( iter1 , iter2 ) :
    """random sampler for equal _ splits func"""

implementation:
pool1 = tuple ( iter1 )
pool2 = tuple ( iter2 )
ind1 = random . sample ( pool1 , 2 )
ind2 = random . sample ( pool2 , 2 )
return tuple ( ind1 + ind2 )


signature:
def n_choose_k ( n , k ) :
    """get the number of quartets as n - choose - k . This is used
 in equal splits to decide whether a split should be exhaustively sampled
 or randomly sampled . Edges near tips can be exhaustive while highly nested
 edges probably have too many quartets"""

implementation: return int ( reduce ( MUL , ( Fraction ( n - i , i + 1 ) for i in range ( k ) ) , 1 ) )


signature:
def count_snps ( mat ) :
    """get dstats from the count array and return as a float tuple"""

implementation:
# # get [ aabb , baba , abba , aaab ]
snps = np . zeros ( 4 , dtype = np . uint32 )
# # get concordant ( aabb ) pis sites
snps [ 0 ] = np . uint32 ( mat [ 0 , 5 ] + mat [ 0 , 10 ] + mat [ 0 , 15 ] + mat [ 5 , 0 ] + mat [ 5 , 10 ] + mat [ 5 , 15 ] + mat [ 10 , 0 ] + mat [ 10 , 5 ] + mat [ 10 , 15 ] + mat [ 15 , 0 ] + mat [ 15 , 5 ] + mat [ 15 , 10 ] )
# # get discordant ( baba ) sites
for i in range ( 16 ) :
    if i % 5 :
        snps [ 1 ] += mat [ i , i ]
# # get discordant ( abba ) sites
snps [ 2 ] = mat [ 1 , 4 ] + mat [ 2 , 8 ] + mat [ 3 , 12 ] + mat [ 4 , 1 ] + mat [ 6 , 9 ] + mat [ 7 , 13 ] + mat [ 8 , 2 ] + mat [ 9 , 6 ] + mat [ 11 , 14 ] + mat [ 12 , 3 ] + mat [ 13 , 7 ] + mat [ 14 , 11 ]
# # get autapomorphy sites
snps [ 3 ] = ( mat . sum ( ) - np . diag ( mat ) . sum ( ) ) - snps [ 2 ]
return snps


signature:
def subsample_snps_map ( seqchunk , nmask , maparr ) :
    """removes ncolumns from snparray prior to matrix calculation , and
 subsamples ' linked ' snps ( those from the same RAD locus ) such that
 for these four samples only 1 SNP per locus is kept . This information
 comes from the ' map ' array ( map file ) ."""

implementation:
# # mask columns that contain Ns
rmask = np . zeros ( seqchunk . shape [ 1 ] , dtype = np . bool_ )
# # apply mask to the mapfile
last_loc = - 1
for idx in xrange ( maparr . shape [ 0 ] ) :
    if maparr [ idx ] != last_loc :
        if not nmask [ idx ] :
            rmask [ idx ] = True
        last_loc = maparr [ idx ]
# # apply mask
# newarr = seqchunk [ : , rmask ]
# # return smaller Nmasked array
return rmask


signature:
def chunk_to_matrices ( narr , mapcol , nmask ) :
    """numba compiled code to get matrix fast .
 arr is a 4 x N seq matrix converted to np . int8
 I convert the numbers for ATGC into their respective index for the MAT
 matrix , and leave all others as high numbers , i . e . , - = = 45 , N = = 78."""

implementation:
# # get seq alignment and create an empty array for filling
mats = np . zeros ( ( 3 , 16 , 16 ) , dtype = np . uint32 )
# # replace ints with small ints that index their place in the
# # 16x16 . This no longer checks for big ints to exclude , so resolve = True
# # is now the default , TODO .
last_loc = - 1
for idx in xrange ( mapcol . shape [ 0 ] ) :
    if not nmask [ idx ] :
        if not mapcol [ idx ] == last_loc :
            i = narr [ : , idx ]
            mats [ 0 , ( 4 * i [ 0 ] ) + i [ 1 ] , ( 4 * i [ 2 ] ) + i [ 3 ] ] += 1
            last_loc = mapcol [ idx ]
# # fill the alternates
x = np . uint8 ( 0 )
for y in np . array ( [ 0 , 4 , 8 , 12 ] , dtype = np . uint8 ) :
    for z in np . array ( [ 0 , 4 , 8 , 12 ] , dtype = np . uint8 ) :
        mats [ 1 , y : y + np . uint8 ( 4 ) , z : z + np . uint8 ( 4 ) ] = mats [ 0 , x ] . reshape ( 4 , 4 )
        mats [ 2 , y : y + np . uint8 ( 4 ) , z : z + np . uint8 ( 4 ) ] = mats [ 0 , x ] . reshape ( 4 , 4 ) . T
        x += np . uint8 ( 1 )
return mats


signature:
def calculate ( seqnon , mapcol , nmask , tests ) :
    """groups together several numba compiled funcs"""

implementation:
# # create empty matrices
# LOGGER . info ( " tests [ 0 ] % s " , tests [ 0 ] )
# LOGGER . info ( ' seqnon [ [ tests [ 0 ] ] ] % s ' , seqnon [ [ tests [ 0 ] ] ] )
mats = chunk_to_matrices ( seqnon , mapcol , nmask )
# # empty svdscores for each arrangement of seqchunk
svds = np . zeros ( ( 3 , 16 ) , dtype = np . float64 )
qscores = np . zeros ( 3 , dtype = np . float64 )
ranks = np . zeros ( 3 , dtype = np . float64 )
for test in range ( 3 ) : # # get svd scores
    svds [ test ] = np . linalg . svd ( mats [ test ] . astype ( np . float64 ) ) [ 1 ]
    ranks [ test ] = np . linalg . matrix_rank ( mats [ test ] . astype ( np . float64 ) )
# # get minrank , or 11
minrank = int ( min ( 11 , ranks . min ( ) ) )
for test in range ( 3 ) :
    qscores [ test ] = np . sqrt ( np . sum ( svds [ test , minrank : ] ** 2 ) )
# # sort to find the best qorder
best = np . where ( qscores == qscores . min ( ) ) [ 0 ]
# best = qscores [ qscores = = qscores . min ( ) ] [ 0]
bidx = tests [ best ] [ 0 ]
qsnps = count_snps ( mats [ best ] [ 0 ] )
return bidx , qsnps


signature:
def nworker ( data , smpchunk , tests ) :
    """The workhorse function . Not numba ."""

implementation:
# # tell engines to limit threads
# numba . config . NUMBA _ DEFAULT _ NUM _ THREADS = 1
# # open the seqarray view , the modified array is in bootsarr
with h5py . File ( data . database . input , 'r' ) as io5 :
    seqview = io5 [ "bootsarr" ] [ : ]
    maparr = io5 [ "bootsmap" ] [ : ]
# # create an N - mask array of all seq cols ( this isn ' t really too slow )
nall_mask = seqview [ : ] == 78
# # tried numba compiling everythign below here , but was not faster
# # than making nmask w / axis arg in numpy
# # get the input arrays ready
rquartets = np . zeros ( ( smpchunk . shape [ 0 ] , 4 ) , dtype = np . uint16 )
rweights = None
# rweights = np . ones ( smpchunk . shape [ 0 ] , dtype = np . float64)
rdstats = np . zeros ( ( smpchunk . shape [ 0 ] , 4 ) , dtype = np . uint32 )
# times = [ ]
# # fill arrays with results using numba funcs
for idx in xrange ( smpchunk . shape [ 0 ] ) : # # get seqchunk for 4 samples ( 4 , ncols )
    sidx = smpchunk [ idx ]
    seqchunk = seqview [ sidx ]
    # # get N - containing columns in 4 - array , and invariant sites .
    nmask = np . any ( nall_mask [ sidx ] , axis = 0 )
    nmask += np . all ( seqchunk == seqchunk [ 0 ] , axis = 0 )
    # # < - do we need this ?
    # # get matrices if there are any shared SNPs
    # # returns best - tree index , qscores , and qstats
    # bidx , qscores , qstats = calculate ( seqchunk , maparr [ : , 0 ] , nmask , tests )
    bidx , qstats = calculate ( seqchunk , maparr [ : , 0 ] , nmask , tests )
    # # get weights from the three scores sorted .
    # # Only save to file if the quartet has information
    rdstats [ idx ] = qstats
    rquartets [ idx ] = smpchunk [ idx ] [ bidx ]
return rquartets , rweights , rdstats


signature:
def shuffle_cols ( seqarr , newarr , cols ) :
    """used in bootstrap resampling without a map file"""

implementation:
for idx in xrange ( cols . shape [ 0 ] ) :
    newarr [ : , idx ] = seqarr [ : , cols [ idx ] ]
return newarr


signature:
def resolve_ambigs ( tmpseq ) :
    """returns a seq array with ' RSKYWM ' randomly replaced with resolved bases"""

implementation:
# # iterate over the bases ' RSKWYM ' : [ 82 , 83 , 75 , 87 , 89 , 77]
for ambig in np . uint8 ( [ 82 , 83 , 75 , 87 , 89 , 77 ] ) : # # get all site in this ambig
    idx , idy = np . where ( tmpseq == ambig )
    # # get the two resolutions of the ambig
    res1 , res2 = AMBIGS [ ambig . view ( "S1" ) ]
    # # randomly sample half those sites
    halfmask = np . random . choice ( [ True , False ] , idx . shape [ 0 ] )
    # # replace ambig bases with their resolutions
    for i in xrange ( halfmask . shape [ 0 ] ) :
        if halfmask [ i ] :
            tmpseq [ idx [ i ] , idy [ i ] ] = np . array ( res1 ) . view ( np . uint8 )
        else :
            tmpseq [ idx [ i ] , idy [ i ] ] = np . array ( res2 ) . view ( np . uint8 )
return tmpseq


signature:
def get_spans ( maparr , spans ) :
    """get span distance for each locus in original seqarray"""

implementation:
# # start at 0 , finds change at 1 - index of map file
bidx = 1
spans = np . zeros ( ( maparr [ - 1 , 0 ] , 2 ) , np . uint64 )
# # read through marr and record when locus id changes
for idx in xrange ( 1 , maparr . shape [ 0 ] ) :
    cur = maparr [ idx , 0 ]
    if cur != bidx :
        idy = idx + 1
        spans [ cur - 2 , 1 ] = idx
        spans [ cur - 1 , 0 ] = idx
        bidx = cur
spans [ - 1 , 1 ] = maparr [ - 1 , - 1 ]
return spans


signature:
def get_shape ( spans , loci ) :
    """get shape of new bootstrap resampled locus array"""

implementation:
width = 0
for idx in xrange ( loci . shape [ 0 ] ) :
    width += spans [ loci [ idx ] , 1 ] - spans [ loci [ idx ] , 0 ]
return width


signature:
def fill_boot ( seqarr , newboot , newmap , spans , loci ) :
    """fills the new bootstrap resampled array"""

implementation:
# # column index
cidx = 0
# # resample each locus
for i in xrange ( loci . shape [ 0 ] ) : # # grab a random locus ' s columns
    x1 = spans [ loci [ i ] ] [ 0 ]
    x2 = spans [ loci [ i ] ] [ 1 ]
    cols = seqarr [ : , x1 : x2 ]
    # # randomize columns within colsq
    cord = np . random . choice ( cols . shape [ 1 ] , cols . shape [ 1 ] , replace = False )
    rcols = cols [ : , cord ]
    # # fill bootarr with n columns from seqarr
    # # the required length was already measured
    newboot [ : , cidx : cidx + cols . shape [ 1 ] ] = rcols
    # # fill bootmap with new map info
    newmap [ cidx : cidx + cols . shape [ 1 ] , 0 ] = i + 1
    # # advance column index
    cidx += cols . shape [ 1 ]
# # return the concatenated cols
return newboot , newmap


signature:
def _byteify ( data , ignore_dicts = False ) :
    """converts unicode to utf - 8 when reading in json files"""

implementation:
if isinstance ( data , unicode ) :
    return data . encode ( "utf-8" )
if isinstance ( data , list ) :
    return [ _byteify ( item , ignore_dicts = True ) for item in data ]
if isinstance ( data , dict ) and not ignore_dicts :
    return { _byteify ( key , ignore_dicts = True ) : _byteify ( value , ignore_dicts = True ) for key , value in data . iteritems ( ) }
return data


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 _filter_clades ( clade_counts , cutoff ) :
    """A subfunc of consensus _ tree ( ) . Removes clades that occur
 with freq < cutoff ."""

implementation:
# # store clades that pass filter
passed = [ ]
clades = np . array ( [ list ( i [ 0 ] ) for i in clade_counts ] , dtype = np . int8 )
counts = np . array ( [ i [ 1 ] for i in clade_counts ] , dtype = np . float64 )
for idx in xrange ( clades . shape [ 0 ] ) :
    conflict = False
    if counts [ idx ] < cutoff :
        continue
    if np . sum ( clades [ idx ] ) > 1 : # check the current clade against all the accepted clades to see if
    # it conflicts . A conflict is defined as :
    # 1 . the clades are not disjoint
    # 2 . neither clade is a subset of the other
    # OR :
    # 1 . it is inverse of clade ( affects only < fake > root state )
    # because at root node it mirror images { 0011 : 95 } , { 1100 : 5 } .
        for aidx in passed : # intersect = clade . intersection ( accepted _ clade )
            summed = clades [ idx ] + clades [ aidx ]
            intersect = np . max ( summed ) > 1
            subset_test0 = np . all ( clades [ idx ] - clades [ aidx ] >= 0 )
            subset_test1 = np . all ( clades [ aidx ] - clades [ idx ] >= 0 )
            invert_test = np . bool_ ( clades [ aidx ] ) != np . bool_ ( clades [ idx ] )
            if np . all ( invert_test ) :
                counts [ aidx ] += counts [ idx ]
                conflict = True
            if intersect :
                if ( not subset_test0 ) and ( not subset_test1 ) :
                    conflict = True
    if conflict == False :
        passed . append ( idx )
# # rebuild the dict
rclades = [ ]
# j for i , j in enumerate ( clade _ counts ) if i in passed ]
# # set the counts to include mirrors
for idx in passed :
    rclades . append ( ( clades [ idx ] , counts [ idx ] ) )
return rclades


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 , guidetreefile = self . files . guidetreefile , resolve = self . _resolve , 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 _parse_names ( self ) :
    """parse sample names from the sequence file"""

implementation:
self . samples = [ ]
with iter ( open ( self . files . data , 'r' ) ) as infile :
    infile . next ( ) . strip ( ) . split ( )
    while 1 :
        try :
            self . samples . append ( infile . next ( ) . split ( ) [ 0 ] )
        except StopIteration :
            break


signature:
def _init_seqarray ( self , quiet = False ) :
    """Fills the seqarr with the full data set , and creates a bootsarr copy
 with the following modifications :
 1 ) converts " - " into " N " s , since they are similarly treated as missing .
 2 ) randomly resolve ambiguities ( RSKWYM )
 3 ) convert to uint8 for smaller memory load and faster computation"""

implementation:
# # read in the data ( seqfile )
try :
    spath = open ( self . files . data , 'r' )
except IOError :
    raise IPyradWarningExit ( NO_SNP_FILE . format ( self . files . data ) )
line = spath . readline ( ) . strip ( ) . split ( )
ntax = int ( line [ 0 ] )
nbp = int ( line [ 1 ] )
# # make a tmp seq array
if not quiet :
    print ( "loading seq array [{} taxa x {} bp]" . format ( ntax , nbp ) )
tmpseq = np . zeros ( ( ntax , nbp ) , dtype = np . uint8 )
# # create array storage for real seq and the tmp bootstrap seqarray
with h5py . File ( self . database . input , 'w' ) as io5 :
    io5 . create_dataset ( "seqarr" , ( ntax , nbp ) , dtype = np . uint8 )
    io5 . create_dataset ( "bootsarr" , ( ntax , nbp ) , dtype = np . uint8 )
    io5 . create_dataset ( "bootsmap" , ( nbp , 2 ) , dtype = np . uint32 )
    # # if there is a map file , load it into the bootsmap
    if self . files . mapfile :
        with open ( self . files . mapfile , 'r' ) as inmap : # # parse the map file from txt and save as dataset
            maparr = np . genfromtxt ( inmap , dtype = np . uint64 )
            io5 [ "bootsmap" ] [ : ] = maparr [ : , [ 0 , 3 ] ]
            # # parse the span info from maparr and save to dataset
            spans = np . zeros ( ( maparr [ - 1 , 0 ] , 2 ) , np . uint64 )
            spans = get_spans ( maparr , spans )
            io5 . create_dataset ( "spans" , data = spans )
            if not quiet :
                print ( "max unlinked SNPs per quartet (nloci): {}" . format ( spans . shape [ 0 ] ) )
    else :
        io5 [ "bootsmap" ] [ : , 0 ] = np . arange ( io5 [ "bootsmap" ] . shape [ 0 ] )
    # # fill the tmp array from the input phy
    for line , seq in enumerate ( spath . readlines ( ) ) :
        tmpseq [ line ] = np . array ( list ( seq . split ( ) [ - 1 ] ) ) . view ( np . uint8 )
    # # convert ' - ' or ' _ ' into ' N '
    tmpseq [ tmpseq == 45 ] = 78
    tmpseq [ tmpseq == 95 ] = 78
    # # save array to disk so it can be easily accessed by slicing
    # # This unmodified array is used again later for sampling boots
    io5 [ "seqarr" ] [ : ] = tmpseq
    # # resolve ambiguous IUPAC codes
    if self . _resolve :
        tmpseq = resolve_ambigs ( tmpseq )
    # # convert CATG bases to matrix indices
    tmpseq [ tmpseq == 65 ] = 0
    tmpseq [ tmpseq == 67 ] = 1
    tmpseq [ tmpseq == 71 ] = 2
    tmpseq [ tmpseq == 84 ] = 3
    # # save modified array to disk
    io5 [ "bootsarr" ] [ : ] = tmpseq
    # # memory cleanup
    # del tmpseq
    # # get initial array
    LOGGER . info ( "original seqarr \n %s" , io5 [ "seqarr" ] [ : , : 20 ] )
    LOGGER . info ( "original bootsarr \n %s" , io5 [ "bootsarr" ] [ : , : 20 ] )
    LOGGER . info ( "original bootsmap \n %s" , io5 [ "bootsmap" ] [ : 20 , : ] )


signature:
def _store_N_samples ( self , ncpus ) :
    """Find all quartets of samples and store in a large array
 Create a chunk size for sampling from the array of quartets .
 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 = 4
if self . params . nquartets > 500000 :
    breaks = 8
# # chunk up the data
self . _chunksize = ( self . params . nquartets // ( breaks * ncpus ) + ( self . params . nquartets % ( breaks * ncpus ) ) )
LOGGER . info ( "nquarts = %s, chunk = %s" , self . params . nquartets , self . _chunksize )
# # ' samples ' stores the indices of the quartet .
# # ` quartets ` stores the correct quartet in the order ( 1,2 | 3,4)
# # ` weights ` stores the weight of the quartet in ' quartets '
# # we gzip this for now , but check later if this has a big speed cost
# # create h5 OUT empty arrays
with h5py . File ( self . database . output , 'w' ) as io5 :
    io5 . create_dataset ( "quartets" , ( self . params . nquartets , 4 ) , dtype = np . uint16 , chunks = ( self . _chunksize , 4 ) )
    io5 . create_dataset ( "qstats" , ( self . params . nquartets , 4 ) , dtype = np . uint32 , chunks = ( self . _chunksize , 4 ) )
    io5 . create_group ( "qboots" )
# # append to h5 IN array ( which also has seqarray ) and fill it
with h5py . File ( self . database . input , 'a' ) as io5 : # # create data sets
    io5 . create_dataset ( "samples" , ( self . params . nquartets , 4 ) , dtype = np . uint16 , chunks = ( self . _chunksize , 4 ) , compression = 'gzip' )
    # # populate array with all possible quartets . This allows us to
    # # sample from the total , and also to continue from a checkpoint
    qiter = itertools . combinations ( xrange ( len ( self . samples ) ) , 4 )
    i = 0
    # # fill chunksize at a time for efficiency
    while i < self . params . nquartets :
        if self . params . method != "all" : # # grab the next random 1000
            qiter = [ ]
            while len ( qiter ) < min ( self . _chunksize , io5 [ "samples" ] . shape [ 0 ] ) :
                qiter . append ( random_combination ( range ( len ( self . samples ) ) , 4 ) )
            dat = np . array ( qiter )
        else : # # grab the next ordered chunksize
            dat = np . array ( list ( itertools . islice ( qiter , self . _chunksize ) ) )
        # # store to h5
        io5 [ "samples" ] [ i : i + self . _chunksize ] = dat [ : io5 [ "samples" ] . shape [ 0 ] - i ]
        i += self . _chunksize


signature:
def _store_equal_samples ( self , ncpus ) :
    """sample quartets evenly across splits of the starting tree , and fills
 in remaining samples with random quartet samples . Uses a hash dict to
 not sample the same quartet twice , so for very large trees this can
 take a few minutes to find millions of possible quartet samples ."""

implementation:
# # choose chunker for h5 arr
breaks = 2
if self . params . nquartets < 5000 :
    breaks = 1
if self . params . nquartets > 100000 :
    breaks = 4
if self . params . nquartets > 500000 :
    breaks = 8
self . _chunksize = ( self . params . nquartets // ( breaks * ncpus ) + ( self . params . nquartets % ( breaks * ncpus ) ) )
LOGGER . info ( "nquarts = %s, chunk = %s" , self . params . nquartets , self . _chunksize )
# # create h5 OUT empty arrays
with h5py . File ( self . database . output , 'w' ) as io5 :
    io5 . create_dataset ( "quartets" , ( self . params . nquartets , 4 ) , dtype = np . uint16 , chunks = ( self . _chunksize , 4 ) )
    io5 . create_dataset ( "qstats" , ( self . params . nquartets , 4 ) , dtype = np . uint32 , chunks = ( self . _chunksize , 4 ) )
    io5 . create_group ( "qboots" )
# # get starting tree , unroot , randomly resolve , ladderize
tre = ete3 . Tree ( self . files . guidetreefile , format = 0 )
# tre = toytree . tree ( self . files . guidetreefile , format = 0)
tre . tree . unroot ( )
tre . tree . resolve_polytomy ( recursive = True )
tre . tree . ladderize ( )
# # randomly sample all splits of tree and convert tip names to indices
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 ( no single tips edges )
# # this seemed to cause problems with unsampled tips
splits = [ i for i in splits if all ( [ len ( j ) > 1 for j in i ] ) ]
# # turn each into an iterable split sampler
# # if the nquartets for that split is small , then sample all of them
# # if it is big , then make it a random sampler from that split
qiters = [ ]
# # how many min quartets are we gonna sample from each split ?
squarts = self . params . nquartets // len ( splits )
# # how many iterators can be sampled to saturation ?
nsaturation = 0
for split in 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 .
    if n_choose_k ( len ( split [ 0 ] ) , 2 ) * n_choose_k ( len ( split [ 1 ] ) , 2 ) < squarts * 2 :
        qiter = ( i + j for ( i , j ) in itertools . product ( itertools . combinations ( split [ 0 ] , 2 ) , itertools . combinations ( split [ 1 ] , 2 ) ) )
        nsaturation += 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 ) )
        nsaturation += 1
    # # store all iterators into a list
    qiters . append ( qiter )
# for split in splits :
# print ( split )
# # make qiters infinitely cycling
qiters = itertools . cycle ( qiters )
cycler = itertools . cycle ( range ( len ( splits ) ) )
# # store visiting quartets
sampled = set ( )
# # iterate over qiters sampling from each , if one runs out , keep
# # sampling from remaining qiters . Keep going until samples is filled
with h5py . File ( self . database . input , 'a' ) as io5 : # # create data sets
    io5 . create_dataset ( "samples" , ( self . params . nquartets , 4 ) , dtype = np . uint16 , chunks = ( self . _chunksize , 4 ) , compression = 'gzip' )
    # # fill chunksize at a time for efficiency
    i = 0
    empty = set ( )
    edge_targeted = 0
    random_target = 0
    # # keep filling quartets until nquartets are sampled
    while i < self . params . nquartets :
        qdat = [ ]
        # # keep filling this chunk until its full
        while len ( qdat ) < self . _chunksize : # # grab the next iterator
            qiter = qiters . next ( )
            cycle = cycler . next ( )
            # # sample from iterator
            try :
                qrtsamp = qiter . next ( )
                if tuple ( qrtsamp ) not in sampled :
                    qdat . append ( qrtsamp )
                    sampled . add ( qrtsamp )
                    edge_targeted += 1
                # else :
                # print ( ' repeat ' )
            # # unless iterator is empty , then skip it
            except StopIteration :
                empty . add ( cycle )
            # # break when all edge samplers are empty
            if len ( empty ) == nsaturation :
                break
        # # if array is not full then add random samples
        while len ( qdat ) < self . _chunksize :
            qrtsamp = random_combination ( range ( len ( self . samples ) ) , 4 )
            if tuple ( qrtsamp ) not in sampled :
                qdat . append ( qrtsamp )
                sampled . add ( qrtsamp )
                random_target += 1
        # # stick chunk into h5 array
        dat = np . array ( qdat , dtype = np . uint16 )
        io5 [ "samples" ] [ i : i + self . _chunksize ] = dat [ : io5 [ "samples" ] . shape [ 0 ] - i ]
        i += self . _chunksize
    print ( "  equal sampling: {} edge quartets, {} random quartets " . format ( edge_targeted , random_target ) )


signature:
def _run_qmc ( self , boot ) :
    """runs quartet max - cut on a quartets file"""

implementation:
# # convert to txt file for wQMC
self . _tmp = os . path . join ( self . dirs , ".tmpwtre" )
cmd = [ ip . bins . qmc , "qrtt=" + self . files . qdump , "otre=" + self . _tmp ]
# # run them
proc = subprocess . Popen ( cmd , stderr = subprocess . STDOUT , stdout = subprocess . PIPE )
res = proc . communicate ( )
if proc . returncode : # LOGGER . error ( " Error in QMC : \ n ( { } ) . " . format ( res ) )
    LOGGER . error ( res )
    raise IPyradWarningExit ( res [ 1 ] )
# # read in the tmp files since qmc does not pipe
with open ( self . _tmp ) as intree : # # convert int names back to str names renamer returns a newick str
# tmp = toytree . tree ( intree . read ( ) . strip ( ) )
    tmp = ete3 . Tree ( intree . read ( ) . strip ( ) )
    tmpwtre = self . _renamer ( tmp )
# . tree )
# # save the tree
if boot :
    self . trees . boots = os . path . join ( self . dirs , self . name + ".boots" )
    with open ( self . trees . boots , 'a' ) as outboot :
        outboot . write ( tmpwtre + "\n" )
else :
    self . trees . tree = os . path . join ( self . dirs , self . name + ".tree" )
    with open ( self . trees . tree , 'w' ) as outtree :
        outtree . write ( tmpwtre )
# # save JSON file checkpoint
self . _save ( )


signature:
def _dump_qmc ( self ) :
    """Makes a reduced array that excludes quartets with no information and
 prints the quartets and weights to a file formatted for wQMC"""

implementation:
# # open the h5 database
io5 = h5py . File ( self . database . output , 'r' )
# # create an output file for writing
self . files . qdump = os . path . join ( self . dirs , self . name + ".quartets.txt" )
LOGGER . info ( "qdump file %s" , self . files . qdump )
outfile = open ( self . files . qdump , 'w' )
# # todo : should pull quarts order in randomly ? or doesn ' t matter ?
for idx in xrange ( 0 , self . params . nquartets , self . _chunksize ) : # # get mask of zero weight quartets
# mask = io5 [ " weights " ] [ idx : idx + self . chunksize ] ! = 0
# weight = io5 [ " weights " ] [ idx : idx + self . chunksize ] [ mask ]
# LOGGER . info ( " exluded = % s , mask shape % s " ,
# self . _ chunksize - mask . shape [ 0 ] , mask . shape )
# LOGGER . info ( ' q shape % s ' , io5 [ " quartets " ] [ idx : idx + self . _ chunksize ] . shape )
    masked_quartets = io5 [ "quartets" ] [ idx : idx + self . _chunksize , : ]
    # [ mask , : ]
    quarts = [ list ( j ) for j in masked_quartets ]
    # # format and print
    # chunk = [ " { } , { } | { } , { } : { } " . format ( * i + [ j ] ) for i , j \
    # in zip ( quarts , weight ) ]
    chunk = [ "{},{}|{},{}" . format ( * i ) for i in quarts ]
    outfile . write ( "\n" . join ( chunk ) + "\n" )
# # close output file and h5 database
outfile . close ( )
io5 . close ( )


signature:
def _renamer ( self , tre ) :
    """renames newick from numbers to sample names"""

implementation:
# # get the tre with numbered tree tip labels
names = tre . get_leaves ( )
# # replace numbered names with snames
for name in names :
    name . name = self . samples [ int ( name . name ) ]
# # return with only topology and leaf labels
return tre . write ( format = 9 )


signature:
def _finalize_stats ( self , ipyclient ) :
    """write final tree files"""

implementation:
# # print stats file location :
# print ( STATSOUT . format ( opr ( self . files . stats ) ) )
# # print finished tree information - - - - -
print ( FINALTREES . format ( opr ( self . trees . tree ) ) )
# # print bootstrap information - - - - -
if self . params . nboots : # # get consensus , map values to tree edges , record stats file
    self . _compute_tree_stats ( ipyclient )
    # # print bootstrap info
    print ( BOOTTREES . format ( opr ( self . trees . cons ) , opr ( self . trees . boots ) ) )
# # print the ASCII tree only if its small
if len ( self . samples ) < 20 :
    if self . params . nboots :
        wctre = ete3 . Tree ( self . trees . cons , format = 0 )
        wctre . ladderize ( )
        print ( wctre . get_ascii ( show_internal = True , attributes = [ "dist" , "name" ] ) )
        print ( "" )
    else :
        qtre = ete3 . Tree ( self . trees . tree , format = 0 )
        qtre . ladderize ( )
        # qtre = toytree . tree ( self . trees . tree , format = 0)
        # qtre . tree . unroot ( )
        print ( qtre . get_ascii ( ) )
        print ( "" )
# # print PDF filename & tips - - - - -
docslink = "https://toytree.readthedocs.io/"
citelink = "https://ipyrad.readthedocs.io/tetrad.html"
print ( LINKS . format ( docslink , citelink ) )


signature:
def _save ( self ) :
    """save a JSON file representation of Tetrad Class for checkpoint"""

implementation:
# # save each attribute as dict
fulldict = copy . deepcopy ( self . __dict__ )
for i , j in fulldict . items ( ) :
    if isinstance ( j , Params ) :
        fulldict [ i ] = j . __dict__
fulldumps = json . dumps ( fulldict , sort_keys = False , indent = 4 , separators = ( "," , ":" ) , )
# # save to file , make dir if it wasn ' t made earlier
assemblypath = os . path . join ( self . dirs , self . name + ".tet.json" )
if not os . path . exists ( self . dirs ) :
    os . mkdir ( self . dirs )
# # protect save from interruption
done = 0
while not done :
    try :
        with open ( assemblypath , 'w' ) as jout :
            jout . write ( fulldumps )
        done = 1
    except ( KeyboardInterrupt , SystemExit ) :
        print ( '.' )
        continue


signature:
def _insert_to_array ( self , start , results ) :
    """inputs results from workers into hdf4 array"""

implementation:
qrts , wgts , qsts = results
# qrts , wgts = results
# print ( qrts )
with h5py . File ( self . database . output , 'r+' ) as out :
    chunk = self . _chunksize
    out [ 'quartets' ] [ start : start + chunk ] = qrts
    # # out [ ' weights ' ] [ start : start + chunk ] = wgts
    # # entered as 0 - indexed !
    if self . checkpoint . boots :
        key = "qboots/b{}" . format ( self . checkpoint . boots - 1 )
        out [ key ] [ start : start + chunk ] = qsts
    else :
        out [ "qstats" ] [ start : start + chunk ] = qsts