buildMatrix.py

   1 #
   2 #  buildMatrix.py
   3 #  ENRAGE
   4 #
   5 #  Created by Ali Mortazavi on 3/6/09.
   6 #
   7 import sys, string, optparse
   8 from commoncode import writeLog
   9
  10 versionString = "%prog: version 1.3"
  11 print versionString
  12
  13
  14 def main(argv=None):
  15     if not argv:
  16         argv = sys.argv
  17
  18     usage = "usage: python %prog matrix.step.N-1 data.part matrix.step.N [--rescale] [--truncate maxRPKM] [--log altlogfile]"
  19
  20     parser = optparse.OptionParser(usage=usage)
  21     parser.add_option("--rescale", action="store_true", dest="rescale")
  22     parser.add_option("--truncate", type="int", dest="maxRPKM")
  23     parser.add_option("--log", dest="logfilename")
  24     parser.set_defaults(rescale=False, maxRPKM=None, logfilename="buildMatrix.log")
  25     (options, args) = parser.parse_args(argv[1:])
  26
  27     if len(args) < 3:
  28         print usage
  29         sys.exit(0)
  30
  31     infile = args[0]
  32     colfilename = args[1]
  33     outfilename = args[2]
  34
  35     if options.maxRPKM is not None:
  36         truncateRPKM = True
  37         maxRPKM = options.maxRPKM
  38     else:
  39         truncateRPKM = False
  40         maxRPKM = 100000000
  41
  42     buildMatrix(infile, colfilename, outfilename, truncateRPKM, maxRPKM,
  43                 options.rescale, options.logfilename)
  44
  45
  46 def buildMatrix(inFileName, colfilename, outfilename, truncateRPKM,
  47                 maxRPKM=100000000, rescale=False, logfilename="buildMatrix.log"):
  48
  49     writeLog(logfilename, versionString, string.join(sys.argv[1:]))
  50
  51     if "/" in colfilename:
  52         colname = colfilename.split("/")[-1]
  53     else:
  54         colname = colfilename
  55
  56     fileParts = colname.split(".")
  57     colID =  fileParts[0]
  58
  59     infile = open(inFileName)
  60     colfile = open(colfilename)
  61     outfile = open(outfilename, "w")
  62     header = infile.readline()[:-1]
  63     if header.strip() == "":
  64         header = "#\t"
  65
  66     outfile.write( "%s\t%s\n" % (header, colID))
  67
  68     values = []
  69     min = 20000000000.
  70     max = -1.
  71     untruncatedMax = -1.
  72     for line in colfile:
  73         if doNotProcessLine(line):
  74             continue
  75
  76         fields = line.strip().split()
  77         val = float(fields[-1])
  78         if truncateRPKM and val > maxRPKM:
  79             if val > untruncatedMax:
  80                 untruncatedMax = val
  81
  82             val = maxRPKM
  83
  84         values.append(val)
  85         if val < min:
  86             min = val
  87
  88         if val > max:
  89             max = val
  90
  91     range = max - min
  92     if rescale:
  93         finalValues = [(val - min)/range for val in values]
  94     else:
  95         finalValues = values
  96
  97     for val in finalValues:
  98         line = infile.readline().strip()
  99         line += "\t%1.3f\n" % val
 100         outfile.write(line)
 101
 102     outfile.close()
 103
 104     if untruncatedMax > 0:
 105         max = untruncatedMax
 106
 107     message = "max value in %s was %.2f" % (colname, max)
 108     if untruncatedMax > 0:
 109         message += " but was truncated to %d" % maxRPKM
 110
 111     print message
 112     writeLog(logfilename, versionString, message)
 113
 114
 115 def doNotProcessLine(line):
 116     return line[0] == "#"
 117
 118
 119 if __name__ == "__main__":
 120     main(sys.argv)