| blob | 182ecdfe71c1ab3fd6ea400440eae6aa04eb57fd |
1 /*\r
2 ******************************************************************************\r
3 *Copyright 2010\r
4 * BGI-SHENZHEN\r
5 *All Rights Reserved\r
6 *ATHUOR : Bill Tang\r
7 *CREATE DATE : 2010-7-14\r
8 *CLASS NAME: \r
9 *FUNCTION : definition of the useful functions\r
10 *FILE NAME : tool.cpp\r
11 *UPDATE DATE : 2010-8-3\r
12 *UPDATE BY : Bill Tang\r
13 UPDATE DATE : 2010-9-9\r
14 *UPDATE BY : Bill Tang\r
15 UPDATE DATE : 2010-9-15\r
16 *UPDATE BY : Zhou Guyue\r
17 *******************************************************************************\r
18 */\r
19 \r
20 #include "tools.h"\r
21 #include <iostream>\r
22 #include <fstream>\r
23 #include <string>\r
24 #include <sstream>\r
25 using namespace std;\r
26 \r
27 /* integer to string*/\r
28 std::string myitoa(int num) {\r
29 char num_str[256] = {0};\r
30 sprintf(num_str, "%d", num);\r
31 return num_str;\r
32 }\r
33 \r
34 // A function to spilt string s into vector vec according to char splitchar\r
35 void StringSplit(std::string s, char splitchar, std::vector<std::string>& vec) {\r
36 // promise the vec is empty\r
37 if(vec.size()>0) \r
38 vec.clear();\r
39 // chomp the string\r
40 while(s.size() > 0 && s[0] == splitchar)\r
41 s = s.substr(1, s.length() - 1);\r
42 while(s.size() > 0 && s[s.length() - 1] == splitchar)\r
43 s = s.substr(0, s.length() - 1);\r
44 int length = s.length();\r
45 int start=0;\r
46 for(int i=0;i<length;i++)\r
47 {\r
48 if(s[i] == splitchar)\r
49 {\r
50 vec.push_back(s.substr(start,i - start));\r
51 while(s[i + 1] == splitchar)\r
52 i ++;\r
53 start = i+1;\r
54 }\r
55 else if(i == length-1)// attach last\r
56 {\r
57 vec.push_back(s.substr(start,i+1 - start));\r
58 }\r
59 }\r
60 }\r
61 \r
62 /* count the indel length from the cigar string. return the indel length and the position*/\r
63 int count_indel_len(const std::string cigar, int &pos) {\r
64 const char *tmp = cigar.c_str();\r
65 int begin = 0;\r
66 int end;\r
67 int len = 0;\r
68 pos = 0;\r
69 for (end = 0; end < cigar.size(); end ++) {\r
70 switch (tmp[end]) {\r
71 case 'M':\r
72 {\r
73 pos += atoi(cigar.substr(begin, end - begin).c_str());\r
74 begin = end + 1;\r
75 break;\r
76 }\r
77 case 'I':\r
78 {\r
79 //len = 200 + atoi(cigar.substr(begin, end - begin).c_str()); \r
80 len = 100 + atoi(cigar.substr(begin, end - begin).c_str()); \r
81 begin = end + 1;\r
82 break;\r
83 }\r
84 case 'D':\r
85 {\r
86 //len = 100 + atoi(cigar.substr(begin, end - begin).c_str()); \r
87 len = 200 + atoi(cigar.substr(begin, end - begin).c_str()); \r
88 begin = end + 1;\r
89 break;\r
90 }\r
91 case 'N':\r
92 {\r
93 pos += atoi(cigar.substr(begin, end - begin).c_str());\r
94 begin = end + 1;\r
95 break;\r
96 }\r
97 case 'S':\r
98 {\r
99 pos += atoi(cigar.substr(begin, end - begin).c_str());\r
100 begin = end + 1;\r
101 break;\r
102 }\r
103 case 'H':\r
104 {\r
105 pos += atoi(cigar.substr(begin, end - begin).c_str());\r
106 begin = end + 1;\r
107 break;\r
108 }\r
109 case 'P':\r
110 {\r
111 pos += atoi(cigar.substr(begin, end - begin).c_str());\r
112 begin = end + 1;\r
113 break;\r
114 }\r
115 default : break;\r
116 }\r
117 if (len != 0) {\r
118 break;\r
119 }\r
120 }\r
121 return len;\r
122 }\r
123 \r
124 /* format the sam text to the soap text*/\r
125 std::string alignment_format(const std::string &sam_ali) {\r
126 if(sam_ali.empty()) {\r
127 return NOUSE_ALIGNMENT;\r
128 }\r
129 std::vector<std::string> vec;\r
130 StringSplit(sam_ali, '\t', vec);\r
131 if ((vec.size() < 11) || (vec[5] == "*"))\r
132 return NOUSE_ALIGNMENT;\r
133 \r
134 std::string format;\r
135 int pos = 0;\r
136 int best_hit = 1;\r
137 int flag = atoi(vec[1].c_str());\r
138 int last_clip_num;\r
139 int clip_num = count_soft_clip(vec[5], last_clip_num);\r
140 \r
141 if (clip_num + last_clip_num > vec[9].size())\r
142 return NOUSE_ALIGNMENT;\r
143 \r
144 if (flag & (0x1 << 6)) {\r
145 format = vec[0] + "/1" + "\t"; // add query name\r
146 } else if (flag & (0x1 << 7)) {\r
147 format = vec[0] + "/2" + "\t"; // add query name\r
148 } else {\r
149 format = vec[0] + "\t"; // add query name\r
150 }\r
151 \r
152 format += vec[9].substr(clip_num, vec[9].size() - last_clip_num - clip_num) + "\t"; // add sequence\r
153 format += vec[10].substr(clip_num, vec[10].size() - last_clip_num - clip_num) + "\t"; // add quality\r
154 \r
155 for (int i = 11; i < vec.size(); i++) {\r
156 if (vec[i].find(HIT_COUNT_H0) != std::string::npos) {\r
157 best_hit = atoi(vec[i].substr(5).c_str());\r
158 break;\r
159 } else if (vec[i].find(HIT_COUNT_X0) != std::string::npos) {\r
160 best_hit = atoi(vec[i].substr(5).c_str());\r
161 break;\r
162 }\r
163 }\r
164 format += myitoa(best_hit) + "\t"; // add number of best hit.\r
165 \r
166 format += ((flag & (0x1 << 6)) ? "a" : "b"); // add a/b\r
167 format += "\t";\r
168 format += myitoa(vec[9].size() - clip_num - last_clip_num); // add length\r
169 format += "\t";\r
170 format += ((flag & (0x1 << 4)) ? "-" : "+"); // add strand\r
171 format += "\t";\r
172 format += vec[2] + "\t"; // add chr\r
173 format += vec[3] + "\t"; // add location\r
174 format += myitoa(count_indel_len(vec[5], pos)); // add number of mismatch\r
175 format += "\t";\r
176 format += myitoa(pos - clip_num); // add indel position\r
177 //format += myitoa(pos) + "\t"; // add indel position. changed at 2010-11-22, to compate to the new version soap.\r
178 //format += vec[5];\r
179 return format;\r
180 }\r
181 \r
182 /**\r
183 * DATE: 2010-9-9\r
184 * FUNCTION: count the soft clip number at the read's beginning.\r
185 * PARAMETER: cigar: the cigar string. last_S: last soft clip's length.\r
186 * RETURN: soft clip number.\r
187 */\r
188 int count_soft_clip(const std::string cigar, int &last_S)\r
189 {\r
190 const char *tmp = cigar.c_str();\r
191 last_S = 0;\r
192 // count the last soft clip number.\r
193 if (tmp[cigar.size() - 1] == 'S')\r
194 {\r
195 int j = cigar.size() - 2;\r
196 while (tmp[j] >= '0' && tmp[j] <= '9')\r
197 {\r
198 j--;\r
199 }\r
200 last_S = atoi(cigar.substr(j + 1, cigar.size() - j - 1).c_str());\r
201 }\r
202 \r
203 int i = 0;\r
204 for (; i < cigar.size(); ++i)\r
205 {\r
206 if (tmp[i] >= '0' && tmp[i] <= '9')\r
207 {\r
208 // jump to the next position.\r
209 continue;\r
210 }\r
211 else if (tmp[i] == 'S')\r
212 {\r
213 // return the soft clip number.\r
214 int num = atoi(cigar.substr(0, i).c_str());\r
215 return num;\r
216 }\r
217 else\r
218 {\r
219 // no soft clip infront of the reads.\r
220 break;\r
221 }\r
222 }\r
223 return 0;\r
224 }\r
225 \r
226 /**\r
227 * DATE: 2010-9-7\r
228 * FUNCTION: change string to bool\r
229 * PARAMETER: str: a string which will change to bool. len: the str length.\r
230 * RETURN: bool type of str\r
231 */\r
232 int str_bool(const char* str)\r
233 {\r
234 if ( *str == '0')\r
235 {\r
236 return 0;\r
237 } \r
238 else\r
239 { \r
240 return 1;\r
241 }\r
242 }\r
243 \r
244 /**\r
245 * DATE: 2010-9-20\r
246 * FUNCTION: intialize sfs parament\r
247 * PARAMETER: sfs_path: sfs configuration file path. sfs: a sfs_para struct.\r
248 * RETURN: SFS_PARA_ERROR:if sfs parament is wrong ,SFS_PARA_SUCC : if sfs parament is right.\r
249 * UPDATE: 2010-10-15\r
250 */\r
251 int init_sfs_para(const string sfs_path, SFS_PARA *sfs )\r
252 {\r
253 string in;\r
254 \r
255 /*set default vaule for sfs_par struct*/\r
256 sfs->allow_PathZ = 0;\r
257 \r
258 sfs->bias = 0.03;\r
259 sfs->doBay = 0;\r
260 sfs->doJoint = 0;\r
261 sfs->eps = 0.005;\r
262 sfs->outfiles = "";\r
263 sfs->pThres = 0.01;\r
264 sfs->pvar = 0.015;\r
265 sfs->sigle_MB = 0;\r
266 sfs->sigle_MJ = 0;\r
267 sfs->start_pos = 0;\r
268 sfs->stop_pos = 0;\r
269 sfs->under_FP = 0;\r
270 sfs->writeFr = 0; \r
271 sfs->alternative = 0;\r
272 sfs->jointFold = 1;\r
273 sfs->qs = 52;\r
274 //update 11-29 for control the column\r
275 sfs->sfs_first_last = 0 ;\r
276 \r
277 my_ifstream sfsfile(sfs_path.c_str()); //open sfs configuration file\r
278 \r
279 if (!sfsfile.is_open()) //can not open sfsfile\r
280 {\r
281 cerr << "Cannot open file:" << sfs_path << endl;\r
282 return SFS_PARA_ERROR;\r
283 }\r
284 \r
285 stringstream temp_ss;\r
286 while (sfsfile >> in)\r
287 {\r
288 if (in.empty())\r
289 {\r
290 continue;\r
291 }\r
292 if (in == "doBay")\r
293 {\r
294 sfsfile >> in ;\r
295 sfs->doBay = str_bool(in.c_str());\r
296 continue;\r
297 }\r
298 else if (in == "doJoint")\r
299 {\r
300 sfsfile >> in ;\r
301 sfs->doJoint = str_bool(in.c_str());\r
302 continue;\r
303 }\r
304 else if (in == "writeFr")\r
305 {\r
306 sfsfile >> in ;\r
307 sfs->writeFr = str_bool(in.c_str());\r
308 continue;\r
309 }\r
310 else if (strcmp(in.c_str() , "outfiles") == 0)\r
311 {\r
312 sfsfile >> sfs->outfiles;\r
313 continue;\r
314 }\r
315 else if (in == "start")\r
316 {\r
317 sfsfile >> in ;\r
318 temp_ss << in;\r
319 temp_ss >> sfs->start_pos;\r
320 temp_ss.clear();\r
321 continue;\r
322 }\r
323 else if (in == "stop")\r
324 {\r
325 sfsfile >> in ;\r
326 temp_ss << in;\r
327 temp_ss >> sfs->stop_pos;\r
328 temp_ss.clear();\r
329 continue;\r
330 }\r
331 else if (in == "underFlowProtect")\r
332 {\r
333 sfsfile >> in ;\r
334 sfs->under_FP = str_bool(in.c_str());\r
335 continue;\r
336 }\r
337 else if (in == "allowPhatZero")\r
338 {\r
339 sfsfile >> in ;\r
340 sfs->allow_PathZ = str_bool(in.c_str());\r
341 continue;\r
342 }\r
343 else if ( in == "bias")\r
344 {\r
345 sfsfile >> in ;\r
346 temp_ss << in;\r
347 temp_ss >> sfs->bias;\r
348 temp_ss.clear();\r
349 continue;\r
350 }\r
351 else if (in == "pvar")\r
352 {\r
353 sfsfile >> in ;\r
354 temp_ss << in;\r
355 temp_ss >> sfs->pvar;\r
356 temp_ss.clear();\r
357 continue;\r
358 }\r
359 else if (in == "eps")\r
360 {\r
361 sfsfile >> in ;\r
362 temp_ss << in;\r
363 temp_ss >> sfs->eps;\r
364 temp_ss.clear();\r
365 continue;\r
366 }\r
367 else if (in == "pThres")\r
368 {\r
369 sfsfile >> in ;\r
370 temp_ss << in;\r
371 temp_ss >> sfs->pThres;\r
372 temp_ss.clear();\r
373 continue;\r
374 }\r
375 else if (in == "sigleMinorJoint")\r
376 {\r
377 sfsfile >> in ;\r
378 sfs->sigle_MJ = str_bool(in.c_str());\r
379 continue;\r
380 }\r
381 else if (in == "sigleMinorBay" )\r
382 {\r
383 sfsfile >> in ;\r
384 sfs->sigle_MB = str_bool(in.c_str());\r
385 continue;\r
386 }\r
387 else if (in == "alt" )\r
388 {\r
389 sfsfile >> in ;\r
390 sfs->alternative = atoi(in.c_str());\r
391 continue;\r
392 }\r
393 else if (in == "qs" )\r
394 {\r
395 sfsfile >> sfs->qs;\r
396 continue;\r
397 }\r
398 //else if (in == "sfsfirstlast")\r
399 else if ( in == "lightoutput")\r
400 {\r
401 sfsfile >> in ;\r
402 sfs->sfs_first_last = str_bool(in.c_str());\r
403 continue;\r
404 }\r
405 else\r
406 {\r
407 sfsfile.close();\r
408 return SFS_PARA_ERROR;\r
409 }\r
410 }\r
411 \r
412 sfsfile.close();\r
413 return SFS_PARA_SUCC;\r
414 }\r
415 \r
416 /*sfs help information*/\r
417 /**\r
418 * DATE: 2010-9-7\r
419 * FUNCTION: return help information\r
420 * PARAMETER: void\r
421 * RETURN: void\r
422 */\r
423 void sfs_info()\r
424 { \r
425 cerr << "\t-> (outfiles) -outfiles \n" << endl;\r
426 cerr << "\t-> (which analysis) -doBay -doJoint\n" << endl;\r
427 cerr << "\t-> (optional variables) -eps -bias -pThres -underFlowProtect\n" << endl;\r
428 cerr << "\t-> (optional variables) -start -stop \t(format is chromo:position)\n" << endl;\r
429 cerr << "\t-> (optional variables) -singleMinor[Joint,Bay] [0 or 1] should we only use a singleminor?\n" << endl;\r
430 cerr << "\t-> (optional variables) -sfsfirstlast [0 or 1] should we only output sfs the first column and the last column \n" <<endl;\r
431 }\r