|
0
|
1 #!/usr/bin/env python
|
|
|
2 # Code by Boris Rebolledo-Jaramillo
|
|
|
3 # (boris-at-bx.psu.edu)
|
|
|
4 # Edited by Nick Stoler
|
|
|
5 # (nick-at-bx.psu.edu)
|
|
|
6 # New in this version:
|
|
|
7 # - Add in proper header line if not present
|
|
|
8
|
|
|
9 import os
|
|
|
10 import sys
|
|
2
|
11 import array
|
|
0
|
12 import numpy
|
|
|
13 from rpy2.robjects import Formula
|
|
|
14 from rpy2.robjects.packages import importr
|
|
|
15 from rpy2 import robjects
|
|
|
16
|
|
|
17 def fail(message):
|
|
|
18 sys.stderr.write(message+'\n')
|
|
|
19 sys.exit(1)
|
|
|
20
|
|
|
21 COLUMN_LABELS = ['SAMPLE', 'CHR', 'POS', 'A', 'C', 'G', 'T', 'CVRG',
|
|
|
22 'ALLELES', 'MAJOR', 'MINOR', 'MINOR.FREQ.PERC.'] #, 'STRAND.BIAS']
|
|
|
23
|
|
|
24 args = sys.argv[1:]
|
|
|
25 if len(args) >= 1:
|
|
|
26 infile = args[0]
|
|
|
27 else:
|
|
|
28 fail('Error: No input filename provided (as argument 1).')
|
|
|
29 if len(args) >= 2:
|
|
|
30 outfile = args[1]
|
|
|
31 else:
|
|
|
32 fail('Error: No output filename provided (as argument 2).')
|
|
|
33 if len(args) >= 3:
|
|
|
34 report = args[2]
|
|
|
35 else:
|
|
|
36 report = ''
|
|
|
37
|
|
|
38 # Check input file
|
|
|
39 add_header = False
|
|
|
40 if not os.path.exists(infile):
|
|
|
41 fail('Error: Input file '+infile+' could not be found.')
|
|
|
42 with open(infile, 'r') as lines:
|
|
|
43 line = lines.readline()
|
|
|
44 if not line:
|
|
|
45 fail('Error: Input file seems to be empty')
|
|
|
46 line = line.strip().lstrip('#') # rm whitespace, comment chars
|
|
|
47 labels = line.split("\t")
|
|
|
48 if 'SAMPLE' not in labels or labels[11] != 'MINOR.FREQ.PERC.':
|
|
|
49 sys.stderr.write("Error: Input file does not seem to have a proper header "
|
|
|
50 +"line.\nAdding an artificial header..")
|
|
|
51 add_header = True
|
|
|
52
|
|
2
|
53 r = robjects.r
|
|
0
|
54 base = importr('base')
|
|
|
55 utils = importr('utils')
|
|
|
56 stats = importr('stats')
|
|
|
57 rprint = robjects.globalenv.get("print")
|
|
|
58 graphics = importr('graphics')
|
|
|
59 grdevices = importr('grDevices')
|
|
|
60 grdevices.png(file=outfile, width=1024, height=768)
|
|
|
61
|
|
|
62 # Read file into a data frame
|
|
|
63 if add_header:
|
|
|
64 # add header line manually if not present
|
|
|
65 DATA = utils.read_delim(infile, header=False)
|
|
|
66 labels = robjects.r.names(DATA)
|
|
|
67 for i in range(len(labels)):
|
|
|
68 try:
|
|
|
69 labels[i] = COLUMN_LABELS[i]
|
|
|
70 except IndexError, e:
|
|
|
71 fail("Error in input file: Too many columns (does not match hardcoded "
|
|
|
72 +"column labels).")
|
|
|
73 else:
|
|
|
74 DATA = utils.read_delim(infile)
|
|
|
75 # Remove comment from header, if present
|
|
|
76 labels = robjects.r.names(DATA)
|
|
|
77 if labels[0][0:2] == 'X.':
|
|
|
78 labels[0] = labels[0][2:]
|
|
|
79
|
|
|
80 # Multiply minor allele frequencies by 100 to get percentage
|
|
|
81 # .rx2() looks up a column by its label and returns it as a vector
|
|
|
82 # .ro turns the returned object into one that can be operated on per-element
|
|
|
83 minor_freq = DATA.rx2('MINOR.FREQ.PERC.').ro * 100
|
|
|
84 samples = DATA.rx2('SAMPLE')
|
|
|
85
|
|
|
86 # Formula() creates a Python object representing the R object returned by x ~ y
|
|
|
87 formula = Formula('minor_freq ~ samples')
|
|
|
88 # The "environment" in .getenvironment() is the entire R workspace in which the
|
|
|
89 # Formula object exists. The R workspace meaning all the defined variables.
|
|
|
90 # Here, the .getenvironment() method is being used to set some variables in the
|
|
2
|
91 # R workspace
|
|
0
|
92
|
|
|
93 formula.getenvironment()['minor_freq'] = minor_freq
|
|
|
94 formula.getenvironment()['samples'] = samples
|
|
|
95
|
|
2
|
96
|
|
|
97 r.par(oma=array.array('i', [0,0,0,0]))
|
|
|
98 r.par(mar=array.array('i', [10,4,4,2]))
|
|
|
99 ylimit = array.array('i',[-5,50])
|
|
|
100
|
|
0
|
101 # create boxplot - fill kwargs1 with the options for the boxplot function
|
|
|
102 kwargs1 = {'ylab':"Minor allele frequency (%)", 'col':"gray", 'xaxt':"n",
|
|
|
103 'outpch':"*",'main':"Distribution of minor allele frequencies",
|
|
|
104 'cex.lab':"1.5"}
|
|
2
|
105 p = graphics.boxplot(formula, axes=0,ylim=ylimit, lty=1,**kwargs1)
|
|
|
106
|
|
|
107 table = base.table(DATA.rx2('SAMPLE'))
|
|
|
108 #graphics.text(-1, 1, 'N:', font=2)
|
|
0
|
109 for i in range(1, base.length(table)[0]+1, 1):
|
|
2
|
110 graphics.text(i, -1, table[i-1], font=2)
|
|
0
|
111
|
|
|
112 graphlabels = base.names(table)
|
|
2
|
113 kwargs3 = {'pos':"-2", 'las':"2", 'cex.axis':"1"}
|
|
|
114 graphics.axis(1, at=range(1, len(graphlabels)+1, 1),labels=graphlabels, **kwargs3)
|
|
|
115 graphics.axis(2,at=(range(0,60,10)),pos=0,font=2)
|
|
0
|
116 grdevices.dev_off()
|
|
|
117
|
|
|
118 if not report:
|
|
|
119 sys.exit(0)
|
|
|
120
|
|
|
121
|
|
|
122 ####################################
|
|
|
123 # GENERATE REPORT
|
|
|
124 # report should be something like:
|
|
|
125 # SAMPLE NoHET MEDIAN MAD TEST
|
|
|
126 # s1 7 10% n p/w/f
|
|
|
127 # n <= 5 pass
|
|
|
128 # 6 <= n <=10 warn
|
|
|
129 # n >= 11 fail
|
|
|
130 # MAD <= 2.0 fail
|
|
|
131 # MAD > 2.0 pass
|
|
|
132 ###################################
|
|
|
133
|
|
|
134 SAMPLES=[]
|
|
|
135 for i in range(len(table)):
|
|
|
136 SAMPLES.append(base.names(table)[i])
|
|
|
137
|
|
|
138 def boxstats(data,sample):
|
|
|
139 VALUES = [100*float(x.strip().split('\t')[11]) for x in list(open(data)) if x.strip().split('\t')[0]==sample]
|
|
|
140 NoHET = len(VALUES)
|
|
|
141 MEDIAN = numpy.median(VALUES)
|
|
|
142 MAD = numpy.median([abs(i - MEDIAN) for i in VALUES]) # Median absolute distance (robust spread statistic)
|
|
|
143 return [NoHET,MEDIAN, MAD]
|
|
|
144
|
|
|
145 boxreport = open(report, "w+")
|
|
|
146 boxreport.write("SAMPLE\tTOTAL.SITES\tMEDIAN.FREQ.\tMAD.FREQ\tEVAL\n")
|
|
|
147 for sample in SAMPLES:
|
|
|
148 ENTRY = [sample] + boxstats(infile,sample)
|
|
|
149 if ENTRY[1] <= 5:
|
|
|
150 ENTRY.append('pass')
|
|
|
151 elif 6 <= ENTRY[1] <=10:
|
|
|
152 ENTRY.append('warn')
|
|
|
153 elif ENTRY[1] >= 11:
|
|
|
154 ENTRY.append('fail')
|
|
|
155 if ENTRY[3] <=2.0:
|
|
|
156 ENTRY.append('fail')
|
|
|
157 elif ENTRY[3] >2.0:
|
|
|
158 ENTRY.append('pass')
|
|
|
159 if len(set(ENTRY[4:6])) == 2:
|
|
|
160 ENTRY.append('warn')
|
|
|
161 else:
|
|
|
162 ENTRY.append(list(set(ENTRY[4:6]))[0])
|
|
|
163 boxreport.write ('%s\t%d\t%.1f\t%.1f\t%s\n' % tuple([ENTRY[i] for i in [0,1,2,3,6]]))
|
|
|
164
|
|
|
165 boxreport.close()
|
|
|
166
|
|
|
167
|
|
|
168
|
|
|
169
|
