tests.py

from unittest import TestCase, main
from io import StringIO

from genomediff import Metadata, GenomeDiff
from genomediff.parser import GenomeDiffParser
from genomediff.records import (Record, TYPE_SPECIFIC_FIELDS,
                                PreservedFloat, PreservedInt)
from genomediff.schema import check_field


class ParserTestCase(TestCase):
    def test_parse(self):
        file = StringIO("""
#=GENOME_DIFF	1.0
#=AUTHOR test
SNP	1	23423	NC_000913	223	A	gene_name=mhpE
RA	2		NC_000913	223	0	G	A	frequency=0.1366
                        """.strip())
        p = GenomeDiffParser(fsock=file)
        self.assertEqual([
                             Metadata('GENOME_DIFF', '1.0'),
                             Metadata('AUTHOR', 'test'),
                             Record('SNP', 1, parent_ids=[23423], new_seq='A', seq_id='NC_000913', position=223, gene_name='mhpE'),
                             Record('RA', 2, new_base='A', frequency=0.1366, position=223, seq_id='NC_000913',
                                    insert_position=0,
                                    ref_base='G')],
                         list(p)
        )

    def test_parse_dot_missing_parent_ids(self):
        file = StringIO("""
#=GENOME_DIFF	1.0
#=AUTHOR test
SNP	1	23423	NC_000913	223	A	gene_name=mhpE
RA	2	.	NC_000913	223	0	G	A	frequency=0.1366
                        """.strip())
        p = GenomeDiffParser(fsock=file)
        self.assertEqual([
                             Metadata('GENOME_DIFF', '1.0'),
                             Metadata('AUTHOR', 'test'),
                             Record('SNP', 1, parent_ids=[23423], new_seq='A', seq_id='NC_000913', position=223, gene_name='mhpE'),
                             Record('RA', 2, new_base='A', frequency=0.1366, position=223, seq_id='NC_000913',
                                    insert_position=0,
                                    ref_base='G')],
                         list(p)
        )


class GenomeDiffTestCase(TestCase):
    def test_document(self):
        file = StringIO("""
#=GENOME_DIFF	1.0
#=AUTHOR test
SNP	1	23423	NC_000913	223	A
RA	2		NC_000913	223	0	G	A
                        """.strip())

        document = GenomeDiff.read(file)

        self.assertEqual({'AUTHOR': 'test', 'GENOME_DIFF': '1.0'}, document.metadata)

        snp_record = Record('SNP', 1, document, [23423], seq_id='NC_000913', new_seq='A', position=223)
        ra_record = Record('RA', 2, document, None, position=223, seq_id='NC_000913', insert_position=0, new_base='A',
                           ref_base='G')

        self.assertEqual([snp_record], document.mutations)
        self.assertEqual([ra_record], document.evidence)
        self.assertEqual(snp_record, document[1])
        self.assertEqual(ra_record, document[2])


class RecordTestCase(TestCase):
    def test_simple(self):
        snp_record = Record('SNP', 1, parent_ids=[23423], seq_id='NC_000913', new_seq='A', position=223, test='more')

        self.assertEqual('SNP', snp_record.type)
        self.assertEqual(1, snp_record.id)
        self.assertEqual('A', snp_record.new_seq)
        self.assertEqual('more', snp_record.test)


class ParentResolveTestCase(TestCase):
    def test_resolve(self):
        file = StringIO("""
#=GENOME_DIFF	1.0
#=AUTHOR test
SNP	1	2	NC_000913	223	A
RA	2		NC_000913	223	0	G	A
                        """.strip())
        document = GenomeDiff.read(file)
        self.assertEqual(document[1].parents, [document[2]])

class RecordComparisonTestCase(TestCase):
    def test_cmp1(self):
        file1 = StringIO("""
#=GENOME_DIFF	1.0
#=CREATED	20:02:17 23 Jan 2019
#=PROGRAM	breseq 0.33.2 
#=COMMAND	breseq -r LCA.gff3 sequence-data/DM0 evolved re-runs (Rohan)/ZDBp889_R1.fastq.gz sequence-data/DM0 evolved re-runs (Rohan)/ZDBp889_R2.fastq.gz sequence-data/ZDBp889_reads.fastq -o consensus/ZDBp889
#=REFSEQ	LCA.gff3
#=READSEQ	sequence-data/DM0 evolved re-runs (Rohan)/ZDBp889_R1.fastq.gz
#=READSEQ	sequence-data/DM0 evolved re-runs (Rohan)/ZDBp889_R2.fastq.gz
#=READSEQ	sequence-data/ZDBp889_reads.fastq
#=CONVERTED-BASES	644779377
#=CONVERTED-READS	14448149
#=INPUT-BASES	645034321
#=INPUT-READS	14455411
#=MAPPED-BASES	602854657
#=MAPPED-READS	13788351
SNP	1	34	REL606	72313	C
        """.strip())

        document1 = GenomeDiff.read(file1)
        
        file2 = StringIO("""
#=GENOME_DIFF	1.0
#=CREATED	16:49:49 23 Jan 2019
#=PROGRAM	breseq 0.33.2 
#=COMMAND	breseq -r LCA.gff3 sequence-data/DM0 evolved re-runs (Rohan)/ZDB67_R1.fastq.gz sequence-data/DM0 evolved re-runs (Rohan)/ZDB67_R2.fastq.gz -o consensus/ZDB67
#=REFSEQ	LCA.gff3
#=READSEQ	sequence-data/DM0 evolved re-runs (Rohan)/ZDB67_R1.fastq.gz
#=READSEQ	sequence-data/DM0 evolved re-runs (Rohan)/ZDB67_R2.fastq.gz
#=CONVERTED-BASES	114566968
#=CONVERTED-READS	419781
#=INPUT-BASES	114567554
#=INPUT-READS	419783
#=MAPPED-BASES	92472620
#=MAPPED-READS	339813
SNP	1	12	REL606	72313	C
        """.strip())

        document2 = GenomeDiff.read(file2)
        self.assertEqual(document1.mutations,document2.mutations)


    def test_cmp2(self):
        file1 = StringIO("""
#=GENOME_DIFF	1.0
SNP	1	12	REL606	72313	C	aa_new_seq=G	aa_position=92	aa_ref_seq=D	codon_new_seq=GGC	codon_number=92	codon_position=2	codon_ref_seq=GAC	gene_name=araA	gene_position=275	gene_product=L-arabinose isomerase	gene_strand=<	genes_overlapping=araA	locus_tag=ECB_00064	locus_tags_overlapping=ECB_00064	mutation_category=snp_nonsynonymous	position_end=72313	position_start=72313	snp_type=nonsynonymous	transl_table=11
        """.strip())

        document1 = GenomeDiff.read(file1)
        
        file2 = StringIO("""
#=GENOME_DIFF	1.0
SNP	1	34	REL606	72313	C	aa_new_seq=G	aa_position=92	aa_ref_seq=D	codon_new_seq=GGC	codon_number=92	codon_position=2	codon_ref_seq=GAC	gene_name=araA	gene_position=275	gene_product=L-arabinose isomerase	gene_strand=<	genes_overlapping=araA	locus_tag=ECB_00064	locus_tags_overlapping=ECB_00064	mutation_category=snp_nonsynonymous	position_end=72313	position_start=72313	snp_type=nonsynonymous	transl_table=11
        """.strip())

        document2 = GenomeDiff.read(file2)
        self.assertEqual(document1.mutations,document2.mutations)


class NewBreseqEntryTypesTestCase(TestCase):
    """Entry types breseq gained after TYPE_SPECIFIC_FIELDS was first written.

    Field names and order come from breseq's own line_specification in
    genome_diff_entry.cpp. Before these were added, reading a current breseq
    output.gd raised KeyError on the first SC line and no part of the file loaded.
    """

    def _read(self, *lines):
        return GenomeDiff.read(iter(('#=GENOME_DIFF\t1.0',) + lines))

    def test_sc(self):
        record = list(self._read('SC\t1\t.\tREL606\t511576\t-1').evidence)[0]
        self.assertEqual(
            {'seq_id': 'REL606', 'position': 511576, 'strand': -1},
            dict(record.attributes))

    def test_pd_has_dps_two_sided_shape(self):
        record = list(self._read(
            'PD\t1\t.\tREL606\t511576\t-1\tREL606\t511730\t1').evidence)[0]
        self.assertEqual('REL606', record.attributes['side_1_seq_id'])
        self.assertEqual(511730, record.attributes['side_2_position'])

    def test_an_unknown_type_still_reads(self):
        """A type this table has never heard of must not fail the whole file: its
        positional fields cannot be named, so they are kept in order."""
        record = list(self._read('ZZZ\t1\t.\tfoo\t42\tkey=value').mutations)[0]
        self.assertEqual(['foo', '42'], record.attributes['unknown_fields'])
        self.assertEqual('value', record.attributes['key'])

    def test_an_unknown_type_round_trips(self):
        record = list(self._read('ZZZ\t1\t.\tfoo\t42\tkey=value').mutations)[0]
        self.assertEqual('ZZZ\t1\t.\tfoo\t42\tkey=value', str(record))


class FieldGuardTestCase(TestCase):
    """Value guards mirrored from breseq's diff_entry_field_variable_types."""

    def test_each_guard_accepts_and_rejects(self):
        for entry_type, key, good, bad in [
                ('SNP', 'position', 5, 0),              # positive integer
                ('RA', 'insert_position', 0, -1),       # non-negative integer
                ('JC', 'overlap', -4, 'x'),             # any integer
                ('MOB', 'strand', -1, 2),               # strand
                ('SNP', 'new_seq', 'ACGTN', 'ACGTX')]:  # base sequence
            self.assertIsNone(check_field(entry_type, key, good),
                              '{} {}={} should be accepted'.format(entry_type, key, good))
            self.assertIsNotNone(check_field(entry_type, key, bad),
                                 '{} {}={} should be rejected'.format(entry_type, key, bad))

    def test_insert_position_floor_depends_on_entry_type(self):
        """The sole entry in breseq's diff_entry_load_field_variable_types."""
        self.assertIsNone(check_field('RA', 'insert_position', 0))
        self.assertIsNotNone(check_field('INS', 'insert_position', 0))
        self.assertIsNone(check_field('INS', 'insert_position', 1))

    def test_size_of_zero_is_allowed(self):
        """breseq lists 'size' twice and the later NonNegativeInteger entry wins."""
        self.assertIsNone(check_field('DEL', 'size', 0))
        self.assertIsNotNone(check_field('DEL', 'size', -1))

    def test_non_integral_values_are_rejected(self):
        self.assertIsNotNone(check_field('SNP', 'position', 1.5))
        self.assertIsNotNone(check_field('MOB', 'strand', True))

    def test_unguarded_and_absent_values_pass(self):
        self.assertIsNone(check_field('SNP', 'gene_name', 'araA'))
        self.assertIsNone(check_field('SNP', 'position', None))

    def test_table_matches_breseq_for_renamed_validation_types(self):
        self.assertEqual(('seq_id', 'enzyme'), TYPE_SPECIFIC_FIELDS['PFGE'])
        self.assertEqual(
            ('seq_id', 'primer_1_start', 'primer_1_end', 'primer_2_start', 'primer_2_end'),
            TYPE_SPECIFIC_FIELDS['RFLP'])


class AccessorAPITestCase(TestCase):
    def _snp(self):
        return Record('SNP', 1, None, None, seq_id='REL606', position=72313, new_seq='C')

    def test_get_and_set(self):
        record = self._snp()
        self.assertEqual(72313, record.get('position'))
        record.set('position', 100)
        self.assertEqual(100, record.get('position'))

    def test_get_missing(self):
        record = self._snp()
        self.assertEqual('fallback', record.get('nope', 'fallback'))
        self.assertRaises(KeyError, record.get, 'nope')

    def test_set_rejects_a_bad_value(self):
        record = self._snp()
        self.assertRaises(ValueError, record.set, 'position', 0)
        self.assertRaises(ValueError, record.set, 'new_seq', 'XYZ')
        self.assertEqual(72313, record.get('position'))  # unchanged

    def test_superseded_field_names_still_resolve(self):
        record = Record('RFLP', 1, None, None, seq_id='REL606', primer_1_start=1,
                        primer_1_end=2, primer_2_start=3, primer_2_end=4)
        self.assertEqual(1, record.primer1_start)
        self.assertEqual(2, record.get('primer1_end'))
        pfge = Record('PFGE', 2, None, None, seq_id='REL606', enzyme='EcoRI')
        self.assertEqual('EcoRI', pfge.restriction_enzyme)
        self.assertEqual('EcoRI', pfge.get('restriction_enzyme'))

    def test_set_stores_under_the_breseq_name(self):
        record = Record('PFGE', 1, None, None, seq_id='REL606')
        record.set('restriction_enzyme', 'EcoRI')
        self.assertEqual('EcoRI', record.attributes['enzyme'])
        self.assertNotIn('restriction_enzyme', record.attributes)

    def test_validate(self):
        self.assertEqual([], self._snp().validate())
        bad = Record('SNP', 1, None, None, seq_id='REL606', position=-1, new_seq='Q')
        self.assertEqual(2, len(bad.validate()))

    def test_unknown_fields_are_not_validated(self):
        record = Record('ZZZ', 1, None, None, unknown_fields=['foo', '42'])
        self.assertEqual([], record.validate())


class ReadLeniencyTestCase(TestCase):
    BAD = ('#=GENOME_DIFF\t1.0',
           'SNP\t1\t.\tREL606\t0\tC',
           'JC\t2\t.\tREL606\t100\t2\tREL606\t200\t1\t0')

    def test_a_guard_violation_does_not_lose_the_file(self):
        document = GenomeDiff.read(iter(self.BAD))
        self.assertEqual(1, len(document.mutations))
        self.assertEqual(1, len(document.evidence))
        self.assertEqual(2, len(document.parse_errors))
        # 1-based, counting the metadata header as line 1.
        self.assertEqual([2, 3], [e.line_number for e in document.parse_errors])

    def test_strict_raises(self):
        self.assertRaises(ValueError, GenomeDiff.read, iter(self.BAD), True)

    def test_a_valid_file_reports_no_errors(self):
        document = GenomeDiff.read(iter(
            ('#=GENOME_DIFF\t1.0', 'SNP\t1\t.\tREL606\t72313\tC')))
        self.assertEqual([], document.parse_errors)
        self.assertEqual([], document.validate())

    def test_a_violating_entry_still_round_trips(self):
        document = GenomeDiff.read(iter(self.BAD))
        self.assertEqual('SNP\t1\t.\tREL606\t0\tC', str(document.mutations[0]))

    def test_a_truncated_line_is_reported_not_raised(self):
        document = GenomeDiff.read(iter(('#=GENOME_DIFF\t1.0', 'DEL\t1\t.\tREL606')))
        record = document.mutations[0]
        self.assertEqual('REL606', record.get('seq_id'))
        self.assertIsNone(record.get('size'))
        self.assertEqual(1, len(document.parse_errors))
        self.assertIn('missing', document.parse_errors[0].message)

    def test_a_non_key_value_column_is_reported_not_raised(self):
        """gdtools VALIDATE reports this and carries on; so must we."""
        document = GenomeDiff.read(iter(
            ('#=GENOME_DIFF\t1.0', 'SNP\t1\t.\tchr\t100\tA\tgarbagefield')))
        self.assertEqual(1, len(document.mutations))
        self.assertEqual('A', document.mutations[0].get('new_seq'))
        self.assertEqual(1, len(document.parse_errors))
        self.assertIn('not a key=value pair', document.parse_errors[0].message)

    def test_a_trailing_tab_is_not_an_error(self):
        """breseq ignores the empty last column a trailing tab leaves behind."""
        document = GenomeDiff.read(iter(
            ('#=GENOME_DIFF\t1.0', 'MOB\t.\t.\tREL606-5\t2345\tIS150\t1\t3\t')))
        self.assertEqual([], document.parse_errors)
        self.assertEqual(3, document.mutations[0].get('duplication_size'))

    def test_document_validate_finds_what_read_collected(self):
        document = GenomeDiff.read(iter(self.BAD))
        self.assertEqual(len(document.parse_errors), len(document.validate()))


class MissingEntryIdTestCase(TestCase):
    """breseq writes '.' for an entry that has no id; requiring a number raised."""

    LINES = ('#=GENOME_DIFF\t1.0',
             'MOB\t.\t.\tREL606-5\t1000\tIS1\t-1\t3',
             'INS\t.\t.\tREL606-5\t1820930\tA',
             'SNP\t7\t.\tREL606-5\t100\tA')

    def test_entries_without_an_id_are_read(self):
        document = GenomeDiff.read(iter(self.LINES))
        self.assertEqual(3, len(document.mutations))
        self.assertEqual([None, None, 7], [r.id for r in document.mutations])
        self.assertEqual([], document.parse_errors)

    def test_an_absent_id_round_trips_as_a_dot(self):
        document = GenomeDiff.read(iter(self.LINES))
        self.assertEqual('MOB\t.\t.\tREL606-5\t1000\tIS1\t-1\t3',
                         str(document.mutations[0]))

    def test_only_real_ids_are_indexed(self):
        document = GenomeDiff.read(iter(self.LINES))
        self.assertEqual('REL606-5', document[7].get('seq_id'))
        self.assertRaises(KeyError, document.__getitem__, None)


class NumberFormattingTestCase(TestCase):
    """A number must write back as it was read, not as Python formats it."""

    def test_scientific_notation_survives_a_round_trip(self):
        line = ('RA\t1\t.\tNC_001416\t54\t0\tG\tA'
                '\tfrequency=8.39314286e-01\tpolymorphism_score=1.00000e+00')
        document = GenomeDiff.read(iter(('#=GENOME_DIFF\t1.0', line)))
        self.assertEqual(line, str(document.evidence[0]))

    def test_a_preserved_number_is_still_a_number(self):
        record = GenomeDiff.read(iter(
            ('#=GENOME_DIFF\t1.0',
             'RA\t1\t.\tX\t5\t0\tG\tA\tfrequency=8.39314286e-01'))).evidence[0]
        frequency = record.get('frequency')
        self.assertIsInstance(frequency, float)
        self.assertAlmostEqual(0.839314286, frequency)
        self.assertTrue(frequency < 1.0)
        self.assertAlmostEqual(1.678628572, frequency * 2)
        self.assertTrue(record.satisfies('frequency>=0.8'))

    def test_ordinary_numbers_are_not_wrapped(self):
        """Only values that would not format back need to carry their text."""
        convert = GenomeDiffParser._convert_value
        self.assertIs(int, type(convert('100')))
        self.assertIs(int, type(convert('-1')))
        self.assertIs(float, type(convert('0.1366')))
        self.assertIs(PreservedFloat, type(convert('8.39314286e-01')))
        self.assertIs(PreservedInt, type(convert('007')))

    def test_an_absent_value_writes_as_a_dot(self):
        """A key=value field of '.' must not come back out as 'None'."""
        line = 'RA\t1\t.\tNC_001416\t54\t0\tG\tA\tmajor_base=.\tminor_base=.'
        document = GenomeDiff.read(iter(('#=GENOME_DIFF\t1.0', line)))
        self.assertIsNone(document.evidence[0].get('major_base'))
        self.assertEqual(line, str(document.evidence[0]))
        self.assertNotIn('None', str(document.evidence[0]))


if __name__ == '__main__':
    main()