allpy: e361a7b7d9aa allpy/base.py

allpy

view allpy/base.py @ 262:e361a7b7d9aa

Moved contents of allpy.sequence into allpy.base

author	Daniil Alexeyevsky <me.dendik@gmail.com>
date	Tue, 14 Dec 2010 21:06:42 +0300
parents	d60628e29b24
children	e3783fca343e

line source

1 import sys

2 import os

3 import os.path

4 from tempfile import NamedTemporaryFile

5 import urllib2

7 import config

8 from graph import Graph

9 from Bio.PDB import Superimposer, CaPPBuilder, PDBIO

10 from Bio.PDB.DSSP import make_dssp_dict

11 from allpy.pdb import std_id, pdb_id_parse, get_structure

12 from fasta import save_fasta

13 import data.codes

15 class MonomerType(object):

16 """Class of monomer types.

18 Each MonomerType object represents a known monomer type, e.g. Valine,

19 and is referenced to by each instance of monomer in a given sequence.

21 - `name`: full name of monomer type

22 - `code1`: one-letter code

23 - `code3`: three-letter code

24 - `is_modified`: either of True or False

26 class atributes:

28 - `by_code1`: a mapping from one-letter code to MonomerType object

29 - `by_code3`: a mapping from three-letter code to MonomerType object

30 - `by_name`: a mapping from monomer name to MonomerType object

31 - `instance_type`: class of Monomer objects to use when creating new

32 objects; this must be redefined in descendent classes

34 All of the class attributes MUST be redefined when subclassing.

35 """

37 by_code1 = {}

38 by_code3 = {}

39 by_name = {}

40 instance_type = None

42 def __init__(self, name="", code1="", code3="", is_modified=False):

43 self.name = name.capitalize()

44 self.code1 = code1.upper()

45 self.code3 = code3.upper()

46 self.is_modified = bool(is_modified)

47 if not is_modified:

48 self.by_code1[self.code1] = self

49 self.by_code3[code3] = self

50 self.by_name[name] = self

51 # We duplicate distinguished long names into MonomerType itself,

52 # so that we can use MonomerType.from_code3 to create the relevant

53 # type of monomer.

54 MonomerType.by_code3[code3] = self

55 MonomerType.by_name[name] = self

57 @classmethod

58 def _initialize(cls, type_letter, codes=data.codes.codes):

59 """Create all relevant instances of MonomerType.

61 `type_letter` is either of:

63 - 'p' for protein

64 - 'd' for DNA

65 - 'r' for RNA

67 `codes` is a table of monomer codes

68 """

69 for type, code1, is_modified, code3, name in codes:

70 if type == type_letter:

71 cls(name, code1, code3, is_modified)

73 @classmethod

74 def from_code1(cls, code1):

75 """Return monomer type by one-letter code."""

76 return cls.by_code1[code1.upper()]

78 @classmethod

79 def from_code3(cls, code3):

80 """Return monomer type by three-letter code."""

81 return cls.by_code3[code3.upper()]

83 @classmethod

84 def from_name(cls, name):

85 """Return monomer type by name."""

86 return cls.by_name[name.capitalize()]

88 def instance(self):

89 """Create a new monomer of given type."""

90 return self.instance_type(self)

92 def __eq__(self, other):

93 if hasattr(other, "type"):

94 return self is other.type

95 return self is other

97 class Monomer(object):

98 """Monomer object.

100 attributes:

101

102 - `type`: type of monomer (a MonomerType object)

103

104 class attribute `monomer_type` is MonomerType or either of it's subclasses,

105 it is used when creating new monomers. It MUST be redefined when subclassing Monomer.

106 """

107 monomer_type = MonomerType

108

109 def __init__(self, type):

110 self.type = type

111

112 @classmethod

113 def from_code1(cls, code1):

114 return cls(cls.monomer_type.by_code1[code1.upper()])

115

116 @classmethod

117 def from_code3(cls, code3):

118 return cls(cls.monomer_type.by_code3[code3.upper()])

119

120 @classmethod

121 def from_name(cls, name):

122 return cls(cls.monomer_type.by_name[name.capitalize()])

123

124 def __eq__(self, other):

125 if hasattr(other, "type"):

126 return self.type is other.type

127 return self.type is other

128

129 class Sequence(list):

130 """ Sequence of Monomers

131

132 list of monomer objects (aminoacids or nucleotides)

133

134 Mandatory data:

135 * name -- str with the name of sequence

136 * description -- str with description of the sequence

137

138 Optional (may be empty):

139 * source -- source of sequence

140 * pdb_chain -- Bio.PDB.Chain

141 * pdb_file -- file object

142

143 * pdb_residues -- {Monomer: Bio.PDB.Residue}

144 * pdb_secstr -- {Monomer: 'Secondary structure'}

145 Code Secondary structure

146 H alpha-helix

147 B Isolated beta-bridge residue

148 E Strand

149 G 3-10 helix

150 I pi-helix

151 T Turn

152 S Bend

153 - Other

154

155

156 ?TODO: global pdb_structures

157 """

158 def __init__(self, monomers=None, name='', description=""):

159 if not monomers:

160 monomers = []

161 self.name = name

162 self.description = description

163 self.monomers = monomers

164 self.pdb_chains = []

165 self.pdb_files = {}

166 self.pdb_residues = {}

167 self.pdb_secstr = {}

168

169 def __len__(self):

170 return len(self.monomers)

171

172 def __str__(self):

173 """ Returns sequence in one-letter code """

174 return ''.join([monomer.type.code1 for monomer in self.monomers])

175

176 def __eq__(self, other):

177 """ Returns if all corresponding monomers of this sequences are equal

178

179 If lengths of sequences are not equal, returns False

180 """

181 return len(self) == len(other) and \

182 all([a==b for a, b in zip(self.monomers, other.monomers)])

183

184 def __ne__(self, other):

185 return not (self == other)

186

187 def set_pdb_chain(self, pdb_file, pdb_id, pdb_chain, pdb_model=0):

188 """ Reads Pdb chain from file

189

190 and align each Monomer with PDB.Residue (TODO)

191 """

192 name = std_id(pdb_id, pdb_chain, pdb_model)

193 structure = get_structure(pdb_file, name)

194 chain = structure[pdb_model][pdb_chain]

195 self.pdb_chains.append(chain)

196 self.pdb_residues[chain] = {}

197 self.pdb_secstr[chain] = {}

198 pdb_sequence = Sequence.from_pdb_chain(chain)

199 a = alignment.Alignment.from_sequences(self, pdb_sequence)

200 a.muscle_align()

201 for monomer, pdb_monomer in a.column(sequence=pdb_sequence, original=self):

202 if pdb_sequence.pdb_has(chain, pdb_monomer):

203 residue = pdb_sequence.pdb_residues[chain][pdb_monomer]

204 self.pdb_residues[chain][monomer] = residue

205 self.pdb_files[chain] = pdb_file

206

207 def pdb_unload(self):

208 """ Delete all pdb-connected links """

209 #~ gc.get_referrers(self.pdb_chains[0])

210 self.pdb_chains = []

211 self.pdb_residues = {}

212 self.pdb_secstr = {} # FIXME

213 self.pdb_files = {} # FIXME

214

215 @staticmethod

216 def from_str(fasta_str, name='', description='', monomer_kind=AminoAcidType):

217 """ Import data from one-letter code

218

219 monomer_kind is class, inherited from MonomerType

220 """

221 monomers = [monomer_kind.from_code1(aa).instance() for aa in fasta_str]

222 return Sequence(monomers, name, description)

223

224 @staticmethod

225 def from_pdb_chain(chain):

226 """ Returns Sequence with Monomers with link to Bio.PDB.Residue

227

228 chain is Bio.PDB.Chain

229 """

230 cappbuilder = CaPPBuilder()

231 peptides = cappbuilder.build_peptides(chain)

232 sequence = Sequence()

233 sequence.pdb_chains = [chain]

234 sequence.pdb_residues[chain] = {}

235 sequence.pdb_secstr[chain] = {}

236 for peptide in peptides:

237 for ca_atom in peptide.get_ca_list():

238 residue = ca_atom.get_parent()

239 monomer = AminoAcidType.from_pdb_residue(residue).instance()

240 sequence.pdb_residues[chain][monomer] = residue

241 sequence.monomers.append(monomer)

242 return sequence

243

244 def pdb_auto_add(self, conformity_info=None, pdb_directory='./tmp'):

245 """ Adds pdb information to each monomer

246

247 Returns if information has been successfully added

248 TODO: conformity_file

249

250 id-format lava flow

251 """

252 if not conformity_info:

253 path = os.path.join(pdb_directory, self.name)

254 if os.path.exists(path) and os.path.getsize(path):

255 match = pdb_id_parse(self.name)

256 self.pdb_chain_add(open(path), match['code'],

257 match['chain'], match['model'])

258 else:

259 match = pdb_id_parse(self.name)

260 if match:

261 code = match['code']

262 pdb_filename = config.pdb_dir % code

263 if not os.path.exists(pdb_filename) or not os.path.getsize(pdb_filename):

264 url = config.pdb_url % code

265 print "Download %s" % url

266 pdb_file = open(pdb_filename, 'w')

267 data = urllib2.urlopen(url).read()

268 pdb_file.write(data)

269 pdb_file.close()

270 print "Save %s" % pdb_filename

271 pdb_file = open(pdb_filename)

272 self.pdb_chain_add(pdb_file, code, match['chain'], match['model'])

273

274 def pdb_save(self, out_filename, pdb_chain):

275 """ Saves pdb_chain to out_file """

276 class GlySelect(Select):

277 def accept_chain(self, chain):

278 if chain == pdb_chain:

279 return 1

280 else:

281 return 0

282 io = PDBIO()

283 structure = chain.get_parent()

284 io.set_structure(structure)

285 io.save(out_filename, GlySelect())

286

287

288 def pdb_add_sec_str(self, pdb_chain):

289 """ Add secondary structure data """

290 tmp_file = NamedTemporaryFile(delete=False)

291 tmp_file.close()

292 pdb_file = self.pdb_files[pdb_chain].name

293 os.system("dsspcmbi %(pdb)s %(tmp)s" % {'pdb': pdb_file, 'tmp': tmp_file.name})

294 dssp, keys = make_dssp_dict(tmp_file.name)

295 for monomer in self.monomers:

296 if self.pdb_has(pdb_chain, monomer):

297 residue = self.pdb_residues[pdb_chain][monomer]

298 try:

299 d = dssp[(pdb_chain.get_id(), residue.get_id())]

300 self.pdb_secstr[pdb_chain][monomer] = d[1]

301 except:

302 print "No dssp information about %s at %s" % (monomer, pdb_chain)

303 os.unlink(tmp_file.name)

304

305 def pdb_has(self, chain, monomer):

306 return chain in self.pdb_residues and monomer in self.pdb_residues[chain]

307

308 def secstr_has(self, chain, monomer):

309 return chain in self.pdb_secstr and monomer in self.pdb_secstr[chain]

310

311 @staticmethod

312 def file_slice(file, n_from, n_to, fasta_name='', name='', description='', monomer_kind=AminoAcidType):

313 """ Build and return sequence, consisting of part of sequence from file

314

315 Does not control gaps

316 """

317 inside = False

318 number_used = 0

319 s = ''

320 for line in file:

321 line = line.split()

322 if not inside:

323 if line.startswith('>%s' % fasta_name):

324 inside = True

325 else:

326 n = len(line)

327 s += line[(n_from - number_user):(n_to - number_user)]

328 return Sequence.from_str(s, name, description, monomer_kind)

329

330 class Alignment(dict):

331 """ Alignment

332

333 {<Sequence object>:[<Monomer object>,None,<Monomer object>]}

334 keys are the Sequence objects, values are the lists, which

335 contain monomers of those sequences or None for gaps in the

336 corresponding sequence of alignment

337 """

338 # _sequences -- list of Sequence objects. Sequences don't contain gaps

339 # - see sequence.py module

340

341 def __init__(self, *args):

342 """overloaded constructor

343

344 Alignment() -> new empty Alignment

345 Alignment(sequences, body) -> new Alignment with sequences and

346 body initialized from arguments

347 Alignment(fasta_file) -> new Alignment, read body and sequences

348 from fasta file

349

350 """

351 if len(args)>1:#overloaded constructor

352 self.sequences=args[0]

353 self.body=args[1]

354 elif len(args)==0:

355 self.sequences=[]

356 self.body={}

357 else:

358 self.sequences, self.body = Alignment.from_fasta(args[0])

359

360 def length(self):

361 """ Returns width, ie length of each sequence with gaps """

362 return max([len(line) for line in self.body.values()])

363

364 def height(self):

365 """ The number of sequences in alignment (it's thickness). """

366 return len(self.body)

367

368 def identity(self):

369 """ Calculate the identity of alignment positions for colouring.

370

371 For every (row, column) in alignment the percentage of the exactly

372 same residue in the same column in the alignment is calculated.

373 The data structure is just like the Alignment.body, but istead of

374 monomers it contains float percentages.

375 """

376 # Oh, God, that's awful! Absolutely not understandable.

377 # First, calculate percentages of amino acids in every column

378 contribution = 1.0 / len(self.sequences)

379 all_columns = []

380 for position in range(len(self)):

381 column_percentage = {}

382 for seq in self.body:

383 if self.body[seq][position] is not None:

384 aa = self.body[seq][position].code

385 else:

386 aa = None

387 if aa in allpy.data.amino_acids:

388 if aa in column_percentage.keys():

389 column_percentage[aa] += contribution

390 else:

391 column_percentage[aa] = contribution

392 all_columns.append(column_percentage)

393 # Second, map these percentages onto the alignment

394 self.identity_percentages = {}

395 for seq in self.sequences:

396 self.identity_percentages[seq] = []

397 for seq in self.identity_percentages:

398 line = self.identity_percentages[seq]

399 for position in range(len(self)):

400 if self.body[seq][position] is not None:

401 aa = self.body[seq][position].code

402 else:

403 aa = None

404 line.append(all_columns[position].get(aa))

405 return self.identity_percentages

406

407 @staticmethod

408 def from_fasta(file, monomer_kind=AminoAcidType):

409 """ Import data from fasta file

410

411 monomer_kind is class, inherited from MonomerType

412

413 >>> import alignment

414 >>> sequences,body=alignment.Alignment.from_fasta(open("test.fasta"))

415 """

416 import re

417

418 sequences = []

419 body = {}

420

421 raw_sequences = file.read().split(">")

422 if len(raw_sequences) <= 1:

423 raise Exception("Wrong format of fasta-file %s" % file.name)

424

425 raw_sequences = raw_sequences[1:] #ignore everything before the first >

426 for raw in raw_sequences:

427 parsed_raw_sequence = raw.split("\n")

428 parsed_raw_sequence = [s.strip() for s in parsed_raw_sequence]

429 name_and_description = parsed_raw_sequence[0]

430 name_and_description = name_and_description.split(" ",1)

431 if len(name_and_description) == 2:

432 name, description = name_and_description

433 elif len(name_and_description) == 1:

434 #if there is description

435 name = name_and_description[0]

436 description = ''

437 else:

438 raise Exception("Wrong name of sequence %(name)$ fasta-file %(file)s" % \

439 {'name': name, 'file': file.name})

440

441 if len(parsed_raw_sequence) <= 1:

442 raise Exception("Wrong format of sequence %(name)$ fasta-file %(file)s" % \

443 {'name': name, 'file': file.name})

444 string = ""

445 for piece in parsed_raw_sequence[1:]:

446 piece_without_whitespace_chars = re.sub("\s", "", piece)

447 string += piece_without_whitespace_chars

448 monomers = [] #convert into Monomer objects

449 body_list = [] #create the respective list in body dict

450 for current_monomer in string:

451 if current_monomer not in ["-", ".", "~"]:

452 monomers.append(monomer_kind.from_code1(current_monomer).instance())

453 body_list.append(monomers[-1])

454 else:

455 body_list.append(None)

456 s = sequence.Sequence(monomers, name, description)

457 sequences.append(s)

458 body[s] = body_list

459 return sequences, body

460

461 @staticmethod

462 def from_sequences(*sequences):

463 """ Constructs new alignment from sequences

464

465 Add None's to right end to make equal lengthes of alignment sequences

466 """

467 alignment = Alignment()

468 alignment.sequences = sequences

469 max_length = max(len(sequence) for sequence in sequences)

470 for sequence in sequences:

471 gaps_count = max_length - len(sequence)

472 alignment.body[sequence] = sequence.monomers + [None] * gaps_count

473 return alignment

474

475 def save_fasta(self, out_file, long_line=70, gap='-'):

476 """ Saves alignment to given file

477

478 Splits long lines to substrings of length=long_line

479 To prevent this, set long_line=None

480 """

481 block.Block(self).save_fasta(out_file, long_line=long_line, gap=gap)

482

483 def muscle_align(self):

484 """ Simple align ths alignment using sequences (muscle)

485

486 uses old Monomers and Sequences objects

487 """

488 tmp_file = NamedTemporaryFile(delete=False)

489 self.save_fasta(tmp_file)

490 tmp_file.close()

491 os.system("muscle -in %(tmp)s -out %(tmp)s" % {'tmp': tmp_file.name})

492 sequences, body = Alignment.from_fasta(open(tmp_file.name))

493 for sequence in self.sequences:

494 try:

495 new_sequence = [i for i in sequences if sequence==i][0]

496 except:

497 raise Exception("Align: Cann't find sequence %s in muscle output" % \

498 sequence.name)

499 old_monomers = iter(sequence.monomers)

500 self.body[sequence] = []

501 for monomer in body[new_sequence]:

502 if not monomer:

503 self.body[sequence].append(monomer)

504 else:

505 old_monomer = old_monomers.next()

506 if monomer != old_monomer:

507 raise Exception("Align: alignment errors")

508 self.body[sequence].append(old_monomer)

509 os.unlink(tmp_file.name)

510

511 def column(self, sequence=None, sequences=None, original=None):

512 """ returns list of columns of alignment

513

514 sequence or sequences:

515 if sequence is given, then column is (original_monomer, monomer)

516 if sequences is given, then column is (original_monomer, {sequence: monomer})

517 if both of them are given, it is an error

518 original (Sequence type):

519 if given, this filters only columns represented by original sequence

520 """

521 if sequence and sequences:

522 raise Exception("Wrong usage. read help")

523 indexes = dict([(v, k) for( k, v) in enumerate(self.sequences)])

524 alignment = self.body.items()

525 alignment.sort(key=lambda i: indexes[i[0]])

526 alignment = [monomers for seq, monomers in alignment]

527 for column in zip(*alignment):

528 if not original or column[indexes[original]]:

529 if sequence:

530 yield (column[indexes[original]], column[indexes[sequence]])

531 else:

532 yield (column[indexes[original]],

533 dict([(s, column[indexes[s]]) for s in sequences]))

534

535 def secstr(self, sequence, pdb_chain, gap='-'):

536 """ Returns string representing secondary structure """

537 return ''.join([

538 (sequence.pdb_secstr[pdb_chain][m] if sequence.secstr_has(pdb_chain, m) else gap)

539 for m in self.body[sequence]])

540

541 class Block(object):

542 """ Block of alignment

543

544 Mandatory data:

545 * self.alignment -- alignment object, which the block belongs to

546 * self.sequences - set of sequence objects that contain monomers

547 and/or gaps, that constitute the block

548 * self.positions -- list of positions of the alignment.body that

549 are included in the block; position[i+1] is always to the right from position[i]

550

551 Don't change self.sequences -- it may be a link to other block.sequences

552

553 How to create a new block:

554 >>> import alignment

555 >>> import block

556 >>> proj = alignment.Alignment(open("test.fasta"))

557 >>> block1 = block.Block(proj)

558 """

559

560 def __init__(self, alignment, sequences=None, positions=None):

561 """ Builds new block from alignment

562

563 if sequences==None, all sequences are used

564 if positions==None, all positions are used

565 """

566 if sequences == None:

567 sequences = set(alignment.sequences) # copy

568 if positions == None:

569 positions = range(len(alignment))

570 self.alignment = alignment

571 self.sequences = sequences

572 self.positions = positions

573

574 def save_fasta(self, out_file, long_line=70, gap='-'):

575 """ Saves alignment to given file in fasta-format

576

577 No changes in the names, descriptions or order of the sequences

578 are made.

579 """

580 for sequence in self.sequences:

581 alignment_monomers = self.alignment.body[sequence]

582 block_monomers = [alignment_monomers[i] for i in self.positions]

583 string = ''.join([m.type.code1 if m else '-' for m in block_monomers])

584 save_fasta(out_file, string, sequence.name, sequence.description, long_line)

585

586 def geometrical_cores(self, max_delta=config.delta,

587 timeout=config.timeout, minsize=config.minsize,

588 ac_new_atoms=config.ac_new_atoms,

589 ac_count=config.ac_count):

590 """ Returns length-sorted list of blocks, representing GCs

591

592 max_delta -- threshold of distance spreading

593 timeout -- Bron-Kerbosh timeout (then fast O(n ln n) algorithm)

594 minsize -- min size of each core

595 ac_new_atoms -- min part or new atoms in new alternative core

596 current GC is compared with each of already selected GCs

597 if difference is less then ac_new_atoms, current GC is skipped

598 difference = part of new atoms in current core

599 ac_count -- max number of cores (including main core)

600 -1 means infinity

601 If more than one pdb chain for some sequence provided, consider all of them

602 cost is calculated as 1 / (delta + 1)

603 delta in [0, +inf) => cost in (0, 1]

604 """

605 nodes = self.positions

606 lines = {}

607 for i in self.positions:

608 for j in self.positions:

609 if i < j:

610 distances = []

611 for sequence in self.sequences:

612 for chain in sequence.pdb_chains:

613 m1 = self.alignment.body[sequence][i]

614 m2 = self.alignment.body[sequence][j]

615 if m1 and m2:

616 r1 = sequence.pdb_residues[chain][m1]

617 r2 = sequence.pdb_residues[chain][m2]

618 ca1 = r1['CA']

619 ca2 = r2['CA']

620 d = ca1 - ca2 # Bio.PDB feature

621 distances.append(d)

622 if len(distances) >= 2:

623 delta = max(distances) - min(distances)

624 if delta <= max_delta:

625 lines[Graph.line(i, j)] = 1.0 / (1.0 + max_delta)

626 graph = Graph(nodes, lines)

627 cliques = graph.cliques(timeout=timeout, minsize=minsize)

628 GCs = []

629 for clique in cliques:

630 for GC in GCs:

631 if len(clique - set(GC.positions)) < ac_new_atoms * len(clique):

632 break

633 else:

634 GCs.append(Block(self.alignment, self.sequences, clique))

635 if ac_count != -1 and len(GCs) >= ac_count:

636 break

637 return GCs

638

639 def xstring(self, x='X', gap='-'):

640 """ Returns string consisting of gap chars and chars x at self.positions

641

642 Length of returning string = length of alignment

643 """

644 monomers = [False] * len(self.alignment)

645 for i in self.positions:

646 monomers[i] = True

647 return ''.join([x if m else gap for m in monomers])

648

649 def save_xstring(self, out_file, name, description='', x='X', gap='-', long_line=70):

650 """ Save xstring and name in fasta format """

651 save_fasta(out_file, self.xstring(x=x, gap=gap), name, description, long_line)

652

653 def monomers(self, sequence):

654 """ Iterates monomers of this sequence from this block """

655 alignment_sequence = self.alignment.body[sequence]

656 return (alignment_sequence[i] for i in self.positions)

657

658 def ca_atoms(self, sequence, pdb_chain):

659 """ Iterates Ca-atom of monomers of this sequence from this block """

660 return (sequence.pdb_residues[pdb_chain][monomer] for monomer in self.monomers())

661

662 def sequences_chains(self):

663 """ Iterates pairs (sequence, chain) """

664 for sequence in self.alignment.sequences:

665 if sequence in self.sequences:

666 for chain in sequence.pdb_chains:

667 yield (sequence, chain)

668

669 def superimpose(self):

670 """ Superimpose all pdb_chains in this block """

671 sequences_chains = list(self.sequences_chains())

672 if len(sequences_chains) >= 1:

673 sup = Superimposer()

674 fixed_sequence, fixed_chain = sequences_chains.pop()

675 fixed_atoms = self.ca_atoms(fixed_sequence, fixed_chain)

676 for sequence, chain in sequences_chains:

677 moving_atoms = self.ca_atoms(sequence, chain)

678 sup.set_atoms(fixed_atoms, moving_atoms)

679 # Apply rotation/translation to the moving atoms

680 sup.apply(moving_atoms)

681

682 def pdb_save(self, out_file):

683 """ Save all sequences

684

685 Returns {(sequence, chain): CHAIN}

686 CHAIN is chain letter in new file

687 """

688 tmp_file = NamedTemporaryFile(delete=False)

689 tmp_file.close()

690

691 for sequence, chain in self.sequences_chains():

692 sequence.pdb_save(tmp_file.name, chain)

693 # TODO: read from tmp_file.name

694 # change CHAIN

695 # add to out_file

696

697 os.unlink(NamedTemporaryFile)

698

699 # vim: set ts=4 sts=4 sw=4 et: