Source code for cnvlib.gary

"""A generic array of genomic positions."""
from __future__ import print_function, absolute_import, division

import logging
import sys
import warnings

import numpy as np
import pandas as pd

from . import core, ngfrills

[docs]class GenomicArray(object): """An array of genomic intervals. Base class for genomic data structures. Can represent most BED-like tabular formats with arbitrary additional columns. """ _required_columns = ("chromosome", "start", "end") _required_dtypes = ("string", "int", "int") def __init__(self, data_table, meta_dict=None): # Validation if len(data_table): if not all(c in data_table.columns for c in self._required_columns): raise ValueError("data table must have at least columns " + repr(self._required_columns)) # Ensure chromosomes are strings to avoid integer conversion of 1, 2... if not isinstance(data_table.chromosome.iat[0], basestring): data_table["chromosome"] = (data_table["chromosome"] .astype("string")) elif not isinstance(data_table, pd.DataFrame): # Empty but conformant table data_table = self._make_blank() = data_table self.meta = (dict(meta_dict) if meta_dict is not None and len(meta_dict) else {}) @staticmethod
[docs] def row2label(row): return "{}:{}-{}".format(row.chromosome, row.start, row.end)
@classmethod def _make_blank(cls): """Create an empty dataframe with the columns required by this class.""" table = pd.DataFrame({key: [] for key in cls._required_columns}) for col, dtype in zip(cls._required_columns, cls._required_dtypes): table[col] = table[col].astype(dtype) return table @classmethod
[docs] def from_columns(cls, columns, meta_dict=None): """Create a new instance from column arrays, given as a dict.""" table = pd.DataFrame.from_dict(columns) ary = cls(table, meta_dict) ary.sort_columns() return ary
[docs] def from_rows(cls, rows, columns=None, meta_dict=None): """Create a new instance from a list of rows, as tuples or arrays.""" if columns is None: columns = cls._required_columns table = pd.DataFrame.from_records(rows, columns=columns) return cls(table, meta_dict)
[docs] def as_columns(self, **columns): """Wrap the named columns in this instance's metadata.""" return self.__class__.from_columns(columns, self.meta)
# return self.__class__([:, columns], self.meta.copy())
[docs] def as_dataframe(self, dframe): """Wrap the given pandas dataframe in this instance's metadata.""" return self.__class__(dframe.reset_index(drop=True), self.meta.copy())
# def as_index(self, index): # """Subset with fancy/boolean indexing; reuse this instance's metadata.""" # """Extract rows by indices, reusing this instance's metadata.""" # if isinstance(index, (int, slice)): # return self.__class__([index], self.meta.copy()) # else: # return self.__class__([index], self.meta.copy())
[docs] def as_rows(self, rows): """Wrap the given rows in this instance's metadata.""" return self.from_rows(rows,, meta_dict=self.meta)
# Container behaviour def __eq__(self, other): return (isinstance(other, self.__class__) and def __len__(self): return len( def __contains__(self, key): return key in def __getitem__(self, index): """Access a portion of the data. Cases: - single integer: a row, as pd.Series - string row name: a column, as pd.Series - a boolean array: masked rows, as_dataframe - tuple of integers: selected rows, as_dataframe """ if isinstance(index, int): # A single row return[index] # return self.as_dataframe([index:index+1]) elif isinstance(index, basestring): # A column, by name return[index] elif (isinstance(index, tuple) and len(index) == 2 and index[1] in # Row index, column index -> cell value return[index] elif isinstance(index, slice): # return self.as_dataframe( return self.as_dataframe([index]) else: # Iterable -- selected row indices or boolean array, probably try: if isinstance(index, type(None)) or len(index) == 0: empty = pd.DataFrame( return self.as_dataframe(empty) except TypeError: raise TypeError("object of type %r " % type(index) + "cannot be used as an index into a " + self.__class__.__name__) return self.as_dataframe([index]) # return self.as_dataframe( def __setitem__(self, index, value): """Assign to a portion of the data. """ #[index] = value if isinstance(index, int):[index] = value elif isinstance(index, basestring):[index] = value elif (isinstance(index, tuple) and len(index) == 2 and index[1] in[index] = value else: assert isinstance(index, slice) or len(index) > 0[index] = value def __delitem__(self, index): return NotImplemented def __iter__(self): return __next__ = next @property def chromosome(self): return['chromosome'] @property def start(self): return['start'] @property def end(self): return['end'] @property def sample_id(self): return self.meta.get('sample_id') # Traversal
[docs] def autosomes(self, also=()): """Select chromosomes w/ integer names, ignoring any 'chr' prefixes.""" with warnings.catch_warnings(): # NB: We're not using the deprecated part of this pandas method warnings.simplefilter("ignore", UserWarning) is_auto = self.chromosome.str.match(r"(chr)?\d+$", as_indexer=True, na=False) if not is_auto.any(): # The autosomes, if any, are not named with plain integers return self if also: if isinstance(also, basestring): also = [also] for a_chrom in also: is_auto |= (self.chromosome == a_chrom) return self[is_auto]
[docs] def by_chromosome(self): """Iterate over bins grouped by chromosome name.""" for chrom, subtable in"chromosome", sort=False): yield chrom, self.as_dataframe(subtable)
[docs] def by_ranges(self, other, mode='inner', keep_empty=True): """Group rows by another GenomicArray's bin coordinate ranges. Returns an iterable of (bin, GenomicArray of overlapping rows)). Usually used for grouping probes or SNVs by CNV segments. `mode` determines what to do with bins that overlap a boundary of the selection. Values are: - ``inner``: Drop the bins on the selection boundary, don't emit them. - ``outer``: Keep/emit those bins as they are. - ``trim``: Emit those bins but alter their boundaries to match the selection; the bin start or end position is replaced with the selection boundary position. [default] Bins in this array that fall outside the other array's bins are skipped. """ chrom_lookup = dict(self.by_chromosome()) for chrom, bin_rows in other.by_chromosome(): if chrom in chrom_lookup: subranges = chrom_lookup[chrom]._iter_ranges( None, bin_rows['start'], bin_rows['end'], mode) for bin_row, subrange in zip(bin_rows, subranges): yield bin_row, subrange else: if keep_empty: for bin_row in bin_rows: yield bin_row, self.as_rows([])
[docs] def coords(self, also=()): """Iterate over plain coordinates of each bin: chromosome, start, end. With `also`, also include those columns. Example, yielding rows in BED format: >>> probes.coords(also=["name", "strand"]) """ cols = list(GenomicArray._required_columns) if also: if isinstance(also, basestring): cols.append(also) else: cols.extend(also) coordframe =[:, cols] return coordframe.itertuples(index=False)
[docs] def labels(self): return, axis=1)
[docs] def in_range(self, chrom=None, start=None, end=None, mode='inner'): """Get the GenomicArray portion within the given genomic range. `mode` works as in `by_ranges`: ``outer`` includes bins straddling the range boundaries, ``trim`` additionally alters the straddling bins' endpoints to match the range boundaries, and ``inner`` excludes those bins. """ if isinstance(start, (int, float, np.float64)): start = [int(start)] if isinstance(end, (int, float, np.float64)): end = [int(end)] results = self._iter_ranges(chrom, start, end, mode) return next(results)
[docs] def in_ranges(self, chrom=None, starts=None, ends=None, mode='inner'): """Get the GenomicArray portion within the specified ranges. Same as `in_ranges` but the `starts` and `ends` are arrays of equal length, and the output concatenates all the selected bins. """ return self.concat(self._iter_ranges(chrom, starts, ends, mode))
def _iter_ranges(self, chrom, starts, ends, mode): """Iterate through sub-ranges.""" assert mode in ('inner', 'outer', 'trim') if chrom: assert isinstance(chrom, basestring) # ENH: accept array? try: table =[['chromosome'] == chrom] except KeyError: raise KeyError("Chromosome %s is not in this probe set" % chrom) else: # Unsafe, but faster if we've already subsetted by chromosome table = # Edge cases if not len(table): yield self.as_rows([]) raise StopIteration if starts is None and ends is None: yield self.as_dataframe(table) raise StopIteration if starts is not None and len(starts): if mode == 'inner': # Only rows entirely after the start point start_idxs = table.start.searchsorted(starts) else: # Include all rows overlapping the start point start_idxs = table.end.searchsorted(starts, 'right') else: starts = np.zeros(len(ends) if ends is not None else 1, dtype=np.int_) start_idxs = starts.copy() if ends is not None and len(ends): if mode == 'inner': end_idxs = table.end.searchsorted(ends, 'right') else: end_idxs = table.start.searchsorted(ends) else: end_idxs = np.repeat(len(table), len(starts)) ends = [None] * len(starts) for start_idx, start_val, end_idx, end_val in zip(start_idxs, starts, end_idxs, ends): subtable = table[start_idx:end_idx] if mode == 'trim': subtable = subtable.copy() # Update 5' endpoints to the boundary if start_val: subtable.start = subtable.start.clip_lower(start_val) # Update 3' endpoints to the boundary if end_val: subtable.end = subtable.end.clip_upper(end_val) yield self.as_dataframe(subtable)
[docs] def match_to_bins(self, other, key, default=0.0, fill=False, summary_func=np.median): """Take values of the other array at each of this array's bins. Assign `default` to indices that fall outside the other array's bins, or chromosomes that appear in `self` but not `other`. Return an array of the `key` column values in `other` corresponding to this array's bin locations, the same length as this array. """ def rows2value(rows): if len(rows) == 0: return default elif len(rows) == 1: return rows[0, key] else: return summary_func(rows[key]) all_out_vals = [rows2value(other_rows) for _bin, other_rows in other.by_ranges(self, mode='outer', keep_empty=True)] return np.asarray(all_out_vals)
# Modification
[docs] def add(self, other): """Combine this array's data with another GenomicArray (in-place). Any optional columns must match between both arrays. """ if not isinstance(other, self.__class__): raise ValueError("Argument (type %s) is not a %s instance" % (type(other), self.__class__)) if not len( return self.copy() = pd.concat([,]) self.sort()
[docs] def concat(self, others): """Concatenate several GenomicArrays, keeping this array's metadata. This array's data table is not implicitly included in the result. """ result = self.as_dataframe(pd.concat([ for otr in others])) result.sort() return result
[docs] def copy(self): """Create an independent copy of this object.""" return self.as_dataframe(
[docs] def add_columns(self, **columns): """Create a new CNA, adding the specified extra columns to this CNA.""" # return self.as_dataframe(**columns)) result = self.copy() for key, values in columns.iteritems(): result[key] = values return result
[docs] def keep_columns(self, columns): """Extract a subset of columns, reusing this instance's metadata.""" return self.__class__([:, columns], self.meta.copy())
[docs] def drop_extra_columns(self): """Remove any optional columns from this GenomicArray. Returns a new copy with only the core columns retained: log2 value, chromosome, start, end, bin name. """ table =[:, self._required_columns] return self.as_dataframe(table)
[docs] def select(self, selector=None, **kwargs): """Take a subset of rows where the given condition is true. Arguments can be a function (lambda expression) returning a bool, which will be used to select True rows, and/or keyword arguments like gene="Background" or chromosome="chr7", which will select rows where the keyed field equals the specified value. """ table = if selector is not None: table = table[table.apply(selector, axis=1)] for key, val in kwargs.items(): assert key in self table = table[table[key] == val] return self.as_dataframe(table)
[docs] def shuffle(self): """Randomize the order of bins in this array (in-place).""" np.random.seed(0xA5EED) # For reproducible results order = np.arange(len( np.random.shuffle(order) =[order] return order
[docs] def sort(self): """Sort this array's bins in-place, with smart chromosome ordering.""" table = table['SORT_KEY'] = self.chromosome.apply(core.sorter_chrom) table.sort_values(by=['SORT_KEY', 'start'], inplace=True) del table['SORT_KEY'] = table.reset_index(drop=True)
[docs] def sort_columns(self): """Sort this array's columns in-place, per class definition.""" extra_cols = [] for col in if col not in self._required_columns: extra_cols.append(col) sorted_colnames = list(self._required_columns) + sorted(extra_cols) assert len(sorted_colnames) == len( =
# I/O @classmethod
[docs] def read(cls, infile, sample_id=None): if sample_id is None: if isinstance(infile, basestring): sample_id = core.fbase(infile) else: sample_id = '<unknown>' # Create a multi-index of genomic coordinates (like GRanges) try: table = pd.read_table(infile, na_filter=False, dtype={'chromosome': 'string'}, # index_col=['chromosome', 'start'] ) except ValueError: # File is blank/empty, most likely"Blank file %s", infile) table = cls._make_blank() # XXX Pending pandas 0.17: # table['chromosome'] = pd.Categorical(table['chromosome'], # table.chromosome.drop_duplicates(), # ordered=True) # table.set_index(['chromosome', 'start'], inplace=True) return cls(table, {"sample_id": sample_id})
[docs] def write(self, outfile=None): """Write the wrapped data table to a file or handle in tabular format. The format is BED-like, but with a header row included and with arbitrary extra columns. To combine multiple samples in one file and/or convert to another format, see the 'export' subcommand. """ with ngfrills.safe_write(outfile or sys.stdout) as handle:, index=False, sep='\t', float_format='%.6g')