Context Navigation

lav.py

リビジョン 3, 17.8 KB (コミッタ: kohda, 14 年前)
Install Unix tools http://hannonlab.cshl.edu/galaxy_unix_tools/galaxy.html

行番号
1	"""
2	Support for reading and writing the LAV format produced by the `blastz`_
3	pairwise aligner.
4
5	.. _blastz: http://www.bx.psu.edu/miller_lab/
6	"""
7
8	from bx.align import *
9	import bx.seq
10
11	import sys,math,StringIO
12	import itertools
13	from bx import interval_index_file
14
15	class Reader(object):
16	"""Iterate over all lav blocks in a file in order"""
17
18	def __init__(self,file,path_subs=None,fail_to_ns=False):
19	self.file = file
20	self.lineNumber = 0
21	self.path_subs = path_subs # list of (prefix,replacement) to allow
22	if (self.path_subs == None): # .. redirection of sequence file paths
23	self.path_subs = [] # .. on different machines
24	self.fail_to_ns = fail_to_ns # True => if sequences fail to open,
25	# create a fake file of all Ns
26
27	self.d_stanza_text = None
28
29	self.seq1_filename = None
30	self.seq1_file = None
31	self.seq1_header = None
32	self.seq1_start = None
33	self.seq1_end = None
34	self.seq1_strand = None
35	self.seq1_contig = None
36	self.seq1_src = None
37	self.seq1_gap = None
38
39	self.seq2_filename = None
40	self.seq2_file = None
41	self.seq2_header = None
42	self.seq2_start = None
43	self.seq2_end = None
44	self.seq2_strand = None
45	self.seq2_contig = None
46	self.seq2_src = None
47	self.seq2_gap = None
48
49	def next(self):
50	while (True):
51	line = self.fetch_line(strip=None,requireLine=False)
52	assert (line), "unexpected end of file (missing #:eof)"
53	line = line.rstrip()
54	if (line == ""): # (allow blank lines between stanzas)
55	continue
56	if (line == "#:eof"):
57	line = self.file.readline().rstrip()
58	assert (not line), "extra line after #:eof (line %d, \"%s\")" \
59	% (self.lineNumber,line)
60	return None
61	if (line == "#:lav"):
62	continue
63	if (line.startswith("d {")):
64	self.d_stanza_text = self.parse_unknown_stanza()
65	continue
66	if (line.startswith("s {")):
67	self.parse_s_stanza()
68	continue
69	if (line.startswith("h {")):
70	self.parse_h_stanza()
71	continue
72	if (line.startswith("a {")):
73	(score,pieces) = self.parse_a_stanza()
74	break
75	if (line.endswith("{")):
76	self.parse_unknown_stanza()
77	continue
78	assert (False), "incomprehensible line (line %d, \"%s\")" \
79	% (self.lineNumber,line)
80	return self.build_alignment(score,pieces)
81
82	def __iter__(self):
83	return ReaderIter(self)
84
85	def close(self):
86	self.file.close()
87
88	def open_seqs(self):
89	if (self.seq1_file != None) and (self.seq2_file != None):
90	return
91
92	if (self.seq1_file == None):
93	if (self.seq1_strand == "+"): revcomp = False
94	else: revcomp = "-5'"
95	if (self.seq1_contig == 1): contig = None
96	else: contig = self.seq1_contig
97	try:
98	f = file(self.seq1_filename,"rb")
99	except:
100	if (self.fail_to_ns):
101	f = StringIO.StringIO(">seq1\n" + ("n" * (self.seq1_end - self.seq1_start)))
102	revcomp = False
103	contig = 1
104	else:
105	assert (False), "failed to open %s" % self.seq1_filename
106	self.seq1_file = bx.seq.seq_file(f,revcomp=revcomp,contig=contig)
107	self.seq1_gap = self.seq1_file.gap
108	try:
109	name1 = self.header_to_src_name(self.seq1_header)
110	except ValueError:
111	try:
112	name1 = self.path_to_src_name(self.seq1_filename)
113	except ValueError:
114	name1 = "seq1"
115	(species1,chrom1) = src_split(name1)
116	self.seq1_src = src_merge(species1,chrom1,contig)
117	if (contig != None): chrom1 += "[%s]" % contig
118
119	if (self.seq2_file == None):
120	if (self.seq2_strand == "+"): revcomp = False
121	else: revcomp = "-5'"
122	if (self.seq2_contig == 1): contig = None
123	else: contig = self.seq2_contig
124	try:
125	f = file(self.seq2_filename,"rb")
126	except:
127	if (self.fail_to_ns):
128	f = StringIO.StringIO(">seq2\n" + ("n" * (self.seq2_end - self.seq2_start)))
129	revcomp = False
130	contig = 1
131	else:
132	assert (False), "failed to open %s" % self.seq1_filename
133	self.seq2_file = bx.seq.seq_file(f,revcomp=revcomp,contig=contig)
134	self.seq2_gap = self.seq2_file.gap
135	try:
136	name2 = self.header_to_src_name(self.seq2_header)
137	except ValueError:
138	try:
139	name2 = self.path_to_src_name(self.seq2_filename)
140	except ValueError:
141	name2 = "seq2"
142	(species2,chrom2) = src_split(name2)
143	self.seq2_src = src_merge(species2,chrom2,contig)
144	if (contig != None): chrom2 += "[%s]" % contig
145
146	length1 = self.seq1_file.length
147	length2 = self.seq2_file.length
148	assert (species1 != species2) or (chrom1 != chrom2) or (length1 == length2), \
149	"conflicting lengths for %s (%d and %d)" % (self.seq1_src,length1,length2)
150
151	self.species_to_lengths = {}
152	self.species_to_lengths[species1] = {}
153	self.species_to_lengths[species2] = {} # (OK if it clobbers line above)
154	self.species_to_lengths[species1][chrom1] = self.seq1_file.length
155	self.species_to_lengths[species2][chrom2] = self.seq2_file.length
156
157	def close_seqs(self):
158	if (self.seq1_file != None):
159	self.seq1_file.close()
160	self.seq1_file = None
161	if (self.seq2_file != None):
162	self.seq2_file.close()
163	self.seq2_file = None
164
165	def parse_s_stanza(self):
166	self.close_seqs()
167	line = self.fetch_line(report=" in s-stanza")
168	(self.seq1_filename,
169	self.seq1_start,
170	self.seq1_end,
171	self.seq1_strand,
172	self.seq1_contig) = self.parse_s_seq(line)
173
174	line = self.fetch_line(report=" in s-stanza")
175	(self.seq2_filename,
176	self.seq2_start,
177	self.seq2_end,
178	self.seq2_strand,
179	self.seq2_contig) = self.parse_s_seq(line)
180
181	line = self.fetch_line(report=" in s-stanza")
182	assert (line == "}"), "improper s-stanza terminator (line %d, \"%s\")" \
183	% (self.lineNumber,line)
184
185	def parse_s_seq(self,line):
186	fields = line.split()
187	filename = fields[0].strip('"')
188	start = int(fields[1]) - 1
189	end = int(fields[2])
190	contig = int(fields[4])
191	if (fields[3] == "1"): strand = "-"
192	else: strand = "+"
193	if (filename.endswith("-")):
194	assert (strand == "-"), "strand mismatch in \"%s\"" % line
195	filename = filename[:-1]
196	filename = do_path_subs(filename,self.path_subs)
197	return (filename,start,end,strand,contig)
198
199
200	def parse_h_stanza(self):
201	line = self.fetch_line(strip='"',report=" in h-stanza")
202	self.seq1_header = line
203	self.seq1_header_prefix = ""
204	if (line.startswith(">")):
205	self.seq1_header = line[1:].strip()
206	self.seq1_header_prefix = ">"
207	self.seq1_header = self.seq1_header.split(None,1)
208	if (len(self.seq1_header) > 0): self.seq1_header = self.seq1_header[0]
209	else: self.seq1_header = "seq1"
210
211	line = self.fetch_line(strip='"',report=" in h-stanza")
212	self.seq2_header = line
213	self.seq2_header_prefix = ""
214	if (line.startswith(">")):
215	self.seq2_header = line[1:].strip()
216	self.seq2_header_prefix = ">"
217	self.seq2_header = self.seq2_header.split(None,1)
218	if (len(self.seq2_header) > 0): self.seq2_header = self.seq2_header[0]
219	else: self.seq2_header = "seq2"
220
221	line = self.fetch_line(report=" in h-stanza")
222	assert (line == "}"), "improper h-stanza terminator (line %d, \"%s\")" \
223	% (self.lineNumber,line)
224
225
226	def parse_a_stanza(self):
227	"""returns the pair (score,pieces)
228	where pieces is a list of ungapped segments (start1,start2,length,pctId)
229	with start1,start2 origin-0"""
230	# 's' line -- score, 1 field
231	line = self.fetch_line(report=" in a-stanza")
232	fields = line.split()
233	assert (fields[0] == "s"), "s line expected in a-stanza (line %d, \"%s\")" \
234	% (self.lineNumber,line)
235	try: score = int(fields[1])
236	except: score = float(fields[1])
237
238	# 'b' line -- begin positions in seqs, 2 fields
239	line = self.fetch_line(report=" in a-stanza")
240	fields = line.split()
241	assert (fields[0] == "b"), "b line expected in a-stanza (line %d, \"%s\")" \
242	% (self.lineNumber,line)
243	beg1 = int(fields[1]) - 1
244	beg2 = int(fields[2]) - 1
245
246	# 'e' line -- end positions in seqs, 2 fields
247	line = self.fetch_line(report=" in a-stanza")
248	fields = line.split()
249	assert (fields[0] == "e"), "e line expected in a-stanza (line %d, \"%s\")" \
250	% (self.lineNumber,line)
251	len1 = int(fields[1]) - beg1
252	len2 = int(fields[2]) - beg2
253
254	# 'l' lines
255	pieces = []
256	while (True):
257	line = self.fetch_line(report=" in a-stanza")
258	fields = line.split()
259	if (fields[0] != "l"):
260	break
261	start1 = int(fields[1]) - 1
262	start2 = int(fields[2]) - 1
263	length = int(fields[3]) - start1
264	length2 = int(fields[4]) - start2
265	try: pctId = int(fields[5])
266	except: pctId = float(fields[5])
267	assert (length2 == length), "length mismatch in a-stanza"
268	pieces.append((start1+self.seq1_start,start2+self.seq2_start,length,pctId))
269	assert (line == "}"), "improper a-stanza terminator (line %d, \"%s\")" \
270	% (self.lineNumber,line)
271	return (score,pieces)
272
273	def parse_unknown_stanza(self):
274	lines = []
275	while (True):
276	line = self.fetch_line()
277	assert (line), "unexpected end of file (missing #:eof)"
278	if (line == "}"): break
279	lines.append(line)
280	return " " + "\n ".join(lines) + "\n"
281
282	def fetch_line(self,strip=True,requireLine=True,report=""):
283	if (strip == None): line = self.file.readline()
284	elif (strip == True): line = self.file.readline().strip()
285	else: line = self.file.readline().strip().strip(strip)
286	self.lineNumber += 1
287	if (requireLine):
288	assert (line), \
289	"unexpected blank line or end of file%s (line %d)" \
290	% (report,self.lineNumber)
291	return line
292
293
294	def d_stanza(self):
295	if (self.d_stanza_text == None): return ""
296	return "d {\n%s}" % self.d_stanza_text
297
298	def s_stanza(self):
299	if (self.seq1_filename == None): return ""
300
301	if (self.seq1_strand == "-"): seq1_strand = "1"
302	else: seq1_strand = "0"
303	if (self.seq2_strand == "-"): seq2_strand = "1"
304	else: seq2_strand = "0"
305
306	s = " \"%s\" %d %d %s %d\n"\
307	% (self.seq1_filename,self.seq2_start+1,self.seq1_end,
308	seq1_strand,self.seq1_contig)
309	s += " \"%s\" %d %d %s %d\n"\
310	% (self.seq2_filename,self.seq2_start+1,self.seq2_end,
311	seq2_strand,self.seq2_contig)
312
313	return "s {\n%s}" % s
314
315	def h_stanza(self):
316	if (self.seq1_header == None): return ""
317	s = " \"%s%s\"\n" % (self.seq1_header_prefix,self.seq1_header)
318	s += " \"%s%s\"\n" % (self.seq2_header_prefix,self.seq2_header)
319	return "h {\n%s}" % s
320
321	def build_alignment(self,score,pieces):
322	"""converts a score and pieces to an alignment"""
323	# build text
324	self.open_seqs()
325	text1 = text2 = ""
326	end1 = end2 = None
327	for (start1,start2,length,pctId) in pieces:
328	if (end1 != None):
329	if (start1 == end1): # insertion in sequence 2
330	text1 += self.seq1_gap * (start2-end2)
331	text2 += self.seq2_file.get(end2,start2-end2)
332	else: # insertion in sequence 1
333	text1 += self.seq1_file.get(end1,start1-end1)
334	text2 += self.seq2_gap * (start1-end1)
335
336	text1 += self.seq1_file.get(start1,length)
337	text2 += self.seq2_file.get(start2,length)
338	end1 = start1 + length
339	end2 = start2 + length
340	# create alignment
341	start1 = pieces[0][0]
342	start2 = pieces[0][1]
343	end1 = pieces[-1][0] + pieces[-1][2]
344	end2 = pieces[-1][1] + pieces[-1][2]
345	size1 = end1 - start1
346	size2 = end2 - start2
347	a = Alignment(score=score,species_to_lengths=self.species_to_lengths)
348	#if (self.seq1_strand == "-"): start1 = self.seq1_file.length - end1
349	a.add_component(Component(self.seq1_src,start1,size1,self.seq1_strand,text=text1))
350	#if (self.seq2_strand == "-"): start2 = self.seq2_file.length - end2
351	a.add_component(Component(self.seq2_src,start2,size2,self.seq2_strand,text=text2))
352	return a
353
354	def path_to_src_name(self,path_name):
355	# converts, e.g. ".../hg18/seq/chr13.nib" to "hg18.chr13"
356	if (path_name == None) or (path_name == ""): raise ValueError
357	if (path_name.endswith(".nib")): path_name = path_name[:-4]
358	if (path_name.endswith(".fa")): path_name = path_name[:-3]
359	if (path_name.endswith(".fasta")): path_name = path_name[:-6]
360	slash = path_name.rfind("/")
361	if (slash == -1): return path_name
362	name = path_name[slash+1:]
363	path_name = path_name[:slash]
364	if (path_name.endswith("/seq")): path_name = path_name[:-4]
365	slash = path_name.rfind("/")
366	if (slash != -1): path_name = path_name[slash+1:]
367	return path_name + "." + name
368
369	def header_to_src_name(self,header):
370	# converts, e.g. "hg18.chr13:115404472-117281897" to "hg18.chr13"
371	if (header == None) or (header == ""): raise ValueError
372	colon = header.rfind(":")
373	if (colon != -1): header = header[:colon]
374	if ("/" in header): raise ValueError
375	if (header.count(".") == 0):
376	return header
377	header = header.split(".")
378	if (header[0] == "") or (header[1] == ""): raise ValueError
379	return ".".join(header)
380
381	class ReaderIter(object):
382	def __init__(self,reader):
383	self.reader = reader
384	def __iter__(self):
385	return self
386	def next(self):
387	v = self.reader.next()
388	if (not v): raise StopIteration
389	return v
390
391
392	class LavAsPiecesReader(Reader):
393	"""Iterate over all lav blocks in a file in order, returning alignments
394	as score and pieces, as returned by Reader.parse_a_stanza"""
395
396	def build_alignment(self,score,pieces):
397	return (score,pieces)
398
399
400	class Writer(object):
401
402	# blockHash is a hash from (src1,strand1,src2,strand2) to a list of blocks;
403	# the blocks are collected on each call to write(), but the actual writing
404	# does not occur until close().
405
406	def __init__(self,file,attributes={}):
407	self.file = file
408	self.fname1 = None
409	self.fname2 = None
410	self.block = 0
411	self.blockHash = {} # (see note above)
412
413	if ("name_format_1" in attributes):
414	self.fname1 = attributes["name_format_1"]
415	if ("name_format_2" in attributes):
416	self.fname2 = attributes["name_format_2"]
417
418	if ("d_stanza" in attributes):
419	write_lav_marker(self)
420	print >>self.file,"d {"
421	print >>self.file,attributes["d_stanza"]
422	print >>self.file,"}"
423
424	def write(self,alignment):
425	if (len(alignment.components) != 2):
426	raise "%d-component alignment is not compatible with lav" % \
427	len(alignment.components)
428
429	c1 = alignment.components[0]
430	c2 = alignment.components[1]
431
432	key = (c1.src,c1.strand,c2.src,c2.strand)
433	if (key not in self.blockHash): self.blockHash[key] = []
434	self.blockHash[key].append(alignment)
435	self.block += 1
436
437	def close(self):
438	keys = [key for key in self.blockHash]
439	keys = sort_keys_by_chrom (keys)
440	for key in keys:
441	(src1,strand1,src2,strand2) = key
442	alignment = self.blockHash[key][0]
443	self.src1 = src1
444	self.strand1 = strand1
445	self.length1 = alignment.src_size(src1)
446	self.src2 = src2
447	self.strand2 = strand2
448	self.length2 = alignment.src_size(src2)
449	self.write_s_stanza()
450	self.write_h_stanza()
451	for alignment in self.blockHash[key]:
452	self.write_a_stanza(alignment)
453	self.write_trailer()
454	if (self.file != sys.stdout): self.file.close()
455
456	def write_s_stanza(self):
457	self.write_lav_marker()
458	(strand1,flag1) = minus_or_nothing(self.strand1)
459	(strand2,flag2) = minus_or_nothing(self.strand2)
460	fname1 = build_filename(self.fname1,self.src1)
461	fname2 = build_filename(self.fname2,self.src2)
462	print >>self.file,"s {"
463	print >>self.file," \"%s%s\" 1 %d %d 1" \
464	% (fname1,strand1,self.length1,flag1)
465	print >>self.file," \"%s%s\" 1 %d %d 1" \
466	% (fname2,strand2,self.length2,flag2)
467	print >>self.file,"}"
468
469	def write_h_stanza(self):
470	strand1 = rc_or_nothing(self.strand1)
471	strand2 = rc_or_nothing(self.strand2)
472	print >>self.file,"h {"
473	print >>self.file," \"> %s%s\"" % (self.src1,strand1)
474	print >>self.file," \"> %s%s\"" % (self.src2,strand2)
475	print >>self.file,"}"
476
477	def write_a_stanza(self,alignment):
478	c1 = alignment.components[0]
479	pos1 = c1.start
480	text1 = c1.text.upper()
481	c2 = alignment.components[1]
482	pos2 = c2.start
483	text2 = c2.text.upper()
484
485	# collect ungapped pieces
486
487	pieces = []
488	piece1 = None
489
490	for ix in range(len(text1)):
491	ch1 = text1[ix]
492	ch2 = text2[ix]
493
494	nonGap = (ch1 != "-") and (ch2 != "-")
495	if (nonGap):
496	if (piece1 == None): # new piece starts
497	(piece1,piece2,idCount) = (pos1,pos2,0)
498	if (ch1 == ch2): idCount += 1
499	elif (piece1 != None): # new gap starts
500	size = pos1 - piece1
501	pctId = (200idCount + size) / (2size)
502	pieces.append((piece1,piece2,size,pctId))
503	piece1 = None
504
505	if (ch1 != "-"): pos1 += 1
506	if (ch2 != "-"): pos2 += 1
507
508	if (piece1 != None):
509	size = pos1 - piece1
510	pctId = (200idCount + size) / (2size)
511	pieces.append((piece1,piece2,size,pctId))
512
513	# write the block
514
515	(start1,start2,size,pctId) = pieces[-1] # get end of final piece
516	end1 = start1 + size
517	end2 = start2 + size
518
519	(start1,start2,size,pctId) = pieces[0] # get start of first piece
520
521	score = int(round(alignment.score))
522
523	print >>self.file,"a {"
524	print >>self.file," s %s" % score
525	print >>self.file," b %d %d" % (start1+1,start2+1)
526	print >>self.file," e %d %d" % (end1, end2)
527	for (start1,start2,size,pctId) in pieces:
528	print >>self.file," l %d %d %d %d %d" \
529	% (start1+1,start2+1,start1+size,start2+size,pctId)
530	print >>self.file,"}"
531
532	def write_lav_marker(self):
533	print >>self.file,"#:lav"
534
535	def write_trailer(self):
536	print >>self.file,"#:eof"
537
538	def sort_keys_by_chrom (keys):
539	decorated = [(chrom_key(src1),strand1,chrom_key(src2),strand2,(src1,strand1,src2,strand2)) \
540	for (src1,strand1,src2,strand2) in keys]
541	decorated.sort()
542	return [key for (src1,strand1,src2,strand2,key) in decorated]
543
544	def chrom_key (src):
545	(species,chrom) = src_split(src)
546	if (chrom.startswith("chr")): chrom = chrom[3:]
547	try: chrom = int(chrom)
548	except ValueError: pass
549	return chrom
550
551	def build_filename(fmt,src):
552	if (fmt == None): return src
553	num = fmt.count("%s")
554	if (num == 0): return fmt
555	(species,chrom) = src_split(src)
556	if (num == 1): return fmt % chrom
557	return fmt % (species,chrom)
558
559	def minus_or_nothing(strand):
560	if (strand == "-"): return ("-",1)
561	else: return ("",0)
562
563	def rc_or_nothing(strand):
564	if (strand == "-"): return " (reverse complement)"
565	else: return ""
566
567	def do_path_subs(path,path_subs):
568	for (prefix,replacement) in path_subs:
569	if (path.startswith(prefix)):
570	return replacement + path[len(prefix):]
571	return path

Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。

Context Navigation

root/galaxy-central/eggs/bx_python-0.5.0_dev_f74aec067563-py2.6-macosx-10.6-universal-ucs2.egg/bx/align/lav.py

異なるフォーマットでダウンロード: