htsworkflow/pipelines/genome_mapper.py

   1 #!/usr/bin/python
   2 import glob
   3 import sys
   4 import os
   5 import re
   6
   7 import logging
   8
   9 from htsworkflow.util.alphanum import alphanum
  10
  11 class DuplicateGenome(Exception): pass
  12
  13
  14 def _has_metainfo(genome_dir):
  15   metapath = os.path.join(genome_dir, '_metainfo_')
  16   if os.path.isfile(metapath):
  17     return True
  18   else:
  19     return False
  20
  21 def getAvailableGenomes(genome_base_dir):
  22   """
  23   raises IOError (on genome_base_dir not found)
  24   raises DuplicateGenome on duplicate genomes found.
  25
  26   returns a double dictionary (i.e. d[species][build] = path)
  27   """
  28
  29   # Need valid directory
  30   if not os.path.exists(genome_base_dir):
  31     msg = "Directory does not exist: %s" % (genome_base_dir)
  32     raise IOError, msg
  33
  34   # Find all subdirectories
  35   filepath_list = glob.glob(os.path.join(genome_base_dir, '*'))
  36   potential_genome_dirs = \
  37     [ filepath for filepath in filepath_list if os.path.isdir(filepath)]
  38
  39   # Get list of metadata files
  40   genome_dir_list = \
  41     [ dirpath \
  42       for dirpath in potential_genome_dirs \
  43       if _has_metainfo(dirpath) ]
  44
  45   # Genome double dictionary
  46   d = {}
  47
  48   for genome_dir in genome_dir_list:
  49     line = open(os.path.join(genome_dir, '_metainfo_'), 'r').readline().strip()
  50
  51     # Get species, build... log and skip on failure
  52     try:
  53       species, build = line.split('|')
  54     except:
  55       logging.warning('Skipping: Invalid metafile (%s) line: %s' \
  56                       % (metafile, line))
  57       continue
  58
  59     build_dict = d.setdefault(species, {})
  60     if build in build_dict:
  61       msg = "Duplicate genome for %s|%s" % (species, build)
  62       raise DuplicateGenome, msg
  63
  64     build_dict[build] = genome_dir
  65
  66   return d
  67
  68
  69 class constructMapperDict(object):
  70     """
  71     Emulate a dictionary to map genome|build names to paths.
  72
  73     It uses the dictionary generated by getAvailableGenomes.
  74     """
  75     def __init__(self, genome_dict):
  76         self.genome_dict = genome_dict
  77
  78     def __getitem__(self, key):
  79         """
  80         Return the best match for key
  81         """
  82         elements = re.split("\|", key)
  83
  84         try:
  85             if len(elements) == 1:
  86                 # we just the species name
  87                 # get the set of builds
  88                 builds = self.genome_dict[elements[0]]
  89
  90                 # sort build names the way humans would
  91                 keys = builds.keys()
  92                 keys.sort(cmp=alphanum)
  93
  94                 # return the path from the 'last' build name
  95                 return builds[keys[-1]]
  96
  97             elif len(elements) == 2:
  98                 # we have species, and build name
  99                 return self.genome_dict[elements[0]][elements[1]]
 100             else:
 101                 raise KeyError("Unrecognized key")
 102         except KeyError, e:
 103             logging.error('Unrecognized genome identifier: %s' % str((elements),))
 104             return "NoGenomeAvailable"
 105
 106     def keys(self):
 107         keys = []
 108         for species in self.genome_dict.keys():
 109             for build in self.genome_dict[species]:
 110                 keys.append([species+'|'+build])
 111         return keys
 112
 113     def values(self):
 114         values = []
 115         for species in self.genome_dict.keys():
 116             for build in self.genome_dict[species]:
 117                 values.append(self.genome_dict[species][build])
 118         return values
 119
 120     def items(self):
 121         items = []
 122         for species in self.genome_dict.keys():
 123             for build in self.genome_dict[species]:
 124                 key = [species+'|'+build]
 125                 value = self.genome_dict[species][build]
 126                 items.append((key, value))
 127         return items
 128
 129 if __name__ == '__main__':
 130
 131   if len(sys.argv) != 2:
 132     print 'useage: %s <base_genome_dir>' % (sys.argv[0])
 133     sys.exit(1)
 134
 135   d = getAvailableGenomes(sys.argv[1])
 136   d2 = constructMapperDict(d)
 137
 138   for k,v in d2.items():
 139     print '%s: %s' % (k,v)
 140
 141