ipyrad/assemble/demultiplex.py
def barmatch2(data, tups, cutters, longbar, matchdict, fnum):
    """
    cleaner barmatch func...
    """

    ## how many reads to store before writing to disk
    waitchunk = int(1e6)
    ## pid name for this engine
    epid = os.getpid()

    ## counters for total reads, those with cutsite, and those that matched
    filestat = np.zeros(3, dtype=np.int)
    ## store observed sample matches
    samplehits = {}
    ## dictionaries to store first and second reads until writing to file
    dsort1 = {} 
    dsort2 = {} 
    ## dictionary for all bars matched in sample
    dbars = {} 

    ## fill for sample names
    for sname in data.barcodes:
        if "-technical-replicate-" in sname:
            sname = sname.rsplit("-technical-replicate", 1)[0]
        samplehits[sname] = 0
        dsort1[sname] = []
        dsort2[sname] = []
        dbars[sname] = set()

    ## store observed bars
    barhits = {}
    for barc in matchdict:
        barhits[barc] = 0

    ## store others
    misses = {}
    misses['_'] = 0

    ## build func for finding barcode
    getbarcode = get_barcode_func(data, longbar)

    ## get quart iterator of reads
    if tups[0].endswith(".gz"):
        ofunc = gzip.open
    else:
        ofunc = open

    ## create iterators 
    ofile1 = ofunc(tups[0], 'r')
    fr1 = iter(ofile1) 
    quart1 = itertools.izip(fr1, fr1, fr1, fr1)
    if tups[1]:
        ofile2 = ofunc(tups[1], 'r')
        fr2 = iter(ofile2)  
        quart2 = itertools.izip(fr2, fr2, fr2, fr2)
        quarts = itertools.izip(quart1, quart2)
    else:
        quarts = itertools.izip(quart1, iter(int, 1))

    ## go until end of the file
    while 1:
        try:
            read1, read2 = quarts.next()
            read1 = list(read1)
            filestat[0] += 1
        except StopIteration:
            break
    
        barcode = ""
        ## Get barcode_R2 and check for matching sample name
        if '3rad' in data.paramsdict["datatype"]:
            ## Here we're just reusing the findbcode function
            ## for R2, and reconfiguring the longbar tuple to have the
            ## maxlen for the R2 barcode
            ## Parse barcode. Use the parsing function selected above.
            barcode1 = find3radbcode(cutters=cutters, 
                                longbar=longbar, read1=read1)
            barcode2 = find3radbcode(cutters=cutters, 
                                longbar=(longbar[2], longbar[1]), read1=read2)
            barcode = barcode1 + "+" + barcode2
        else:
            ## Parse barcode. Uses the parsing function selected above.
            barcode = getbarcode(cutters, read1, longbar)
   
        ## find if it matches 
        sname_match = matchdict.get(barcode)

        if sname_match:
            #sample_index[filestat[0]-1] = snames.index(sname_match) + 1
            ## record who matched
            dbars[sname_match].add(barcode)
            filestat[1] += 1
            filestat[2] += 1
            samplehits[sname_match] += 1
            barhits[barcode] += 1
            if barcode in barhits:
                barhits[barcode] += 1
            else:
                barhits[barcode] = 1
    
            ## trim off barcode
            lenbar = len(barcode)
            if '3rad' in data.paramsdict["datatype"]:
                ## Iff 3rad trim the len of the first barcode
                lenbar = len(barcode1)
    
            if data.paramsdict["datatype"] == '2brad':
                overlen = len(cutters[0][0]) + lenbar + 1
                read1[1] = read1[1][:-overlen] + "\n"
                read1[3] = read1[3][:-overlen] + "\n"
            else:
                read1[1] = read1[1][lenbar:]
                read1[3] = read1[3][lenbar:]
    
            ## Trim barcode off R2 and append. Only 3rad datatype
            ## pays the cpu cost of splitting R2
            if '3rad' in data.paramsdict["datatype"]:
                read2 = list(read2)
                read2[1] = read2[1][len(barcode2):]
                read2[3] = read2[3][len(barcode2):]
    
            ## append to dsort
            dsort1[sname_match].append("".join(read1))
            if 'pair' in data.paramsdict["datatype"]:
                dsort2[sname_match].append("".join(read2))

        else:
            misses["_"] += 1
            if barcode:
                filestat[1] += 1

        ## how can we make it so all of the engines aren't trying to write to
        ## ~100-200 files all at the same time? This is the I/O limit we hit..
        ## write out at 100K to keep memory low. It is fine on HPC which can 
        ## write parallel, but regular systems might crash
        if not filestat[0] % waitchunk:
            ## write the remaining reads to file"
            writetofile(data, dsort1, 1, epid)
            if 'pair' in data.paramsdict["datatype"]:
                writetofile(data, dsort2, 2, epid)
            ## clear out dsorts
            for sample in data.barcodes:
                if "-technical-replicate-" in sname:
                    sname = sname.rsplit("-technical-replicate", 1)[0]
                dsort1[sname] = []
                dsort2[sname] = []
            ## reset longlist
            #longlist = np.zeros(waitchunk, dtype=np.uint32)                

    ## close open files
    ofile1.close()
    if tups[1]:
        ofile2.close()

    ## write the remaining reads to file
    writetofile(data, dsort1, 1, epid)
    if 'pair' in data.paramsdict["datatype"]:
        writetofile(data, dsort2, 2, epid)

    ## return stats in saved pickle b/c return_queue is too small
    ## and the size of the match dictionary can become quite large
    samplestats = [samplehits, barhits, misses, dbars]
    outname = os.path.join(data.dirs.fastqs, "tmp_{}_{}.p".format(epid, fnum))
    with open(outname, 'w') as wout:
        pickle.dump([filestat, samplestats], wout)

    return outname
def barmatch2(data, tups, cutters, longbar, matchdict, fnum):
    """
    cleaner barmatch func...
    """

    ## how many reads to store before writing to disk
    waitchunk = int(1e6)
    ## pid name for this engine
    epid = os.getpid()

    ## counters for total reads, those with cutsite, and those that matched
    filestat = np.zeros(3, dtype=np.int)
    ## store observed sample matches
    samplehits = {}
    ## dictionaries to store first and second reads until writing to file
    dsort1 = {} 
    dsort2 = {} 
    ## dictionary for all bars matched in sample
    dbars = {} 

    ## fill for sample names
    for sname in data.barcodes:
        if "-technical-replicate-" in sname:
            sname = sname.rsplit("-technical-replicate", 1)[0]
        samplehits[sname] = 0
        dsort1[sname] = []
        dsort2[sname] = []
        dbars[sname] = set()

    ## store observed bars
    barhits = {}
    for barc in matchdict:
        barhits[barc] = 0

    ## store others
    misses = {}
    misses['_'] = 0

    ## build func for finding barcode
    getbarcode = get_barcode_func(data, longbar)

    ## get quart iterator of reads
    if tups[0].endswith(".gz"):
        ofunc = gzip.open
    else:
        ofunc = open

    ## create iterators 
    ofile1 = ofunc(tups[0], 'r')
    fr1 = iter(ofile1) 
    quart1 = itertools.izip(fr1, fr1, fr1, fr1)
    if tups[1]:
        ofile2 = ofunc(tups[1], 'r')
        fr2 = iter(ofile2)  
        quart2 = itertools.izip(fr2, fr2, fr2, fr2)
        quarts = itertools.izip(quart1, quart2)
    else:
        quarts = itertools.izip(quart1, iter(int, 1))

    ## go until end of the file
    while 1:
        try:
            read1, read2 = quarts.next()
            read1 = list(read1)
            filestat[0] += 1
        except StopIteration:
            break
    
        barcode = ""
        ## Get barcode_R2 and check for matching sample name
        if '3rad' in data.paramsdict["datatype"]:
            ## Here we're just reusing the findbcode function
            ## for R2, and reconfiguring the longbar tuple to have the
            ## maxlen for the R2 barcode
            ## Parse barcode. Use the parsing function selected above.
            barcode1 = find3radbcode(cutters=cutters, 
                                longbar=longbar, read1=read1)
            barcode2 = find3radbcode(cutters=cutters, 
                                longbar=(longbar[2], longbar[1]), read1=read2)
            barcode = barcode1 + "+" + barcode2
        else:
            ## Parse barcode. Uses the parsing function selected above.
            barcode = getbarcode(cutters, read1, longbar)
   
        ## find if it matches 
        sname_match = matchdict.get(barcode)

        if sname_match:
            #sample_index[filestat[0]-1] = snames.index(sname_match) + 1
            ## record who matched
            dbars[sname_match].add(barcode)
            filestat[1] += 1
            filestat[2] += 1
            samplehits[sname_match] += 1
            barhits[barcode] += 1
            if barcode in barhits:
                barhits[barcode] += 1
            else:
                barhits[barcode] = 1
    
            ## trim off barcode
            lenbar = len(barcode)
            if '3rad' in data.paramsdict["datatype"]:
                ## Iff 3rad trim the len of the first barcode
                lenbar = len(barcode1)
    
            if data.paramsdict["datatype"] == '2brad':
                overlen = len(cutters[0][0]) + lenbar + 1
                read1[1] = read1[1][:-overlen] + "\n"
                read1[3] = read1[3][:-overlen] + "\n"
            else:
                read1[1] = read1[1][lenbar:]
                read1[3] = read1[3][lenbar:]
    
            ## Trim barcode off R2 and append. Only 3rad datatype
            ## pays the cpu cost of splitting R2
            if '3rad' in data.paramsdict["datatype"]:
                read2 = list(read2)
                read2[1] = read2[1][len(barcode2):]
                read2[3] = read2[3][len(barcode2):]
    
            ## append to dsort
            dsort1[sname_match].append("".join(read1))
            if 'pair' in data.paramsdict["datatype"]:
                dsort2[sname_match].append("".join(read2))

        else:
            misses["_"] += 1
            if barcode:
                filestat[1] += 1

        ## how can we make it so all of the engines aren't trying to write to
        ## ~100-200 files all at the same time? This is the I/O limit we hit..
        ## write out at 100K to keep memory low. It is fine on HPC which can 
        ## write parallel, but regular systems might crash
        if not filestat[0] % waitchunk:
            ## write the remaining reads to file"
            writetofile(data, dsort1, 1, epid)
            if 'pair' in data.paramsdict["datatype"]:
                writetofile(data, dsort2, 2, epid)
            ## clear out dsorts
            for sample in data.barcodes:
                if "-technical-replicate-" in sname:
                    sname = sname.rsplit("-technical-replicate", 1)[0]
                dsort1[sname] = []
                dsort2[sname] = []
            ## reset longlist
            #longlist = np.zeros(waitchunk, dtype=np.uint32)                

    ## close open files
    ofile1.close()
    if tups[1]:
        ofile2.close()

    ## write the remaining reads to file
    writetofile(data, dsort1, 1, epid)
    if 'pair' in data.paramsdict["datatype"]:
        writetofile(data, dsort2, 2, epid)

    ## return stats in saved pickle b/c return_queue is too small
    ## and the size of the match dictionary can become quite large
    samplestats = [samplehits, barhits, misses, dbars]
    outname = os.path.join(data.dirs.fastqs, "tmp_{}_{}.p".format(epid, fnum))
    with open(outname, 'w') as wout:
        pickle.dump([filestat, samplestats], wout)

    return outname
def barmatch2 ( data , tups , cutters , longbar , matchdict , fnum ) : ## how many reads to store before writing to disk waitchunk = int ( 1e6 ) ## pid name for this engine epid = os . getpid ( ) ## counters for total reads, those with cutsite, and those that matched filestat = np . zeros ( 3 , dtype = np . int ) ## store observed sample matches samplehits = { } ## dictionaries to store first and second reads until writing to file dsort1 = { } dsort2 = { } ## dictionary for all bars matched in sample dbars = { } ## fill for sample names for sname in data . barcodes : if "-technical-replicate-" in sname : sname = sname . rsplit ( "-technical-replicate" , 1 ) [ 0 ] samplehits [ sname ] = 0 dsort1 [ sname ] = [ ] dsort2 [ sname ] = [ ] dbars [ sname ] = set ( ) ## store observed bars barhits = { } for barc in matchdict : barhits [ barc ] = 0 ## store others misses = { } misses [ '_' ] = 0 ## build func for finding barcode getbarcode = get_barcode_func ( data , longbar ) ## get quart iterator of reads if tups [ 0 ] . endswith ( ".gz" ) : ofunc = gzip . open else : ofunc = open ## create iterators  ofile1 = ofunc ( tups [ 0 ] , 'r' ) fr1 = iter ( ofile1 ) quart1 = itertools . izip ( fr1 , fr1 , fr1 , fr1 ) if tups [ 1 ] : ofile2 = ofunc ( tups [ 1 ] , 'r' ) fr2 = iter ( ofile2 ) quart2 = itertools . izip ( fr2 , fr2 , fr2 , fr2 ) quarts = itertools . izip ( quart1 , quart2 ) else : quarts = itertools . izip ( quart1 , iter ( int , 1 ) ) ## go until end of the file while 1 : try : read1 , read2 = quarts . next ( ) read1 = list ( read1 ) filestat [ 0 ] += 1 except StopIteration : break barcode = "" ## Get barcode_R2 and check for matching sample name if '3rad' in data . paramsdict [ "datatype" ] : ## Here we're just reusing the findbcode function ## for R2, and reconfiguring the longbar tuple to have the ## maxlen for the R2 barcode ## Parse barcode. Use the parsing function selected above. barcode1 = find3radbcode ( cutters = cutters , longbar = longbar , read1 = read1 ) barcode2 = find3radbcode ( cutters = cutters , longbar = ( longbar [ 2 ] , longbar [ 1 ] ) , read1 = read2 ) barcode = barcode1 + "+" + barcode2 else : ## Parse barcode. Uses the parsing function selected above. barcode = getbarcode ( cutters , read1 , longbar ) ## find if it matches  sname_match = matchdict . get ( barcode ) if sname_match : #sample_index[filestat[0]-1] = snames.index(sname_match) + 1 ## record who matched dbars [ sname_match ] . add ( barcode ) filestat [ 1 ] += 1 filestat [ 2 ] += 1 samplehits [ sname_match ] += 1 barhits [ barcode ] += 1 if barcode in barhits : barhits [ barcode ] += 1 else : barhits [ barcode ] = 1 ## trim off barcode lenbar = len ( barcode ) if '3rad' in data . paramsdict [ "datatype" ] : ## Iff 3rad trim the len of the first barcode lenbar = len ( barcode1 ) if data . paramsdict [ "datatype" ] == '2brad' : overlen = len ( cutters [ 0 ] [ 0 ] ) + lenbar + 1 read1 [ 1 ] = read1 [ 1 ] [ : - overlen ] + "\n" read1 [ 3 ] = read1 [ 3 ] [ : - overlen ] + "\n" else : read1 [ 1 ] = read1 [ 1 ] [ lenbar : ] read1 [ 3 ] = read1 [ 3 ] [ lenbar : ] ## Trim barcode off R2 and append. Only 3rad datatype ## pays the cpu cost of splitting R2 if '3rad' in data . paramsdict [ "datatype" ] : read2 = list ( read2 ) read2 [ 1 ] = read2 [ 1 ] [ len ( barcode2 ) : ] read2 [ 3 ] = read2 [ 3 ] [ len ( barcode2 ) : ] ## append to dsort dsort1 [ sname_match ] . append ( "" . join ( read1 ) ) if 'pair' in data . paramsdict [ "datatype" ] : dsort2 [ sname_match ] . append ( "" . join ( read2 ) ) else : misses [ "_" ] += 1 if barcode : filestat [ 1 ] += 1 ## how can we make it so all of the engines aren't trying to write to ## ~100-200 files all at the same time? This is the I/O limit we hit.. ## write out at 100K to keep memory low. It is fine on HPC which can  ## write parallel, but regular systems might crash if not filestat [ 0 ] % waitchunk : ## write the remaining reads to file" writetofile ( data , dsort1 , 1 , epid ) if 'pair' in data . paramsdict [ "datatype" ] : writetofile ( data , dsort2 , 2 , epid ) ## clear out dsorts for sample in data . barcodes : if "-technical-replicate-" in sname : sname = sname . rsplit ( "-technical-replicate" , 1 ) [ 0 ] dsort1 [ sname ] = [ ] dsort2 [ sname ] = [ ] ## reset longlist #longlist = np.zeros(waitchunk, dtype=np.uint32)                 ## close open files ofile1 . close ( ) if tups [ 1 ] : ofile2 . close ( ) ## write the remaining reads to file writetofile ( data , dsort1 , 1 , epid ) if 'pair' in data . paramsdict [ "datatype" ] : writetofile ( data , dsort2 , 2 , epid ) ## return stats in saved pickle b/c return_queue is too small ## and the size of the match dictionary can become quite large samplestats = [ samplehits , barhits , misses , dbars ] outname = os . path . join ( data . dirs . fastqs , "tmp_{}_{}.p" . format ( epid , fnum ) ) with open ( outname , 'w' ) as wout : pickle . dump ( [ filestat , samplestats ] , wout ) return outname
cleaner barmatch func...
cleaner barmatch func ...