6
+ − 1 #!/usr/bin/env python
+ − 2 # encoding: utf-8
+ − 3 '''
+ − 4 Module to query a set of accurate mass values detected by high-resolution mass spectrometers
+ − 5 against various repositories/services such as METabolomics EXPlorer database or the
+ − 6 MFSearcher service (http://webs2.kazusa.or.jp/mfsearcher/).
+ − 7
+ − 8 It will take the input file and for each record it will query the
+ − 9 molecular mass in the selected repository/service. If one or more compounds are found
+ − 10 then extra information regarding these compounds is added to the output file.
+ − 11
+ − 12 The output file is thus the input file enriched with information about
+ − 13 related items found in the selected repository/service.
+ − 14
+ − 15 The service should implement the following interface:
+ − 16
+ − 17 http://service_url/mass?targetMs=500&margin=1&marginUnit=ppm&output=txth (txth means there is guaranteed to be a header line before the data)
+ − 18
+ − 19 The output should be tab separated and should contain the following columns (in this order)
+ − 20 db-name molecular-formula dbe formula-weight id description
+ − 21
+ − 22
+ − 23 '''
+ − 24 import csv
+ − 25 import sys
+ − 26 import fileinput
+ − 27 import urllib2
+ − 28 import time
+ − 29 from collections import OrderedDict
+ − 30
+ − 31 __author__ = "Pieter Lukasse"
+ − 32 __contact__ = "pieter.lukasse@wur.nl"
+ − 33 __copyright__ = "Copyright, 2014, Plant Research International, WUR"
+ − 34 __license__ = "Apache v2"
+ − 35
+ − 36 def _process_file(in_xsv, delim='\t'):
+ − 37 '''
+ − 38 Generic method to parse a tab-separated file returning a dictionary with named columns
+ − 39 @param in_csv: input filename to be parsed
+ − 40 '''
+ − 41 data = list(csv.reader(open(in_xsv, 'rU'), delimiter=delim))
+ − 42 return _process_data(data)
+ − 43
+ − 44 def _process_data(data):
+ − 45
+ − 46 header = data.pop(0)
+ − 47 # Create dictionary with column name as key
+ − 48 output = OrderedDict()
+ − 49 for index in xrange(len(header)):
+ − 50 output[header[index]] = [row[index] for row in data]
+ − 51 return output
+ − 52
+ − 53
+ − 54 def _query_and_add_data(input_data, molecular_mass_col, repository_dblink, error_margin, margin_unit):
+ − 55
+ − 56 '''
+ − 57 This method will iterate over the record in the input_data and
+ − 58 will enrich them with the related information found (if any) in the
+ − 59 chosen repository/service
+ − 60
+ − 61 # TODO : could optimize this with multi-threading, see also nice example at http://stackoverflow.com/questions/2846653/python-multithreading-for-dummies
+ − 62 '''
+ − 63 merged = []
+ − 64
+ − 65 for i in xrange(len(input_data[input_data.keys()[0]])):
+ − 66 # Get the record in same dictionary format as input_data, but containing
+ − 67 # a value at each column instead of a list of all values of all records:
+ − 68 input_data_record = OrderedDict(zip(input_data.keys(), [input_data[key][i] for key in input_data.keys()]))
+ − 69
+ − 70 # read the molecular mass :
+ − 71 molecular_mass = input_data_record[molecular_mass_col]
+ − 72
+ − 73
+ − 74 # search for related records in repository/service:
+ − 75 data_found = None
+ − 76 if molecular_mass != "":
+ − 77 molecular_mass = float(molecular_mass)
+ − 78
+ − 79 # 1- search for data around this MM:
+ − 80 query_link = repository_dblink + "/mass?targetMs=" + str(molecular_mass) + "&margin=" + str(error_margin) + "&marginUnit=" + margin_unit + "&output=txth"
+ − 81
+ − 82 data_found = _fire_query_and_return_dict(query_link + "&_format_result=tsv")
+ − 83 data_type_found = "MM"
+ − 84
+ − 85
+ − 86 if data_found == None:
+ − 87 # If still nothing found, just add empty columns
+ − 88 extra_cols = ['', '','','','','']
+ − 89 else:
+ − 90 # Add info found:
+ − 91 extra_cols = _get_extra_info_and_link_cols(data_found, data_type_found, query_link)
+ − 92
+ − 93 # Take all data and merge it into a "flat"/simple array of values:
+ − 94 field_values_list = _merge_data(input_data_record, extra_cols)
+ − 95
+ − 96 merged.append(field_values_list)
+ − 97
+ − 98 # return the merged/enriched records:
+ − 99 return merged
+ − 100
+ − 101
+ − 102 def _get_extra_info_and_link_cols(data_found, data_type_found, query_link):
+ − 103 '''
+ − 104 This method will go over the data found and will return a
+ − 105 list with the following items:
+ − 106 - details of hits found :
+ − 107 db-name molecular-formula dbe formula-weight id description
+ − 108 - Link that executes same query
+ − 109
+ − 110 '''
+ − 111
+ − 112 # set() makes a unique list:
+ − 113 db_name_set = []
+ − 114 molecular_formula_set = []
+ − 115 id_set = []
+ − 116 description_set = []
+ − 117
+ − 118
+ − 119 if 'db-name' in data_found:
+ − 120 db_name_set = set(data_found['db-name'])
+ − 121 elif '# db-name' in data_found:
+ − 122 db_name_set = set(data_found['# db-name'])
+ − 123 if 'molecular-formula' in data_found:
+ − 124 molecular_formula_set = set(data_found['molecular-formula'])
+ − 125 if 'id' in data_found:
+ − 126 id_set = set(data_found['id'])
+ − 127 if 'description' in data_found:
+ − 128 description_set = set(data_found['description'])
+ − 129
+ − 130 result = [data_type_found,
+ − 131 _to_xsv(db_name_set),
+ − 132 _to_xsv(molecular_formula_set),
+ − 133 _to_xsv(id_set),
+ − 134 _to_xsv(description_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 # remove comment lines if any (only leave the one that has "molecular-formula" word in it...compatible with kazusa service):
+ − 161 data_rows_to_remove = []
+ − 162 for data_row in data_rows:
+ − 163 if data_row == "" or (data_row[0] == '#' and "molecular-formula" not in data_row):
+ − 164 data_rows_to_remove.append(data_row)
+ − 165
+ − 166 for data_row in data_rows_to_remove:
+ − 167 data_rows.remove(data_row)
+ − 168
+ − 169 # check if there is any data in the response:
+ − 170 if len(data_rows) <= 1 or data_rows[1].strip() == '':
+ − 171 # means there is only the header row...so no hits:
+ − 172 return None
+ − 173
+ − 174 for data_row in data_rows:
+ − 175 if not data_row.strip() == '':
+ − 176 row_as_list = _str_to_list(data_row, delimiter='\t')
+ − 177 result.append(row_as_list)
+ − 178
+ − 179 # return result processed into a dict:
+ − 180 return _process_data(result)
+ − 181
+ − 182 except urllib2.HTTPError, e:
+ − 183 raise Exception( "HTTP error for URL: " + url + " : %s - " % e.code + e.reason)
+ − 184 except urllib2.URLError, e:
+ − 185 raise Exception( "Network error: %s" % e.reason.args[1] + ". Administrator: please check if service [" + url + "] is accessible from your Galaxy server. ")
+ − 186
+ − 187 def _str_to_list(data_row, delimiter='\t'):
+ − 188 result = []
+ − 189 for column in data_row.split(delimiter):
+ − 190 result.append(column)
+ − 191 return result
+ − 192
+ − 193
+ − 194 # alternative: ?
+ − 195 # s = requests.Session()
+ − 196 # s.verify = False
+ − 197 # #s.auth = (token01, token02)
+ − 198 # resp = s.get(url, params={'name': 'anonymous'}, stream=True)
+ − 199 # content = resp.content
+ − 200 # # transform to dictionary:
+ − 201
+ − 202
+ − 203
+ − 204
+ − 205 def _merge_data(input_data_record, extra_cols):
+ − 206 '''
+ − 207 Adds the extra information to the existing data record and returns
+ − 208 the combined new record.
+ − 209 '''
+ − 210 record = []
+ − 211 for column in input_data_record:
+ − 212 record.append(input_data_record[column])
+ − 213
+ − 214
+ − 215 # add extra columns
+ − 216 for column in extra_cols:
+ − 217 record.append(column)
+ − 218
+ − 219 return record
+ − 220
+ − 221
+ − 222 def _save_data(data_rows, headers, out_csv):
+ − 223 '''
+ − 224 Writes tab-separated data to file
+ − 225 @param data_rows: dictionary containing merged/enriched dataset
+ − 226 @param out_csv: output csv file
+ − 227 '''
+ − 228
+ − 229 # Open output file for writing
+ − 230 outfile_single_handle = open(out_csv, 'wb')
+ − 231 output_single_handle = csv.writer(outfile_single_handle, delimiter="\t")
+ − 232
+ − 233 # Write headers
+ − 234 output_single_handle.writerow(headers)
+ − 235
+ − 236 # Write one line for each row
+ − 237 for data_row in data_rows:
+ − 238 output_single_handle.writerow(data_row)
+ − 239
+ − 240 def _get_repository_URL(repository_file):
+ − 241 '''
+ − 242 Read out and return the URL stored in the given file.
+ − 243 '''
+ − 244 file_input = fileinput.input(repository_file)
+ − 245 try:
+ − 246 for line in file_input:
+ − 247 if line[0] != '#':
+ − 248 # just return the first line that is not a comment line:
+ − 249 return line
+ − 250 finally:
+ − 251 file_input.close()
+ − 252
+ − 253
+ − 254 def main():
+ − 255 '''
+ − 256 Query main function
+ − 257
+ − 258 The input file can be any tabular file, as long as it contains a column for the molecular mass.
+ − 259 This column is then used to query against the chosen repository/service Database.
+ − 260 '''
+ − 261 seconds_start = int(round(time.time()))
+ − 262
+ − 263 input_file = sys.argv[1]
+ − 264 molecular_mass_col = sys.argv[2]
+ − 265 repository_file = sys.argv[3]
+ − 266 error_margin = float(sys.argv[4])
+ − 267 margin_unit = sys.argv[5]
+ − 268 output_result = sys.argv[6]
+ − 269
+ − 270 # Parse repository_file to find the URL to the service:
+ − 271 repository_dblink = _get_repository_URL(repository_file)
+ − 272
+ − 273 # Parse tabular input file into dictionary/array:
+ − 274 input_data = _process_file(input_file)
+ − 275
+ − 276 # Query data against repository :
+ − 277 enriched_data = _query_and_add_data(input_data, molecular_mass_col, repository_dblink, error_margin, margin_unit)
+ − 278 headers = input_data.keys() + ['SEARCH hits for ','SEARCH hits: db-names', 'SEARCH hits: molecular-formulas ',
+ − 279 'SEARCH hits: ids','SEARCH hits: descriptions', 'Link to SEARCH hits'] #TODO - add min and max formula weigth columns
+ − 280
+ − 281 _save_data(enriched_data, headers, output_result)
+ − 282
+ − 283 seconds_end = int(round(time.time()))
+ − 284 print "Took " + str(seconds_end - seconds_start) + " seconds"
+ − 285
+ − 286
+ − 287
+ − 288 if __name__ == '__main__':
+ − 289 main()
