--- a
+++ b/tests/test_taps.py
@@ -0,0 +1,270 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+import unittest
+import pysam
+import os
+
+from singlecellmultiomics.molecule import TAPSNlaIIIMolecule, TAPS
+from singlecellmultiomics.fragment import NlaIIIFragment
+from singlecellmultiomics.utils import create_MD_tag
+
+from singlecellmultiomics.utils import complement
+
+
+class TestTAPs(unittest.TestCase):
+
+    def test_all(self):
+        temp_folder = 'data'
+
+        enable_ref_write=True
+
+
+        ref_path = f'{temp_folder}/ref.fa'
+        alignments_path = f'{temp_folder}/alignments.bam'
+
+        if not os.path.exists(temp_folder):
+            os.makedirs(temp_folder)
+        # Create reference bam file
+
+        refseq = 'TTAATCATGAAACCGTGGAGGCAAATCGGAGTGTAAGGCTTGACTGGATTCCTACGTTGCGTAGGTTCATGGGGGG'
+        if enable_ref_write:
+            with open(ref_path, 'w') as f:
+                f.write(f">chr1\n{refseq}\n>chr2\n{complement(refseq)}\n""")
+
+            # This command needs to finish, which is not working properly during testing
+            pysam.faidx(ref_path)
+
+        # CATG at base 5
+        # Create BAM file with NLA fragment:
+
+        alignments_path_unsorted = f'{alignments_path}.unsorted.bam'
+        with pysam.AlignmentFile(alignments_path_unsorted,'wb',reference_names=['chr1'],reference_lengths=[len(refseq)]) as bam:
+
+
+            ### Nla III mate pair example, containing 2 CpGs and 1 call on the wrong strand
+            read_A = pysam.AlignedSegment(bam.header)
+            read_A.reference_name = 'chr1'
+            read_A.reference_start = 5
+            # Before last A is a bogus G>A conversion to test strandness:
+            read_A.query_sequence = 'CATGAAACCGTGGAGGCAAATTGGAGTAT'
+            read_A.cigarstring = f'{len(read_A.query_sequence)}M'
+            read_A.qual = 'A'*len(read_A.query_sequence)
+            read_A.mapping_quality = 60
+            read_A.query_name = 'EX1_GA_CONV_2x_CpG_TAPS'
+            read_A.set_tag('SM', 'Cell_A')
+            read_A.is_read1 = True
+            read_A.is_read2 = False
+            read_A.set_tag('lh','TG')
+            # Set substitution tag:
+            read_A.set_tag('MD',
+                           create_MD_tag(
+                                   refseq[read_A.reference_start:read_A.reference_end], read_A.query_sequence))
+            read_A.is_paired = True
+            read_A.is_proper_pair = True
+
+
+            # Create a second read which is a mate of the previous
+            read_B = pysam.AlignedSegment(bam.header)
+            read_B.reference_name = 'chr1'
+            read_B.reference_start = 25
+            read_B.query_sequence = refseq[25:60].replace('TGT','TAT').replace('CG', 'TG')
+            read_B.cigarstring = f'{len(read_B.query_sequence)}M'
+            read_B.qual = 'A'*len(read_B.query_sequence)
+            read_B.mapping_quality = 60
+            read_B.is_read2 = True
+            read_B.is_read1 = False
+            read_B.is_reverse = True
+            read_B.query_name = 'EX1_GA_CONV_2x_CpG_TAPS'
+            read_B.set_tag('SM', 'Cell_A')
+            read_B.set_tag('lh','TG')
+            read_B.set_tag('MD',
+                       create_MD_tag(refseq[read_B.reference_start:read_B.reference_end],
+                                     read_B.query_sequence,
+                               ))
+            read_B.is_paired = True
+            read_B.is_proper_pair = True
+
+            read_A.next_reference_id = read_B.reference_id
+            read_A.next_reference_start = read_B.reference_start
+            read_B.next_reference_id = read_A.reference_id
+            read_B.next_reference_start = read_A.reference_start
+
+            read_A.mate_is_reverse = read_B.is_reverse
+            read_B.mate_is_reverse = read_A.is_reverse
+
+            bam.write(read_A)
+            bam.write(read_B)
+
+            ### Nla III mate pair example, dove tailed over random primer
+            # , containing 1 CpGs and one 1 call in the dove tail which should not be called
+            read_C = pysam.AlignedSegment(bam.header)
+            read_C.reference_name = 'chr1'
+            read_C.reference_start = 5
+            read_C.query_sequence = 'CATGAAACCGTGGAGGC'.replace('ACC','ATC').replace('AGGC','CGGT')
+            read_C.cigarstring = f'{len(read_C.query_sequence)}M'
+            read_C.qual = 'A'*len(read_C.query_sequence)
+            read_C.mapping_quality = 60
+            read_C.query_name = 'EX2_GA_DOVE'
+            read_C.set_tag('SM', 'Cell_A')
+            read_C.is_read1 = True
+            read_C.set_tag('lh','TG')
+            # Set substitution tag:
+            read_C.set_tag('MD',
+                           create_MD_tag(
+                                   refseq[read_C.reference_start:read_C.reference_end],
+                               read_C.query_sequence))
+
+            read_C.is_paired = True
+            read_C.is_proper_pair = True
+
+
+            # Create a second read which is a mate of the previous
+            read_D = pysam.AlignedSegment(bam.header)
+            read_D.reference_name = 'chr1'
+            read_D.reference_start = 10
+            read_D.query_sequence = refseq[10:15].replace('ACC','GTC')
+            read_D.cigarstring = f'{len(read_D.query_sequence)}M'
+            read_D.qual = 'A'*len(read_D.query_sequence)
+            read_D.mapping_quality = 60
+            read_D.is_read2 = True
+            read_D.is_read1 = False
+            read_D.is_reverse = True
+            read_D.query_name = 'EX2_GA_DOVE'
+            read_D.set_tag('SM', 'Cell_A')
+            read_D.set_tag('lh','TG')
+            read_D.set_tag('MD',
+                       create_MD_tag(refseq[read_D.reference_start:read_D.reference_end],
+                                     read_D.query_sequence,
+                               ))
+            read_D.is_paired = True
+            read_D.is_proper_pair = True
+
+            read_C.next_reference_id = read_D.reference_id
+            read_C.next_reference_start = read_D.reference_start
+            read_D.next_reference_id = read_C.reference_id
+            read_D.next_reference_start = read_C.reference_start
+
+            read_C.mate_is_reverse = read_D.is_reverse
+            read_D.mate_is_reverse = read_C.is_reverse
+
+            bam.write(read_C)
+            bam.write(read_D)
+
+            ########################################
+            # Reverse dovetailed (2 way) alignment #
+            ########################################
+
+            read_E = pysam.AlignedSegment(bam.header)
+            read_E.reference_name = 'chr1'
+            read_E.query_sequence = refseq[2:71].replace('CATGAA','CATAAA').replace('CGG','CAG')
+            read_E.reference_start = 71 - len(read_E.query_sequence)
+            read_E.cigarstring = f'{len(read_E.query_sequence)}M'
+            read_E.qual = 'A'*len(read_E.query_sequence)
+            read_E.mapping_quality = 60
+            read_E.query_name = 'EX2_GA_2xDOVE_rev'
+            read_E.set_tag('SM', 'Cell_A')
+            read_E.is_read2 = False
+            read_E.is_read1 = True
+            read_E.set_tag('lh','TG')
+            read_E.is_reverse = True
+            # Set substitution tag:
+            read_E.set_tag('MD',
+                           create_MD_tag(
+                                   refseq[read_E.reference_start:read_E.reference_end],
+                               read_E.query_sequence))
+            read_E.set_tag('ri','read_E')
+            read_E.is_paired = True
+            read_E.is_proper_pair = True
+
+
+            # Create a second read which is a mate of the previous
+            read_F = pysam.AlignedSegment(bam.header)
+            read_F.reference_name = 'chr1'
+            read_F.reference_start = 10
+            read_F.query_sequence = refseq[10:74].replace('CGG','CAG').replace('GGGG','GAGG')
+            read_F.cigarstring = f'{len(read_F.query_sequence)}M'
+            read_F.qual = 'A'*len(read_F.query_sequence)
+            read_F.mapping_quality = 60
+            read_F.is_read1 = False
+            read_F.is_read2 = True
+            read_F.is_reverse = False
+            read_F.query_name = 'EX2_GA_2xDOVE_rev'
+            read_F.set_tag('ri','read_F')
+            read_F.set_tag('SM', 'Cell_A')
+            read_F.set_tag('lh','TG')
+            read_F.set_tag('MD',
+                       create_MD_tag(refseq[read_F.reference_start:read_F.reference_end],
+                                     read_F.query_sequence,
+                               ))
+            read_F.is_paired = True
+            read_F.is_proper_pair = True
+
+
+            read_F.mate_is_reverse = read_E.is_reverse
+            read_E.mate_is_reverse = read_F.is_reverse
+
+            read_E.next_reference_id = read_F.reference_id
+            read_E.next_reference_start = read_F.reference_start
+            read_F.next_reference_id = read_E.reference_id
+            read_F.next_reference_start = read_E.reference_start
+
+            bam.write(read_E)
+            bam.write(read_F)
+
+
+        pysam.sort(alignments_path_unsorted, '-o', alignments_path)
+        pysam.index(alignments_path)
+
+        taps = TAPS()
+        with pysam.FastaFile(ref_path) as reference:
+
+            self.assertEqual(reference.fetch('chr1', 26, 26 + 3),'CGG')
+            molecule = TAPSNlaIIIMolecule(
+                NlaIIIFragment([read_A, read_B]),
+                reference=reference,
+                taps=taps,
+                taps_strand='F'
+            )
+            molecule.__finalise__()
+
+            calls = molecule.methylation_call_dict
+            print(calls)
+            print(calls[('chr1', 54)])
+            self.assertEqual( calls['chr1', 54]['context'], 'Z')
+            self.assertEqual( calls['chr1', 26]['context'], 'Z')
+            self.assertNotIn(  ('chr1', 26 + 6), calls)
+
+
+            molecule = TAPSNlaIIIMolecule(
+                NlaIIIFragment([read_E, read_F]),
+                reference =reference,
+                taps = taps,
+                taps_strand='F'
+            )
+            molecule.__finalise__()
+
+            # Test dove-tail detection:
+            self.assertNotIn( ('chr1', 71) , molecule.methylation_call_dict)
+            self.assertNotIn(('chr1', 8) , molecule.methylation_call_dict)
+
+
+
+            molecule = TAPSNlaIIIMolecule(
+                NlaIIIFragment([read_C, read_D]),
+                reference=reference,
+                taps=taps,
+                taps_strand='F'
+            )
+            molecule.__finalise__()
+
+            calls = molecule.methylation_call_dict
+
+            self.assertEqual(calls['chr1', 12]['context'], 'X')
+
+            # Check that dove tail is not included:
+            self.assertNotIn(('chr1', 21), calls)
+
+
+
+if __name__ == '__main__':
+    unittest.main()