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     print "getting uniq reads"
  89     uniqDict = RDS.getReadsDict(withChrom=True, withPairID=True, doUniqs=True, readIDDict=True)
  90     print "got uniqs"
  91
  92     if doSplices:
  93         addSplicesToUniqReads(RDS, uniqDict)
  94
  95     if doVerbose:
  96         print len(uniqDict), time.ctime()
  97
  98     outfile = open(outfilename,"w")
  99     diffChrom = 0
 100     distal = 0
 101     total = 0
 102     for readID in uniqDict:
 103         readList = uniqDict[readID]
 104         if len(readList) == 2:
 105             total += 1
 106             start1 = readList[0]["start"]
 107             sense1 = readList[0]["sense"]
 108             chrom1 = readList[0]["chrom"]
 109             start2 = readList[1]["start"]
 110             sense2 = readList[1]["sense"]
 111             chrom2 = readList[1]["chrom"]
 112
 113             if chrom1 != chrom2:
 114                 diffChrom += 1
 115                 if sameChromOnly:
 116                     continue
 117                 else:
 118                     outline = "%s\t%s\t%d\t%s\t%s\t%d\t%s" % (readID, chrom1, start1, sense1, chrom2, start2, sense2)
 119                     print >> outfile, outline
 120                     if doVerbose:
 121                         print diffChrom, outline
 122             else:
 123                 dist = abs(start1 - start2)
 124                 if minDist < dist < maxDist:
 125                     distal += 1
 126                     outline = "%s\t%s\t%d\t%s\t%d\t%s\t%d" % (readID, chrom1, start1, sense1, start2, sense2, dist)
 127                     print >> outfile, outline
 128                     if doVerbose:
 129                         print distal, outline
 130
 131     outfile.write("#distal: %d\tdiffChrom: %d\tpossible: %d\n" % (distal, diffChrom, total))
 132     total = float(total)
 133     if total < 1:
 134         total = 1.
 135
 136     outfile.write("#distal %2.2f pct\tdiffChrom %2.2f pct\n" % ((100. * distal/total), (100. * diffChrom/total)))
 137     outfile.close()
 138     print "distal: %d\tdiffChrom: %d\tpossible: %d" % (distal, diffChrom, int(total))
 139     print "distal: %2.2f pct\tdiffChrom: %2.2f pct\n" % ((100. * distal/total), (100. * diffChrom/total))
 140     print time.ctime()
 141
 142
 143 def addSplicesToUniqReads(RDS, uniqDict):
 144     print "getting splices"
 145     splicesDict = RDS.getSplicesDict(withChrom=True, withPairID=True, readIDDict=True, splitRead=True)
 146     print "got splices"
 147     for readID in splicesDict:
 148         theRead = splicesDict[readID]
 149         read0 = theRead[0]
 150         read0["start"] = read0["startL"]
 151         del read0["stopL"]
 152         try:
 153             uniqDict[readID].append(read0)
 154         except:
 155             if len(theRead) == 4:
 156                 read2 = theRead[2]
 157                 read2["start"] = read2["startL"]
 158                 del read2["stopL"]
 159                 uniqDict[readID] = [read0,read2]
 160
 161
 162 if __name__ == "__main__":
 163     main(sys.argv)