Handle paired-end eland files.
[htsworkflow.git] / htsworkflow / pipelines / test / test_runfolder110.py
1 #!/usr/bin/env python
2
3 from datetime import datetime, date
4 import os
5 import tempfile
6 import shutil
7 import unittest
8
9 from htsworkflow.pipelines import firecrest
10 from htsworkflow.pipelines import bustard
11 from htsworkflow.pipelines import gerald
12 from htsworkflow.pipelines import runfolder
13 from htsworkflow.pipelines.runfolder import ElementTree
14
15 from htsworkflow.pipelines.test.simulate_runfolder import *
16
17
18 def make_runfolder(obj=None):
19     """
20     Make a fake runfolder, attach all the directories to obj if defined
21     """
22     # make a fake runfolder directory
23     temp_dir = tempfile.mkdtemp(prefix='tmp_runfolder_')
24
25     runfolder_dir = os.path.join(temp_dir,
26                                  '081017_HWI-EAS229_0062_30J55AAXX')
27     os.mkdir(runfolder_dir)
28
29     data_dir = os.path.join(runfolder_dir, 'Data')
30     os.mkdir(data_dir)
31
32     firecrest_dir = os.path.join(data_dir,
33                                  'C1-37_Firecrest1.9.6_20-10-2008_diane')
34     os.mkdir(firecrest_dir)
35
36     matrix_dir = os.path.join(firecrest_dir, 'Matrix')
37     os.mkdir(matrix_dir)
38     make_matrix(matrix_dir)
39
40     bustard_dir = os.path.join(firecrest_dir,
41                                'Bustard1.9.6_20-10-2008_diane')
42     os.mkdir(bustard_dir)
43     make_phasing_params(bustard_dir)
44
45     gerald_dir = os.path.join(bustard_dir,
46                               'GERALD_20-10-2008_diane')
47     os.mkdir(gerald_dir)
48     make_gerald_config(gerald_dir)
49     make_summary_htm_110(gerald_dir)
50     make_eland_multi(gerald_dir)
51
52     if obj is not None:
53         obj.temp_dir = temp_dir
54         obj.runfolder_dir = runfolder_dir
55         obj.data_dir = data_dir
56         obj.image_analysis_dir = firecrest_dir
57         obj.matrix_dir = matrix_dir
58         obj.bustard_dir = bustard_dir
59         obj.gerald_dir = gerald_dir
60
61
62 class RunfolderTests(unittest.TestCase):
63     """
64     Test components of the runfolder processing code
65     which includes firecrest, bustard, and gerald
66     """
67     def setUp(self):
68         # attaches all the directories to the object passed in
69         make_runfolder(self)
70
71     def tearDown(self):
72         shutil.rmtree(self.temp_dir)
73
74     def test_firecrest(self):
75         """
76         Construct a firecrest object
77         """
78         f = firecrest.firecrest(self.image_analysis_dir)
79         self.failUnlessEqual(f.version, '1.9.6')
80         self.failUnlessEqual(f.start, 1)
81         self.failUnlessEqual(f.stop, 37)
82         self.failUnlessEqual(f.user, 'diane')
83         self.failUnlessEqual(f.date, date(2008,10,20))
84
85         xml = f.get_elements()
86         # just make sure that element tree can serialize the tree
87         xml_str = ElementTree.tostring(xml)
88
89         f2 = firecrest.Firecrest(xml=xml)
90         self.failUnlessEqual(f.version, f2.version)
91         self.failUnlessEqual(f.start,   f2.start)
92         self.failUnlessEqual(f.stop,    f2.stop)
93         self.failUnlessEqual(f.user,    f2.user)
94
95     def test_bustard(self):
96         """
97         construct a bustard object
98         """
99         b = bustard.bustard(self.bustard_dir)
100         self.failUnlessEqual(b.version, '1.9.6')
101         self.failUnlessEqual(b.date,    date(2008,10,20))
102         self.failUnlessEqual(b.user,    'diane')
103         self.failUnlessEqual(len(b.phasing), 8)
104         self.failUnlessAlmostEqual(b.phasing[8].phasing, 0.0099)
105
106         xml = b.get_elements()
107         b2 = bustard.Bustard(xml=xml)
108         self.failUnlessEqual(b.version, b2.version)
109         self.failUnlessEqual(b.date,    b2.date )
110         self.failUnlessEqual(b.user,    b2.user)
111         self.failUnlessEqual(len(b.phasing), len(b2.phasing))
112         for key in b.phasing.keys():
113             self.failUnlessEqual(b.phasing[key].lane,
114                                  b2.phasing[key].lane)
115             self.failUnlessEqual(b.phasing[key].phasing,
116                                  b2.phasing[key].phasing)
117             self.failUnlessEqual(b.phasing[key].prephasing,
118                                  b2.phasing[key].prephasing)
119
120     def test_gerald(self):
121         # need to update gerald and make tests for it
122         g = gerald.gerald(self.gerald_dir)
123
124         self.failUnlessEqual(g.version,
125             '@(#) Id: GERALD.pl,v 1.68.2.2 2007/06/13 11:08:49 km Exp')
126         self.failUnlessEqual(g.date, datetime(2008,4,19,19,8,30))
127         self.failUnlessEqual(len(g.lanes), len(g.lanes.keys()))
128         self.failUnlessEqual(len(g.lanes), len(g.lanes.items()))
129
130
131         # list of genomes, matches what was defined up in
132         # make_gerald_config.
133         # the first None is to offset the genomes list to be 1..9
134         # instead of pythons default 0..8
135         genomes = [None, '/g/dm3', '/g/equcab1', '/g/equcab1', '/g/canfam2',
136                          '/g/hg18', '/g/hg18', '/g/hg18', '/g/hg18', ]
137
138         # test lane specific parameters from gerald config file
139         for i in range(1,9):
140             cur_lane = g.lanes[i]
141             self.failUnlessEqual(cur_lane.analysis, 'eland')
142             self.failUnlessEqual(cur_lane.eland_genome, genomes[i])
143             self.failUnlessEqual(cur_lane.read_length, '32')
144             self.failUnlessEqual(cur_lane.use_bases, 'Y'*32)
145
146         # I want to be able to use a simple iterator
147         for l in g.lanes.values():
148           self.failUnlessEqual(l.analysis, 'eland')
149           self.failUnlessEqual(l.read_length, '32')
150           self.failUnlessEqual(l.use_bases, 'Y'*32)
151
152         # raw cluster numbers extracted from summary file
153         # its the first +/- value in the lane results summary
154         # section
155         clusters = [None,
156                     (190220, 15118), (190560, 14399),
157                     (187597, 12369), (204142, 16877),
158                     (247308, 11600), (204298, 15640),
159                     (202707, 15404), (198075, 14702),]
160
161         self.failUnlessEqual(len(g.summary), 1)
162         for i in range(1,9):
163             summary_lane = g.summary[0][i]
164             self.failUnlessEqual(summary_lane.cluster, clusters[i])
165             self.failUnlessEqual(summary_lane.lane, i)
166
167         xml = g.get_elements()
168         # just make sure that element tree can serialize the tree
169         xml_str = ElementTree.tostring(xml)
170         g2 = gerald.Gerald(xml=xml)
171
172         # do it all again after extracting from the xml file
173         self.failUnlessEqual(g.version, g2.version)
174         self.failUnlessEqual(g.date, g2.date)
175         self.failUnlessEqual(len(g.lanes.keys()), len(g2.lanes.keys()))
176         self.failUnlessEqual(len(g.lanes.items()), len(g2.lanes.items()))
177
178         # test lane specific parameters from gerald config file
179         for i in range(1,9):
180             g_lane = g.lanes[i]
181             g2_lane = g2.lanes[i]
182             self.failUnlessEqual(g_lane.analysis, g2_lane.analysis)
183             self.failUnlessEqual(g_lane.eland_genome, g2_lane.eland_genome)
184             self.failUnlessEqual(g_lane.read_length, g2_lane.read_length)
185             self.failUnlessEqual(g_lane.use_bases, g2_lane.use_bases)
186
187         self.failUnlessEqual(len(g.summary), 1)
188         # test (some) summary elements
189         for i in range(1,9):
190             g_summary = g.summary[0][i]
191             g2_summary = g2.summary[0][i]
192             self.failUnlessEqual(g_summary.cluster, g2_summary.cluster)
193             self.failUnlessEqual(g_summary.lane, g2_summary.lane)
194
195             g_eland = g.eland_results
196             g2_eland = g2.eland_results
197             for lane in g_eland.results[0].keys():
198                 g_results = g_eland.results[0][lane]
199                 g2_results = g2_eland.results[0][lane]
200                 self.failUnlessEqual(g_results.reads,
201                                      g2_results.reads)
202                 self.failUnlessEqual(len(g_results.mapped_reads),
203                                      len(g2_results.mapped_reads))
204                 for k in g_results.mapped_reads.keys():
205                     self.failUnlessEqual(g_results.mapped_reads[k],
206                                          g2_results.mapped_reads[k])
207
208                 self.failUnlessEqual(len(g_results.match_codes),
209                                      len(g2_results.match_codes))
210                 for k in g_results.match_codes.keys():
211                     self.failUnlessEqual(g_results.match_codes[k],
212                                          g2_results.match_codes[k])
213
214
215     def test_eland(self):
216         hg_map = {'Lambda.fa': 'Lambda.fa'}
217         for i in range(1,22):
218           short_name = 'chr%d.fa' % (i,)
219           long_name = 'hg18/chr%d.fa' % (i,)
220           hg_map[short_name] = long_name
221
222         genome_maps = { 1:hg_map, 2:hg_map, 3:hg_map, 4:hg_map,
223                         5:hg_map, 6:hg_map, 7:hg_map, 8:hg_map }
224         eland = gerald.eland(self.gerald_dir, genome_maps=genome_maps)
225
226         for i in range(1,9):
227             lane = eland.results[0][i]
228             self.failUnlessEqual(lane.reads, 4)
229             self.failUnlessEqual(lane.sample_name, "s")
230             self.failUnlessEqual(lane.lane_id, i)
231             self.failUnlessEqual(len(lane.mapped_reads), 15)
232             self.failUnlessEqual(lane.mapped_reads['hg18/chr5.fa'], 4)
233             self.failUnlessEqual(lane.match_codes['U0'], 1)
234             self.failUnlessEqual(lane.match_codes['R0'], 2)
235             self.failUnlessEqual(lane.match_codes['U1'], 1)
236             self.failUnlessEqual(lane.match_codes['R1'], 9)
237             self.failUnlessEqual(lane.match_codes['U2'], 0)
238             self.failUnlessEqual(lane.match_codes['R2'], 12)
239             self.failUnlessEqual(lane.match_codes['NM'], 1)
240             self.failUnlessEqual(lane.match_codes['QC'], 0)
241
242         xml = eland.get_elements()
243         # just make sure that element tree can serialize the tree
244         xml_str = ElementTree.tostring(xml)
245         e2 = gerald.ELAND(xml=xml)
246
247         for i in range(1,9):
248             l1 = eland.results[0][i]
249             l2 = e2.results[0][i]
250             self.failUnlessEqual(l1.reads, l2.reads)
251             self.failUnlessEqual(l1.sample_name, l2.sample_name)
252             self.failUnlessEqual(l1.lane_id, l2.lane_id)
253             self.failUnlessEqual(len(l1.mapped_reads), len(l2.mapped_reads))
254             self.failUnlessEqual(len(l1.mapped_reads), 15)
255             for k in l1.mapped_reads.keys():
256                 self.failUnlessEqual(l1.mapped_reads[k],
257                                      l2.mapped_reads[k])
258
259             self.failUnlessEqual(len(l1.match_codes), 9)
260             self.failUnlessEqual(len(l1.match_codes), len(l2.match_codes))
261             for k in l1.match_codes.keys():
262                 self.failUnlessEqual(l1.match_codes[k],
263                                      l2.match_codes[k])
264
265     def test_runfolder(self):
266         runs = runfolder.get_runs(self.runfolder_dir)
267
268         # do we get the flowcell id from the filename?
269         self.failUnlessEqual(len(runs), 1)
270         name = 'run_30J55AAXX_2008-10-20.xml'
271         self.failUnlessEqual(runs[0].name, name)
272
273         # do we get the flowcell id from the FlowcellId.xml file
274         make_flowcell_id(self.runfolder_dir, '30J55AAXX')
275         runs = runfolder.get_runs(self.runfolder_dir)
276         self.failUnlessEqual(len(runs), 1)
277         name = 'run_30J55AAXX_2008-10-20.xml'
278         self.failUnlessEqual(runs[0].name, name)
279
280         r1 = runs[0]
281         xml = r1.get_elements()
282         xml_str = ElementTree.tostring(xml)
283
284         r2 = runfolder.PipelineRun(xml=xml)
285         self.failUnlessEqual(r1.name, r2.name)
286         self.failIfEqual(r2.image_analysis, None)
287         self.failIfEqual(r2.bustard, None)
288         self.failIfEqual(r2.gerald, None)
289
290
291 def suite():
292     return unittest.makeSuite(RunfolderTests,'test')
293
294 if __name__ == "__main__":
295     unittest.main(defaultTest="suite")
296