allpy: 670ae384aa37 test/usecase2.py

allpy

view test/usecase2.py @ 372:670ae384aa37

Fixed bugs & refactored base.Alignment.process & it's helpers - there were a few name errors - process would previously behave bad if new alignment was shorter - refactored and simplified code for copying sequence contents

author	Daniil Alexeyevsky <me.dendik@gmail.com>
date	Mon, 31 Jan 2011 18:44:23 +0300
parents	cc9f4972ea46
children	bed32775625a

line source

1 # Fragments are in pair_repeat.fasta

2 from allpy import dna

3 width = 15

4 treshold = 14

6 def my_column_mark(column, threshold):

7 """Helper to mark column (given as dict) by identity."""

8 count = {}

9 for sequence, monomer in column:

10 code = monomer.code1

11 count[code] = count.get(code, 0) + 1

12 for code in count:

13 if count[code] > threshold:

14 return "+"

15 return "-"

17 def my_pair_mark(column):

18 """Helper to mark column of 2 sequences (given as list) by identity."""

19 if column[0] is None or column[1] is None:

20 return "-"

21 if column[0].code1 == column[0].code1:

22 return "+"

23 return "-"

25 def find_runs(markup):

26 """Fund long positive runs.

28 This obscure and probably broken function has nothing to do with allpy,

29 so it's presence in the example is unnecessary.

30 """

31 position = 0

32 count = 0

33 plus_positions=[]

34 for i in range(len(markup)):

35 position += 1

36 if position < width :

37 if markup[i]=="+":

38 count += 1

39 continue

40 if position > width:

41 if markup[i-width] == "+":

42 count -=1

43 if count >= treshold:

44 plus_positions.append(position)

45 if len(plus_positions)==0:

46 raise Exception("No blocks in alignment")

48 blocks=[]

49 start = plus_positions[0]-width + 1

50 stop = plus_positions[0]

51 for p in plus_positions[1:]:

52 if p == stop +1:

53 stop = p

54 continue

55 blocks.append((start,stop))

56 start = p - width + 1

57 stop = p

58 return blocks

60 def main():

61 alignment = dna.Alignment.from_file(open("pair_repeat.fasta"))

62 if len(alignment.sequences) != 2:

63 raise Exception("Input must have exactly 2 sequences!")

64 alignment.realign("needle", gap_open = 0)

65 markup = []

66 for column in alignment.columns_as_lists():

67 markup.append(my_pair_mark(column))

68 markup = alignment.map_columns(my_pair_mark)

69 print find_runs(markup)

71 try:

72 main()

73 except Exception, e:

74 print "An error has occured:", e