Документ взят из кэша поисковой машины. Адрес оригинального документа : http://kodomo.fbb.msu.ru/hg/allpy/raw-rev/b8ddd35fac65 Дата изменения: Unknown Дата индексирования: Tue Oct 2 08:01:14 2012 Кодировка: Поисковые слова: ccd

# HG changeset patch
# User Daniil Alexeyevsky
# Date 1320788965 -10800
# Node ID b8ddd35fac65cd22dc3d7fb92b0ea05bbf72371d
# Parent 0615fe6265cd87522266c6f07f2900380201ad2f
Added internal support for new file format: FastQ

diff -r 0615fe6265cd -r b8ddd35fac65 allpy/fileio/auto.py
--- a/allpy/fileio/auto.py Tue Nov 08 19:12:37 2011 +0300
+++ b/allpy/fileio/auto.py Wed Nov 09 00:49:25 2011 +0300
@@ -1,4 +1,5 @@
import fasta
+import fastq
import markup
import bio
import emboss
@@ -8,6 +9,8 @@
def __new__(cls, file, format="fasta", **kw):
if format == "fasta":
return fasta.FastaFile(file, **kw)
+ elif format == "fastq":
+ return fastq.FastqFile(file, **kw)
elif format == 'markup':
return markup.MarkupFile(file, **kw)
elif format.startswith('markup:'):
diff -r 0615fe6265cd -r b8ddd35fac65 allpy/fileio/fastq.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/allpy/fileio/fastq.py Wed Nov 09 00:49:25 2011 +0300
@@ -0,0 +1,47 @@
+class FastqFile(object):
+ """FastQ parser & writer."""
+
+ def __init__(self, file, format='fasta', gaps=None):
+ # NOTE: 'gaps' argument is ignored!
+ self.file = file
+ self.format = format
+ self.phred_base = 33
+ try:
+ self.source = file.name
+ except Exception:
+ self.source = ''
+
+ def write_alignment(self, alignment):
+ """Append alignment to the file."""
+ write = self.file.write
+ for row in alignment.rows_as_lists():
+ header = '%s %s' % (row.sequence.name, row.sequence.description)
+ chars = "".join(monomer.code1 for monomer in row)
+ phred = "".join(
+ chr(monomer.quality + self.phred_base)
+ for monomer in row
+ )
+ write('@%s\n%s\n+\n%s\n' % (header, row, phred))
+
+ def read_alignment(self, alignment):
+ """Append file contents to the alignment."""
+ for n, line in enumerate(self.file):
+ n = n % 4
+ line = line.rstrip('\n')
+ if n == 0:
+ assert line[0] == '@'
+ name, description = line.split(' ', 1)
+ elif n == 1:
+ # XXX: this is ugly!
+ # XXX append_row_from_string should return the new Sequence
+ alignment.append_row_from_string(line, name, description, self.source, [])
+ sequence = alignment.sequences[-1]
+ elif n == 2:
+ assert line[0] == '+'
+ elif n == 3:
+ assert len(sequence) == len(line)
+ quality = sequence.add_markup('quality')
+ for monomer, phred in zip(sequence, line):
+ monomer.quality = ord(phred) - self.phred_base
+
+# vim: set et ts=4 sts=4 sw=4:
diff -r 0615fe6265cd -r b8ddd35fac65 allpy/markups.py
--- a/allpy/markups.py Tue Nov 08 19:12:37 2011 +0300
+++ b/allpy/markups.py Wed Nov 09 00:49:25 2011 +0300
@@ -152,6 +152,9 @@
return monomer.code1.lower()
raise AssertionError("Unknown monomer case")

+class SequenceQualityMarkup(base.SequenceMarkup, IntMarkupMixin):
+ name = 'quality'
+
class SequencePdbResiMarkup(base.SequenceMarkup, IntMarkupMixin):
name = 'pdb_resi'

diff -r 0615fe6265cd -r b8ddd35fac65 test/test_fastq.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/test/test_fastq.py Wed Nov 09 00:49:25 2011 +0300
@@ -0,0 +1,27 @@
+from StringIO import StringIO
+from allpy import dna
+
+example = """
+@HWI-ST992:108:D06T0ACXX:6:1107:13562:94150 1:N:0:
+TCAATGGTGCTGTATAAAACACGTAGAATTACTATTATGGATTATATTCGTCTTCATTTGTTAAAATCAGTGGCTTAACTGGCTAAATTTTGTGGACTATA
++
+#8::4AD,A+CEDEIEDFIE:FF;(>AAAA>
+@HWI-ST992:108:D06T0ACXX:6:1107:13630:94150 1:N:0:
+GGTGAAGCGCTTTGAGTCCCTCTCCAGCCAGCGCTTCCTGCAGGGCGTGATCGTGGCGCTCATCACGGGCATGTTCTGGTGGCAGCGCGGCCGCGGCGGCT
++
+<@@DDFDBHHGFHIAHHIIIIIIIIIGHGIGAGHIHIIIIHIIII?G@BHHHHABDDDDDDDDDDDD?B@BCD@C@CDDACDDDDDDDDDBDBDBDD@B<>
+""".strip()
+
+def test_fastq_io():
+ file = StringIO(example)
+ aln = dna.Alignment().append_file(file, format="fastq")
+ assert len(aln.sequences) == 2
+ a, b = aln.sequences
+ assert str(a)[:10] == "TCAATGGTGC"
+ assert str(b)[:10] == "GGTGAAGCGC"
+ qa = list(a.markups['quality'].sorted_values())
+ assert qa[:10] == [2, 23, 25, 25, 19, 32, 35, 11, 32, 10]
+ qb = list(b.markups['quality'].sorted_values())
+ assert qb[:10] == [27, 31, 31, 35, 35, 37, 35, 33, 39, 39]
+
+# vim: set et ts=4 sts=4 sw=4: