3 from datetime import datetime, date
9 from gaworkflow.pipeline import firecrest
10 from gaworkflow.pipeline import bustard
11 from gaworkflow.pipeline import gerald
12 from gaworkflow.pipeline import runfolder
13 from gaworkflow.pipeline.runfolder import ElementTree
16 def make_flowcell_id(runfolder_dir, flowcell_id=None):
17 if flowcell_id is None:
18 flowcell_id = '207BTAAXY'
20 config = """<?xml version="1.0"?>
23 </FlowcellId>""" % (flowcell_id,)
24 config_dir = os.path.join(runfolder_dir, 'Config')
26 if not os.path.exists(config_dir):
28 pathname = os.path.join(config_dir, 'FlowcellId.xml')
29 f = open(pathname,'w')
33 def make_matrix(matrix_dir):
34 contents = """# Auto-generated frequency response matrix
41 -0.10 -0.10 1.17 -0.03
44 s_matrix = os.path.join(matrix_dir, 's_matrix.txt')
45 f = open(s_matrix, 'w')
49 def make_phasing_params(bustard_dir):
50 for lane in range(1,9):
51 pathname = os.path.join(bustard_dir, 'params%d.xml' % (lane))
52 f = open(pathname, 'w')
53 f.write("""<Parameters>
54 <Phasing>0.009900</Phasing>
55 <Prephasing>0.003500</Prephasing>
60 def make_gerald_config(gerald_dir):
61 config_xml = """<RunParameters>
62 <ChipWideRunParameters>
63 <ANALYSIS>default</ANALYSIS>
64 <BAD_LANES></BAD_LANES>
65 <BAD_TILES></BAD_TILES>
66 <CONTAM_DIR></CONTAM_DIR>
67 <CONTAM_FILE></CONTAM_FILE>
68 <ELAND_GENOME>Need_to_specify_ELAND_genome_directory</ELAND_GENOME>
69 <ELAND_MULTIPLE_INSTANCES>8</ELAND_MULTIPLE_INSTANCES>
70 <ELAND_REPEAT></ELAND_REPEAT>
71 <EMAIL_DOMAIN>domain.com</EMAIL_DOMAIN>
72 <EMAIL_LIST>diane</EMAIL_LIST>
73 <EMAIL_SERVER>localhost:25</EMAIL_SERVER>
74 <EXPT_DIR>/home/diane/gec/080416_HWI-EAS229_0024_207BTAAXX/Data/C1-33_Firecrest1.8.28_19-04-2008_diane/Bustard1.8.28_19-04-2008_diane</EXPT_DIR>
75 <EXPT_DIR_ROOT>/home/diane/gec</EXPT_DIR_ROOT>
77 <GENOME_DIR>/home/diane/proj/SolexaPipeline-0.2.2.6/Goat/../Gerald/../../Genomes</GENOME_DIR>
78 <GENOME_FILE>Need_to_specify_genome_file_name</GENOME_FILE>
79 <HAMSTER_FLAG>genome</HAMSTER_FLAG>
80 <OUT_DIR>/home/diane/gec/080416_HWI-EAS229_0024_207BTAAXX/Data/C1-33_Firecrest1.8.28_19-04-2008_diane/Bustard1.8.28_19-04-2008_diane/GERALD_19-04-2008_diane</OUT_DIR>
81 <POST_RUN_COMMAND></POST_RUN_COMMAND>
82 <PRB_FILE_SUFFIX>_prb.txt</PRB_FILE_SUFFIX>
83 <PURE_BASES>12</PURE_BASES>
84 <QF_PARAMS>'((CHASTITY>=0.6))'</QF_PARAMS>
85 <QHG_FILE_SUFFIX>_qhg.txt</QHG_FILE_SUFFIX>
86 <QUALITY_FORMAT>--symbolic</QUALITY_FORMAT>
87 <READ_LENGTH>32</READ_LENGTH>
88 <SEQUENCE_FORMAT>--scarf</SEQUENCE_FORMAT>
89 <SEQ_FILE_SUFFIX>_seq.txt</SEQ_FILE_SUFFIX>
90 <SIG_FILE_SUFFIX_DEPHASED>_sig2.txt</SIG_FILE_SUFFIX_DEPHASED>
91 <SIG_FILE_SUFFIX_NOT_DEPHASED>_sig.txt</SIG_FILE_SUFFIX_NOT_DEPHASED>
92 <SOFTWARE_VERSION>@(#) Id: GERALD.pl,v 1.68.2.2 2007/06/13 11:08:49 km Exp</SOFTWARE_VERSION>
93 <TILE_REGEX>s_[1-8]_[0-9][0-9][0-9][0-9]</TILE_REGEX>
94 <TILE_ROOT>s</TILE_ROOT>
95 <TIME_STAMP>Sat Apr 19 19:08:30 2008</TIME_STAMP>
96 <TOOLS_DIR>/home/diane/proj/SolexaPipeline-0.2.2.6/Goat/../Gerald</TOOLS_DIR>
97 <USE_BASES>all</USE_BASES>
98 <WEB_DIR_ROOT>http://host.domain.com/yourshare/</WEB_DIR_ROOT>
99 </ChipWideRunParameters>
100 <LaneSpecificRunParameters>
113 <s_2>/g/equcab1</s_2>
114 <s_3>/g/equcab1</s_3>
115 <s_4>/g/canfam2</s_4>
132 <s_1>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_1>
133 <s_2>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_2>
134 <s_3>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_3>
135 <s_4>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_4>
136 <s_5>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_5>
137 <s_6>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_6>
138 <s_7>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_7>
139 <s_8>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_8>
141 </LaneSpecificRunParameters>
144 pathname = os.path.join(gerald_dir, 'config.xml')
145 f = open(pathname,'w')
150 def make_summary_htm(gerald_dir):
151 summary_htm = """<!--RUN_TIME Mon Apr 21 11:52:25 2008 -->
152 <!--SOFTWARE_VERSION @(#) $Id: jerboa.pl,v 1.31 2007/03/05 17:52:15 km Exp $-->
156 <a name="Top"><h2><title>080416_HWI-EAS229_0024_207BTAAXX Summary</title></h2></a>
157 <h1>Summary Information For Experiment 080416_HWI-EAS229_0024_207BTAAXX on Machine HWI-EAS229</h1>
158 <h2><br></br>Chip Summary<br></br></h2>
159 <table border="1" cellpadding="5">
160 <tr><td>Machine</td><td>HWI-EAS229</td></tr>
161 <tr><td>Run Folder</td><td>080416_HWI-EAS229_0024_207BTAAXX</td></tr>
162 <tr><td>Chip ID</td><td>unknown</td></tr>
164 <h2><br></br>Lane Parameter Summary<br></br></h2>
165 <table border="1" cellpadding="5">
169 <td>Sample Target</td>
181 <td>'((CHASTITY>=0.6))'</td>
182 <td><a href="#Lane1">Lane 1</a></td>
190 <td>'((CHASTITY>=0.6))'</td>
191 <td><a href="#Lane2">Lane 2</a></td>
199 <td>'((CHASTITY>=0.6))'</td>
200 <td><a href="#Lane3">Lane 3</a></td>
208 <td>'((CHASTITY>=0.6))'</td>
209 <td><a href="#Lane4">Lane 4</a></td>
217 <td>'((CHASTITY>=0.6))'</td>
218 <td><a href="#Lane5">Lane 5</a></td>
226 <td>'((CHASTITY>=0.6))'</td>
227 <td><a href="#Lane6">Lane 6</a></td>
235 <td>'((CHASTITY>=0.6))'</td>
236 <td><a href="#Lane7">Lane 7</a></td>
244 <td>'((CHASTITY>=0.6))'</td>
245 <td><a href="#Lane8">Lane 8</a></td>
248 <h2><br></br>Lane Results Summary<br></br></h2>
249 <table border="1" cellpadding="5">
254 <td>Av 1st Cycle Int </td>
255 <td>% intensity after 20 cycles </td>
256 <td>% PF Clusters </td>
257 <td>% Align (PF) </td>
258 <td>Av Alignment Score (PF) </td>
259 <td> % Error Rate (PF) </td>
263 <td>17421 +/- 2139</td>
264 <td>7230 +/- 801</td>
265 <td>23.73 +/- 10.79</td>
266 <td>13.00 +/- 22.91</td>
267 <td>32.03 +/- 18.45</td>
268 <td>6703.57 +/- 3753.85</td>
269 <td>4.55 +/- 4.81</td>
273 <td>20311 +/- 2402</td>
274 <td>7660 +/- 678</td>
275 <td>17.03 +/- 4.40</td>
276 <td>40.74 +/- 30.33</td>
277 <td>29.54 +/- 9.03</td>
278 <td>5184.02 +/- 1631.54</td>
279 <td>3.27 +/- 3.94</td>
283 <td>20193 +/- 2399</td>
284 <td>7700 +/- 797</td>
285 <td>15.75 +/- 3.30</td>
286 <td>56.56 +/- 17.16</td>
287 <td>27.33 +/- 7.48</td>
288 <td>4803.49 +/- 1313.31</td>
289 <td>3.07 +/- 2.86</td>
293 <td>15537 +/- 2531</td>
294 <td>7620 +/- 1392</td>
295 <td>15.37 +/- 3.79</td>
296 <td>63.05 +/- 18.30</td>
297 <td>15.88 +/- 4.99</td>
298 <td>3162.13 +/- 962.59</td>
299 <td>3.11 +/- 2.22</td>
303 <td>32047 +/- 3356</td>
304 <td>8093 +/- 831</td>
305 <td>23.79 +/- 6.18</td>
306 <td>53.36 +/- 18.06</td>
307 <td>48.04 +/- 13.77</td>
308 <td>9866.23 +/- 2877.30</td>
309 <td>2.26 +/- 1.16</td>
313 <td>32946 +/- 4753</td>
314 <td>8227 +/- 736</td>
315 <td>24.07 +/- 4.69</td>
316 <td>54.65 +/- 12.57</td>
317 <td>50.98 +/- 10.54</td>
318 <td>10468.86 +/- 2228.53</td>
319 <td>2.21 +/- 2.33</td>
323 <td>39504 +/- 4171</td>
324 <td>8401 +/- 785</td>
325 <td>22.55 +/- 4.56</td>
326 <td>45.22 +/- 10.34</td>
327 <td>48.41 +/- 9.67</td>
328 <td>9829.40 +/- 1993.20</td>
329 <td>2.26 +/- 1.11</td>
333 <td>37998 +/- 3792</td>
334 <td>8443 +/- 1211</td>
335 <td>39.03 +/- 7.52</td>
336 <td>42.16 +/- 12.35</td>
337 <td>40.98 +/- 14.89</td>
338 <td>8128.87 +/- 3055.34</td>
339 <td>3.57 +/- 2.77</td>
345 pathname = os.path.join(gerald_dir, 'Summary.htm')
346 f = open(pathname, 'w')
350 def make_eland_results(gerald_dir):
351 eland_result = """>HWI-EAS229_24_207BTAAXX:1:7:599:759 ACATAGNCACAGACATAAACATAGACATAGAC U0 1 1 3 chrUextra.fa 28189829 R D.
352 >HWI-EAS229_24_207BTAAXX:1:7:205:842 AAACAANNCTCCCAAACACGTAAACTGGAAAA U1 0 1 0 chr2L.fa 8796855 R DD 24T
353 >HWI-EAS229_24_207BTAAXX:1:7:776:582 AGCTCANCCGATCGAAAACCTCNCCAAGCAAT NM 0 0 0
354 >HWI-EAS229_24_207BTAAXX:1:7:205:842 AAACAANNCTCCCAAACACGTAAACTGGAAAA U1 0 1 0 Lambda.fa 8796855 R DD 24T
357 pathname = os.path.join(gerald_dir,
358 's_%d_eland_result.txt' % (i,))
359 f = open(pathname, 'w')
360 f.write(eland_result)
363 class RunfolderTests(unittest.TestCase):
365 Test components of the runfolder processing code
366 which includes firecrest, bustard, and gerald
369 # make a fake runfolder directory
370 self.temp_dir = tempfile.mkdtemp(prefix='tmp_runfolder_')
372 self.runfolder_dir = os.path.join(self.temp_dir,
373 '080102_HWI-EAS229_0010_207BTAAXX')
374 os.mkdir(self.runfolder_dir)
376 self.data_dir = os.path.join(self.runfolder_dir, 'Data')
377 os.mkdir(self.data_dir)
379 self.firecrest_dir = os.path.join(self.data_dir,
380 'C1-33_Firecrest1.8.28_12-04-2008_diane'
382 os.mkdir(self.firecrest_dir)
383 self.matrix_dir = os.path.join(self.firecrest_dir, 'Matrix')
384 os.mkdir(self.matrix_dir)
385 make_matrix(self.matrix_dir)
387 self.bustard_dir = os.path.join(self.firecrest_dir,
388 'Bustard1.8.28_12-04-2008_diane')
389 os.mkdir(self.bustard_dir)
390 make_phasing_params(self.bustard_dir)
392 self.gerald_dir = os.path.join(self.bustard_dir,
393 'GERALD_12-04-2008_diane')
394 os.mkdir(self.gerald_dir)
395 make_gerald_config(self.gerald_dir)
396 make_summary_htm(self.gerald_dir)
397 make_eland_results(self.gerald_dir)
400 shutil.rmtree(self.temp_dir)
402 def test_firecrest(self):
404 Construct a firecrest object
406 f = firecrest.firecrest(self.firecrest_dir)
407 self.failUnlessEqual(f.version, '1.8.28')
408 self.failUnlessEqual(f.start, 1)
409 self.failUnlessEqual(f.stop, 33)
410 self.failUnlessEqual(f.user, 'diane')
411 self.failUnlessEqual(f.date, date(2008,4,12))
413 xml = f.get_elements()
414 # just make sure that element tree can serialize the tree
415 xml_str = ElementTree.tostring(xml)
417 f2 = firecrest.Firecrest(xml=xml)
418 self.failUnlessEqual(f.version, f2.version)
419 self.failUnlessEqual(f.start, f2.start)
420 self.failUnlessEqual(f.stop, f2.stop)
421 self.failUnlessEqual(f.user, f2.user)
422 self.failUnlessEqual(f.date, f2.date)
424 def test_bustard(self):
426 construct a bustard object
428 b = bustard.bustard(self.bustard_dir)
429 self.failUnlessEqual(b.version, '1.8.28')
430 self.failUnlessEqual(b.date, date(2008,4,12))
431 self.failUnlessEqual(b.user, 'diane')
432 self.failUnlessEqual(len(b.phasing), 8)
433 self.failUnlessAlmostEqual(b.phasing[8].phasing, 0.0099)
435 xml = b.get_elements()
436 b2 = bustard.Bustard(xml=xml)
437 self.failUnlessEqual(b.version, b2.version)
438 self.failUnlessEqual(b.date, b2.date )
439 self.failUnlessEqual(b.user, b2.user)
440 self.failUnlessEqual(len(b.phasing), len(b2.phasing))
441 for key in b.phasing.keys():
442 self.failUnlessEqual(b.phasing[key].lane,
443 b2.phasing[key].lane)
444 self.failUnlessEqual(b.phasing[key].phasing,
445 b2.phasing[key].phasing)
446 self.failUnlessEqual(b.phasing[key].prephasing,
447 b2.phasing[key].prephasing)
449 def test_gerald(self):
450 # need to update gerald and make tests for it
451 g = gerald.gerald(self.gerald_dir)
453 self.failUnlessEqual(g.version,
454 '@(#) Id: GERALD.pl,v 1.68.2.2 2007/06/13 11:08:49 km Exp')
455 self.failUnlessEqual(g.date, datetime(2008,4,19,19,8,30))
456 self.failUnlessEqual(len(g.lanes), len(g.lanes.keys()))
457 self.failUnlessEqual(len(g.lanes), len(g.lanes.items()))
460 # list of genomes, matches what was defined up in
461 # make_gerald_config.
462 # the first None is to offset the genomes list to be 1..9
463 # instead of pythons default 0..8
464 genomes = [None, '/g/dm3', '/g/equcab1', '/g/equcab1', '/g/canfam2',
465 '/g/hg18', '/g/hg18', '/g/hg18', '/g/hg18', ]
467 # test lane specific parameters from gerald config file
469 cur_lane = g.lanes[str(i)]
470 self.failUnlessEqual(cur_lane.analysis, 'eland')
471 self.failUnlessEqual(cur_lane.eland_genome, genomes[i])
472 self.failUnlessEqual(cur_lane.read_length, '32')
473 self.failUnlessEqual(cur_lane.use_bases, 'Y'*32)
475 # test data extracted from summary file
477 (17421, 2139), (20311, 2402), (20193, 2399), (15537, 2531),
478 (32047, 3356), (32946, 4753), (39504, 4171), (37998, 3792)]
481 summary_lane = g.summary[str(i)]
482 self.failUnlessEqual(summary_lane.cluster, clusters[i])
483 self.failUnlessEqual(summary_lane.lane, str(i))
485 xml = g.get_elements()
486 # just make sure that element tree can serialize the tree
487 xml_str = ElementTree.tostring(xml)
488 g2 = gerald.Gerald(xml=xml)
490 # do it all again after extracting from the xml file
491 self.failUnlessEqual(g.version, g2.version)
492 self.failUnlessEqual(g.date, g2.date)
493 self.failUnlessEqual(len(g.lanes.keys()), len(g2.lanes.keys()))
494 self.failUnlessEqual(len(g.lanes.items()), len(g2.lanes.items()))
496 # test lane specific parameters from gerald config file
498 g_lane = g.lanes[str(i)]
499 g2_lane = g2.lanes[str(i)]
500 self.failUnlessEqual(g_lane.analysis, g2_lane.analysis)
501 self.failUnlessEqual(g_lane.eland_genome, g2_lane.eland_genome)
502 self.failUnlessEqual(g_lane.read_length, g2_lane.read_length)
503 self.failUnlessEqual(g_lane.use_bases, g2_lane.use_bases)
505 # test (some) summary elements
507 g_summary = g.summary[str(i)]
508 g2_summary = g2.summary[str(i)]
509 self.failUnlessEqual(g_summary.cluster, g2_summary.cluster)
510 self.failUnlessEqual(g_summary.lane, g2_summary.lane)
512 g_eland = g.eland_results
513 g2_eland = g2.eland_results
514 for lane in g_eland.keys():
515 self.failUnlessEqual(g_eland[lane].reads,
516 g2_eland[lane].reads)
517 self.failUnlessEqual(len(g_eland[lane].mapped_reads),
518 len(g2_eland[lane].mapped_reads))
519 for k in g_eland[lane].mapped_reads.keys():
520 self.failUnlessEqual(g_eland[lane].mapped_reads[k],
521 g2_eland[lane].mapped_reads[k])
523 self.failUnlessEqual(len(g_eland[lane].match_codes),
524 len(g2_eland[lane].match_codes))
525 for k in g_eland[lane].match_codes.keys():
526 self.failUnlessEqual(g_eland[lane].match_codes[k],
527 g2_eland[lane].match_codes[k])
530 def test_eland(self):
531 dm3_map = { 'chrUextra.fa' : 'dm3/chrUextra.fa',
532 'chr2L.fa': 'dm3/chr2L.fa',
533 'Lambda.fa': 'Lambda.fa'}
534 genome_maps = { '1':dm3_map, '2':dm3_map, '3':dm3_map, '4':dm3_map,
535 '5':dm3_map, '6':dm3_map, '7':dm3_map, '8':dm3_map }
536 eland = gerald.eland(self.gerald_dir, genome_maps=genome_maps)
540 self.failUnlessEqual(lane.reads, 4)
541 self.failUnlessEqual(len(lane.mapped_reads), 3)
542 self.failUnlessEqual(lane.mapped_reads['Lambda.fa'], 1)
543 self.failUnlessEqual(lane.mapped_reads['dm3/chr2L.fa'], 1)
544 self.failUnlessEqual(lane.match_codes['U1'], 2)
545 self.failUnlessEqual(lane.match_codes['NM'], 1)
547 xml = eland.get_elements()
548 # just make sure that element tree can serialize the tree
549 xml_str = ElementTree.tostring(xml)
550 e2 = gerald.ELAND(xml=xml)
555 self.failUnlessEqual(l1.reads, l2.reads)
556 self.failUnlessEqual(len(l1.mapped_reads), len(l2.mapped_reads))
557 self.failUnlessEqual(len(l1.mapped_reads), 3)
558 for k in l1.mapped_reads.keys():
559 self.failUnlessEqual(l1.mapped_reads[k],
562 self.failUnlessEqual(len(l1.match_codes), 9)
563 self.failUnlessEqual(len(l1.match_codes), len(l2.match_codes))
564 for k in l1.match_codes.keys():
565 self.failUnlessEqual(l1.match_codes[k],
568 def test_runfolder(self):
569 runs = runfolder.get_runs(self.runfolder_dir)
571 # do we get the flowcell id from the filename?
572 self.failUnlessEqual(len(runs), 1)
573 self.failUnlessEqual(runs[0].name, 'run_207BTAAXX_2008-04-19.xml')
575 # do we get the flowcell id from the FlowcellId.xml file
576 make_flowcell_id(self.runfolder_dir, '207BTAAXY')
577 runs = runfolder.get_runs(self.runfolder_dir)
578 self.failUnlessEqual(len(runs), 1)
579 self.failUnlessEqual(runs[0].name, 'run_207BTAAXY_2008-04-19.xml')
582 xml = r1.get_elements()
583 xml_str = ElementTree.tostring(xml)
585 r2 = runfolder.PipelineRun(xml=xml)
586 self.failUnlessEqual(r1.name, r2.name)
587 self.failIfEqual(r2.firecrest, None)
588 self.failIfEqual(r2.bustard, None)
589 self.failIfEqual(r2.gerald, None)
593 return unittest.makeSuite(RunfolderTests,'test')
595 if __name__ == "__main__":
596 unittest.main(defaultTest="suite")