distalPairs.py

   1 #
   2 #  distalPairs.py
   3 #  ENRAGE
   4 #
   5 #  Created by Ali Mortazavi on 10/14/08.
   6 #
   7
   8
   9 try:
  10     import psyco
  11     psyco.full()
  12 except:
  13     pass
  14
  15 import sys
  16 import time
  17 import optparse
  18 import ReadDataset
  19 from commoncode import getConfigParser, getConfigOption, getConfigIntOption, getConfigBoolOption
  20
  21
  22 def main(argv=None):
  23     if not argv:
  24         argv = sys.argv
  25
  26     print "distalPairs: version 3.4"
  27     print "looks at all chromosomes simultaneously: is both slow and takes up large amount of RAM"
  28     usage = "usage: python %prog minDist rdsfile outfile [--sameChrom] [--splices] [--maxDist bp] [--verbose] [--cache cachepages]"
  29
  30     parser = optparse.OptionParser(usage=usage)
  31     parser.add_option("--sameChrom", action="store_true", dest="sameChromOnly")
  32     parser.add_option("--splices", action="store_true", dest="doSplices")
  33     parser.add_option("--verbose", action="store_true", dest="doVerbose")
  34     parser.add_option("--maxDist", type="int", dest="maxDist")
  35     parser.add_option("--cache", type="int", dest="cachePages")
  36     parser.set_defaults(sameChromOnly=False, doSplices=False, doVerbose=False, maxDist=1000000000, cachePages=None)
  37     (options, args) = parser.parse_args(argv[1:])
  38
  39     if len(args) < 3:
  40         print usage
  41         sys.exit(1)
  42
  43     minDist = int(args[0])
  44     rdsfile = args[1]
  45     outfilename = args[2]
  46
  47     distalPairs(minDist, rdsfile, outfilename, options.sameChromOnly, options.doSplices, options.doVerbose, options.maxDist, options.cachePages)
  48
  49
  50 def makeParser(usage=""):
  51     parser = optparse.OptionParser(usage=usage)
  52     parser.add_option("--sameChrom", action="store_true", dest="sameChromOnly")
  53     parser.add_option("--splices", action="store_true", dest="doSplices")
  54     parser.add_option("--verbose", action="store_true", dest="doVerbose")
  55     parser.add_option("--maxDist", type="int", dest="maxDist")
  56     parser.add_option("--cache", type="int", dest="cachePages")
  57
  58     configParser = getConfigParser()
  59     section = "distalPairs"
  60     sameChromOnly = getConfigBoolOption(configParser, section, "sameChromOnly", False)
  61     doSplices = getConfigBoolOption(configParser, section, "doSplices", False)
  62     doVerbose = getConfigBoolOption(configParser, section, "doVerbose", False)
  63     maxDist = getConfigIntOption(configParser, section, "maxDist", 1000000000)
  64     cachePages = getConfigOption(configParser, section, "cachePages", None)
  65
  66     parser.set_defaults(sameChromOnly=sameChromOnly, doSplices=doSplices, doVerbose=doVerbose, maxDist=maxDist, cachePages=cachePages)
  67
  68     return parser
  69
  70
  71 def distalPairs(minDist, rdsfile, outfilename, sameChromOnly=False, doSplices=False, doVerbose=False, maxDist=1000000000, cachePages=None):
  72     if cachePages is not None:
  73         doCache = True
  74     else:
  75         doCache = False
  76         cachePages = -1
  77
  78     RDS = ReadDataset.ReadDataset(rdsfile, verbose = True, cache=doCache)
  79     if not RDS.hasIndex():
  80         print "Will not attempt to run on unIndexed dataset - please index with rdsmetadata.py and rerun"
  81         sys.exit(1)
  82
  83     if cachePages > RDS.getDefaultCacheSize():
  84         RDS.setDBcache(cachePages)
  85
  86     print time.ctime()
  87
  88     if doSplices:
  89         print "getting splices"
  90         splicesDict = RDS.getSplicesDict(withChrom=True, withPairID=True, readIDDict=True, splitRead=True)
  91         print "got splices"
  92
  93     print "getting uniq reads"
  94     uniqDict = RDS.getReadsDict(withChrom=True, withPairID=True, doUniqs=True, readIDDict=True)
  95     print "got uniqs"
  96
  97     if doSplices:
  98         for readID in splicesDict:
  99             theRead = splicesDict[readID]
 100             read0 = theRead[0]
 101             del read0[1]
 102             try:
 103                 uniqDict[readID].append(read0)
 104             except:
 105                 if len(theRead) == 4:
 106                     read2 = theRead[2]
 107                     del read2[1]
 108                     uniqDict[readID] = [read0,read2]
 109
 110     if doVerbose:
 111         print len(uniqDict), time.ctime()
 112
 113     outfile = open(outfilename,"w")
 114
 115     diffChrom = 0
 116     distal = 0
 117     total = 0
 118     for readID in uniqDict:
 119         readList = uniqDict[readID]
 120         if len(readList) == 2:
 121             total += 1
 122             start1 = readList[0]["start"]
 123             sense1 = readList[0]["sense"]
 124             chrom1 = readList[0]["chrom"]
 125             start2 = readList[1]["start"]
 126             sense2 = readList[1]["sense"]
 127             chrom2 = readList[1]["chrom"]
 128
 129             if chrom1 != chrom2:
 130                 diffChrom += 1
 131                 if sameChromOnly:
 132                     continue
 133                 else:
 134                     outline = "%s\t%s\t%d\t%s\t%s\t%d\t%s" % (readID, chrom1, start1, sense1, chrom2, start2, sense2)
 135                     outfile.write(outline + "\n")
 136                     if doVerbose:
 137                         print diffChrom, outline
 138             else:
 139                 dist = abs(start1 - start2)
 140
 141                 if minDist < dist < maxDist:
 142                     distal += 1
 143                     outline = "%s\t%s\t%d\t%s\t%d\t%s\t%d" % (readID, chrom1, start1, sense1, start2, sense2, dist)
 144                     outfile.write(outline + "\n")
 145                     if doVerbose:
 146                         print distal, outline
 147
 148     outfile.write("#distal: %d\tdiffChrom: %d\tpossible: %d\n" % (distal, diffChrom, total))
 149     total = float(total)
 150     if total < 1:
 151         total = 1.
 152
 153     outfile.write("#distal %2.2f pct\tdiffChrom %2.2f pct\n" % ((100. * distal/total), (100. * diffChrom/total)))
 154     outfile.close()
 155     print "distal: %d\tdiffChrom: %d\tpossible: %d" % (distal, diffChrom, int(total))
 156     print "distal: %2.2f pct\tdiffChrom: %2.2f pct\n" % ((100. * distal/total), (100. * diffChrom/total))
 157     print time.ctime()
 158
 159
 160 if __name__ == "__main__":
 161     main(sys.argv)