ipyrad/assemble/refmap.py
def fetch_cluster_pairs(data, samfile, chrom, rstart, rend):
    """ 
    Builds a paired cluster from the refmapped data.
    """
    ## 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
def fetch_cluster_pairs(data, samfile, chrom, rstart, rend):
    """ 
    Builds a paired cluster from the refmapped data.
    """
    ## 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
def fetch_cluster_pairs ( data , samfile , chrom , rstart , rend ) : ## 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
Builds a paired cluster from the refmapped data.
Builds a paired cluster from the refmapped data .