signature:
def add_axes_and_nodes ( self ) :
    """Adds the axes ( i . e . 2 or 3 axes , not to be confused with matplotlib
 axes ) and the nodes that belong to each axis ."""

implementation:
for i , ( group , nodelist ) in enumerate ( self . nodes . items ( ) ) :
    theta = self . group_theta ( group )
    if self . has_edge_within_group ( group ) :
        theta = theta - self . minor_angle
        self . plot_nodes ( nodelist , theta , group )
        theta = theta + 2 * self . minor_angle
        self . plot_nodes ( nodelist , theta , group )
    else :
        self . plot_nodes ( nodelist , theta , group )


signature:
def find_node_group_membership ( self , node ) :
    """Identifies the group for which a node belongs to ."""

implementation:
for group , nodelist in self . nodes . items ( ) :
    if node in nodelist :
        return group


signature:
def get_idx ( self , node ) :
    """Finds the index of the node in the sorted list ."""

implementation:
group = self . find_node_group_membership ( node )
return self . nodes [ group ] . index ( node )


signature:
def node_radius ( self , node ) :
    """Computes the radial position of the node ."""

implementation: return self . get_idx ( node ) * self . scale + self . internal_radius


signature:
def node_theta ( self , node ) :
    """Convenience function to find the node ' s theta angle ."""

implementation:
group = self . find_node_group_membership ( node )
return self . group_theta ( group )


signature:
def draw_edge ( self , n1 , n2 , d , group ) :
    """Renders the given edge ( n1 , n2 ) to the plot ."""

implementation:
start_radius = self . node_radius ( n1 )
start_theta = self . node_theta ( n1 )
end_radius = self . node_radius ( n2 )
end_theta = self . node_theta ( n2 )
start_theta , end_theta = self . correct_angles ( start_theta , end_theta )
start_theta , end_theta = self . adjust_angles ( n1 , start_theta , n2 , end_theta )
middle1_radius = np . min ( [ start_radius , end_radius ] )
middle2_radius = np . max ( [ start_radius , end_radius ] )
if start_radius > end_radius :
    middle1_radius , middle2_radius = middle2_radius , middle1_radius
middle1_theta = np . mean ( [ start_theta , end_theta ] )
middle2_theta = np . mean ( [ start_theta , end_theta ] )
startx , starty = get_cartesian ( start_radius , start_theta )
middle1x , middle1y = get_cartesian ( middle1_radius , middle1_theta )
middle2x , middle2y = get_cartesian ( middle2_radius , middle2_theta )
# middlex , middley = get _ cartesian ( middle _ radius , middle _ theta )
endx , endy = get_cartesian ( end_radius , end_theta )
verts = [ ( startx , starty ) , ( middle1x , middle1y ) , ( middle2x , middle2y ) , ( endx , endy ) ]
codes = [ Path . MOVETO , Path . CURVE4 , Path . CURVE4 , Path . CURVE4 ]
path = Path ( verts , codes )
if self . edge_colormap is None :
    edgecolor = 'black'
else :
    edgecolor = self . edge_colormap [ group ]
patch = patches . PathPatch ( path , lw = self . linewidth , facecolor = 'none' , edgecolor = edgecolor , alpha = 0.3 )
self . ax . add_patch ( patch )


signature:
def add_edges ( self ) :
    """Draws all of the edges in the graph ."""

implementation:
for group , edgelist in self . edges . items ( ) :
    for ( u , v , d ) in edgelist :
        self . draw_edge ( u , v , d , group )


signature:
def draw ( self ) :
    """The master function that is called that draws everything ."""

implementation:
self . ax . set_xlim ( - self . plot_radius ( ) , self . plot_radius ( ) )
self . ax . set_ylim ( - self . plot_radius ( ) , self . plot_radius ( ) )
self . add_axes_and_nodes ( )
self . add_edges ( )
self . ax . axis ( 'off' )


signature:
def adjust_angles ( self , start_node , start_angle , end_node , end_angle ) :
    """This function adjusts the start and end angles to correct for
 duplicated axes ."""

implementation:
start_group = self . find_node_group_membership ( start_node )
end_group = self . find_node_group_membership ( end_node )
if start_group == 0 and end_group == len ( self . nodes . keys ( ) ) - 1 :
    if self . has_edge_within_group ( start_group ) :
        start_angle = correct_negative_angle ( start_angle - self . minor_angle )
    if self . has_edge_within_group ( end_group ) :
        end_angle = correct_negative_angle ( end_angle + self . minor_angle )
elif start_group == len ( self . nodes . keys ( ) ) - 1 and end_group == 0 :
    if self . has_edge_within_group ( start_group ) :
        start_angle = correct_negative_angle ( start_angle + self . minor_angle )
    if self . has_edge_within_group ( end_group ) :
        end_angle = correct_negative_angle ( end_angle - self . minor_angle )
elif start_group < end_group :
    if self . has_edge_within_group ( end_group ) :
        end_angle = correct_negative_angle ( end_angle - self . minor_angle )
    if self . has_edge_within_group ( start_group ) :
        start_angle = correct_negative_angle ( start_angle + self . minor_angle )
elif end_group < start_group :
    if self . has_edge_within_group ( start_group ) :
        start_angle = correct_negative_angle ( start_angle - self . minor_angle )
    if self . has_edge_within_group ( end_group ) :
        end_angle = correct_negative_angle ( end_angle + self . minor_angle )
return start_angle , end_angle


signature:
def correct_angles ( self , start_angle , end_angle ) :
    """This function corrects for the following problems in the edges :"""

implementation:
# Edges going the anti - clockwise direction involves angle = 0.
if start_angle == 0 and ( end_angle - start_angle > np . pi ) :
    start_angle = np . pi * 2
if end_angle == 0 and ( end_angle - start_angle < - np . pi ) :
    end_angle = np . pi * 2
# Case when start _ angle = = end _ angle :
if start_angle == end_angle :
    start_angle = start_angle - self . minor_angle
    end_angle = end_angle + self . minor_angle
return start_angle , end_angle


signature:
def mods_genre ( self ) :
    """Guesses an appropriate MODS XML genre type ."""

implementation:
type2genre = { 'conference' : 'conference publication' , 'book chapter' : 'bibliography' , 'unpublished' : 'article' }
tp = str ( self . type ) . lower ( )
return type2genre . get ( tp , tp )


signature:
def _produce_author_lists ( self ) :
    """Parse authors string to create lists of authors ."""

implementation:
# post - process author names
self . authors = self . authors . replace ( ', and ' , ', ' )
self . authors = self . authors . replace ( ',and ' , ', ' )
self . authors = self . authors . replace ( ' and ' , ', ' )
self . authors = self . authors . replace ( ';' , ',' )
# list of authors
self . authors_list = [ author . strip ( ) for author in self . authors . split ( ',' ) ]
# simplified representation of author names
self . authors_list_simple = [ ]
# author names represented as a tuple of given and family name
self . authors_list_split = [ ]
# tests if title already ends with a punctuation mark
self . title_ends_with_punct = self . title [ - 1 ] in [ '.' , '!' , '?' ] if len ( self . title ) > 0 else False
suffixes = [ 'I' , 'II' , 'III' , 'IV' , 'V' , 'VI' , 'VII' , 'VIII' , "Jr." , "Sr." ]
prefixes = [ 'Dr.' ]
prepositions = [ 'van' , 'von' , 'der' , 'de' , 'den' ]
# further post - process author names
for i , author in enumerate ( self . authors_list ) :
    if author == '' :
        continue
    names = author . split ( ' ' )
    # check if last string contains initials
    if ( len ( names [ - 1 ] ) <= 3 ) and names [ - 1 ] not in suffixes and all ( c in ascii_uppercase for c in names [ - 1 ] ) : # turn " Gauss CF " into " C . F . Gauss "
        names = [ c + '.' for c in names [ - 1 ] ] + names [ : - 1 ]
    # number of suffixes
    num_suffixes = 0
    for name in names [ : : - 1 ] :
        if name in suffixes :
            num_suffixes += 1
        else :
            break
    # abbreviate names
    for j , name in enumerate ( names [ : - 1 - num_suffixes ] ) : # don ' t try to abbreviate these
        if j == 0 and name in prefixes :
            continue
        if j > 0 and name in prepositions :
            continue
        if ( len ( name ) > 2 ) or ( len ( name ) and ( name [ - 1 ] != '.' ) ) :
            k = name . find ( '-' )
            if 0 < k + 1 < len ( name ) : # take care of dash
                names [ j ] = name [ 0 ] + '.-' + name [ k + 1 ] + '.'
            else :
                names [ j ] = name [ 0 ] + '.'
    if len ( names ) :
        self . authors_list [ i ] = ' ' . join ( names )
        # create simplified / normalized representation of author name
        if len ( names ) > 1 :
            for name in names [ 0 ] . split ( '-' ) :
                name_simple = self . simplify_name ( ' ' . join ( [ name , names [ - 1 ] ] ) )
                self . authors_list_simple . append ( name_simple )
        else :
            self . authors_list_simple . append ( self . simplify_name ( names [ 0 ] ) )
        # number of prepositions
        num_prepositions = 0
        for name in names :
            if name in prepositions :
                num_prepositions += 1
        # splitting point
        sp = 1 + num_suffixes + num_prepositions
        self . authors_list_split . append ( ( ' ' . join ( names [ : - sp ] ) , ' ' . join ( names [ - sp : ] ) ) )
# list of authors in BibTex format
self . authors_bibtex = ' and ' . join ( self . authors_list )
# overwrite authors string
if len ( self . authors_list ) > 2 :
    self . authors = ', and ' . join ( [ ', ' . join ( self . authors_list [ : - 1 ] ) , self . authors_list [ - 1 ] ] )
elif len ( self . authors_list ) > 1 :
    self . authors = ' and ' . join ( self . authors_list )
else :
    self . authors = self . authors_list [ 0 ]


signature:
def get_publications ( context , template = 'publications/publications.html' ) :
    """Get all publications ."""

implementation:
types = Type . objects . filter ( hidden = False )
publications = Publication . objects . select_related ( )
publications = publications . filter ( external = False , type__in = types )
publications = publications . order_by ( '-year' , '-month' , '-id' )
if not publications :
    return ''
# load custom links and files
populate ( publications )
return render_template ( template , context [ 'request' ] , { 'publications' : publications } )


signature:
def get_publication ( context , id ) :
    """Get a single publication ."""

implementation:
pbl = Publication . objects . filter ( pk = int ( id ) )
if len ( pbl ) < 1 :
    return ''
pbl [ 0 ] . links = pbl [ 0 ] . customlink_set . all ( )
pbl [ 0 ] . files = pbl [ 0 ] . customfile_set . all ( )
return render_template ( 'publications/publication.html' , context [ 'request' ] , { 'publication' : pbl [ 0 ] } )


signature:
def get_publication_list ( context , list , template = 'publications/publications.html' ) :
    """Get a publication list ."""

implementation:
list = List . objects . filter ( list__iexact = list )
if not list :
    return ''
list = list [ 0 ]
publications = list . publication_set . all ( )
publications = publications . order_by ( '-year' , '-month' , '-id' )
if not publications :
    return ''
# load custom links and files
populate ( publications )
return render_template ( template , context [ 'request' ] , { 'list' : list , 'publications' : publications } )


signature:
def tex_parse ( string ) :
    """Renders some basic TeX math to HTML ."""


signature:
def parse ( string ) :
    """Takes a string in BibTex format and returns a list of BibTex entries , where
 each entry is a dictionary containing the entries ' key - value pairs .
 @ type string : string
 @ param string : bibliography in BibTex format
 @ rtype : list
 @ return : a list of dictionaries representing a bibliography"""

implementation:
# bibliography
bib = [ ]
# make sure we are dealing with unicode strings
if not isinstance ( string , six . text_type ) :
    string = string . decode ( 'utf-8' )
# replace special characters
for key , value in special_chars :
    string = string . replace ( key , value )
string = re . sub ( r'\\[cuHvs]{?([a-zA-Z])}?' , r'\1' , string )
# split into BibTex entries
entries = re . findall ( r'(?u)@(\w+)[ \t]?{[ \t]*([^,\s]*)[ \t]*,?\s*((?:[^=,\s]+\s*\=\s*(?:"[^"]*"|{(?:[^{}]*|{[^{}]*})*}|[^,}]*),?\s*?)+)\s*}' , string )
for entry in entries : # parse entry
    pairs = re . findall ( r'(?u)([^=,\s]+)\s*\=\s*("[^"]*"|{(?:[^{}]*|{[^{}]*})*}|[^,]*)' , entry [ 2 ] )
    # add to bibliography
    bib . append ( { 'type' : entry [ 0 ] . lower ( ) , 'key' : entry [ 1 ] } )
    for key , value in pairs : # post - process key and value
        key = key . lower ( )
        if value and value [ 0 ] == '"' and value [ - 1 ] == '"' :
            value = value [ 1 : - 1 ]
        if value and value [ 0 ] == '{' and value [ - 1 ] == '}' :
            value = value [ 1 : - 1 ]
        if key not in [ 'booktitle' , 'title' ] :
            value = value . replace ( '}' , '' ) . replace ( '{' , '' )
        else :
            if value . startswith ( '{' ) and value . endswith ( '}' ) :
                value = value [ 1 : ]
                value = value [ : - 1 ]
        value = value . strip ( )
        value = re . sub ( r'\s+' , ' ' , value )
        # store pair in bibliography
        bib [ - 1 ] [ key ] = value
return bib


signature:
def swap ( self , qs ) :
    """Swap the positions of this object with a reference object ."""

implementation:
try :
    replacement = qs [ 0 ]
except IndexError : # already first / last
    return
if not self . _valid_ordering_reference ( replacement ) :
    raise ValueError ( "%r can only be swapped with instances of %r which %s equals %r." % ( self , self . __class__ , self . order_with_respect_to , self . _get_order_with_respect_to ( ) ) )
self . order , replacement . order = replacement . order , self . order
self . save ( )
replacement . save ( )


signature:
def up ( self ) :
    """Move this object up one position ."""

implementation: self . swap ( self . get_ordering_queryset ( ) . filter ( order__lt = self . order ) . order_by ( '-order' ) )


signature:
def down ( self ) :
    """Move this object down one position ."""

implementation: self . swap ( self . get_ordering_queryset ( ) . filter ( order__gt = self . order ) )


signature:
def to ( self , order ) :
    """Move object to a certain position , updating all affected objects to move accordingly up or down ."""

implementation:
if order is None or self . order == order : # object is already at desired position
    return
qs = self . get_ordering_queryset ( )
if self . order > order :
    qs . filter ( order__lt = self . order , order__gte = order ) . update ( order = F ( 'order' ) + 1 )
else :
    qs . filter ( order__gt = self . order , order__lte = order ) . update ( order = F ( 'order' ) - 1 )
self . order = order
self . save ( )


signature:
def above ( self , ref ) :
    """Move this object above the referenced object ."""

implementation:
if not self . _valid_ordering_reference ( ref ) :
    raise ValueError ( "%r can only be moved above instances of %r which %s equals %r." % ( self , self . __class__ , self . order_with_respect_to , self . _get_order_with_respect_to ( ) ) )
if self . order == ref . order :
    return
if self . order > ref . order :
    o = ref . order
else :
    o = self . get_ordering_queryset ( ) . filter ( order__lt = ref . order ) . aggregate ( Max ( 'order' ) ) . get ( 'order__max' ) or 0
self . to ( o )


signature:
def below ( self , ref ) :
    """Move this object below the referenced object ."""

implementation:
if not self . _valid_ordering_reference ( ref ) :
    raise ValueError ( "%r can only be moved below instances of %r which %s equals %r." % ( self , self . __class__ , self . order_with_respect_to , self . _get_order_with_respect_to ( ) ) )
if self . order == ref . order :
    return
if self . order > ref . order :
    o = self . get_ordering_queryset ( ) . filter ( order__gt = ref . order ) . aggregate ( Min ( 'order' ) ) . get ( 'order__min' ) or 0
else :
    o = ref . order
self . to ( o )


signature:
def top ( self ) :
    """Move this object to the top of the ordered stack ."""

implementation:
o = self . get_ordering_queryset ( ) . aggregate ( Min ( 'order' ) ) . get ( 'order__min' )
self . to ( o )


signature:
def bottom ( self ) :
    """Move this object to the bottom of the ordered stack ."""

implementation:
o = self . get_ordering_queryset ( ) . aggregate ( Max ( 'order' ) ) . get ( 'order__max' )
self . to ( o )


signature:
def unapi ( request ) :
    """This view implements unAPI 1.0 ( see http : / / unapi . info ) ."""

implementation:
id = request . GET . get ( 'id' )
format = request . GET . get ( 'format' )
if format is not None :
    try :
        publications = Publication . objects . filter ( pk = int ( id ) )
        if not publications :
            raise ValueError
    except ValueError : # invalid id
        return HttpResponse ( '\n' . join ( [ '<?xml version="1.0" encoding="UTF-8"?>' , '<error>Invalid ID.</error>' ] ) , content_type = "application/xml" , status = 404 )
    if format == 'bibtex' : # return BibTex encoded publication
        return render ( request , 'publications/publication.bib' , { 'publication' : publications [ 0 ] } , content_type = 'text/x-bibtex; charset=UTF-8' )
    if format == 'mods' : # return MODS encoded publication
        return render ( request , 'publications/publications.mods' , { 'publications' : publications } , content_type = 'application/xml; charset=UTF-8' )
    if format == 'ris' : # return MODS encoded publication
        return render ( request , 'publications/publications.ris' , { 'publications' : publications } , content_type = 'application/x-research-info-systems; charset=UTF-8' )
    # invalid format
    return HttpResponse ( '\n' . join ( [ '<?xml version="1.0" encoding="UTF-8"?>' , '<error>Invalid format.</error>' ] ) , content_type = "application/xml" , status = 406 )
if id is not None :
    return HttpResponse ( '\n' . join ( [ '<?xml version="1.0" encoding="UTF-8"?>' , '<formats id="{0}">' . format ( id ) , '<format name="bibtex" type="text/x-bibtex" />' , '<format name="ris" type="application/x-research-info-systems" />' , '<format name="mods" type="application/xml" />' , '</formats>' ] ) , content_type = "application/xml" )
return HttpResponse ( '\n' . join ( [ '<?xml version="1.0" encoding="UTF-8"?>' , '<formats>' , '<format name="bibtex" type="text/x-bibtex" />' , '<format name="ris" type="application/x-research-info-systems" />' , '<format name="mods" type="application/xml" />' , '</formats>' ] ) , content_type = "application/xml" )


signature:
def populate ( publications ) :
    """Load custom links and files from database and attach to publications ."""

implementation:
customlinks = CustomLink . objects . filter ( publication__in = publications )
customfiles = CustomFile . objects . filter ( publication__in = publications )
publications_ = { }
for publication in publications :
    publication . links = [ ]
    publication . files = [ ]
    publications_ [ publication . id ] = publication
for link in customlinks :
    publications_ [ link . publication_id ] . links . append ( link )
for file in customfiles :
    publications_ [ file . publication_id ] . files . append ( file )


signature:
def worker ( self ) :
    """Calculates the quartet weights for the test at a random
 subsampled chunk of loci ."""

implementation:
# # subsample loci
fullseqs = self . sample_loci ( )
# # find all iterations of samples for this quartet
liters = itertools . product ( * self . imap . values ( ) )
# # run tree inference for each iteration of sampledict
hashval = uuid . uuid4 ( ) . hex
weights = [ ]
for ridx , lidx in enumerate ( liters ) : # # get subalignment for this iteration and make to nex
    a , b , c , d = lidx
    sub = { }
    for i in lidx :
        if self . rmap [ i ] == "p1" :
            sub [ "A" ] = fullseqs [ i ]
        elif self . rmap [ i ] == "p2" :
            sub [ "B" ] = fullseqs [ i ]
        elif self . rmap [ i ] == "p3" :
            sub [ "C" ] = fullseqs [ i ]
        else :
            sub [ "D" ] = fullseqs [ i ]
    # # write as nexus file
    nex = [ ]
    for tax in list ( "ABCD" ) :
        nex . append ( ">{}         {}" . format ( tax , sub [ tax ] ) )
    # # check for too much missing or lack of variants
    nsites , nvar = count_var ( nex )
    # # only run test if there ' s variation present
    if nvar > self . minsnps : # # format as nexus file
        nexus = "{} {}\n" . format ( 4 , len ( fullseqs [ a ] ) ) + "\n" . join ( nex )
        # # infer ML tree
        treeorder = self . run_tree_inference ( nexus , "{}.{}" . format ( hashval , ridx ) )
        # # add to list
        weights . append ( treeorder )
# # cleanup - remove all files with the hash val
rfiles = glob . glob ( os . path . join ( tempfile . tempdir , "*{}*" . format ( hashval ) ) )
for rfile in rfiles :
    if os . path . exists ( rfile ) :
        os . remove ( rfile )
# # return result as weights for the set topologies .
trees = [ "ABCD" , "ACBD" , "ADBC" ]
wdict = { i : float ( weights . count ( i ) ) / len ( weights ) for i in trees }
return wdict


signature:
def get_order ( tre ) :
    """return tree order"""

implementation:
anode = tre . tree & ">A"
sister = anode . get_sisters ( ) [ 0 ]
sisters = ( anode . name [ 1 : ] , sister . name [ 1 : ] )
others = [ i for i in list ( "ABCD" ) if i not in sisters ]
return sorted ( sisters ) + sorted ( others )


signature:
def count_var ( nex ) :
    """count number of sites with cov = 4 , and number of variable sites ."""

implementation:
arr = np . array ( [ list ( i . split ( ) [ - 1 ] ) for i in nex ] )
miss = np . any ( arr == "N" , axis = 0 )
nomiss = arr [ : , ~ miss ]
nsnps = np . invert ( np . all ( nomiss == nomiss [ 0 , : ] , axis = 0 ) ) . sum ( )
return nomiss . shape [ 1 ] , nsnps


signature:
def sample_loci ( self ) :
    """finds loci with sufficient sampling for this test"""

implementation:
# # store idx of passing loci
idxs = np . random . choice ( self . idxs , self . ntests )
# # open handle , make a proper generator to reduce mem
with open ( self . data ) as indata :
    liter = ( indata . read ( ) . strip ( ) . split ( "|\n" ) )
# # store data as dict
seqdata = { i : "" for i in self . samples }
# # put chunks into a list
for idx , loc in enumerate ( liter ) :
    if idx in idxs : # # parse chunk
        lines = loc . split ( "\n" ) [ : - 1 ]
        names = [ i . split ( ) [ 0 ] for i in lines ]
        seqs = [ i . split ( ) [ 1 ] for i in lines ]
        dd = { i : j for i , j in zip ( names , seqs ) }
        # # add data to concatenated seqdict
        for name in seqdata :
            if name in names :
                seqdata [ name ] += dd [ name ]
            else :
                seqdata [ name ] += "N" * len ( seqs [ 0 ] )
# # concatenate into a phylip file
return seqdata


signature:
def run_tree_inference ( self , nexus , idx ) :
    """Write nexus to tmpfile , runs phyml tree inference , and parses
 and returns the resulting tree ."""

implementation:
# # create a tmpdir for this test
tmpdir = tempfile . tempdir
tmpfile = os . path . join ( tempfile . NamedTemporaryFile ( delete = False , prefix = str ( idx ) , dir = tmpdir , ) )
# # write nexus to tmpfile
tmpfile . write ( nexus )
tmpfile . flush ( )
# # infer the tree
rax = raxml ( name = str ( idx ) , data = tmpfile . name , workdir = tmpdir , N = 1 , T = 2 )
rax . run ( force = True , block = True , quiet = True )
# # clean up
tmpfile . close ( )
# # return tree order
order = get_order ( toytree . tree ( rax . trees . bestTree ) )
return "" . join ( order )


signature:
def run ( self , ipyclient ) :
    """parallelize calls to worker function ."""

implementation:
# # connect to parallel client
lbview = ipyclient . load_balanced_view ( )
# # iterate over tests
asyncs = [ ]
for test in xrange ( self . ntests ) : # # submit jobs to run
    async = lbview . apply ( worker , self )
    asyncs . append ( async )
# # wait for jobs to finish
ipyclient . wait ( )
# # check for errors
for async in asyncs :
    if not async . successful ( ) :
        raise Exception ( "Error: {}" . format ( async . result ( ) ) )
# # return results as df
results = [ i . result ( ) for i in asyncs ]
self . results_table = pd . DataFrame ( results )


signature:
def plot ( self ) :
    """return a toyplot barplot of the results table ."""

implementation:
if self . results_table == None :
    return "no results found"
else :
    bb = self . results_table . sort_values ( by = [ "ABCD" , "ACBD" ] , ascending = [ False , True ] , )
    # # make a barplot
    import toyplot
    c = toyplot . Canvas ( width = 600 , height = 200 )
    a = c . cartesian ( )
    m = a . bars ( bb )
    return c , a , m


signature:
def plot ( self , pcs = [ 1 , 2 ] , ax = None , cmap = None , cdict = None , legend = True , title = None , outfile = None ) :
    """Do the PCA and plot it .
 Parameters
 pcs : list of ints
 ax : matplotlib axis
 cmap : matplotlib colormap
 cdict : dictionary mapping pop names to colors
 legend : boolean , whether or not to show the legend"""

implementation:
# # Specify which 2 pcs to plot , default is pc1 and pc2
pc1 = pcs [ 0 ] - 1
pc2 = pcs [ 1 ] - 1
if pc1 < 0 or pc2 > self . ncomponents - 1 :
    raise IPyradError ( "PCs are 1-indexed. 1 is min & {} is max" . format ( self . ncomponents ) )
# # Convert genotype data to allele count data
# # We do this here because we might want to try different ways
# # of accounting for missing data and " alt " allele counts treat
# # missing data as " ref "
allele_counts = self . genotypes . to_n_alt ( )
# # Actually do the pca
if self . ncomponents > len ( self . samples_vcforder ) :
    self . ncomponents = len ( self . samples_vcforder )
    print ( "  INFO: # PCs < # samples. Forcing # PCs = {}" . format ( self . ncomponents ) )
coords , model = allel . pca ( allele_counts , n_components = self . ncomponents , scaler = 'patterson' )
self . pcs = pd . DataFrame ( coords , index = self . samples_vcforder , columns = [ "PC{}" . format ( x ) for x in range ( 1 , self . ncomponents + 1 ) ] )
# # Just allow folks to pass in the name of the cmap they want to use
if isinstance ( cmap , str ) :
    try :
        cmap = cm . get_cmap ( cmap )
    except :
        raise IPyradError ( "  Bad cmap value: {}" . format ( cmap ) )
if not cmap and not cdict :
    if not self . quiet :
        print ( "  Using default cmap: Spectral" )
    cmap = cm . get_cmap ( 'Spectral' )
if cmap :
    if cdict :
        print ( "  Passing in both cmap and cdict defaults to using the cmap value." )
    popcolors = cmap ( np . arange ( len ( self . pops ) ) / len ( self . pops ) )
    cdict = { i : j for i , j in zip ( self . pops . keys ( ) , popcolors ) }
fig = ""
if not ax :
    fig = plt . figure ( figsize = ( 6 , 5 ) )
    ax = fig . add_subplot ( 1 , 1 , 1 )
x = coords [ : , pc1 ]
y = coords [ : , pc2 ]
for pop in self . pops : # # Don ' t include pops with no samples , it makes the legend look stupid
# # TODO : This doesn ' t prevent empty pops from showing up in the legend for some reason .
    if len ( self . pops [ pop ] ) > 0 :
        mask = np . isin ( self . samples_vcforder , self . pops [ pop ] )
        ax . plot ( x [ mask ] , y [ mask ] , marker = 'o' , linestyle = ' ' , color = cdict [ pop ] , label = pop , markersize = 6 , mec = 'k' , mew = .5 )
ax . set_xlabel ( 'PC%s (%.1f%%)' % ( pc1 + 1 , model . explained_variance_ratio_ [ pc1 ] * 100 ) )
ax . set_ylabel ( 'PC%s (%.1f%%)' % ( pc2 + 1 , model . explained_variance_ratio_ [ pc2 ] * 100 ) )
if legend :
    ax . legend ( bbox_to_anchor = ( 1 , 1 ) , loc = 'upper left' )
if fig :
    fig . tight_layout ( )
if title :
    ax . set_title ( title )
if outfile :
    try :
        plt . savefig ( outfile , format = "png" , bbox_inches = "tight" )
    except :
        print ( "  Saving pca.plot() failed to save figure to {}" . format ( outfile ) )
return ax


signature:
def plot_pairwise_dist ( self , labels = None , ax = None , cmap = None , cdict = None , metric = "euclidean" ) :
    """Plot pairwise distances between all samples
 labels : bool or list
 by default labels aren ' t included . If labels = = True , then labels are read in
 from the vcf file . Alternatively , labels can be passed in as a list , should
 be same length as the number of samples ."""

implementation:
allele_counts = self . genotypes . to_n_alt ( )
dist = allel . pairwise_distance ( allele_counts , metric = metric )
if not ax :
    fig = plt . figure ( figsize = ( 5 , 5 ) )
    ax = fig . add_subplot ( 1 , 1 , 1 )
if isinstance ( labels , bool ) :
    if labels :
        labels = list ( self . samples_vcforder )
elif isinstance ( labels , type ( None ) ) :
    pass
else : # # If not bool or None ( default ) , then check to make sure the list passed in
# # is the right length
    if not len ( labels ) == len ( self . samples_vcforder ) :
        raise IPyradError ( LABELS_LENGTH_ERROR . format ( len ( labels ) , len ( self . samples_vcforder ) ) )
allel . plot . pairwise_distance ( dist , labels = labels , ax = ax , colorbar = False )


signature:
def copy ( self ) :
    """returns a copy of the pca analysis object"""

implementation:
cp = copy . deepcopy ( self )
cp . genotypes = allel . GenotypeArray ( self . genotypes , copy = True )
return cp


signature:
def loci2cf ( name , locifile , popdict , wdir = None , ipyclient = None ) :
    """Convert ipyrad . loci file to an iqtree - pomo ' counts ' file
 Parameters :
 name :
 A prefix name for output files that will be produced
 locifile :
 A . loci file produced by ipyrad .
 popdict :
 A python dictionary grouping Clade names to Sample names .
 Example : { " A " : [ ' a ' , ' b ' , ' c ' ] , " B " : [ ' d ' , ' e ' , ' f ' ] }
 ipyclient :
 If you pass it an ipyclient it will distribute work over
 remote engines , otherwise we use multiprocessing ( todo ) ."""

implementation:
# # working directory , make sure it exists
if wdir :
    wdir = os . path . abspath ( wdir )
    if not os . path . exists ( wdir ) :
        raise IPyradWarningExit ( " working directory (wdir) does not exist" )
else :
    wdir = os . path . curdir
# # output file path
name = name . rsplit ( ".cf" ) [ 0 ]
outfile = os . path . join ( wdir , "{}.cf" . format ( name ) )
out = open ( outfile , 'w' )
# # parse loci file
with open ( locifile ) as inloc :
    loci = inloc . read ( ) . strip ( ) . split ( "|\n" )
# # get all names
names = list ( itertools . chain ( * popdict . values ( ) ) )
popkeys = sorted ( popdict . keys ( ) )
# # count nsites
nsites = sum ( len ( loc . split ( "\n" ) [ 0 ] . split ( ) [ 1 ] ) for loc in loci [ : ] )
# # print the header
out . write ( HEADER . format ( ** { "NPOP" : len ( popdict ) , "NSITES" : nsites , "VTAXA" : "\t" . join ( popkeys ) } ) )
# # build print string
outstr = "chr{:<8}  {:<4}  "
for cidx in xrange ( len ( popkeys ) ) :
    outstr += "{:<8}  "
toprint = [ ]
for idx in xrange ( len ( loci ) ) :
    dat = loci [ idx ] . split ( "\n" )
    seqs = np . array ( [ list ( i . split ( ) [ 1 ] ) for i in dat [ : - 1 ] ] )
    names = [ i . split ( ) [ 0 ] for i in dat [ : - 1 ] ]
    data = np . zeros ( ( seqs . shape [ 1 ] , len ( popkeys ) , 4 ) , dtype = np . uint16 )
    for sidx in xrange ( seqs . shape [ 1 ] ) :
        for cidx in xrange ( len ( popkeys ) ) :
            for name in popdict [ popkeys [ cidx ] ] :
                if name in names :
                    base = seqs [ names . index ( name ) , sidx ]
                    if base in list ( "ACGT" ) :
                        data [ sidx , cidx , BASE2IDX [ base ] ] += 2
                    elif base in list ( "RSYMKW" ) :
                        base1 , base2 = AMBIGS [ base ]
                        data [ sidx , cidx , BASE2IDX [ base1 ] ] += 1
                        data [ sidx , cidx , BASE2IDX [ base2 ] ] += 1
        # # print string for one locus
        sdat = [ "," . join ( [ str ( i ) for i in i . tolist ( ) ] ) for i in data [ sidx ] ]
        # print outstr . format ( idx + 1 , sidx + 1 , * sdat )
        toprint . append ( outstr . format ( idx + 1 , sidx + 1 , * sdat ) )
    # # if 10K loci , then print and clear
    if not idx % 10000 :
        out . write ( "\n" . join ( toprint ) + "\n" )
        toprint = [ ]
# # close handle
out . write ( "\n" . join ( toprint ) + "\n" )
out . close ( )


signature:
def loci2migrate ( name , locifile , popdict , mindict = 1 ) :
    """A function to build an input file for the program migrate from an ipyrad
 . loci file , and a dictionary grouping Samples into populations .
 Parameters :
 name : ( str )
 The name prefix for the migrate formatted output file .
 locifile : ( str )
 The path to the . loci file produced by ipyrad .
 popdict : ( dict )
 A Python dictionary grouping Samples into Populations .
 Examples :
 You can create the population dictionary by hand , and pass in the path
 to your . loci file as a string .
 > > popdict = { ' A ' : [ ' a ' , ' b ' , ' c ' ] , ' B ' : [ ' d ' , ' e ' , ' f ' ] }
 > > loci2migrate ( " outfile . migrate " , " . / mydata . loci " , popdict )
 Or , if you load your ipyrad . Assembly object from it ' s JSON file , you can
 access the loci file path and population information from there directly .
 > > data = ip . load _ json ( " mydata . json " )
 > > loci2migrate ( " outfile . migrate " , data . outfiles . loci , data . populations )"""

implementation:
# # I / O
outfile = open ( name + ".migrate" , 'w' )
infile = open ( locifile , 'r' )
# # minhits dictionary can be an int ( all same ) or a dictionary ( set each )
if isinstance ( mindict , int ) :
    mindict = { pop : mindict for pop in popdict }
else :
    mindict = mindict
# # filter data to only the loci that have data for mindict setting
keep = [ ]
MINS = zip ( taxa . keys ( ) , minhits )
# # read in data to sample names
loci = infile . read ( ) . strip ( ) . split ( "|" ) [ : - 1 ]
for loc in loci :
    samps = [ i . split ( ) [ 0 ] . replace ( ">" , "" ) for i in loc . split ( "\n" ) if ">" in i ]
    # # filter for coverage
    GG = [ ]
    for group , mins in MINS :
        GG . append ( sum ( [ i in samps for i in taxa [ group ] ] ) >= int ( mins ) )
    if all ( GG ) :
        keep . append ( loc )
# # print data to file
print >> outfile , len ( taxa ) , len ( keep ) , "( npops nloci for data set" , data . name + ".loci" , ")"
# # print all data for each population at a time
done = 0
for group in taxa : # # print a list of lengths of each locus
    if not done :
        loclens = [ len ( loc . split ( "\n" ) [ 1 ] . split ( ) [ - 1 ] . replace ( "x" , "n" ) . replace ( "n" , "" ) ) for loc in keep ]
        print >> outfile , " " . join ( map ( str , loclens ) )
        done += 1
    # # print a list of number of individuals in each locus
    indslist = [ ]
    for loc in keep :
        samps = [ i . split ( ) [ 0 ] . replace ( ">" , "" ) for i in loc . split ( "\n" ) if ">" in i ]
        inds = sum ( [ i in samps for i in taxa [ group ] ] )
        indslist . append ( inds )
    print >> outfile , " " . join ( map ( str , indslist ) ) , group
    # # print sample id , spaces , and sequence data
    # for loc in range ( len ( keep ) ) :
    for loc in range ( len ( keep ) ) :
        seqs = [ i . split ( ) [ - 1 ] for i in keep [ loc ] . split ( "\n" ) if i . split ( ) [ 0 ] . replace ( ">" , "" ) in taxa [ group ] ]
        for i in range ( len ( seqs ) ) :
            print >> outfile , group [ 0 : 8 ] + "_" + str ( i ) + ( " " * ( 10 - len ( group [ 0 : 8 ] + "_" + str ( i ) ) ) ) + seqs [ i ] . replace ( "x" , "n" ) . replace ( "n" , "" )
outfile . close ( )


signature:
def update ( assembly , idict , count ) :
    """updates dictionary with the next . 5M reads from the super long string
 phylip file . Makes for faster reading ."""

implementation:
data = iter ( open ( os . path . join ( assembly . dirs . outfiles , assembly . name + ".phy" ) , 'r' ) )
ntax , nchar = data . next ( ) . strip ( ) . split ( )
# # read in max N bp at a time
for line in data :
    tax , seq = line . strip ( ) . split ( )
    idict [ tax ] = idict [ tax ] [ 100000 : ]
    idict [ tax ] += seq [ count : count + 100000 ]
del line
return idict


signature:
def makephy ( data , samples , longname ) :
    """builds phy output . If large files writes 50000 loci at a time to tmp
 files and rebuilds at the end"""

implementation:
# # order names
names = [ i . name for i in samples ]
names . sort ( )
# # read in loci file
locifile = os . path . join ( data . dirs . outfiles , data . name + ".loci" )
locus = iter ( open ( locifile , 'rb' ) )
# # dict for saving the full matrix
fdict = { name : [ ] for name in names }
# # list for saving locus number and locus range for partitions
partitions = [ ]
initial_pos = 1
# # remove empty column sites and append edited seqs to dict F
done = 0
nloci = 0
nbases = 0
# # TODO : This should be fixed . it cycles through reading each locus
# # until nloci is less than this large number . It should really just
# # read to the end of the file , so it ' ll do all loci no matter how
# # many there are .
while nloci < 5000000 :
    seqs = [ ]
    # arrayed = np . array ( [ ] )
    anames = [ ]
    while 1 : # # get next locus
        try :
            samp = locus . next ( )
        except StopIteration :
            done = 1
            break
        if "//" in samp :
            nloci += 1
            break
        else :
            try :
                name , seq = samp . split ( )
            except ValueError :
                print samp
            anames . append ( name [ 1 : ] )
            seqs . append ( seq . strip ( ) )
    # # reset
    arrayed = np . array ( [ list ( i ) for i in seqs ] )
    if done :
        break
    # # create mask for columns that are empty or
    # # that are paired - end separators ( compatible w / pyrad v2 and v3)
    # mask = [ i for i in range ( len ( arrayed . T ) ) if np . any ( [
    # # still surely a better way to vectorize this . . .
    mask = [ i for i in arrayed . T if any ( [ j not in list ( "-Nn" ) for j in i ] ) ]
    masked = np . dstack ( mask ) [ 0 ]
    # # partition information
    loc_name = "p" + str ( nloci )
    loc_range = str ( initial_pos ) + "-" + str ( len ( masked [ 0 ] ) + initial_pos - 1 )
    initial_pos += len ( masked [ 0 ] )
    partitions . append ( loc_name + "=" + loc_range )
    # # uncomment to print block info ( used to partition by locus )
    # blockend + = minray
    # print blockend ,
    # print loc
    # print arrayed
    # # append data to dict
    for name in names :
        if name in anames : # fdict [ name ] . append ( arrayed [ anames . index ( name ) , mask ] . tostring ( ) )
            fdict [ name ] . append ( masked [ anames . index ( name ) , : ] . tostring ( ) )
        else :
            fdict [ name ] . append ( "N" * masked . shape [ 1 ] )
            # fdict [ name ] . append ( " N " * len ( arrayed [ 0 , mask ] ) )
    # # add len to total length
    nbases += len ( fdict [ name ] [ - 1 ] )
    # # after x iterations tmp pickle fdict ?
    if not nloci % 1e4 : # # concat strings
        for name in fdict :
            with open ( os . path . join ( assembly . dirs . outfiles , "tmp" , "{}_{}.phy.tmp" . format ( name , nloci ) ) , 'wb' ) as wout :
                wout . write ( "" . join ( fdict [ name ] ) )
        del fdict
        fdict = { name : [ ] for name in names }
# # print out . PHY file , if really big , pull form multiple tmp pickle
superout = open ( os . path . join ( assembly . dirs . outfiles , assembly . name + ".phy" ) , 'wb' )
print >> superout , len ( names ) , nbases
if nloci < 1e4 :
    for name in names :
        print >> superout , name + ( " " * ( ( longname + 3 ) - len ( name ) ) ) + "" . join ( fdict [ name ] )
else :
    for name in names :
        superout . write ( "{}{}{}" . format ( name , " " * ( ( longname + 3 ) - len ( name ) ) , "" . join ( fdict [ name ] ) ) )
        tmpfiles = glob . glob ( os . path . join ( assembly . dirs . outfiles , "tmp" , name + "*.phy.tmp" ) )
        tmpfiles . sort ( )
        for tmpf in tmpfiles :
            with open ( tmpf , 'rb' ) as tmpin :
                superout . write ( tmpin . read ( ) )
            os . remove ( tmpf )
        superout . write ( "\n" )
superout . close ( )
raxml_part_out = open ( os . path . join ( assembly . dirs . outfiles , assembly . name + ".phy.partitions" ) , 'w' )
for partition in partitions :
    print >> raxml_part_out , "DNA, %s" % ( partition )
raxml_part_out . close ( )
return partitions


signature:
def makenex ( assembly , names , longname , partitions ) :
    """PRINT NEXUS"""

implementation:
# # make nexus output
data = iter ( open ( os . path . join ( assembly . dirs . outfiles , assembly . name + ".phy" ) , 'r' ) )
nexout = open ( os . path . join ( assembly . dirs . outfiles , assembly . name + ".nex" ) , 'wb' )
ntax , nchar = data . next ( ) . strip ( ) . split ( " " )
print >> nexout , "#NEXUS"
print >> nexout , "BEGIN DATA;"
print >> nexout , "  DIMENSIONS NTAX=%s NCHAR=%s;" % ( ntax , nchar )
print >> nexout , "  FORMAT DATATYPE=DNA MISSING=N GAP=- INTERLEAVE=YES;"
print >> nexout , "  MATRIX"
idict = { }
# # read in max 1M bp at a time
for line in data :
    tax , seq = line . strip ( ) . split ( )
    idict [ tax ] = seq [ 0 : 100000 ]
del line
nameorder = idict . keys ( )
nameorder . sort ( )
n = 0
tempn = 0
sz = 100
while n < len ( seq ) :
    for tax in nameorder :
        print >> nexout , "  " + tax + " " * ( ( longname - len ( tax ) ) + 3 ) + idict [ tax ] [ tempn : tempn + sz ]
    n += sz
    tempn += sz
    print >> nexout , ""
    if not n % 100000 : # print idict [ tax ] [ tempn : tempn + sz ]
        idict = update ( assembly , idict , n )
        tempn -= 100000
print >> nexout , ';'
print >> nexout , 'END;'
# # # partitions info
print >> nexout , "BEGIN SETS;"
for partition in partitions :
    print >> nexout , "  CHARSET %s;" % ( partition )
print >> nexout , "END;"
nexout . close ( )


signature:
def make ( assembly , samples ) :
    """Make phylip and nexus formats . This is hackish since I ' m recycling the
 code whole - hog from pyrad V3 . Probably could be good to go back through
 and clean up the conversion code some time ."""

implementation:
# # get the longest name
longname = max ( [ len ( i ) for i in assembly . samples . keys ( ) ] )
names = [ i . name for i in samples ]
partitions = makephy ( assembly , samples , longname )
makenex ( assembly , names , longname , partitions )


signature:
def sample_cleanup ( data , sample ) :
    """Clean up a bunch of loose files ."""

implementation:
umap1file = os . path . join ( data . dirs . edits , sample . name + "-tmp-umap1.fastq" )
umap2file = os . path . join ( data . dirs . edits , sample . name + "-tmp-umap2.fastq" )
unmapped = os . path . join ( data . dirs . refmapping , sample . name + "-unmapped.bam" )
samplesam = os . path . join ( data . dirs . refmapping , sample . name + ".sam" )
split1 = os . path . join ( data . dirs . edits , sample . name + "-split1.fastq" )
split2 = os . path . join ( data . dirs . edits , sample . name + "-split2.fastq" )
refmap_derep = os . path . join ( data . dirs . edits , sample . name + "-refmap_derep.fastq" )
for f in [ umap1file , umap2file , unmapped , samplesam , split1 , split2 , refmap_derep ] :
    try :
        os . remove ( f )
    except :
        pass


signature:
def index_reference_sequence ( data , force = False ) :
    """Index the reference sequence , unless it already exists . Also make a mapping
 of scaffolds to index numbers for later user in steps 5-6."""

implementation:
# # get ref file from params
refseq_file = data . paramsdict [ 'reference_sequence' ]
index_files = [ ]
# # Check for existence of index files . Default to bwa unless you specify smalt
if "smalt" in data . _hackersonly [ "aligner" ] : # These are smalt index files . Only referenced here to ensure they exist
    index_files . extend ( [ ".sma" , ".smi" ] )
else :
    index_files . extend ( [ ".amb" , ".ann" , ".bwt" , ".pac" , ".sa" ] )
# # samtools specific index
index_files . extend ( [ ".fai" ] )
# # If reference sequence already exists then bail out of this func
if not force :
    if all ( [ os . path . isfile ( refseq_file + i ) for i in index_files ] ) :
        return
# if data . _ headers :
# print ( INDEX _ MSG . format ( data . _ hackersonly [ " aligner " ] ) )
if "smalt" in data . _hackersonly [ "aligner" ] : # # Create smalt index for mapping
# # smalt index [ - k < wordlen > ] [ - s < stepsiz > ] < index _ name > < reference _ file >
    cmd1 = [ ipyrad . bins . smalt , "index" , "-k" , str ( data . _hackersonly [ "smalt_index_wordlen" ] ) , refseq_file , refseq_file ]
else : # # bwa index < reference _ file >
    cmd1 = [ ipyrad . bins . bwa , "index" , refseq_file ]
# # call the command
LOGGER . info ( " " . join ( cmd1 ) )
proc1 = sps . Popen ( cmd1 , stderr = sps . STDOUT , stdout = sps . PIPE )
error1 = proc1 . communicate ( ) [ 0 ]
# # simple samtools index for grabbing ref seqs
cmd2 = [ ipyrad . bins . samtools , "faidx" , refseq_file ]
LOGGER . info ( " " . join ( cmd2 ) )
proc2 = sps . Popen ( cmd2 , stderr = sps . STDOUT , stdout = sps . PIPE )
error2 = proc2 . communicate ( ) [ 0 ]
# # error handling
if proc1 . returncode :
    raise IPyradWarningExit ( error1 )
if error2 :
    if "please use bgzip" in error2 :
        raise IPyradWarningExit ( NO_ZIP_BINS . format ( refseq_file ) )
    else :
        raise IPyradWarningExit ( error2 )


signature:
def mapreads ( data , sample , nthreads , force ) :
    """Attempt to map reads to reference sequence . This reads in the fasta files
 ( samples . files . edits ) , and maps each read to the reference . Unmapped reads
 are dropped right back in the de novo pipeline . Reads that map successfully
 are processed and pushed downstream and joined with the rest of the data
 post muscle _ align .
 Mapped reads end up in a sam file ."""

implementation:
LOGGER . info ( "Entering mapreads(): %s %s" , sample . name , nthreads )
# # This is the input derep file , for paired data we need to split the data ,
# # and so we will make sample . files . dereps = = [ derep1 , derep2 ] , but for
# # SE data we can simply use sample . files . derep = = [ derepfile ] .
derepfile = os . path . join ( data . dirs . edits , sample . name + "_derep.fastq" )
sample . files . dereps = [ derepfile ]
# # This is the final output files containing merged / concat derep ' d refmap ' d
# # reads that did not match to the reference . They will be back in
# # merge / concat ( - - nnnnn - - ) format ready to be input to vsearch , if needed .
mumapfile = sample . files . unmapped_reads
umap1file = os . path . join ( data . dirs . edits , sample . name + "-tmp-umap1.fastq" )
umap2file = os . path . join ( data . dirs . edits , sample . name + "-tmp-umap2.fastq" )
# # split the derepfile into the two handles we designate
if "pair" in data . paramsdict [ "datatype" ] :
    sample . files . split1 = os . path . join ( data . dirs . edits , sample . name + "-split1.fastq" )
    sample . files . split2 = os . path . join ( data . dirs . edits , sample . name + "-split2.fastq" )
    sample . files . dereps = [ sample . files . split1 , sample . files . split2 ]
    split_merged_reads ( sample . files . dereps , derepfile )
# # ( cmd1 ) smalt < task > [ TASK _ OPTIONS ] [ < index _ name > < file _ name _ A > [ < file _ name _ B > ] ]
# # - f sam : Output as sam format , tried : clip : to hard mask output
# # but it shreds the unmapped reads ( outputs empty fq )
# # - l [ pe , mp , pp ] : If paired end select the orientation of each read
# # - n # : Number of threads to use
# # - x : Perform a more exhaustive search
# # - y # : proportion matched to reference ( sequence similarity )
# # - o : output file
# # : Reference sequence
# # : Input file ( s ) , in a list . One for R1 and one for R2
# # - c # : proportion of the query read length that must be covered
# # ( cmd1 ) bwa mem [ OPTIONS ] < index _ name > < file _ name _ A > [ < file _ name _ B > ] > < output _ file >
# # - t # : Number of threads
# # - M : Mark split alignments as secondary .
# # ( cmd2 ) samtools view [ options ] < in . bam > | < in . sam > | < in . cram > [ region . . . ]
# # - b = write to . bam
# # - q = Only keep reads with mapq score > = 30 ( seems to be pretty standard )
# # - F = Select all reads that DON ' T have these flags .
# # 0x4 ( segment unmapped )
# # 0x100 ( Secondary alignment )
# # 0x800 ( supplementary alignment )
# # - U = Write out all reads that don ' t pass the - F filter
# # ( all unmapped reads go to this file ) .
# # TODO : Should eventually add ` - q 13 ` to filter low confidence mapping .
# # If you do this it will throw away some fraction of reads . Ideally you ' d
# # catch these and throw them in with the rest of the unmapped reads , but
# # I can ' t think of a straightforward way of doing that . There should be
# # a ` - Q ` flag to only keep reads below the threshold , but i realize that
# # would be of limited use besides for me .
# # ( cmd3 ) samtools sort [ options . . . ] [ in . bam ]
# # - T = Temporary file name , this is required by samtools , ignore it
# # Here we hack it to be samhandle . tmp cuz samtools cleans it up
# # - O = Output file format , in this case bam
# # - o = Output file name
if "smalt" in data . _hackersonly [ "aligner" ] : # # The output SAM data is written to file ( - o )
# # input is either ( derep ) or ( derep - split1 , derep - split2)
    cmd1 = [ ipyrad . bins . smalt , "map" , "-f" , "sam" , "-n" , str ( max ( 1 , nthreads ) ) , "-y" , str ( data . paramsdict [ 'clust_threshold' ] ) , "-o" , os . path . join ( data . dirs . refmapping , sample . name + ".sam" ) , "-x" , data . paramsdict [ 'reference_sequence' ] ] + sample . files . dereps
    cmd1_stdout = sps . PIPE
    cmd1_stderr = sps . STDOUT
else :
    cmd1 = [ ipyrad . bins . bwa , "mem" , "-t" , str ( max ( 1 , nthreads ) ) , "-M" , data . paramsdict [ 'reference_sequence' ] ] + sample . files . dereps
    # # Insert optional flags for bwa
    try :
        bwa_args = data . _hackersonly [ "bwa_args" ] . split ( )
        bwa_args . reverse ( )
        for arg in bwa_args :
            cmd1 . insert ( 2 , arg )
    except KeyError : # # Do nothing
        pass
    cmd1_stdout = open ( os . path . join ( data . dirs . refmapping , sample . name + ".sam" ) , 'w' )
    cmd1_stderr = None
# # Reads in the SAM file from cmd1 . It writes the unmapped data to file
# # and it pipes the mapped data to be used in cmd3
cmd2 = [ ipyrad . bins . samtools , "view" , "-b" , # # TODO : This introduces a bug with PE right now . Think about the case where
# # R1 has low qual mapping and R2 has high . You get different numbers
# # of reads in the unmapped tmp files . FML .
# " - q " , " 30 " ,
"-F" , "0x904" , "-U" , os . path . join ( data . dirs . refmapping , sample . name + "-unmapped.bam" ) , os . path . join ( data . dirs . refmapping , sample . name + ".sam" ) ]
# # this is gonna catch mapped bam output from cmd2 and write to file
cmd3 = [ ipyrad . bins . samtools , "sort" , "-T" , os . path . join ( data . dirs . refmapping , sample . name + ".sam.tmp" ) , "-O" , "bam" , "-o" , sample . files . mapped_reads ]
# # TODO : Unnecessary ?
# # this is gonna read the sorted BAM file and index it . only for pileup ?
cmd4 = [ ipyrad . bins . samtools , "index" , sample . files . mapped_reads ]
# # this is gonna read in the unmapped files , args are added below ,
# # and it will output fastq formatted unmapped reads for merging .
# # - v 45 sets the default qscore arbitrarily high
cmd5 = [ ipyrad . bins . samtools , "bam2fq" , "-v 45" , os . path . join ( data . dirs . refmapping , sample . name + "-unmapped.bam" ) ]
# # Insert additional arguments for paired data to the commands .
# # We assume Illumina paired end reads for the orientation
# # of mate pairs ( orientation : - - - > < - - - - ) .
if 'pair' in data . paramsdict [ "datatype" ] :
    if "smalt" in data . _hackersonly [ "aligner" ] : # # add paired flag ( - l pe ) to cmd1 right after ( smalt map . . . )
        cmd1 . insert ( 2 , "pe" )
        cmd1 . insert ( 2 , "-l" )
    else : # # No special PE flags for bwa
        pass
    # # add samtools filter for only keep if both pairs hit
    # # 0x1 - Read is paired
    # # 0x2 - Each read properly aligned
    cmd2 . insert ( 2 , "0x3" )
    cmd2 . insert ( 2 , "-f" )
    # # tell bam2fq that there are output files for each read pair
    cmd5 . insert ( 2 , umap1file )
    cmd5 . insert ( 2 , "-1" )
    cmd5 . insert ( 2 , umap2file )
    cmd5 . insert ( 2 , "-2" )
else :
    cmd5 . insert ( 2 , mumapfile )
    cmd5 . insert ( 2 , "-0" )
# # Running cmd1 creates ref _ mapping / sname . sam ,
LOGGER . debug ( " " . join ( cmd1 ) )
proc1 = sps . Popen ( cmd1 , stderr = cmd1_stderr , stdout = cmd1_stdout )
# # This is really long running job so we wrap it to ensure it dies .
try :
    error1 = proc1 . communicate ( ) [ 0 ]
except KeyboardInterrupt :
    proc1 . kill ( )
# # raise error if one occurred in smalt
if proc1 . returncode :
    raise IPyradWarningExit ( error1 )
# # Running cmd2 writes to ref _ mapping / sname . unmapped . bam , and
# # fills the pipe with mapped BAM data
LOGGER . debug ( " " . join ( cmd2 ) )
proc2 = sps . Popen ( cmd2 , stderr = sps . STDOUT , stdout = sps . PIPE )
# # Running cmd3 pulls mapped BAM from pipe and writes to
# # ref _ mapping / sname . mapped - sorted . bam .
# # Because proc2 pipes to proc3 we just communicate this to run both .
LOGGER . debug ( " " . join ( cmd3 ) )
proc3 = sps . Popen ( cmd3 , stderr = sps . STDOUT , stdout = sps . PIPE , stdin = proc2 . stdout )
error3 = proc3 . communicate ( ) [ 0 ]
if proc3 . returncode :
    raise IPyradWarningExit ( error3 )
proc2 . stdout . close ( )
# # Later we ' re gonna use samtools to grab out regions using ' view ' , and to
# # do that we need it to be indexed . Let ' s index it now .
LOGGER . debug ( " " . join ( cmd4 ) )
proc4 = sps . Popen ( cmd4 , stderr = sps . STDOUT , stdout = sps . PIPE )
error4 = proc4 . communicate ( ) [ 0 ]
if proc4 . returncode :
    raise IPyradWarningExit ( error4 )
# # Running cmd5 writes to either edits / sname - refmap _ derep . fastq for SE
# # or it makes edits / sname - tmp - umap { 12 } . fastq for paired data , which
# # will then need to be merged .
LOGGER . debug ( " " . join ( cmd5 ) )
proc5 = sps . Popen ( cmd5 , stderr = sps . STDOUT , stdout = sps . PIPE )
error5 = proc5 . communicate ( ) [ 0 ]
if proc5 . returncode :
    raise IPyradWarningExit ( error5 )
# # Finally , merge the unmapped reads , which is what cluster ( )
# # expects . If SE , just rename the outfile . In the end
# # < sample > - refmap _ derep . fq will be the final output
if 'pair' in data . paramsdict [ "datatype" ] :
    LOGGER . info ( "Merging unmapped reads {} {}" . format ( umap1file , umap2file ) )
    merge_pairs_after_refmapping ( data , [ ( umap1file , umap2file ) ] , mumapfile )


signature:
def fetch_cluster_se ( data , samfile , chrom , rstart , rend ) :
    """Builds a single end cluster from the refmapped data ."""

implementation:
# # If SE then we enforce the minimum overlap distance to avoid the
# # staircase syndrome of multiple reads overlapping just a little .
overlap_buffer = data . _hackersonly [ "min_SE_refmap_overlap" ]
# # the * _ buff variables here are because we have to play patty
# # cake here with the rstart / rend vals because we want pysam to
# # enforce the buffer for SE , but we want the reference sequence
# # start and end positions to print correctly for downstream .
rstart_buff = rstart + overlap_buffer
rend_buff = rend - overlap_buffer
# # Reads that map to only very short segements of the reference
# # sequence will return buffer end values that are before the
# # start values causing pysam to complain . Very short mappings .
if rstart_buff > rend_buff :
    tmp = rstart_buff
    rstart_buff = rend_buff
    rend_buff = tmp
# # Buffering can ' t make start and end equal or pysam returns nothing .
if rstart_buff == rend_buff :
    rend_buff += 1
# # store pairs
rdict = { }
clust = [ ]
iterreg = [ ]
iterreg = samfile . fetch ( chrom , rstart_buff , rend_buff )
# # use dict to match up read pairs
for read in iterreg :
    if read . qname not in rdict :
        rdict [ read . qname ] = read
# # sort dict keys so highest derep is first ( ' seed ' )
sfunc = lambda x : int ( x . split ( ";size=" ) [ 1 ] . split ( ";" ) [ 0 ] )
rkeys = sorted ( rdict . keys ( ) , key = sfunc , reverse = True )
# # get blocks from the seed for filtering , bail out if seed is not paired
try :
    read1 = rdict [ rkeys [ 0 ] ]
except ValueError :
    LOGGER . error ( "Found bad cluster, skipping - key:{} rdict:{}" . format ( rkeys [ 0 ] , rdict ) )
    return ""
# # the starting blocks for the seed
poss = read1 . get_reference_positions ( full_length = True )
seed_r1start = min ( poss )
seed_r1end = max ( poss )
# # store the seed - - - - -
if read1 . is_reverse :
    seq = revcomp ( read1 . seq )
else :
    seq = read1 . seq
# # store , could write orient but just + for now .
size = sfunc ( rkeys [ 0 ] )
clust . append ( ">{}:{}:{};size={};*\n{}" . format ( chrom , seed_r1start , seed_r1end , size , seq ) )
# # If there ' s only one hit in this region then rkeys will only have
# # one element and the call to ` rkeys [ 1 : ] ` will raise . Test for this .
if len ( rkeys ) > 1 : # # store the hits to the seed - - - - -
    for key in rkeys [ 1 : ] :
        skip = False
        try :
            read1 = rdict [ key ]
        except ValueError : # # enter values that will make this read get skipped
            read1 = rdict [ key ] [ 0 ]
            skip = True
        # # orient reads only if not skipping
        if not skip :
            poss = read1 . get_reference_positions ( full_length = True )
            minpos = min ( poss )
            maxpos = max ( poss )
            # # store the seq
            if read1 . is_reverse :
                seq = revcomp ( read1 . seq )
            else :
                seq = read1 . seq
            # # store , could write orient but just + for now .
            size = sfunc ( key )
            clust . append ( ">{}:{}:{};size={};+\n{}" . format ( chrom , minpos , maxpos , size , seq ) )
        else : # # seq is excluded , though , we could save it and return
        # # it as a separate cluster that will be aligned separately .
            pass
return clust


signature:
def fetch_cluster_pairs ( data , samfile , chrom , rstart , rend ) :
    """Builds a paired cluster from the refmapped data ."""

implementation:
# # store pairs
rdict = { }
clust = [ ]
# # grab the region and make tuples of info
iterreg = samfile . fetch ( chrom , rstart , rend )
# # use dict to match up read pairs
for read in iterreg :
    if read . qname not in rdict :
        rdict [ read . qname ] = [ read ]
    else :
        rdict [ read . qname ] . append ( read )
# # sort dict keys so highest derep is first ( ' seed ' )
sfunc = lambda x : int ( x . split ( ";size=" ) [ 1 ] . split ( ";" ) [ 0 ] )
rkeys = sorted ( rdict . keys ( ) , key = sfunc , reverse = True )
# # get blocks from the seed for filtering , bail out if seed is not paired
try :
    read1 , read2 = rdict [ rkeys [ 0 ] ]
except ValueError :
    return 0
# # the starting blocks for the seed
poss = read1 . get_reference_positions ( ) + read2 . get_reference_positions ( )
seed_r1start = min ( poss )
seed_r2end = max ( poss )
# # store the seed - - - - -
# # Simplify . R1 and R2 are always on opposite strands , but the
# # orientation is variable . We revcomp and order the reads to
# # preserve genomic order .
reads_overlap = False
if read1 . is_reverse :
    if read2 . aend > read1 . get_blocks ( ) [ 0 ] [ 0 ] :
        reads_overlap = True
        seq = read2 . seq + "nnnn" + revcomp ( read1 . seq )
    else :
        seq = read2 . seq + "nnnn" + read1 . seq
else :
    if read1 . aend > read2 . get_blocks ( ) [ 0 ] [ 0 ] :
        reads_overlap = True
        seq = read1 . seq + "nnnn" + revcomp ( read2 . seq )
    else :
        seq = read1 . seq + "nnnn" + read2 . seq
# # store , could write orient but just + for now .
size = sfunc ( rkeys [ 0 ] )
clust . append ( ">{}:{}:{};size={};*\n{}" . format ( chrom , seed_r1start , seed_r2end , size , seq ) )
# # If there ' s only one hit in this region then rkeys will only have
# # one element and the call to ` rkeys [ 1 : ] ` will raise . Test for this .
if len ( rkeys ) > 1 : # # store the hits to the seed - - - - -
    for key in rkeys [ 1 : ] :
        skip = False
        try :
            read1 , read2 = rdict [ key ]
        except ValueError : # # enter values that will make this read get skipped
            read1 = rdict [ key ] [ 0 ]
            read2 = read1
            skip = True
        # # orient reads and filter out ones that will not align well b / c
        # # they do not overlap enough with the seed
        poss = read1 . get_reference_positions ( ) + read2 . get_reference_positions ( )
        minpos = min ( poss )
        maxpos = max ( poss )
        # # skip if more than one hit location
        if read1 . has_tag ( "SA" ) or read2 . has_tag ( "SA" ) :
            skip = True
        # # store if read passes
        if ( abs ( minpos - seed_r1start ) < 50 ) and ( abs ( maxpos - seed_r2end ) < 50 ) and ( not skip ) : # # store the seq
            if read1 . is_reverse : # # do reads overlap
                if read2 . aend > read1 . get_blocks ( ) [ 0 ] [ 0 ] :
                    reads_overlap = True
                    seq = read2 . seq + "nnnn" + revcomp ( read1 . seq )
                else :
                    seq = read2 . seq + "nnnn" + read1 . seq
            else :
                if read1 . aend > read2 . get_blocks ( ) [ 0 ] [ 0 ] :
                    reads_overlap = True
                    seq = read1 . seq + "nnnn" + revcomp ( read2 . seq )
                else :
                    seq = read1 . seq + "nnnn" + read2 . seq
            # # store , could write orient but just + for now .
            size = sfunc ( key )
            clust . append ( ">{}:{}:{};size={};+\n{}" . format ( chrom , minpos , maxpos , size , seq ) )
        else : # # seq is excluded , though , we could save it and return
        # # it as a separate cluster that will be aligned separately .
            pass
# # merge the pairs prior to returning them
# # Remember , we already tested for quality scores , so
# # merge _ after _ pysam will generate arbitrarily high scores
# # It would be nice to do something here like test if
# # the average insert length + 2 stdv is > 2 * read len
# # so you can switch off merging for mostly non - overlapping data
if reads_overlap :
    if data . _hackersonly [ "refmap_merge_PE" ] :
        clust = merge_after_pysam ( data , clust )
        # clust = merge _ pair _ pipes ( data , clust )
return clust


signature:
def ref_build_and_muscle_chunk ( data , sample ) :
    """1 . Run bedtools to get all overlapping regions
 2 . Parse out reads from regions using pysam and dump into chunk files .
 We measure it out to create 10 chunk files per sample .
 3 . If we really wanted to speed this up , though it is pretty fast already ,
 we could parallelize it since we can easily break the regions into
 a list of chunks ."""

implementation:
# # get regions using bedtools
regions = bedtools_merge ( data , sample ) . strip ( ) . split ( "\n" )
nregions = len ( regions )
chunksize = ( nregions / 10 ) + ( nregions % 10 )
LOGGER . debug ( "nregions {} chunksize {}" . format ( nregions , chunksize ) )
# # create an output file to write clusters to
idx = 0
tmpfile = os . path . join ( data . tmpdir , sample . name + "_chunk_{}.ali" )
# # remove old files if they exist to avoid append errors
for i in range ( 11 ) :
    if os . path . exists ( tmpfile . format ( i ) ) :
        os . remove ( tmpfile . format ( i ) )
fopen = open
# # If reference + denovo we drop the reads back into clust . gz
# # and let the muscle _ chunker do it ' s thing back in cluster _ within
if data . paramsdict [ "assembly_method" ] == "denovo+reference" :
    tmpfile = os . path . join ( data . dirs . clusts , sample . name + ".clust.gz" )
    fopen = gzip . open
# # build clusters for aligning with muscle from the sorted bam file
samfile = pysam . AlignmentFile ( sample . files . mapped_reads , 'rb' )
# " . / tortas _ refmapping / PZ70 - mapped - sorted . bam " , " rb " )
# # fill clusts list and dump periodically
clusts = [ ]
nclusts = 0
for region in regions :
    chrom , pos1 , pos2 = region . split ( )
    try : # # fetches pairs quickly but then goes slow to merge them .
        if "pair" in data . paramsdict [ "datatype" ] :
            clust = fetch_cluster_pairs ( data , samfile , chrom , int ( pos1 ) , int ( pos2 ) )
        # # fetch but no need to merge
        else :
            clust = fetch_cluster_se ( data , samfile , chrom , int ( pos1 ) , int ( pos2 ) )
    except IndexError as inst :
        LOGGER . error ( "Bad region chrom:start-end {}:{}-{}" . format ( chrom , pos1 , pos2 ) )
        continue
    if clust :
        clusts . append ( "\n" . join ( clust ) )
        nclusts += 1
        if nclusts == chunksize : # # write to file
            tmphandle = tmpfile . format ( idx )
            with fopen ( tmphandle , 'a' ) as tmp : # LOGGER . debug ( " Writing tmpfile - { } " . format ( tmpfile . format ( idx ) ) )
            # if data . paramsdict [ " assembly _ method " ] = = " denovo + reference " :
            # # # This is dumb , but for this method you need to prepend the
            # # # separator to maintain proper formatting of clust . gz
                tmp . write ( "\n//\n//\n" . join ( clusts ) + "\n//\n//\n" )
            idx += 1
            nclusts = 0
            clusts = [ ]
if clusts : # # write remaining to file
    with fopen ( tmpfile . format ( idx ) , 'a' ) as tmp : # tmp . write ( " \ n / / \ n / / \ n " + ( " \ n / / \ n / / \ n " . join ( clusts ) ) )
        tmp . write ( "\n//\n//\n" . join ( clusts ) + "\n//\n//\n" )
    clusts = [ ]
if not data . paramsdict [ "assembly_method" ] == "denovo+reference" :
    chunkfiles = glob . glob ( os . path . join ( data . tmpdir , sample . name + "_chunk_*.ali" ) )
    LOGGER . info ( "created chunks %s" , chunkfiles )
# # cleanup
samfile . close ( )


signature:
def ref_muscle_chunker ( data , sample ) :
    """Run bedtools to get all overlapping regions . Pass this list into the func
 ' get _ overlapping _ reads ' which will write fastq chunks to the clust . gz file .
 1 ) Run bedtools merge to get a list of all contiguous blocks of bases
 in the reference seqeunce where one or more of our reads overlap .
 The output will look like this :
 1 45230754 45230783
 1 74956568 74956596
 1 116202035 116202060"""

implementation:
LOGGER . info ( 'entering ref_muscle_chunker' )
# # Get regions , which will be a giant list of 5 - tuples , of which we ' re
# # only really interested in the first three : ( chrom , start , end ) position .
regions = bedtools_merge ( data , sample )
if len ( regions ) > 0 : # # this calls bam _ region _ to _ fasta a billion times
    get_overlapping_reads ( data , sample , regions )
else :
    msg = "No reads mapped to reference sequence - {}" . format ( sample . name )
    LOGGER . warn ( msg )


signature:
def get_overlapping_reads ( data , sample , regions ) :
    """For SE data , this pulls mapped reads out of sorted mapped bam files and
 appends them to the clust . gz file so they fall into downstream
 ( muscle alignment ) analysis .
 For PE data , this pulls mapped reads out of sorted mapped bam files , splits
 R1s from R2s and writes them to separate files . Once all reads are written ,
 it calls merge _ reads ( vsearch ) to find merged and non - merged reads . These
 are then put into clust . gz with either an nnnn separator or as merged .
 The main func being called here is ' bam _ region _ to _ fasta ' , which calls
 samtools to pull out the mapped reads .
 1 ) Coming into this function we have sample . files . mapped _ reads
 as a sorted bam file , and a passed in list of regions to evaluate .
 2 ) Get all reads overlapping with each individual region .
 3 ) Pipe to vsearch for clustering .
 4 ) Append to the clust . gz file ."""

implementation:
# # storage and counter
locus_list = [ ]
reads_merged = 0
# # Set the write mode for opening clusters file .
# # 1 ) if " reference " then only keep refmapped , so use ' wb ' to overwrite
# # 2 ) if " denovo + reference " then ' ab ' adds to end of denovo clust file
write_flag = 'wb'
if data . paramsdict [ "assembly_method" ] == "denovo+reference" :
    write_flag = 'ab'
# # file handle for writing clusters
sample . files . clusters = os . path . join ( data . dirs . clusts , sample . name + ".clust.gz" )
outfile = gzip . open ( sample . files . clusters , write_flag )
# # write a separator if appending to clust . gz
if data . paramsdict [ "assembly_method" ] == "denovo+reference" :
    outfile . write ( "\n//\n//\n" )
# # Make a process to pass in to bam _ region _ to _ fasta so we can just reuse
# # it rather than recreating a bunch of subprocesses . Saves hella time .
proc1 = sps . Popen ( "sh" , stdin = sps . PIPE , stdout = sps . PIPE , universal_newlines = True )
# Wrap this in a try so we can easily locate errors
try : # # For each identified region , build the pileup and write out the fasta
    for line in regions . strip ( ) . split ( "\n" ) : # Blank lines returned from bedtools screw things up . Filter them .
        if line == "" :
            continue
        # # get elements from bedtools region
        chrom , region_start , region_end = line . strip ( ) . split ( ) [ 0 : 3 ]
        # # bam _ region _ to _ fasta returns a chunk of fasta sequence
        args = [ data , sample , proc1 , chrom , region_start , region_end ]
        clust = bam_region_to_fasta ( * args )
        # # If bam _ region _ to _ fasta fails for some reason it ' ll return [ ] ,
        # # in which case skip the rest of this . Normally happens if reads
        # # map successfully , but too far apart .
        if not clust :
            continue
        # # Store locus in a list
        # LOGGER . info ( " clust from bam - region - to - fasta \ n % s " , clust )
        locus_list . append ( clust )
        # # write chunk of 1000 loci and clear list to minimize memory
        if not len ( locus_list ) % 1000 :
            outfile . write ( "\n//\n//\n" . join ( locus_list ) + "\n//\n//\n" )
            locus_list = [ ]
    # # write remaining
    if any ( locus_list ) :
        outfile . write ( "\n//\n//\n" . join ( locus_list ) )
    else : # # If it ' s empty , strip off the last \ n / / \ n / / \ n from the outfile .
        pass
    # # close handle
    outfile . close ( )
except Exception as inst :
    LOGGER . error ( "Exception inside get_overlapping_reads - {}" . format ( inst ) )
    raise
finally :
    if "pair" in data . paramsdict [ "datatype" ] :
        LOGGER . info ( "Total merged reads for {} - {}" . format ( sample . name , reads_merged ) )
        sample . stats . reads_merged = reads_merged


signature:
def split_merged_reads ( outhandles , input_derep ) :
    """Takes merged / concat derep file from vsearch derep and split it back into
 separate R1 and R2 parts .
 - sample _ fastq : a list of the two file paths to write out to .
 - input _ reads : the path to the input merged reads"""

implementation:
handle1 , handle2 = outhandles
splitderep1 = open ( handle1 , 'w' )
splitderep2 = open ( handle2 , 'w' )
with open ( input_derep , 'r' ) as infile : # # Read in the infile two lines at a time : ( seqname , sequence )
    duo = itertools . izip ( * [ iter ( infile ) ] * 2 )
    # # lists for storing results until ready to write
    split1s = [ ]
    split2s = [ ]
    # # iterate over input splitting , saving , and writing .
    idx = 0
    while 1 :
        try :
            itera = duo . next ( )
        except StopIteration :
            break
        # # split the duo into separate parts and inc counter
        part1 , part2 = itera [ 1 ] . split ( "nnnn" )
        idx += 1
        # # R1 needs a newline , but R2 inherits it from the original file
        # # store parts in lists until ready to write
        split1s . append ( "{}{}\n" . format ( itera [ 0 ] , part1 ) )
        split2s . append ( "{}{}" . format ( itera [ 0 ] , part2 ) )
        # # if large enough then write to file
        if not idx % 10000 :
            splitderep1 . write ( "" . join ( split1s ) )
            splitderep2 . write ( "" . join ( split2s ) )
            split1s = [ ]
            split2s = [ ]
# # write final chunk if there is any
if any ( split1s ) :
    splitderep1 . write ( "" . join ( split1s ) )
    splitderep2 . write ( "" . join ( split2s ) )
# # close handles
splitderep1 . close ( )
splitderep2 . close ( )


signature:
def check_insert_size ( data , sample ) :
    """check mean insert size for this sample and update
 hackersonly . max _ inner _ mate _ distance if need be . This value controls how
 far apart mate pairs can be to still be considered for bedtools merging
 downstream ."""

implementation:
# # pipe stats output to grep
cmd1 = [ ipyrad . bins . samtools , "stats" , sample . files . mapped_reads ]
cmd2 = [ "grep" , "SN" ]
proc1 = sps . Popen ( cmd1 , stderr = sps . STDOUT , stdout = sps . PIPE )
proc2 = sps . Popen ( cmd2 , stderr = sps . STDOUT , stdout = sps . PIPE , stdin = proc1 . stdout )
# # get piped result
res = proc2 . communicate ( ) [ 0 ]
# # raise exception on failure and do cleanup
if proc2 . returncode :
    raise IPyradWarningExit ( "error in %s: %s" , cmd2 , res )
# # starting vals
avg_insert = 0
stdv_insert = 0
avg_len = 0
# # iterate over results
for line in res . split ( "\n" ) :
    if "insert size average" in line :
        avg_insert = float ( line . split ( ":" ) [ - 1 ] . strip ( ) )
    elif "insert size standard deviation" in line : # # hack to fix sim data when stdv is 0.0 . Shouldn ' t
    # # impact real data bcz stdv gets rounded up below
        stdv_insert = float ( line . split ( ":" ) [ - 1 ] . strip ( ) ) + 0.1
    elif "average length" in line :
        avg_len = float ( line . split ( ":" ) [ - 1 ] . strip ( ) )
LOGGER . debug ( "avg {} stdv {} avg_len {}" . format ( avg_insert , stdv_insert , avg_len ) )
# # If all values return successfully set the max inner mate distance .
# # This is tricky . avg _ insert is the average length of R1 + R2 + inner mate
# # distance . avg _ len is the average length of a read . If there are lots
# # of reads that overlap then avg _ insert will be close to but bigger than
# # avg _ len . We are looking for the right value for ` bedtools merge - d `
# # which wants to know the max distance between reads .
if all ( [ avg_insert , stdv_insert , avg_len ] ) : # # If 2 * the average length of a read is less than the average
# # insert size then most reads DO NOT overlap
    if stdv_insert < 5 :
        stdv_insert = 5.
    if ( 2 * avg_len ) < avg_insert :
        hack = avg_insert + ( 3 * np . math . ceil ( stdv_insert ) ) - ( 2 * avg_len )
    # # If it is > than the average insert size then most reads DO
    # # overlap , so we have to calculate inner mate distance a little
    # # differently .
    else :
        hack = ( avg_insert - avg_len ) + ( 3 * np . math . ceil ( stdv_insert ) )
    # # set the hackerdict value
    LOGGER . info ( "stdv: hacked insert size is %s" , hack )
    data . _hackersonly [ "max_inner_mate_distance" ] = int ( np . math . ceil ( hack ) )
else : # # If something fsck then set a relatively conservative distance
    data . _hackersonly [ "max_inner_mate_distance" ] = 300
    LOGGER . debug ( "inner mate distance for {} - {}" . format ( sample . name , data . _hackersonly [ "max_inner_mate_distance" ] ) )


signature:
def bedtools_merge ( data , sample ) :
    """Get all contiguous genomic regions with one or more overlapping
 reads . This is the shell command we ' ll eventually run
 bedtools bamtobed - i 1A _ 0 . sorted . bam | bedtools merge [ - d 100]
 - i < input _ bam > : specifies the input file to bed ' ize
 - d < int > : For PE set max distance between reads"""

implementation:
LOGGER . info ( "Entering bedtools_merge: %s" , sample . name )
mappedreads = os . path . join ( data . dirs . refmapping , sample . name + "-mapped-sorted.bam" )
# # command to call ` bedtools bamtobed ` , and pipe output to stdout
# # Usage : bedtools bamtobed [ OPTIONS ] - i < bam >
# # Usage : bedtools merge [ OPTIONS ] - i < bam >
cmd1 = [ ipyrad . bins . bedtools , "bamtobed" , "-i" , mappedreads ]
cmd2 = [ ipyrad . bins . bedtools , "merge" , "-i" , "-" ]
# # If PE the - d flag to tell bedtools how far apart to allow mate pairs .
# # If SE the - d flag is negative , specifying that SE reads need to
# # overlap by at least a specific number of bp . This prevents the
# # stairstep syndrome when a + and - read are both extending from
# # the same cutsite . Passing a negative number to ` merge - d ` gets this done .
if 'pair' in data . paramsdict [ "datatype" ] :
    check_insert_size ( data , sample )
    # cmd2 . insert ( 2 , str ( data . _ hackersonly [ " max _ inner _ mate _ distance " ] ) )
    cmd2 . insert ( 2 , str ( data . _hackersonly [ "max_inner_mate_distance" ] ) )
    cmd2 . insert ( 2 , "-d" )
else :
    cmd2 . insert ( 2 , str ( - 1 * data . _hackersonly [ "min_SE_refmap_overlap" ] ) )
    cmd2 . insert ( 2 , "-d" )
# # pipe output from bamtobed into merge
LOGGER . info ( "stdv: bedtools merge cmds: %s %s" , cmd1 , cmd2 )
proc1 = sps . Popen ( cmd1 , stderr = sps . STDOUT , stdout = sps . PIPE )
proc2 = sps . Popen ( cmd2 , stderr = sps . STDOUT , stdout = sps . PIPE , stdin = proc1 . stdout )
result = proc2 . communicate ( ) [ 0 ]
proc1 . stdout . close ( )
# # check for errors and do cleanup
if proc2 . returncode :
    raise IPyradWarningExit ( "error in %s: %s" , cmd2 , result )
# # Write the bedfile out , because it ' s useful sometimes .
if os . path . exists ( ipyrad . __debugflag__ ) :
    with open ( os . path . join ( data . dirs . refmapping , sample . name + ".bed" ) , 'w' ) as outfile :
        outfile . write ( result )
# # Report the number of regions we ' re returning
nregions = len ( result . strip ( ) . split ( "\n" ) )
LOGGER . info ( "bedtools_merge: Got # regions: %s" , nregions )
return result


signature:
def trim_reference_sequence ( fasta ) :
    """If doing PE and R1 / R2 don ' t overlap then the reference sequence
 will be quite long and will cause indel hell during the
 alignment stage . Here trim the reference sequence to the length
 of the merged reads . Input is a list of alternating locus labels
 and sequence data . The first locus label is the reference
 sequence label and the first seq is the reference seq . Returns
 the same list except with the reference sequence trimmed to
 the length of the rad tags"""

implementation:
LOGGER . debug ( "pre - {}" . format ( fasta [ 0 ] ) )
# # If the reads are merged then the reference sequence should be the
# # same length as the merged pair . If unmerged then we have to fix it .
if "nnnn" in fasta [ 1 ] :
    r1_len = len ( fasta [ 1 ] . split ( "\n" ) [ 1 ] . split ( "nnnn" ) [ 0 ] )
    r2_len = len ( fasta [ 1 ] . split ( "\n" ) [ 1 ] . split ( "nnnn" ) [ 1 ] )
    new_seq = fasta [ 0 ] . split ( "\n" ) [ 1 ] [ : r1_len ] + ( "nnnn" ) + revcomp ( fasta [ 0 ] . split ( "\n" ) [ 1 ] [ - r2_len : ] )
    fasta [ 0 ] = fasta [ 0 ] . split ( "\n" ) [ 0 ] + "\n" + new_seq
LOGGER . debug ( "post - {}" . format ( fasta [ 0 ] ) )
return fasta


signature:
def bam_region_to_fasta ( data , sample , proc1 , chrom , region_start , region_end ) :
    """Take the chromosome position , and start and end bases and return sequences
 of all reads that overlap these sites . This is the command we ' re building :
 samtools view - b 1A _ sorted . bam 1:116202035-116202060 | samtools bam2fq < options > -
 - b : output bam format
 -0 : For SE , output all reads to this file
 -1 / - 2 : For PE , output first and second reads to different files
 - : Tell samtools to read in from the pipe
 Write out the sam output and parse it to return as fasta for clust . gz file .
 We also grab the reference sequence with a @ REF header to aid in alignment
 for single - end data . This will be removed post - alignment ."""

implementation:
# # output bam file handle for storing genome regions
bamf = sample . files . mapped_reads
if not os . path . exists ( bamf ) :
    raise IPyradWarningExit ( "  file not found - %s" , bamf )
# chrom = re . escape ( repr ( chrom ) ) [ 1 : - 1 ] . replace ( ' \ \ \ \ ' , ' \ \ ' )
# LOGGER . info ( " before : % s " , chrom )
chrom . replace ( "|" , r"\|" )
chrom . replace ( "(" , r"\(" )
chrom . replace ( ")" , r"\)" )
# LOGGER . info ( " after : % s " , chrom )
# # What we want to do is have the num - chrom dict as an arg , then build this
# # string as three ints [ chrom - int , pos - start , pos - end ]
# cint = cdict [ chrom ]
# cpstring = " _ _ { } _ { } _ { } _ _ " . format ( cint , int ( region _ start ) + 1 , region _ end )
# # a string argument as input to commands , indexed at either 0 or 1,
# # and with pipe characters removed from chromo names
# # rstring _ id1 is for fetching the reference sequence bcz faidx is
# # 1 indexed
rstring_id1 = "{}:{}-{}" . format ( chrom , str ( int ( region_start ) + 1 ) , region_end )
# # rstring _ id0 is just for printing out the reference CHROM / POS
# # in the read name
rstring_id0 = "{}:{}-{}" . format ( chrom , region_start , region_end )
# # If SE then we enforce the minimum overlap distance to avoid the
# # staircase syndrome of multiple reads overlapping just a little .
overlap_buffer = 0
if not "pair" in data . paramsdict [ "datatype" ] :
    overlap_buffer = data . _hackersonly [ "min_SE_refmap_overlap" ]
# # rstring _ id0 _ buffered is for samtools view . We have to play patty
# # cake here with the two rstring _ id0s because we want ` view ` to
# # enforce the buffer for SE , but we want the reference sequence
# # start and end positions to print correctly for downstream .
rstring_id0_buffered = "{}:{}-{}" . format ( chrom , int ( region_start ) + overlap_buffer , int ( region_end ) - overlap_buffer )
# # The " samtools faidx " command will grab this region from reference
# # which we ' ll paste in at the top of each stack to aid alignment .
cmd1 = [ ipyrad . bins . samtools , "faidx" , data . paramsdict [ "reference_sequence" ] , rstring_id1 , " ; echo __done__" ]
# # Call the command , I found that it doesn ' t work with shell = False if
# # the refstring is ' MT ' : 100-200 ' , but it works if it is MT : 100-200.
LOGGER . info ( "Grabbing bam_region_to_fasta:\n {}" . format ( cmd1 ) )
# proc1 = sps . Popen ( cmd1 , stderr = sps . STDOUT , stdout = sps . PIPE )
# ref = proc1 . communicate ( ) [ 0]
# if proc1 . returncode :
# raise IPyradWarningExit ( " error in % s : % s " , cmd1 , ref )
# # push the samtools faidx command to our subprocess , then accumulate
# # the results from stdout
print ( " " . join ( cmd1 ) , file = proc1 . stdin )
ref = ""
for line in iter ( proc1 . stdout . readline , "//\n" ) :
    if "__done__" in line :
        break
    ref += line
# # initialize the fasta list .
fasta = [ ]
# # parse sam to fasta . Save ref location to name .
# # Set size = an improbably large value so the REF sequence
# # sorts to the top for muscle aligning .
try :
    name , seq = ref . strip ( ) . split ( "\n" , 1 )
    seq = "" . join ( seq . split ( "\n" ) )
    fasta = [ "{}_REF;size={};+\n{}" . format ( name , 1000000 , seq ) ]
except ValueError as inst :
    LOGGER . error ( "ref failed to parse - {}" . format ( ref ) )
    LOGGER . error ( " " . join ( cmd1 ) )
# # if PE then you have to merge the reads here
if "pair" in data . paramsdict [ "datatype" ] : # # PE R1 can either be on the forward or the reverse strand .
# # Samtools view always outputs reads with respect to the
# # forward strand . This means that reads with R1 on reverse
# # end up with the R1 and R2 reference sequences swapped
# # in the clust . gz file . There is a way to fix it but it ' s
# # very annoying and i ' m not sure if it ' s worth it . . .
# # Drop the reference sequence for now . . .
# # If you ever fix this be sure to remove the reference sequence
# # from each cluster post alignment in cluster _ within / align _ and _ parse ( )
    fasta = [ ]
    # # Create temporary files for R1 , R2 and merged , which we will pass to
    # # the function merge _ pairs ( ) which calls vsearch to test merging .
    # # If you are on linux then creating the temp files in / dev / shm
    # # should improve performance
    if os . path . exists ( "/dev/shm" ) :
        prefix = os . path . join ( "/dev/shm" , "{}-{}" . format ( sample . name , rstring_id0 ) )
    else :
        prefix = os . path . join ( data . dirs . refmapping , "{}-{}" . format ( sample . name , rstring_id0 ) )
    read1 = "{}-R1" . format ( prefix )
    read2 = "{}-R2" . format ( prefix )
    merged = "{}-merged" . format ( prefix )
    # # Grab all the reads that map to this genomic location and dump
    # # fastq to R1 and R2 files .
    # # ` - v 45 ` sets the default qscore to something high
    cmd1 = " " . join ( [ ipyrad . bins . samtools , "view" , "-b" , bamf , rstring_id0 ] )
    cmd2 = " " . join ( [ ipyrad . bins . samtools , "bam2fq" , "-v" , "45" , "-1" , read1 , "-2" , read2 , "-" , "; echo __done__" ] )
    cmd = " | " . join ( [ cmd1 , cmd2 ] )
    print ( cmd , file = proc1 . stdin )
    for line in iter ( proc1 . stdout . readline , "//\n" ) :
        if "__done__" in line :
            break
    # # run commands , pipe 1 - > 2 , then cleanup
    # # proc1 = sps . Popen ( cmd1 , stderr = sps . STDOUT , stdout = sps . PIPE )
    # # proc2 = sps . Popen ( cmd2 , stderr = sps . STDOUT , stdout = sps . PIPE , stdin = proc1 . stdout )
    # # res = proc2 . communicate ( ) [ 0]
    # # if proc2 . returncode :
    # # raise IPyradWarningExit ( " error { } : { } " . format ( cmd2 , res ) )
    # # proc1 . stdout . close ( )
    # # merge the pairs . 0 means don ' t revcomp bcz samtools already
    # # did it for us . 1 means " actually merge " .
    try : # # return number of merged reads , writes merged data to ' merged '
    # # we don ' t yet do anything with the returned number of merged
        _ = merge_pairs ( data , [ ( read1 , read2 ) ] , merged , 0 , 1 )
        with open ( merged , 'r' ) as infile :
            quatro = itertools . izip ( * [ iter ( infile ) ] * 4 )
            while 1 :
                try :
                    bits = quatro . next ( )
                except StopIteration :
                    break
                # # TODO : figure out a real way to get orientation for PE
                orient = "+"
                fullfast = ">{a};{b};{c};{d}\n{e}" . format ( a = bits [ 0 ] . split ( ";" ) [ 0 ] , b = rstring_id1 , c = bits [ 0 ] . split ( ";" ) [ 1 ] , d = orient , e = bits [ 1 ] . strip ( ) )
                # , e = bits [ 9 ] )
                fasta . append ( fullfast )
            # # TODO : If you ever figure out a good way to get the reference
            # # sequence included w / PE then this commented call is useful
            # # for trimming the reference sequence to be the right length .
            # # If doing PE and R1 / R2 don ' t overlap then the reference sequence
            # # will be quite long and will cause indel hell during the
            # # alignment stage . Here trim the reference sequence to the length
            # # of the merged reads .
            # # This is commented out because we aren ' t currently including the
            # # ref seq for PE alignment .
            # fasta = trim _ reference _ sequence ( fasta )
    except ( OSError , ValueError , IPyradError ) as inst : # # ValueError raised inside merge _ pairs ( ) if it can ' t open one
    # # or both of the files . Write this out , but ignore for now .
    # # Failed merging , probably unequal number of reads in R1 and R2
    # # IPyradError raised if merge _ pairs can ' t read either R1 or R2
    # # file .
    # # Skip this locus ?
        LOGGER . debug ( "Failed to merge reads, continuing; %s" , inst )
        LOGGER . error ( "cmd - {}" . format ( cmd ) )
        return ""
    finally : # # Only clean up the files if they exist otherwise it ' ll raise .
        if os . path . exists ( merged ) :
            os . remove ( merged )
        if os . path . exists ( read1 ) :
            os . remove ( read1 )
        if os . path . exists ( read2 ) :
            os . remove ( read2 )
else :
    try : # # SE if faster than PE cuz it skips writing intermedidate files
    # # rstring _ id0 _ buffered is going to enforce the required
    # # min _ SE _ refmap _ overlap on either end of this region .
        cmd2 = [ ipyrad . bins . samtools , "view" , bamf , rstring_id0_buffered ]
        proc2 = sps . Popen ( cmd2 , stderr = sps . STDOUT , stdout = sps . PIPE )
        # # run and check outputs
        res = proc2 . communicate ( ) [ 0 ]
        if proc2 . returncode :
            raise IPyradWarningExit ( "{} {}" . format ( cmd2 , res ) )
        # # if the region string is malformated you ' ll get back a warning
        # # from samtools
        if "[main_samview]" in res :
            raise IPyradError ( "Bad reference region {}" . format ( rstring_id0_buffered ) )
        # # do not join seqs that
        for line in res . strip ( ) . split ( "\n" ) :
            bits = line . split ( "\t" )
            # # Read in the 2nd field ( FLAGS ) , convert to binary
            # # and test if the 7th bit is set which indicates revcomp
            orient = "+"
            if int ( '{0:012b}' . format ( int ( bits [ 1 ] ) ) [ 7 ] ) :
                orient = "-"
                # # Don ' t actually revcomp the sequence because samtools
                # # writes out reference sequence on the forward strand
                # # as well as reverse strand hits from the bam file .
                # bits [ 9 ] = revcomp ( bits [ 9 ] )
            # # Rip insert the mapping position between the seq label and
            # # the vsearch derep size .
            fullfast = ">{a};{b};{c};{d}\n{e}" . format ( a = bits [ 0 ] . split ( ";" ) [ 0 ] , b = rstring_id0 , c = bits [ 0 ] . split ( ";" ) [ 1 ] , d = orient , e = bits [ 9 ] )
            fasta . append ( fullfast )
    except IPyradError as inst : # # If the mapped fragment is too short then the you ' ll get
    # # regions that look like this : scaffold262:299039-299036
    # # Just carry on , it ' s not a big deal .
        LOGGER . debug ( "Got a bad region string: {}" . format ( inst ) )
        return ""
    except ( OSError , ValueError , Exception ) as inst : # # Once in a blue moon something fsck and it breaks the
    # # assembly . No reason to give up if . 001 % of reads fail
    # # so just skip this locus .
        LOGGER . error ( "Failed get reads at a locus, continuing; %s" , inst )
        LOGGER . error ( "cmd - {}" . format ( cmd2 ) )
        return ""
return "\n" . join ( fasta )


signature:
def refmap_stats ( data , sample ) :
    """Get the number of mapped and unmapped reads for a sample
 and update sample . stats"""

implementation:
# # shorter names
mapf = os . path . join ( data . dirs . refmapping , sample . name + "-mapped-sorted.bam" )
umapf = os . path . join ( data . dirs . refmapping , sample . name + "-unmapped.bam" )
# # get from unmapped
cmd1 = [ ipyrad . bins . samtools , "flagstat" , umapf ]
proc1 = sps . Popen ( cmd1 , stderr = sps . STDOUT , stdout = sps . PIPE )
result1 = proc1 . communicate ( ) [ 0 ]
# # get from mapped
cmd2 = [ ipyrad . bins . samtools , "flagstat" , mapf ]
proc2 = sps . Popen ( cmd2 , stderr = sps . STDOUT , stdout = sps . PIPE )
result2 = proc2 . communicate ( ) [ 0 ]
# # store results
# # If PE , samtools reports the _ actual _ number of reads mapped , both
# # R1 and R2 , so here if PE divide the results by 2 to stay consistent
# # with how we ' ve been reporting R1 and R2 as one " read pair "
if "pair" in data . paramsdict [ "datatype" ] :
    sample . stats [ "refseq_unmapped_reads" ] = int ( result1 . split ( ) [ 0 ] ) / 2
    sample . stats [ "refseq_mapped_reads" ] = int ( result2 . split ( ) [ 0 ] ) / 2
else :
    sample . stats [ "refseq_unmapped_reads" ] = int ( result1 . split ( ) [ 0 ] )
    sample . stats [ "refseq_mapped_reads" ] = int ( result2 . split ( ) [ 0 ] )
sample_cleanup ( data , sample )


signature:
def refmap_init ( data , sample , force ) :
    """create some file handles for refmapping"""

implementation:
# # make some persistent file handles for the refmap reads files
sample . files . unmapped_reads = os . path . join ( data . dirs . edits , "{}-refmap_derep.fastq" . format ( sample . name ) )
sample . files . mapped_reads = os . path . join ( data . dirs . refmapping , "{}-mapped-sorted.bam" . format ( sample . name ) )


signature:
def parse_command_line ( ) :
    """Parse CLI args . Only three options now ."""

implementation:
# # create the parser
parser = argparse . ArgumentParser ( formatter_class = argparse . RawDescriptionHelpFormatter , epilog = """
  * Example command-line usage ----------------------------------------------

  * Read in sequence/SNP data file, provide linkage, output name, ambig option.
     tetrad -s data.snps.phy -n test             ## input phylip and give name
     tetrad -s data.snps.phy -l data.snps.map    ## sample one SNP per locus
     tetrad -s data.snps.phy -n noambigs -r 0    ## do not use hetero sites

  * Load saved/checkpointed analysis from '.tet.json' file, or force restart.
     tetrad -j test.tet.json -b 100         ## continue 'test' until 100 boots
     tetrad -j test.tet.json -b 100 -f      ## force restart of 'test'

  * Sampling modes: 'equal' uses guide tree to sample quartets more efficiently
     tetrad -s data.snps.phy -m all                       ## sample all quartets
     tetrad -s data.snps.phy -m random -q 1e6 -x 123      ## sample 1M randomly
     tetrad -s data.snps.phy -m equal -q 1e6 -t guide.tre ## sample 1M across tree

  * Connect to N cores on a computer (default without -c arg is to use all avail.)
     tetrad -s data.snps.phy -c 20

  * Start an MPI cluster to connect to nodes across multiple available hosts.
     tetrad -s data.snps.phy --MPI

  * Connect to a manually started ipcluster instance with default or named profile
     tetrad -s data.snps.phy --ipcluster        ## connects to default profile
     tetrad -s data.snps.phy --ipcluster pname  ## connects to profile='pname'

  * Further documentation: http://ipyrad.readthedocs.io/analysis.html
    """ )
# # get version from ipyrad
ipyversion = str ( pkg_resources . get_distribution ( 'ipyrad' ) )
parser . add_argument ( '-v' , '--version' , action = 'version' , version = "tetrad " + ipyversion . split ( ) [ 1 ] )
parser . add_argument ( '-f' , "--force" , action = 'store_true' , help = "force overwrite of existing data" )
parser . add_argument ( '-s' , metavar = "seq" , dest = "seq" , type = str , default = None , help = "path to input phylip file (only SNPs)" )
parser . add_argument ( '-j' , metavar = 'json' , dest = "json" , type = str , default = None , help = "load checkpointed/saved analysis from JSON file." )
parser . add_argument ( '-m' , metavar = "method" , dest = "method" , type = str , default = "all" , help = "method for sampling quartets (all, random, or equal)" )
parser . add_argument ( '-q' , metavar = "nquartets" , dest = "nquartets" , type = int , default = 0 , help = "number of quartets to sample (if not -m all)" )
parser . add_argument ( '-b' , metavar = "boots" , dest = "boots" , type = int , default = 0 , help = "number of non-parametric bootstrap replicates" )
parser . add_argument ( '-l' , metavar = "map_file" , dest = "map" , type = str , default = None , help = "map file of snp linkages (e.g., ipyrad .snps.map)" )
parser . add_argument ( '-r' , metavar = "resolve" , dest = 'resolve' , type = int , default = 1 , help = "randomly resolve heterozygous sites (default=1)" )
parser . add_argument ( '-n' , metavar = "name" , dest = "name" , type = str , default = "test" , help = "output name prefix (default: 'test')" )
parser . add_argument ( '-o' , metavar = "workdir" , dest = "workdir" , type = str , default = "./analysis-tetrad" , help = "output directory (default: creates ./analysis-tetrad)" )
parser . add_argument ( '-t' , metavar = "starting_tree" , dest = "tree" , type = str , default = None , help = "newick file starting tree for equal splits sampling" )
parser . add_argument ( "-c" , metavar = "CPUs/cores" , dest = "cores" , type = int , default = 0 , help = "setting -c improves parallel efficiency with --MPI" )
parser . add_argument ( "-x" , metavar = "random_seed" , dest = "rseed" , type = int , default = None , help = "random seed for quartet sampling and/or bootstrapping" )
parser . add_argument ( '-d' , "--debug" , action = 'store_true' , help = "print lots more info to debugger: ipyrad_log.txt." )
parser . add_argument ( "--MPI" , action = 'store_true' , help = "connect to parallel CPUs across multiple nodes" )
parser . add_argument ( "--invariants" , action = 'store_true' , help = "save a (large) database of all invariants" )
parser . add_argument ( "--ipcluster" , metavar = "ipcluster" , dest = "ipcluster" , type = str , nargs = "?" , const = "default" , help = "connect to ipcluster profile (default: 'default')" )
# # if no args then return help message
if len ( sys . argv ) == 1 :
    parser . print_help ( )
    sys . exit ( 1 )
# # parse args
args = parser . parse_args ( )
# # RAISE errors right away for some bad argument combinations :
if args . method not in [ "random" , "equal" , "all" ] :
    raise IPyradWarningExit ( "  method argument (-m) must be one of" + """ "all", "random", or "equal.\n""" )
# # if ' random ' require nquarts argument
# if args . method = = ' random ' :
# if not args . nquartets :
# raise IPyradWarningExit ( \
# " Number of quartets ( - q ) is required with method = random \ n " )
# # if ' equal ' method require starting tree and nquarts
# if args . method = = ' equal ' :
# raise IPyradWarningExit ( \
# " The equal sampling method is currently for developers only . \ n " )
# if not args . nquartets :
# raise IPyradWarningExit ( \
# " Number of quartets ( - q ) is required with method = equal \ n " )
# if not args . tree :
# raise IPyradWarningExit ( \
# " Input guide tree ( - t ) is required with method = equal \ n " )
# # required args
if not any ( x in [ "seq" , "json" ] for x in vars ( args ) . keys ( ) ) :
    print ( """
    Bad arguments: tetrad command must include at least one of (-s or -j)
    """ )
    parser . print_help ( )
    sys . exit ( 1 )
return args