test/test_PairedEndTrack.py

   1 #!/usr/bin/env python
   2 # Time-stamp: <2024-10-15 16:07:27 Tao Liu>
   3
   4 import unittest
   5 from MACS3.Signal.PairedEndTrack import PETrackI, PETrackII
   6 import numpy as np
   7
   8
   9 class Test_PETrackI(unittest.TestCase):
  10     def setUp(self):
  11         self.input_regions = [(b"chrY", 0, 100),
  12                               (b"chrY", 70, 270),
  13                               (b"chrY", 70, 100),
  14                               (b"chrY", 80, 160),
  15                               (b"chrY", 80, 160),
  16                               (b"chrY", 80, 180),
  17                               (b"chrY", 80, 180),
  18                               (b"chrY", 85, 185),
  19                               (b"chrY", 85, 285),
  20                               (b"chrY", 85, 285),
  21                               (b"chrY", 85, 285),
  22                               (b"chrY", 85, 385),
  23                               (b"chrY", 90, 190),
  24                               (b"chrY", 90, 190),
  25                               (b"chrY", 90, 191),
  26                               (b"chrY", 150, 190),
  27                               (b"chrY", 150, 250),
  28                               ]
  29         self.t = sum([x[2]-x[1] for x in self.input_regions])
  30
  31     def test_add_loc(self):
  32         pe = PETrackI()
  33         for (c, l, r) in self.input_regions:
  34             pe.add_loc(c, l, r)
  35         pe.finalize()
  36         # roughly check the numbers...
  37         self.assertEqual(pe.total, 17)
  38         self.assertEqual(pe.length, self.t)
  39
  40     def test_filter_dup(self):
  41         pe = PETrackI()
  42         for (c, l, r) in self.input_regions:
  43             pe.add_loc(c, l, r)
  44         pe.finalize()
  45         # roughly check the numbers...
  46         self.assertEqual(pe.total, 17)
  47         self.assertEqual(pe.length, self.t)
  48
  49         # filter out more than 3 tags
  50         pe.filter_dup(3)
  51         self.assertEqual(pe.total, 17)
  52
  53         # filter out more than 2 tags
  54         pe.filter_dup(2)
  55         self.assertEqual(pe.total, 16)
  56
  57         # filter out more than 1 tag
  58         pe.filter_dup(1)
  59         self.assertEqual(pe.total, 12)
  60
  61     def test_sample_num(self):
  62         pe = PETrackI()
  63         for (c, l, r) in self.input_regions:
  64             pe.add_loc(c, l, r)
  65         pe.finalize()
  66         pe.sample_num(10)
  67         self.assertEqual(pe.total, 10)
  68
  69     def test_sample_percent(self):
  70         pe = PETrackI()
  71         for (c, l, r) in self.input_regions:
  72             pe.add_loc(c, l, r)
  73         pe.finalize()
  74         pe.sample_percent(0.5)
  75         self.assertEqual(pe.total, 8)
  76
  77
  78 class Test_PETrackII(unittest.TestCase):
  79     def setUp(self):
  80         self.input_regions = [(b"chrY", 0, 100, b"0w#AAACGAAAGACTCGGA", 2),
  81                               (b"chrY", 70, 170, b"0w#AAACGAAAGACTCGGA", 1),
  82                               (b"chrY", 80, 190, b"0w#AAACGAAAGACTCGGA", 1),
  83                               (b"chrY", 85, 180, b"0w#AAACGAAAGACTCGGA", 3),
  84                               (b"chrY", 100, 190, b"0w#AAACGAAAGACTCGGA", 1),
  85                               (b"chrY", 0, 100, b"0w#AAACGAACAAGTAACA", 1),
  86                               (b"chrY", 70, 170, b"0w#AAACGAACAAGTAACA", 2),
  87                               (b"chrY", 80, 190, b"0w#AAACGAACAAGTAACA", 1),
  88                               (b"chrY", 85, 180, b"0w#AAACGAACAAGTAACA", 1),
  89                               (b"chrY", 100, 190, b"0w#AAACGAACAAGTAACA", 3),
  90                               (b"chrY", 10, 110, b"0w#AAACGAACAAGTAAGA", 1),
  91                               (b"chrY", 50, 160, b"0w#AAACGAACAAGTAAGA", 2),
  92                               (b"chrY", 100, 170, b"0w#AAACGAACAAGTAAGA", 3)
  93                               ]
  94         self.pileup_p = np.array([10, 50, 70, 80, 85, 100, 110, 160, 170, 180, 190], dtype="i4")
  95         self.pileup_v = np.array([3.0, 4.0, 6.0, 9.0, 11.0, 15.0, 19.0, 18.0, 16.0, 10.0, 6.0], dtype="f4")
  96         self.peak_str = "chrom:chrY     start:80        end:180 name:peak_1     score:19        summit:105\n"
  97         self.subset_barcodes = {b'0w#AAACGAACAAGTAACA', b"0w#AAACGAACAAGTAAGA"}
  98         self.subset_pileup_p = np.array([10, 50, 70, 80, 85, 100, 110, 160, 170, 180, 190], dtype="i4")
  99         self.subset_pileup_v = np.array([1.0, 2.0, 4.0, 6.0, 7.0, 8.0, 13.0, 12.0, 10.0, 5.0, 4.0], dtype="f4")
 100         self.subset_peak_str = "chrom:chrY      start:100       end:170 name:peak_1     score:13        summit:105\n"
 101
 102         self.t = sum([(x[2]-x[1]) * x[4] for x in self.input_regions])
 103
 104     def test_add_frag(self):
 105         pe = PETrackII()
 106         for (c, l, r, b, C) in self.input_regions:
 107             pe.add_loc(c, l, r, b, C)
 108         pe.finalize()
 109         # roughly check the numbers...
 110         self.assertEqual(pe.total, 22)
 111         self.assertEqual(pe.length, self.t)
 112
 113         # subset
 114         pe_subset = pe.subset(self.subset_barcodes)
 115         # roughly check the numbers...
 116         self.assertEqual(pe_subset.total, 14)
 117         self.assertEqual(pe_subset.length, 1305)
 118
 119     def test_pileup(self):
 120         pe = PETrackII()
 121         for (c, l, r, b, C) in self.input_regions:
 122             pe.add_loc(c, l, r, b, C)
 123         pe.finalize()
 124         bdg = pe.pileup_bdg()
 125         d = bdg.get_data_by_chr(b'chrY')  # (p, v) of ndarray
 126         np.testing.assert_array_equal(d[0], self.pileup_p)
 127         np.testing.assert_array_equal(d[1], self.pileup_v)
 128
 129         pe_subset = pe.subset(self.subset_barcodes)
 130         bdg = pe_subset.pileup_bdg()
 131         d = bdg.get_data_by_chr(b'chrY')  # (p, v) of ndarray
 132         np.testing.assert_array_equal(d[0], self.subset_pileup_p)
 133         np.testing.assert_array_equal(d[1], self.subset_pileup_v)
 134
 135     def test_pileup2(self):
 136         pe = PETrackII()
 137         for (c, l, r, b, C) in self.input_regions:
 138             pe.add_loc(c, l, r, b, C)
 139         pe.finalize()
 140         bdg = pe.pileup_bdg2()
 141         d = bdg.get_data_by_chr(b'chrY')  # (p, v) of ndarray
 142         np.testing.assert_array_equal(d['p'], self.pileup_p)
 143         np.testing.assert_array_equal(d['v'], self.pileup_v)
 144
 145         pe_subset = pe.subset(self.subset_barcodes)
 146         bdg = pe_subset.pileup_bdg2()
 147         d = bdg.get_data_by_chr(b'chrY')  # (p, v) of ndarray
 148         np.testing.assert_array_equal(d['p'], self.subset_pileup_p)
 149         np.testing.assert_array_equal(d['v'], self.subset_pileup_v)
 150
 151     def test_callpeak(self):
 152         pe = PETrackII()
 153         for (c, l, r, b, C) in self.input_regions:
 154             pe.add_loc(c, l, r, b, C)
 155         pe.finalize()
 156         bdg = pe.pileup_bdg()  # bedGraphTrackI object
 157         peaks = bdg.call_peaks(cutoff=10, min_length=20, max_gap=10)
 158         self.assertEqual(str(peaks), self.peak_str)
 159
 160         pe_subset = pe.subset(self.subset_barcodes)
 161         bdg = pe_subset.pileup_bdg()
 162         peaks = bdg.call_peaks(cutoff=10, min_length=20, max_gap=10)
 163         self.assertEqual(str(peaks), self.subset_peak_str)
 164
 165     def test_callpeak2(self):
 166         pe = PETrackII()
 167         for (c, l, r, b, C) in self.input_regions:
 168             pe.add_loc(c, l, r, b, C)
 169         pe.finalize()
 170         bdg = pe.pileup_bdg2()  # bedGraphTrackII object
 171         peaks = bdg.call_peaks(cutoff=10, min_length=20, max_gap=10)
 172         self.assertEqual(str(peaks), self.peak_str)
 173
 174         pe_subset = pe.subset(self.subset_barcodes)
 175         bdg = pe_subset.pileup_bdg2()
 176         peaks = bdg.call_peaks(cutoff=10, min_length=20, max_gap=10)
 177         self.assertEqual(str(peaks), self.subset_peak_str)