source: trunk/brs2010p35/src/jag/test/test_jag.py @ 2

Last change on this file since 2 was 2, checked in by tim.te.beek@…, 8 years ago

Add revision 338 snapshot of previous SVN instance, omitting PLINK binaries

File size: 20.1 KB
Line 
1'''
2Test cases for JAG project relating to anything.
3'''
4
5from jag.association_analysis import AssociationAnalysis
6from jag.file_fetch_and_write import InAndOut
7from jag.plink import Plink
8import hashlib
9import glob
10import os
11import unittest
12import logging as log
13
14class TestJag(unittest.TestCase):
15    '''
16    Test class for JAG project, relating to larger functions call multiple sub routines.
17   
18    Extending unittest.TestCase to take advantage of the self.assertXYZ methods.
19    '''
20    def test_step1(self):
21        """
22        test association analysis (aka step1)
23        """
24        log.info('I\'m a log message, with magic values: %s & %d', 'like this', 12)
25
26        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
27
28        pathcheck = prefix + "hapmap_CEU_r23_60"
29        phenofile = prefix + "hapmap_pheno.txt"
30        group = prefix + "hapmap_group.txt"
31
32        plink = Plink()
33        plink.bfile = (pathcheck)
34        plink.pheno_file = phenofile
35
36        inoutput = InAndOut()
37
38
39        aa = AssociationAnalysis(inoutput)
40        results = aa.run_step1(plink, group)
41        results = results["1"]["clusterresults"]
42
43        round_to_decimals = 2
44        self.assertAlmostEqual(455.56, results.get_sum_neglog_10_score("FG1"), round_to_decimals)
45        self.assertEqual(53, results.get_n_genes_of_group("FG1"))
46        self.assertEqual(1130, results.get_n_snps("FG1"))
47        self.assertAlmostEqual(274.45, results.get_sum_neglog_10_score("FG18"), round_to_decimals)
48        self.assertEqual(49, results.get_n_genes_of_group("FG18"))
49        self.assertEqual(629, results.get_n_snps("FG18"))
50
51        #Clean up by removing old files
52        remove_old_files()
53
54    def test_step1_without_pheno(self):
55        """
56        test association analysis  without phenotype file(aka step1)
57        """
58        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
59
60        pathcheck = prefix + "hapmap_CEU_r23_60"
61        #phenofile = prefix + "pheno60.txt"
62        group = prefix + "hapmap_group.txt"
63
64        plink2 = Plink()
65        plink2.bfile = (pathcheck)
66        #plink.set_pheno_file(phenofile)
67
68
69        inoutput = InAndOut()
70
71
72        aa = AssociationAnalysis(inoutput)
73        results2 = aa.run_step1(plink2, group)
74
75        results2 = results2["1"]["clusterresults"]
76
77
78
79        round_to_decimals = 2
80        self.assertAlmostEqual(455.56, results2.get_sum_neglog_10_score("FG1"), round_to_decimals)
81        self.assertEqual(53, results2.get_n_genes_of_group("FG1"))
82        self.assertEqual(1130, results2.get_n_snps("FG1"))
83        self.assertAlmostEqual(274.45, results2.get_sum_neglog_10_score("FG18"), round_to_decimals)
84        self.assertEqual(49, results2.get_n_genes_of_group("FG18"))
85        self.assertEqual(629, results2.get_n_snps("FG18"))
86        remove_old_files()
87
88    def test_step1_with_cont_pheno_age(self):
89        """
90     
91               test association analysis (aka step1)
92        """
93        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
94
95        pathcheck = prefix + "hapmap_CEU_r23_60"
96        phenofile = prefix + "hapmap_pheno_qt.txt"
97        group = prefix + "hapmap_group.txt"
98
99        age = prefix + "hapmap_age.txt"
100        plink = Plink()
101        plink.bfile = (pathcheck)
102        plink.covar_file = age
103        plink.pheno_file = phenofile
104        plink.switches += "--linear "
105
106        inoutput = InAndOut()
107
108
109        aa = AssociationAnalysis(inoutput)
110        results = aa.run_step1(plink, group)
111
112        results = results["1"]["clusterresults"]
113
114        round_to_decimals = 2
115        self.assertAlmostEqual(513.89, results.get_sum_neglog_10_score("FG1"), round_to_decimals)
116        self.assertEqual(53, results.get_n_genes_of_group("FG1"))
117        self.assertEqual(1130, results.get_n_snps("FG1"))
118        self.assertAlmostEqual(288.30, results.get_sum_neglog_10_score("FG18"), round_to_decimals)
119        self.assertEqual(49, results.get_n_genes_of_group("FG18"))
120        self.assertEqual(629, results.get_n_snps("FG18"))
121
122        #Clean up by removing old files
123        remove_old_files()
124
125    def test_step1_with_cont_pheno_age_sex(self):
126        """
127     
128               test association analysis (aka step1)
129        """
130        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
131
132        pathcheck = prefix + "hapmap_CEU_r23_60"
133        phenofile = prefix + "hapmap_pheno_qt.txt"
134        group = prefix + "hapmap_group.txt"
135
136        age = prefix + "hapmap_age.txt"
137        plink = Plink()
138        plink.bfile = (pathcheck)
139        plink.covar_file = age
140        plink.pheno_file = phenofile
141        plink.switches += "--linear --sex"
142
143        inoutput = InAndOut()
144
145
146        aa = AssociationAnalysis(inoutput)
147        results = aa.run_step1(plink, group)
148
149        results = results["1"]["clusterresults"]
150
151        round_to_decimals = 2
152        self.assertAlmostEqual(513.80, results.get_sum_neglog_10_score("FG1"), round_to_decimals)
153        self.assertEqual(53, results.get_n_genes_of_group("FG1"))
154        self.assertEqual(1130, results.get_n_snps("FG1"))
155        self.assertAlmostEqual(292.46, results.get_sum_neglog_10_score("FG18"), round_to_decimals)
156        self.assertEqual(49, results.get_n_genes_of_group("FG18"))
157        self.assertEqual(629, results.get_n_snps("FG18"))
158
159        #Clean up by removing old files
160        remove_old_files()
161
162    def test_step1_with_cont_pheno(self):
163        """
164     
165               test association analysis (aka step1)
166        """
167        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
168
169        pathcheck = prefix + "hapmap_CEU_r23_60"
170        phenofile = prefix + "hapmap_pheno_qt.txt"
171        group = prefix + "hapmap_group.txt"
172
173        plink = Plink()
174        plink.bfile = (pathcheck)
175        plink.pheno_file = phenofile
176
177        inoutput = InAndOut()
178
179
180        aa = AssociationAnalysis(inoutput)
181        results = aa.run_step1(plink, group)
182
183        results = results["1"]["clusterresults"]
184
185        round_to_decimals = 2
186        self.assertAlmostEqual(513.00, results.get_sum_neglog_10_score("FG1"), round_to_decimals)
187        self.assertEqual(53, results.get_n_genes_of_group("FG1"))
188        self.assertEqual(1130, results.get_n_snps("FG1"))
189        self.assertAlmostEqual(284.32, results.get_sum_neglog_10_score("FG18"), round_to_decimals)
190        self.assertEqual(49, results.get_n_genes_of_group("FG18"))
191        self.assertEqual(629, results.get_n_snps("FG18"))
192
193        #Clean up by removing old files
194        remove_old_files()
195
196    def test_step1_with_covar_age_logistic(self):
197        """
198     
199               test association analysis (aka step1)
200        """
201        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
202
203        pathcheck = prefix + "hapmap_CEU_r23_60"
204        group = prefix + "hapmap_group.txt"
205        age = prefix + "hapmap_age.txt"
206        plink = Plink()
207        plink.bfile = (pathcheck)
208        plink.covar_file = age
209        plink.switches += "--logistic "
210
211        inoutput = InAndOut()
212
213        aa = AssociationAnalysis(inoutput)
214        results = aa.run_step1(plink, group)
215        print(str(results))
216        results = results["1"]["clusterresults"]
217
218        round_to_decimals = 2
219        self.assertAlmostEqual(456.10, results.get_sum_neglog_10_score("FG1"), round_to_decimals)
220        self.assertEqual(53, results.get_n_genes_of_group("FG1"))
221        self.assertEqual(1130, results.get_n_snps("FG1"))
222        self.assertAlmostEqual(271.72, results.get_sum_neglog_10_score("FG18"), round_to_decimals)
223        self.assertEqual(49, results.get_n_genes_of_group("FG18"))
224        self.assertEqual(629, results.get_n_snps("FG18"))
225
226        #Clean up by removing old files
227        remove_old_files()
228
229    def test_step1_with_covar_age_gender_logistic(self):
230        """
231     
232               test association analysis (aka step1)
233        """
234        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
235
236        pathcheck = prefix + "hapmap_CEU_r23_60"
237        group = prefix + "hapmap_group.txt"
238        age = prefix + "hapmap_age.txt"
239        plink = Plink()
240        plink.bfile = (pathcheck)
241        plink.covar_file = age
242        plink.switches += "--logistic "
243        plink.switches += "--sex "
244
245        inoutput = InAndOut()
246
247        aa = AssociationAnalysis(inoutput)
248        results = aa.run_step1(plink, group)
249        print(str(results))
250        results = results["1"]["clusterresults"]
251
252        round_to_decimals = 2
253        self.assertAlmostEqual(461.07, results.get_sum_neglog_10_score("FG1"), round_to_decimals)
254        self.assertEqual(53, results.get_n_genes_of_group("FG1"))
255        self.assertEqual(1130, results.get_n_snps("FG1"))
256        self.assertAlmostEqual(271.97, results.get_sum_neglog_10_score("FG18"), round_to_decimals)
257        self.assertEqual(49, results.get_n_genes_of_group("FG18"))
258        self.assertEqual(629, results.get_n_snps("FG18"))
259
260        #Clean up by removing old files
261        remove_old_files()
262
263    def test_step1_with_logistic(self):
264        """
265     
266               test association analysis (aka step1)
267        """
268        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
269
270        pathcheck = prefix + "hapmap_CEU_r23_60"
271        group = prefix + "hapmap_group.txt"
272        plink = Plink()
273        plink.bfile = pathcheck
274        plink.switches += "--logistic "
275
276
277        inoutput = InAndOut()
278
279        aa = AssociationAnalysis(inoutput)
280        results = aa.run_step1(plink, group)
281        print(str(results))
282        results = results["1"]["clusterresults"]
283
284        round_to_decimals = 2
285        self.assertAlmostEqual(450.00, results.get_sum_neglog_10_score("FG1"), round_to_decimals)
286        self.assertEqual(53, results.get_n_genes_of_group("FG1"))
287        self.assertEqual(1130, results.get_n_snps("FG1"))
288        self.assertAlmostEqual(356.36, results.get_sum_neglog_10_score("FG18"), round_to_decimals)
289        self.assertEqual(49, results.get_n_genes_of_group("FG18"))
290        self.assertEqual(629, results.get_n_snps("FG18"))
291
292        #Clean up by removing old files
293        remove_old_files()
294
295    def test_step1_with_linear(self):
296        """
297     
298               test association analysis (aka step1)
299        """
300        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
301
302        pathcheck = prefix + "hapmap_CEU_r23_60"
303        group = prefix + "hapmap_group.txt"
304        plink = Plink()
305        plink.bfile = pathcheck
306        plink.switches += "--linear "
307
308
309        inoutput = InAndOut()
310
311        aa = AssociationAnalysis(inoutput)
312        results = aa.run_step1(plink, group)
313        print(str(results))
314        results = results["1"]["clusterresults"]
315
316        round_to_decimals = 2
317        self.assertAlmostEqual(450.00, results.get_sum_neglog_10_score("FG1"), round_to_decimals)
318        self.assertEqual(53, results.get_n_genes_of_group("FG1"))
319        self.assertEqual(1130, results.get_n_snps("FG1"))
320        self.assertAlmostEqual(356.36, results.get_sum_neglog_10_score("FG18"), round_to_decimals)
321        self.assertEqual(49, results.get_n_genes_of_group("FG18"))
322        self.assertEqual(629, results.get_n_snps("FG18"))
323
324        #Clean up by removing old files
325        remove_old_files()
326
327
328
329    def test_step2_(self):
330        """
331     
332              test plain permutations
333        """
334        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
335
336        pathcheck = prefix + "hapmap_CEU_r23_60"
337        group = prefix + "hapmap_group.txt"
338        plink = Plink()
339        plink.bfile = pathcheck
340
341        inoutput = InAndOut()
342
343        from jag.permutation import Permutation
344        permutation = Permutation(plink, group, inoutput)
345        permutation.no_emp = False
346        seed = 11
347        gene_group = False
348        permutation.localpermutation(2, seed, gene_group)
349
350        perm_hash = hashlib.sha256(open(permutation.files["1"]["perm"], 'rb').read()).hexdigest()
351        perm_hash_orig = '8e141f0af513ba562eb59266591e510c4d6bf4334c6e385f2d2f810e70b2faca'
352
353        empp_hash = hashlib.sha256(open(permutation.files["1"]["empp"], 'rb').read()).hexdigest()
354        empp_hash_orig = '8fd205a3931bbf5c4afa29734111bf410ee8728820a650043e63f57d01fd6df7'
355
356        self.assertEqual(perm_hash_orig, perm_hash)
357
358        self.assertEqual(empp_hash_orig, empp_hash)
359
360        permutation.localpermutation(2, seed, gene_group)
361
362
363        #Clean up by removing old files
364        remove_old_files()
365
366    def test_step2_with_prefix(self):
367        """
368     
369              test plain permutations with strange bug when having 2 numbers in the end of outfile
370        """
371        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
372        pathcheck = prefix + "hapmap_CEU_r23_60"
373        phenofile = prefix + "multiple_pheno.txt"
374        group = prefix + "hapmap_group.txt"
375
376        plink = Plink()
377        plink.bfile = (pathcheck)
378        plink.pheno_file = phenofile
379
380        inoutput = InAndOut()
381        inoutput.set_outfile("testing12")
382        inoutput.run_rproject = False
383        from jag.permutation import Permutation
384        permutation = Permutation(plink, group, inoutput)
385        permutation.no_emp = False
386        seed = 11
387        gene_group = False
388        permutation.localpermutation(2, seed, gene_group)
389
390        perm_hash = hashlib.sha256(open(permutation.files["1"]["perm"], 'rb').read()).hexdigest()
391        perm_hash_orig = '8e141f0af513ba562eb59266591e510c4d6bf4334c6e385f2d2f810e70b2faca'
392
393        empp_hash = hashlib.sha256(open(permutation.files["1"]["empp"], 'rb').read()).hexdigest()
394        empp_hash_orig = '8fd205a3931bbf5c4afa29734111bf410ee8728820a650043e63f57d01fd6df7'
395
396        self.assertEqual(perm_hash_orig, perm_hash,"first perm")
397        self.assertEqual(empp_hash_orig, empp_hash,"first emp")
398       
399
400        perm_hash = hashlib.sha256(open(permutation.files["2"]["perm"], 'rb').read()).hexdigest()
401        perm_hash_orig = 'fa501ca70c32d54ced9b272afff1134b320a2fc64ad8d43f4395444d8f5381c5'
402
403        empp_hash = hashlib.sha256(open(permutation.files["2"]["empp"], 'rb').read()).hexdigest()
404        empp_hash_orig = '84da2bbbeb881164e53ff4655ce5f80c2e4422c402b4296912b2b676a96c935d'
405        self.assertEqual(perm_hash_orig, perm_hash,"second perm")
406        self.assertEqual(empp_hash_orig, empp_hash,"second emp")
407        remove_old_files()
408#        inoutput = InAndOut()
409#        inoutput.set_outfile("testing12")
410#        from jag.permutation import Permutation
411#        permutation = Permutation(plink, group, inoutput)
412#        permutation.no_emp = False
413#        seed = 12
414#        gene_group = False
415#        permutation.localpermutation(2, seed, gene_group)
416#        perm_hash = hashlib.sha256(open(permutation.files["1"]["perm"], 'rb').read()).hexdigest()
417#        perm_hash_orig = '8e141f0af513ba562eb59266591e510c4d6bf4334c6e385f2d2f810e70b2faca'
418#
419#        empp_hash = hashlib.sha256(open(permutation.files["1"]["empp"], 'rb').read()).hexdigest()
420#        empp_hash_orig = '8fd205a3931bbf5c4afa29734111bf410ee8728820a650043e63f57d01fd6df7'
421#        self.assertEqual(perm_hash_orig, perm_hash)
422#        self.assertEqual(empp_hash_orig, empp_hash)
423#        #Clean up by removing old files
424       
425
426    def test_step3_(self):
427        """
428     
429               test step 3 merging permutation files
430
431        """
432        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
433
434        pathcheck = prefix + "hapmap_CEU_r23_60"
435        group = prefix + "hapmap_group.txt"
436        plink = Plink()
437        plink.bfile = pathcheck
438
439        inoutput = InAndOut()
440
441        from jag.permutation import Permutation
442        permutation = Permutation(plink, group, inoutput)
443        permutation.no_emp = False
444        seed = 11
445        gene_group = False
446        permutation.localpermutation(2, seed, gene_group)
447
448        perm_hash = hashlib.sha256(open(permutation.files["1"]["perm"], 'rb').read()).hexdigest()
449        perm_hash_orig = '8e141f0af513ba562eb59266591e510c4d6bf4334c6e385f2d2f810e70b2faca'
450
451        empp_hash = hashlib.sha256(open(permutation.files["1"]["empp"], 'rb').read()).hexdigest()
452        empp_hash_orig = '8fd205a3931bbf5c4afa29734111bf410ee8728820a650043e63f57d01fd6df7'
453
454        self.assertEqual(perm_hash_orig, perm_hash)
455
456        self.assertEqual(empp_hash_orig, empp_hash)
457        #do some extra perm to merge
458        seed = 12
459        permutation.localpermutation(3, seed, gene_group)
460
461        import jag.mergepermutations as mergepermutations
462        merger = mergepermutations.MergePermutation(inoutput)
463        merger.mergeresults()
464        print(str(merger.files))
465        print(merger.files["P1"]["perm"])
466        print(merger.files["P1"]["empp"])
467
468        #Clean up by removing old files
469        remove_old_files()
470
471    def test_step3_with_prefix(self):
472        """
473     
474               test step 3 merging permutation files with prefix
475        """
476        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
477
478        pathcheck = prefix + "hapmap_CEU_r23_60"
479        group = prefix + "hapmap_group.txt"
480        plink = Plink()
481        plink.bfile = pathcheck
482
483        inoutput = InAndOut()
484        inoutput.set_outfile("nice_prefix")
485        from jag.permutation import Permutation
486        permutation = Permutation(plink, group, inoutput)
487        permutation.no_emp = False
488        seed = 11
489        gene_group = False
490        permutation.localpermutation(2, seed, gene_group)
491
492        perm_hash = hashlib.sha256(open(permutation.files["1"]["perm"], 'rb').read()).hexdigest()
493        perm_hash_orig = '8e141f0af513ba562eb59266591e510c4d6bf4334c6e385f2d2f810e70b2faca'
494
495        empp_hash = hashlib.sha256(open(permutation.files["1"]["empp"], 'rb').read()).hexdigest()
496        empp_hash_orig = '8fd205a3931bbf5c4afa29734111bf410ee8728820a650043e63f57d01fd6df7'
497
498        self.assertEqual(perm_hash_orig, perm_hash)
499
500        self.assertEqual(empp_hash_orig, empp_hash)
501        #do some extra perm to merge
502        seed = 12
503        permutation.localpermutation(3, seed, gene_group)
504
505        import jag.mergepermutations as mergepermutations
506        merger = mergepermutations.MergePermutation(inoutput)
507        merger.mergeresults()
508        print(str(merger.files))
509        print(merger.files["P1"]["perm"])
510        print(merger.files["P1"]["empp"])
511
512        #Clean up by removing old files
513        remove_old_files()
514
515    def test_step2_no_emp(self):
516        """
517     
518               test step 3 merging permutation files with prefix
519        """
520        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
521
522        pathcheck = prefix + "hapmap_CEU_r23_60"
523        group = prefix + "hapmap_group.txt"
524        plink = Plink()
525        plink.bfile = pathcheck
526
527        inoutput = InAndOut()
528        from jag.permutation import Permutation
529        permutation = Permutation(plink, group, inoutput)
530        #no_emp is set to True
531        permutation.no_emp = True
532        seed = 11
533        gene_group = False
534        permutation.localpermutation(2, seed, gene_group)
535
536        perm_hash = hashlib.sha256(open(permutation.files["1"]["perm"], 'rb').read()).hexdigest()
537        perm_hash_orig = '8e141f0af513ba562eb59266591e510c4d6bf4334c6e385f2d2f810e70b2faca'
538
539
540
541        self.assertEqual(perm_hash_orig, perm_hash)
542
543       
544
545        #Clean up by removing old files
546        remove_old_files()
547#    def test_step4a_draw_genes(self):
548#        prefix = os.path.dirname(os.path.abspath(__file__)) + "/../../../hapmap/"
549#
550#        import jag.drawrandom as drawrandom
551#        inoutput = InAndOut()
552#        complete_gene_snp_mapping = prefix + "all_snps_in_gene_NCBI.txt"
553#        group = prefix + "hapmap_group.txt"
554#
555#
556#        draw = drawrandom.DrawRandom(complete_gene_snp_mapping, group, inoutput)
557#        draw.setseed(11)
558#
559#        import jag.drawrandom as drawrandom
560#        genes = drawrandom.Genes(draw)
561#        amount = 100
562#        set_name = "FG1"
563#        genes.genes(set_name, amount)
564
565
566
567
568def remove_old_files():
569    """
570    remove files made by jag runs
571    """
572    currentdir = os.path.abspath(os.path.curdir) + "/"
573    files = glob.glob(currentdir + "*sumlog*.out")
574    files.extend(glob.glob(currentdir + "*.*assoc*"))
575    files.extend(glob.glob(currentdir + "*empp*.out"))
576    files.extend(glob.glob(currentdir + "*perm*.out"))
577    files.extend(glob.glob(currentdir + "*png"))
578    for filename in files:
579        if os.path.exists(filename):
580            os.remove(filename)
581
Note: See TracBrowser for help on using the repository browser.