-Fixing off-by-one bugs #20

Author: kwade4

poplars/sequence_locator.py

@@ -75,7 +75,7 @@ def set_pos_from_cds(self, region_coords):
         if self.ncoords is not None:
             if self.region_name != '5\'LTR':
                 self.rel_pos['CDS'].append((self.ncoords[0] + 1 - region_coords[0]))
-                self.rel_pos['CDS'].append((self.ncoords[1] + 1 - region_coords[0]) + 1)
+                self.rel_pos['CDS'].append((self.ncoords[1] + 1 - region_coords[0]))
 
     def set_pos_from_gstart(self):
         self.rel_pos['gstart'] = self.ncoords
@@ -161,7 +161,7 @@ def get_overlap(self, region, coord_pair, base):
             end = min(local_pair[1] + 1, len(seq))
             if base == 'nucl':
                 overlap.append(seq[start - 1: end + 1])
-                overlap.append(self.local_to_global_index(region, [start, end], base))
+                overlap.append(self.local_to_global_index(region, [start, end + 1], base))
             else:
                 overlap.append(seq[start - 1: end - 1])
                 overlap.append(self.local_to_global_index(region, [start, end - 1], base))
@@ -455,7 +455,7 @@ def find_matches(base, ref_regions, match_coordinates):
 
     for coord in match_coordinates:
         start_aln = coord[0]
-        end_aln = coord[1] - 1  # -1 to account for match.end() in regex
+        end_aln = coord[1]
 
         for ref_region in ref_regions:
             if ref_region.region_name != "Complete":
@@ -504,9 +504,9 @@ def set_protein_equivalents(query_reg, ref_regions):
     for ref_reg in ref_regions:
         if ref_reg.region_name == query_reg.region_name and ref_reg.region_name not in non_coding:
             if ref_reg.codon_aln is not None and query_reg.pcoords is not None:
-                prot_equiv = ref_reg.codon_aln[query_reg.ncoords[0]: query_reg.ncoords[1]]
+                prot_equiv = ref_reg.codon_aln[query_reg.pcoords[0]: query_reg.pcoords[1]]
                 prot_equiv = re.sub('[-]', '', prot_equiv)
-                query_reg.set_aa_seq(prot_equiv)
+                query_reg.set_sequence('prot', prot_equiv)
 
     return prot_equiv
 
@@ -561,7 +561,7 @@ def output_retrieved_region(region, outfile=None):
 
         if region.rel_pos['gstart']:
             print("\t\tNucleotide position relative to genome start:\t{} --> {}"
-                  .format(region.rel_pos['gstart'][0] + 1, region.rel_pos['gstart'][1] + 1))
+                  .format(region.rel_pos['gstart'][0], region.rel_pos['gstart'][1]))
 
         if region.rel_pos['pstart']:
             print("\t\tAmino acid position relative to protein start:\t{} --> {}"
@@ -642,7 +642,7 @@ def output_overlap(overlap_regions, outfile=None):
 
             if region.rel_pos['gstart']:
                 print("\t\tNucleotide position relative to genome start:\t{} --> {}"
-                      .format(region.rel_pos['gstart'][0] + 1, region.rel_pos['gstart'][1] + 1))
+                      .format(region.rel_pos['gstart'][0], region.rel_pos['gstart'][1]))
 
             if region.rel_pos['pstart']:
                 print("\t\tAmino acid position relative to protein start:\t{} --> {}"
@@ -725,7 +725,6 @@ def retrieve(base, ref_regions, region, qstart=1, qend='end'):
             query_region = GenomeRegion(region)
 
             # Set local and global coordinates
-            # query_region.set_coords([qstart, qend], base)
             global_coords = GenomeRegion.local_to_global_index(ref_region, [qstart, qend], base)
             query_region.set_coords(global_coords, base)
 

poplars/tests/test_sequence_locator.py

@@ -205,13 +205,13 @@ def testGagFromAAStart(self):
 class TestMakeCodonAln(unittest.TestCase):
 
     def testSimpleUse(self):
-        region = GenomeRegion('p6', [2134, 2298], 'CTTCAGAGCAGACCAGAGCCAACAGCCCCACCAGAAGAGAGCTTCAGGTCTGGGGTAGAGACAACAAC'
+        region = GenomeRegion('p6', [2134, 2292], 'CTTCAGAGCAGACCAGAGCCAACAGCCCCACCAGAAGAGAGCTTCAGGTCTGGGGTAGAGACAACAAC'
                                                   'TCCCCCTCAGAAGCAGGAGCCGATAGACAAGGAACTGTATCCTTTAACTTCCCTCAGGTCACTCTTTG'
                                                   'GCAACGACCCCTCGTCACAATAA',
-                              [449, 500], 'LQSRPEPTAPPEESFRSGVETTTPPQKQEPIDKELYPLTSLRSLFGNDPSSQ')
+                                    [449, 500],   'LQSRPEPTAPPEESFRSGVETTTPPQKQEPIDKELYPLTSLRSLFGNDPSSQ*')
         region.make_codon_aln()
         expected = '-L--Q--S--R--P--E--P--T--A--P--P--E--E--S--F--R--S--G--V--E--T--T--T--P--P--Q-' \
-                   '-K--Q--E--P--I--D--K--E--L--Y--P--L--T--S--L--R--S--L--F--G--N--D--P--S--S--Q-'
+                   '-K--Q--E--P--I--D--K--E--L--Y--P--L--T--S--L--R--S--L--F--G--N--D--P--S--S--Q--*-'
         result = region.codon_aln
         self.assertEqual(expected, result)
 
@@ -1060,5 +1060,62 @@ def testInputAllSIV(self):
         self.assertEqual(expected_ref_prot, result_ref_prot)
 
 
+class TestSetNucleotideEquivalent(InputTestCase):
+
+    def testSimpleNuclEquiv(self):
+        configs = handle_args('hiv', 'prot', self.hiv_nt_seq_path, self.hiv_ncoords_path,
+                              self.hiv_aa_seq_path, self.hiv_pcoords_path)
+
+        ref_nt_seq = configs[0][0][1]
+        ref_aa_seq = configs[1]
+        nt_coords = configs[2]
+        aa_coords = configs[3]
+
+        with open(nt_coords) as ncoords, open(aa_coords) as pcoords:
+            ref_regions = set_regions('prot', ncoords, ref_nt_seq, pcoords, ref_aa_seq)
+        for region in ref_regions:
+            region.make_codon_aln()
+
+        query_region = GenomeRegion('Gag', [2133, 2292], 'GCTGAAGCAATGAGCCAAGTAACAAATTCAGCTACCATAATG',
+                                    [364, 377], 'AEAMSQVTNSATIM')
+
+        expected = 'GCTGAAGCAATGAGCCAAGTAACAAATTCAGCTACCATAATG'
+        result = set_nucleotide_equivalents(query_region, ref_regions)
+        self.assertEqual(expected, result)
+
+    def testWithGaps(self):
+        ref_regions = [GenomeRegion('test1', [1, 31], 'GGGGGCCCGGGTTTAAACCCGGGTTTAAATTTC', [1, 11], 'GGPGLNPGLNF')]
+        for region in ref_regions:
+            region.make_codon_aln()
+        q_region = GenomeRegion('test1', [12, 27],  'TTAAACCCGGGTTTA', [1, 5], 'LNPGL')
+        expected = 'TTAAACCCGGGTTTA'
+        result = set_nucleotide_equivalents(q_region, ref_regions)
+        self.assertEqual(expected, result)
+
+
+class TestSetProteinEquivalent(InputTestCase):
+
+    def testHIV(self):
+        configs = handle_args('hiv', 'nucl', self.hiv_nt_seq_path, self.hiv_ncoords_path,
+                              self.hiv_aa_seq_path, self.hiv_pcoords_path)
+
+        ref_nt_seq = configs[0][0][1]
+        ref_aa_seq = configs[1]
+        nt_coords = configs[2]
+        aa_coords = configs[3]
+
+        with open(nt_coords) as ncoords, open(aa_coords) as pcoords:
+            ref_regions = set_regions('nucl', ncoords, ref_nt_seq, pcoords, ref_aa_seq)
+
+        for reg in ref_regions:
+            reg.make_codon_aln()
+
+        query_region = GenomeRegion('Gag', [2133, 2292], 'GCTGAAGCAATGAGCCAAGTAACAAATTCAGCTACCATAATG',
+                                    [364, 377], 'AEAMSQVTNSATIM')
+        expected = 'AEAMSQVTNSATIM'
+        result = set_protein_equivalents(query_region, ref_regions)
+        self.assertEqual(expected, result)
+
+
 if __name__ == '__main__':
     unittest.main()