#!/usr/bin/env python """Code for handling multiple sequence alignments. In particular: - SequenceCollection handles both aligned and unaligned sequences. - Alignment and its subclasses handle multiple sequence alignments, storing the raw sequences and a gap map. Useful for very long alignments, e.g. genomics data. - DenseAlignment and its subclasses handle multiple sequence alignments as arrays of characters. Especially useful for short alignments that contain many sequences. WARNING: The various alignment objects try to guess the input type from the input, but this behavior has a few quirks. In particular, if the input is a sequence of two-item sequences (e.g. a list of two-character strings), each sequence will be unpacked and the first item will be used as the label, the second as the sequence. For example, Alignment(['AA','CC','AA']) produces an alignment of three 1-character strings labeled A, C and A respectively. The reason for this is that the common case is that you have passed in a stream of two-item label, sequence pairs. However, this can cause confusion when testing. """ from __future__ import division from types import GeneratorType from cogent.core.annotation import Map, _Annotatable import cogent #will use to get at cogent.parse.fasta.MinimalFastaParser, #which is a circular import otherwise. from cogent.format.alignment import save_to_filename from cogent.core.info import Info as InfoClass from cogent.core.sequence import frac_same, ModelSequence from cogent.maths.stats.util import Freqs from cogent.format.fasta import fasta_from_alignment from cogent.format.phylip import phylip_from_alignment from cogent.format.nexus import nexus_from_alignment from cogent.parse.gff import GffParser, parse_attributes from numpy import nonzero, array, logical_or, logical_and, logical_not, \ transpose, arange, zeros, ones, take, put, uint8, ndarray from numpy.random import randint, permutation from cogent.util.dict2d import Dict2D import logging LOG = logging.getLogger('cogent.data') from copy import copy from cogent.core.profile import Profile __author__ = "Peter Maxwell and Rob Knight" __copyright__ = "Copyright 2007-2009, The Cogent Project" __credits__ = ["Peter Maxwell", "Rob Knight", "Gavin Huttley", "Jeremy Widmann", "Catherine Lozupone", "Matthew Wakefield", "Micah Hamady", "Daniel McDonald"] __license__ = "GPL" __version__ = "1.4.1" __maintainer__ = "Rob Knight" __email__ = "rob@spot.colorado.edu" __status__ = "Production" class DataError(Exception): pass eps = 1e-6 #small number: 1-eps is almost 1, and is used for things like the #default number of gaps to allow in a column. def assign_sequential_names(ignored, num_seqs, base_name='seq', start_at=0): """Returns list of num_seqs sequential, unique names. First argument is ignored; expect this to be set as a class attribute. """ return ['%s_%s' % (base_name,i) for i in range(start_at,start_at+num_seqs)] class SeqLabeler(object): """Allows flexible seq labeling in toFasta().""" def __init__(self, aln, label_f=assign_sequential_names, **kwargs): """Initializes a new seq labeler.""" self._aln = aln self._label_f = label_f self._map = dict(zip(aln.Names, label_f(len(aln.Names, **kwargs)))) def __call__(self, s): """Returns seq name from seq id""" return self._map[s.Name] def coerce_to_string(s): """Converts an arbitrary sequence into a string.""" if isinstance(s, str): #if it's a string, OK as is return s if isinstance(s, Aligned): #if it's an Aligned object, convert to string return str(s) curr = str(s) #if its string is the same length, return that if len(curr) == len(s): return curr try: return ''.join(s) #assume it's a seq of chars except(TypeError, ValueError): return ''.join(map(str, s)) #general case (slow, might not be correct) def seqs_from_array(a, Alphabet=None): """SequenceCollection from array of pos x seq: names are integers. This is an InputHandler for SequenceCollection. It converts an arbitrary array of numbers into Sequence objects using seq_constructor, and leaves the sequences unlabeled. """ return list(transpose(a)), None def seqs_from_model_seqs(seqs, Alphabet=None): """Alignment from ModelSequence objects: seqs -> array, names from seqs. This is an InputHandler for SequenceCollection. It converts a list of Sequence objects with _data and Name properties into a SequenceCollection that uses those sequences. """ return seqs, [s.Name for s in seqs] def seqs_from_generic(seqs, Alphabet=None): """SequenceCollection from generic seq x pos data: seq of seqs of chars. This is an InputHandler for SequenceCollection. It converts a generic list (each item in the list will be mapped onto an object using seq_constructor and assigns sequential integers (0-based) as names. """ names = [] for s in seqs: if hasattr(s, 'Name'): names.append(s.Name) else: names.append(None) return seqs, names def seqs_from_fasta(seqs, Alphabet=None): """SequenceCollection from FASTA-format string or lines. This is an InputHandler for SequenceCollection. It converts a FASTA-format string or collection of lines into a SequenceCollection object, preserving order.. """ if isinstance(seqs, str): seqs = seqs.splitlines() names, seqs = zip(*list(cogent.parse.fasta.MinimalFastaParser(seqs))) return list(seqs), list(names) def seqs_from_dict(seqs, Alphabet=None): """SequenceCollection from dict of {label:seq_as_str}. This is an InputHandler for SequenceCollection. It converts a dict in which the keys are the names and the values are the sequences (sequence only, no whitespace or other formatting) into a SequenceCollection. Because the dict doesn't preserve order, the result will not necessarily be in alphabetical order.""" names, seqs = map(list, zip(*seqs.items())) return seqs, names def seqs_from_kv_pairs(seqs, Alphabet=None): """SequenceCollection from list of (key, val) pairs. This is an InputHandler for SequenceCollection. It converts a dict in which the keys are the names and the values are the sequences (sequence only, no whitespace or other formatting) into a SequenceCollection. Because the dict doesn't preserve order, the result will be in arbitrary order.""" names, seqs = map(list, zip(*seqs)) return seqs, names def seqs_from_aln(seqs, Alphabet=None): """SequenceCollection from existing SequenceCollection object: copies data. This is relatively inefficient: you should really use the copy() method instead, which duplicates the internal data structures. """ return seqs.Seqs, seqs.Names def seqs_from_empty(obj, *args, **kwargs): """SequenceCollection from empty data: raise exception.""" raise ValueError, "Cannot create empty SequenceCollection." class SequenceCollection(object): """Base class for Alignment, but also just stores unaligned seqs. Handles shared functionality: detecting the input type, writing out the sequences as different formats, translating the sequences, chopping off stop codons, looking up sequences by name, etc. A SequenceCollection must support: - input handlers for different data types - SeqData: behaves like list of lists of chars, holds seq data - Seqs: behaves like list of Sequence objects, iterable in name order - Names: behaves like list of names for the sequence objects - NamedSeqs: behaves like dict of {name:seq} - MolType: specifies what kind of sequences are in the collection """ InputHandlers = { 'array': seqs_from_array, 'model_seqs': seqs_from_model_seqs, 'generic': seqs_from_generic, 'fasta': seqs_from_fasta, 'collection': seqs_from_aln, 'aln': seqs_from_aln, 'dense_aln': seqs_from_aln, 'dict': seqs_from_dict, 'empty': seqs_from_empty, 'kv_pairs':seqs_from_kv_pairs, } IsArray = set(['array', 'model_seqs']) DefaultNameFunction = assign_sequential_names def __init__(self, data, Names=None, Alphabet=None, MolType=None, \ Name=None, Info=None, conversion_f=None, is_array=False, \ force_same_data=False, \ remove_duplicate_names=False, label_to_name=None, suppress_named_seqs=False): """Initialize self with data and optionally Info. We are always going to convert to characters, so Sequence objects in the collection will lose additional special attributes they have. This is somewhat inefficient, so it might be worth revisiting this decision later. The handling of sequence names requires special attention. Depending on the input data, we might get the names from the sequences themselves, or we might add them from Names that are passed in. However, the Names attribute controls the order that we examine the sequences in, so if it is passed in it should override the order that we got from the input data (e.g. you might pass in unlabeled sequences with the names ['b','a'], so that you want the first sequence to be called 'b' and the second to be called 'a', or you might pass in labeled sequences, e.g. as a dict, and the names ['b','a'], indicating that you want the sequence called b to be first and the sequence called a to be second despite the fact that they are in arbitrary order in the original input. In this second situation, it is imortant that the sequences not be relabeled. This is handled as followed. If the sequences are passed in using a method that does not carry the names with it, the Names that are passed in will be handed out to successive sequences. If the sequences are passed in using a method that does carry the names with it, the Names that are passed in will be used to order the sequences, but they will not be relabeled. Note that if you're passing in a data type that is already labeled (e.g. a list of Sequence objects) you _must_ have unique names beforehand. It's possible that this additional handling should be moved to a separate object; the motivation for having it on Alignment __init__ is that it's easy for users to construct Alignment objects directly. Parameters: data: Data to convert into a SequenceCollection Names: Order of Names in the alignment. Should match the names of the sequences (after processing by label_to_name if present). Alphabet: Alphabet to use for the alignment (primarily important for DenseAlignment) MolType: MolType to be applied to the Alignment and to each seq. Name: Name of the SequenceCollection. Info: Info object to be attached to the alignment itself. conversion_f: Function to convert string into sequence. is_array: True if input is an array, False otherwise. force_same_data: True if data will be used as the same object. remove_duplicate_names: True if duplicate names are to be silently deleted instead of raising errors. label_to_name: if present, converts name into f(name). """ #read all the data in if we were passed a generator if isinstance(data, GeneratorType): data = list(data) #set the Name self.Name = Name #figure out alphabet and moltype self.Alphabet, self.MolType = \ self._get_alphabet_and_moltype(Alphabet, MolType, data) if not isinstance(Info, InfoClass): if Info: Info = InfoClass(Info) else: Info = InfoClass() self.Info = Info #if we're forcing the same data, skip the validation if force_same_data: self._force_same_data(data, Names) curr_seqs = data #otherwise, figure out what we got and coerce it into the right type else: per_seq_names, curr_seqs, name_order = \ self._names_seqs_order(conversion_f, data, Names, is_array, \ label_to_name, remove_duplicate_names, \ Alphabet=self.Alphabet) self.Names = name_order #will take only the seqs and names that are in name_order if per_seq_names != name_order: good_indices = [] for n in name_order: good_indices.append(per_seq_names.index(n)) if hasattr(curr_seqs, 'astype'): #it's an array #much faster to check than to raise exception in this case curr_seqs = take(curr_seqs, good_indices, axis=0) else: curr_seqs = [curr_seqs[i] for i in good_indices] per_seq_names = name_order #create NamedSeqs dict for fast lookups if not suppress_named_seqs: self.NamedSeqs = self._make_named_seqs(self.Names, curr_seqs) #Sequence objects behave like sequences of chars, so no difference #between Seqs and SeqData. Note that this differs for Alignments, #so be careful which you use if writing methods that should work for #both SequenceCollections and Alignments. self._set_additional_attributes(curr_seqs) def __str__(self): """Returns self in FASTA-format, respecting name order.""" return ''.join(['>%s\n%s\n' % (name, self.getGappedSeq(name)) for name in self.Names]) def _make_named_seqs(self, names, seqs): """Returns NamedSeqs: dict of name:seq.""" name_seq_tuples = zip(names, seqs) for n, s in name_seq_tuples: s.Name = n return dict(name_seq_tuples) def _set_additional_attributes(self, curr_seqs): """Sets additional attributes based on current seqs: class-specific.""" self.SeqData = curr_seqs self._seqs = curr_seqs try: self.SeqLen = max(map(len, curr_seqs)) except ValueError: #got empty sequence, for some reason? self.SeqLen = 0 def _force_same_data(self, data, Names): """Forces dict that was passed in to be used as self.NamedSeqs""" self.NamedSeqs = data self.Names = Names or data.keys() def copy(self): """Returns deep copy of self.""" result = self.__class__(self, MolType=self.MolType) def _get_alphabet_and_moltype(self, Alphabet, MolType, data): """Returns Alphabet and MolType, giving MolType precedence.""" if Alphabet is None and MolType is None: if hasattr(data, 'MolType'): MolType = data.MolType elif hasattr(data, 'Alphabet'): Alphabet = data.Alphabet #check for containers else: curr_item = self._get_container_item(data) if hasattr(curr_item, 'MolType'): MolType = curr_item.MolType elif hasattr(curr_item, 'Alphabet'): Alphabet = curr_item.Alphabet else: MolType = self.MolType #will be BYTES by default if Alphabet is not None and MolType is None: MolType = Alphabet.MolType if MolType is not None and Alphabet is None: try: Alphabet = MolType.Alphabets.DegenGapped except AttributeError: Alphabet = MolType.Alphabet return Alphabet, MolType def _get_container_item(self, data): """Checks container for item with Alphabet or MolType""" curr_item = None if hasattr(data, 'itervalues'): curr_item = data.itervalues().next() else: try: curr_item = iter(data).next() except: pass return curr_item def _strip_duplicates(self, names, seqs): """Internal function to strip duplicates from list of names""" if len(set(names)) == len(names): return set(), names, seqs #if we got here, there are duplicates unique_names = {} duplicates = {} fixed_names = [] fixed_seqs = [] for n, s in zip(names, seqs): if n in unique_names: duplicates[n] = 1 else: unique_names[n] = 1 fixed_names.append(n) fixed_seqs.append(s) if type(seqs) is ndarray: fixed_seqs = array(fixed_seqs, seqs.dtype) return duplicates, fixed_names, fixed_seqs def _names_seqs_order(self, conversion_f, data, Names, is_array, \ label_to_name, remove_duplicate_names, Alphabet=None): """Internal function to figure out names, seqs, and name_order.""" #figure out conversion function and whether it's an array if not conversion_f: input_type = self._guess_input_type(data) is_array = input_type in self.IsArray conversion_f = self.InputHandlers[input_type] #set seqs and names as properties if Alphabet: seqs, names = conversion_f(data, Alphabet=Alphabet) else: seqs, names = conversion_f(data) if names and label_to_name: names = map(label_to_name, names) curr_seqs = self._coerce_seqs(seqs, is_array) #if no names were passed in as Names, if we obtained them from #the seqs we should use them, but otherwise we should use the #default names if Names is None: if (names is None) or (None in names): per_seq_names = name_order = \ self.DefaultNameFunction(len(curr_seqs)) else: #got names from seqs per_seq_names = name_order = names else: #otherwise, names were passed in as Names: use this as the order #if we got names from the sequences, but otherwise assign the #names to successive sequences in order if (names is None) or (None in names): per_seq_names = name_order = Names else: #got names from seqs, so assume name_order is in Names per_seq_names = names name_order = Names #check for duplicate names duplicates, fixed_names, fixed_seqs = \ self._strip_duplicates(per_seq_names, curr_seqs) if duplicates: if remove_duplicate_names: per_seq_names, curr_seqs = fixed_names, fixed_seqs #if name_order doesn't have the same names as per_seq_names, #replace it with per_seq_names if (set(name_order) != set(per_seq_names)) or\ (len(name_order) != len(per_seq_names)): name_order = per_seq_names else: raise ValueError, \ "Some names were not unique. Duplicates are:\n" + \ str(sorted(duplicates.keys())) return per_seq_names, curr_seqs, name_order def _coerce_seqs(self, seqs, is_array): """Controls how seqs are coerced in _names_seqs_order. Override in subclasses where this behavior should differ. """ if is_array: seqs = map(str, map(self.MolType.ModelSeq, seqs)) return map(self.MolType.Sequence, seqs) def _guess_input_type(self, data): """Guesses input type of data; returns result as key of InputHandlers. First checks whether data is an Alignment, then checks for some common string formats, then tries to do it based on string or array properties. Returns 'empty' if check fails, i.e. if it can't recognize the sequence as a specific type. Note that bad sequences are not guaranteed to return 'empty', and may be recognized as another type incorrectly. """ if isinstance(data, DenseAlignment): return 'dense_aln' if isinstance(data, Alignment): return 'aln' if isinstance(data, SequenceCollection): return 'collection' if isinstance(data, dict): return 'dict' if isinstance(data, str): if data.startswith('>'): return 'fasta' else: return 'generic' first = None try: first = data[0] except (IndexError, TypeError): pass try: first = iter(data).next() except (IndexError, TypeError, StopIteration): pass if first is None: return 'empty' try: if isinstance(first, ModelSequence): #model sequence base type return 'model_seqs' elif hasattr(first, 'dtype'): #array object return 'array' elif isinstance(first, str) and first.startswith('>'): return 'fasta' else: try: dict(data) return 'kv_pairs' except (TypeError, ValueError): pass return 'generic' except (IndexError, TypeError), e: return 'empty' def __cmp__(self, other): """cmp first tests as dict, then as str.""" c = cmp(self.NamedSeqs, other) if not c: return 0 else: return cmp(str(self), str(other)) def keys(self): """keys uses self.Names, which defaults to known keys if None. Note: returns copy, not original. """ return self.Names[:] def values(self): """values returns values corresponding to self.Names.""" return [self.NamedSeqs[n] for n in self.Names] def items(self): """items returns (name, value) pairs.""" return [(n, self.NamedSeqs[n]) for n in self.Names] def iterSeqs(self, seq_order=None): """Iterates over values (sequences) in the alignment, in order. seq_order: list of keys giving the order in which seqs will be returned. Defaults to self.Names. Note that only these sequences will be returned, and that KeyError will be raised if there are sequences in order that have been deleted from the Alignment. If self.Names is None, returns the sequences in the same order as self.NamedSeqs.values(). Use map(f, self.seqs()) to apply the constructor f to each seq. f must accept a single list as an argument. Always returns references to the same objects that are values of the alignment. """ ns = self.NamedSeqs get = ns.__getitem__ for key in seq_order or self.Names: yield get(key) def _take_seqs(self): return list(self.iterSeqs()) Seqs = property(_take_seqs) #access as attribute if using default order. def takeSeqs(self, seqs, negate=False, **kwargs): """Returns new Alignment containing only specified seqs. Note that the seqs in the new alignment will be references to the same objects as the seqs in the old alignment. """ get = self.NamedSeqs.__getitem__ result = {} if negate: #copy everything except the specified seqs negated_names = [] row_lookup = dict.fromkeys(seqs) for r, row in self.NamedSeqs.items(): if r not in row_lookup: result[r] = row negated_names.append(r) seqs = negated_names #remember to invert the list of names else: #copy only the specified seqs for r in seqs: result[r] = get(r) if result: return self.__class__(result, Names=seqs, **kwargs) else: return {} #safe value; can't construct empty alignment def getSeqIndices(self, f, negate=False): """Returns list of keys of seqs where f(row) is True. List will be in the same order as self.Names, if present. """ get = self.NamedSeqs.__getitem__ #negate function if necessary if negate: new_f = lambda x: not f(x) else: new_f = f #get all the seqs where the function is True return [key for key in self.Names if new_f(get(key))] def takeSeqsIf(self, f, negate=False, **kwargs): """Returns new Alignment containing seqs where f(row) is True. Note that the seqs in the new Alignment are the same objects as the seqs in the old Alignment, not copies. """ #pass negate to get SeqIndices return self.takeSeqs(self.getSeqIndices(f, negate), **kwargs) def iterItems(self, seq_order=None, pos_order=None): """Iterates over elements in the alignment. seq_order (names) can be used to select a subset of seqs. pos_order (positions) can be used to select a subset of positions. Always iterates along a seq first, then down a position (transposes normal order of a[i][j]; possibly, this should change).. WARNING: Alignment.iterItems() is not the same as alignment.iteritems() (which is the built-in dict iteritems that iterates over key-value pairs). """ if pos_order: for row in self.iterSeqs(seq_order): for i in pos_order: yield row[i] else: for row in self.iterSeqs(seq_order): for i in row: yield i Items = property(iterItems) def getItems(self, items, negate=False): """Returns list containing only specified items. items should be a list of (row_key, col_key) tuples. """ get = self.NamedSeqs.__getitem__ if negate: #have to cycle through every item and check that it's not in #the list of items to return item_lookup = dict.fromkeys(map(tuple, items)) result = [] for r in self.Names: curr_row = get(r) for c in range(len(curr_row)): if (r, c) not in items: result.append(curr_row[c]) return result #otherwise, just pick the selected items out of the list else: return [get(row)[col] for row, col in items] def getItemIndices(self, f, negate=False): """Returns list of (key,val) tuples where f(self.NamedSeqs[key][val]).""" get = self.NamedSeqs.__getitem__ if negate: new_f = lambda x: not f(x) else: new_f = f result = [] for row_label in self.Names: curr_row = get(row_label) for col_idx, item in enumerate(curr_row): if new_f(item): result.append((row_label, col_idx)) return result def getItemsIf(self, f, negate=False): """Returns list of items where f(self.NamedSeqs[row][col]) is True.""" return self.getItems(self.getItemIndices(f, negate)) def getSimilar(self, target, min_similarity=0.0, max_similarity=1.0, \ metric=frac_same, transform=None): """Returns new Alignment containing sequences similar to target. target: sequence object to compare to. Can be in the alignment. min_similarity: minimum similarity that will be kept. Default 0.0. max_similarity: maximum similarity that will be kept. Default 1.0. (Note that both min_similarity and max_similarity are inclusive.) metric: similarity function to use. Must be f(first_seq, second_seq). The default metric is fraction similarity, ranging from 0.0 (0% identical) to 1.0 (100% identical). The Sequence classes have lots of methods that can be passed in as unbound methods to act as the metric, e.g. fracSameGaps. transform: transformation function to use on the sequences before the metric is calculated. If None, uses the whole sequences in each case. A frequent transformation is a function that returns a specified range of a sequence, e.g. eliminating the ends. Note that the transform applies to both the real sequence and the target sequence. WARNING: if the transformation changes the type of the sequence (e.g. extracting a string from an RnaSequence object), distance metrics that depend on instance data of the original class may fail. """ if transform: target = transform(target) m = lambda x: metric(target, x) if transform: def f(x): result = m(transform(x)) return min_similarity <= result <= max_similarity else: def f(x): result = m(x) return min_similarity <= result <= max_similarity return self.takeSeqsIf(f) def distanceMatrix(self, f): """Returns Matrix containing pairwise distances between sequences. f is the distance function f(x,y) -> distance between x and y. It's often useful to pass an unbound method in as f. Does not assume that f(x,y) == f(y,x) or that f(x,x) == 0. """ get = self.NamedSeqs.__getitem__ seqs = self.NamedSeqs.keys() result = Dict2D() for i in seqs: for j in seqs: d = f(get(i), get(j)) if i not in result: result[i] = {} if j not in result: result[j] = {} result[i][j] = d result[j][i] = d return result def isRagged(self): """Returns True if alignment has sequences of different lengths.""" seqs = self.Seqs #Get all sequences in alignment length = len(seqs[0]) #Get length of first sequence for seq in seqs: #If lengths differ if length != len(seq): return True #lengths were all equal return False def toPhylip(self, generic_label=True, make_seqlabel=None): """ Return alignment in PHYLIP format and mapping to sequence ids raises exception if invalid alignment Arguments: - make_seqlabel: callback function that takes the seq object and returns a label str """ return phylip_from_alignment(self, generic_label=generic_label, make_seqlabel=make_seqlabel) def toFasta(self, make_seqlabel=None): """Return alignment in Fasta format Arguments: - make_seqlabel: callback function that takes the seq object and returns a label str """ return fasta_from_alignment(self, make_seqlabel=make_seqlabel) def toNexus(self, seq_type, interleave_len=50): """ Return alignment in NEXUS format and mapping to sequence ids **NOTE** Not that every sequence in the alignment MUST come from a different species!! (You can concatenate multiple sequences from same species together before building tree) seq_type: dna, rna, or protein Raises exception if invalid alignment """ return nexus_from_alignment(self, seq_type, interleave_len=interleave_len) def getIntMap(self,prefix='seq_'): """Returns a dict with names mapped to enumerates integer names. - prefix: prefix for sequence label. Default = 'seq_' - int_keys is a dict mapping int names to sorted original names. """ get = self.NamedSeqs.__getitem__ int_keys = dict([(prefix+str(i),k) for i,k in \ enumerate(sorted(self.NamedSeqs.keys()))]) int_map = dict([(k, copy(get(v))) for k,v in int_keys.items()]) return int_map, int_keys def getNumSeqs(self): """Returns the number of sequences in the alignment.""" return len(self.NamedSeqs) def copyAnnotations(self, unaligned): """Copies annotations from seqs in unaligned to self, matching by name. Alignment programs like ClustalW don't preserve annotations, so this method is available to copy annotations off the unaligned sequences. unaligned should be a dictionary of Sequence instances. Ignores sequences that are not in self, so safe to use on larger dict of seqs that are not in the current collection/alignment. """ for name, seq in unaligned.items(): if name in self.NamedSeqs: self.NamedSeqs[name].copyAnnotations(seq) def annotateFromGff(self, f): """Copies annotations from gff-format file to self. Matches by name of sequence. This method expects a file handle, not the name of a file. Skips sequences in the file that are not in self. """ for (name, source, feature, start, end, score, strand, frame, attributes, comments) in GffParser(f): if name in self.NamedSeqs: self.NamedSeqs[name].addFeature( feature, parse_attributes(attributes), [(start,end)]) ''' self.NamedSeqs[seqname].data.addFeature( feature, parse_attributes(attributes), [(start, end)]) ''' def replaceSeqs(self, seqs): """Returns new alignment with same shape but with data taken from seqs. Primary use is for aligning codons from protein alignment, or, more generally, substituting in codons from a set of protein sequences (not necessarily aligned). For this reason, it takes characters from seqs three at a time rather than one at a time (i.e. 3 characters in seqs are put in place of 1 character in self). If seqs is an alignment, any gaps in it will be ignored. """ if hasattr(seqs, 'NamedSeqs'): seqs = seqs.NamedSeqs else: seqs = SequenceCollection(seqs).NamedSeqs new_seqs = [] for label in self.Names: aligned = self.NamedSeqs[label] seq = seqs[label] if isinstance(seq, Aligned): seq = seq.data new_seqs.append((label, Aligned(aligned.map * 3, seq))) return self.__class__(new_seqs) def getGappedSeq(self, seq_name, recode_gaps=False): """Return a gapped Sequence object for the specified seqname. Note: return type may depend on what data was loaded into the SequenceCollection or Alignment. """ return self.NamedSeqs[seq_name] def __add__(self, other): """Concatenates sequence data for same names""" aligned = isinstance(self, Alignment) if len(self.NamedSeqs) != len(other.NamedSeqs): raise ValueError("Alignments don't have same number of sequences") concatenated = [] for name in self.Names: if name not in other.Names: raise ValueError("Right alignment doesn't have a '%s'" % name) if aligned: new_seq = self.NamedSeqs[name].getGappedSeq() + \ other.getGappedSeq(name) else: new_seq = self.NamedSeqs[name] + other.getGappedSeq(name) concatenated.append(new_seq) new = self.__class__(MolType=self.MolType, data=zip(self.Names, concatenated)) if aligned: left = [a for a in self._shiftedAnnotations(new, 0) \ if a.map.End <= len(self)] right = [a for a in other._shiftedAnnotations(new, len(self)) \ if a.map.Start >= len(self)] new.annotations = left + right return new def addSeqs(self, other): """Adds sequences from other to self. Returns a new object. other must be of same class as self or coerceable to that class.. """ assert not isinstance(other, str), "Must provide a series of seqs "+\ "or an alignment" self_seq_class = self.Seqs[0].__class__ try: combined = self.Seqs + other.Seqs except AttributeError: combined = self.Seqs + list(other) for seq in combined: assert seq.__class__ == self_seq_class,\ "Seq classes different: Expected %s, Got %s" % \ (seq.__class__, self_seq_class) return self.__class__(data=combined) def writeToFile(self, filename=None, format=None, **kwargs): """Write the alignment to a file, preserving order of sequences. Arguments: - filename: name of the sequence file - format: format of the sequence file If format is None, will attempt to infer format from the filename suffix. """ if filename is None: raise DataError('no filename specified') # need to turn the alignment into a dictionary align_dict = {} for seq_name in self.Names: align_dict[seq_name] = str(self.NamedSeqs[seq_name]) if format is None and '.' in filename: # allow extension to work if provided format = filename[filename.rfind(".")+1:] if 'order' not in kwargs: kwargs['order'] = self.Names save_to_filename(align_dict, filename, format, **kwargs) def __len__(self): """len of SequenceCollection returns length of longest sequence.""" return self.SeqLen def getTranslation(self, gc=None, **kwargs): """Returns a new alignment object with the DNA sequences translated, using the current codon moltype, into an amino acid sequence. """ translated = [] aligned = isinstance(self, Alignment) # do the translation try: for seqname in self.Names: if aligned: seq = self.getGappedSeq(seqname) else: seq = self.NamedSeqs[seqname] pep = seq.getTranslation(gc) translated.append((seqname, pep)) return self.__class__(translated, **kwargs) except AttributeError, msg: raise AttributeError, "%s -- %s" % (msg, "Did you set a DNA MolType?") def getSeq(self, seqname): """Return a sequence object for the specified seqname. """ return self.NamedSeqs[seqname] def todict(self): """Returns the alignment as dict of names -> strings. Note: returns strings, NOT Sequence objects. """ align_dict = {} for seq_name in self.Names: align_dict[seq_name] = str(self.NamedSeqs[seq_name]) return align_dict def getPerSequenceAmbiguousPositions(self): """Returns dict of seq:{position:char} for ambiguous chars. Used in likelihood calculations. """ result = {} for name in self.Names: result[name] = ambig = {} for (i, motif) in enumerate(self.getGappedSeq(name)): if self.MolType.isAmbiguity(motif): ambig[i] = motif return result def degap(self, **kwargs): """Returns copy in which sequences have no gaps.""" new_seqs = [] aligned = isinstance(self, Alignment) for seq_name in self.Names: if aligned: seq = self.NamedSeqs[seq_name].data else: seq = self.NamedSeqs[seq_name] new_seqs.append((seq_name, seq.degap())) return SequenceCollection(MolType=self.MolType, data=new_seqs, **kwargs) def withModifiedTermini(self): """Changes the termini to include termini char instead of gapmotif. Useful to correct the standard gap char output by most alignment programs when aligned sequences have different ends. """ seqs = [] for name in self.Names: seq = self.NamedSeqs[name].withTerminiUnknown() seqs.append((name, seq)) return self.__class__(MolType=self.MolType, data=seqs) def hasTerminalStops(self, gc=None): """Returns True if any sequence has a terminal stop codon.""" stops = [] aligned = isinstance(self, Alignment) for seq_name in self.Names: if aligned: seq = self.NamedSeqs[seq_name].data else: seq = self.NamedSeqs[seq_name] stops.append(seq.hasTerminalStop(gc=gc)) return max(stops) def withoutTerminalStopCodons(self, gc=None, **kwargs): """Removes any terminal stop codons from the sequences""" new_seqs = [] aligned = isinstance(self, Alignment) for seq_name in self.Names: old_seq = self.NamedSeqs[seq_name] if aligned: new_seq = old_seq.data.withoutTerminalStopCodon(gc=gc) new_seq = Aligned(old_seq.map, new_seq) else: new_seq = old_seq.withoutTerminalStopCodon(gc=gc) new_seqs.append((seq_name, new_seq)) return self.__class__(MolType=self.MolType, data=new_seqs, **kwargs) def getSeqNames(self): """Return a list of sequence names.""" return self.Names[:] def getMotifProbs(self, alphabet=None, include_ambiguity=False, exclude_unobserved=False, allow_gap=False, pseudocount=0): """Return a dictionary of motif probs. Arguments: - include_ambiguity: if True resolved ambiguous codes are included in estimation of frequencies, default is False. - exclude_unobserved: if True, motifs that are not present in the alignment are excluded from the returned dictionary, default is False. - allow_gap: allow gap motif """ if alphabet is None: alphabet = self.MolType.Alphabet if allow_gap: alphabet = alphabet.Gapped counts = {} for seq_name in self.Names: sequence = self.NamedSeqs[seq_name] motif_len = alphabet.getMotifLen() if motif_len > 1: posns = range(0, len(sequence)+1-motif_len, motif_len) sequence = [sequence[i:i+motif_len] for i in posns] for motif in sequence: if not allow_gap: if self.MolType.Gap in motif: continue if motif in counts: counts[motif] += 1 else: counts[motif] = 1 probs = {} if not exclude_unobserved: for motif in alphabet: probs[motif] = pseudocount for (motif, count) in counts.items(): motif_set = alphabet.resolveAmbiguity(motif) if len(motif_set) > 1: if include_ambiguity: count = float(count) / len(motif_set) else: continue for motif in motif_set: probs[motif] = probs.get(motif, pseudocount) + count total = float(sum(probs.values())) for motif in probs: probs[motif] /= total return probs def getGapCount(self, seq_name): return len(self.NamedSeqs[seq_name].map.gaps()) def getSeqFreqs(self): """Returns Profile of counts: seq by character. See documentation for _get_freqs: this just wraps it and converts the result into a Profile object organized per-sequence (i.e. per row). """ return Profile(self._get_freqs(0), self.Alphabet) def _make_gaps_ok(self, allowed_gap_frac): """Makes the gaps_ok function used by omitGapPositions and omitGapSeqs. Need to make the function because if it's a method of Alignment, it has unwanted 'self' and 'allowed_gap_frac' parameters that impede the use of map() in takeSeqsIf. WARNING: may not work correctly if component sequences have gaps that are not the Alignment gap character. This is because the gaps are checked at the position level (and the positions are lists), rather than at the sequence level. Working around this issue would probably cause a significant speed penalty. """ def gaps_ok(seq): seq_len = len(seq) try: num_gaps = seq.countGaps() except AttributeError: num_gaps = len(filter(self.MolType.Gaps.__contains__, seq)) return num_gaps / seq_len <= allowed_gap_frac return gaps_ok def omitGapPositions(self, allowed_gap_frac=1-eps, del_seqs=False, \ allowed_frac_bad_cols=0, seq_constructor=None): """Returns new alignment where all cols have <= allowed_gap_frac gaps. allowed_gap_frac says what proportion of gaps is allowed in each column (default is 1-eps, i.e. all cols with at least one non-gap character are preserved). If del_seqs is True (default:False), deletes the sequences that don't have gaps where everything else does. Otherwise, just deletes the corresponding column from all sequences, in which case real data as well as gaps can be removed. Uses seq_constructor(seq) to make each new sequence object. Note: a sequence that is all gaps will not be deleted by del_seqs (even if all the positions have been deleted), since it has no non-gaps in positions that are being deleted for their gap content. Possibly, this decision should be revisited since it may be a surprising result (and there are more convenient ways to return the sequences that consist wholly of gaps). """ if seq_constructor is None: seq_constructor = self.MolType.Sequence gaps_ok = self._make_gaps_ok(allowed_gap_frac) #if we're not deleting the 'naughty' seqs that contribute to the #gaps, it's easy... if not del_seqs: return self.takePositionsIf(f=gaps_ok, \ seq_constructor=seq_constructor) #otherwise, we have to figure out which seqs to delete. #if we get here, we're doing del_seqs. cols_to_delete = dict.fromkeys(self.getPositionIndices(gaps_ok, \ negate=True)) default_gap_f = self.MolType.Gaps.__contains__ bad_cols_per_row = {} for key, row in self.NamedSeqs.items(): try: is_gap = row.Alphabet.Gaps.__contains__ except AttributeError: is_gap = default_gap_f for col in cols_to_delete: if not is_gap(str(row)[col]): if key not in bad_cols_per_row: bad_cols_per_row[key] = 1 else: bad_cols_per_row[key] += 1 #figure out which of the seqs we're deleting get = self.NamedSeqs.__getitem__ seqs_to_delete = {} for key, count in bad_cols_per_row.items(): if float(count)/len(get(key)) >= allowed_frac_bad_cols: seqs_to_delete[key] = True #It's _much_ more efficient to delete the seqs before the cols. good_seqs = self.takeSeqs(seqs_to_delete, negate=True) cols_to_keep = dict.fromkeys(range(self.SeqLen)) for c in cols_to_delete: del cols_to_keep[c] if good_seqs: return good_seqs.takePositions(cols=cols_to_keep.keys(), \ seq_constructor=seq_constructor) else: return {} def omitGapSeqs(self, allowed_gap_frac=0): """Returns new alignment with seqs that have <= allowed_gap_frac. allowed_gap_frac should be a fraction between 0 and 1 inclusive. Default is 0. """ gaps_ok = self._make_gaps_ok(allowed_gap_frac) return self.takeSeqsIf(gaps_ok) def omitGapRuns(self, allowed_run=1): """Returns new alignment where all seqs have runs of gaps <=allowed_run. Note that seqs with exactly allowed_run gaps are not deleted. Default is for allowed_run to be 1 (i.e. no consecutive gaps allowed). Because the test for whether the current gap run exceeds the maximum allowed gap run is only triggered when there is at least one gap, even negative values for allowed_run will still let sequences with no gaps through. """ def ok_gap_run(x): try: is_gap = x.Alphabet.Gaps.__contains__ except AttributeError: is_gap = self.MolType.Gaps.__contains__ curr_run = max_run = 0 for i in x: if is_gap(i): curr_run += 1 if curr_run > allowed_run: return False else: curr_run = 0 #can only get here if max_run was never exceeded (although this #does include the case where the sequence is empty) return True return self.takeSeqsIf(ok_gap_run) def omitSeqsTemplate(self, template_name, gap_fraction, gap_run): """Returns new alignment where all seqs are well aligned with template. gap_fraction = fraction of positions that either have a gap in the template but not in the seq or in the seq but not in the template gap_run = number of consecutive gaps tolerated in query relative to sequence or sequence relative to query """ template = self.NamedSeqs[template_name] gap_filter = make_gap_filter(template, gap_fraction, gap_run) return self.takeSeqsIf(gap_filter) def toDna(self): """Returns the alignment as DNA.""" new = {} aligned = isinstance(self, Alignment) for name, seq in self.NamedSeqs.items(): if aligned: seq = seq.getGappedSeq() new[name] = seq.toDna() return self.__class__(data=new, Name = self.Name, Info = self.Info) def toRna(self): """Returns the alignment as RNA""" new = {} aligned = isinstance(self, Alignment) for name, seq in self.NamedSeqs.items(): if aligned: seq = seq.getGappedSeq() new[name] = seq.toRna() return self.__class__(data=new, Name = self.Name, Info = self.Info) def rc(self): """Returns the reverse complement alignment""" new = {} aligned = isinstance(self, Alignment) for name, seq in self.NamedSeqs.items(): if aligned: seq = seq.getGappedSeq() new[name] = seq.rc() rc = self.__class__(data=new, Name = self.Name, Info = self.Info) if isinstance(self, _Annotatable): self._annotations_nucleic_reversed_on(rc) return rc def reversecomplement(self): """Returns the reverse complement alignment. A synonymn for rc.""" return self.rc() def padSeqs(self, pad_length=None, **kwargs): """Returns copy in which sequences are padded to same length. pad_length: Length all sequences are to be padded to. Will pad to max sequence length if pad_length is None or less than max length. """ #get max length max_len = max([len(s) for s in self.Seqs]) #If a pad_length was passed in, make sure it is valid if pad_length is not None: pad_length = int(pad_length) if pad_length < max_len: raise ValueError, \ "pad_length must be at greater or equal to maximum sequence length: %s"\ %(str(max_len)) #pad_length is max sequence length. else: pad_length = max_len #Get new sequence list new_seqs = [] aligned = isinstance(self, Alignment) #for each sequence, pad gaps to end for seq_name in self.Names: if aligned: seq = self.NamedSeqs[seq_name].data else: seq = self.NamedSeqs[seq_name] padded_seq = seq + '-'*(pad_length-len(seq)) new_seqs.append((seq_name, padded_seq)) #return new SequenceCollection object return SequenceCollection(MolType=self.MolType, data=new_seqs, **kwargs) class Aligned(object): """One sequence in an alignment, a map between alignment coordinates and sequence coordinates""" def __init__(self, map, data, length=None): #Unlike the normal map constructor, here we take a list of pairs of #alignment coordinates, NOT a list of pairs of sequence coordinates if isinstance(map, list): map = Map(map, parent_length=length).inverse() self.map = map self.data = data if hasattr(data, 'Info'): self.Info = data.Info if hasattr(data, 'Name'): self.Name = data.Name def copy(self, memo=None, _nil=[], constructor='ignored'): """Returns a shallow copy of self WARNING: cogent.core.sequence.Sequence does NOT implement a copy method, as such, the data member variable of the copied object will maintain reference to the original object. WARNING: cogent.core.location.Map does NOT implement a copy method, as such, the data member variable of the copied object will maintain reference to the original object. """ return self.__class__(self.map, self.data) def __repr__(self): return '%s of %s' % (repr(self.map), repr(self.data)) def withTerminiUnknown(self): return self.__class__(self.map.withTerminiUnknown(), self.data) def copyAnnotations(self, other): self.data.copyAnnotations(other) def annotateFromGff(self, f): self.data.annotate_from_gff(f) def addFeature(self, *args, **kwargs): self.data.addFeature(*args, **kwargs) def __str__(self): """Returns string representation of aligned sequence, incl. gaps.""" return str(self.getGappedSeq()) def __cmp__(self, other): """Compares based on string representations.""" return cmp(str(self), str(other)) def __iter__(self): """Iterates over sequence one motif (e.g. char) at a time, incl. gaps""" return self.data.gappedByMapMotifIter(self.map) def getGappedSeq(self, recode_gaps=False): """Returns sequence as an object, including gaps.""" return self.data.gappedByMap(self.map, recode_gaps) def __len__(self): # these make it look like Aligned should be a subclass of Map, # but then you have to be careful with __getitem__, __init__ and inverse. return len(self.map) def __getitem__(self, slice): return Aligned(self.map[slice], self.data) def getTracks(self, policy): policy = policy.at(self.map.inverse()) return self.data.getTracks(policy) def remappedTo(self, map): #assert map is self.parent_map or ... ? #print 'REMAP', self.map, self #print 'ONTO', map, map.inverse() result = Aligned(map[self.map.inverse()].inverse(), self.data) #print 'GIVES', result.map, result #print return result def getAnnotationsMatching(self, alignment, *args): for annot in self.data.getAnnotationsMatching(*args): yield annot.remappedTo(alignment, self.map.inverse()) def gapVector(self): """Returns gapVector of GappedSeq, for omitGapPositions.""" return self.getGappedSeq().gapVector() def _masked_annotations(self, annot_types, mask_char, shadow): """returns a new aligned sequence with regions defined by align_spans and shadow masked.""" new_data = self.data.withMaskedAnnotations(annot_types, mask_char, shadow) # we remove the mask annotations from self and new_data return self.__class__(self.map, new_data) class AlignmentI(object): """Alignment interface object. Contains methods shared by implementations. Note that subclasses should inherit both from AlignmentI and from SequenceCollection (typically). Alignments are expected to be immutable once created. No mechanism is provided for maintaining reference consistency if data in the alignment are modified. An Alignment is expected to be able to generate the following: - Seqs: Sequence objects in the alignment, can turn themselves into strings. These are usually thought of as "rows" in an alignment. - Positions: Vectors representing data in each position in the alignment These are usually thought of as "columns" in an alignment. - SeqData: Vectors representing data in each sequence in the alignment, not necessarily guaranteed to turn themselves into a string - Items: Iterator over the characters in the alignment - Names: List of names of sequences in the alignment. Used for display order. A cheap way to omit or reorder sequences is to modify the list of names. - NamedSeqs: Dict of name -> seq object, used for lookup. - MolType: MolType of the alignment. """ DefaultGap = '-' #default gap character for padding GapChars = dict.fromkeys('-?') #default gap chars for comparisons def iterPositions(self, pos_order=None): """Iterates over positions in the alignment, in order. pos_order refers to a list of indices (ints) specifying the column order. This lets you rearrange positions if you want to (e.g. to pull out individual codon positions). Note that self.iterPositions() always returns new objects, by default lists of elements. Use map(f, self.iterPositions) to apply the constructor or function f to the resulting lists (f must take a single list as a parameter). Note that some sequences (e.g. ViennaStructures) have rules that prevent arbitrary strings of their symbols from being valid objects. Will raise IndexError if one of the indices in order exceeds the sequence length. This will always happen on ragged alignments: assign to self.SeqLen to set all sequences to the same length. """ get = self.NamedSeqs.__getitem__ pos_order = pos_order or xrange(self.SeqLen) seq_order = self.Names for pos in pos_order: yield [get(seq)[pos] for seq in seq_order] Positions = property(iterPositions) def takePositions(self, cols, negate=False, seq_constructor=None): """Returns new Alignment containing only specified positions. By default, the seqs will be lists, but an alternative constructor can be specified. Note that takePositions will fail on ragged positions. """ if seq_constructor is None: seq_constructor = self.MolType.Sequence result = {} #if we're negating, pick out all the positions except specified indices if negate: col_lookup = dict.fromkeys(cols) for key, row in self.NamedSeqs.items(): result[key] = seq_constructor([row[i] for i in range(len(row)) \ if i not in col_lookup]) #otherwise, just get the requested indices else: for key, row in self.NamedSeqs.items(): result[key] = seq_constructor([row[i] for i in cols]) return self.__class__(result, Names=self.Names) def getPositionIndices(self, f, negate=False): """Returns list of column indices for which f(col) is True.""" #negate f if necessary if negate: new_f = lambda x: not f(x) else: new_f = f return [i for i, col in enumerate(self.Positions) if new_f(col)] def takePositionsIf(self, f, negate=False, seq_constructor=None): """Returns new Alignment containing cols where f(col) is True. Note that the seqs in the new Alignment are always new objects. Default constructor is list(), but an alternative can be passed in. """ if seq_constructor is None: seq_constructor = self.MolType.Sequence return self.takePositions(self.getPositionIndices(f, negate), \ seq_constructor=seq_constructor) def IUPACConsensus(self, alphabet=None): """Returns string containing IUPAC consensus sequence of the alignment. """ if alphabet is None: alphabet = self.MolType consensus = [] degen = alphabet.degenerateFromSequence for col in self.Positions: consensus.append(degen(coerce_to_string(col))) return coerce_to_string(consensus) def columnFreqs(self, constructor=Freqs): """Returns list of Freqs with item counts for each column. """ return map(constructor, self.Positions) def columnProbs(self, constructor=Freqs): """Returns FrequencyDistribuutions w/ prob. of each item per column. Implemented as a list of normalized Freqs objects. """ freqs = self.columnFreqs(constructor) for fd in freqs: fd.normalize() return freqs def majorityConsensus(self, transform=None, constructor=Freqs): """Returns list containing most frequent item at each position. Optional parameter transform gives constructor for type to which result will be converted (useful when consensus should be same type as originals). """ col_freqs = self.columnFreqs(constructor) consensus = [freq.Mode for freq in col_freqs] if transform == str: return coerce_to_string(consensus) elif transform: return transform(consensus) else: return consensus def uncertainties(self, good_items=None): """Returns Shannon uncertainty at each position. Usage: information_list = alignment.information(good_items=None) If good_items is supplied, deletes any symbols that are not in good_items. """ uncertainties = [] #calculate column probabilities if necessary if hasattr(self, 'PositionumnProbs'): probs = self.PositionumnProbs else: probs = self.columnProbs() #calculate uncertainty for each column for prob in probs: #if there's a list of valid symbols, need to delete everything else if good_items: prob = prob.copy() #do not change original #get rid of any symbols not in good_items for symbol in prob.keys(): if symbol not in good_items: del prob[symbol] #normalize the probabilities and add to the list prob.normalize() uncertainties.append(prob.Uncertainty) return uncertainties def scoreMatrix(self): """Returns a position specific score matrix for the alignment.""" return Dict2D(dict([(i,Freqs(col)) for i, col in enumerate(self.Positions)])) def _get_freqs(self, index=None): """Gets array of freqs along index 0 (= positions) or 1 (= seqs). index: if 0, will calculate the frequency of each symbol in each position (=column) in the alignment. Will return 2D array where the first index is the position, and the second index is the index of the symbol in the alphabet. For example, for the TCAG DNA Alphabet, result[3][0] would store the count of T at position 3 (i.e. the 4th position in the alignment. if 1, does the same thing except that the calculation is performed for each sequence, so the 2D array has the sequence index as the first index, and the symbol index as the second index. For example, for the TCAG DNA Alphabet, result[3][0] would store the count of T in the sequence at index 3 (i.e. the 4th sequence). First an DenseAligment object is created, next the calculation is done on this object. It is important that the DenseAlignment is initialized with the same MolType and Alphabet as the original Alignment. """ da = DenseAlignment(self, MolType=self.MolType, Alphabet=self.Alphabet) return da._get_freqs(index) def getPosFreqs(self): """Returns Profile of counts: position by character. See documentation for _get_freqs: this just wraps it and converts the result into a Profile object organized per-position (i.e. per column). """ return Profile(self._get_freqs(1), self.Alphabet) def sample(self, n=None, with_replacement=False, motif_length=1, \ randint=randint, permutation=permutation): """Returns random sample of positions from self, e.g. to bootstrap. Arguments: - n: the number of positions to sample from the alignment. Default is alignment length - with_replacement: boolean flag for determining if sampled positions - random_series: a random number generator with .randint(min,max) .random() methods Notes: By default (resampling all positions without replacement), generates a permutation of the positions of the alignment. Setting with_replacement to True and otherwise leaving parameters as defaults generates a standard bootstrap resampling of the alignment. """ population_size = len(self) // motif_length if not n: n = population_size if with_replacement: locations = randint(0, population_size, n) else: assert n <= population_size, (n, population_size, motif_length) locations = permutation(population_size)[:n] positions = [(loc*motif_length, (loc+1)*motif_length) for loc in locations] sample = Map(positions, parent_length=len(self)) return self.gappedByMap(sample, Info=self.Info) def slidingWindows(self, window, step): """Generator yielding new Alignments of given length and interval. Arguments: - window: The length of each returned alignment. - step: The interval between the start of the successive alignment objects returned. """ for pos in range(0, len(self)-window+1,step): yield self[pos:pos+window] def aln_from_array(a, array_type=None, Alphabet=None): """Alignment from array of pos x seq: no change, names are integers. This is an InputHandler for Alignment. It converts an arbitrary array of numbers without change, but adds successive integer names (0-based) to each sequence (i.e. column) in the input a. Data type of input is unchanged. """ if array_type is None: result = a.copy() else: result = a.astype(array_type) return transpose(result), None def aln_from_model_seqs(seqs, array_type=None, Alphabet=None): """Alignment from ModelSequence objects: seqs -> array, names from seqs. This is an InputHandler for Alignment. It converts a list of Sequence objects with _data and Label properties into the character array Alignment needs. All sequences must be the same length. WARNING: Assumes that the ModelSeqs are already in the right alphabet. If this is not the case, e.g. if you are putting sequences on a degenerate alphabet into a non-degenerate alignment or you are putting protein sequences into a DNA alignment, there will be problems with the alphabet mapping (i.e. the resulting sequences may be meaningless). WARNING: Data type of return array is not guaranteed -- check in caller! """ data, names = [], [] for s in seqs: data.append(s._data) names.append(s.Name) result = array(data) if array_type: result = result.astype(array_type) return result, names def aln_from_generic(data, array_type=None, Alphabet=None): """Alignment from generic seq x pos data: sequence of sequences of chars. This is an InputHandler for Alignment. It converts a generic list (each item in the list will be mapped onto an Array object, with character transformations, all items must be the same length) into a numpy array, and assigns sequential integers (0-based) as names. WARNING: Data type of return array is not guaranteed -- check in caller! """ result = array(map(Alphabet.toIndices, data)) names = [] for d in data: if hasattr(d, 'Name'): names.append(d.Name) else: names.append(None) if array_type: result = result.astype(array_type) return result, names def aln_from_collection(seqs, array_type=None, Alphabet=None): """Alignment from SequenceCollection object, or its subclasses.""" names = seqs.Names data = [seqs.NamedSeqs[i] for i in names] result = array(map(Alphabet.toIndices, data)) if array_type: result = result.astype(array_type) return result, names def aln_from_fasta(seqs, array_type=None, Alphabet=None): """Alignment from FASTA-format string or lines. This is an InputHandler for Alignment. It converts a FASTA-format string or collection of lines into an Alignment object. All sequences must be the same length. WARNING: Data type of return array is not guaranteed -- check in caller! """ if isinstance(seqs, str): seqs = seqs.splitlines() return aln_from_model_seqs([ModelSequence(s, Name=l, Alphabet=Alphabet)\ for l, s in cogent.parse.fasta.MinimalFastaParser(seqs)], array_type) def aln_from_dict(aln, array_type=None, Alphabet=None): """Alignment from dict of {label:seq_as_str}. This is an InputHandler for Alignment. It converts a dict in which the keys are the names and the values are the sequences (sequence only, no whitespace or other formatting) into an alignment. Because the dict doesn't preserve order, the result will be in alphabetical order.""" names, seqs = zip(*sorted(aln.items())) result = array(map(Alphabet.toIndices, seqs), array_type) return result, list(names) def aln_from_kv_pairs(aln, array_type=None, Alphabet=None): """Alignment from sequence of (key, value) pairs. This is an InputHandler for Alignment. It converts a list in which the first item of each pair is the label and the second item is the sequence (sequence only, no whitespace or other formatting) into an alignment. Because the dict doesn't preserve order, the result will be in arbitrary order.""" names, seqs = zip(*aln) result = array(map(Alphabet.toIndices, seqs), array_type) return result, list(names) def aln_from_dense_aln(aln, array_type=None, Alphabet=None): """Alignment from existing DenseAlignment object: copies data. Retrieves data from Positions field. Uses copy(), so array data type should be unchanged. """ if array_type is None: result = aln.ArrayPositions.copy() else: result = aln.ArrayPositions.astype(array_type) return transpose(result), aln.Names[:] def aln_from_empty(obj, *args, **kwargs): """Alignment from empty data: raise exception.""" raise ValueError, "Cannot create empty alignment." #Implementation of Alignment base class class DenseAlignment(AlignmentI, SequenceCollection): """Holds a dense array representing a multiple sequence alignment. An Alignment is _often_, but not necessarily, an array of chars. You might want to use some other data type for the alignment if you have a large number of symbols. For example, codons on an ungapped DNA alphabet has 4*4*4=64 entries so can fit in a standard char data type, but tripeptides on the 20-letter ungapped protein alphabet has 20*20*20=8000 entries so can _not_ fit in a char and values will wrap around (i.e. you will get an unpredictable, wrong value for any item whose index is greater than the max value, e.g. 255 for uint8), so in this case you would need to use UInt16, which can hold 65536 values. DO NOT USE SIGNED DATA TYPES FOR YOUR ALIGNMENT ARRAY UNLESS YOU LOVE MISERY AND HARD-TO-DEBUG PROBLEMS. Implementation: aln[i] returns position i in the alignment. aln.Positions[i] returns the same as aln[i] -- usually, users think of this as a 'column', because alignment editors such as Clustal typically display each sequence as a row so a position that cuts across sequences is a column. aln.Seqs[i] returns a sequence, or 'row' of the alignment in standard terminology. WARNING: aln.Seqs and aln.Positions are different views of the same array, so if you change one you will change the other. This will no longer be true if you assign to Seqs or Positions directly, so don't do it. If you want to change the data in the whole array, always assign to a slice so that both views update: aln.Seqs[:] = x instead of aln.Seqs = x. If you get the two views out of sync, you will get all sorts of exceptions. No validation is performed on aln.Seqs and aln.Positions for performance reasons, so this can really get you into trouble. Alignments are immutable, though this is not enforced. If you change the data after the alignment is created, all sorts of bad things might happen. Class properties: Alphabet: should be an Alphabet object. Must provide mapping between items (possibly, but not necessarily, characters) in the alignment and indices of those characters in the resulting Alignment object. SequenceType: Constructor to use when building sequences. Default: Sequence InputHandlers: dict of {input_type:input_handler} where input_handler is from the InputHandlers above and input_type is a result of the method self._guess_input_type (should always be a string). Creating a new array will always result in a new object unless you use the force_same_object=True parameter. WARNING: Rebinding the Names attribute in a DenseAlignment is not recommended because not all methods will use the updated name order. This is because the original sequence and name order are used to produce data structures that are cached for efficiency, and are not updated if you change the Names attribute. WARNING: DenseAlignment strips off Info objects from sequences that have them, primarily for efficiency. """ MolType = None #will be set to BYTES on moltype import Alphabet = None #will be set to BYTES.Alphabet on moltype import InputHandlers = { 'array':aln_from_array, 'model_seqs':aln_from_model_seqs, 'generic':aln_from_generic, 'fasta':aln_from_fasta, 'dense_aln':aln_from_dense_aln, 'aln':aln_from_collection, 'collection':aln_from_collection, 'dict':aln_from_dict, 'kv_pairs':aln_from_kv_pairs, 'empty':aln_from_empty, } def __init__(self, *args, **kwargs): """Returns new DenseAlignment object. Inherits from SequenceCollection. """ kwargs['suppress_named_seqs'] = True super(DenseAlignment, self).__init__(*args, **kwargs) self.ArrayPositions = transpose(\ self.SeqData.astype(self.Alphabet.ArrayType)) self.ArraySeqs = transpose(self.ArrayPositions) self.SeqData = self.ArraySeqs self.SeqLen = len(self.ArrayPositions) def _force_same_data(self, data, Names): """Forces array that was passed in to be used as self.ArrayPositions""" if isinstance(data, DenseAlignment): data = data._positions self.ArrayPositions = data self.Names = Names or self.DefaultNameFunction(len(data[0])) def _get_positions(self): """Override superclass Positions to return positions as symbols.""" return map(self.Alphabet.fromIndices, self.ArrayPositions) Positions = property(_get_positions) def _get_named_seqs(self): if not hasattr(self, '_named_seqs'): seqs = map(self.Alphabet.toString, self.ArraySeqs) if self.MolType: seqs = map(self.MolType.Sequence, seqs) self._named_seqs = self._make_named_seqs(self.Names, seqs) return self._named_seqs NamedSeqs = property(_get_named_seqs) def keys(self): """Supports dict-like interface: returns names as keys.""" return self.Names def values(self): """Supports dict-like interface: returns seqs as Sequence objects.""" return [self.Alphabet.MolType.ModelSeq(i, Alphabet=self.Alphabet) \ for i in self.ArraySeqs] def items(self): """Supports dict-like interface; returns (name, seq) pairs.""" return zip(self.keys(), self.values()) def __iter__(self): """iter(aln) iterates over positions, returning array slices. Each item in the result is be a position ('column' in standard terminology) within the alignment, with the sequneces in the same order as in the names. The result shares data with the original array, so if you change the result you change the Alignment. """ return iter(self.Positions) def __getitem__(self, item): """getitem delegates to self.Positions., returning array slices. The result is a column or slice of columns, supporting full slice functionality (including stride). Use this to get a selection of positions from the alignment. Result shares data with the original array, so if you change the result you change the Alignment. """ return self.Positions[item] def _coerce_seqs(self, seqs, is_array): """Controls how seqs are coerced in _names_seqs_order. Override in subclasses where this behavior should differ. """ return seqs def getSubAlignment(self, seqs=None, pos=None, invert_seqs=False, \ invert_pos=False): """Returns subalignment of specified sequences and positions. seqs and pos can be passed in as lists of sequence indices to keep or positions to keep. invert_seqs: if True (default False), gets everything _except_ the specified sequences. invert_pos: if True (default False), gets everything _except_ the specified positions. Unlike most of the other code that gets things out of an alignment, this method returns a new alignment that does NOT share data with the original alignment. """ #figure out which positions to keep, and keep them if pos is not None: if invert_pos: pos_mask = ones(len(self.ArrayPositions)) put(pos_mask, pos, 0) pos = nonzero(pos_mask)[0] data = take(self.ArrayPositions, pos, axis=0) else: data = self.ArrayPositions #figure out which sequences to keep, and keep them if seqs is not None: if invert_seqs: seq_mask = ones(len(self.ArraySeqs)) put(seq_mask, seqs, 0) seqs = nonzero(seq_mask)[0] data = take(data, seqs, 1) names = [self.Names[i] for i in seqs] else: names = self.Names return self.__class__(data, map(str,names), self.Alphabet, \ conversion_f=aln_from_array) def __str__(self): """Returns FASTA-format string. Should be able to handle joint alphabets, e.g. codons. """ result = [] names = map(str, self.Names) max_label_length = max(map(len, names)) + 1 seq2str = self.Alphabet.fromIndices for l, s in zip(self.Names, self.ArraySeqs): result.append('>'+str(l)+'\n'+''.join(seq2str(s))) return '\n'.join(result) + '\n' def _get_freqs(self, index=None): """Gets array of freqs along index 0 (= positions) or 1 (= seqs). index: if 0, will calculate the frequency of each symbol in each position (=column) in the alignment. Will return 2D array where the first index is the position, and the second index is the index of the symbol in the alphabet. For example, for the TCAG DNA Alphabet, result[3][0] would store the count of T at position 3 (i.e. the 4th position in the alignment. if 1, does the same thing except that the calculation is performed for each sequence, so the 2D array has the sequence index as the first index, and the symbol index as the second index. For example, for the TCAG DNA Alphabet, result[3][0] would store the count of T in the sequence at index 3 (i.e. the 4th sequence). """ if index: a = self.ArrayPositions else: a = self.ArraySeqs count_f = self.Alphabet.counts return array(map(count_f, a)) def getPosFreqs(self): """Returns Profile of counts: position by character. See documentation for _get_freqs: this just wraps it and converts the result into a Profile object organized per-position (i.e. per column). """ return Profile(self._get_freqs(1), self.Alphabet) def getSeqEntropy(self): """Returns array containing Shannon entropy for each seq in self. Uses the profile object from getSeqFreqs (see docstring) to calculate the per-symbol entropy in each sequence in the alignment, i.e. the uncertainty about each symbol in each sequence (or row). This can be used to, for instance, filter low-complexity sequences. """ p = self.getSeqFreqs() p.normalizePositions() return p.rowUncertainty() def getPosEntropy(self): """Returns array containing Shannon entropy for each pos in self. Uses the profile object from getPosFreqs (see docstring) to calculate the per-symbol entropy in each position in the alignment, i.e. the uncertainty about each symbol at each position (or column). This can be used to, for instance, detect the level of conservation at each position in an alignment. """ p = self.getPosFreqs() p.normalizePositions() return p.rowUncertainty() def IUPACConsensus(self, alphabet=None): """Returns string containing IUPAC consensus sequence of the alignment. """ if alphabet is None: alphabet = self.MolType consensus = [] degen = alphabet.degenerateFromSequence for col in self.Positions: consensus.append(degen(str(alphabet.ModelSeq(col, \ Alphabet=alphabet.Alphabets.DegenGapped)))) return coerce_to_string(consensus) def _make_gaps_ok(self, allowed_gap_frac): """Makes the gaps_ok function used by omitGapPositions and omitGapSeqs. Need to make the function because if it's a method of Alignment, it has unwanted 'self' and 'allowed_gap_frac' parameters that impede the use of map() in takeSeqsIf. WARNING: may not work correctly if component sequences have gaps that are not the Alignment gap character. This is because the gaps are checked at the column level (and the positions are lists), rather than at the row level. Working around this issue would probably cause a significant speed penalty. """ def gaps_ok(seq): seq_len = len(seq) if hasattr(seq, 'countGaps'): num_gaps = seq.countGaps() elif hasattr(seq, 'count'): num_gaps = seq.count(self.Alphabet.Gap) else: num_gaps = sum(seq==self.Alphabet.GapIndex) return num_gaps / seq_len <= allowed_gap_frac return gaps_ok def columnFreqs(self, constructor=Freqs): """Returns list of Freqs with item counts for each column. """ return map(constructor, self.Positions) def sample(self, n=None, with_replacement=False, motif_length=1, \ randint=randint, permutation=permutation): """Returns random sample of positions from self, e.g. to bootstrap. Arguments: - n: the number of positions to sample from the alignment. Default is alignment length - with_replacement: boolean flag for determining if sampled positions - randint and permutation: functions for random integer in a specified range, and permutation, respectively. Notes: By default (resampling all positions without replacement), generates a permutation of the positions of the alignment. Setting with_replacement to True and otherwise leaving parameters as defaults generates a standard bootstrap resampling of the alignment. """ population_size = len(self) // motif_length if not n: n = population_size if with_replacement: locations = randint(0, population_size, n) else: assert n <= population_size, (n, population_size, motif_length) locations = permutation(population_size)[:n] #check if we need to convert coords for multi-width motifs if motif_length > 1: locations = (locations*motif_length).repeat(motif_length) wrapped_locations =locations.reshape((n,motif_length)) wrapped_locations += arange(motif_length) positions = take(self.ArrayPositions, locations, 0) result = self.__class__(positions.T,force_same_data=True, \ Info=self.Info, Names=self.Names) return result def aln_from_fasta_codons(seqs, array_type=None, Alphabet=None): """Codon alignment from FASTA-format string or lines. This is an InputHandler for taking a FASTA-format string of individual bases and converting it into an array by way of a CodonSequence object that groups triples of bases together and converts them into symbols on the codon alphabet (i.e. each group of 3 bases together is coded by a single symbol). This needs to override the normal aln_from_fasta InputHandler, which asssumes that it can convert the string into the array directly without this grouping step. """ if isinstance(seqs, str): seqs = seqs.split('\n') return aln_from_model_seqs([CodonSequenceGap(s, Label=l) for l, s \ in cogent.parse.fasta.MinimalFastaParser(seqs)]) def xsample(self, n=None, with_replacement=False, motif_length=1, \ random_series=random): """Returns random sample of positions from self, e.g. to bootstrap. Arguments: - n: the number of positions to sample from the alignment. Default is alignment length - with_replacement: boolean flag for determining if sampled positions - random_series: a random number generator with .randint(min,max) .random() methods Notes: By default (resampling all positions without replacement), generates a permutation of the positions of the alignment. Setting with_replacement to True and otherwise leaving parameters as defaults generates a standard bootstrap resampling of the alignment. """ population_size = len(self) // motif_length if not n: n = population_size if with_replacement: locations = [random_series.randint(0, population_size) for samp in xrange(n)] else: assert n <= population_size, (n, population_size, motif_length) locations = random_series.sample(xrange(population_size), n) positions = [(loc*motif_length, (loc+1)*motif_length) for loc in locations] sample = Map(positions, parent_length=len(self)) return self.gappedByMap(sample, Info=self.Info) class CodonDenseAlignment(DenseAlignment): """Stores alignment of gapped codons, no degenerate symbols.""" InputHandlers = { 'array':aln_from_array, 'seqs':aln_from_model_seqs, 'generic':aln_from_generic, 'fasta':aln_from_fasta_codons, 'dense_aln':aln_from_dense_aln, 'aln': aln_from_collection, 'collection':aln_from_collection, 'dict':aln_from_dict, 'empty':aln_from_empty, } def make_gap_filter(template, gap_fraction, gap_run): """Returns f(seq) -> True if no gap runs and acceptable gap fraction. Calculations relative to template. gap_run = number of consecutive gaps allowed in either the template or seq gap_fraction = fraction of positions that either have a gap in the template but not in the seq or in the seq but not in the template NOTE: template and seq must both be ModelSequence objects. """ template_gaps = array(template.gapVector()) def result(seq): """Returns True if seq adhers to the gap threshold and gap fraction.""" seq_gaps = array(seq.gapVector()) #check if gap amount bad if sum(seq_gaps!=template_gaps)/float(len(seq)) > gap_fraction: return False #check if gap runs bad if '\x01'*gap_run in logical_and(seq_gaps, \ logical_not(template_gaps)).astype(uint8).tostring(): return False #check if insertion runs bad elif '\x01'*gap_run in logical_and(template_gaps, \ logical_not(seq_gaps)).astype(uint8).tostring(): return False return True return result class Alignment(_Annotatable, AlignmentI, SequenceCollection): MolType = None #note: this is reset to ASCII in moltype module def __init__(self, *args, **kwargs): """Returns new Alignment object: see SequenceCollection.""" SequenceCollection.__init__(self, *args, **kwargs) #need to convert seqs to Aligned objects seqs = self.SeqData names = self.Names self._motif_probs = {} self._type = self.MolType.gettype() lengths = map(len, self.SeqData) if lengths and (max(lengths) != min(lengths)): raise DataError, "Not all sequences are the same length:\n" + \ "max is %s, min is %s" % (max(lengths), min(lengths)) aligned_seqs = [] for s, n in zip(seqs, names): if isinstance(s, Aligned): s.Name = n #ensure consistency aligned_seqs.append(s) else: aligned_seqs.append(self._seq_to_aligned(s, n)) self.NamedSeqs = self.AlignedSeqs = dict(zip(names, aligned_seqs)) self.SeqData = self._seqs = aligned_seqs def _coerce_seqs(self, seqs, is_array): if not min([isinstance(seq, _Annotatable) or isinstance(seq, Aligned) for seq in seqs]): seqs = map(self.MolType.Sequence, seqs) return seqs def _seq_to_aligned(self, seq, key): """Converts seq to Aligned object -- override in subclasses""" new_seq = self.MolType.Sequence(seq, key) aligned = Aligned(*new_seq.parseOutGaps()) if hasattr(seq, "annotations"): aligned.data.copyAnnotations(seq) return aligned def getTracks(self, policy): # drawing code related # same as sequence but annotations go below sequence tracks return policy.tracksForAlignment(self) def getChildTracks(self, policy): """The only Alignment method required for cogent.draw""" tracks = [] for label in self.Names: seq = self.NamedSeqs[label] tracks += seq.getTracks(policy.copy(seqname=label)) return tracks def __repr__(self): seqs = [] limit = 10 delimiter = '' for (count, name) in enumerate(self.Names): if count == 3: seqs.append('...') break elts = list(self.getGappedSeq(name)[:limit+1]) if len(elts) > limit: elts.append('...') seqs.append("%s[%s]" % (name, delimiter.join(elts))) seqs = ', '.join(seqs) return "%s x %s %s alignment: %s" % (len(self.Names), self.SeqLen, self._type, seqs) def _mapped(self, slicemap): align = [] for name in self.Names: align.append((name, self.NamedSeqs[name][slicemap])) return self.__class__(MolType=self.MolType, data=align) def gappedByMap(self, keep, **kwargs): # keep is a Map seqs = [] for seq_name in self.Names: aligned = self.NamedSeqs[seq_name] seqmap = aligned.map[keep] seq = aligned.data.gappedByMap(seqmap) seqs.append((seq_name, seq)) return self.__class__(MolType=self.MolType, data=seqs, **kwargs) def getProjectedAnnotations(self, seq_name, *args): aligned = self.NamedSeqs[seq_name] seq_annots = self.getAnnotationsMatching(*args) return [a.remappedTo(aligned.data, aligned.map) for a in seq_annots] def getAnnotationsFromSequence(self, seq_name, *args): aligned = self.NamedSeqs[seq_name] return aligned.getAnnotationsMatching(self, *args) def getAnnotationsFromAnySequence(self, *args): result = [] for seq_name in self.Names: result.extend(self.getAnnotationsFromSequence(seq_name, *args)) return result def getBySequenceAnnotation(self, seq_name, *args): result = [] for feature in self.getAnnotationsFromSequence(seq_name, *args): segment = self[feature.map.Start:feature.map.End] segment.Name = '%s "%s" %s to %s of %s' % ( feature.type, feature.Name, feature.map.Start, feature.map.End, self.Name or '') result.append(segment) return result def withMaskedAnnotations(self, annot_types, mask_char=None, shadow=False): """returns an alignment with annot_types regions replaced by mask_char if shadow is False, otherwise all other regions are masked. Arguments: - annot_types: annotation type(s) - mask_char: must be a character valid for the seq MolType. The default value is the most ambiguous character, eg. '?' for DNA - shadow: whether to mask the annotated regions, or everything but the annotated regions""" masked_seqs = [] for seq in self.Seqs: # we mask each sequence using these spans masked_seqs += [seq._masked_annotations(annot_types,mask_char,shadow)] new = self.__class__(data=masked_seqs, Info=self.Info, Name=self.Name) return new def variablePositions(self, include_gap_motif = True): """Return a list of variable position indexes. Arguments: - include_gap_motif: if False, sequences with a gap motif in a column are ignored.""" seqs = [self.getGappedSeq(n) for n in self.Names] seq1 = seqs[0] positions = zip(*seqs[1:]) result = [] for (position, (motif1, column)) in enumerate(zip(seq1,positions)): for motif in column: if motif != motif1: if include_gap_motif: result.append(position) break elif motif != '-' and motif1 != '-': result.append(position) break return result def filtered(self, predicate, motif_length=1, **kwargs): """The alignment positions where predicate(column) is true. Arguments: - predicate: a callback function that takes an tuple of motifs and returns True/False - motif_length: length of the motifs the sequences should be split into, eg. 3 for filtering aligned codons.""" gv = [] kept = False seqs = [self.getGappedSeq(n).getInMotifSize(motif_length, **kwargs) for n in self.Names] positions = zip(*seqs) for (position, column) in enumerate(positions): keep = predicate(column) if kept != keep: gv.append(position*motif_length) kept = keep if kept: gv.append(len(positions)*motif_length) locations = [(gv[i], gv[i+1]) for i in range(0, len(gv), 2)] keep = Map(locations, parent_length=len(self)) return self.gappedByMap(keep, Info=self.Info) def getSeq(self, seqname): """Return a ungapped Sequence object for the specified seqname. Note: always returns Sequence object, not ModelSequence. """ return self.NamedSeqs[seqname].data def getGappedSeq(self, seq_name, recode_gaps=False): """Return a gapped Sequence object for the specified seqname. Note: always returns Sequence object, not ModelSequence. """ return self.NamedSeqs[seq_name].getGappedSeq(recode_gaps) def iterPositions(self, pos_order=None): """Iterates over positions in the alignment, in order. pos_order refers to a list of indices (ints) specifying the column order. This lets you rearrange positions if you want to (e.g. to pull out individual codon positions). Note that self.iterPositions() always returns new objects, by default lists of elements. Use map(f, self.iterPositions) to apply the constructor or function f to the resulting lists (f must take a single list as a parameter). Note that some sequences (e.g. ViennaStructures) have rules that prevent arbitrary strings of their symbols from being valid objects. Will raise IndexError if one of the indices in order exceeds the sequence length. This will always happen on ragged alignments: assign to self.SeqLen to set all sequences to the same length. """ get = self.NamedSeqs.__getitem__ pos_order = pos_order or xrange(self.SeqLen) seq_order = self.Names aligned_objs = [get(seq) for seq in seq_order] seqs = map(str, aligned_objs) for pos in pos_order: yield [seq[pos] for seq in seqs] Positions = property(iterPositions) def withGapsFrom(self, template): """Same alignment but overwritten with the gaps from 'template'""" if len(self) != len(template): raise ValueError("Template alignment must be same length") gap = self.Alphabet.Gap tgp = template.Alphabet.Gap result = {} for name in self.Names: seq = self.getGappedSeq(name) if name not in template.Names: raise ValueError("Template alignment doesn't have a '%s'" % name) gsq = template.getGappedSeq(name) assert len(gsq) == len(seq) combo = [] for (s,g) in zip(seq, gsq): if g == tgp: combo.append(gap) else: combo.append(s) result[name] = combo return Alignment(result, Alphabet=self.Alphabet.withGapMotif())