ipyrad/assemble/consens_se.py
def newconsensus(data, sample, tmpchunk, optim):
    """ 
    new faster replacement to consensus 
    """
    ## do reference map funcs?
    isref = "reference" in data.paramsdict["assembly_method"]

    ## temporarily store the mean estimates to Assembly
    data._este = data.stats.error_est.mean()
    data._esth = data.stats.hetero_est.mean()

    ## get number relative to tmp file
    tmpnum = int(tmpchunk.split(".")[-1])

    ## prepare data for reading
    clusters = open(tmpchunk, 'rb')
    pairdealer = itertools.izip(*[iter(clusters)]*2)    
    maxlen = data._hackersonly["max_fragment_length"]

    ## write to tmp cons to file to be combined later
    consenshandle = os.path.join(
        data.dirs.consens, sample.name+"_tmpcons."+str(tmpnum))
    tmp5 = consenshandle.replace("_tmpcons.", "_tmpcats.")
    with h5py.File(tmp5, 'w') as io5:
        io5.create_dataset("cats", (optim, maxlen, 4), dtype=np.uint32)
        io5.create_dataset("alls", (optim, ), dtype=np.uint8)
        io5.create_dataset("chroms", (optim, 3), dtype=np.int64)

        ## local copies to use to fill the arrays
        catarr = io5["cats"][:]
        nallel = io5["alls"][:]
        refarr = io5["chroms"][:]

    ## if reference-mapped then parse the fai to get index number of chroms
    if isref:
        fai = pd.read_csv(data.paramsdict["reference_sequence"] + ".fai", 
                names=['scaffold', 'size', 'sumsize', 'a', 'b'],
                sep="\t")
        faidict = {j:i for i,j in enumerate(fai.scaffold)}

    ## store data for stats counters
    counters = {"name" : tmpnum,
                "heteros": 0,
                "nsites" : 0,
                "nconsens" : 0}

    ## store data for what got filtered
    filters = {"depth" : 0,
               "maxh" : 0,
               "maxn" : 0}

    ## store data for writing
    storeseq = {}

    ## set max limits
    if 'pair' in data.paramsdict["datatype"]:
        maxhet = sum(data.paramsdict["max_Hs_consens"])
        maxn = sum(data.paramsdict["max_Ns_consens"])
    else:
        maxhet = data.paramsdict["max_Hs_consens"][0]
        maxn = data.paramsdict["max_Ns_consens"][0]

    ## load the refmap dictionary if refmapping
    done = 0
    while not done:
        try:
            done, chunk = clustdealer(pairdealer, 1)
        except IndexError:
            raise IPyradError("clustfile formatting error in %s", chunk)

        if chunk:
            ## get names and seqs
            piece = chunk[0].strip().split("\n")
            names = piece[0::2]
            seqs = piece[1::2]

            ## pull replicate read info from seqs
            reps = [int(sname.split(";")[-2][5:]) for sname in names]

            ## IF this is a reference mapped read store the chrom and pos info
            ## -1 defaults to indicating an anonymous locus, since we are using
            ## the faidict as 0 indexed. If chrompos fails it defaults to -1
            ref_position = (-1, 0, 0)
            if isref:
                try:
                    ## parse position from name string
                    name, _, _ = names[0].rsplit(";", 2)
                    chrom, pos0, pos1 = name.rsplit(":", 2)
                    
                    ## pull idx from .fai reference dict 
                    chromint = faidict[chrom] + 1
                    ref_position = (int(chromint), int(pos0), int(pos1))
                    
                except Exception as inst:
                    LOGGER.debug("Reference sequence chrom/pos failed for {}".format(names[0]))
                    LOGGER.debug(inst)
                    
            ## apply read depth filter
            if nfilter1(data, reps):

                ## get stacks of base counts
                sseqs = [list(seq) for seq in seqs]
                arrayed = np.concatenate(
                    [[seq]*rep for seq, rep in zip(sseqs, reps)])
                arrayed = arrayed[:, :maxlen]
                
                ## get consens call for each site, applies paralog-x-site filter
                #consens = np.apply_along_axis(basecall, 0, arrayed, data)
                consens = basecaller(
                    arrayed, 
                    data.paramsdict["mindepth_majrule"], 
                    data.paramsdict["mindepth_statistical"],
                    data._esth, 
                    data._este,
                    )

                ## apply a filter to remove low coverage sites/Ns that
                ## are likely sequence repeat errors. This is only applied to
                ## clusters that already passed the read-depth filter (1)
                if "N" in consens:
                    try:
                        consens, arrayed = removerepeats(consens, arrayed)

                    except ValueError as _:
                        LOGGER.info("Caught a bad chunk w/ all Ns. Skip it.")
                        continue

                ## get hetero sites
                hidx = [i for (i, j) in enumerate(consens) \
                            if j in list("RKSYWM")]
                nheteros = len(hidx)
                
                ## filter for max number of hetero sites
                if nfilter2(nheteros, maxhet):
                    ## filter for maxN, & minlen
                    if nfilter3(consens, maxn):
                        ## counter right now
                        current = counters["nconsens"]
                        ## get N alleles and get lower case in consens
                        consens, nhaps = nfilter4(consens, hidx, arrayed)
                        ## store the number of alleles observed
                        nallel[current] = nhaps

                        ## store a reduced array with only CATG
                        catg = np.array(\
                            [np.sum(arrayed == i, axis=0)  \
                            for i in list("CATG")],
                            dtype='uint32').T
                        catarr[current, :catg.shape[0], :] = catg
                        refarr[current] = ref_position

                        ## store the seqdata for tmpchunk
                        storeseq[counters["name"]] = "".join(list(consens))
                        counters["name"] += 1
                        counters["nconsens"] += 1
                        counters["heteros"] += nheteros
                    else:
                        #LOGGER.debug("@haplo")
                        filters['maxn'] += 1
                else:
                    #LOGGER.debug("@hetero")
                    filters['maxh'] += 1
            else:
                #LOGGER.debug("@depth")
                filters['depth'] += 1
                
    ## close infile io
    clusters.close()

    ## write final consens string chunk
    if storeseq:
        with open(consenshandle, 'wb') as outfile:
            outfile.write("\n".join([">"+sample.name+"_"+str(key)+"\n"+\
                                   str(storeseq[key]) for key in storeseq]))

    ## write to h5 array, this can be a bit slow on big data sets and is not 
    ## currently convered by progressbar movement.
    with h5py.File(tmp5, 'a') as io5:
        io5["cats"][:] = catarr
        io5["alls"][:] = nallel
        io5["chroms"][:] = refarr
    del catarr
    del nallel
    del refarr

    ## return stats
    counters['nsites'] = sum([len(i) for i in storeseq.itervalues()])
    return counters, filters
def newconsensus(data, sample, tmpchunk, optim):
    """ 
    new faster replacement to consensus 
    """
    ## do reference map funcs?
    isref = "reference" in data.paramsdict["assembly_method"]

    ## temporarily store the mean estimates to Assembly
    data._este = data.stats.error_est.mean()
    data._esth = data.stats.hetero_est.mean()

    ## get number relative to tmp file
    tmpnum = int(tmpchunk.split(".")[-1])

    ## prepare data for reading
    clusters = open(tmpchunk, 'rb')
    pairdealer = itertools.izip(*[iter(clusters)]*2)    
    maxlen = data._hackersonly["max_fragment_length"]

    ## write to tmp cons to file to be combined later
    consenshandle = os.path.join(
        data.dirs.consens, sample.name+"_tmpcons."+str(tmpnum))
    tmp5 = consenshandle.replace("_tmpcons.", "_tmpcats.")
    with h5py.File(tmp5, 'w') as io5:
        io5.create_dataset("cats", (optim, maxlen, 4), dtype=np.uint32)
        io5.create_dataset("alls", (optim, ), dtype=np.uint8)
        io5.create_dataset("chroms", (optim, 3), dtype=np.int64)

        ## local copies to use to fill the arrays
        catarr = io5["cats"][:]
        nallel = io5["alls"][:]
        refarr = io5["chroms"][:]

    ## if reference-mapped then parse the fai to get index number of chroms
    if isref:
        fai = pd.read_csv(data.paramsdict["reference_sequence"] + ".fai", 
                names=['scaffold', 'size', 'sumsize', 'a', 'b'],
                sep="\t")
        faidict = {j:i for i,j in enumerate(fai.scaffold)}

    ## store data for stats counters
    counters = {"name" : tmpnum,
                "heteros": 0,
                "nsites" : 0,
                "nconsens" : 0}

    ## store data for what got filtered
    filters = {"depth" : 0,
               "maxh" : 0,
               "maxn" : 0}

    ## store data for writing
    storeseq = {}

    ## set max limits
    if 'pair' in data.paramsdict["datatype"]:
        maxhet = sum(data.paramsdict["max_Hs_consens"])
        maxn = sum(data.paramsdict["max_Ns_consens"])
    else:
        maxhet = data.paramsdict["max_Hs_consens"][0]
        maxn = data.paramsdict["max_Ns_consens"][0]

    ## load the refmap dictionary if refmapping
    done = 0
    while not done:
        try:
            done, chunk = clustdealer(pairdealer, 1)
        except IndexError:
            raise IPyradError("clustfile formatting error in %s", chunk)

        if chunk:
            ## get names and seqs
            piece = chunk[0].strip().split("\n")
            names = piece[0::2]
            seqs = piece[1::2]

            ## pull replicate read info from seqs
            reps = [int(sname.split(";")[-2][5:]) for sname in names]

            ## IF this is a reference mapped read store the chrom and pos info
            ## -1 defaults to indicating an anonymous locus, since we are using
            ## the faidict as 0 indexed. If chrompos fails it defaults to -1
            ref_position = (-1, 0, 0)
            if isref:
                try:
                    ## parse position from name string
                    name, _, _ = names[0].rsplit(";", 2)
                    chrom, pos0, pos1 = name.rsplit(":", 2)
                    
                    ## pull idx from .fai reference dict 
                    chromint = faidict[chrom] + 1
                    ref_position = (int(chromint), int(pos0), int(pos1))
                    
                except Exception as inst:
                    LOGGER.debug("Reference sequence chrom/pos failed for {}".format(names[0]))
                    LOGGER.debug(inst)
                    
            ## apply read depth filter
            if nfilter1(data, reps):

                ## get stacks of base counts
                sseqs = [list(seq) for seq in seqs]
                arrayed = np.concatenate(
                    [[seq]*rep for seq, rep in zip(sseqs, reps)])
                arrayed = arrayed[:, :maxlen]
                
                ## get consens call for each site, applies paralog-x-site filter
                #consens = np.apply_along_axis(basecall, 0, arrayed, data)
                consens = basecaller(
                    arrayed, 
                    data.paramsdict["mindepth_majrule"], 
                    data.paramsdict["mindepth_statistical"],
                    data._esth, 
                    data._este,
                    )

                ## apply a filter to remove low coverage sites/Ns that
                ## are likely sequence repeat errors. This is only applied to
                ## clusters that already passed the read-depth filter (1)
                if "N" in consens:
                    try:
                        consens, arrayed = removerepeats(consens, arrayed)

                    except ValueError as _:
                        LOGGER.info("Caught a bad chunk w/ all Ns. Skip it.")
                        continue

                ## get hetero sites
                hidx = [i for (i, j) in enumerate(consens) \
                            if j in list("RKSYWM")]
                nheteros = len(hidx)
                
                ## filter for max number of hetero sites
                if nfilter2(nheteros, maxhet):
                    ## filter for maxN, & minlen
                    if nfilter3(consens, maxn):
                        ## counter right now
                        current = counters["nconsens"]
                        ## get N alleles and get lower case in consens
                        consens, nhaps = nfilter4(consens, hidx, arrayed)
                        ## store the number of alleles observed
                        nallel[current] = nhaps

                        ## store a reduced array with only CATG
                        catg = np.array(\
                            [np.sum(arrayed == i, axis=0)  \
                            for i in list("CATG")],
                            dtype='uint32').T
                        catarr[current, :catg.shape[0], :] = catg
                        refarr[current] = ref_position

                        ## store the seqdata for tmpchunk
                        storeseq[counters["name"]] = "".join(list(consens))
                        counters["name"] += 1
                        counters["nconsens"] += 1
                        counters["heteros"] += nheteros
                    else:
                        #LOGGER.debug("@haplo")
                        filters['maxn'] += 1
                else:
                    #LOGGER.debug("@hetero")
                    filters['maxh'] += 1
            else:
                #LOGGER.debug("@depth")
                filters['depth'] += 1
                
    ## close infile io
    clusters.close()

    ## write final consens string chunk
    if storeseq:
        with open(consenshandle, 'wb') as outfile:
            outfile.write("\n".join([">"+sample.name+"_"+str(key)+"\n"+\
                                   str(storeseq[key]) for key in storeseq]))

    ## write to h5 array, this can be a bit slow on big data sets and is not 
    ## currently convered by progressbar movement.
    with h5py.File(tmp5, 'a') as io5:
        io5["cats"][:] = catarr
        io5["alls"][:] = nallel
        io5["chroms"][:] = refarr
    del catarr
    del nallel
    del refarr

    ## return stats
    counters['nsites'] = sum([len(i) for i in storeseq.itervalues()])
    return counters, filters
new faster replacement to consensus