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