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

リビジョン 3, 17.8 KB (コミッタ: kohda, 14 年 前)

Install Unix tools  http://hannonlab.cshl.edu/galaxy_unix_tools/galaxy.html
行番号 
1"""
2Support for reading and writing the LAV format produced by the `blastz`_
3pairwise aligner.
4
5.. _blastz: http://www.bx.psu.edu/miller_lab/
6"""
7
8from bx.align import *
9import bx.seq
10
11import sys,math,StringIO
12import itertools
13from bx import interval_index_file
14
15class 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
381class 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
392class 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
400class 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 = (200*idCount + size) / (2*size)
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 = (200*idCount + size) / (2*size)
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
538def     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
544def     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
551def 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
559def minus_or_nothing(strand):
560        if (strand == "-"): return ("-",1)
561        else:               return ("",0)
562
563def rc_or_nothing(strand):
564        if (strand == "-"): return " (reverse complement)"
565        else:               return ""
566
567def 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 を参照してください。