solexa2srf likes to produce output, so my trick of watching the
[htsworkflow.git] / gaworkflow / pipeline / 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 gaworkflow.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         if len(elements) == 1:
85             # we just the species name
86             # get the set of builds
87             builds = self.genome_dict[elements[0]]
88             
89             # sort build names the way humans would
90             keys = builds.keys()
91             keys.sort(cmp=alphanum)
92             
93             # return the path from the 'last' build name
94             return builds[keys[-1]]
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 keys(self):
103         keys = []
104         for species in self.genome_dict.keys():
105             for build in self.genome_dict[species]:
106                 keys.append([species+'|'+build])
107         return keys
108             
109     def values(self):
110         values = []
111         for species in self.genome_dict.keys():
112             for build in self.genome_dict[species]:
113                 values.append(self.genome_dict[species][build])
114         return values
115        
116     def items(self):
117         items = []
118         for species in self.genome_dict.keys():
119             for build in self.genome_dict[species]:
120                 key = [species+'|'+build]
121                 value = self.genome_dict[species][build]
122                 items.append((key, value))
123         return items
124             
125 if __name__ == '__main__':
126
127   if len(sys.argv) != 2:
128     print 'useage: %s <base_genome_dir>' % (sys.argv[0])
129     sys.exit(1)
130
131   d = getAvailableGenomes(sys.argv[1])
132   d2 = constructMapperDict(d)
133
134   for k,v in d2.items():
135     print '%s: %s' % (k,v)
136   
137