Mercurial > repos > guerler > springsuite

comparison spring_roc.py @ 34:b300ddbbf9d0 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit 0410e2fadc4e9fc1df6010de7b3678154cbdfe62-dirty"
author guerler
date Tue, 24 Nov 2020 17:55:07 +0000
parents 3071750405c9
children 0bcc0a269916
comparison
equal deleted inserted replaced
33:f115fbf3ac63 34:b300ddbbf9d0
1 #! /usr/bin/env python 1 #! /usr/bin/env python
2 import argparse 2 import argparse
3 import math 3 import math
4 import random 4 import random
5 from os.path import isfile 5 from os.path import isfile
6 from datetime import datetime
7 6
8 from matplotlib import pyplot as plt 7 from matplotlib import pyplot as plt
9 8
10 9
11 def getIds(rawIds): 10 def getIds(rawIds):
205 204
206 # process prediction file 205 # process prediction file
207 print("Loading prediction file...") 206 print("Loading prediction file...")
208 prediction, _ = getReference(args.input, scoreCol=2) 207 prediction, _ = getReference(args.input, scoreCol=2)
209 208
210 # get subcellular locations from UniProt export 209 # determine negative set
211 locations = dict() 210 print("Identifying non-interacting pairs...")
212 if isfile(args.locations):
213 regions = list()
214 if args.regions:
215 regions = args.regions.split(",")
216 with open(args.locations) as locFile:
217 for line in locFile:
218 searchKey = "SUBCELLULAR LOCATION"
219 searchPos = line.find(searchKey)
220 if searchPos != -1:
221 uniId = line.split()[0]
222 locStart = searchPos + len(searchKey) + 1
223 locId = line[locStart:].split()[0]
224 if regions:
225 if locId not in regions:
226 continue
227 if uniId in filterA or uniId in filterB:
228 locations[uniId] = locId
229 print("Found %d subcellular locations." % (len(list(locations.keys()))))
230
231 # estimate background noise
232 print("Estimating background noise...")
233 negative = set() 211 negative = set()
234 filterAList = sorted(list(filterA)) 212 if isfile(args.negative):
235 filterBList = sorted(list(filterB)) 213 # load from explicit file
236 negativeRequired = positiveCount 214 with open(args.negative) as file:
237 random.seed(0) 215 for line in file:
238 totalAttempts = int(len(filterAList) * len(filterBList) / 2) 216 cols = line.split()
239 while totalAttempts > 0: 217 nameA = cols[0]
240 totalAttempts = totalAttempts - 1 218 nameB = cols[1]
241 nameA = random.choice(filterAList) 219 key = getKey(nameA, nameB)
242 nameB = random.choice(filterBList) 220 if key not in putative and key not in negative:
243 if locations: 221 negative.add(key)
244 if nameA not in locations or nameB not in locations: 222 else:
245 continue 223 # get subcellular locations from UniProt export
246 if locations[nameA] == locations[nameB]: 224 locations = dict()
247 continue 225 if isfile(args.locations):
248 key = getKey(nameA, nameB) 226 regions = list()
249 if key not in putative and key not in negative: 227 if args.regions:
250 negative.add(key) 228 regions = args.regions.split(",")
251 negativeRequired = negativeRequired - 1 229 with open(args.locations) as locFile:
252 if negativeRequired == 0: 230 for line in locFile:
253 break 231 searchKey = "SUBCELLULAR LOCATION"
232 searchPos = line.find(searchKey)
233 if searchPos != -1:
234 uniId = line.split()[0]
235 locStart = searchPos + len(searchKey) + 1
236 locId = line[locStart:].split()[0]
237 if regions:
238 if locId not in regions:
239 continue
240 if uniId in filterA or uniId in filterB:
241 locations[uniId] = locId
242 print("Found %d subcellular locations." % (len(list(locations.keys()))))
243
244 # randomly sample non-interacting pairs
245 filterAList = sorted(list(filterA))
246 filterBList = sorted(list(filterB))
247 negativeRequired = positiveCount
248 random.seed(0)
249 totalAttempts = int(len(filterAList) * len(filterBList) / 2)
250 while totalAttempts > 0:
251 totalAttempts = totalAttempts - 1
252 nameA = random.choice(filterAList)
253 nameB = random.choice(filterBList)
254 if locations:
255 if nameA not in locations or nameB not in locations:
256 continue
257 if locations[nameA] == locations[nameB]:
258 continue
259 key = getKey(nameA, nameB)
260 if key not in putative and key not in negative:
261 negative.add(key)
262 negativeRequired = negativeRequired - 1
263 if negativeRequired == 0:
264 break
254 265
255 # create plot 266 # create plot
256 print("Producing plot data...") 267 print("Producing plot data...")
257 print("Total count in prediction file: %d." % len(prediction)) 268 print("Total count in prediction file: %d." % len(prediction))
258 print("Total count in positive file: %d." % len(positive)) 269 print("Total count in positive file: %d." % len(positive))
269 plt.savefig(args.output, format="png") 280 plt.savefig(args.output, format="png")
270 281
271 282
272 if __name__ == "__main__": 283 if __name__ == "__main__":
273 parser = argparse.ArgumentParser(description='Create ROC plot.') 284 parser = argparse.ArgumentParser(description='Create ROC plot.')
274 parser.add_argument('-i', '--input', help='Input prediction file.', required=True) 285 parser.add_argument('-i', '--input', help='Input prediction file (2-columns).', required=True)
275 parser.add_argument('-b', '--biogrid', help='BioGRID interaction database file', required=True) 286 parser.add_argument('-b', '--biogrid', help='BioGRID interaction database file', required=True)
276 parser.add_argument('-l', '--locations', help='UniProt export table with subcellular locations', required=False) 287 parser.add_argument('-l', '--locations', help='UniProt export table with subcellular locations', default="", required=False)
277 parser.add_argument('-r', '--regions', help='Comma-separated regions', required=False) 288 parser.add_argument('-r', '--regions', help='Comma-separated regions', required=False)
289 parser.add_argument('-n', '--negative', help='Negative set (2-columns)', default="", required=False)
278 parser.add_argument('-e', '--experiment', help='Type (physical/genetic)', default="", required=False) 290 parser.add_argument('-e', '--experiment', help='Type (physical/genetic)', default="", required=False)
279 parser.add_argument('-t', '--throughput', help='Throughput (low/high)', default="", required=False) 291 parser.add_argument('-t', '--throughput', help='Throughput (low/high)', default="", required=False)
280 parser.add_argument('-m', '--method', help='Method e.g. Two-hybrid', default="", required=False) 292 parser.add_argument('-m', '--method', help='Method e.g. Two-hybrid', default="", required=False)
281 parser.add_argument('-o', '--output', help='Output (png)', required=True) 293 parser.add_argument('-o', '--output', help='Output (png)', required=True)
282 args = parser.parse_args() 294 args = parser.parse_args()