rnaEditing.py

   1 import sys
   2 import optparse
   3 import chksnp, getSNPs, getSNPGeneInfo, analyzego, getNovelSNPs, makeSNPtrack, rnaAToIFilter
   4 from commoncode import countDuplicatesInList, getConfigParser, getConfigOption
   5
   6
   7 def main(argv=None):
   8     if not argv:
   9         argv = sys.argv
  10
  11     usage = "usage: python %prog dbfile snpsfile genome rpkmfile [options]"
  12
  13     parser = getParser(usage)
  14     (options, args) = parser.parse_args(argv[1:])
  15
  16     if len(args) < 4:
  17         print usage
  18         sys.exit(1)
  19
  20     dbfile = args[0]
  21     hitfile = args[1]
  22     genome = args[2]
  23     rpkmfilename = args[3]
  24
  25     rnaEditing(dbfile, hitfile, genome, rpkmfilename, options)
  26
  27
  28 def getParser(usage):
  29     parser = optparse.OptionParser(usage=usage)
  30     parser.add_option("--goprefix", dest="prefix")
  31     parser.add_option("--novelsnp", dest="novelsnpoutfilename")
  32     parser.add_option("--bedfile", dest="bedoutfilename")
  33     parser.add_option("--cache", type="int", dest="cachePages")
  34     parser.add_option("--snpDB", action="append", dest="snpDBList",
  35                       help="additional snp db files to check will be searched in order given")
  36
  37     configParser = getConfigParser()
  38     section = "rnaEditing"
  39     prefix = getConfigOption(configParser, section, "prefix", None)
  40     novelsnpoutfilename = getConfigOption(configParser, section, "novelsnpoutfilename", None)
  41     bedoutfilename = getConfigOption(configParser, section, "bedoutfilename", None)
  42     cachePages = getConfigOption(configParser, section, "cachePages", None)
  43
  44     parser.set_defaults(prefix=prefix, novelsnpoutfilename=novelsnpoutfilename, bedoutfilename=bedoutfilename,
  45                         cachePages=cachePages, snpDBList=[])
  46
  47     return parser
  48
  49
  50 def rnaEditing(dbfile, hitfile, genome, rpkmfilename, options):
  51     if options.cachePages is not None:
  52         doCache = True
  53     else:
  54         doCache = False
  55
  56     # get the SNPs
  57     snpList = getSNPs.getSNPs(hitfile, 3, 0.25, doCache, options.cachePages, forceChr=True)
  58
  59     # check for existing SNPs
  60     dbList = [dbfile]
  61     for dbFileName in options.snpDBList:
  62         dbList.append(dbFileName)
  63
  64     snpPropertiesList = chksnp.chkSNP(dbList, snpList, options.cachePages)
  65
  66     # get the neighboring genes
  67     geneInfoList = getSNPGeneInfo.getSNPGeneInfo(genome, snpPropertiesList, rpkmfilename, doCache, flankBP=10000)
  68
  69     # filter out for the A-to-I events in the same direction as the genes
  70     filteredSNPs = rnaAToIFilter.rnaAToIFilter(geneInfoList)
  71
  72     # count the number of different bases that have been called for each gene
  73     # pick a set of genes with a high number of sites (here 5)
  74     geneList = getGenesWithMultipleSNPs(filteredSNPs, minCount=5)
  75
  76     if options.prefix is not None:
  77         analyzego.analyzeGO(genome, geneList, options.prefix, translateGene=True, fieldID=1)
  78
  79     if options.novelsnpoutfilename is not None:
  80         getNovelSNPs.writeNovelSNPFile(genome, filteredSNPs, options.novelsnpoutfilename)
  81
  82     if options.bedoutfilename is not None:
  83         makeSNPtrack.writeSNPsBedfile(filteredSNPs, "rnaEdit_sample", options.bedoutfilename)
  84
  85
  86 def getGenesWithMultipleSNPs(snpList, minCount=1):
  87     geneList = []
  88     for snpEntry in snpList:
  89         geneList.append(snpEntry[11])
  90
  91     duplicateCountList = countDuplicatesInList(geneList)
  92
  93     geneList = []
  94     for (gene, count) in duplicateCountList:
  95         if count >= minCount:
  96             geneList.append(gene)
  97
  98     return geneList
  99
 100
 101 if __name__ == "__main__":
 102     main(sys.argv)