test/test_ScoreTrack.py

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