14 from commoncode import getConfigParser, getConfigOption, getConfigBoolOption, getConfigFloatOption
16 print "profilebins: version 2.3"
23 usage = "usage: python %prog label infile1 [--upstream infile2] [--downstream infile3] [--uplength kb] [--downlength kb] [--gene geneName] [--genes genefile] [--append] outfile"
25 parser = makeParser(usage)
26 (options, args) = parser.parse_args(argv[1:])
36 profilebins(label, infilename, outfilename, options.upfilename, options.downfilename,
37 options.uplength, options.downlength, options.gene, options.genefile,
41 def makeParser(usage=""):
42 parser = optparse.OptionParser(usage=usage)
43 parser.add_option("--upstream", dest="upfilename")
44 parser.add_option("--downstream", dest="downfilename")
45 parser.add_option("--uplength", type="float", dest="uplength")
46 parser.add_option("--downlength", type="int", dest="")
47 parser.add_option("--gene", dest="gene")
48 parser.add_option("--genes", dest="genefile")
49 parser.add_option("--append", action="store_true", dest="doAppend")
51 configParser = getConfigParser()
52 section = "profilebins"
53 upfilename = getConfigOption(configParser, section, "upfilename", None)
54 downfilename = getConfigOption(configParser, section, "downfilename", None)
55 uplength = getConfigFloatOption(configParser, section, "uplength", 0.0)
56 downlength = getConfigFloatOption(configParser, section, "downlength", 0.0)
57 gene = getConfigOption(configParser, section, "gene", None)
58 genefile = getConfigOption(configParser, section, "genefile", None)
59 doAppend = getConfigBoolOption(configParser, section, "doAppend", False)
61 parser.set_defaults(upfilename=upfilename, downfilename=downfilename, uplength=uplength, downlength=downlength,
62 gene=gene, genefile=genefile, doAppend=doAppend)
67 def profilebins(label, infilename, outfilename, upfilename=None, downfilename=None,
68 uplength=0.0, downlength=0.0, gene=None, genefile=None, doAppend=False):
70 fileList = [infilename]
77 if genefile is not None:
79 fields = line.strip().split()
81 geneList.append(fields[0])
83 geneList.append(line.strip())
87 if upfilename is not None:
88 fileList = [upfilename, infilename]
90 if downfilename is not None:
91 fileList.append(downfilename)
97 partLength = [uplength, 10.]
98 partOffset = [-1. * uplength, 0.]
101 partLength.append(downlength)
102 partOffset.append(10.)
106 for afile in fileList:
109 line = infile.readline()
110 fields = line.strip().split()
111 numBins = len(fields) - 4
114 weight = float(fields[2])
115 if restrictGenes and geneName in geneList:
116 totalWeight += weight
119 for myBin in fields[4:]:
120 if not restrictGenes or (restrictGenes and geneName in geneList):
121 totalBins[-1].append(weight * float(myBin))
123 totalBins[-1].append(0.)
126 fields = line.strip().split()
128 if restrictGenes and geneName not in geneList:
131 weight = float(fields[2])
133 for myBin in fields[4:]:
134 totalBins[-1][index] += weight * float(myBin)
137 totalWeight += weight
142 outfile = open(outfilename, "a")
144 outfile = open(outfilename, "w")
145 outfile.write("x-axis")
147 for partBins in totalBins:
148 partLen = partLength[partIndex]
149 numBins = len(partBins)
150 for binIndex in range(numBins):
151 outfile.write("\t%.2f" % (partOffset[partIndex] + (binIndex * partLen/numBins)))
155 outfile.write("\tweight\n")
158 for partBins in totalBins:
159 for aBin in partBins:
160 percent = aBin / totalWeight
161 outfile.write("\t%.1f" % percent)
163 totalPercent += percent
165 outfile.write("\t%.1f\n" % totalWeight)
172 if __name__ == "__main__":