SOX12 and DELTA mediate ecdysone-driven adult neural cell proliferation during metamorphosis

1. #include "process_argv.h"
2. #include <stdlib.h>
3. #include <iostream>
4. #include <stdio.h>
5. #include <getopt.h>
6. #include <map>
7. #include <fstream>
8. #include "global_parameter.h"
9. #include "sys/sysinfo.h"
10. // #include <string>
11. using namespace ::std;
12. #define ADA_RATIO 0.4
14. map<string, vector<string>> wrong_paras;
16. void check_module(int argc, char *argv[])
17. {
18. if (argc < 2)
19. {
20. printModule();
21. }
22. else
23. {
24. set<string> modules;
25. modules.insert("filter");
26. modules.insert("filtersRNA");
27. modules.insert("filterMeta");
28. #ifdef _PROCESSHTS
29. modules.insert("filterHts");
30. #endif
31. modules.insert("filterStLFR");
32. string module = argv[1];
33. if (modules.find(module) == modules.end())
34. {
35. if (module == "-h" || module == "--help")
36. {
37. printModule();
38. }
39. else if (module == "-v" || module == "--version")
40. {
41. printVersion();
42. }
43. else
44. {
45. if (module == "filterHts")
46. {
47. cout
48. << "Error:filterHts module not installed, please re-make after set USEHTS true in Makefile, details see Install part in Readme"
49. << endl;
50. }
51. cerr << "Error:no such module,type -h/--help for help" << endl;
52. exit(1);
53. }
54. }
55. else
56. {
57. if (argc == 2)
58. {
59. if (module == "filterHts")
60. {
61. printHtsUsage();
62. }
63. else
64. {
65. printUsage(module);
66. }
67. }
68. }
69. }
70. }
72. int global_parameter_initial(int argc, char *argv[], C_global_parameter &gp)
73. {
74. // const char *shortOptions = "f:r:1:2:K:M:A:l:T:q:n:m:p:d3in:N:t:e:c:SO:P:Q:L:I:G:a:o:C:D:R:W:5:6:7:8:9:Eb:x:y:z:hv";
75. string c_module(argv[1]);
76. // const char *shortOptions ="j1:2:R:W:C:D:o:5:8:E:Jaf:r:Z:z:c:d:k:Y:K:F:iQ:G:l:q:m:x:y:n:p:g:X:t:B:O:P:7e:T:3:4:L:w:M:A:9:S:s:U:u:b:0:hv";
77. const char *shortOptions = "j1:2:C:D:o:c:E:Jf:r:l:q:m:x:y:n:p:g:X:t:T:3:4:L:w:hv";
78. // move some uncommonly used parameters to config file
79. // include RW58aZzYcdkKFiBOP7eMA9SsUub0QG
80. const struct option longOptions[] = {
81. // common parameter
82. // input and output file
83. //{"mode",1,NULL,'E'},
84. {"streaming", 0, NULL, 'j'},
85. {"fq1", 1, NULL, '1'},
86. {"fq2", 1, NULL, '2'},
88. // { "trimFq1" , 1, NULL, 'R' },//uncommonly used parameters
89. // { "trimFq2" , 1, NULL, 'W' },//uncommonly used parameters
91. {"cleanFq1", 1, NULL, 'C'},
92. {"cleanFq2", 1, NULL, 'D'},
93. {"outDir", 1, NULL, 'o'},
94. {"configFile", 1, NULL, 'c'},
96. // input and output file type
97. // { "seqType" , 1, NULL, '5' },//uncommonly used parameters
98. // {"outFileType",1,NULL,'8'},//uncommonly used parameters
100. // reference for cram input
101. {"ref", 1, NULL, 'E'},
102. {"ada_trim", 0, NULL, 'J'},
103. // {"contam_trim",0,NULL,'a'},//uncommonly used parameters
104. {"adapter1", 1, NULL, 'f'},
105. {"adapter2", 1, NULL, 'r'},
106. // { "contam1" , 1, NULL, 'Z' },//uncommonly used parameters
107. // { "contam2" , 1, NULL, 'z' },//uncommonly used parameters
109. // { "ctMatchR" , 1,NULL,'Y'},//uncommonly used parameters
110. //
111. // { "global_contams",1,NULL,'c'},//uncommonly used parameters
112. // //uncommonly used parameters
113. // { "glob_cotm_mR",1,NULL,'d'},//uncommonly used parameters
114. // { "glob_cotm_mM",1,NULL,'k'},//uncommonly used parameters
115. // { "tile" , 1, NULL, 'K' },//uncommonly used parameters
116. // { "fov" , 1, NULL, 'F' },//uncommonly used parameters
117. // //index remove
118. // { "index" , 0, NULL, 'i' },//uncommonly used parameters
119. //
120. // //base quality
121. // { "qualSys" , 1, NULL, 'Q' },//uncommonly used parameters
122. // { "outQualSys" , 1, NULL, 'G' },//uncommonly used parameters
123. {"lowQual", 1, NULL, 'l'},
124. {"qualRate", 1, NULL, 'q'},
125. {"mean", 1, NULL, 'm'},
126. {"trimBadHead", 1, NULL, 'x'},
127. {"trimBadTail", 1, NULL, 'y'},
128. // base content
129. {"nRate", 1, NULL, 'n'},
130. {"highA", 1, NULL, 'p'},
131. {"polyG_tail", 1, NULL, 'g'},
132. {"polyX", 1, NULL, 'X'},
133. {"trim", 1, NULL, 't'},
134. // { "TtoU" , 0, &TtoU, 1 },
135. // { "UtoT" , 0, &UtoT, 1 },
137. // {"baseConvert",1,NULL,'B'},//uncommonly used parameters
138. // //{ "small" , 0, NULL, 'S' },
139. // //PE reads
140. // { "overlap" , 1, NULL, 'O' },//uncommonly used parameters
141. // { "mis" , 1, NULL, 'P' },//uncommonly used parameters
142. // { "pe_info" , 0, NULL, '7' },//uncommonly used parameters
143. // //computer resource
144. // { "patch" , 1, NULL, 'e' },//uncommonly used parameters
146. {"thread", 1, NULL, 'T'},
147. //{ "split_line",1,NULL, '6'},
148. // read length limit
149. // {"maxReadLen",1,NULL,'3'},//uncommonly used parameters
150. {"minReadLen", 1, NULL, '4'},
151. // reads number limit
152. // {"'totalReadsNum'" , 1, NULL, 'c' },
153. // {"totalReadsNum",1,NULL,'L'},
154. {"output_clean", 1, NULL, 'w'},
156. // {"adaMis",1,NULL,'M'},//uncommonly used parameters
157. // {"adaMR",1,NULL,'A'},//uncommonly used parameters
158. // {"adaEdge",1,NULL,'9'},//uncommonly used parameters
159. //
160. // {"adaRCtg",1,NULL,'S'},//uncommonly used parameters
161. // {"adaRAr",1,NULL,'s'},//uncommonly used parameters
162. // {"adaRMa",1,NULL,'U'},//uncommonly used parameters
163. // {"adaREr",1,NULL,'u'},//uncommonly used parameters
164. // {"adaRMm",1,NULL,'b'},//uncommonly used parameters
166. // { "append" , 1, NULL, 'a' },
167. // { "log",1,NULL,'0'},//uncommonly used parameters
168. {"help", 0, NULL, 'h'},
169. {"version", 0, NULL, 'v'},
170. {0, 0, 0, 0}};
172. int nextOpt;
173. gp.module_name = argv[1];
174. if (gp.module_name == "filtersRNA")
175. {
176. gp.min_read_length = 18;
177. gp.max_read_length = 49;
178. }
179. gp.log = "log";
180. // int error=0;
181. while (-1 != (nextOpt = getopt_long(argc, argv, shortOptions, longOptions, NULL)))
182. {
183. switch (nextOpt)
184. {
185. case 'E':
186. gp.reference.assign(optarg);
187. break;
188. case 'j':
189. gp.is_streaming = true;
190. break;
191. case '1':
192. {
193. gp.fq1_path.assign(optarg);
194. if (gp.fq1_path.rfind(".gz") == gp.fq1_path.size() - 3)
195. {
196. gp.inputGzformat = true;
197. }
198. else
199. {
200. gp.inputGzformat = false;
201. }
202. break;
203. }
204. case '2':
205. gp.fq2_path.assign(optarg);
206. break;
207. // case 'R':{
208. // gp.trim_fq1.assign(optarg);
209. // if(gp.trim_fq1.rfind(".gz")==gp.trim_fq1.size()-3){
210. // gp.trimOutGzformat=true;
211. // }else{
212. // gp.trimOutGzformat=false;
213. // }
214. // break;
215. // }
216. // case 'W':gp.trim_fq2.assign(optarg);break;
217. case 'C':
218. {
219. gp.clean_fq1.assign(optarg);
220. if (gp.clean_fq1.rfind(".gz") == gp.clean_fq1.size() - 3)
221. {
222. gp.cleanOutGzFormat = true;
223. }
224. else
225. {
226. gp.cleanOutGzFormat = false;
227. }
228. break;
229. }
230. case 'D':
231. gp.clean_fq2.assign(optarg);
232. break;
233. case 'o':
234. gp.output_dir.assign(optarg);
235. break;
236. // case '5':gp.seq_type.assign(optarg);break;
237. // case '8':gp.output_file_type.assign(optarg);break;
238. case 'J':
239. gp.adapter_discard_or_trim = "trim";
240. break;
241. // case 'a':gp.contam_discard_or_trim="trim";break;
242. case 'f':
243. {
244. string validSeq = "ACGTacgtNn";
245. ifstream ifAda1(optarg);
246. if (!ifAda1)
247. {
248. gp.adapter1_seq.assign(optarg);
249. gp.ada1s.push_back(gp.adapter1_seq);
250. for (int i = 0; i < gp.adapter1_seq.size(); i++)
251. {
252. if (validSeq.find(gp.adapter1_seq[i]) == string::npos)
253. {
254. cerr << "Error:invalid character found in adapter:" << gp.adapter1_seq[i]
255. << ". Only ACGTacgtNn are supported" << endl;
256. exit(1);
257. }
258. }
259. break;
260. }
261. else
262. {
263. cout << "input adapter1 list file:" << optarg << endl;
264. string adaLine;
265. while (getline(ifAda1, adaLine))
266. {
267. gp.ada1s.push_back(adaLine);
268. }
269. }
270. ifAda1.close();
271. break;
272. }
273. case 'r':
274. {
275. // gp.adapter2_seq.assign(optarg);break;
276. string validSeq = "ACGTacgtNn";
277. ifstream ifAda2(optarg);
278. if (!ifAda2)
279. {
280. gp.adapter2_seq.assign(optarg);
281. gp.ada2s.push_back(gp.adapter2_seq);
282. for (int i = 0; i < gp.adapter2_seq.size(); i++)
283. {
284. if (validSeq.find(gp.adapter2_seq[i]) == string::npos)
285. {
286. cerr << "Error:invalid character found in adapter:" << gp.adapter2_seq[i]
287. << ". Only ACGTacgtNn are supported" << endl;
288. exit(1);
289. }
290. }
291. break;
292. }
293. else
294. {
295. cout << "input adapter2 list file:" << optarg << endl;
296. string adaLine;
297. while (getline(ifAda2, adaLine))
298. {
299. gp.ada2s.push_back(adaLine);
300. }
301. }
302. ifAda2.close();
303. break;
304. }
305. // case 'Z':gp.contam1_seq.assign(optarg);break;
306. case 'z':
307. gp.contam2_seq.assign(optarg);
308. break;
309. // case 'c':gp.global_contams.assign(optarg);break;
310. case 'c':
311. initFromConfigFile(gp, optarg);
312. break; //-c means config file path ,no longer for global contams
313. // case 'd':gp.g_mrs.assign(optarg);break;
314. // case 'k':gp.g_mms.assign(optarg);break;
315. // case 'Y':gp.ctMatchR.assign(optarg);break;
316. // case 'K':gp.tile.assign(optarg);break;
317. // case 'F':gp.fov.assign(optarg);break;
318. // case 'i':gp.index_remove = true;break;
319. // case 'Q':{
320. // gp.qualityPhred=atoi(optarg);
321. // if(gp.qualityPhred==1){
322. // gp.qualityPhred=64;
323. // }else if(gp.qualityPhred==2){
324. // gp.qualityPhred=33;
325. // }
326. // break;
327. // }
328. // case 'G':{
329. // gp.outputQualityPhred=atoi(optarg);
330. // if(gp.outputQualityPhred==1){
331. // gp.outputQualityPhred=64;
332. // }else if(gp.outputQualityPhred==2){
333. // gp.outputQualityPhred=33;
334. // }
335. // break;
336. // }
337. case 'l':
338. gp.lowQual = atoi(optarg);
339. break;
340. case 'q':
341. gp.lowQualityBaseRatio = atof(optarg);
342. break;
343. case 'm':
344. gp.meanQuality = atoi(optarg);
345. break;
346. case 'x':
347. gp.trimBadHead.assign(optarg);
348. break;
349. case 'y':
350. gp.trimBadTail.assign(optarg);
351. break;
352. case 'n':
353. gp.n_ratio = atof(optarg);
354. break;
355. case 'p':
356. gp.highA_ratio = atof(optarg);
357. break;
358. case 'g':
359. gp.polyG_tail = atof(optarg);
360. break;
361. case 'X':
362. gp.polyX_num = atof(optarg);
363. break;
364. case 't':
365. gp.trim.assign(optarg);
366. break;
367. // case 'B':gp.base_convert.assign(optarg);break;
368. // case 'O':gp.overlap_length=atoi(optarg);break;
369. // case 'P':gp.peMismatchRatio=atof(optarg);break;
370. // case '7':gp.whether_add_pe_info=true;break;
371. // case 'e':gp.patchSize=atoi(optarg);break;
372. case 'T':
373. gp.threads_num = atoi(optarg);
374. break;
375. // case '6':gp.split_line=atoi(optarg);break;
376. // case '3':gp.max_read_length=atoi(optarg);break;
377. case '4':
378. gp.min_read_length = atoi(optarg);
379. break;
380. // case 'L':{
381. // string tmp_str;
382. // tmp_str.assign(optarg);
383. // if(tmp_str.find("head")==string::npos){
384. // gp.total_reads_num_random=true;
385. // //gp.catWhenrunning=false;
386. // for(int i=0;i!=tmp_str.size();i++){
387. // if(!isdigit(tmp_str[i]) && tmp_str[i]!='.'){
388. // cerr<<"Error:-L value should be a positive integer or float"<<endl;
389. // exit(1);
390. // }
391. // }
392. // }else{
393. // gp.total_reads_num_random=false;
394. // tmp_str.erase(tmp_str.find("head"),4);
395. // if(tmp_str.find(".")!=string::npos){
396. // cerr<<"Error:-L value should be a integer when with head suffix"<<endl;
397. // exit(1);
398. // }else{
399. // for(int i=0;i!=tmp_str.size();i++){
400. // if(!isdigit(tmp_str[i])){
401. // cerr<<"Error:-L value should be an integer when with head suffix"<<endl;
402. // exit(1);
403. // }
404. // }
405. // }
406. // }
407. //
408. // float tmp_val=atof(optarg);
409. // if(tmp_val==0){
410. // cerr<<"Error:-L value should be a positive integer or float"<<endl;
411. // exit(1);
412. // }
413. // gp.total_reads_num=tmp_val;
414. // if(tmp_val<1){
415. // gp.f_total_reads_ratio=tmp_val;
416. // }else{
417. // istringstream is_str(tmp_str);
418. // is_str>>gp.l_total_reads_num;
419. // }
420. // if(gp.f_total_reads_ratio>0 && gp.l_total_reads_num>0){
421. // cerr<<"Error:reads number and ratio should not be both assigned at the same time"<<endl;
422. // exit(1);
423. // }
424. // break;
425. // }
426. case 'w':
427. {
428. string paraCheck(optarg);
429. for (int i = 0; i != paraCheck.size(); i++)
430. {
431. if (!isdigit(paraCheck[i]))
432. {
433. cerr << "Error:-w value should be a positive integer" << endl;
434. exit(1);
435. }
436. }
437. gp.cleanOutSplit = atoi(paraCheck.c_str());
438. if (gp.cleanOutSplit == 0)
439. {
440. cerr << "Error:-w value should be a positive integer" << endl;
441. exit(1);
442. }
443. break;
444. }
445. // case 'M':{
446. // wrong_paras["filtersRNA"].emplace_back("-M|--adaMis");
447. // string tmp_str;
448. // tmp_str.assign(optarg);
449. // if(tmp_str.find(",")==string::npos){
450. // gp.adaMis=atoi(optarg);
451. // gp.adaMis2=gp.adaMis;
452. // }else{
453. // vector<string> values;
454. // line_split(tmp_str,',',values);
455. // if(values.size()<2){
456. // cerr<<"Error:expected two values in -M parameter"<<endl;
457. // exit(1);
458. // }
459. // gp.adaMis=atoi(values[0].c_str());
460. // gp.adaMis2=atoi(values[1].c_str());
461. // }
462. // break;
463. // }
464. // case 'A':{
465. // wrong_paras["filtersRNA"].emplace_back("-A|adaMR");
466. // string tmp_str;
467. // tmp_str.assign(optarg);
468. // if(tmp_str.find(",")==string::npos){
469. // gp.adaMR=atof(optarg);
470. // gp.adaMR2=gp.adaMR;
471. // }else{
472. // vector<string> values;
473. // line_split(tmp_str,',',values);
474. // if(values.size()<2){
475. // cerr<<"Error:expected two values in -A parameter"<<endl;
476. // exit(1);
477. // }
478. // gp.adaMR=atof(values[0].c_str());
479. // gp.adaMR2=atof(values[1].c_str());
480. // }
481. // break;
482. // }
483. // case '9':{
484. // wrong_paras["filtersRNA"].emplace_back("-9|--adaEdge");
485. // string tmp_str;
486. // tmp_str.assign(optarg);
487. // if(tmp_str.find(",")==string::npos){
488. // gp.adaEdge=atoi(optarg);
489. // gp.adaEdge2=gp.adaEdge;
490. // }else{
491. // vector<string> values;
492. // line_split(tmp_str,',',values);
493. // if(values.size()<2){
494. // cerr<<"Error:expected two values in -9 parameter"<<endl;
495. // exit(1);
496. // }
497. // gp.adaEdge=atoi(values[0].c_str());
498. // gp.adaEdge2=atoi(values[1].c_str());
499. // }
500. // break;
501. // }
502. // case 'S':gp.adaRCtg=atoi(optarg);wrong_paras["filter"].emplace_back("-S|--adaRCtg");break;
503. // case 's':gp.adaRAr=atof(optarg);wrong_paras["filter"].emplace_back("-s|--adaRAr");break;
504. // case 'U':gp.adaRMa=atoi(optarg);wrong_paras["filter"].emplace_back("-U|--adaRMa");break;
505. // case 'u':gp.adaREr=atof(optarg);wrong_paras["filter"].emplace_back("-u|--adaREr");break;
506. // case 'b':gp.adaRMm=atoi(optarg);wrong_paras["filter"].emplace_back("-b|--adaRMm");break;
507. // // case 'd':gp.rmdup = true;break;
508. // case '0':gp.log.assign(optarg);break;
509. case 'v':
510. printVersion();
511. return 1;
512. case 'h':
513. printUsage(c_module);
514. return 1;
515. default:
516. {
517. exit(1);
518. }
519. }
520. }
521. if (argc != optind + 1)
522. {
523. cerr << "Error:please check the options" << endl;
524. exit(1);
525. }
526. if (!gp.rmdup || gp.cleanOutSplit <= 0)
527. {
528. }
529. else
530. {
531. cerr << "Warning:generating split files(-w was set) would become slower when rmdup function was on" << endl;
532. }
533. if (gp.log.find("/") == string::npos)
534. {
535. gp.log = gp.output_dir + "/" + gp.log;
536. }
537. if (gp.fq1_path.rfind(".gz") != gp.fq1_path.size() - 3)
538. {
539. gp.mode = "ssd";
540. }
541. if (gp.patchSize == 0)
542. {
543. gp.patchSize = gp.threads_num * 20000 / 8;
544. }
545. /*
546. int min_adapter_length=gp.adapter1_seq.size()>gp.adapter2_seq.size()?gp.adapter2_seq.size():gp.adapter1_seq.size();
547. if(min_adapter_length>20){
548. gp.adaEdge=min_adapter_length*ADA_RATIO;
549. }
550. */
551. return 0;
552. }
554. bool check_parameter(int argc, char *argv[], C_global_parameter &gp)
555. {
556. // 对于filterHts模块不再做检查
557. bool pe_data = false;
558. if (gp.module_name != "filterHts")
559. {
560. if (!gp.fq1_path.empty())
561. {
562. if (file_exist_and_not_empty(gp.fq1_path) == 0)
563. {
564. cerr << "Error:input fastq1 is required" << endl;
565. exit(1);
566. }
567. }
568. else
569. {
570. cerr << "Error:input fastq1 is required" << endl;
571. exit(1);
572. }
573. if (gp.output_dir.empty())
574. {
575. cerr << "Error:output directory is required" << endl;
576. exit(1);
577. }
578. if (!gp.fq2_path.empty())
579. {
580. pe_data = true;
581. if (file_exist_and_not_empty(gp.fq2_path) == 0)
582. {
583. cerr << "Error:input fastq2 is required" << endl;
584. exit(1);
585. }
586. if (gp.fq1_path == gp.fq2_path)
587. {
588. cerr << "Error:input fq1 and fq2 are the same,please check the parameters" << endl;
589. exit(1);
590. }
591. }
592. if (gp.clean_fq1.empty())
593. {
594. cerr << "Error:output clean fastq is required" << endl;
595. exit(1);
596. }
597. else
598. {
599. if (pe_data)
600. {
601. if (gp.clean_fq2.empty())
602. {
603. cerr << "Error:output clean fastq2 is required" << endl;
604. exit(1);
605. }
606. if (!(
607. gp.clean_fq1.rfind(".gz") == gp.clean_fq1.size() - 3 && gp.clean_fq2.rfind(".gz") == gp.clean_fq2.size() - 3) &&
608. !(
609. gp.clean_fq1.rfind(".gz") != gp.clean_fq1.size() - 3 && gp.clean_fq2.rfind(".gz") != gp.clean_fq2.size() - 3))
610. {
611. cerr << "Error:the format of clean fastq1 is inconsistent with fastq2" << endl;
612. exit(1);
613. }
614. if (gp.cleanOutSplit > 0 || gp.total_reads_num > 0)
615. {
616. if (gp.clean_fq1.rfind(".gz") != gp.clean_fq1.size() - 3 && gp.clean_fq2.rfind(".gz") != gp.clean_fq2.size() - 3)
617. {
618. cerr << "Error:the clean out fastq should be non-gz format when clean output reads are limited"
619. << endl;
620. exit(1);
621. }
622. }
623. }
624. }
625. if (!pe_data && gp.module_name != "filtersRNA")
626. {
627. if (!gp.adapter2_seq.empty())
628. {
629. cerr << "Error:no need adapter2" << endl;
630. exit(1);
631. }
632. }
633. if (!pe_data)
634. {
635. if (!gp.trim_fq2.empty() || !gp.clean_fq2.empty())
636. {
637. cerr << "Error:input file is not pe data" << endl;
638. exit(1);
639. }
640. }
641. else
642. {
643. if (!(
644. gp.fq1_path.rfind(".gz") == gp.fq1_path.size() - 3 && gp.fq2_path.rfind(".gz") == gp.fq2_path.size() - 3) &&
645. !(
646. gp.fq1_path.rfind(".gz") != gp.fq1_path.size() - 3 && gp.fq2_path.rfind(".gz") != gp.fq2_path.size() - 3))
647. {
648. cerr << "Error:the format of input fastq1 is inconsistent with fastq2" << endl;
649. exit(1);
650. }
651. }
653. if (gp.seq_type != "0" && gp.seq_type != "1")
654. {
655. cerr << "Error:seq_type value should be 0 or 1" << endl;
656. exit(1);
657. }
658. if (gp.output_file_type != "fastq" && gp.output_file_type != "fasta")
659. {
660. cerr << "Error:output_file_type value should be fastq or fasta" << endl;
661. exit(1);
662. }
663. }
664. if (gp.module_name == "filterStLFR")
665. {
666. // check barcode list file
667. if (gp.barcodeListPath == "")
668. {
669. cerr << "Error:barcode list not assigned" << endl;
670. exit(1);
671. }
672. else
673. {
674. ifstream bl(gp.barcodeListPath);
675. if (!bl)
676. {
677. cerr << "Error:cannot open such file," << gp.barcodeListPath << endl;
678. exit(1);
679. }
680. bl.close();
681. }
682. // check barcode region format
683. if (gp.barcodeRegionStr.find("_") == string::npos)
684. {
685. cerr << "Error:barcode region format error, it should be set as 101_10,117_10,133_10" << endl;
686. exit(1);
687. }
688. else
689. {
690. vector<string> eles;
691. line_split(gp.barcodeRegionStr, ',', eles);
692. if (eles.size() != 3)
693. {
694. cerr << "Error:barcode region format error, it should be set as 101_10,117_10,133_10" << endl;
695. exit(1);
696. }
697. for (vector<string>::iterator ix = eles.begin(); ix != eles.end(); ix++)
698. {
699. vector<string> eles2;
700. line_split(*ix, '_', eles2);
701. if (eles2.size() != 2)
702. {
703. cerr << "Error:barcode region format error, it should be set as 101_10,117_10,133_10" << endl;
704. exit(1);
705. }
706. }
707. }
708. }
709. if (!gp.adapter_discard_or_trim.empty())
710. {
711. if (gp.adapter_discard_or_trim != "trim" && gp.adapter_discard_or_trim != "discard")
712. {
713. cerr << "Error:adapter_discard_or_trim value should be trim or discard" << endl;
714. exit(1);
715. }
716. }
717. if (!gp.tile.empty())
718. {
719. if (gp.tile.find("-") != string::npos)
720. {
721. vector<string> tmp_eles;
722. line_split(gp.tile, '-', tmp_eles);
723. string start, end;
724. for (string::size_type ix = tmp_eles[0].size() - 1; ix >= 0; ix--)
725. {
726. if (!isalpha(tmp_eles[0][ix]))
727. {
728. start.insert(start.begin(), tmp_eles[0][ix]);
729. }
730. }
731. for (string::size_type ix = tmp_eles[1].size() - 1; ix >= 0; ix--)
732. {
733. if (!isalpha(tmp_eles[1][ix]))
734. {
735. end.insert(end.begin(), tmp_eles[1][ix]);
736. }
737. }
738. if (atoi(start.c_str()) > atoi(end.c_str()))
739. {
740. cerr << "Error:tile value format error" << endl;
741. exit(1);
742. }
743. }
744. for (string::size_type ix = 0; ix != gp.tile.size(); ix++)
745. {
746. if (!isalnum(gp.tile[ix]) && gp.tile[ix] != '-' && gp.tile[ix] != ',')
747. {
748. cerr << "Error:tile value format error" << endl;
749. exit(1);
750. }
751. }
752. }
753. if (gp.qualityPhred != 64 && gp.qualityPhred != 33)
754. {
755. cerr << "Error:qualityPhred value error" << endl;
756. exit(1);
757. }
758. if (gp.outputQualityPhred != 64 && gp.outputQualityPhred != 33)
759. {
760. cerr << "Error:outputQualityPhred value error" << endl;
761. exit(1);
762. }
763. if (gp.module_name == "filter" || gp.module_name == "filterMeta")
764. {
765. if (wrong_paras["filter"].size() > 0)
766. {
767. cerr << "Error:these parameters should not appear in the module," << join_vector(wrong_paras["filter"], ',')
768. << endl;
769. exit(1);
770. }
771. }
772. else
773. {
774. if (gp.module_name == "filtersRNA")
775. {
776. // adaMis,adaMR;
777. if (wrong_paras["filtersRNA"].size() > 0)
778. {
779. cerr << "Error:these parameters should not appear in the module,"
780. << join_vector(wrong_paras["filtersRNA"], ',') << endl;
781. exit(1);
782. }
783. }
784. }
785. if (gp.cleanOutSplit != 0 && gp.cleanOutSplit < gp.patchSize)
786. {
787. cerr << "Error: output reads in each clean fastq file(-w) should be more than patch size(-e)" << endl;
788. exit(1);
789. }
790. if (gp.module_name != "filterHts")
791. {
792. if (!gp.trim.empty())
793. {
794. vector<string> tmp_eles;
795. line_split(gp.trim, ',', tmp_eles);
796. if (pe_data)
797. {
798. if (tmp_eles.size() != 4)
799. {
800. cerr << "Error:trim value format error" << endl;
801. exit(1);
802. }
803. }
804. else
805. {
806. if (tmp_eles.size() != 2)
807. {
808. cerr << "Error:trim value format error" << endl;
809. exit(1);
810. }
811. }
812. for (int i = 0; i < gp.trim.size(); i++)
813. {
814. if (!isdigit(gp.trim[i]) && gp.trim[i] != ',')
815. {
816. cerr << "Error:trim value format error:" << gp.trim << endl;
817. cerr << "e.g.: -t 10 2 10 2" << endl;
818. exit(1);
819. }
820. }
821. }
822. if (!gp.trimBadHead.empty())
823. {
824. vector<string> tmp_eles;
825. line_split(gp.trimBadHead, ',', tmp_eles);
826. if (pe_data)
827. {
828. if (tmp_eles.size() != 2)
829. {
830. cerr << "Error:trimBadHead value format error" << endl;
831. exit(1);
832. }
833. }
834. else
835. {
836. if (tmp_eles.size() != 1)
837. {
838. cerr << "Error:trimBadHead value format error" << endl;
839. exit(1);
840. }
841. }
842. }
843. if (!gp.trimBadTail.empty())
844. {
845. vector<string> tmp_eles;
846. line_split(gp.trimBadTail, ',', tmp_eles);
847. if (pe_data)
848. {
849. if (tmp_eles.size() != 2)
850. {
851. cerr << "Error:trimBadTail value format error" << endl;
852. exit(1);
853. }
854. }
855. else
856. {
857. if (tmp_eles.size() != 1)
858. {
859. cerr << "Error:trimBadTail value format error" << endl;
860. exit(1);
861. }
862. }
863. }
864. }
865. if (!gp.base_convert.empty())
866. {
867. set<string> acgt_s;
868. string tmp_s = "ACGTacgt";
869. string tmp_single;
870. for (string::size_type ix = 0; ix != tmp_s.size(); ix++)
871. {
872. tmp_single.insert(tmp_single.end(), tmp_s[ix]);
873. acgt_s.insert(tmp_single);
874. tmp_single = "";
875. }
876. if (gp.base_convert.find("TO") == string::npos && gp.base_convert.find("2") == string::npos)
877. {
878. cerr << "Error:base_convert value format error" << endl;
879. exit(1);
880. }
881. else
882. {
883. string firstC = gp.base_convert.substr(0, 1);
884. string lastC = gp.base_convert.substr(gp.base_convert.size() - 1, 1);
885. if (acgt_s.find(firstC) == acgt_s.end() || acgt_s.find(lastC) == acgt_s.end())
886. {
887. cerr << "Error:base_convert value format error" << endl;
888. exit(1);
889. }
890. }
891. }
892. if (gp.cleanOutSplit > 0 && gp.total_reads_num > 0)
893. {
894. cerr << "Error:-w and -L cannot be both assigned" << endl;
895. exit(1);
896. }
897. if (gp.cleanOutSplit > 0)
898. {
899. gp.clean_file_reads = gp.cleanOutSplit;
900. }
901. else if (gp.l_total_reads_num > 0)
902. {
903. gp.clean_file_reads = gp.l_total_reads_num;
904. }
905. if (gp.threads_num > get_nprocs())
906. {
907. gp.threads_num = get_nprocs();
908. cerr << "Warning:threads number exceeds the system cpu number" << endl;
909. // exit(1);
910. }
911. if (gp.patchSize > 5000000)
912. {
913. cerr << "Error:patchSize cannot exceed 5M considering memory usage" << endl;
914. exit(1);
915. }
916. return 0;
917. }
919. void printModule()
920. {
921. cout << endl;
922. cout << "Program: SOAPnuke\n";
923. cout << "Version: " << PACKAGEVERSION << "." << MINORVERSION << endl;
924. cout << "Contact: GongChun<[email hidden]> ChenYuXin<[email hidden]>" << endl;
925. cout << "Command:\n";
926. cout << " filter preprocessing normal Fastq files\n";
927. #ifdef _PROCESSHTS
928. cout << " filterHts preprocessing BAM/CRAM files\n";
929. #endif
930. cout << " filterStLFR preprocessing stLFR Fastq files\n";
931. cout << " filtersRNA preprocessing sRNA Fastq files\n";
932. // cout << " filterDGE preprocessing DGE sequences\n"; //not include filterDGE in thie version
933. cout << " filterMeta preprocessing Meta Fastq files\n";
934. cout << endl;
935. exit(1);
936. }
938. void printHtsUsage()
939. {
940. cout << "Usage: "
941. << "filterHts"
942. << " [OPTION]... \n";
943. cout << "\ncommonly used parameters\n";
944. cout << "\t-E, --ref\t\tFILE\t\treference file(required when process cram format)\n";
945. cout << "\t-1, \t\t\tFILE\t\tinput bam/cram file(required)\n";
946. cout << "\t-2, \t\t\tFILE\t\toutput bam/cram file(required)\n";
947. cout << "\t-o, --outDir\tSTR\t\toutput directory, directory must exists\n";
948. cout << "\n";
949. cout << "\t-f, --adapter1\t\tSTR\t\tadapter sequence or list file of read1\n";
950. cout << "\t-r, --adapter2\t\tSTR\t\tadapter sequence or list file of read2 (if PE)\n";
951. cout << "\t-J, --ada_trim\t\t\t\ttrim read when find adapter[default:discard]\n";
952. cout << "\t-T, --thread\t\tINT\t\tthreads number used in process[6]\n";
953. cout
954. << "\t-c, --configFile\tSTR\t\tconfig file which include uncommonly used parameters. Each line contains a parameter, for value needed parameter: adaMis=2, for bool parameter:contam_trim, which means change mode from discard to trim\n";
955. cout << "\t-l, --lowQual\t\tINT\t\tlow quality threshold [default:5]\n";
956. cout << "\t-q, --qualRate\t\tFLOAT\t\tlow quality rate [default:0.5]\n";
957. cout << "\t-n, --nRate\t\tFLOAT\t\tN rate threshold [default:0.05]\n";
958. cout << "\t-m, --mean\t\tFLOAT\t\tfilter reads with low average quality\n";
959. cout << "\t-p, --highA\t\tFLOAT\t\tfilter reads if ratio of A in a read exceed [FLOAT]\n";
960. cout << "\t-g, --polyG_tail\tFLOAT\t\tfilter reads if found polyG in tail [INT]\n";
961. cout << "\t-X, --polyX\t\tINT\t\tfilter reads if a read contains polyX [INT]\n";
962. cout << "\t-i, --index\t\t\t\tremove index\n";
963. cout
964. << "\t-L, --totalReadsNum\tINT/FLOAT\tnumber/fraction of reads you want to keep in the output clean fq file(cannot be assigned when -w is given).\n";
965. cout
966. << "\t \t\t\t\t\tIt will extract reads randomly through the total clean fq file by default, you also can get the head reads\n";
967. cout << "\t \t\t\t\t\tfor save time by add head suffix to the integer(e.g. -L 10000000head)\n";
968. cout << "\t-4, --minReadLen\tINT\t\tread min length,default 18 for filtersRNA,30 for other modules\n";
969. cout << "\t-h, --help\t\t\t\thelp" << endl;
970. cout << "\t-v, --version\t\t\t\tshow version" << endl;
971. cout
972. << "\tExample: ./SOAPnuke filterHts -l 10 -q 0.1 -n 0.01 -f AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA -r AAGTCGGATCGTAGCCATGTCGTTCTGTGAGCCAAGGAGTTG --ref chr21.fa -1 input.bam -2 output.cram -o result"
973. << endl;
974. cout << "\n";
975. cout << "uncommonly used parameters, which should be set in config file" << endl;
976. // cout << "\tcontam1\t\tSTR\t\tcontaminant sequence(s) for read1, split by comma\n";
977. // cout << "\tcontam2\t\tSTR\t\tcontaminant sequence(s) for read2, split by comma\n";
978. cout << "\tctMatchR\tFLOAT/STR\tcontam's shortest consistent matching ratio [default:0.2]\n";
979. // cout << "\tglobal_contams\tSTR\t\tglobal contaminant sequences which need to be detected, split by comma if more than 1\n";
980. // cout << "\tglob_cotm_mR\tSTR\t\tminimum match ratio in global contaminant sequences detection\n";
981. // cout << "\tglob_cotm_mM\tSTR\t\tmaximum mismatch number in global contaminant sequences detection\n";
982. cout << "\tseqType\t\tINT\t\tSequence fq name type, 0->old fastq name, 1->new fastq name [0]\n";
983. cout << "\t\t\t\t\t\t\t\told fastq name: @FCD1PB1ACXX:4:1101:1799:2201#GAAGCACG/2\n";
984. cout << "\t\t\t\t\t\t\t\tnew fastq name: @HISEQ:310:C5MH9ANXX:1:1101:3517:2043 2:N:0:TCGGTCAC\n";
985. cout << "\ttile\t\tSTR\t\ttile number to ignore reads, such as [1101-1104,1205]\n";
986. cout << "\tfov\t\tSTR\t\tfov number to ignore reads (only for zebra-platform data), such as [C001R003,C003R004]\n";
987. cout << "\tAdapter related:\n";
989. // cout << "\tcontam_trim\t\t\ttrim read when find contam[default:discard]\n";
990. cout << "\tadaMis\t\tINT,[INT]\tthe max mismatch number when match the adapter (depend on -f/-r).If different values are required for fq1 and fq2, you can set two values seperated by comma,e.g. 1,2 (same as -A/-9)[1]\n";
991. cout << "\tadaMR\t\tFLOAT,[FLOAT]\tadapter's shortest match ratio (depend on -f/-r) [0.5]\n";
992. cout << "\tadaEdge\t\tINT,[INT]\tthe min length for segmental alignment [6]\n";
993. cout << endl;
994. cout << "\toverlap\t\tINT\t\tfilter the small insert size.Not filter until the value exceed 1[-1]\n";
995. cout << "\tmis\t\tFLOAT\t\tthe maximum mismatch ratio when find overlap between PE reads(depend on -O)[0.1]\n";
996. cout << "\n";
997. cout << "\tpatch\t\tINT\t\treads number of a patch processed[400000]\n";
998. cout << "\n";
999. cout << "\tqualSys\t\tINT\t\tquality system 1:64, 2:33[default:2]\n";
1000. cout << "\toutQualSys\tINT\t\tout quality system 1:64, 2:33[default:2]\n";
1001. cout << "\tmaxReadLen\tINT\t\tread max length,default 49 for filtersRNA\n";
1003. cout << "\n";
1005. cout << "\tpe_info\t\t\t\tAdd /1, /2 at the end of fastq name.[default:not add]\n";
1006. cout << "\tbaseConvert\tSTR\t\tconvert base when write data,example:TtoU, means convert base T to base U in the output\n";
1007. cout << "\tlog\t\tSTR\t\tlog file\n";
1008. cout << "\n";
1010. exit(1);
1011. }
1013. void printUsage(string c_module)
1014. {
1015. cout << "Usage: " << c_module << " [OPTION]... \n";
1016. cout << "\ncommonly used parameters\n";
1017. // cout << "\t-E, --mode\t\tSTR\t\tif pigz software is available,you can assign ssd mode for accelerate(-E ssd). Default is non SSD\n";
1018. cout << "\t-1, --fq1\t\tFILE\t\tfq1 file(required), .gz or normal text format are both supported(required)\n";
1019. cout << "\t-2, --fq2\t\tFILE\t\tfq2 file(used when process PE data), format should be same as fq1 file, both are gz or both are normal text\n";
1020. cout << "\t-C, --cleanFq1\t\tSTR\t\treads which passed QC from fq1 file would output to this file\n";
1021. cout << "\t-D, --cleanFq2\t\tSTR\t\treads which passed QC from fq2 file would output to this file\n";
1022. cout << "\t-o, --outDir\t\tSTR\t\tOutput directory. Processed fq files and statistical results would be output to here\n";
1023. cout << "\t-f, --adapter1\t\tSTR\t\tadapter sequence or list file of read1\n";
1024. cout << "\t-r, --adapter2\t\tSTR\t\tadapter sequence or list file of read2 (if PE)\n";
1025. cout << "\t-J, --ada_trim\t\t\t\ttrim read when find adapter[default:discard]\n";
1026. cout << "\t-T, --thread\t\tINT\t\tthreads number used in process[6]\n";
1027. cout << "\t-c, --configFile\tSTR\t\tconfig file which include uncommonly used parameters. Each line contains a parameter, for value needed parameter: adaMis=2, for bool parameter:contam_trim, which means change mode from discard to trim\n";
1028. cout << "\t-l, --lowQual\t\tINT\t\tlow quality threshold [default:5]\n";
1029. cout << "\t-q, --qualRate\t\tFLOAT\t\tlow quality rate [default:0.5]\n";
1030. cout << "\t-n, --nRate\t\tFLOAT\t\tN rate threshold [default:0.05]\n";
1031. cout << "\t-m, --mean\t\tFLOAT\t\tfilter reads with low average quality\n";
1032. cout << "\t-p, --highA\t\tFLOAT\t\tfilter reads if ratio of A in a read exceed [FLOAT]\n";
1033. cout << "\t-g, --polyG_tail\tFLOAT\t\tfilter reads if found polyG in tail [INT]\n";
1034. cout << "\t-X, --polyX\t\tINT\t\tfilter reads if a read contains polyX [INT]\n";
1035. cout << "\t-4, --minReadLen\tINT\t\tread min length,default 18 for filtersRNA,30 for other modules\n";
1036. cout << "\t-w, --cleanOutSplit\tINT\t\tmax reads number in each output clean fastq file\n";
1037. cout << "\t-h, --help\t\t\t\thelp" << endl;
1038. cout << "\t-v, --version\t\t\t\tshow version" << endl;
1039. if (c_module == "filter")
1040. {
1041. cout
1042. << "\tExample: ./SOAPnuke filter -l 10 -q 0.1 -n 0.01 -f AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA -r AAGTCGGATCGTAGCCATGTCGTTCTGTGAGCCAAGGAGTTG -1 test.r1.fq.gz -2 test.r2.fq.gz -C clean_1.fq.gz -D clean_2.fq.gz -o result -T 8"
1043. << endl;
1044. }
1045. cout << "\n";
1046. cout << "uncommonly used parameters, which should be set in config file" << endl;
1047. // cout << "\tcontam1\t\tSTR\t\tcontaminant sequence(s) for fq1 file, split by comma\n";
1048. // cout << "\tcontam2\t\tSTR\t\tcontaminant sequence(s) for fq2 file, split by comma. The two parameters is similar to -f/-r, but support a list of sequences\n";
1049. cout << "\tAdapter related:\n";
1050. cout << "\tctMatchR\tFLOAT/STR\tcontam's shortest consistent matching ratio [default:0.2]\n";
1051. // cout << "\tglobal_contams\tSTR\t\tglobal contaminant sequences which need to be detected, split by comma if more than 1\n";
1052. // cout << "\tglob_cotm_mR\tSTR\t\tminimum match ratio in global contaminant sequences detection\n";
1053. // cout << "\tglob_cotm_mM\tSTR\t\tmaximum mismatch number in global contaminant sequences detection\n";
1054. cout << "\tseqType\t\tINT\t\tSequence fq name type, 0->old fastq name, 1->new fastq name [0]\n";
1055. cout << "\t\t\t\t\t\t\t\told fastq name: @FCD1PB1ACXX:4:1101:1799:2201#GAAGCACG/2\n";
1056. cout << "\t\t\t\t\t\t\t\tnew fastq name: @HISEQ:310:C5MH9ANXX:1:1101:3517:2043 2:N:0:TCGGTCAC\n";
1057. cout << "\ttrimFq1\t\tSTR\t\ttrim fq1 file name(gz format)[optional]\n";
1058. cout
1059. << "\ttrimFq2\t\tSTR\t\ttrim fq2 file name[optional]. If trim related parameters were set on, these output files would include the total reads which only do trimming. For example, if read A failed QC after trimming, it will still output to -R/-W, but not to -C/-D\n";
1060. cout << "\ttile\t\tSTR\t\ttile number to ignore reads, such as [1101-1104,1205]\n";
1061. cout << "\tfov\t\tSTR\t\tfov number to ignore reads (only for zebra-platform data), such as [C001R003,C003R004]\n";
1062. // cout << "\tcontam_trim\t\t\ttrim read when find contam[default:discard]\n";
1063. if (c_module == "filter" || c_module == "filtermeta")
1064. {
1065. // cout << "filter and filtermeta module adapter related parameter\n";
1066. cout
1067. << "\tadaMis\t\tINT,[INT]\tthe max mismatch number when match the adapter (depend on -f/-r).If different values are required for fq1 and fq2, you can set two values seperated by comma,e.g. 1,2 (same as -A/-9)[1]\n";
1068. cout << "\tadaMR\t\tFLOAT,[FLOAT]\tadapter's shortest match ratio (depend on -f/-r) [0.5]\n";
1069. cout << "\tadaEdge\t\tINT,[INT]\tthe min length for segmental alignment [6]\n";
1070. }
1071. else if (c_module == "filtersRNA")
1072. {
1073. // cout << "filtersRNA module adapter related parameter\n";
1074. cout << " find 5' adapter\n";
1075. cout << "\tadaRCtg\tINT\t\tmini 5' adapter continuous alignment length (default: 6)\n";
1076. cout << "\tadaRAr\tFLOAT\t\tmini alignment rate when find 5' adapter: alignment/tag (default: 0.8)\n";
1077. cout << " find 3' adapter\n";
1078. cout << "\tadaRMa\tINT\t\tmini alignment length when find 3' adapter (default: 5)\n";
1079. cout << "\tadaREr\tFLOAT\t\tMax error rate when find 3' adapter (mismatch/match) (default: 0.4)\n";
1080. cout << "\tadaRMm\tINT\t\tMax mismatch number when find 3' adapter (default: 4)" << endl;
1081. }
1082. else if (c_module == "filterStLFR")
1083. {
1084. /*if(para=="barcodeListPath"){
1085. gp.barcodeListPath.assign(value.c_str());
1086. }
1087. if(para=="barcodeRegionStr"){
1088. gp.barcodeRegionStr.assign(value.c_str());
1089. }
1090. if(para=="notCutNoLFR"){
1091. gp.notCutNoLFR=true;
1092. }*/
1093. cout << "\tbarcodeListPath\tSTR\t\tbarcode list of two columns:sequence and barcodeID\n";
1094. cout
1095. << "\tbarcodeRegionStr\tSTR\t\tbarcode regions, such as: 101_10,117_10,145_10 or 101_10,117_10,133_10\n";
1096. cout << "\tnotCutNoLFR\t\t\t\tdo not cut sequence when fail found barcode\n";
1097. cout << "\tinputAsList\t\t\t\tinput file list not a file\n";
1098. cout << "\ttenX\t\t\t\t\toutput tenX format\n";
1099. }
1100. cout << "\toutFileType\tSTR\t\toutput file format: fastq or fasta[fastq]\n";
1102. cout << endl;
1103. if (c_module == "filterStLFR" || c_module == "filter")
1104. {
1105. cout
1106. << "\trmdup\t\t\tremove duplicate reads. The function contains a certain false positive rate which means a small number of non-duplicate reads would be marked as duplication, and limited in large reads number\n";
1107. // cout<<"\tapproximateReadsNum\tSTR/LONG\tapproximate reads number. We suggest you set this parameter if you know. e.g. 100m"<<endl;
1108. // cout<<"\tmemSizeUsedInRmdup\tSTR/LONG\tmaximum memory size used in rmdup. e.g. 1g\n";
1109. // cout<<"\texpectedFalsePositive\tFLOAT\texpected false positive rate which means maybe a <expectedFalsePositive> of reads would be marked as duplication which actually not. e.g. 1e-7"<<endl;
1110. }
1111. cout << "\tindex\t\t\t\tremove index\n";
1112. cout
1113. << "\ttotalReadsNum\tINT/FLOAT\tnumber/fraction of reads you want to keep in the output clean fq file(cannot be assigned when -w is given).\n";
1114. cout
1115. << "\t \t\t\t\t\tIt will extract reads randomly through the total clean fq file by default, you also can get the head reads\n";
1116. cout << "\t \t\t\t\t\tfor save time by add head suffix to the integer(e.g. totalReadsNum=10000000head)\n";
1117. cout << "\ttrim\t\tINT,INT,INT,INT\n";
1118. cout
1119. << "\t \t\t\t\t\ttrim some bp of the read's head and tail, they means: (PE type:read1's head and tail and read2's head and tail [0,0,0,0]; SE type:read head and tail [0,0])\n";
1121. cout
1122. << "\ttrimBadHead\tINT,INT\t\tTrim from head ends until meeting high-quality base or reach the length threshold, set (quality threshold,MaxLengthForTrim) [0,0]\n";
1123. cout
1124. << "\ttrimBadTail\tINT,INT\t\tTrim from tail ends until meeting high-quality base or reach the length threshold, set (quality threshold,MaxLengthForTrim) [0,0]\n";
1125. cout << "\n";
1127. cout << "\toverlap\t\tINT\t\tfilter the small insert size.Not filter until the value exceed 1[-1]\n";
1128. cout << "\tmis\t\tFLOAT\t\tthe maximum mismatch ratio when find overlap between PE reads(depend on -O)[0.1]\n";
1129. cout << "\n";
1130. // cout << "\tsplit_line\tINT\t\tsplit raw fastq by <split_line>, default 10M reads per file (if ssd mode is open)\n";
1131. cout << "\tpatch\t\tINT\t\treads number of a patch processed[400000]\n";
1132. cout << "\n";
1133. cout << "\tqualSys\t\tINT\t\tquality system 1:64, 2:33[default:2]\n";
1134. cout << "\toutQualSys\tINT\t\tout quality system 1:64, 2:33[default:2]\n";
1135. cout << "\tmaxReadLen\tINT\t\tread max length,default 49 for filtersRNA\n";
1137. // cout << "\tcleanOutSplit\tINT\t\tmax reads number in each output clean fastq file\n";
1138. // cout << "\tappend STR the log's output place : console or file [console]\n";
1140. cout << "\n";
1142. cout << "\tpe_info\t\t\t\tAdd /1, /2 at the end of fastq name.[default:not add]\n";
1143. cout
1144. << "\tbaseConvert\tSTR\t\tconvert base when write data,example:TtoU, means convert base T to base U in the output\n";
1145. cout << "\tmaxBaseQuality\tINT\t\tmaximum base quality[42]\n";
1146. cout << "\tlog\t\tSTR\t\tlog file\n";
1147. cout << "\n";
1149. exit(1);
1150. }
1152. void printVersion()
1153. {
1154. cerr << "SOAPnuke filter tools version " << PACKAGEVERSION << "." << MINORVERSION << "\n";
1155. exit(1);
1156. }
1158. void initFromConfigFile(C_global_parameter &gp, char *configFile)
1159. {
1160. set<string> legalParas, boolParas;
1161. legalParas.insert("trimFq1");
1162. legalParas.insert("trimFq2");
1163. legalParas.insert("seqType");
1164. legalParas.insert("outFileType");
1165. legalParas.insert("contam_trim");
1166. legalParas.insert("contam1");
1167. legalParas.insert("contam2");
1168. legalParas.insert("ctMatchR");
1169. legalParas.insert("global_contams");
1170. legalParas.insert("glob_cotm_mR");
1171. legalParas.insert("glob_cotm_mM");
1172. legalParas.insert("tile");
1173. legalParas.insert("fov");
1174. legalParas.insert("index");
1175. legalParas.insert("qualSys");
1176. legalParas.insert("outQualSys");
1177. legalParas.insert("baseConvert");
1178. legalParas.insert("maxBaseQuality");
1179. legalParas.insert("overlap");
1180. legalParas.insert("mis");
1181. legalParas.insert("pe_info");
1182. legalParas.insert("patch");
1183. legalParas.insert("maxReadLen");
1184. legalParas.insert("adaMis");
1185. legalParas.insert("adaMR");
1186. legalParas.insert("adaEdge");
1187. legalParas.insert("adaRCtg");
1188. legalParas.insert("adaRAr");
1189. legalParas.insert("adaRMa");
1190. legalParas.insert("adaREr");
1191. legalParas.insert("adaRMm");
1192. legalParas.insert("log");
1193. legalParas.insert("totalReadsNum");
1194. legalParas.insert("cleanOutSplit");
1195. legalParas.insert("trim");
1196. legalParas.insert("trimBadHead");
1197. legalParas.insert("trimBadTail");
1198. // stLFR
1199. legalParas.insert("barcodeListPath");
1200. legalParas.insert("barcodeRegionStr");
1201. legalParas.insert("notCutNoLFR");
1202. legalParas.insert("inputAsList");
1203. legalParas.insert("tenX");
1204. legalParas.insert("rmdup");
1205. // legalParas.insert("approximateReadsNum");
1206. // legalParas.insert("memSizeUsedInRmdup");
1207. // legalParas.insert("expectedFalsePositive");
1209. boolParas.insert("index");
1210. boolParas.insert("pe_info");
1211. boolParas.insert("contam_trim");
1212. // stLFR
1213. boolParas.insert("notCutNoLFR");
1214. boolParas.insert("inputAsList");
1215. boolParas.insert("tenX");
1217. boolParas.insert("rmdup");
1218. string configFilePath(configFile);
1219. ifstream ifConfig;
1220. ifConfig.open(configFilePath.c_str());
1221. if (!ifConfig)
1222. {
1223. cerr << "Error:cannot open such file," << configFile << endl;
1224. exit(1);
1225. }
1226. string line;
1227. while (getline(ifConfig, line))
1228. {
1229. string para, value;
1230. if (line.find("#") == 0)
1231. {
1232. continue;
1233. }
1234. if (line.find("=") != string::npos)
1235. {
1236. vector<string> eles;
1237. line_split(line, '=', eles);
1238. if (eles.size() != 2)
1239. {
1240. cerr << "Error:unrecgonized format parameter," << line << endl;
1241. exit(1);
1242. }
1243. para = eles[0];
1244. value = eles[1];
1245. chomp_space(para, "all");
1246. chomp_space(value, "all");
1247. }
1248. else
1249. {
1250. para = line;
1251. if (boolParas.find(para) == boolParas.end())
1252. {
1253. cerr << "Error:this parameter should set a value," << para << endl;
1254. exit(1);
1255. }
1256. }
1257. if (legalParas.find(para) == legalParas.end())
1258. {
1259. cerr << "Error:no such parameter," << para << endl;
1260. exit(1);
1261. }
1262. if (para == "trimFq1")
1263. {
1264. gp.trim_fq1.assign(value);
1265. if (gp.trim_fq1.rfind(".gz") == gp.trim_fq1.size() - 3)
1266. {
1267. gp.trimOutGzformat = true;
1268. }
1269. else
1270. {
1271. gp.trimOutGzformat = false;
1272. }
1273. }
1274. if (para == "trimFq2")
1275. {
1276. gp.trim_fq2.assign(value);
1277. }
1278. if (para == "seqType")
1279. {
1280. gp.seq_type.assign(value);
1281. }
1282. if (para == "outFileType")
1283. {
1284. gp.output_file_type.assign(value);
1285. }
1286. if (para == "contam_trim")
1287. {
1288. gp.contam_discard_or_trim = "trim";
1289. }
1290. if (para == "contam1")
1291. {
1292. gp.contam1_seq.assign(value);
1293. }
1294. if (para == "contam2")
1295. {
1296. gp.contam2_seq.assign(value);
1297. }
1298. if (para == "ctMatchR")
1299. {
1300. gp.ctMatchR.assign(value.c_str());
1301. }
1302. if (para == "global_contams")
1303. {
1304. gp.global_contams.assign(value.c_str());
1305. }
1306. if (para == "glob_cotm_mR")
1307. {
1308. gp.g_mrs.assign(value.c_str());
1309. }
1310. if (para == "glob_cotm_mM")
1311. {
1312. gp.g_mms.assign(value.c_str());
1313. }
1314. if (para == "tile")
1315. {
1316. gp.tile.assign(value.c_str());
1317. }
1318. if (para == "fov")
1319. {
1320. gp.fov.assign(value.c_str());
1321. }
1322. if (para == "index")
1323. {
1324. gp.index_remove = true;
1325. }
1326. if (para == "qualSys")
1327. {
1328. gp.qualityPhred = atoi(value.c_str());
1329. if (gp.qualityPhred == 1)
1330. {
1331. gp.qualityPhred = 64;
1332. }
1333. else if (gp.qualityPhred == 2)
1334. {
1335. gp.qualityPhred = 33;
1336. }
1337. }
1338. if (para == "outQualSys")
1339. {
1340. gp.outputQualityPhred = atoi(value.c_str());
1341. if (gp.outputQualityPhred == 1)
1342. {
1343. gp.outputQualityPhred = 64;
1344. }
1345. else if (gp.outputQualityPhred == 2)
1346. {
1347. gp.outputQualityPhred = 33;
1348. }
1349. }
1350. if (para == "baseConvert")
1351. {
1352. gp.base_convert.assign(value.c_str());
1353. }
1354. if (para == "maxBaseQuality")
1355. {
1356. gp.maxBaseQuality = atoi(value.c_str());
1357. }
1358. if (para == "overlap")
1359. {
1360. gp.overlap_length = atoi(value.c_str());
1361. }
1362. if (para == "mis")
1363. {
1364. gp.peMismatchRatio = atof(value.c_str());
1365. }
1366. if (para == "pe_info")
1367. {
1368. gp.whether_add_pe_info = true;
1369. }
1370. if (para == "patch")
1371. {
1372. gp.patchSize = atoi(value.c_str());
1373. }
1374. if (para == "maxReadLen")
1375. {
1376. gp.max_read_length = atoi(value.c_str());
1377. }
1378. if (para == "adaMis")
1379. {
1380. wrong_paras["filtersRNA"].emplace_back("-M|--adaMis");
1381. string tmp_str;
1382. tmp_str.assign(value.c_str());
1383. if (tmp_str.find(",") == string::npos)
1384. {
1385. gp.adaMis = atoi(value.c_str());
1386. gp.adaMis2 = gp.adaMis;
1387. }
1388. else
1389. {
1390. vector<string> values;
1391. line_split(tmp_str, ',', values);
1392. if (values.size() < 2)
1393. {
1394. cerr << "Error:expected two values in -M parameter" << endl;
1395. exit(1);
1396. }
1397. gp.adaMis = atoi(values[0].c_str());
1398. gp.adaMis2 = atoi(values[1].c_str());
1399. }
1400. }
1401. if (para == "adaMR")
1402. {
1403. wrong_paras["filtersRNA"].emplace_back("-A|adaMR");
1404. string tmp_str;
1405. tmp_str.assign(value.c_str());
1406. if (tmp_str.find(",") == string::npos)
1407. {
1408. gp.adaMR = atof(value.c_str());
1409. gp.adaMR2 = gp.adaMR;
1410. }
1411. else
1412. {
1413. vector<string> values;
1414. line_split(tmp_str, ',', values);
1415. if (values.size() < 2)
1416. {
1417. cerr << "Error:expected two values in -A parameter" << endl;
1418. exit(1);
1419. }
1420. gp.adaMR = atof(values[0].c_str());
1421. gp.adaMR2 = atof(values[1].c_str());
1422. }
1423. }
1424. if (para == "adaEdge")
1425. {
1426. wrong_paras["filtersRNA"].emplace_back("-9|--adaEdge");
1427. string tmp_str;
1428. tmp_str.assign(value.c_str());
1429. if (tmp_str.find(",") == string::npos)
1430. {
1431. gp.adaEdge = atoi(value.c_str());
1432. gp.adaEdge2 = gp.adaEdge;
1433. }
1434. else
1435. {
1436. vector<string> values;
1437. line_split(tmp_str, ',', values);
1438. if (values.size() < 2)
1439. {
1440. cerr << "Error:expected two values in -9 parameter" << endl;
1441. exit(1);
1442. }
1443. gp.adaEdge = atoi(values[0].c_str());
1444. gp.adaEdge2 = atoi(values[1].c_str());
1445. }
1446. }
1447. if (para == "adaRCtg")
1448. {
1449. gp.adaRCtg = atoi(value.c_str());
1450. wrong_paras["filter"].emplace_back("-S|--adaRCtg");
1451. }
1452. if (para == "adaRAr")
1453. {
1454. gp.adaRAr = atof(value.c_str());
1455. wrong_paras["filter"].emplace_back("-s|--adaRAr");
1456. }
1457. if (para == "adaRMa")
1458. {
1459. gp.adaRMa = atoi(value.c_str());
1460. wrong_paras["filter"].emplace_back("-U|--adaRMa");
1461. }
1462. if (para == "adaREr")
1463. {
1464. gp.adaREr = atof(value.c_str());
1465. wrong_paras["filter"].emplace_back("-u|--adaREr");
1466. }
1467. if (para == "adaRMm")
1468. {
1469. gp.adaRMm = atoi(value.c_str());
1470. wrong_paras["filter"].emplace_back("-b|--adaRMm");
1471. }
1472. if (para == "log")
1473. {
1474. gp.log.assign(value.c_str());
1475. }
1476. if (para == "totalReadsNum")
1477. {
1478. string tmp_str;
1479. tmp_str.assign(value.c_str());
1480. if (tmp_str.find("head") == string::npos)
1481. {
1482. gp.total_reads_num_random = true;
1483. // gp.catWhenrunning=false;
1484. for (int i = 0; i != tmp_str.size(); i++)
1485. {
1486. if (!isdigit(tmp_str[i]) && tmp_str[i] != '.')
1487. {
1488. cerr << "Error:-L value should be a positive integer or float" << endl;
1489. exit(1);
1490. }
1491. }
1492. }
1493. else
1494. {
1495. gp.total_reads_num_random = false;
1496. tmp_str.erase(tmp_str.find("head"), 4);
1497. if (tmp_str.find(".") != string::npos)
1498. {
1499. cerr << "Error:-L value should be a integer when with head suffix" << endl;
1500. exit(1);
1501. }
1502. else
1503. {
1504. for (int i = 0; i != tmp_str.size(); i++)
1505. {
1506. if (!isdigit(tmp_str[i]))
1507. {
1508. cerr << "Error:-L value should be an integer when with head suffix" << endl;
1509. exit(1);
1510. }
1511. }
1512. }
1513. }
1515. float tmp_val = atof(value.c_str());
1516. if (tmp_val == 0)
1517. {
1518. cerr << "Error:-L value should be a positive integer or float" << endl;
1519. exit(1);
1520. }
1521. gp.total_reads_num = tmp_val;
1522. if (tmp_val < 1)
1523. {
1524. gp.f_total_reads_ratio = tmp_val;
1525. }
1526. else
1527. {
1528. istringstream is_str(tmp_str);
1529. is_str >> gp.l_total_reads_num;
1530. }
1531. if (gp.f_total_reads_ratio > 0 && gp.l_total_reads_num > 0)
1532. {
1533. cerr << "Error:reads number and ratio should not be both assigned at the same time" << endl;
1534. exit(1);
1535. }
1536. }
1537. if (para == "cleanOutSplit")
1538. {
1539. string paraCheck = value;
1540. for (int i = 0; i != paraCheck.size(); i++)
1541. {
1542. if (!isdigit(paraCheck[i]))
1543. {
1544. cerr << "Error:-w value should be a positive integer" << endl;
1545. exit(1);
1546. }
1547. }
1548. gp.cleanOutSplit = atoi(paraCheck.c_str());
1549. if (gp.cleanOutSplit == 0)
1550. {
1551. cerr << "Error:-w value should be a positive integer" << endl;
1552. exit(1);
1553. }
1554. }
1555. if (para == "trim")
1556. {
1557. gp.trim = value;
1558. }
1559. if (para == "trimBadHead")
1560. {
1561. gp.trimBadHead = value;
1562. }
1563. if (para == "trimBadTail")
1564. {
1565. gp.trimBadTail = value;
1566. }
1567. /*legalParas.insert("cleanOutSplit");
1568. legalParas.insert("trim");
1569. legalParas.insert("trimBadHead");
1570. legalParas.insert("trimBadTail");
1571. */
1572. // stLFR
1573. /*string barcodeListPath;
1574. string barcodeRegionStr;
1575. bool notCutNoLFR;
1576. */
1577. if (para == "barcodeListPath")
1578. {
1579. gp.barcodeListPath.assign(value.c_str());
1580. }
1581. if (para == "barcodeRegionStr")
1582. {
1583. gp.barcodeRegionStr.assign(value.c_str());
1584. }
1585. if (para == "notCutNoLFR")
1586. {
1587. gp.notCutNoLFR = true;
1588. }
1589. if (para == "inputAsList")
1590. {
1591. gp.inputAsList = true;
1592. }
1593. if (para == "tenX")
1594. {
1595. gp.tenX = true;
1596. }
1598. if (para == "rmdup")
1599. {
1600. gp.rmdup = true;
1601. }
1603. // if(para=="approximateReadsNum"){
1604. // if(isdigit(value[value.size()-1])){
1605. // gp.approximateReadsNum=atol(value.c_str());
1606. // }else if(value[value.size()-1]=='M' || value[value.size()-1]=='m'){
1607. // string realValue=value.substr(0,value.size()-1);
1608. // gp.approximateReadsNum=atol(realValue.c_str())*1024*1024;
1609. // }else if(value[value.size()-1]=='G' || value[value.size()-1]=='g') {
1610. // string realValue = value.substr(0, value.size() - 1);
1611. // gp.approximateReadsNum = atol(realValue.c_str()) * 1024 * 1024 * 1024;
1612. // }else{
1613. // cerr<<"Error:cannot recognize value,"<<value<<endl;
1614. // exit(1);
1615. // }
1616. // }
1617. // if(para=="memSizeUsedInRmdup"){
1618. // if(isdigit(value[value.size()-1])){
1619. // gp.memSizeUsedInRmdup=atol(value.c_str());
1620. // }else if(value[value.size()-1]=='M' || value[value.size()-1]=='m'){
1621. // string realValue=value.substr(0,value.size()-1);
1622. // gp.memSizeUsedInRmdup=atol(realValue.c_str())*1024*1024;
1623. // }else if(value[value.size()-1]=='G' || value[value.size()-1]=='g') {
1624. // string realValue = value.substr(0, value.size() - 1);
1625. // gp.memSizeUsedInRmdup = atol(realValue.c_str()) * 1024 * 1024 * 1024;
1626. // }else{
1627. // cerr<<"Error:cannot recognize value,"<<value<<endl;
1628. // exit(1);
1629. // }
1630. // }
1631. // if(para=="expectedFalsePositive"){
1632. // gp.expectedFalsePositive=stof(value.c_str());
1633. // if(gp.expectedFalsePositive==0 || gp.expectedFalsePositive==1){
1634. // cout<<"Warning:expected FP you set is "<<gp.expectedFalsePositive<<", please check the parameter"<<endl;
1635. // }
1636. // }
1637. }
1638. }

process_argv.cpp at commit 2d5b727, under GPL-3.0 · at the source