Mercurial > repos > iuc > reprof

comparison reprof.py @ 1:141da185be70 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/reprof commit 3aca39eed712c5b3ab8faac62accab588704ffd9
author iuc
date Wed, 02 Dec 2015 15:01:36 -0500
parents 7c33ed152672
children fb0936cf5bef
comparison
equal deleted inserted replaced
0:7c33ed152672 1:141da185be70
11 from Bio.SeqFeature import SeqFeature, FeatureLocation 11 from Bio.SeqFeature import SeqFeature, FeatureLocation
12 12
13 def run_reprof(query_path, modeldir): 13 def run_reprof(query_path, modeldir):
14 outtmp = tempfile.NamedTemporaryFile(delete=False) 14 outtmp = tempfile.NamedTemporaryFile(delete=False)
15 cmd = [ 15 cmd = [
16 './reprof/scripts/reprof', 16 'reprof',
17 '-i', query_path, 17 '-i', query_path,
18 '--modeldir=%s' % modeldir, 18 '--modeldir=%s' % modeldir,
19 '-o', outtmp.name 19 '-o', outtmp.name
20 ] 20 ]
21 subprocess.check_call(cmd) 21 subprocess.check_call(cmd)
73 rec.features.append(region_feat) 73 rec.features.append(region_feat)
74 74
75 with open(path, 'a') as handle: 75 with open(path, 'a') as handle:
76 GFF.write([rec], handle) 76 GFF.write([rec], handle)
77 77
78 def storeClassData(path, id, phel):
79 h = float(phel.count('H')) / float(len(phel))
80 e = float(phel.count('E')) / float(len(phel))
81
82 if h > .45 and e < .05:
83 classification = 'all-alpha'
84 elif h < .05 and e > .45:
85 classification = 'all-beta'
86 elif h > .3 and e > .2:
87 classification = 'alpha-beta'
88 else:
89 classification = 'mixed'
90
91 with open(path, 'a') as handle:
92 handle.write("{0}\t{1}\n".format(id, classification))
93
78 def main(fasta, modeldir): 94 def main(fasta, modeldir):
79 for record in SeqIO.parse(fasta, 'fasta'): 95 for record in SeqIO.parse(fasta, 'fasta'):
80 tmp = tempfile.NamedTemporaryFile(delete=False) 96 tmp = tempfile.NamedTemporaryFile(delete=False)
81 SeqIO.write([record], tmp, 'fasta') 97 SeqIO.write([record], tmp, 'fasta')
82 tmp.close() 98 tmp.close()
84 # Run reprof 100 # Run reprof
85 data = process_reprof_report(run_reprof(tmp.name, modeldir)) 101 data = process_reprof_report(run_reprof(tmp.name, modeldir))
86 for col in ('RI_S', 'P10', 'RI_A', 'PACC', 'PREL', 'pH', 'pE', 'pL'): 102 for col in ('RI_S', 'P10', 'RI_A', 'PACC', 'PREL', 'pH', 'pE', 'pL'):
87 storeWigData(data['idx'], data[col], record.id, col + '.wig') 103 storeWigData(data['idx'], data[col], record.id, col + '.wig')
88 104
105 eco = ['{ECO:0000255|reprof_1.0.1}']
89 storeGff3Data( 106 storeGff3Data(
90 'secondary_structure.gff3', record.id, data['idx'], data['PHEL'], 107 'secondary_structure.gff3', record.id, data['idx'], data['PHEL'],
91 { 108 {
92 'H': { 109 'H': {
93 'type': 'peptide_helix', 110 'type': 'peptide_helix',
94 'label': ['Helix'], 111 'description': ['Helix'],
95 'evidence': ['ECO:0000255'] 112 'evidence': eco,
96 }, 113 },
97 'E': { 114 'E': {
98 'type': 'beta_strand', 115 'type': 'beta_strand',
99 'label': ['Extended/Sheet'], 116 'description': ['Extended/Sheet'],
100 'evidence': ['ECO:0000255'] 117 'evidence': eco,
101 }, 118 },
102 'L': { 119 'L': {
103 'type': 'loop', 120 'type': 'loop',
104 'label': ['Loop'], 121 'description': ['Loop'],
105 'evidence': ['ECO:0000255'] 122 'evidence': eco,
106 } 123 }
107 } 124 }
108 ) 125 )
109 126
110 storeGff3Data( 127 storeGff3Data(
111 'solvent_accessibility.gff3', record.id, data['idx'], data['Pbe'], 128 'solvent_accessibility.gff3', record.id, data['idx'], data['Pbe'],
112 { 129 {
113 'B': { 130 'B': {
114 'type': 'experimental_result_region', 131 'type': 'experimental_result_region',
115 'label': ['Buried'], 132 'description': ['Buried'],
116 'evidence': ['ECO:0000255'] 133 'evidence': eco,
117 }, 134 },
118 'E': { 135 'E': {
119 'type': 'experimental_result_region', 136 'type': 'experimental_result_region',
120 'label': ['Exposed'], 137 'description': ['Exposed'],
121 'evidence': ['ECO:0000255'] 138 'evidence': eco,
122 }, 139 },
123 } 140 }
124 ) 141 )
142
143 storeClassData('report.tsv', record.id, data['PHEL'])
125 144
126 145
127 if __name__ == '__main__': 146 if __name__ == '__main__':
128 # Grab all of the filters from our plugin loader 147 # Grab all of the filters from our plugin loader
129 parser = argparse.ArgumentParser(description='Wrapper for reprof') 148 parser = argparse.ArgumentParser(description='Wrapper for reprof')