Package Bio :: Package Sequencing :: Module Phd
[hide private]
[frames] | no frames]

Source Code for Module Bio.Sequencing.Phd

  1  # Copyright 2004 by Cymon J. Cox and Frank Kauff.  All rights reserved. 
  2  # Copyright 2008 by Michiel de Hoon.  All rights reserved. 
  3  # This code is part of the Biopython distribution and governed by its 
  4  # license.  Please see the LICENSE file that should have been included 
  5  # as part of this package. 
  6  """ 
  7  Parser for PHD files output by PHRED and used by PHRAP and CONSED. 
  8   
  9  This module can be used used directly which will return Record objects 
 10  which should contain all the original data in the file. 
 11   
 12  Alternatively, using Bio.SeqIO with the "phd" format will call this module 
 13  internally.  This will give SeqRecord objects for each contig sequence. 
 14  """ 
 15   
 16  from Bio import Seq 
 17  from Bio.Alphabet import generic_dna 
 18   
 19  CKEYWORDS=['CHROMAT_FILE','ABI_THUMBPRINT','PHRED_VERSION','CALL_METHOD',\ 
 20          'QUALITY_LEVELS','TIME','TRACE_ARRAY_MIN_INDEX','TRACE_ARRAY_MAX_INDEX',\ 
 21          'TRIM','TRACE_PEAK_AREA_RATIO','CHEM','DYE'] 
 22   
23 -class Record:
24 """Hold information from a PHD file."""
25 - def __init__(self):
26 self.file_name = '' 27 self.comments={} 28 for kw in CKEYWORDS: 29 self.comments[kw.lower()]=None 30 self.sites = [] 31 self.seq = '' 32 self.seq_trimmed = ''
33 34
35 -def read(handle):
36 """Reads the next PHD record from the file, returning it as a Record object. 37 38 This function reads PHD file data line by line from the handle, 39 and returns a single Record object. 40 """ 41 for line in handle: 42 if line.startswith("BEGIN_SEQUENCE"): 43 record = Record() 44 record.file_name = line[15:].rstrip() 45 break 46 else: 47 return # No record found 48 49 for line in handle: 50 if line.startswith("BEGIN_COMMENT"): 51 break 52 else: 53 raise ValueError("Failed to find BEGIN_COMMENT line") 54 55 for line in handle: 56 line = line.strip() 57 if not line: 58 continue 59 if line=="END_COMMENT": 60 break 61 keyword, value = line.split(":", 1) 62 keyword = keyword.lower() 63 value = value.strip() 64 if keyword in ('chromat_file', 65 'phred_version', 66 'call_method', 67 'chem', 68 'dye', 69 'time', 70 'basecaller_version', 71 'trace_processor_version'): 72 record.comments[keyword] = value 73 elif keyword in ('abi_thumbprint', 74 'quality_levels', 75 'trace_array_min_index', 76 'trace_array_max_index'): 77 record.comments[keyword] = int(value) 78 elif keyword=='trace_peak_area_ratio': 79 record.comments[keyword] = float(value) 80 elif keyword=='trim': 81 first, last, prob = value.split() 82 record.comments[keyword] = (int(first), int(last), float(prob)) 83 else: 84 raise ValueError("Failed to find END_COMMENT line") 85 86 for line in handle: 87 if line.startswith('BEGIN_DNA'): 88 break 89 else: 90 raise ValueError("Failed to find BEGIN_DNA line") 91 92 for line in handle: 93 if line.startswith('END_DNA'): 94 break 95 else: 96 base, quality, location = line.split() 97 record.sites.append((base, quality, location)) 98 99 for line in handle: 100 if line.startswith("END_SEQUENCE"): 101 break 102 else: 103 raise ValueError("Failed to find END_SEQUENCE line") 104 105 record.seq = Seq.Seq(''.join([n[0] for n in record.sites]), generic_dna) 106 if record.comments['trim'] is not None: 107 first, last = record.comments['trim'][:2] 108 record.seq_trimmed = record.seq[first:last] 109 110 return record
111
112 -def parse(handle):
113 """Iterates over a file returning multiple PHD records. 114 115 The data is read line by line from the handle. The handle can be a list 116 of lines, an open file, or similar; the only requirement is that we can 117 iterate over the handle to retrieve lines from it. 118 119 Typical usage: 120 121 records = parse(handle) 122 for record in records: 123 # do something with the record object 124 """ 125 while True: 126 record = read(handle) 127 if not record: 128 return 129 yield record
130 131 132 # ---------- Everything below is deprecated 133 134 from Bio import File 135 from Bio.ParserSupport import * 136 137
138 -class Iterator:
139 """Iterates over a file of multiple PHD records (DEPRECATED). 140 141 Methods: 142 next Return the next record from the stream, or None. 143 """ 144
145 - def __init__(self, handle, parser=None):
146 """Create a new iterator. 147 148 Create a new iterator. handle is a file-like object. parser 149 is an optional Parser object to change the results into another form. 150 If set to None, then the raw contents of the file will be returned. 151 """ 152 import warnings 153 warnings.warn("Bio.Sequencing.Phd.Iterator is deprecated. Please use Bio.Sequencing.Phd.parse(handle) instead of Bio.Sequencing.Phd.Iterator(handle, RecordParser())", DeprecationWarning) 154 self._uhandle = File.UndoHandle(handle) 155 self._parser = parser
156
157 - def next(self):
158 lines = [] 159 while 1: 160 line = self._uhandle.readline() 161 if not line: 162 break 163 # If a new record, then put the line back and stop. 164 if lines and line[:14] == 'BEGIN_SEQUENCE': 165 self._uhandle.saveline(line) 166 break 167 lines.append(line) 168 169 if not lines: 170 return None 171 172 data = ''.join(lines) 173 if self._parser is not None: 174 return self._parser.parse(File.StringHandle(data)) 175 return data
176
177 - def __iter__(self):
178 """Iterate over the PHY file record by record.""" 179 return iter(self.next, None)
180
181 -class RecordParser(AbstractParser):
182 """Parses PHD file data into a Record object (DEPRECATED)."""
183 - def __init__(self):
184 import warnings 185 warnings.warn("Bio.Sequencing.Phd.RecordParser is deprecated. Please use Bio.Sequencing.Phd.read(handle) instead of Bio.Sequencing.Phd.RecordParser().parse(handle)", DeprecationWarning) 186 self._scanner = _Scanner() 187 self._consumer = _RecordConsumer()
188
189 - def parse(self, handle):
190 if isinstance(handle, File.UndoHandle): 191 uhandle = handle 192 else: 193 uhandle = File.UndoHandle(handle) 194 self._scanner.feed(uhandle, self._consumer) 195 return self._consumer.data
196 197
198 -class _Scanner:
199 """Scans a PHD-formatted file (DEPRECATED). 200 201 Methods: 202 feed - Feed one PHD record. 203 """
204 - def feed(self, handle, consumer):
205 """Reads in PDH data from the handle for scanning. 206 207 Feed in PHD data for scanning. handle is a file-like object 208 containing PHD data. consumer is a Consumer object that will 209 receive events as the PHD data is scanned. 210 """ 211 assert isinstance(handle, File.UndoHandle), \ 212 "handle must be an UndoHandle" 213 if handle.peekline(): 214 self._scan_record(handle, consumer)
215
216 - def _scan_record(self, uhandle, consumer):
217 self._scan_begin_sequence(uhandle, consumer) 218 self._scan_comments(uhandle, consumer) 219 self._scan_dna(uhandle, consumer) 220 consumer.end_sequence()
221
222 - def _scan_begin_sequence(self, uhandle, consumer):
223 read_and_call(uhandle, consumer.begin_sequence, start = 'BEGIN_SEQUENCE')
224
225 - def _scan_comments(self, uhandle, consumer):
226 227 read_and_call_while(uhandle, consumer.noevent, blank=1) 228 read_and_call(uhandle, consumer.noevent, start = 'BEGIN_COMMENT') 229 read_and_call_while(uhandle, consumer.noevent, blank=1) 230 231 while 1: 232 for kw in CKEYWORDS: 233 if attempt_read_and_call(uhandle,getattr(consumer,kw.lower()),start=kw+':'): 234 break # recognized keyword: end for loop and do another while 235 else: 236 break # no keywords found: end while loop 237 238 read_and_call_while(uhandle, consumer.noevent, blank=1) 239 read_and_call(uhandle, consumer.noevent, start = 'END_COMMENT')
240
241 - def _scan_dna(self, uhandle, consumer):
242 while 1: 243 line = uhandle.readline() 244 if is_blank_line(line) or line == 'BEGIN_DNA\n': 245 continue 246 elif line == 'END_DNA\n': 247 break 248 consumer.read_dna(line)
249 250
251 -class _RecordConsumer(AbstractConsumer):
252 """Consumer that converts a PHD record to a Record object (DEPRECATED)."""
253 - def __init__(self):
254 self.data = None
255
256 - def begin_sequence(self, line):
257 self.data = Record() 258 self.data.file_name = line[15:].rstrip()
259
260 - def end_sequence(self):
261 self.data.seq = Seq.Seq(''.join([n[0] for n in self.data.sites]), generic_dna) 262 if self.data.comments['trim'] is not None : 263 first = self.data.comments['trim'][0] 264 last = self.data.comments['trim'][1] 265 self.data.seq_trimmed = self.data.seq[first:last]
266
267 - def chromat_file(self, line):
268 self.data.comments['chromat_file'] = line[13:-1].strip()
269
270 - def abi_thumbprint(self, line):
271 self.data.comments['abi_thumbprint'] = int(line[15:-1].strip())
272
273 - def phred_version(self, line):
274 self.data.comments['phred_version'] = line[14:-1].strip()
275
276 - def call_method(self, line):
277 self.data.comments['call_method'] = line[12:-1].strip()
278
279 - def quality_levels(self, line):
280 self.data.comments['quality_levels'] = int(line[15:-1].strip())
281
282 - def time(self, line):
283 self.data.comments['time'] = line[5:-1].strip()
284
285 - def trace_array_min_index(self, line):
286 self.data.comments['trace_array_min_index'] = int(line[22:-1].strip())
287
288 - def trace_array_max_index(self, line):
289 self.data.comments['trace_array_max_index'] = int(line[22:-1].strip())
290
291 - def trim(self, line):
292 first, last, prob = line[5:-1].split() 293 self.data.comments['trim'] = (int(first), int(last), float(prob))
294
295 - def trace_peak_area_ratio(self, line):
296 self.data.comments['trace_peak_area_ratio'] = float(line[22:-1].strip())
297
298 - def chem(self, line):
299 self.data.comments['chem'] = line[5:-1].strip()
300
301 - def dye(self, line):
302 self.data.comments['dye'] = line[4:-1].strip()
303
304 - def read_dna(self, line):
305 base, quality, location = line.split() 306 self.data.sites.append((base, quality, location))
307 308 if __name__ == "__main__" : 309 print "Quick self test" 310 #Test the iterator, 311 handle = open("../../Tests/Phd/phd1") 312 records = parse(handle) 313 for record in records: 314 print record.file_name, len(record.seq) 315 handle.close() 316 print "Done" 317