allpy: 01547d8d5c36 test/usecase2.py

allpy

view test/usecase2.py @ 326:01547d8d5c36

save_fasta: always prepend space before description, even if no name is given

author	Daniil Alexeyevsky <me.dendik@gmail.com>
date	Fri, 17 Dec 2010 00:57:35 +0300
parents
children	cc9f4972ea46

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_as_list):

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(column_marks)):

35 position += 1

36 if position < width :

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

38 count += 1

39 column_window_marks += "-"

40 continue

41 if position > width:

42 if column_marks[i-width] == "+":

43 count -=1

44 if count >= treshold:

45 plus_positions.append(position)

46 if len(positions)==0:

47 raise Exception("No blocks in alignment")

48 blocks=[]

49 start = positions[0]-width + 1

50 stop = positions[0]

52 for p in positions[1:]:

53 if p == stop +1:

54 stop = p

55 cuntinue

56 blocks.append((start,stop))

57 start = p - width + 1

58 stop = p

60 def main():

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

62 if alignment.height != 2:

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

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

65 markup = alignment.map_columns(my_pair_mark)

66 print find_runs(markup)

68 try:

69 main()

70 except Exception, e:

71 print "An error has occured:", e