test/test_ScoreTrack.py

   1 #!/usr/bin/env python
   2 # Time-stamp: <2019-09-25 14:44:07 taoliu>
   3
   4 import io
   5 import unittest
   6 from numpy.testing import assert_equal,  assert_almost_equal, assert_array_equal
   7
   8 from MACS2.IO.ScoreTrack import *
   9 from MACS2.IO.BedGraph import bedGraphTrackI
  10
  11 class Test_TwoConditionScores(unittest.TestCase):
  12     def setUp(self):
  13         self.t1bdg = bedGraphTrackI()
  14         self.t2bdg = bedGraphTrackI()
  15         self.c1bdg = bedGraphTrackI()
  16         self.c2bdg = bedGraphTrackI()
  17         self.test_regions1 = [("chrY",0,70,0.00,0.01),
  18                               ("chrY",70,80,7.00,0.5),
  19                               ("chrY",80,150,0.00,0.02)]
  20         self.test_regions2 = [("chrY",0,75,20.0,4.00),
  21                               ("chrY",75,90,35.0,6.00),
  22                               ("chrY",90,150,10.0,15.00)]
  23         for a in self.test_regions1:
  24             self.t1bdg.safe_add_loc(a[0],a[1],a[2],a[3])
  25             self.c1bdg.safe_add_loc(a[0],a[1],a[2],a[4])
  26
  27         for a in self.test_regions2:
  28             self.t2bdg.safe_add_loc(a[0],a[1],a[2],a[3])
  29             self.c2bdg.safe_add_loc(a[0],a[1],a[2],a[4])
  30
  31         self.twoconditionscore = ScoreTrack.TwoConditionScores( self.t1bdg,
  32                                                                 self.c1bdg,
  33                                                                 self.t2bdg,
  34                                                                 self.c2bdg,
  35                                                                 1.0,
  36                                                                 1.0 )
  37         self.twoconditionscore.build()
  38         self.twoconditionscore.finalize()
  39         (self.cat1,self.cat2,self.cat3) = twoconditionscore.call_peaks(min_length=10, max_gap=10, cutoff=3)
  40
  41
  42 class Test_ScoreTrackII(unittest.TestCase):
  43
  44     def setUp(self):
  45         # for initiate scoretrack
  46         self.test_regions1 = [(b"chrY",10,100,10),
  47                               (b"chrY",60,10,10),
  48                               (b"chrY",110,15,20),
  49                               (b"chrY",160,5,20),
  50                               (b"chrY",210,20,5)]
  51         self.treat_edm = 10
  52         self.ctrl_edm = 5
  53         # for different scoring method
  54         self.p_result = [60.49, 0.38, 0.08, 0.0, 6.41] # -log10(p-value), pseudo count 1 added
  55         self.q_result = [58.17, 0.0, 0.0, 0.0, 5.13] # -log10(q-value) from BH, pseudo count 1 added
  56         self.l_result = [58.17, 0.0, -0.28, -3.25, 4.91] # log10 likelihood ratio, pseudo count 1 added
  57         self.f_result = [0.96, 0.00, -0.12, -0.54, 0.54] # note, pseudo count 1 would be introduced.
  58         self.d_result = [90.00, 0, -5.00, -15.00, 15.00]
  59         self.m_result = [10.00, 1.00, 1.50, 0.50, 2.00]
  60         # for norm
  61         self.norm_T = np.array([[ 10, 100,  20,   0],
  62                                 [ 60,  10,  20,   0],
  63                                 [110,  15,  40,   0],
  64                                 [160,   5,  40,   0],
  65                                 [210,  20,  10,   0]]).transpose()
  66         self.norm_C = np.array([[ 10,  50,  10,   0],
  67                                 [ 60,   5,  10,   0],
  68                                 [110,   7.5,  20,   0],
  69                                 [160,   2.5,  20,   0],
  70                                 [210,  10,   5,   0]]).transpose()
  71         self.norm_M = np.array([[ 10,  10,   2,   0],
  72                                 [ 60,   1,   2,   0],
  73                                 [110,   1.5,   4,   0],
  74                                 [160,   0.5,   4,   0],
  75                                 [210,   2,   1,   0]]).transpose()
  76         self.norm_N = np.array([[ 10, 100,  10,   0],  # note precision lost
  77                                 [ 60,  10,  10,   0],
  78                                 [110,  15,  20,   0],
  79                                 [160,   5,  20,   0],
  80                                 [210,  20,   5,   0]]).transpose()
  81
  82         # for write_bedGraph
  83         self.bdg1 = """chrY     0       10      100.00000
  84 chrY    10      60      10.00000
  85 chrY    60      110     15.00000
  86 chrY    110     160     5.00000
  87 chrY    160     210     20.00000
  88 """
  89         self.bdg2 = """chrY     0       60      10.00000
  90 chrY    60      160     20.00000
  91 chrY    160     210     5.00000
  92 """
  93         self.bdg3 = """chrY     0       10      60.48912
  94 chrY    10      60      0.37599
  95 chrY    60      110     0.07723
  96 chrY    110     160     0.00006
  97 chrY    160     210     6.40804
  98 """
  99         # for peak calls
 100         self.peak1 = """chrY    0       60      peak_1  60.48912
 101 chrY    160     210     peak_2  6.40804
 102 """
 103         self.summit1 = """chrY  5       6       peak_1  60.48912
 104 chrY    185     186     peak_2  6.40804
 105 """
 106         self.xls1    ="""chr    start   end     length  abs_summit      pileup  -log10(pvalue)  fold_enrichment -log10(qvalue)  name
 107 chrY    1       60      60      6       100.00  63.27251        9.18182 -1.00000        MACS_peak_1
 108 chrY    161     210     50      186     20.00   7.09102 3.50000 -1.00000        MACS_peak_2
 109 """
 110
 111     def assertListAlmostEqual ( self, a, b, places =2 ):
 112         return all( [self.assertAlmostEqual(x, y, places=places) for (x, y) in zip( a, b)] )
 113
 114     def test_compute_scores(self):
 115         s1 = scoreTrackII( self.treat_edm, self.ctrl_edm )
 116         s1.add_chromosome( b"chrY", 5 )
 117         for a in self.test_regions1:
 118             s1.add( a[0],a[1],a[2],a[3] )
 119
 120         s1.set_pseudocount ( 1.0 )
 121
 122         s1.change_score_method( ord('p') )
 123         r = s1.get_data_by_chr(b"chrY")
 124         self.assertListAlmostEqual( [round(x,2) for x in r[3]], self.p_result )
 125
 126         s1.change_score_method( ord('q') )
 127         r = s1.get_data_by_chr(b"chrY")
 128         self.assertListAlmostEqual( [round(x,2) for x in list(r[3])], self.q_result )
 129
 130         s1.change_score_method( ord('l') )
 131         r = s1.get_data_by_chr(b"chrY")
 132         self.assertListAlmostEqual( [round(x,2) for x in list(r[3])], self.l_result )
 133
 134         s1.change_score_method( ord('f') )
 135         r = s1.get_data_by_chr(b"chrY")
 136         self.assertListAlmostEqual( [round(x,2) for x in list(r[3])], self.f_result )
 137
 138         s1.change_score_method( ord('d') )
 139         r = s1.get_data_by_chr(b"chrY")
 140         self.assertListAlmostEqual( [round(x,2) for x in list(r[3])], self.d_result )
 141
 142         s1.change_score_method( ord('m') )
 143         r = s1.get_data_by_chr(b"chrY")
 144         self.assertListAlmostEqual( [round(x,2) for x in list(r[3])], self.m_result )
 145
 146     def test_normalize(self):
 147         s1 = scoreTrackII( self.treat_edm, self.ctrl_edm )
 148         s1.add_chromosome( b"chrY", 5 )
 149         for a in self.test_regions1:
 150             s1.add( a[0],a[1],a[2],a[3] )
 151
 152         s1.change_normalization_method( ord('T') )
 153         r = s1.get_data_by_chr(b"chrY")
 154         assert_array_equal( r, self.norm_T )
 155
 156         s1.change_normalization_method( ord('C') )
 157         r = s1.get_data_by_chr(b"chrY")
 158         assert_array_equal( r, self.norm_C )
 159
 160         s1.change_normalization_method( ord('M') )
 161         r = s1.get_data_by_chr(b"chrY")
 162         assert_array_equal( r, self.norm_M )
 163
 164         s1.change_normalization_method( ord('N') )
 165         r = s1.get_data_by_chr(b"chrY")
 166         assert_array_equal( r, self.norm_N )
 167
 168     def test_writebedgraph ( self ):
 169         s1 = scoreTrackII( self.treat_edm, self.ctrl_edm )
 170         s1.add_chromosome( b"chrY", 5 )
 171         for a in self.test_regions1:
 172             s1.add( a[0],a[1],a[2],a[3] )
 173
 174         s1.change_score_method( ord('p') )
 175
 176         strio = io.StringIO()
 177         s1.write_bedGraph( strio, "NAME", "DESC", 1 )
 178         self.assertEqual( strio.getvalue(), self.bdg1 )
 179         strio = io.StringIO()
 180         s1.write_bedGraph( strio, "NAME", "DESC", 2 )
 181         self.assertEqual( strio.getvalue(), self.bdg2 )
 182         strio = io.StringIO()
 183         s1.write_bedGraph( strio, "NAME", "DESC", 3 )
 184         self.assertEqual( strio.getvalue(), self.bdg3 )
 185
 186     def test_callpeak ( self ):
 187         s1 = scoreTrackII( self.treat_edm, self.ctrl_edm )
 188         s1.add_chromosome( b"chrY", 5 )
 189         for a in self.test_regions1:
 190             s1.add( a[0],a[1],a[2],a[3] )
 191
 192         s1.change_score_method( ord('p') )
 193         p = s1.call_peaks( cutoff = 0.10, min_length=10, max_gap=10 )
 194         strio = io.StringIO()
 195         p.write_to_bed( strio, trackline = False )
 196         self.assertEqual( strio.getvalue(), self.peak1 )
 197
 198         strio = io.StringIO()
 199         p.write_to_summit_bed( strio, trackline = False )
 200         self.assertEqual( strio.getvalue(), self.summit1 )
 201
 202         strio = io.StringIO()
 203         p.write_to_xls( strio )
 204         self.assertEqual( strio.getvalue(), self.xls1 )
 205