6
+ − 1 #!/usr/bin/env python
+ − 2 # encoding: utf-8
+ − 3 '''
+ − 4 Module to query a set of identifications against the METabolomics EXPlorer database.
+ − 5
+ − 6 It will take the input file and for each record it will query the
+ − 7 molecular mass in the selected MetExp DB. If one or more compounds are found in the
+ − 8 MetExp DB then extra information regarding these compounds is added to the output file.
+ − 9
+ − 10 The output file is thus the input file enriched with information about
+ − 11 related items found in the selected MetExp DB.
+ − 12 '''
+ − 13 import csv
+ − 14 import sys
+ − 15 import fileinput
+ − 16 import urllib2
+ − 17 import time
+ − 18 from collections import OrderedDict
+ − 19
+ − 20 __author__ = "Pieter Lukasse"
+ − 21 __contact__ = "pieter.lukasse@wur.nl"
+ − 22 __copyright__ = "Copyright, 2014, Plant Research International, WUR"
+ − 23 __license__ = "Apache v2"
+ − 24
+ − 25 def _process_file(in_xsv, delim='\t'):
+ − 26 '''
+ − 27 Generic method to parse a tab-separated file returning a dictionary with named columns
+ − 28 @param in_csv: input filename to be parsed
+ − 29 '''
+ − 30 data = list(csv.reader(open(in_xsv, 'rU'), delimiter=delim))
+ − 31 return _process_data(data)
+ − 32
+ − 33 def _process_data(data):
+ − 34
+ − 35 header = data.pop(0)
+ − 36 # Create dictionary with column name as key
+ − 37 output = OrderedDict()
+ − 38 for index in xrange(len(header)):
+ − 39 output[header[index]] = [row[index] for row in data]
+ − 40 return output
+ − 41
+ − 42
+ − 43 def _query_and_add_data(input_data, casid_col, formula_col, molecular_mass_col, metexp_dblink, separation_method):
+ − 44 '''
+ − 45 This method will iterate over the record in the input_data and
+ − 46 will enrich them with the related information found (if any) in the
+ − 47 MetExp Database.
+ − 48
+ − 49 # TODO : could optimize this with multi-threading, see also nice example at http://stackoverflow.com/questions/2846653/python-multithreading-for-dummies
+ − 50 '''
+ − 51 merged = []
+ − 52
+ − 53 for i in xrange(len(input_data[input_data.keys()[0]])):
+ − 54 # Get the record in same dictionary format as input_data, but containing
+ − 55 # a value at each column instead of a list of all values of all records:
+ − 56 input_data_record = OrderedDict(zip(input_data.keys(), [input_data[key][i] for key in input_data.keys()]))
+ − 57
+ − 58 # read the molecular mass and formula:
+ − 59 cas_id = input_data_record[casid_col]
+ − 60 formula = input_data_record[formula_col]
+ − 61 molecular_mass = input_data_record[molecular_mass_col]
+ − 62
+ − 63 # search for related records in MetExp:
+ − 64 data_found = None
+ − 65 if cas_id != "undef":
+ − 66 # 1- search for other experiments where this CAS id has been found:
+ − 67 query_link = metexp_dblink + "/find_entries/query?cas_nr="+ cas_id + "&method=" + separation_method
+ − 68 data_found = _fire_query_and_return_dict(query_link + "&_format_result=tsv")
+ − 69 data_type_found = "CAS"
+ − 70 if data_found == None:
+ − 71 # 2- search for other experiments where this FORMULA has been found:
+ − 72 query_link = metexp_dblink + "/find_entries/query?molecule_formula="+ formula + "&method=" + separation_method
+ − 73 data_found = _fire_query_and_return_dict(query_link + "&_format_result=tsv")
+ − 74 data_type_found = "FORMULA"
+ − 75 if data_found == None:
+ − 76 # 3- search for other experiments where this MM has been found:
+ − 77 query_link = metexp_dblink + "/find_entries/query?molecule_mass="+ molecular_mass + "&method=" + separation_method
+ − 78 data_found = _fire_query_and_return_dict(query_link + "&_format_result=tsv")
+ − 79 data_type_found = "MM"
+ − 80
+ − 81 if data_found == None:
+ − 82 # If still nothing found, just add empty columns
+ − 83 extra_cols = ['', '','','','','','','']
+ − 84 else:
+ − 85 # Add info found:
+ − 86 extra_cols = _get_extra_info_and_link_cols(data_found, data_type_found, query_link)
+ − 87
+ − 88 # Take all data and merge it into a "flat"/simple array of values:
+ − 89 field_values_list = _merge_data(input_data_record, extra_cols)
+ − 90
+ − 91 merged.append(field_values_list)
+ − 92
+ − 93 # return the merged/enriched records:
+ − 94 return merged
+ − 95
+ − 96
+ − 97 def _get_extra_info_and_link_cols(data_found, data_type_found, query_link):
+ − 98 '''
+ − 99 This method will go over the data found and will return a
+ − 100 list with the following items:
+ − 101 - Experiment details where hits have been found :
+ − 102 'organism', 'tissue','experiment_name','user_name','column_type'
+ − 103 - Link that executes same query
+ − 104
+ − 105 '''
+ − 106 # set() makes a unique list:
+ − 107 organism_set = []
+ − 108 tissue_set = []
+ − 109 experiment_name_set = []
+ − 110 user_name_set = []
+ − 111 column_type_set = []
+ − 112 cas_nr_set = []
+ − 113
+ − 114 if 'organism' in data_found:
+ − 115 organism_set = set(data_found['organism'])
+ − 116 if 'tissue' in data_found:
+ − 117 tissue_set = set(data_found['tissue'])
+ − 118 if 'experiment_name' in data_found:
+ − 119 experiment_name_set = set(data_found['experiment_name'])
+ − 120 if 'user_name' in data_found:
+ − 121 user_name_set = set(data_found['user_name'])
+ − 122 if 'column_type' in data_found:
+ − 123 column_type_set = set(data_found['column_type'])
+ − 124 if 'CAS' in data_found:
+ − 125 cas_nr_set = set(data_found['CAS'])
+ − 126
+ − 127
+ − 128 result = [data_type_found,
+ − 129 _to_xsv(organism_set),
+ − 130 _to_xsv(tissue_set),
+ − 131 _to_xsv(experiment_name_set),
+ − 132 _to_xsv(user_name_set),
+ − 133 _to_xsv(column_type_set),
+ − 134 _to_xsv(cas_nr_set),
+ − 135 #To let Excel interpret as link, use e.g. =HYPERLINK("http://stackoverflow.com", "friendly name"):
+ − 136 "=HYPERLINK(\""+ query_link + "\", \"Link to entries found in DB \")"]
+ − 137 return result
+ − 138
+ − 139
+ − 140 def _to_xsv(data_set):
+ − 141 result = ""
+ − 142 for item in data_set:
+ − 143 result = result + str(item) + "|"
+ − 144 return result
+ − 145
+ − 146
+ − 147 def _fire_query_and_return_dict(url):
+ − 148 '''
+ − 149 This method will fire the query as a web-service call and
+ − 150 return the results as a list of dictionary objects
+ − 151 '''
+ − 152
+ − 153 try:
+ − 154 data = urllib2.urlopen(url).read()
+ − 155
+ − 156 # transform to dictionary:
+ − 157 result = []
+ − 158 data_rows = data.split("\n")
+ − 159
+ − 160 # check if there is any data in the response:
+ − 161 if len(data_rows) <= 1 or data_rows[1].strip() == '':
+ − 162 # means there is only the header row...so no hits:
+ − 163 return None
+ − 164
+ − 165 for data_row in data_rows:
+ − 166 if not data_row.strip() == '':
+ − 167 row_as_list = _str_to_list(data_row, delimiter='\t')
+ − 168 result.append(row_as_list)
+ − 169
+ − 170 # return result processed into a dict:
+ − 171 return _process_data(result)
+ − 172
+ − 173 except urllib2.HTTPError, e:
+ − 174 raise Exception( "HTTP error for URL: " + url + " : %s - " % e.code + e.reason)
+ − 175 except urllib2.URLError, e:
+ − 176 raise Exception( "Network error: %s" % e.reason.args[1] + ". Administrator: please check if MetExp service [" + url + "] is accessible from your Galaxy server. ")
+ − 177
+ − 178 def _str_to_list(data_row, delimiter='\t'):
+ − 179 result = []
+ − 180 for column in data_row.split(delimiter):
+ − 181 result.append(column)
+ − 182 return result
+ − 183
+ − 184
+ − 185 # alternative: ?
+ − 186 # s = requests.Session()
+ − 187 # s.verify = False
+ − 188 # #s.auth = (token01, token02)
+ − 189 # resp = s.get(url, params={'name': 'anonymous'}, stream=True)
+ − 190 # content = resp.content
+ − 191 # # transform to dictionary:
+ − 192
+ − 193
+ − 194
+ − 195
+ − 196 def _merge_data(input_data_record, extra_cols):
+ − 197 '''
+ − 198 Adds the extra information to the existing data record and returns
+ − 199 the combined new record.
+ − 200 '''
+ − 201 record = []
+ − 202 for column in input_data_record:
+ − 203 record.append(input_data_record[column])
+ − 204
+ − 205
+ − 206 # add extra columns
+ − 207 for column in extra_cols:
+ − 208 record.append(column)
+ − 209
+ − 210 return record
+ − 211
+ − 212
+ − 213 def _save_data(data_rows, headers, out_csv):
+ − 214 '''
+ − 215 Writes tab-separated data to file
+ − 216 @param data_rows: dictionary containing merged/enriched dataset
+ − 217 @param out_csv: output csv file
+ − 218 '''
+ − 219
+ − 220 # Open output file for writing
+ − 221 outfile_single_handle = open(out_csv, 'wb')
+ − 222 output_single_handle = csv.writer(outfile_single_handle, delimiter="\t")
+ − 223
+ − 224 # Write headers
+ − 225 output_single_handle.writerow(headers)
+ − 226
+ − 227 # Write one line for each row
+ − 228 for data_row in data_rows:
+ − 229 output_single_handle.writerow(data_row)
+ − 230
+ − 231 def _get_metexp_URL(metexp_dblink_file):
+ − 232 '''
+ − 233 Read out and return the URL stored in the given file.
+ − 234 '''
+ − 235 file_input = fileinput.input(metexp_dblink_file)
+ − 236 try:
+ − 237 for line in file_input:
+ − 238 if line[0] != '#':
+ − 239 # just return the first line that is not a comment line:
+ − 240 return line
+ − 241 finally:
+ − 242 file_input.close()
+ − 243
+ − 244
+ − 245 def main():
+ − 246 '''
+ − 247 MetExp Query main function
+ − 248
+ − 249 The input file can be any tabular file, as long as it contains a column for the molecular mass
+ − 250 and one for the formula of the respective identification. These two columns are then
+ − 251 used to query against MetExp Database.
+ − 252 '''
+ − 253 seconds_start = int(round(time.time()))
+ − 254
+ − 255 input_file = sys.argv[1]
+ − 256 casid_col = sys.argv[2]
+ − 257 formula_col = sys.argv[3]
+ − 258 molecular_mass_col = sys.argv[4]
+ − 259 metexp_dblink_file = sys.argv[5]
+ − 260 separation_method = sys.argv[6]
+ − 261 output_result = sys.argv[7]
+ − 262
+ − 263 # Parse metexp_dblink_file to find the URL to the MetExp service:
+ − 264 metexp_dblink = _get_metexp_URL(metexp_dblink_file)
+ − 265
+ − 266 # Parse tabular input file into dictionary/array:
+ − 267 input_data = _process_file(input_file)
+ − 268
+ − 269 # Query data against MetExp DB :
+ − 270 enriched_data = _query_and_add_data(input_data, casid_col, formula_col, molecular_mass_col, metexp_dblink, separation_method)
+ − 271 headers = input_data.keys() + ['METEXP hits for ','METEXP hits: organisms', 'METEXP hits: tissues',
+ − 272 'METEXP hits: experiments','METEXP hits: user names','METEXP hits: column types', 'METEXP hits: CAS nrs', 'Link to METEXP hits']
+ − 273
+ − 274 _save_data(enriched_data, headers, output_result)
+ − 275
+ − 276 seconds_end = int(round(time.time()))
+ − 277 print "Took " + str(seconds_end - seconds_start) + " seconds"
+ − 278
+ − 279
+ − 280
+ − 281 if __name__ == '__main__':
+ − 282 main()
