remove poorly aligned positions and divergent regions with Gblocks
173560
def run_gblocks(align_fasta_file, **kwargs):
    
    cl = GblocksCommandline(aln_file=align_fasta_file, **kwargs)
    r, e = cl.run()

    print("Gblocks:", cl, file=sys.stderr)

    if e:
        print("***Gblocks could not run", file=sys.stderr)
        return None
    else:
        print(r, file=sys.stderr)
        alignp = re.sub(r'.*Gblocks alignment:.*\(([0-9]{1,3}) %\).*', \
            r'\1', r, flags=re.DOTALL)
        alignp = int(alignp)
        if alignp <= 10:
            print("** WARNING ** Only %s %% positions retained by Gblocks. " \
                "Results aborted. Using original alignment instead.\n" % alignp, file=sys.stderr)
            return None
        else:
            return align_fasta_file+"-gb"
784921
def msgblock(key, text, side='|'):
    
    blocked_text = ''.join(
        [' + --- ', key, ' ---\n'] +
        [' ' + side + ' ' + line + '\n' for line in text.split('\n')] +
        [' L ___ ', key, ' ___\n']
    )
    return blocked_text
1168577
def block_sep1(self, Y):
        
        """

        return Y[(slice(None),)*self.blkaxis + (slice(self.blkidx, None),)]
1097381
def _tiles_from_cleaned_bounds(tp, bounds, zoom):
    """"""
    lb = _tile_from_xy(tp, bounds.left, bounds.bottom, zoom, on_edge_use="rt")
    rt = _tile_from_xy(tp, bounds.right, bounds.top, zoom, on_edge_use="lb")
    for tile_id in product([zoom], range(rt.row, lb.row + 1), range(lb.col, rt.col + 1)):
        yield tp.tile(*tile_id)
1350092
def _coarsen_reshape(self, windows, boundary, side):
        
        if not utils.is_dict_like(boundary):
            boundary = {d: boundary for d in windows.keys()}

        if not utils.is_dict_like(side):
            side = {d: side for d in windows.keys()}

        # remove unrelated dimensions
        boundary = {k: v for k, v in boundary.items() if k in windows}
        side = {k: v for k, v in side.items() if k in windows}

        for d, window in windows.items():
            if window <= 0:
                raise ValueError('window must be > 0. Given {}'.format(window))

        variable = self
        for d, window in windows.items():
            # trim or pad the object
            size = variable.shape[self._get_axis_num(d)]
            n = int(size / window)
            if boundary[d] == 'exact':
                if n * window != size:
                    raise ValueError(
                        'Could not coarsen a dimension of size {} with '
                        'window {}'.format(size, window))
            elif boundary[d] == 'trim':
                if side[d] == 'left':
                    variable = variable.isel({d: slice(0, window * n)})
                else:
                    excess = size - window * n
                    variable = variable.isel({d: slice(excess, None)})
            elif boundary[d] == 'pad':  # pad
                pad = window * n - size
                if pad < 0:
                    pad += window
                if side[d] == 'left':
                    pad_widths = {d: (0, pad)}
                else:
                    pad_widths = {d: (pad, 0)}
                variable = variable.pad_with_fill_value(pad_widths)
            else:
                raise TypeError(
                    "{} is invalid for boundary. Valid option is 'exact', "
                    "'trim' and 'pad'".format(boundary[d]))

        shape = []
        axes = []
        axis_count = 0
        for i, d in enumerate(variable.dims):
            if d in windows:
                size = variable.shape[i]
                shape.append(int(size / windows[d]))
                shape.append(windows[d])
                axis_count += 1
                axes.append(i + axis_count)
            else:
                shape.append(variable.shape[i])

        return variable.data.reshape(shape), tuple(axes)
237484
def gblocks(self,
                new_path = None,
                seq_type = 'nucl' or 'prot'):
        """"""
        # Temporary path #
        if new_path is None: final = self.__class__(new_temp_path())
        else:                final = self.__class__(new_path)
        # Mapping every sequence name with a random name #
        orig_name_to_temp = {seq.description: 'name' + str(i) for i,seq in enumerate(self)}
        temp_name_to_orig = {v: k for k, v in orig_name_to_temp.items()}
        # Rename every sequence with a random name #
        temp_fasta = self.rename_sequences(orig_name_to_temp)
        # Options #
        if seq_type == 'nucl': t_option = "-t=d"
        if seq_type == 'prot': t_option = "-t=p"
        # Run it #
        result = sh.gblocks91(temp_fasta.path, t_option, '-p=n', "-b4=3", "-b3=20", "-b5=a", _ok_code=[0,1])
        created_file = temp_fasta.path + '-gb'
        assert os.path.exists(created_file)
        # Check errors #
        if "Execution terminated" in result.stdout: raise Exception("gblocks crashed again.")
        # Back #
        temp_fasta.rename_sequences(temp_name_to_orig, final)
        # Return #
        return final
1282821
def remove_exclude_regions(orig_bed, base_file, items, remove_entire_feature=False):
    
    from bcbio.structural import shared as sshared
    out_bed = os.path.join("%s-noexclude.bed" % (utils.splitext_plus(base_file)[0]))
    if not utils.file_uptodate(out_bed, orig_bed):
        exclude_bed = sshared.prepare_exclude_file(items, base_file)
        with file_transaction(items[0], out_bed) as tx_out_bed:
            pybedtools.BedTool(orig_bed).subtract(pybedtools.BedTool(exclude_bed),
                                                  A=remove_entire_feature, nonamecheck=True).saveas(tx_out_bed)
    if utils.file_exists(out_bed):
        return out_bed
    else:
        return orig_bed
966079
def _replace_bm(self):
        """"""
        self._block_matcher = cv2.StereoBM(preset=self._bm_preset,
                                          ndisparities=self._search_range,
                                          SADWindowSize=self._window_size)
26149
def clean_locks(root=None):
    '''
    
    '''
    LCK = "removed"
    out = {LCK: 0}
    locks = os.path.join(root, os.path.relpath(LOCKS, os.path.sep)) if root else LOCKS
    if not os.path.exists(locks):
        return out

    for node in __zypper__(root=root).xml.call('cl').getElementsByTagName("message"):
        text = node.childNodes[0].nodeValue.lower()
        if text.startswith(LCK):
            out[LCK] = text.split(" ")[1]
            break

    return out
620613
def _vertical_segment(old_offs, new_offs, spacing, radii):
    '''
    '''
    return np.array(((new_offs[0] - radii[0], old_offs[1] + spacing[1]),
                     (new_offs[0] - radii[1], new_offs[1]),
                     (new_offs[0] + radii[1], new_offs[1]),
                     (new_offs[0] + radii[0], old_offs[1] + spacing[1])))
68530
def align(self):
        """"""

        bits = self._bits
        self._buffer = 0
        self._bits = 0
        return bits
98851
public void removeMask(long[] mask) {
    int dest = BitsUtil.nextSetBit(mask, 0);
    if(dest < 0) {
      return;
    }
    int src = BitsUtil.nextSetBit(mask, dest);
    while(src < numEntries) {
      if(!BitsUtil.get(mask, src)) {
        entries[dest] = entries[src];
        dest++;
      }
      src++;
    }
    int rm = src - dest;
    while(dest < numEntries) {
      entries[dest] = null;
      dest++;
    }
    numEntries -= rm;
  }
928679
public static void alignDatanodes(DatanodeInfo[] dstLocs,
      DatanodeInfo[] srcLocs) {
    for (int i = 0; i < dstLocs.length; i++) {
      for (int j = 0; j < srcLocs.length; j++) {
        if (i == j)
          continue;
        if (dstLocs[i].equals(srcLocs[j])) {
          if (i < j) {
            swap(i, j, srcLocs);
          } else {
            swap(i, j, dstLocs);
          }
          break;
        }
      }
    }
  }
1149938
function cleanBogus( block ) {
			var bogus = [];
			var last = getLast( block ), node, previous;

			if ( last ) {
				// Check for bogus at the end of this block.
				// e.g. <p>foo<br /></p>
				maybeBogus( last, 1 ) && bogus.push( last );

				while ( last ) {
					// Check for bogus at the end of any pseudo block contained.
					if ( isBlockBoundary( last ) && ( node = getPrevious( last ) ) && maybeBogus( node ) ) {
						// Bogus must have inline proceeding, instead single BR between two blocks,
						// is considered as filler, e.g. <hr /><br /><hr />
						if ( ( previous = getPrevious( node ) ) && !isBlockBoundary( previous ) )
							bogus.push( node );
						// Convert the filler into appropriate form.
						else {
							createFiller( isOutput ).insertAfter( node );
							node.remove();
						}
					}

					last = last.previous;
				}
			}

			// Now remove all bogus collected from above.
			for ( var i = 0 ; i < bogus.length ; i++ )
				bogus[ i ].remove();
		}
1020755
def correctGrid(self, img, grid):
        '''
        
        '''

        self.img = imread(img)
        h = self.homography  # TODO: cleanup only needed to get newBorder attr.

        if self.opts['do_correctIntensity']:
            self.img = self.img / self._getTiltFactor(self.img.shape)

        s0, s1 = grid.shape[:2]
        n0, n1 = s0 - 1, s1 - 1

        snew = self._newBorders
        b = self.opts['border']

        sx, sy = (snew[0] - 2 * b) // n0, (snew[1] - 2 * b) // n1

        out = np.empty(snew[::-1], dtype=self.img.dtype)

        def warp(ix, iy, objP, outcut):
            shape = outcut.shape[::-1]
            quad = grid[ix:ix + 2,
                        iy:iy + 2].reshape(4, 2)[np.array([0, 2, 3, 1])]
            hcell = cv2.getPerspectiveTransform(
                quad.astype(np.float32), objP)
            cv2.warpPerspective(self.img, hcell, shape, outcut,
                                flags=cv2.INTER_LANCZOS4,
                                **self.opts['cv2_opts'])
            return quad

        objP = np.array([[0, 0],
                         [sx, 0],
                         [sx, sy],
                         [0, sy]], dtype=np.float32)
        # INNER CELLS
        for ix in range(1, n0 - 1):
            for iy in range(1, n1 - 1):
                sub = out[iy * sy + b: (iy + 1) * sy + b,
                          ix * sx + b: (ix + 1) * sx + b]
#                 warp(ix, iy, objP, sub)

                shape = sub.shape[::-1]
                quad = grid[ix:ix + 2,
                            iy:iy + 2].reshape(4, 2)[np.array([0, 2, 3, 1])]
#                 print(quad, objP)

                hcell = cv2.getPerspectiveTransform(
                    quad.astype(np.float32), objP)
                cv2.warpPerspective(self.img, hcell, shape, sub,
                                    flags=cv2.INTER_LANCZOS4,
                                    **self.opts['cv2_opts'])

#         return out
        # TOP CELLS
        objP[:, 1] += b
        for ix in range(1, n0 - 1):
            warp(ix, 0, objP, out[: sy + b,
                                  ix * sx + b: (ix + 1) * sx + b])
        # BOTTOM CELLS
        objP[:, 1] -= b
        for ix in range(1, n0 - 1):
            iy = (n1 - 1)
            y = iy * sy + b
            x = ix * sx + b
            warp(ix, iy, objP, out[y: y + sy + b, x: x + sx])
        # LEFT CELLS
        objP[:, 0] += b
        for iy in range(1, n1 - 1):
            y = iy * sy + b
            warp(0, iy, objP, out[y: y + sy, : sx + b])
        # RIGHT CELLS
        objP[:, 0] -= b
        ix = (n0 - 1)
        x = ix * sx + b
        for iy in range(1, n1 - 1):
            y = iy * sy + b
            warp(ix, iy, objP, out[y: y + sy, x: x + sx + b])
        # BOTTOM RIGHT CORNER
        warp(n0 - 1, n1 - 1, objP, out[-sy - b - 1:, x: x + sx + b])
#         #TOP LEFT CORNER
        objP += (b, b)
        warp(0, 0, objP, out[0: sy + b, 0: sx + b])
        # TOP RIGHT CORNER
        objP[:, 0] -= b
#         x = (n0-1)*sx+b
        warp(n0 - 1, 0, objP, out[: sy + b, x: x + sx + b])
#         #BOTTOM LEFT CORNER
        objP += (b, -b)
        warp(0, n1 - 1, objP, out[-sy - b - 1:, : sx + b])
        return out
102850
def erase_line (self): # <ESC>[2K
        ''''''

        self.fill_region (self.cur_r, 1, self.cur_r, self.cols)
1158922
def has_aligned_reads(align_bam, region=None):
    
    import pybedtools
    if region is not None:
        if isinstance(region, six.string_types) and os.path.isfile(region):
            regions = [tuple(r) for r in pybedtools.BedTool(region)]
        else:
            regions = [region]
    with pysam.Samfile(align_bam, "rb") as cur_bam:
        if region is not None:
            for region in regions:
                if isinstance(region, six.string_types):
                    for item in cur_bam.fetch(str(region)):
                        return True
                else:
                    for item in cur_bam.fetch(str(region[0]), int(region[1]), int(region[2])):
                        return True
        else:
            for item in cur_bam:
                if not item.is_unmapped:
                    return True
    return False