#!/usr/bin/env python """Provides Profile and ProfileError object and CharMeaningProfile Owner: Sandra Smit (Sandra Smit) """ from __future__ import division #SUPPORT2425 #from __future__ import with_statement from string import maketrans, translate from numpy import array, sum, transpose, reshape, ones, zeros,\ take, float64, ravel, nonzero, log, put, concatenate, argmax, cumsum,\ sort, argsort, searchsorted, logical_and, asarray, uint8, add, subtract,\ multiply, divide, newaxis, alltrue, max, all, isfinite #from numpy.oldnumeric import sum from numpy.random import random from cogent.util.array import euclidean_distance, row_degeneracy,\ column_degeneracy, row_uncertainty, column_uncertainty, safe_log from cogent.format.table import formattedCells ##SUPPORT2425 import numpy #from cogent.util.unit_test import numpy_err __author__ = "Sandra Smit" __copyright__ = "Copyright 2007-2009, The Cogent Project" __credits__ = ["Sandra Smit", "Gavin Huttley", "Rob Knight" "Peter Maxwell"] __license__ = "GPL" __version__ = "1.4.1" __maintainer__ = "Sandra Smit" __email__ = "sandra.smit@colorado.edu" __status__ = "Production" class ProfileError(Exception): """Error raised for exceptions occuring in the Profile object""" pass class Profile(object): """Profile class """ def __init__(self, Data, Alphabet, CharOrder=None): """Initializes a new Profile object. Data: numpy 2D array with the Profile data in it. Specifically, each row of the array corresponds to a position in the Alignment. Each column of the array corresponds to a character in the Alphabet. Alphabet: an Alphabet object or anything that can act as a list of characters CharOrder: optional list of characters to which the columns in the Data correspond. """ self.Data = Data self.Alphabet = Alphabet if CharOrder is None: self.CharOrder = list(self.Alphabet) else: self.CharOrder = CharOrder #the translation table is needed for making consensus sequences, #but will fail if the alphabet isn't made of chars (in which case, #we'll just skip the translation table, and certain downstream #operations may fail). try: self._translation_table = self._make_translation_table() except: pass def __str__(self): """Returns string representation of self.Data""" return str(self.Data) def _make_translation_table(self): """Makes a translation tables between the CharOrder and indices """ indices = ''.join(map(chr, range(len(self.CharOrder)))) chars = ''.join(map(str,self.CharOrder)) return maketrans(chars, indices) def hasValidData(self, err=1e-16): """Returns True if all rows in self.Data add up to one err -- float, maximum deviation from 1 allowed, default is 1e-16 Rounding errors might occur, so a small deviation from 1 is allowed. The default tolerance is 1e-16. """ obs_sums = sums = sum(self.Data,1) lower_bound = ones(len(self.Data)) - err upper_bound = ones(len(self.Data)) + err tfs = sum(self.Data,1) == ones(len(self.Data)) if (lower_bound <= obs_sums).all() and (obs_sums <= upper_bound).all(): return True return False def hasValidAttributes(self): """Checks Alphabet, CharOrder, and size of self.Data""" if not reduce(logical_and, [c in self.Alphabet\ for c in self.CharOrder]): return False elif self.Data.shape[1] != len(self.CharOrder): return False return True def isValid(self): """Check whether everything in the Profile is valid""" vd = self.hasValidData() va = self.hasValidAttributes() return vd and va def dataAt(self, pos, character=None): """Return data for a certain position (row!) and character (column) pos -- int, position (row) in the profile character -- str, character from the CharacterOrder If character is None, all data for the position is returned. """ if not 0 <= pos < len(self.Data): raise ProfileError(\ "Position %s is not present in the profile"%(pos)) if character is None: return self.Data[pos,:] else: if character not in self.CharOrder: raise ProfileError(\ "Character %s is not present in the profile's CharacterOrder"\ %(character)) return self.Data[pos, self.CharOrder.index(character)] def copy(self): """Returns a copy of the Profile object WARNING: the data, alphabet and the character order are the same object in the original and the copy. This means you can rebind the attributes, but modifying them will change them in both the original and the copy. """ return self.__class__(self.Data, self.Alphabet, self.CharOrder) def normalizePositions(self): """Normalizes the data by position (the rows!) to one It does not make sense to normalize anything with negative numbers in there. However, the method does NOT check for that, because it would slow down the calculations too much. It will work, but you might get very unexpected results. The method will raise an error when one or more rows add up to one. It checks explicitly for that to avoid OverflowErrors, ZeroDivisionErrors, and infinities in the results. WARNING: this method works in place with respect to the Profile object, not with respect to the Data attribute. Normalization rebinds self.Data to a new array. """ row_sums = sum(self.Data,1) if (row_sums == 0).any(): zero_indices = nonzero(row_sums==0)[0].tolist() raise ProfileError,\ "Can't normalize profile, rows at indices %s add up to zero"\ %(zero_indices) else: self.Data = self.Data/row_sums[:,newaxis] def normalizeSequences(self): """Normalized the data by sequences (the columns) to one It does not make sense to normalize anything with negative numbers in there. However, the method does NOT check for that, because it would slow down the calculations too much. It will work, but you might get very unexpected results. The method will raise an error when one or more columns add up to one. It checks explicitly for that to avoid OverflowErrors, ZeroDivisionErrors, and infinities in the results. WARNING: this method works in place with respect to the Profile object, not with respect to the Data attribute. Normalization rebinds self.Data to a new array. """ col_sums = sum(self.Data, axis=0) if (col_sums == 0).any(): zero_indices = nonzero(col_sums==0)[0].tolist() raise ProfileError,\ "Can't normalize profile, columns at indices %s add up to zero"\ %(zero_indices) else: self.Data = self.Data/col_sums def prettyPrint(self, include_header=False, transpose_data=False,\ column_limit=None, col_sep='\t'): """Returns a string method of the data and character order. include_header: include charcter order or not transpose_data: data as is (rows are positions) or transposed (rows are characters) to line it up with an alignment column_limit = int, maximum number of columns displayed col_sep = string, column separator """ h = self.CharOrder d = self.Data if column_limit is None: max_col_idx = d.shape[1] else: max_col_idx = column_limit if include_header and not transpose_data: r = [h]+d.tolist() elif include_header and transpose_data: r =[[x] + y for x,y in zip(h,transpose(d).tolist())] elif transpose_data: r = transpose(d).tolist() else: r = d.tolist() # resize the result based on the column limit if column_limit is not None: r = [row[:column_limit] for row in r] # nicely format the table content, discard the header (already included) if r: new_header, formatted_res = formattedCells(r) else: formatted_res = r return '\n'.join([col_sep.join(map(str,i)) for i in formatted_res]) def reduce(self,other,op=add,normalize_input=True,normalize_output=True): """Reduces two profiles with some operator and returns a new Profile other: Profile object op: operator (e.g. add, subtract, multiply, divide) normalize_input: whether the input profiles will be normalized before collapsing. The default is True. normalize_output: whether the output profile will be normalized. The default is True This function is intented for use on normalized profiles. For safety it'll try to normalize the data before collapsing them. If you do not normalize your data and set normalize_input to False, you might get unexpected results. It does check whether self.Data and other.Data have the same shape It does not check whether self and other have the same CharOrder. The resulting Profile gets the alphabet and char order from self. """ if self.Data.shape != other.Data.shape: raise ProfileError,\ "Cannot collapse profiles of different size: %s, %s"\ %(self.Data.shape,other.Data.shape) if normalize_input: self.normalizePositions() other.normalizePositions() try: ##SUPPORT2425 ori_err = numpy.geterr() numpy.seterr(divide='raise') try: new_data = op(self.Data, other.Data) finally: numpy.seterr(**ori_err) #with numpy_err(divide='raise'): #new_data = op(self.Data, other.Data) except (OverflowError, ZeroDivisionError, FloatingPointError): raise ProfileError, "Can't do operation on input profiles" result = Profile(new_data, self.Alphabet, self.CharOrder) if normalize_output: result.normalizePositions() return result def __add__(self,other): """Binary + operator: adds two profiles element-wise. Input and output are NOT normalized. """ return self.reduce(other, op=add, normalize_input=False,\ normalize_output=False) def __sub__(self,other): """Binary - operator: subtracts two profiles element-wise Input and output are NOT normalized. """ return self.reduce(other, op=subtract, normalize_input=False,\ normalize_output=False) def __mul__(self,other): """* operator: multiplies two profiles element-wise Input and output are NOT normalized. """ return self.reduce(other, op=multiply, normalize_input=False,\ normalize_output=False) def __div__(self,other): """/ operator for old-style division: divides 2 profiles element-wise. Used when __future__.divsion not imported Input and output are NOT normalized. """ return self.reduce(other, op=divide, normalize_input=False,\ normalize_output=False) def __truediv__(self,other): """/ operator for new-style division: divides 2 profiles element-wise. Used when __future__.division is in action. Input and output are NOT normalized. """ return self.reduce(other, op=divide, normalize_input=False,\ normalize_output=False) def distance(self, other, method=euclidean_distance): """Returns the distance between two profiles other: Profile object method: function used to calculated the distance between two arrays. WARNING: In principle works only on profiles of the same size. However, when one of the two profiles is 1D (which shouldn't happen) and can be aligned with the other profile the distance is still calculated and may give unexpected results. """ try: return method(self.Data, other.Data) except ValueError: #frames not aligned raise ProfileError,\ "Profiles have different size (and are not aligned): %s %s"\ %(self.Data.shape,other.Data.shape) def toOddsMatrix(self, symbol_freqs=None): """Returns the OddsMatrix of a profile as a new Profile. symbol_freqs: per character array of background frequencies e.g. [.25,.25,.25,.25] for equal frequencies for each of the four bases. If no symbol frequencies are provided, all symbols will get equal freqs. The length of symbol freqs should match the number of columns in the profile! If symbol freqs contains a zero entry, a ProfileError is raised. This is done to prevent either a ZeroDivisionError (raised when zero is an int) or 'inf' in the resulting matrix (which happens when zero is a float). """ pl = self.Data.shape[1] #profile length #if symbol_freqs is None, create an array with equal frequencies if symbol_freqs is None: symbol_freqs = ones(pl)/pl else: symbol_freqs = array(symbol_freqs) #raise error when symbol_freqs has wrong length if len(symbol_freqs) != pl: raise ProfileError,\ "Length of symbol freqs should be %s, but is %s"\ %(pl,len(symbol_freqs)) #raise error when symbol freqs contains zero (to prevent #ZeroDivisionError or 'inf' in the resulting matrix) if sum(symbol_freqs != 0, 0) != len(symbol_freqs): raise ProfileError,\ "Symbol frequency is not allowed to be zero: %s"\ %(symbol_freqs) #calculate the OddsMatrix log_odds = self.Data/symbol_freqs return Profile(log_odds, self.Alphabet, self.CharOrder) def toLogOddsMatrix(self, symbol_freqs=None): """Returns the LogOddsMatrix of a profile as a new Profile/ symbol_freqs: per character array of background frequencies e.g. [.25,.25,.25,.25] for equal frequencies for each of the four bases. See toOddsMatrix for more information. """ odds = self.toOddsMatrix(symbol_freqs) log_odds = safe_log(odds.Data) return Profile(log_odds, self.Alphabet, self.CharOrder) def _score_indices(self, seq_indices, offset=0): """Returns score of the profile for each slice of the seq_indices seq_indices: translation of sequence into indices that match the characters in the CharOrder of the profile offset: where to start the matching procedure This function doesn't do any input validation. That is done in 'score' See method 'score' for more information. """ data = self.Data pl = len(data) #profile length (number of positions) sl = len(seq_indices) r = range(pl) #fixed range result = [] for starting_pos in range(offset, len(seq_indices)-pl+1): slice = seq_indices[starting_pos:starting_pos+pl] result.append(sum(array([data[i] for i in zip(r,slice)]), axis=0)) return array(result) def _score_profile(self, profile, offset=0): """Returns score of the profile against the input_profile. profile: Profile of a sequence or alignment that has to be scored offset: where to start the matching procedure This function doesn't do any input validation. That is done in 'score' See method 'score' for more information. """ data = self.Data self_l = len(data) #profile length other_l = len(profile.Data) #other profile length result = [] for start in range(offset,other_l-self_l+1): stop = start + self_l slice = profile.Data[start:stop,:] result.append(sum(self.Data*slice)) return array(result) def score(self, input_data, offset=0): """Returns a score of the profile against input_data (Profile or Seq). seq: Profile or Sequence object (or string) offset: starting index for searching in seq/profile Returns the score of the profile against all possible subsequences/ subprofiles of the input_data. This method determines how well a profile fits at different places in the sequence. This is very useful when the profile is a motif and you want to find the position in the sequence that best matches the profile/motif. Sequence Example: ================= T C A G 0 .2 .4 .4 0 1 .1 0 .9 0 2 .1 .2 .3 .4 Sequence: TCAAGT pos 0: TCA -> 0.5 pos 1: CAA -> 1.6 pos 2: AAG -> 1.7 pos 3: AGT -> 0.5 So the subsequence starting at index 2 in the sequence has the best match with the motif Profile Example: ================ Profile: same as above Profile to score: T C A G 0 1 0 0 0 1 0 1 0 0 2 0 0 .5 .5 3 0 0 0 1 4 .25 .25 .25 .25 pos 0: rows 0,1,2 -> 0.55 pos 1: rows 1,2,3 -> 1.25 pos 2: rows 2,3,4 -> 0.45 """ #set up some local variables data = self.Data pl = len(data) #profile length is_profile = False #raise error if profile is empty if not data.any(): raise ProfileError,"Can't score an empty profile" #figure out what the input_data type is if isinstance(input_data,Profile): is_profile = True to_score_length = len(input_data.Data) #raise error if CharOrders don't match if self.CharOrder != input_data.CharOrder: raise ProfileError, "Profiles must have same character order" else: #assumes it get a sequence to_score_length = len(input_data) #Profile should fit at least once in the sequence/profile_to_score if to_score_length < pl: raise ProfileError,\ "Sequence or Profile to score should be at least %s "%(pl)+\ "characters long, but is %s."%(to_score_length) #offset should be valid if not offset <= (to_score_length - pl): raise ProfileError, "Offset must be <= %s, but is %s"\ %((to_score_length-pl), offset) #call the apropriate scoring function if is_profile: return self._score_profile(input_data, offset) else: #translate seq to indices if hasattr(self, '_translation_table'): seq_indices = array(map(ord,translate(str(input_data),\ self._translation_table))) else: #need to figure out where each item is in the charorder idx = self.CharOrder.index seq_indices = array(map(idx, input_data)) #raise error if some sequence characters are not in the CharOrder if (seq_indices > len(self.CharOrder)).any(): raise ProfileError,\ "Sequence contains characters that are not in the "+\ "CharOrder" #now the profile is scored against the list of indices return self._score_indices(seq_indices,offset) def rowUncertainty(self): """Returns the uncertainty (Shannon's entropy) for each row in profile Entropy is returned in BITS (not in NATS). """ if not self.Data.any(): return array([]) try: return row_uncertainty(self.Data) except ValueError: raise ProfileError,\ "Profile has to be two dimensional to calculate rowUncertainty" def columnUncertainty(self): """Returns uncertainty (Shannon's entropy) for each column in profile Uncertainty is returned in BITS (not in NATS). """ if not self.Data.any(): return array([]) try: return column_uncertainty(self.Data) except ValueError: raise ProfileError,\ "Profile has to be two dimensional to calculate columnUncertainty" def rowDegeneracy(self, cutoff=0.5): """Returns how many chars are needed to cover the cutoff value. cutoff: value that should be covered in each row For example: pos 0: [.1,.2,.3,.4] char order=TCAG. If cutoff=0.75 -> degeneracy = 3 (degenearate char for CAG) If cutoff=0.25 -> degeneracy = 1 (G alone covers this cutoff) If cutoff=0.5 -> degeneracy = 2 (degenerate char for AG) If the cutoff value is not reached in the row, the returned value will be clipped to the length of the character order (=the number of columns in the Profile). """ try: return row_degeneracy(self.Data,cutoff) except ValueError: raise ProfileError,\ "Profile has to be two dimensional to calculate rowDegeneracy" def columnDegeneracy(self, cutoff=0.5): """Returns how many chars are neede to cover the cutoff value See rowDegeneracy for more information. """ try: return column_degeneracy(self.Data,cutoff) except ValueError: raise ProfileError,\ "Profile has to be two dimensional to calculate columnDegeneracy" def rowMax(self): """Returns ara containing most frequent element in each row of the profile.""" return max(self.Data, 1) def toConsensus(self, cutoff=None, fully_degenerate=False,\ include_all=False): """Returns the consensus sequence from a profile. cutoff: cutoff value, determines how much should be covered in a position (row) of the profile. Example: pos 0 [.2,.1,.3,.4] (CharOrder: TCAG). To cover .65 (=cutoff) we need two characters: A and G, which results in the degenerate character R. fully_degenerate: determines whether the fully degenerate character is returned at a position. For the example above an 'N' would be returned. inlcude_all: all possibilities are included in the degenerate character. Example: row = UCAG = [.1,.3,.3,.3] cutoff = .4, consensus = 'V' (even though only 2 chars would be enough to reach the cutoff value). The Alphabet of the Profile should implement degenerateFromSequence. Note that cutoff has priority over fully_degenerate. In other words, if you specify a cutoff value and set fully_degenerate to true, the calculation will be done with the cutoff value. If nothing gets passed in, the maximum argument is chosen. In the first example above G will be returned. """ #set up some local variables co = array(self.CharOrder, 'c') alpha = self.Alphabet data = self.Data #determine the action. Cutoff takes priority over fully_degenerate if cutoff: result = [] degen = self.rowDegeneracy(cutoff) sorted = argsort(data) if include_all: #if include_all include all possiblilities in the degen char for row_idx, (num_to_keep, row) in enumerate(zip(degen,sorted)): to_take = [item for item in row[-num_to_keep:]\ if item in nonzero(data[row_idx])[0]] +\ [item for item in nonzero(data[row_idx] ==\ data[row_idx,row[-num_to_keep]])[0] if item in\ nonzero(data[row_idx])[0]] result.append(alpha.degenerateFromSequence(\ map(str,take(co, to_take, axis=0)))) else: for row_idx, (num_to_keep, row) in enumerate(zip(degen,sorted)): result.append(alpha.degenerateFromSequence(\ map(str,take(co, [item for item in row[-num_to_keep:]\ if item in nonzero(data[row_idx])[0]])))) elif not fully_degenerate: result = take(co, argmax(self.Data, axis=-1), axis=0) else: result = [] for row in self.Data: result.append(alpha.degenerateFromSequence(\ map(str,take(co, nonzero(row)[0], axis=0)))) return ''.join(map(str,result)) def randomIndices(self, force_accumulate=False, random_f = random): """Returns random indices matching current probability matrix. Stores cumulative sum (sort of) of probability matrix in self._accumulated; Use force_accumulate to reset if you change the matrix in place (which you shouldn't do anyway). The returned indices correspond to the characters in the CharOrder of the Profile. """ if force_accumulate or not hasattr(self, '_accumulated'): self._accumulated = cumsum(self.Data, 1) choices = random_f(len(self.Data)) return array([searchsorted(v, c) for v, c in\ zip(self._accumulated, choices)]) def randomSequence(self, force_accumulate=False, random_f = random): """Returns random sequence matching current probability matrix. Stores cumulative sum (sort of) of probability matrix in self._accumulated; Use force_accumulate to reset if you change the matrix in place (which you shouldn't do anyway). """ co = self.CharOrder random_indices = self.randomIndices(force_accumulate,random_f) return ''.join(map(str,take(co,random_indices))) def CharMeaningProfile(alphabet, char_order=None, split_degenerates=False): """Returns a Profile with the meaning of each character in the alphabet alphabet: Alphabet object (should have 'Degenerates'if split_degenerates is set to True) char_order: string indicating the order of the characters in the profile split_degenerates: whether the meaning of degenerate symbols in the alphabet should be split up among the characters in the char order, or ignored. The returned profile has 255 rows (one for each ascii character) and one column for each character in the character order. The profile specifies the meaning of each character in the alphabet. Chars in the character order will count as a full character by themselves, degenerate characters might split their 'value' over several other charcters in the character order. Splitting up degenerates: only degenerate characters of which the full set of symbols it maps onto are in the character order are split up, others are ignored. E.g. in the DnaAlphabet, if the char order is TACG, ? (which maps to TCAG-) wouldn't be split up, 'R' (which maps to 'AG') would. Any degenerate characters IN the character order will NOT be split up. It doesn't make sense to split up a character that is in the char order because it would create an empty column in the profile, so it might as well be left out alltogether. Example 1: Alphabet = DnaAlphabet Character order = "TCAG" Split degenerates = False All the nonzero rows in the resulting profile are: 65: [0,0,1,0] (A) 67: [0,1,0,0] (C) 71: [0,0,0,1] (G) 84: [1,0,0,0] (T) All other rows will be [0,0,0,0]. Example 2: Alphabet = DnaAlphabet Character order = "AGN" Split degenerates = True All the nonzero rows in the resulting profile are: 65: [1,0,0] (A) 71: [0,1,0] (G) 78: [0,0,1] (N) 82: [.5,.5,0] (R) All other rows will be [0,0,0]. Errors are raised when the character order is empty or when there's a character in the character order that is not in the alphabet. """ if not char_order: #both testing for None and for empty string char_order = list(alphabet) char_order = array(char_order, 'c') lc = len(char_order) #length char_order #initialize the profile. 255 rows (one for each ascii char), one column #for each character in the character order result = zeros([255,lc],float64) if split_degenerates: degen = alphabet.Degenerates for degen_char in degen: #if all characters that the degenerate character maps onto are #in the character order, split its value up according to the #alphabet curr_degens = degen[degen_char] if all(map(char_order.__contains__, curr_degens)): contains = map(curr_degens.__contains__, char_order) result[ord(degen_char)] = \ array(contains, float)/len(curr_degens) #for each character in the character order, make an entry of ones and #zeros, matching the character order for c in char_order: c = str(c) if c not in alphabet: raise ValueError, "Found character in the character order "+\ "that is not in the specified alphabet: %s"%(c) result[ord(c)] = array(c*lc, 'c') == char_order return Profile(Data=result,Alphabet=alphabet,CharOrder=char_order)