annotate mzID_process2.py @ 65:a551998e1068 draft

Uploaded
author bornea
date Sat, 27 Aug 2016 22:13:31 -0400
parents
children b8169cc0b860
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
65
a551998e1068 Uploaded
bornea
parents:
diff changeset
1 # -*- coding: utf-8 -*-
a551998e1068 Uploaded
bornea
parents:
diff changeset
2 """
a551998e1068 Uploaded
bornea
parents:
diff changeset
3 Python-code: Preprocess mzIdentML
a551998e1068 Uploaded
bornea
parents:
diff changeset
4 @author = Brent Kuenzi
a551998e1068 Uploaded
bornea
parents:
diff changeset
5 @email = Brent.Kuenzi@moffitt.org
a551998e1068 Uploaded
bornea
parents:
diff changeset
6 """
a551998e1068 Uploaded
bornea
parents:
diff changeset
7 #######################################################################################
a551998e1068 Uploaded
bornea
parents:
diff changeset
8 ## Description: ##
a551998e1068 Uploaded
bornea
parents:
diff changeset
9 # This program will create inter, prey, and bait files from mzIdentML files
a551998e1068 Uploaded
bornea
parents:
diff changeset
10 ## Required input: ##
a551998e1068 Uploaded
bornea
parents:
diff changeset
11 # 1) mzIdentML file to be reformatted
a551998e1068 Uploaded
bornea
parents:
diff changeset
12 # 2) minimum PSM for quantification
a551998e1068 Uploaded
bornea
parents:
diff changeset
13
a551998e1068 Uploaded
bornea
parents:
diff changeset
14
a551998e1068 Uploaded
bornea
parents:
diff changeset
15 import sys
a551998e1068 Uploaded
bornea
parents:
diff changeset
16 import os
a551998e1068 Uploaded
bornea
parents:
diff changeset
17
a551998e1068 Uploaded
bornea
parents:
diff changeset
18 ins_path = sys.argv[5]
a551998e1068 Uploaded
bornea
parents:
diff changeset
19
a551998e1068 Uploaded
bornea
parents:
diff changeset
20 class ReturnValue1(object):
a551998e1068 Uploaded
bornea
parents:
diff changeset
21 def __init__(self, sequence, gene):
a551998e1068 Uploaded
bornea
parents:
diff changeset
22 self.seqlength = sequence
a551998e1068 Uploaded
bornea
parents:
diff changeset
23 self.genename = gene
a551998e1068 Uploaded
bornea
parents:
diff changeset
24 class ReturnValue2(object):
a551998e1068 Uploaded
bornea
parents:
diff changeset
25 def __init__(self, inter, accessions):
a551998e1068 Uploaded
bornea
parents:
diff changeset
26 self.inter = inter
a551998e1068 Uploaded
bornea
parents:
diff changeset
27 self.accessions = accessions
a551998e1068 Uploaded
bornea
parents:
diff changeset
28 def read_tab(infile):
a551998e1068 Uploaded
bornea
parents:
diff changeset
29 with open(infile,'r') as x:
a551998e1068 Uploaded
bornea
parents:
diff changeset
30 output = []
a551998e1068 Uploaded
bornea
parents:
diff changeset
31 for line in x:
a551998e1068 Uploaded
bornea
parents:
diff changeset
32 line = line.strip()
a551998e1068 Uploaded
bornea
parents:
diff changeset
33 temp = line.split('\t')
a551998e1068 Uploaded
bornea
parents:
diff changeset
34 output.append(temp)
a551998e1068 Uploaded
bornea
parents:
diff changeset
35 return output
a551998e1068 Uploaded
bornea
parents:
diff changeset
36 def printProgress (iteration, total, prefix = '', suffix = '', decimals = 1, barLength = 100):
a551998e1068 Uploaded
bornea
parents:
diff changeset
37 """
a551998e1068 Uploaded
bornea
parents:
diff changeset
38 Call in a loop to create terminal progress bar
a551998e1068 Uploaded
bornea
parents:
diff changeset
39 @params:
a551998e1068 Uploaded
bornea
parents:
diff changeset
40 iteration - Required : current iteration (Int)
a551998e1068 Uploaded
bornea
parents:
diff changeset
41 total - Required : total iterations (Int)
a551998e1068 Uploaded
bornea
parents:
diff changeset
42 prefix - Optional : prefix string (Str)
a551998e1068 Uploaded
bornea
parents:
diff changeset
43 suffix - Optional : suffix string (Str)
a551998e1068 Uploaded
bornea
parents:
diff changeset
44 decimals - Optional : positive number of decimals in percent complete (Int)
a551998e1068 Uploaded
bornea
parents:
diff changeset
45 barLength - Optional : character length of bar (Int)
a551998e1068 Uploaded
bornea
parents:
diff changeset
46 """
a551998e1068 Uploaded
bornea
parents:
diff changeset
47 formatStr = "{0:." + str(decimals) + "f}"
a551998e1068 Uploaded
bornea
parents:
diff changeset
48 percents = formatStr.format(100 * (iteration / float(total)))
a551998e1068 Uploaded
bornea
parents:
diff changeset
49 filledLength = int(round(barLength * iteration / float(total)))
a551998e1068 Uploaded
bornea
parents:
diff changeset
50 bar = '=' * filledLength + '-' * (barLength - filledLength)
a551998e1068 Uploaded
bornea
parents:
diff changeset
51 sys.stdout.write('\r%s |%s| %s%s %s' % (prefix, bar, percents, '%', suffix)),
a551998e1068 Uploaded
bornea
parents:
diff changeset
52 sys.stdout.flush()
a551998e1068 Uploaded
bornea
parents:
diff changeset
53 if iteration == total:
a551998e1068 Uploaded
bornea
parents:
diff changeset
54 sys.stdout.write('\n')
a551998e1068 Uploaded
bornea
parents:
diff changeset
55 sys.stdout.flush()
a551998e1068 Uploaded
bornea
parents:
diff changeset
56 def get_info(uniprot_accession_in,fasta_db):
a551998e1068 Uploaded
bornea
parents:
diff changeset
57 # Get aminoacid lengths and gene name.
a551998e1068 Uploaded
bornea
parents:
diff changeset
58 error = open('error proteins.txt', 'a+')
a551998e1068 Uploaded
bornea
parents:
diff changeset
59 data = open(fasta_db, 'r')
a551998e1068 Uploaded
bornea
parents:
diff changeset
60 data_lines = data.readlines()
a551998e1068 Uploaded
bornea
parents:
diff changeset
61 db_len = len(data_lines)
a551998e1068 Uploaded
bornea
parents:
diff changeset
62 seqlength = 0
a551998e1068 Uploaded
bornea
parents:
diff changeset
63 count = 0
a551998e1068 Uploaded
bornea
parents:
diff changeset
64 last_line = data_lines[-1]
a551998e1068 Uploaded
bornea
parents:
diff changeset
65 for data_line in data_lines:
a551998e1068 Uploaded
bornea
parents:
diff changeset
66 if ">sp" in data_line:
a551998e1068 Uploaded
bornea
parents:
diff changeset
67 namer = data_line.split("|")[2]
a551998e1068 Uploaded
bornea
parents:
diff changeset
68 if uniprot_accession_in == data_line.split("|")[1]:
a551998e1068 Uploaded
bornea
parents:
diff changeset
69 match = count + 1
a551998e1068 Uploaded
bornea
parents:
diff changeset
70 if 'GN=' in data_line:
a551998e1068 Uploaded
bornea
parents:
diff changeset
71 lst = data_line.split('GN=')
a551998e1068 Uploaded
bornea
parents:
diff changeset
72 lst2 = lst[1].split(' ')
a551998e1068 Uploaded
bornea
parents:
diff changeset
73 genename = lst2[0]
a551998e1068 Uploaded
bornea
parents:
diff changeset
74 if 'GN=' not in data_line:
a551998e1068 Uploaded
bornea
parents:
diff changeset
75 genename = 'NA'
a551998e1068 Uploaded
bornea
parents:
diff changeset
76 while ">sp" not in data_lines[match]:
a551998e1068 Uploaded
bornea
parents:
diff changeset
77 if match <= db_len:
a551998e1068 Uploaded
bornea
parents:
diff changeset
78 seqlength = seqlength + len(data_lines[match].strip())
a551998e1068 Uploaded
bornea
parents:
diff changeset
79 if data_lines[match] == last_line:
a551998e1068 Uploaded
bornea
parents:
diff changeset
80 break
a551998e1068 Uploaded
bornea
parents:
diff changeset
81 match = match + 1
a551998e1068 Uploaded
bornea
parents:
diff changeset
82 else:
a551998e1068 Uploaded
bornea
parents:
diff changeset
83 break
a551998e1068 Uploaded
bornea
parents:
diff changeset
84 return ReturnValue1(seqlength, genename)
a551998e1068 Uploaded
bornea
parents:
diff changeset
85 if uniprot_accession_in == namer.split(" ")[0]:
a551998e1068 Uploaded
bornea
parents:
diff changeset
86 match = count + 1
a551998e1068 Uploaded
bornea
parents:
diff changeset
87 # Ensures consistent spacing throughout.
a551998e1068 Uploaded
bornea
parents:
diff changeset
88 if 'GN=' in data_line:
a551998e1068 Uploaded
bornea
parents:
diff changeset
89 lst = data_line.split('GN=')
a551998e1068 Uploaded
bornea
parents:
diff changeset
90 lst2 = lst[1].split(' ')
a551998e1068 Uploaded
bornea
parents:
diff changeset
91 genename = lst2[0]
a551998e1068 Uploaded
bornea
parents:
diff changeset
92 if 'GN=' not in data_line:
a551998e1068 Uploaded
bornea
parents:
diff changeset
93 genename = 'NA'
a551998e1068 Uploaded
bornea
parents:
diff changeset
94 while ">sp" not in data_lines[match]:
a551998e1068 Uploaded
bornea
parents:
diff changeset
95 if match <= db_len:
a551998e1068 Uploaded
bornea
parents:
diff changeset
96 seqlength = seqlength + len(data_lines[match].strip())
a551998e1068 Uploaded
bornea
parents:
diff changeset
97 if data_lines[match] == last_line:
a551998e1068 Uploaded
bornea
parents:
diff changeset
98 break
a551998e1068 Uploaded
bornea
parents:
diff changeset
99 match = match + 1
a551998e1068 Uploaded
bornea
parents:
diff changeset
100 else:
a551998e1068 Uploaded
bornea
parents:
diff changeset
101 break
a551998e1068 Uploaded
bornea
parents:
diff changeset
102 return ReturnValue1(seqlength, genename)
a551998e1068 Uploaded
bornea
parents:
diff changeset
103 count = count + 1
a551998e1068 Uploaded
bornea
parents:
diff changeset
104 if seqlength == 0:
a551998e1068 Uploaded
bornea
parents:
diff changeset
105 error.write(uniprot_accession_in + '\t' + "Uniprot not in Fasta" + '\n')
a551998e1068 Uploaded
bornea
parents:
diff changeset
106 error.close
a551998e1068 Uploaded
bornea
parents:
diff changeset
107 seqlength = 'NA'
a551998e1068 Uploaded
bornea
parents:
diff changeset
108 genename = 'NA'
a551998e1068 Uploaded
bornea
parents:
diff changeset
109 return ReturnValue1(seqlength, genename)
a551998e1068 Uploaded
bornea
parents:
diff changeset
110 def make_inter(mzIdentML,replicate,grouping):
a551998e1068 Uploaded
bornea
parents:
diff changeset
111 accession_index = mzIdentML[0].index("accession")
a551998e1068 Uploaded
bornea
parents:
diff changeset
112 PSMs = {}
a551998e1068 Uploaded
bornea
parents:
diff changeset
113 accessions = []
a551998e1068 Uploaded
bornea
parents:
diff changeset
114 cnt = 0
a551998e1068 Uploaded
bornea
parents:
diff changeset
115 unique_lines = [mzIdentML[1:]]
a551998e1068 Uploaded
bornea
parents:
diff changeset
116 for i in mzIdentML[1:]:
a551998e1068 Uploaded
bornea
parents:
diff changeset
117 PSMs[i[accession_index]] = 0
a551998e1068 Uploaded
bornea
parents:
diff changeset
118 if i[accession_index] not in accessions:
a551998e1068 Uploaded
bornea
parents:
diff changeset
119 accessions.append(i[accession_index])
a551998e1068 Uploaded
bornea
parents:
diff changeset
120 if i not in unique_lines:
a551998e1068 Uploaded
bornea
parents:
diff changeset
121 unique_lines.append(i)
a551998e1068 Uploaded
bornea
parents:
diff changeset
122 for i in accessions:
a551998e1068 Uploaded
bornea
parents:
diff changeset
123 for j in unique_lines[1:]:
a551998e1068 Uploaded
bornea
parents:
diff changeset
124 if j[accession_index] == i:
a551998e1068 Uploaded
bornea
parents:
diff changeset
125 PSMs[j[accession_index]] +=1
a551998e1068 Uploaded
bornea
parents:
diff changeset
126 inter = ""
a551998e1068 Uploaded
bornea
parents:
diff changeset
127 for i in accessions:
a551998e1068 Uploaded
bornea
parents:
diff changeset
128 inter = inter + replicate + "\t" + grouping + "\t" + i + "\t" + str(PSMs[i]) + "\n"
a551998e1068 Uploaded
bornea
parents:
diff changeset
129 return ReturnValue2(inter,accessions)
a551998e1068 Uploaded
bornea
parents:
diff changeset
130
a551998e1068 Uploaded
bornea
parents:
diff changeset
131
a551998e1068 Uploaded
bornea
parents:
diff changeset
132 files = sys.argv[1]
a551998e1068 Uploaded
bornea
parents:
diff changeset
133 file_list = files.split(",")
a551998e1068 Uploaded
bornea
parents:
diff changeset
134 bait = read_tab(sys.argv[2])
a551998e1068 Uploaded
bornea
parents:
diff changeset
135 make_prey = sys.argv[3]
a551998e1068 Uploaded
bornea
parents:
diff changeset
136 db = sys.argv[4]
a551998e1068 Uploaded
bornea
parents:
diff changeset
137 if db == "None":
a551998e1068 Uploaded
bornea
parents:
diff changeset
138 db = str(ins_path) + "/SwissProt_HUMAN_2015_12.fasta"
a551998e1068 Uploaded
bornea
parents:
diff changeset
139 make_bait = sys.argv[6]
a551998e1068 Uploaded
bornea
parents:
diff changeset
140 bait_bool = sys.argv[7]
a551998e1068 Uploaded
bornea
parents:
diff changeset
141 prey_file = sys.argv[8]
a551998e1068 Uploaded
bornea
parents:
diff changeset
142 bait_out = sys.argv[9]
a551998e1068 Uploaded
bornea
parents:
diff changeset
143 inter_out = sys.argv[10]
a551998e1068 Uploaded
bornea
parents:
diff changeset
144
a551998e1068 Uploaded
bornea
parents:
diff changeset
145 def bait_create(baits, infile):
a551998e1068 Uploaded
bornea
parents:
diff changeset
146 # Verifies the Baits are valid in the Scaffold file and writes the Bait.txt.
a551998e1068 Uploaded
bornea
parents:
diff changeset
147 baits = make_bait.split()
a551998e1068 Uploaded
bornea
parents:
diff changeset
148 i = 0
a551998e1068 Uploaded
bornea
parents:
diff changeset
149 bait_file_tmp = open("bait.txt", "w")
a551998e1068 Uploaded
bornea
parents:
diff changeset
150 order = []
a551998e1068 Uploaded
bornea
parents:
diff changeset
151 bait_cache = []
a551998e1068 Uploaded
bornea
parents:
diff changeset
152 while i < len(baits):
a551998e1068 Uploaded
bornea
parents:
diff changeset
153 if baits[i+2] == "true":
a551998e1068 Uploaded
bornea
parents:
diff changeset
154 T_C = "C"
a551998e1068 Uploaded
bornea
parents:
diff changeset
155 else:
a551998e1068 Uploaded
bornea
parents:
diff changeset
156 T_C = "T"
a551998e1068 Uploaded
bornea
parents:
diff changeset
157 bait_line = baits[i] + "\t" + baits[i+1] + "\t" + T_C + "\n"
a551998e1068 Uploaded
bornea
parents:
diff changeset
158 bait_cache.append(str(bait_line))
a551998e1068 Uploaded
bornea
parents:
diff changeset
159 i = i + 3
a551998e1068 Uploaded
bornea
parents:
diff changeset
160
a551998e1068 Uploaded
bornea
parents:
diff changeset
161 for cache_line in bait_cache:
a551998e1068 Uploaded
bornea
parents:
diff changeset
162 bait_file_tmp.write(cache_line)
a551998e1068 Uploaded
bornea
parents:
diff changeset
163
a551998e1068 Uploaded
bornea
parents:
diff changeset
164 bait_file_tmp.close()
a551998e1068 Uploaded
bornea
parents:
diff changeset
165
a551998e1068 Uploaded
bornea
parents:
diff changeset
166 if bait_bool == 'false':
a551998e1068 Uploaded
bornea
parents:
diff changeset
167 bait_create(make_bait, infile)
a551998e1068 Uploaded
bornea
parents:
diff changeset
168 bait = "bait.txt"
a551998e1068 Uploaded
bornea
parents:
diff changeset
169 else:
a551998e1068 Uploaded
bornea
parents:
diff changeset
170 bait_temp_file = open(sys.argv[9], 'r')
a551998e1068 Uploaded
bornea
parents:
diff changeset
171 bait_cache = bait_temp_file.readlines()
a551998e1068 Uploaded
bornea
parents:
diff changeset
172 bait_file_tmp = open("bait.txt", "wr")
a551998e1068 Uploaded
bornea
parents:
diff changeset
173 for cache_line in bait_cache:
a551998e1068 Uploaded
bornea
parents:
diff changeset
174 bait_file_tmp.write(cache_line)
a551998e1068 Uploaded
bornea
parents:
diff changeset
175 bait_file_tmp.close()
a551998e1068 Uploaded
bornea
parents:
diff changeset
176 bait = "bait.txt"
a551998e1068 Uploaded
bornea
parents:
diff changeset
177
a551998e1068 Uploaded
bornea
parents:
diff changeset
178 inter = ""
a551998e1068 Uploaded
bornea
parents:
diff changeset
179 cnt = 0
a551998e1068 Uploaded
bornea
parents:
diff changeset
180 accessions = []
a551998e1068 Uploaded
bornea
parents:
diff changeset
181 for i in file_list:
a551998e1068 Uploaded
bornea
parents:
diff changeset
182 cmd = (r"Rscript "+ str(ins_path) +"flatten_mzIdentML.R " + i)
a551998e1068 Uploaded
bornea
parents:
diff changeset
183 os.system(cmd)
a551998e1068 Uploaded
bornea
parents:
diff changeset
184 mzIdentML = read_tab("flat_mzIdentML.txt")
a551998e1068 Uploaded
bornea
parents:
diff changeset
185 inter = inter + make_inter(mzIdentML,bait[cnt][0],bait[cnt][1]).inter
a551998e1068 Uploaded
bornea
parents:
diff changeset
186 accessions.append(make_inter(mzIdentML,bait[cnt][0],bait[cnt][1]).accessions)
a551998e1068 Uploaded
bornea
parents:
diff changeset
187 cnt+=1
a551998e1068 Uploaded
bornea
parents:
diff changeset
188
a551998e1068 Uploaded
bornea
parents:
diff changeset
189 with open("inter.txt","w") as x:
a551998e1068 Uploaded
bornea
parents:
diff changeset
190 x.write(inter)
a551998e1068 Uploaded
bornea
parents:
diff changeset
191 if make_prey == "Y":
a551998e1068 Uploaded
bornea
parents:
diff changeset
192 unique_accessions = []
a551998e1068 Uploaded
bornea
parents:
diff changeset
193 prey = ""
a551998e1068 Uploaded
bornea
parents:
diff changeset
194 for i in accessions:
a551998e1068 Uploaded
bornea
parents:
diff changeset
195 for j in i:
a551998e1068 Uploaded
bornea
parents:
diff changeset
196 if j not in unique_accessions:
a551998e1068 Uploaded
bornea
parents:
diff changeset
197 unique_accessions.append(j)
a551998e1068 Uploaded
bornea
parents:
diff changeset
198 start = 0
a551998e1068 Uploaded
bornea
parents:
diff changeset
199 end = len(unique_accessions)
a551998e1068 Uploaded
bornea
parents:
diff changeset
200 printProgress(start,end,prefix = "Making Prey File:",suffix = "Complete",barLength=50)
a551998e1068 Uploaded
bornea
parents:
diff changeset
201
a551998e1068 Uploaded
bornea
parents:
diff changeset
202 for i in unique_accessions:
a551998e1068 Uploaded
bornea
parents:
diff changeset
203 prey = prey + i + "\t" + str(get_info(i,db).seqlength) + "\t" + get_info(i,db).genename + "\n"
a551998e1068 Uploaded
bornea
parents:
diff changeset
204 start+=1
a551998e1068 Uploaded
bornea
parents:
diff changeset
205 printProgress(start, end)
a551998e1068 Uploaded
bornea
parents:
diff changeset
206 with open("prey.txt","w") as x:
a551998e1068 Uploaded
bornea
parents:
diff changeset
207 x.write(prey)
a551998e1068 Uploaded
bornea
parents:
diff changeset
208
a551998e1068 Uploaded
bornea
parents:
diff changeset
209 os.rename("bait.txt", sys.argv[2])
a551998e1068 Uploaded
bornea
parents:
diff changeset
210 os.rename("inter.txt", sys.argv[10])
a551998e1068 Uploaded
bornea
parents:
diff changeset
211 if str(prey) != "None":
a551998e1068 Uploaded
bornea
parents:
diff changeset
212 os.rename("prey.txt", sys.argv[11])