source: orange-bioinformatics/Orange/bioinformatics/obiGEO.py @ 1625:cefeb35cbfc9

Revision 1625:cefeb35cbfc9, 14.3 KB checked in by mitar, 2 years ago (diff)

Moving files around.

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