|
4
|
1 #######################################################################################
|
|
|
2 # Python-code: SAINT pre-processing from Scaffold "Samples Report" output
|
|
|
3 # Author: Brent Kuenzi
|
|
|
4 #######################################################################################
|
|
|
5 # This program reads in a raw Scaffold "Samples Report" output and a user generated
|
|
|
6 # bait file and autoformats it into prey and interaction files for SAINTexpress
|
|
|
7 # analysis
|
|
|
8 #######################################################################################
|
|
|
9 # Copyright (C) Brent Kuenzi.
|
|
|
10 # Permission is granted to copy, distribute and/or modify this document
|
|
|
11 # under the terms of the GNU Free Documentation License, Version 1.3
|
|
|
12 # or any later version published by the Free Software Foundation;
|
|
|
13 # with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
|
|
|
14 # A copy of the license is included in the section entitled "GNU
|
|
|
15 # Free Documentation License".
|
|
|
16 #######################################################################################
|
|
|
17 ## REQUIRED INPUT ##
|
|
|
18
|
|
|
19 # 1) infile: Scaffold "Samples Report" output
|
|
|
20 # 2) baitfile: SAINT formatted bait file generated in Galaxy
|
|
|
21 # 3) fasta_db: fasta database for use (defaults to SwissProt_HUMAN_2014_08.fasta)
|
|
|
22 # 4) prey: Y or N for generating a prey file
|
|
|
23 # 5) make_bait: String of bait names, assignment, and test or control boolean
|
|
|
24 #######################################################################################
|
|
|
25
|
|
|
26 import sys
|
|
|
27 import os.path
|
|
11
|
28 import re
|
|
4
|
29
|
|
|
30
|
|
|
31 infile = sys.argv[1]
|
|
|
32 #Scaffold "Samples Report" output.
|
|
|
33 prey = sys.argv[2]
|
|
|
34 # Y or N boolean from Galaxy.
|
|
|
35 fasta_db = sys.argv[3]
|
|
|
36 tool_path = sys.argv[8]
|
|
|
37 if fasta_db == "None":
|
|
|
38 fasta_db = str(tool_path) + "/SwissProt_HUMAN_2014_08.fasta"
|
|
|
39 make_bait = sys.argv[6]
|
|
|
40 bait_bool = sys.argv[9]
|
|
|
41
|
|
|
42
|
|
|
43 def bait_create(baits, infile):
|
|
|
44 # Verifies the Baits are valid in the Scaffold file and writes the Bait.txt.
|
|
|
45 baits = make_bait.split()
|
|
|
46 i = 0
|
|
|
47 bait_file_tmp = open("bait.txt", "w")
|
|
|
48 order = []
|
|
|
49 bait_cache = []
|
|
|
50 while i < len(baits):
|
|
|
51 if baits[i+2] == "true":
|
|
|
52 T_C = "C"
|
|
|
53 else:
|
|
|
54 T_C = "T"
|
|
|
55 bait_line = baits[i] + "\t" + baits[i+1] + "\t" + T_C + "\n"
|
|
16
|
56 bait_cache.append(str(bait_line))
|
|
4
|
57 i = i + 3
|
|
|
58
|
|
|
59 for cache_line in bait_cache:
|
|
17
|
60 bait_file_tmp.write(cache_line)
|
|
4
|
61
|
|
|
62 bait_file_tmp.close()
|
|
|
63
|
|
|
64 if bait_bool == 'false':
|
|
|
65 bait_create(make_bait, infile)
|
|
|
66 baitfile = "bait.txt"
|
|
|
67 else:
|
|
|
68 bait_temp_file = open(sys.argv[10], 'r')
|
|
|
69 bait_cache = bait_temp_file.readlines()
|
|
|
70 bait_file_tmp = open("bait.txt", "wr")
|
|
|
71 for cache_line in bait_cache:
|
|
|
72 bait_file_tmp.write(cache_line)
|
|
|
73 bait_file_tmp.close()
|
|
|
74 baitfile = "bait.txt"
|
|
|
75
|
|
|
76
|
|
|
77 class ReturnValue1(object):
|
|
|
78 def __init__(self, sequence, gene):
|
|
|
79 self.seqlength = sequence
|
|
|
80 self.genename = gene
|
|
|
81
|
|
|
82
|
|
|
83 class ReturnValue2(object):
|
|
|
84 def __init__(self, getdata, getproteins, getheader):
|
|
|
85 self.data = getdata
|
|
|
86 self.proteins = getproteins
|
|
|
87 self.header = getheader
|
|
|
88
|
|
|
89
|
|
|
90 def main(Scaffold_input, baits):
|
|
|
91 bait_check(baitfile, Scaffold_input)
|
|
|
92 make_inter(Scaffold_input)
|
|
|
93 if prey == 'true':
|
|
|
94 make_prey(Scaffold_input)
|
|
|
95 no_error_inter(Scaffold_input)
|
|
|
96 os.rename('prey.txt', sys.argv[5])
|
|
|
97 elif prey == 'false':
|
|
|
98 if os.path.isfile('error proteins.txt') == True:
|
|
|
99 no_error_inter(Scaffold_input)
|
|
|
100 pass
|
|
|
101 elif prey != 'true' or 'false':
|
|
|
102 sys.exit("Invalid Prey Argument: Y or N")
|
|
|
103
|
|
|
104
|
|
|
105 def get_info(uniprot_accession_in):
|
|
|
106 # Get aminoacid lengths and gene name.
|
|
|
107 error = open('error proteins.txt', 'a+')
|
|
|
108 data = open(fasta_db, 'r')
|
|
|
109 data_lines = data.readlines()
|
|
|
110 db_len = len(data_lines)
|
|
|
111 seqlength = 0
|
|
|
112 count = 0
|
|
53
|
113 last_line = data_lines[-1]
|
|
4
|
114 for data_line in data_lines:
|
|
|
115 if ">sp" in data_line:
|
|
|
116 namer = data_line.split("|")[2]
|
|
|
117 if uniprot_accession_in == data_line.split("|")[1]:
|
|
|
118 match = count + 1
|
|
|
119 if 'GN=' in data_line:
|
|
10
|
120 lst = data_line.split('GN=')
|
|
4
|
121 lst2 = lst[1].split(' ')
|
|
|
122 genename = lst2[0]
|
|
|
123 if 'GN=' not in data_line:
|
|
|
124 genename = 'NA'
|
|
|
125 while ">sp" not in data_lines[match]:
|
|
|
126 if match <= db_len:
|
|
|
127 seqlength = seqlength + len(data_lines[match].strip())
|
|
53
|
128 if data_lines[match] == last_line:
|
|
|
129 break
|
|
4
|
130 match = match + 1
|
|
|
131 else:
|
|
|
132 break
|
|
|
133 return ReturnValue1(seqlength, genename)
|
|
10
|
134 if uniprot_accession_in == namer.split(" ")[0]:
|
|
4
|
135 match = count + 1
|
|
|
136 # Ensures consistent spacing throughout.
|
|
|
137 if 'GN=' in data_line:
|
|
|
138 lst = data_line.split('GN=')
|
|
|
139 lst2 = lst[1].split(' ')
|
|
|
140 genename = lst2[0]
|
|
|
141 if 'GN=' not in data_line:
|
|
|
142 genename = 'NA'
|
|
|
143 while ">sp" not in data_lines[match]:
|
|
|
144 if match <= db_len:
|
|
|
145 seqlength = seqlength + len(data_lines[match].strip())
|
|
53
|
146 if data_lines[match] == last_line:
|
|
|
147 break
|
|
4
|
148 match = match + 1
|
|
|
149 else:
|
|
|
150 break
|
|
|
151 return ReturnValue1(seqlength, genename)
|
|
|
152 count = count + 1
|
|
|
153 if seqlength == 0:
|
|
|
154 error.write(uniprot_accession_in + '\t' + "Uniprot not in Fasta" + '\n')
|
|
|
155 error.close
|
|
|
156 seqlength = 'NA'
|
|
|
157 genename = 'NA'
|
|
|
158 return ReturnValue1(seqlength, genename)
|
|
|
159
|
|
|
160
|
|
|
161 def readtab(infile):
|
|
|
162 with open(infile, 'r') as input_file:
|
|
|
163 # read in tab-delim text
|
|
|
164 output = []
|
|
|
165 for input_line in input_file:
|
|
|
166 input_line = input_line.strip()
|
|
|
167 temp = input_line.split('\t')
|
|
|
168 output.append(temp)
|
|
|
169 return output
|
|
|
170
|
|
|
171
|
|
|
172 def read_Scaffold(Scaffold_input):
|
|
|
173 # Get data, proteins and header from Scaffold output
|
|
|
174 dupes = readtab(Scaffold_input)
|
|
|
175 cnt = 0
|
|
|
176 for Scaffold_line in dupes:
|
|
|
177 cnt += 1
|
|
|
178 if Scaffold_line[0] == '#':
|
|
|
179 # Finds the start of second header.
|
|
|
180 header_start = cnt-1
|
|
|
181 header = dupes[header_start]
|
|
|
182 prot_start = header.index("Accession Number")
|
|
|
183 data = dupes[header_start+1:len(dupes)-2]
|
|
|
184 # Cut off blank line and END OF FILE.
|
|
|
185 proteins = []
|
|
|
186 for Scaffold_line in data:
|
|
|
187 Scaffold_line[4] = Scaffold_line[4].split()[0]
|
|
|
188 # Removes the (+##) that sometimes is attached.
|
|
11
|
189 uniprot_re = re.compile("[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}")
|
|
12
|
190 for protein in data:
|
|
11
|
191 prot_id = uniprot_re.match(protein[prot_start])
|
|
12
|
192 if prot_id:
|
|
|
193 proteins.append(prot_id.group())
|
|
|
194 else:
|
|
|
195 prot_ids = protein[prot_start].split("|")
|
|
|
196 for prot_id in prot_ids:
|
|
|
197 if "_HUMAN" in prot_id:
|
|
|
198 proteins.append(prot_id)
|
|
|
199 elif "_YEAST" in prot_id:
|
|
|
200 proteins.append(prot_id)
|
|
|
201 elif "_MOUSE" in prot_id:
|
|
|
202 proteins.append(prot_id)
|
|
13
|
203 else:
|
|
14
|
204 print "Accession must be uniprot ID or gene"
|
|
4
|
205 return ReturnValue2(data, proteins, header)
|
|
|
206
|
|
|
207
|
|
|
208 def make_inter(Scaffold_input):
|
|
|
209 bait = readtab(baitfile)
|
|
|
210 data = read_Scaffold(Scaffold_input).data
|
|
|
211 header = read_Scaffold(Scaffold_input).header
|
|
|
212 proteins = read_Scaffold(Scaffold_input).proteins
|
|
|
213 bait_index = []
|
|
|
214 for bait_line in bait:
|
|
|
215 bait_index.append(header.index(bait_line[0]))
|
|
|
216 # Find just the baits defined in bait file.
|
|
|
217 with open('inter.txt', 'w') as inter_file:
|
|
|
218 a = 0; l = 0
|
|
|
219 for bb in bait:
|
|
|
220 for lst in data:
|
|
|
221 inter_file.write(header[bait_index[l]] + '\t' + bb[1] + '\t' + proteins[a] + '\t'
|
|
|
222 + lst[bait_index[l]] + '\n')
|
|
|
223 a += 1
|
|
|
224 if a == len(proteins):
|
|
|
225 a = 0; l += 1
|
|
|
226
|
|
|
227
|
|
|
228 def make_prey(Scaffold_input):
|
|
|
229 proteins = read_Scaffold(Scaffold_input).proteins
|
|
|
230 output_file = open("prey.txt", 'w')
|
|
|
231 for protein in proteins:
|
|
|
232 protein = protein.replace("\n", "")
|
|
|
233 # Remove \n for input into function.
|
|
|
234 protein = protein.replace("\r", "")
|
|
|
235 # Ditto for \r.
|
|
|
236 seq = get_info(protein).seqlength
|
|
|
237 GN = get_info(protein).genename
|
|
|
238 if seq != 'NA':
|
|
33
|
239 if GN != 'NA':
|
|
|
240 output_file.write(protein + "\t" + str(seq) + "\t" + str(GN) + "\n")
|
|
4
|
241 output_file.close()
|
|
|
242
|
|
|
243
|
|
|
244 def no_error_inter(Scaffold_input):
|
|
|
245 # Remake inter file without protein errors from Uniprot.
|
|
49
|
246 err = readtab("error proteins.txt")
|
|
4
|
247 bait = readtab(baitfile)
|
|
|
248 data = read_Scaffold(Scaffold_input).data
|
|
|
249 header = read_Scaffold(Scaffold_input).header
|
|
|
250 bait_index = []
|
|
51
|
251 for bait_line in bait:
|
|
|
252 bait_index.append(header.index(bait_line[0]))
|
|
4
|
253 proteins = read_Scaffold(Scaffold_input).proteins
|
|
|
254 errors = []
|
|
|
255 for e in err:
|
|
|
256 errors.append(e[0])
|
|
51
|
257 with open('inter.txt', 'w') as y:
|
|
4
|
258 l = 0; a = 0
|
|
|
259 for bb in bait:
|
|
|
260 for lst in data:
|
|
51
|
261 if proteins[a] not in errors:
|
|
|
262 y.write(header[bait_index[l]] + '\t' + bb[1] + '\t' + proteins[a] + '\t'
|
|
|
263 + lst[bait_index[l]] + '\n')
|
|
4
|
264 a += 1
|
|
|
265 if a == len(proteins):
|
|
|
266 l += 1; a = 0
|
|
|
267
|
|
|
268
|
|
|
269 def bait_check(bait, Scaffold_input):
|
|
|
270 # Check that bait names share Scaffold header titles.
|
|
|
271 bait_in = readtab(bait)
|
|
|
272 header = read_Scaffold(Scaffold_input).header
|
|
|
273 for i in bait_in:
|
|
|
274 if i[0] not in header:
|
|
|
275 sys.exit("Bait must share header titles with Scaffold output")
|
|
|
276
|
|
|
277 if __name__ == '__main__':
|
|
|
278 main(infile, baitfile)
|
|
|
279
|
|
|
280 os.rename("inter.txt", sys.argv[4])
|
|
|
281 os.rename("bait.txt", sys.argv[7])
|
