#!/usr/bin/env python3

# -*- coding: utf-8 -*-

import unittest

import singlecellmultiomics.molecule

import singlecellmultiomics.fragment

import pysam

import pysamiterators.iterators

import os

from singlecellmultiomics.molecule import MoleculeIterator, CHICMolecule

from singlecellmultiomics.fragment import CHICFragment

"""

These tests check if the Molecule module is working correctly

"""

class TestMolecule(unittest.TestCase):

    def test_chic_cigar_dedup(self):

        i = 0

        with pysam.AlignmentFile('./data/chic_test_region.bam') as alignments:

            for molecule in MoleculeIterator(alignments,CHICMolecule, CHICFragment):

                i+=1

        self.assertEqual(i,1)

    def test_eq(self):

        def get_chic_read(header,  qname, contig ='chr1', start=100, sequence='ATCGGG', cigar=None, umi='CAT', sample='CELL_1', is_reverse=False, read1=True, paired=False,proper_pair=True ):

                read = pysam.AlignedSegment(header)

                read.set_tag('SM',sample) # The sample to which the sample belongs is extracted from the SM tag

                read.set_tag('RX',umi) # The UMI is extracted from the RX tag

                read.set_tag('MX','scCHIC')

                # By default the molecule assignment is done based on the mapping location of read 1:

                read.reference_name = contig

                read.reference_start = start

                read.query_name = qname

                read.query_sequence = sequence

                read.is_reverse = is_reverse

                read.cigarstring = f'{len(sequence)}M' if cigar is None else cigar

                if read1:

                    read.is_read1 = True

                    read.is_read2 = False

                else:

                    read.is_read1 = False

                    read.is_read2 = True

                if paired:

                    read.is_paired = True

                    read.is_proper_pair = proper_pair

                return read

        from singlecellmultiomics.molecule import Molecule

        from singlecellmultiomics.fragment import Fragment

        import pysam

        # Create sam file to write some reads to:

        with pysam.AlignmentFile('test.sam','w',reference_names=['chr1','chr2'],reference_lengths=[1000,1000]) as test_sam:

            read_A = get_chic_read(test_sam.header, 'read_A')

            read_B = get_chic_read(test_sam.header, 'read_B', start=102)

            read_C = get_chic_read(test_sam.header, 'read_C', start=110)

            # A reverse read

            read_D = get_chic_read(test_sam.header, 'read_D', start=110-5, is_reverse=True)

            # A different sample read

            read_E = get_chic_read(test_sam.header, 'read_E', start=110-5, is_reverse=True, sample='CELL_2')

            # A umi error read

            read_F = get_chic_read(test_sam.header, 'read_F', start=110-5, is_reverse=True, sample='CELL_2', umi='CAG')

            # A proper pair

            read_G_a = get_chic_read(test_sam.header, 'read_G', start=102,  paired=True)

            read_G_b = get_chic_read(test_sam.header, 'read_G', start=130, read1=False, paired=True, is_reverse=True)

            # A stale R2

            read_H = get_chic_read(test_sam.header, 'read_H', start=130,  paired=True, proper_pair=False, read1=False)

            # Read softclipped from start

            read_I = get_chic_read(test_sam.header, 'soft_clipped_read', start=104, cigar='2S4M')

            frag_A = Fragment([read_A],umi_hamming_distance=0,assignment_radius=2)

            frag_B = Fragment([read_B],umi_hamming_distance=0,assignment_radius=2)

            frag_C = Fragment([read_C],umi_hamming_distance=0,assignment_radius=2)

            frag_D = Fragment([read_D],umi_hamming_distance=0,assignment_radius=2)

            frag_E = Fragment([read_E],umi_hamming_distance=0,assignment_radius=2)

            frag_F = Fragment([read_F],umi_hamming_distance=0,assignment_radius=2)

            frag_G = Fragment([read_G_a, read_G_b],umi_hamming_distance=0,assignment_radius=2)

            self.assertTrue(frag_A==frag_B)

            self.assertTrue(frag_A==frag_G)

            self.assertFalse(frag_A==frag_C)

            self.assertFalse(frag_A==frag_D)

            self.assertFalse(frag_C==frag_D)

            self.assertFalse(frag_E==frag_D)

            frag_A = CHICFragment([read_A],umi_hamming_distance=0,assignment_radius=2)

            frag_B = CHICFragment([read_B],umi_hamming_distance=0,assignment_radius=2)

            frag_C = CHICFragment([read_C],umi_hamming_distance=0,assignment_radius=2)

            frag_D = CHICFragment([read_D],umi_hamming_distance=0,assignment_radius=2)

            frag_E = CHICFragment([read_E],umi_hamming_distance=1,assignment_radius=2)

            frag_F = CHICFragment([read_F],umi_hamming_distance=1,assignment_radius=2)

            frag_G = CHICFragment([read_G_a, read_G_b],umi_hamming_distance=0,assignment_radius=2)

            frag_I = CHICFragment([read_I],umi_hamming_distance=1,assignment_radius=2)

            frag_H = CHICFragment([None, read_H],umi_hamming_distance=1,assignment_radius=2)

            self.assertFalse(frag_H.is_valid())

            self.assertTrue(frag_A.is_valid())

            self.assertTrue(frag_B.is_valid())

            self.assertTrue(frag_A==frag_B)

            self.assertEqual(frag_G, frag_I)

            self.assertFalse(frag_A==frag_C)

            self.assertFalse(frag_A==frag_D)

            self.assertFalse(frag_C==frag_D)

            self.assertFalse(frag_E==frag_D)

            self.assertTrue(frag_A==frag_B)

            self.assertFalse(frag_A==frag_C)

            molecule_A = CHICMolecule([frag_A])

            print(molecule_A.site_location[1])

            self.assertEqual(molecule_A.site_location[1],98)

            self.assertTrue(molecule_A==frag_B)

            self.assertFalse(molecule_A==frag_C)

            # Test construction of molecule:

            self.assertTrue( molecule_A.add_fragment(frag_B) )

            self.assertEqual(len(molecule_A),2)

            # Frag C cannot be added:

            self.assertFalse( molecule_A.add_fragment(frag_C) )

            self.assertEqual(len(molecule_A),2)

            # Frag D cannot be added:

            self.assertFalse( molecule_A.add_fragment(frag_D) )

            self.assertEqual(len(molecule_A),2)

            # Test moving of site location by read which is aligned more to left:

            molecule = CHICMolecule([frag_B])

            self.assertEqual(molecule.site_location[1],100)

            self.assertTrue( molecule.add_fragment(frag_A) )

            self.assertEqual(molecule.site_location[1],98)

            # Test moving of site location by read which is aligned more to right: (Reverse strand)

            molecule = CHICMolecule([frag_B])

            self.assertEqual(molecule.site_location[1],100)

            self.assertTrue( molecule.add_fragment(frag_A) )

            self.assertEqual(molecule.site_location[1],98)

            # Test umi distance:

            molecule = CHICMolecule([frag_E])

            self.assertTrue( molecule.add_fragment(frag_F) )

            read_F.set_tag('RX','GGG')

            frag_G = CHICFragment([read_F],umi_hamming_distance=1,assignment_radius=2)

            self.assertFalse( molecule.add_fragment(frag_G) )

            # Perform iteration

            reads = [read_A, read_B, read_C, read_D, read_E, read_F, (read_G_a, read_G_b), (None,read_H), (read_I,None)]

            molecules = list(MoleculeIterator( reads,

                            molecule_class = CHICMolecule,

                            fragment_class = CHICFragment,

                            yield_invalid=True,

                            fragment_class_args={

                                'assignment_radius':2,

                                'umi_hamming_distance':1

                            }))

            recovered_reads = []

            for molecule in molecules:

                for read in molecule.iter_reads():

                    recovered_reads.append(read)

                    print(read.query_name)

            self.assertEqual(len(recovered_reads),10)

        os.remove('test.sam')

    def test_chic_nocigar_dedup(self):

        i = 0

        with pysam.AlignmentFile('./data/chic_test_region.bam') as alignments:

            for molecule in MoleculeIterator(alignments,CHICMolecule, CHICFragment,fragment_class_args={'no_umi_cigar_processing':True}):

                i+=1

        self.assertEqual(i,2)

    def test_a_pysam_iterators(self):

        """Test if the pysamiterators package yields the proper amount or mate pairs"""

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            for i,(R1,R2) in enumerate(pysamiterators.iterators.MatePairIterator(f)):

                pass

            self.assertEqual(i,224)

    def test_molecule_iterator_stability(self):

        """Test if the simplest molecule iterator works"""

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.Molecule,

            fragment_class=singlecellmultiomics.fragment.Fragment

            )

            pass

    def test_every_fragment_as_molecule(self):

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            for i,m in enumerate(singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.Molecule,

            fragment_class=singlecellmultiomics.fragment.Fragment,

            every_fragment_as_molecule=True

            )):

                pass

            self.assertEqual(i,340)

    def test_every_fragment_as_molecule_np_iterator(self):

        to_be_found = set()

        amount_of_r1s_to_be_found=0

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            for read in f:

                if read.is_read1 and not read.is_read2:

                    to_be_found.add(read.query_name)

                    amount_of_r1s_to_be_found+=1

        found=set()

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            frag_count=0

            for i,m in enumerate(singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.Molecule,

            fragment_class=singlecellmultiomics.fragment.Fragment,

            every_fragment_as_molecule=True,

            yield_invalid=True,

            iterator_class=pysamiterators.MatePairIteratorIncludingNonProper

            )):

                if m[0].has_R1():

                    frag_count+=1

                    found.add( m[0][0].query_name )

            diff = to_be_found.difference( found )

            tm = found.difference(to_be_found)

            print('missed:')

            print(diff)

            print('too much:')

            print(tm)

            self.assertEqual( len(diff), 0)

            self.assertEqual(frag_count,amount_of_r1s_to_be_found)

            self.assertEqual(frag_count,293)

    def test_Molecule_repr_stability(self):

        """Test if the molecule representation function is stable"""

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.Molecule,

            fragment_class=singlecellmultiomics.fragment.Fragment

            )

            for molecule in it:

                str(molecule)

    def test_iterable_molecule_iter(self):

        from singlecellmultiomics.molecule import MoleculeIterator

        from singlecellmultiomics.fragment import Fragment

        with pysam.AlignmentFile('test.sam','w',reference_names=['chr1','chr2'],reference_lengths=[1000,1000]) as test_sam:

            read_A = pysam.AlignedSegment(test_sam.header)

            read_A.set_tag('SM','CELL_1')

            read_A.set_tag('RX','CAT')

            read_A.reference_name = 'chr1'

            read_A.reference_start = 100

            read_A.query_sequence = 'ATCGGG'

            read_A.cigarstring = '6M'

            read_A.mapping_quality = 60

            read_B = pysam.AlignedSegment(test_sam.header)

            read_B.set_tag('SM','CELL_1')

            read_B.set_tag('RX','CAT')

            read_B.reference_name = 'chr1'

            read_B.reference_start = 100

            read_B.query_sequence = 'ATCGG'

            read_B.cigarstring = '5M'

            read_B.mapping_quality = 60

            read_C = pysam.AlignedSegment(test_sam.header)

            read_C.set_tag('SM','CELL_2')

            read_C.set_tag('RX','CAT')

            read_C.reference_name = 'chr1'

            read_C.reference_start = 100

            read_C.query_sequence = 'ATCGG'

            read_C.cigarstring = '5M'

            read_C.mapping_quality = 60

            reads = [  read_A,read_B,read_C ]

            mi = MoleculeIterator( reads , yield_invalid=True)

            molecules=[]

            for molecule in mi:

                molecules.append(molecule)

            self.assertEqual(len(molecules),2)

            self.assertEqual(max( (len(m) for m in molecules) ),2)

            self.assertEqual(min( (len(m) for m in molecules) ),1)

            # Test tags:

            a = molecules[0]

            a.write_tags()

            self.assertEqual(a[0][0].get_tag('TF') , 2)

        os.remove('test.sam')

    def test_NLA_Molecule_repr_stability(self):

        """Test if the molecule representation function is stable"""

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.NlaIIIMolecule,

            fragment_class=singlecellmultiomics.fragment.NlaIIIFragment

            )

            for molecule in it:

                str(molecule)

    def test_NLA_Molecule_ref_id(self):

        """Test if the molecule representation function is stable"""

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.NlaIIIMolecule,

            fragment_class=singlecellmultiomics.fragment.NlaIIIFragment

            )

            for molecule in it:

                self.assertEqual( molecule.get_a_reference_id(),0)

    def test_NLA_molecule_iterator_stability(self):

        """Test if the simplest molecule iterator works"""

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.NlaIIIMolecule,

            fragment_class=singlecellmultiomics.fragment.NlaIIIFragment

            )

            pass

    def test_consensus(self):

        """Test if a right consensus sequence can be produced from a noisy molecule"""

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.Molecule,

            fragment_class=singlecellmultiomics.fragment.NlaIIIFragment,

            fragment_class_args={

                'R1_primer_length':0,

                'R2_primer_length':6,

            }

            )

            for molecule in iter(it):

                #print(molecule.get_sample())

                if  molecule.get_sample()=='APKS3-P19-1-1_91':

                    break

            self.assertEqual(''.join( list(molecule.get_consensus().values()) ), 'CATGAGTTAGATATGGACTCTTCTTCAGACACTTTGTTTAAATTTTAAATTTTTTTCTGATTGCATATTACTAAAAATGTGTTATGAATATTTTCCATATCATTAAACATTCTTCTCAAGCATAACTTTAAATAACTGCATTATAGAAAATTTACGCTACTTTTGTTTTTGTTTTTTTTTTTTTTTTTTTACTATTATTAATAACACGGTGG')

    def test_fragment_sizes(self):

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            for molecule in singlecellmultiomics.molecule.MoleculeIterator(

                alignments=f,

                molecule_class=singlecellmultiomics.molecule.Molecule,

                fragment_class_args={

                    'R1_primer_length':4,

                    'R2_primer_length':6,

                },

                fragment_class=singlecellmultiomics.fragment.NlaIIIFragment):

                # This is a dovetailed molecule, both R1 and R2 overshoot into the adapter

                if molecule.sample=='APKS2-P18-2-1_66' and molecule.umi=='CGC':

                    # The non-dovetailed region of this molecule is

                    # 164834866-164834975

                    # This means a fragment size of 109,

                    # the 4 bases of the CATG are not counted.

                    # the 6 bases of the random primer are also not counted

                    # Resulting in a fragment size of 109 - 10 = 101

                    start, end = molecule[0].get_safe_span()

                    self.assertEqual(

                        abs(end-start)

                        , 107)

                    self.assertEqual(

                        molecule.get_safely_aligned_length()

                        , 107)

    def test_get_match_mismatch_frequency(self):

        """Test if the matches and mismatches of reads in a molecule are counted properly"""

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.Molecule,

            fragment_class_args={

                'R1_primer_length':0,

                'R2_primer_length':0,

            },

            fragment_class=singlecellmultiomics.fragment.NlaIIIFragment)

            for molecule in iter(it):

                #print(molecule.get_sample())

                if  molecule.get_sample()=='APKS3-P19-1-1_91':

                    break

        #print(molecule.get_match_mismatch_frequency())

        self.assertEqual( (628, 13), molecule.get_match_mismatch_frequency() )

    def test_get_consensus_gc_ratio(self):

        """Test if the gc ratio of a molecule is properly determined"""

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

            it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.Molecule,

            fragment_class=singlecellmultiomics.fragment.Fragment,

            fragment_class_args={

                'R1_primer_length':0,

                'R2_primer_length':0,

            }

            )

            for molecule in iter(it):

                #print(molecule.get_sample())

                if  molecule.get_sample()=='APKS3-P19-1-1_91':

                    break

        self.assertAlmostEqual(0.23113207547169812, molecule.get_consensus_gc_ratio() )

    def test_deduplication(self):

        f= pysam.AlignmentFile('./data/mini_nla_test.bam')

        total_frag_count = 0

        for molecule in singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.NlaIIIMolecule,

            fragment_class=singlecellmultiomics.fragment.NlaIIIFragment,

            fragment_class_args={'umi_hamming_distance':0},

            pooling_method=0,

            yield_invalid=True

        ):

            if 'APKS2-P14-1-1' in molecule.sample:

                for fragment in molecule:

                    is_dup = False

                    for read in fragment.reads:

                        if read is not None and read.is_duplicate:

                            is_dup = True

                    if not is_dup:

                        total_frag_count+=1

        self.assertEqual(total_frag_count,13)

    def _pool_test(self, pooling_method=0,hd=0):

        for sample in ['AP1-P22-1-1_318','APKS2-P8-2-2_52']:

            f= pysam.AlignmentFile('./data/mini_nla_test.bam')

            it = singlecellmultiomics.molecule.MoleculeIterator(

                alignments=f,

                molecule_class=singlecellmultiomics.molecule.Molecule,

                fragment_class=singlecellmultiomics.fragment.NlaIIIFragment,

                fragment_class_args={'umi_hamming_distance':hd},

                pooling_method=pooling_method

            )

            molecule_count = 0

            for molecule in iter(it):

                if molecule.get_sample()==sample:

                    molecule_count+=1

            if hd==0:

                # it has one fragment with a sequencing error in the UMI

                self.assertEqual(molecule_count,2)

            else:

                self.assertEqual(molecule_count,1)

    def test_molecule_pooling_vanilla_exact_umi(self):

        self._pool_test(0,0)

    def test_molecule_pooling_nlaIIIoptim_exact_umi(self):

        self._pool_test(1,0)

    def test_molecule_pooling_vanilla_umi_mismatch(self):

        self._pool_test(0,1)

    def test_molecule_pooling_nlaIIIoptim_umi_mismatch(self):

        self._pool_test(1,1)

    def test_max_associated_fragments(self):

        for i in range(1,3):

            with pysam.AlignmentFile('./data/mini_nla_test.bam') as f:

                for molecule in singlecellmultiomics.molecule.MoleculeIterator(

                    alignments=f,

                    molecule_class=singlecellmultiomics.molecule.Molecule,

                    fragment_class=singlecellmultiomics.fragment.NlaIIIFragment,

                    molecule_class_args={'max_associated_fragments':i},

                    fragment_class_args={'umi_hamming_distance':1}

                ):

                    self.assertTrue(len(molecule)<= i , f'{i}, {len(molecule)}')

    def test_rt_reaction_counting_HAMMING_1(self):

        # This is a dictionary containing molecules and the amount of rt reactions for location chr1:164834865

        f= pysam.AlignmentFile('./data/mini_nla_test.bam')

        it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.Molecule,

            fragment_class=singlecellmultiomics.fragment.NlaIIIFragment,

            fragment_class_args={'umi_hamming_distance':1}

        )

        hand_curated_truth = {

            # hd: 1:

            # hard example with N in random primer and N in UMI:

            'APKS2-P8-2-2_52':{'rt_count':3},

            # Simple examples:

            'APKS2-P18-1-1_318':{'rt_count':1},

            'APKS2-P18-1-1_369':{'rt_count':2},

            'APKS2-P18-2-1_66':{'rt_count':1},

            'APKS2-P18-2-1_76':{'rt_count':1},

            'APKS2-P18-2-1_76':{'rt_count':1},

        }

        obtained_rt_count = {}

        for molecule in it:

            site = molecule.get_cut_site()

            if site is not None and site[1]==164834865:

                obtained_rt_count[molecule.get_sample()]  = len( molecule.get_rt_reaction_fragment_sizes(max_N_distance=1) )

        # Validate:

        for sample, truth in hand_curated_truth.items():

            self.assertEqual( obtained_rt_count.get(sample,0),truth.get('rt_count',-1) )

    def test_rt_reaction_counting_HAMMING_0(self):

        # This is a dictionary containing molecules and the amount of rt reactions for location chr1:164834865

        f= pysam.AlignmentFile('./data/mini_nla_test.bam')

        it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.Molecule,

            fragment_class=singlecellmultiomics.fragment.NlaIIIFragment,

            fragment_class_args={'umi_hamming_distance':1}

        )

        hand_curated_truth = {

            # hd: 0:

            # hard example with N in random primer and N in UMI:

            'APKS2-P8-2-2_52':{'rt_count':3},

            # Simple examples:

            'APKS2-P18-1-1_318':{'rt_count':1},

            'APKS2-P18-1-1_369':{'rt_count':2},

            'APKS2-P18-2-1_66':{'rt_count':1},

            'APKS2-P18-2-1_76':{'rt_count':1},

            'APKS2-P18-2-1_76':{'rt_count':1},

        }

        obtained_rt_count = {}

        for molecule in it:

            site = molecule.get_cut_site()

            if site is not None and site[1]==164834865:

                obtained_rt_count[molecule.get_sample()]  = len( molecule.get_rt_reaction_fragment_sizes(max_N_distance=0) )

        # Validate:

        for sample, truth in hand_curated_truth.items():

            self.assertEqual( obtained_rt_count.get(sample,0),truth.get('rt_count',-1) )

    def test_feature_molecule(self):

        import singlecellmultiomics.features

        features = singlecellmultiomics.features.FeatureContainer()

        features.addFeature(chromosome='chr1',start=164835366, end=164835535, data=('test4'),name='test4')

        features.sort()

        f= pysam.AlignmentFile('./data/mini_nla_test.bam')

        it = singlecellmultiomics.molecule.MoleculeIterator(

            alignments=f,

            molecule_class=singlecellmultiomics.molecule.AnnotatedNLAIIIMolecule,

            molecule_class_args={'features':features},

            fragment_class=singlecellmultiomics.fragment.NlaIIIFragment,

            fragment_class_args={'umi_hamming_distance':1}

        )

        hit_count = 0

        for molecule in it:

            molecule.annotate()

            if len(molecule.hits):

                hit_count+=1

        self.assertEqual(hit_count,2)

    """

    def test_classification_consensus(self):

        with pysam.AlignmentFile('./data/mini_nla_test.bam') as f, pysam.AlignmentFile('./data/consensus_write_test.bam','wb',header=f.header) as target_bam:

            molecule_iterator =  singlecellmultiomics.molecule.MoleculeIterator(

                alignments=f,

                molecule_class=singlecellmultiomics.molecule.NlaIIIMolecule,

                fragment_class=singlecellmultiomics.fragment.NlaIIIFragment

            )

            classifier = singlecellmultiomics.molecule.train_consensus_model(

                        molecule_iterator,

                        mask_variants=None,

                        skip_already_covered_bases=False,

                        n_train=100)

            molecule_iterator =  singlecellmultiomics.molecule.MoleculeIterator(

                alignments=f,

                molecule_class=singlecellmultiomics.molecule.NlaIIIMolecule,

                fragment_class=singlecellmultiomics.fragment.NlaIIIFragment

            )

            for i,molecule in enumerate(molecule_iterator):

                read_name=f'consensus_{i}'

                reads = molecule.deduplicate_to_single_CIGAR_spaced(target_bam,

                read_name, classifier,reference=None)

                read = molecule.deduplicate_to_single(target_bam)

                consensus = molecule.get_consensus(classifier)

                """

if __name__ == '__main__':

    unittest.main()