Mercurial > repos > nick > minor_variant_boxplot

comparison hetbox.py @ 4:dfa2e75da6aa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fixed script typo
author nick
date Thu, 30 May 2013 15:41:27 -0400
parents 3a1ce69571e5
children 99cda2ff12b8
comparison
equal deleted inserted replaced
3:e9c34ef993d6 4:dfa2e75da6aa
1 #!/usr/bin/env python 1 #!/usr/bin/env python
2 # New in this version: 2 # New in this version:
3 # - handle commented-out header lines 3 # - option to generate report is triggered simply by including report filename
4 # - did everything through the Rpy2 interface instead of inline R code 4 # as third argument
5 5
6 import os,sys,numpy 6 import os,sys,numpy
7 from rpy2.robjects import Formula 7 from rpy2.robjects import Formula
8 from rpy2.robjects.packages import importr 8 from rpy2.robjects.packages import importr
9 from rpy2 import robjects 9 from rpy2 import robjects
11 def fail(message): 11 def fail(message):
12 sys.stderr.write(message+'\n') 12 sys.stderr.write(message+'\n')
13 sys.exit(1) 13 sys.exit(1)
14 14
15 args = sys.argv[1:] 15 args = sys.argv[1:]
16 if '-r' in args:
17 make_report = True
18 else:
19 make_report = False
20 if len(args) >= 1: 16 if len(args) >= 1:
21 infile = args[0] 17 infile = args[0]
22 else: 18 else:
23 fail('Error: No input filename provided (as argument 1).') 19 fail('Error: No input filename provided (as argument 1).')
24 if len(args) >= 2: 20 if len(args) >= 2:
25 outfile = args[1] 21 outfile = args[1]
26 else: 22 else:
27 fail('Error: No output filename provided (as argument 2).') 23 fail('Error: No output filename provided (as argument 2).')
28 if len(args) >= 3: 24 if len(args) >= 3:
29 report = args[2] 25 report = args[2]
30 elif make_report: 26 else:
31 fail('Error: No output report filename provided (as argument 3).') 27 report = ''
32 28
33 # Check input file 29 # Check input file
34 if not os.path.exists(infile): 30 if not os.path.exists(infile):
35 fail('Error: Input file '+infile+' could not be found.') 31 fail('Error: Input file '+infile+' could not be found.')
36 with open(infile, 'r') as lines: 32 with open(infile, 'r') as lines:
87 graphlabels = base.names(table) 83 graphlabels = base.names(table)
88 kwargs3 = {'pos':"0", 'las':"2", 'cex.axis':"1"} 84 kwargs3 = {'pos':"0", 'las':"2", 'cex.axis':"1"}
89 graphics.axis(1, at=range(1, len(graphlabels)+1, 1), labels=graphlabels, **kwargs3) 85 graphics.axis(1, at=range(1, len(graphlabels)+1, 1), labels=graphlabels, **kwargs3)
90 grdevices.dev_off() 86 grdevices.dev_off()
91 87
92 if not make_report: 88 if not report:
93 sys.exit(0) 89 sys.exit(0)
94 90
95 91
96 #################################### 92 ####################################
97 # GENERATE REPORT 93 # GENERATE REPORT
104 # MAD <= 2.0 fail 100 # MAD <= 2.0 fail
105 # MAD > 2.0 pass 101 # MAD > 2.0 pass
106 ################################### 102 ###################################
107 103
108 SAMPLES=[] 104 SAMPLES=[]
109 for i in range(len(tab)): 105 for i in range(len(table)):
110 SAMPLES.append(base.names(tab)[i]) 106 SAMPLES.append(base.names(table)[i])
111 107
112 def boxstats(data,sample): 108 def boxstats(data,sample):
113 VALUES = [100*float(x.strip().split('\t')[11]) for x in list(open(data)) if x.strip().split('\t')[0]==sample] 109 VALUES = [100*float(x.strip().split('\t')[11]) for x in list(open(data)) if x.strip().split('\t')[0]==sample]
114 NoHET = len(VALUES) 110 NoHET = len(VALUES)
115 MEDIAN = numpy.median(VALUES) 111 MEDIAN = numpy.median(VALUES)