Mercurial > repos > bornea > saint_preprocessing_peptideshaker

comparison Protein_report_processing.py @ 0:5a67d55bb091 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author bornea
date Tue, 23 Aug 2016 18:55:17 -0400
parents
children 39e772aa6e93
comparison
equal deleted inserted replaced
-1:000000000000 0:5a67d55bb091
1 import sys
2 import os
3 from time import sleep
4
5 files = sys.argv[1] # read in a string of file names seperated by ", "
6 # e.g. "Default_Protein_Report.txt, Default_Protein_Report_2.txt"
7 #bait = sys.argv[2] # SAINT formatted bait file
8 # still need a way to match files to bait identifiers
9 # or they can just be required to be put in the order of the bait file
10 quant_type = sys.argv[3] # what metric to use for quantification
11 # "#Validated Peptides", "#Peptides", "#Unique", "#Validated PSMs", "#PSMs"
12 db = sys.argv[4] # fasta database used in SearchGUI and PeptideShaker
13 prey = sys.argv[5]
14 tool_path = sys.argv[7]
15 if fasta_db == "None":
16 fasta_db = str(tool_path) + "/SwissProt_HUMAN_2015_12.fasta"
17 make_bait = sys.argv[6]
18 bait_bool = sys.argv[8]
19
20 def bait_create(baits, infile):
21 # Verifies the Baits are valid in the Scaffold file and writes the Bait.txt.
22 baits = make_bait.split()
23 i = 0
24 bait_file_tmp = open("bait.txt", "w")
25 order = []
26 bait_cache = []
27 while i < len(baits):
28 if baits[i+2] == "true":
29 T_C = "C"
30 else:
31 T_C = "T"
32 bait_line = baits[i] + "\t" + baits[i+1] + "\t" + T_C + "\n"
33 bait_cache.append(str(bait_line))
34 i = i + 3
35
36 for cache_line in bait_cache:
37 bait_file_tmp.write(cache_line)
38
39 bait_file_tmp.close()
40
41 if bait_bool == 'false':
42 bait_create(make_bait, infile)
43 bait = "bait.txt"
44 else:
45 bait_temp_file = open(sys.argv[9], 'r')
46 bait_cache = bait_temp_file.readlines()
47 bait_file_tmp = open("bait.txt", "wr")
48 for cache_line in bait_cache:
49 bait_file_tmp.write(cache_line)
50 bait_file_tmp.close()
51 bait = "bait.txt"
52
53 class ReturnValue1(object):
54 def __init__(self, sequence, gene):
55 self.seqlength = sequence
56 self.genename = gene
57
58 def read_tab(infile):
59 with open(infile,'r') as x:
60 output = []
61 for line in x:
62 line = line.strip()
63 temp = line.split('\t')
64 output.append(temp)
65 return output
66 def printProgress (iteration, total, prefix = '', suffix = '', decimals = 1, barLength = 100):
67 """
68 Call in a loop to create terminal progress bar
69 @params:
70 iteration - Required : current iteration (Int)
71 total - Required : total iterations (Int)
72 prefix - Optional : prefix string (Str)
73 suffix - Optional : suffix string (Str)
74 decimals - Optional : positive number of decimals in percent complete (Int)
75 barLength - Optional : character length of bar (Int)
76 """
77 formatStr = "{0:." + str(decimals) + "f}"
78 percents = formatStr.format(100 * (iteration / float(total)))
79 filledLength = int(round(barLength * iteration / float(total)))
80 bar = '=' * filledLength + '-' * (barLength - filledLength)
81 sys.stdout.write('\r%s |%s| %s%s %s' % (prefix, bar, percents, '%', suffix)),
82 sys.stdout.flush()
83 if iteration == total:
84 sys.stdout.write('\n')
85 sys.stdout.flush()
86 def get_info(uniprot_accession_in,fasta_db):
87 # Get aminoacid lengths and gene name.
88 error = open('error proteins.txt', 'a+')
89 data = open(fasta_db, 'r')
90 data_lines = data.readlines()
91 db_len = len(data_lines)
92 seqlength = 0
93 count = 0
94 for data_line in data_lines:
95 if ">sp" in data_line:
96 namer = data_line.split("|")[2]
97 if uniprot_accession_in == data_line.split("|")[1]:
98 match = count + 1
99 if 'GN=' in data_line:
100 lst = data_line.split('GN=')
101 lst2 = lst[1].split(' ')
102 genename = lst2[0]
103 if 'GN=' not in data_line:
104 genename = 'NA'
105 while ">sp" not in data_lines[match]:
106 if match <= db_len:
107 seqlength = seqlength + len(data_lines[match].strip())
108 match = match + 1
109 else:
110 break
111 return ReturnValue1(seqlength, genename)
112 if uniprot_accession_in == namer.split(" ")[0]:
113 match = count + 1
114 # Ensures consistent spacing throughout.
115 if 'GN=' in data_line:
116 lst = data_line.split('GN=')
117 lst2 = lst[1].split(' ')
118 genename = lst2[0]
119 if 'GN=' not in data_line:
120 genename = 'NA'
121 while ">sp" not in data_lines[match]:
122 if match <= db_len:
123 seqlength = seqlength + len(data_lines[match].strip())
124 match = match + 1
125 else:
126 break
127 return ReturnValue1(seqlength, genename)
128 count = count + 1
129 if seqlength == 0:
130 error.write(uniprot_accession_in + '\t' + "Uniprot not in Fasta" + '\n')
131 error.close
132 seqlength = 'NA'
133 genename = 'NA'
134 return ReturnValue1(seqlength, genename)
135 def concatenate_files(file_list_string, bait_file):
136 file_list = file_list_string.split(", ")
137 bait = read_tab(bait_file)
138 master_table = []
139 header_check = 0
140 file_cnt = 0
141 table_cnt = 0
142 for i in file_list:
143 table = read_tab(i)
144 for j in table:
145 if table_cnt == 0:
146 if header_check == 0:
147 header_check +=1
148 j.append("Replicate")
149 j.append("Bait_Grouping")
150 master_table.append(j)
151 if table_cnt > 0:
152 j.append(bait[file_cnt][0])
153 j.append(bait[file_cnt][1])
154 master_table.append(j)
155 table_cnt +=1
156 file_cnt+=1
157 table_cnt = 0
158 if len(master_table[0]) < len(master_table[1]):
159 master_table[0] = ["#"] + master_table[0]
160 with open("merged_PeptideShaker.txt","w") as x:
161 for i in master_table:
162 x.write("\t".join(i))
163 x.write("\n")
164 return master_table
165 def make_inter(master_table,quant_type):
166 if len(master_table[0]) < len(master_table[1]):
167 master_table[0] = ["#"] + master_table[0]
168 replicate_index = master_table[0].index("Replicate")
169 grouping_index = master_table[0].index("Bait_Grouping")
170 accession_index = master_table[0].index("Main Accession")
171 Quant_index = master_table[0].index(quant_type)
172 inter_file = ""
173 for i in master_table[1:]:
174 line = []
175 line.append(i[replicate_index])
176 line.append(i[grouping_index])
177 line.append(i[accession_index])
178 line.append(i[Quant_index])
179 inter_file = inter_file + "\t".join(line) + "\n"
180 with open("inter.txt","w") as x:
181 x.write(inter_file)
182
183 def make_prey(concat_table,fasta_db):
184 input_data = concat_table
185 if len(input_data[0]) < len(input_data[1]):
186 input_data[0] = ["#"] + input_data[0]
187 accession_index = input_data[0].index("Main Accession")
188 proteins = []
189 for i in input_data[1:]:
190 proteins.append(i[accession_index])
191 output_file = open("prey.txt", 'w')
192 start = 0
193 end = len(proteins)
194
195 # Initial call to print 0% progress
196 printProgress(start, end, prefix = 'Progress:', suffix = 'Complete', barLength = 50)
197
198 for protein in proteins:
199 seq = get_info(protein,fasta_db).seqlength
200 GN = get_info(protein,fasta_db).genename
201 if seq != 'NA':
202 output_file.write(protein + "\t" + str(seq) + "\t" + str(GN) + "\n")
203 start+=1
204 printProgress(start, end, prefix = 'Progress:', suffix = 'Complete', barLength = 50)
205 output_file.close()
206 data = concatenate_files(files,bait)
207 make_inter(data, quant_type)
208 if prey == "true":
209 make_prey(data,db)
210
211 os.rename("bait.txt", sys.argv[2])
212 os.rename("inter.txt", sys.argv[10])
213 os.rename("prey.txt", sys.argv[11])