Mercurial > repos > davidmurphy > codonlogo
diff corebio/seq_io/fasta_io.py @ 0:c55bdc2fb9fa
Uploaded
author | davidmurphy |
---|---|
date | Thu, 27 Oct 2011 12:09:09 -0400 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/corebio/seq_io/fasta_io.py Thu Oct 27 12:09:09 2011 -0400 @@ -0,0 +1,230 @@ +#!/usr/bin/env python + +# Copyright (c) 2005 Gavin E. Crooks <gec@threeplusone.com> +# +# This software is distributed under the MIT Open Source License. +# <http://www.opensource.org/licenses/mit-license.html> +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the "Software"), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# + +"""Read and write sequence information in FASTA format. + +This is a very common format for unannotated biological sequence data, +accepted by many multiple sequence alignment programs. Each sequence +consists of a single-line description, followed by lines of sequence data. +The first character of the description line is a greater-than (">") symbol +in the first column. The first word of the description is often the name or +ID of the sequence. Fasta files containing multiple sequences have one +sequence listed right after another. + + +Example Fasta File :: + +>Lamprey GLOBIN V - SEA LAMPREY +PIVDTGSVA-P------------------LSAAEKTKIRSAWAPVYSTY---ETSGVDILVKFFTSTPAAQEFFPKFKGL +TT-----ADQLKKSA---DVRWHA-ERIINAVNDAVASMDDTEKMS--MKL-RDLSGKH----AKSFQV-----DPQYFK +VLAAVI-AD-TVAAGD--AGFEKLMSM------I---CILLR----S-----A-----Y------------ +>Hagfish GLOBIN III - ATLANTIC HAGFISH +PITDHGQPP-T------------------LSEGDKKAIRESWPQIYKNF---EQNSLAVLLEFLKKFPKAQDSFPKFSAK +KS-------HLEQDP---AVKLQA-EVIINAVNHTIGLMDKEAAMK--KYL-KDLSTKH----STEFQV-----NPDMFK +ELSAVF-VS-TMG-GK--AAYEKLFSI------I---ATLLR----S-----T-----YDA---------- +>Frog HEMOGLOBIN BETA CHAIN - EDIBLE FROG +----------GS-----------------------DLVSGFWGKV--DA---HKIGGEALARLLVVYPWTQRYFTTFGNL +GSADAIC-----HNA---KVLAHG-EKVLAAIGEGLKHPENLKAHY--AKL-SEYHSNK----LHVDPANFRLLGNVFIT +VLARHF-QH-EFTPELQ-HALEAHFCA------V---GDALA----K-----A-----YH----------- + + +""" +import re +from corebio.utils import * +from corebio.seq import * +from corebio.seq_io import * + + +names = ( 'fasta', 'pearson', 'fa') +extensions = ('fa', 'fasta', 'fast', 'seq', 'fsa', 'fst', 'nt', 'aa','fna','mpfa', 'faa', 'fnn','mfasta') + + +example = """ +>Lamprey GLOBIN V - SEA LAMPREY +PIVDTGSVA-P------------------LSAAEKTKIRSAWAPVYSTY---ETSGVDILVKFFTSTPAAQEFFPKFKGL +TT-----ADQLKKSA---DVRWHA-ERIINAVNDAVASMDDTEKMS--MKL-RDLSGKH----AKSFQV-----DPQYFK +VLAAVI-AD-TVAAGD--AGFEKLMSM------I---CILLR----S-----A-----Y------------ + +>Hagfish GLOBIN III - ATLANTIC HAGFISH +PITDHGQPP-T------------------LSEGDKKAIRESWPQIYKNF---EQNSLAVLLEFLKKFPKAQDSFPKFSAK +KS-------HLEQDP---AVKLQA-EVIINAVNHTIGLMDKEAAMK--KYL-KDLSTKH----STEFQV-----NPDMFK +ELSAVF-VS-TMG-GK--AAYEKLFSI------I---ATLLR----S-----T-----YDA---------- + +>Frog HEMOGLOBIN BETA CHAIN - EDIBLE FROG +----------GS-----------------------DLVSGFWGKV--DA---HKIGGEALARLLVVYPWTQRYFTTFGNL +GSADAIC-----HNA---KVLAHG-EKVLAAIGEGLKHPENLKAHY--AKL-SEYHSNK----LHVDPANFRLLGNVFIT +VLARHF-QH-EFTPELQ-HALEAHFCA------V---GDALA----K-----A-----YH----------- + +""" + + +def read(fin, alphabet=None): + """Read and parse a fasta file. + + Args: + fin -- A stream or file to read + alphabet -- The expected alphabet of the data, if given + Returns: + SeqList -- A list of sequences + Raises: + ValueError -- If the file is unparsable + """ + seqs = [ s for s in iterseq(fin, alphabet)] + name = names[0] + if hasattr(fin, "name") : name = fin.name + return SeqList(seqs, name=name) + + +def readseq(fin, alphabet=None) : + """Read one sequence from the file, starting + from the current file position.""" + return iterseq(fin, alphabet).next() + + +def iterseq(fin, alphabet=None): + """ Parse a fasta file and generate sequences. + + Args: + fin -- A stream or file to read + alphabet -- The expected alphabet of the data, if given + Yeilds: + Seq -- One alphabetic sequence at a time. + Raises: + ValueError -- If the file is unparsable + """ + alphabet = Alphabet(alphabet) + + seqs = [] + comments = [] # FIXME: comments before first sequence are lost. + header = None + header_lineno = -1 + + def build_seq(seqs,alphabet, header, header_lineno,comments) : + try : + name = header.split(' ',1)[0] + if comments : + header += '\n' + '\n'.join(comments) + s = Seq( "".join(seqs), alphabet, name=name, description=header) + except ValueError: + raise ValueError( + "Parsed failed with sequence starting at line %d: " + "Character not in alphabet: %s" % (header_lineno, alphabet) ) + return s + + for lineno, line in enumerate(fin) : + line = line.strip() + if line == '' : continue + if line.startswith('>') : + if header is not None : + yield build_seq(seqs,alphabet, header, header_lineno, comments) + header = None + seqs = [] + header = line[1:] + header_lineno = lineno + comments = [] + elif line.startswith(';') : + # Optional (and unusual) comment line + comments.append(line[1:]) + else : + if header is None : + raise ValueError ( + "Parse failed on line %d: sequence before header" + % (lineno) ) + seqs.append(line) + + if not seqs: return + yield build_seq(seqs,alphabet, header, header_lineno, comments) + + +def write(fout, seqs): + """Write a fasta file. + + Args: + fout -- A writable stream. + seqs -- A list of Seq's + """ + if seqs.description : + for line in seqs.description.splitlines(): + print >>fout, ';'+ line + for s in seqs : + writeseq(fout, s) + + +def writeseq(afile, seq): + """ Write a single sequence in fasta format. + + Args: + afile -- A writable stream. + seq -- A Seq instance + """ + + header = seq.description or seq.name or '' + + # We prepend '>' to the first header line + # Additional lines start with ';' to indicate comment lines + if header : + header = header.splitlines() + print >>afile, '>'+header[0] + if len(header) > 1 : + for h in header[1:] : + print >>afile, ';' +h + else : + print >>afile, '>' + + L = len(seq) + line_length = 80 + for n in range (1+ L/line_length) : + print >>afile, seq[n * line_length: (n+1) * line_length] + print >>afile + + +def index(afile, alphabet=None) : + """Return a FileIndex for the fasta file. Sequences can be retrieved + by item number or name. + """ + def parser( afile) : + return readseq(afile, alphabet) + + key = re.compile(r"^>\s*(\S*)") + def linekey( line): + k = key.search(line) + if k is None : return None + return k.group(1) + + return FileIndex(afile, linekey, parser) + + + + + + + + + + + + + \ No newline at end of file