allpy: e83572fff43f test/usecase2.py

allpy

view test/usecase2.py @ 746:e83572fff43f

Roll-back a bug introduces by dirty hand-merge in [723]. (closes #74) (see #76) Boris! Please do not do dirty hand merges! If you did hg fetch here, this bug would not appear! Please, be extremely careful when you do hand merges and double-check your changes. Do a diff with each parent and see what you remove related to the parent! If someone else's code is involved in the merge (which is almost always the case), do that diff twice just to make sure you have not missed anything!

author	Daniil Alexeyevsky <dendik@kodomo.fbb.msu.ru>
date	Mon, 11 Jul 2011 14:29:54 +0400
parents	bed32775625a
children	ddf85d0a8924

line source

1 import sys

2 from allpy import dna

3 from allpy.processors import Needle, Left

4 from allpy.fileio import FastaFile

5 from collections import deque

7 width = 15

8 threshold = 10

10 def has_identity(column):

11 as_list = column.values()

12 return len(column) == 2 and as_list[0] == as_list[1]

14 def is_good_window(window):

15 sum_id = sum(int(has_identity(column)) for column in window)

16 return len(window) == width and sum_id >= threshold

18 def find_runs(alignment):

19 window = deque([], width)

20 blocks = []

21 in_block = False

22 for column in alignment.columns:

23 window.append(column)

24 in_block, was_in_block = is_good_window(window), in_block

25 if in_block and not was_in_block:

26 block = dna.Block.from_alignment(alignment, columns=list(window))

27 blocks.append(block)

28 elif in_block:

29 block.columns.append(column)

30 return blocks

32 def blocks_markup(alignment, blocks):

33 for column in alignment.columns:

34 column.in_block = "-"

35 for block in blocks:

36 for column in block.columns:

37 column.in_block = "+"

38 return "".join(column.in_block for column in alignment.columns)

40 def main():

41 alignment = dna.Alignment().append_file(sys.stdin)

42 assert len(alignment.sequences) == 2, "Input must have TWO sequences!"

43 alignment.realign(Left())

44 alignment.realign(Needle())

45 blocks = find_runs(alignment)

47 for n, block in enumerate(blocks, 1):

48 block.to_file(open("block_%02d.fasta" % n, "w"))

50 alignment.to_file(sys.stdout)

51 FastaFile(sys.stdout).write_string(

52 blocks_markup(alignment, blocks),

53 "markup",

54 "In run with window %s and threshold %s" % (width, threshold)

55 )

57 try:

58 main()

59 except Exception, e:

60 print "An error has occured:", e