[3] | 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 = (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 | |
---|
| 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 |
---|