allpy: b556c96c6719 utils/rasmol

allpy

view utils/rasmol_homology.py @ 1168:b556c96c6719

blocks3d/www Makefile: never check certificates of github, they are too confusing for wget

author	Daniil Alexeyevsky <dendik@kodomo.fbb.msu.ru>
date	Mon, 26 May 2014 17:20:29 +0400
parents	41f7241c2aeb
children

line source

1 #!/usr/bin/python

2 import sys

3 import optparse

4 from allpy import protein, structure, markups, homology

6 flank_size = 5 #flanks' size

8 class Alignment(protein.Alignment, structure.AlignmentMixin):

9 pass

10 class Block(protein.Block, structure.BlockMixin):

11 pass

12 class Sequence(protein.Sequence, structure.SequenceMixin):

13 pass

15 def write(line):

16 options.output_rasmol.write(line + '\n')

17 def warn(text):

18 sys.stderr.write(text + '\n')

20 parser = optparse.OptionParser()

21 parser.add_option('-m', '--markup',

22 help='Markup file to use')

23 parser.add_option('-H', '--homology',

24 help='Homology file to use')

25 parser.add_option('-w', '--min-width', type=int,

26 help='Minimal block width [optional]')

27 parser.add_option('-d', '--min-depth', type=int,

28 help='Minimal block depth [optional]')

29 parser.add_option('-o', '--output-pdb',

30 help='Output PDB file name')

31 parser.add_option('-O', '--output-rasmol',

32 help='Output RasMol script file name [default: stdout]')

33 parser.add_option('-c', '--pdb-cache',

34 help='PDB cache directory [optional]')

35 options, args = parser.parse_args()

37 assert options.markup and options.homology and options.output_pdb

39 if options.output_rasmol:

40 options.output_rasmol = open(options.output_rasmol, 'w')

41 else:

42 options.output_rasmol = sys.stdout

44 download_pdb = structure.CachedDownloadPdb(cache_dir=options.pdb_cache)

46 alignment = Alignment().append_file(open(options.markup), format='markup')

47 alignment.add_markup('number')

49 def pdb_loader(sequence):

50 sequence.__class__ = Sequence

51 sequence.markups['pdb_resi'].add_pdb(download_pdb=download_pdb)

52 sequence.add_markup('number')

54 hmg = homology.MonomerHomology()

55 hmg.read(options.homology, columns=True)

56 hmg.highest_blocks()

57 blocks = hmg.alignment_blocks(alignment)

59 columns = set()

60 for block in blocks:

61 columns |= set(block.columns)

62 columns = list(columns)

64 chains = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"

65 pdb_file = open(options.output_pdb, 'w')

66 monomer2id = {}

67 idmap = {}

68 def process_sequence(sequence, i):

69 for monomer in sequence:

70 monomer2id[monomer] = monomer.pdb_residue.get_id()[1]

71 chain = chains[i % len(chains)]

72 model = i // len(chains)

73 idmap[sequence] = (chain, model)

74 sequence.save_pdb(pdb_file, chain, model)

75 main_sequence = alignment.sequences[-1]

76 pdb_loader(main_sequence)

77 for i, sequence in enumerate(alignment.sequences[:-1]):

78 pdb_loader(sequence)

79 sequence.superimpose(main_sequence, columns, extra_columns=True)

80 process_sequence(sequence, i)

81 sequence.pdb_unload()

82 process_sequence(main_sequence, len(alignment.sequences)-1)

83 main_sequence.pdb_unload()

85 def sequence2name(sequence):

86 return 'str_%s' % sequence.name.replace('/','_').replace(':','_')

88 sequence_names = []

89 for sequence, chain_model in idmap.items():

90 chain, model = chain_model

91 write('define %s *:%s' % (sequence2name(sequence), chain))

92 sequence_names.append(sequence2name(sequence))

94 write('select all')

95 write('color chain')

96 write('restrict none')

97 write('background [230, 240, 250]')

98 write('select none')

99 write('select %s' % ', '.join(sequence_names))

100 write('backbone')

101 write('pause')

102

103 for num, block in enumerate(blocks, 1):

104 name = 'block_%s' % num

105

106 if options.min_width and len(block.columns) < options.min_width:

107 warn("Skipping %s of width %s" % (name, len(block.columns)))

108 continue

109 if options.min_depth and len(block.sequences) < options.min_depth:

110 warn("Skipping %s of depth %s" % (name, len(block.sequences)))

111 continue

112

113 with_flanks = []

114 parts = []

115 for sequence in block.sequences:

116 chain, model = idmap[sequence]

117 resi_from = monomer2id[block.columns[0][sequence]]

118 resi_to = monomer2id[block.columns[-1][sequence]]

119 with_flanks.append('%i-%i:%s' % (resi_from-flank_size, resi_to+flank_size, chain))

120 parts.append('%i-%i:%s' % (resi_from, resi_to, chain))

121

122 selection = ','.join(parts)

123 sequences = ', '.join(sequence2name(seq) for seq in block.sequences)

124 column_from = alignment.markups['number'][block.columns[0]]

125 column_to = alignment.markups['number'][block.columns[-1]]

126

127 write('')

128 write('define with_flanks_%i %s' % (num, ','.join(with_flanks)))

129 write('define %s %s' % (name, selection))

130 write('select with_flanks_%i' % num)

131 write('backbone 60')

132 write('select %s' % (name))

133 write('echo %s %s %i-%i' % (name, sequences, column_from, column_to))

134 write('backbone 150')

135 write('pause')

136 write('select with_flanks_%i' % num)

137 write('backbone ')

138

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