source: orange-bioinformatics/orangecontrib/bio/obiGEO.py @ 1873:0810c5708cc5

Revision 1873:0810c5708cc5, 14.3 KB checked in by Ales Erjavec <ales.erjavec@…>, 6 months ago (diff)

Moved '_bioinformatics' into orangecontrib namespace.

Line 
1from __future__ import absolute_import
2
3import cPickle, gzip, os.path, re
4
5import orange
6from Orange.orng import orngMisc, orngServerFiles
7
8from . import obiData, obiTaxonomy
9
10def spots_mean(x):
11    vs = [v for v in x if v and v<>"?"]
12    if len(vs) == 0: return "?"
13    return sum(vs)/len(vs)
14def spots_median(x):
15    vs = [v for v in x if v and v<>"?"]
16    if len(vs) == 0: return "?"
17    if len(vs) % 2:
18        return sorted(vs)/(len(vs)/2)
19    else:
20        z = sorted(x)
21        return (z[len(vs)/2-1] + z[len(vs)/2]) / 2. 
22    return sum(vs)/len(vs)
23def spots_min(x):
24    vs = [v for v in x if v and v<>"?"]
25    if len(vs) == 0: return "?"
26    return min(vs)/len(vs)
27def spots_max(x):
28    vs = [v for v in x if v and v<>"?"]
29    if len(vs) == 0: return "?"
30    return max(vs)/len(vs)
31
32p_assign = re.compile(" = (.*$)")
33p_tagvalue = re.compile("![a-z]*_([a-z_]*) = (.*)$")   
34tagvalue = lambda x: p_tagvalue.search(x).groups()
35
36DOMAIN = "GEO"
37GDS_INFO_FILENAME = "gds_info.pickled"
38FTP_NCBI = "ftp.ncbi.nih.gov"
39FTP_DIR = "pub/geo/DATA/SOFT/GDS/"
40
41class GDSInfo:
42    def __init__(self, force_update=False):
43        path = orngServerFiles.localpath(DOMAIN, GDS_INFO_FILENAME)
44        if not os.path.exists(path) or force_update:
45            orngServerFiles.download(DOMAIN, GDS_INFO_FILENAME)
46        f = file(path, "rb")
47        self.info, self.excluded = cPickle.load(f)
48    def keys(self): return self.info.keys()
49    def items(self): return self.info.items()
50    def values(self): return self.info.values()
51    def clear(self): return self.info.clear()
52    def __getitem__(self, key): return self.info[key]
53    def __setitem__(self, key, item): self.info[key] = item
54    def __len__(self): return len(self.info)
55    def __iter__(self): return iter(self.info)
56    def __contains__(self, key): return key in self.info
57   
58
59class GeneData:
60    """Store mapping between spot id and gene."""
61    def __init__(self, spot_id, gene_name, d):
62        self.spot_id = spot_id
63        self.gene_name = gene_name
64        self.data = d
65
66class GDS():
67    """GEO DataSet class: read GEO datasets and convert them to ExampleTable."""
68    def __init__(self, gdsname, verbose=False, force_download=False):
69        self.gdsname = gdsname
70        self.verbose = verbose
71        self.force_download = force_download
72        self.filename = orngServerFiles.localpath(DOMAIN, self.gdsname + ".soft.gz")
73        self._getinfo() # to get the info
74        taxid = obiTaxonomy.search(self.info["sample_organism"], exact=True)
75        self.info["taxid"] = taxid[0] if len(taxid)==1 else None
76        self._getspotmap() # to get gene->spot and spot->gene mapping
77        self.genes = self.gene2spots.keys()
78        self.info["gene_count"] = len(self.genes)
79        self.gdsdata = None
80        self.data = None
81       
82    def _download(self):
83        """Download GDS data file if not in local cache or forced download requested."""
84        localpath = orngServerFiles.localpath(DOMAIN)
85        if self.force_download or not os.path.exists(self.filename):
86            ftp = obiData.FtpDownloader(FTP_NCBI, localpath, FTP_DIR)
87            ftp.retrieve(self.gdsname + ".soft.gz", progressCallback=orngMisc.ConsoleProgressBar()
88                         if self.verbose else None)
89            if self.verbose:
90                # print "Downloading %s" % self.gdsname
91                print
92
93    def _getinfo(self):
94        """Parse GDS data file and return a dictionary with info."""
95        getstate = lambda x: x.split(" ")[0][1:] 
96        getid = lambda x: x.rstrip().split(" ")[2]
97        self._download()
98        f = gzip.open(self.filename)
99        state = None; previous_state = None
100   
101        info = {"subsets" : []}
102        subset = None
103   
104        # GDS information part
105        for line in f:
106            if line[0] == "^":
107                previous_state = state; state = getstate(line)
108                if state == "SUBSET":
109                    if subset:
110                        info["subsets"] += [subset]
111                    subset = {"id" : getid(line)}
112                if state == "DATASET":
113                    info["dataset_id"] = getid(line)
114                continue
115            if state == "DATASET":
116                if previous_state == "DATABASE":
117                    tag, value = tagvalue(line)
118                    info[tag] = value
119                else:
120                    if subset:
121                        info["subsets"] += [subset]
122                    break
123            if state == "SUBSET":
124                tag, value = tagvalue(line)
125                if tag == "description" or tag == "type":
126                    subset[tag] = value
127                if tag == "sample_id":
128                    subset[tag] = value.split(",")
129        for t,v in info.items():
130            if "count" in t:
131                info[t] = int(v)
132               
133        # sample information
134        state = None
135        for line in f:
136            if state == "header":
137                info["samples"] = line.rstrip().split("\t")[2:]
138                break
139            if line.startswith("!dataset_table_begin"):
140                state = "header"
141
142        self.info = info
143
144    def _getspotmap(self, include_spots=None):
145        """Return gene to spot and spot to genes mapings."""
146        f = gzip.open(self.filename)
147        for line in f:
148            if line.startswith("!dataset_table_begin"):
149                break
150        f.readline() # skip header
151        spot2gene = {}
152        gene2spots = {}
153        for line in f:
154            if line.startswith("!dataset_table_end"):
155                break
156            spot, gene = line.rstrip().split("\t")[0:2]
157            if include_spots and (spot not in include_spots):
158                continue 
159            spot2gene[spot] = gene
160            gene2spots[gene] = gene2spots.get(gene, []) + [spot]
161   
162        self.spot2gene = spot2gene
163        self.gene2spots = gene2spots
164       
165    def sample_annotations(self, sample_type=None):
166        """Return a dictionary with sample annotation."""
167        annotation = {}
168        for info in self.info["subsets"]:
169            if sample_type and info["type"]<>sample_type:
170                continue
171            for id in info["sample_id"]:
172                annotation.setdefault(id, {})[info["type"]]=info["description"]
173        return annotation
174
175    def sample_to_class(self, sample_type=None, missing_class_value=None):
176        """Return class values for GDS samples."""
177        annotations = self.sample_annotations(sample_type)
178        return dict([(sample, "|".join([a for t,a in ann.items()])) for sample, ann in annotations.items()])
179       
180    def sample_types(self):
181        """Return a set of sample types."""
182        return set([info["type"] for info in self.info["subsets"]])
183   
184    def _parse_soft(self, remove_unknown=None):
185        """Parse GDS data, returns data dictionary."""
186        f = gzip.open(self.filename)
187        mfloat = lambda x: float(x) if x<>'null' else '?'
188   
189        data = {}
190        # find the start of the data part
191        for line in f:
192            if line.startswith("!dataset_table_begin"):
193                break
194        f.readline()
195       
196        # load data
197        for line in f:
198            if line.startswith("!dataset_table_end"):
199                break
200            d = line.rstrip().split("\t")
201            if remove_unknown and (float(d[2:].count('null')) / len(d[2:]) > remove_unknown):
202#                print "ZZZ", float(d[2:].count('null')) / len(d[2:]), d[2:]
203                continue 
204            data[d[0]] = GeneData(d[0], d[1], [mfloat(v) for v in d[2:]])
205        self.gdsdata = data
206   
207    def _to_ExampleTable(self, report_genes=True, merge_function=spots_mean,
208                                sample_type=None, missing_class_value=None, transpose=False):
209        """Convert parsed GEO format to orange, save by genes or by spots."""
210        orng_data = []
211        if transpose: # samples in rows
212            sample2class = self.sample_to_class(sample_type, missing_class_value)
213            cvalues = list(set(sample2class.values()))
214            if None in cvalues:
215                cvalues.remove(None)
216            classvar = orange.EnumVariable(name=sample_type or "class", values=cvalues)
217            if report_genes: # save by genes
218                atts = [orange.FloatVariable(name=gene) for gene in self.gene2spots.keys()]
219                domain = orange.Domain(atts, classvar)
220                for (i, sampleid) in enumerate(self.info["samples"]):
221                    vals = [merge_function([self.gdsdata[spot].data[i] \
222                            for spot in self.gene2spots[gene]]) for gene in self.gene2spots.keys()]
223                    orng_data.append(vals + [sample2class[sampleid]])
224               
225            else: # save by spots
226                spots = self.spot2gene.keys()
227                atts = [orange.FloatVariable(name=id) for id in spots]
228                domain = orange.Domain(atts, classvar)
229                for (i, sampleid) in enumerate(self.info["samples"]):
230                    orng_data.append([self.gdsdata[spot].data[i] for spot in spots] + [sample2class[sampleid]])
231            if missing_class_value == None:
232                orng_data = [example for example in orng_data if example[-1] != None]
233   
234            return orange.ExampleTable(domain, orng_data)
235   
236        else: # genes in rows
237            annotations = self.sample_annotations(sample_type)
238            atts = [orange.FloatVariable(name=ss) for ss in self.info["samples"]]
239            for i, a in enumerate(atts):
240                a.setattr("attributes", annotations[self.info["samples"][i]])
241            domain  = orange.Domain(atts, False)
242   
243            if report_genes: # save by genes
244                domain.addmeta(orange.newmetaid(), orange.StringVariable("gene"))
245                for g in self.gene2spots.keys():
246                    orng_data.append(map(lambda *x: merge_function(x),
247                                         *[self.gdsdata[spot].data for spot in self.gene2spots[g]]))
248            else: # save by spots
249                domain.addmeta(orange.newmetaid(), orange.StringVariable("spot"))
250                spots = self.spot2gene.keys()
251                orng_data = [self.gdsdata[spot].data for spot in spots]
252   
253            data = orange.ExampleTable(domain, orng_data)
254
255            if report_genes:
256                for i, g in enumerate(self.gene2spots.keys()):
257                    data[i]["gene"] = g
258            else:
259                for i, s in enumerate(spots):
260                    data[i]["spot"] = s
261            return data
262       
263    def getdata(self, report_genes=True, merge_function=spots_mean,
264                 sample_type=None, missing_class_value=None,
265                 transpose=False, remove_unknown=None):
266        """Load GDS data and returns a corresponding orange data set,
267        spot<->gene mappings and subset info."""
268        if self.verbose: print "Reading data ..."
269#        if not self.gdsdata:
270        self._parse_soft(remove_unknown = remove_unknown)
271#        if remove_unknown:
272            # some spots were filtered out, need to revise spot<>gene mappings
273        self._getspotmap(include_spots=set(self.gdsdata.keys()))
274        if self.verbose: print "Converting to example table ..."
275        self.data = self._to_ExampleTable(merge_function=merge_function,
276                                          sample_type=sample_type, transpose=transpose,
277                                          report_genes=report_genes)
278        return self.data
279
280    def __str__(self):
281        return "%s (%s), samples=%s, features=%s, genes=%s, subsets=%d" % \
282              (self.info["dataset_id"],
283               self.info["sample_organism"],
284               self.info['sample_count'],
285               self.info['feature_count'],
286               self.info['gene_count'],
287               len(self.info['subsets']),
288               )
289
290
291def _float_or_na(x):
292    if x.isSpecial():
293        return "?"
294    return float(x)
295
296def transpose_class_to_labels(data, attcol="sample"):
297    """Converts data with genes as attributes to data with genes in rows."""
298    if attcol in [v.name for v in data.domain.getmetas().values()]:
299        atts = [orange.FloatVariable(str(d[attcol])) for d in data]
300    else:
301        atts = [orange.FloatVariable("S%d" % i) for i in range(len(data))]
302    for i, d in enumerate(data):
303        atts[i].setattr("class", str(d.getclass()))
304    domain = orange.Domain(atts, False)
305   
306    newdata = []
307    for a in data.domain.attributes:
308        newdata.append([_float_or_na(d[a]) for d in data])
309
310    gene = orange.StringVariable("gene")
311    id = orange.newmetaid()
312    new = orange.ExampleTable(domain, newdata)
313    new.domain.addmeta(id, gene)
314    for i, d in enumerate(new):
315        d[gene] = data.domain.attributes[i].name
316
317    return new
318
319def transpose_labels_to_class(data, class_label=None, gene_label="gene"):
320    """Converts data with genes in rows to data with genes as attributes."""
321    # if no class_label (attribute type) given, guess it from the data
322    if not class_label:
323        l = []
324        for a in data.domain.attributes:
325            l.extend(a.attributes.keys())
326        l = list(set(l))
327        class_label = l[0]
328        if len(set(l)) > 1:
329            import warnings
330            warnings.warn("More than single attribute label types (%s), took %s"
331                          % (", ".join(l), class_label))
332
333    if gene_label in [v.name for v in data.domain.getmetas().values()]:
334        atts = [orange.FloatVariable(str(d[gene_label])) for d in data]
335    else:
336        atts = [orange.FloatVariable("A%d" % i) for i in range(len(data))]
337       
338    classvalues = list(set([a.attributes[class_label] for a in data.domain.attributes]))
339   
340    if all(map(lambda x: isinstance(x, (int, long, float, complex)), classvalues)):
341        classvar = orange.FloatVariable(class_label)
342    else:
343        classvar = orange.EnumVariable(class_label, values=classvalues)
344       
345    domain = orange.Domain(atts + [classvar])
346   
347    newdata = []
348    for a in data.domain.attributes:
349        newdata.append([_float_or_na(d[a]) for d in data] + [a.attributes[class_label]])
350
351    sample = orange.StringVariable("sample")
352    id = orange.newmetaid()
353    new = orange.ExampleTable(domain, newdata)
354    new.domain.addmeta(id, sample)
355    for i, d in enumerate(new):
356        d[sample] = data.domain.attributes[i].name
357
358    return new
359
360def transpose(data):
361    """Transposes data matrix, converts class information to attribute label and back"""
362    if data.domain.classVar:
363        return transpose_class_to_labels(data)
364    else:
365        return transpose_labels_to_class(data)
Note: See TracBrowser for help on using the repository browser.