Merge pull request #200 from dib-lab/fix/query-rc

Fix problem with toggling query sequence strand

kevlar/call.py

@@ -376,6 +376,10 @@ def call(targetlist, querylist, match=1, mismatch=2, gapopen=5, gapextend=0,
                 query.sequence = kevlar.revcom(query.sequence)
             for varcall in make_call(besttarget, query, bestcigar, ksize):
                 yield varcall
+            if bestorientation == -1:
+                # Change it back!
+                # There's a better way to do this, but this works for now.
+                query.sequence = kevlar.revcom(query.sequence)
 
 
 def main(args):

kevlar/tests/data/multibestrc.contig.fa

@@ -0,0 +1,2 @@
+>contig1
+TCAGGCGGGTCCCCAAAGGACACACAAGCAAAACCCTCTGCCCAAGAGGGGTCATCCCAGGGCAATGGCTGGGGCTCAGGCCCAGCCTCATGGGCAGACCGGGCCAGGACCCGACTTGAGAGGGCTCAGGGAAGCCTCAAGCCCTGGGCAAGCCCCTCTCTCCAGGAGCCACATCCCCACTCAAATGAGTGCC

kevlar/tests/data/multibestrc.gdna.fa

@@ -0,0 +1,22 @@
+>17_5040570-5040785
+GCCCACAGGAGGCTCAGGCGGGTCCCCGAAGGACACACAAGCAAAACCCTCTGCCCAAGAGGGGTCATCCCAGGGCAATGGCTGGGACTCAGGTCCAGCCTCATGGGCAGACTGGGTCAGGACCCAACTTGAGAGGGCTCAGGGAAGCCTCAAGCCCTGGGCAAGCCCCTCTCTCCAGGAGCCACATACCCACTCAAATGAGTACCCCCCCATGA
+>17_18303281-18303424
+AGCCCACAGGAGGCTCAGGCGGGTCCCCGAAGGACACACAAGCAAAACCCTCTGCCCAAGAGGGGTCATCCCAGGGCAATGGCTGGGGCTCAGGCCCAGCCTCATGGGCAGACTGGGCCAGGACCCAACTTGAGAGGGCTCAG
+>17_20447649-20447781
+CCCTCTCAAGTTGGGTCCTGGCCCAGTCTGCCCATGAGGCTGGGCCTGAGCCCCAGCCATTGCCCTGGGATGACCCCTCTTGGGCAGAGGGTTTTGCTTGTGTGTCCTTCGGGGACCCGCCTGAGCCTCCTG
+>17_34495668-34495912
+AAGGTCTTGAAGCTCCTCATGGGGGGCACTCATTTGAGTGGGGATGTGGCTCCTGGAGAGAGGGGCTTGCCCAGGGCTTGAGGCTTCCCTGAGCCCTCTCAAGTCGGGTCCTGGCCCAGTCTGCCCATGAGGCTGGGCCTGAGCCCCAGCCATTGCCCTGGGATGACCCCTCTTGGGCAGAGGGTTTTGCTTGTGTGTCCTTTGGGGACCCGCCTGAGCCTCCTGTGGGCTGGGAGTGAGCCAG
+>17_34583711-34583957
+CAAGGTCTTGAAGCTCCTCATGGGGGGCACTCATTTGAGTGGGGATGTGGCTCCTGGAGAGAGGGGCTTGCCCAGGGCTTGAGGCTTCCCTGAGCCCTCTCAAGTCGGGTCCTGGCCCAGTCTGCCCATGAGGCTGGGCCTGAGCCCCAGCCATTGCCCTGGGATGACCCCTCTTGGGCAGAGGGTTTTGCTTGTGTGTCCTTTGGGGACCCGCCTGAGCCTTCTGTGGGCTGGGAGTGAGCCAGA
+>17_34748730-34748868
+ACAAGGTCTTGAAGCTCCTCATGGGGGGCACTCATTTGAGTGGGGATGTGGCTCCTGGAGAGAGGGGCTTGCCCAGGGCTTGAGGCTTCCCTGAGCCCTCTCAAGTCGGGTCCTGGCCCAGTCTGCCCATGAGGCTGG
+>17_34799780-34799928
+GACAAGGTCTTGAAGCTCCTCATGGGGGGCACTCATTTGAGTGGGGATGTGGCTCCTGGAGAGAGGGGCTTGCCCAGGGCTTGAGGCTTCCCTGAGCCCTCTCAAGTCGGGTCCTGGCCCAGTCTGCCCATGAGGCTGGGCCTGAGCC
+>17_36292137-36292277
+GGGCTCAGGCCCAGCCTCATGGGCAGACTGGGCCAGGACCCGACTTGAGAGGGCTCAGGGAAGCCTCAAGCCCTGGGCAAGCCCCTCTCTCCAGGAGCCACATCCCCACTCAAATGAGTGCCCCCCATGAGGAGCTTCAA
+>17_36340338-36340486
+AGACAAGGTCTTGAAGCTCCTCATGGGGGGCACTCATTTGAGTGGGGATGTGGCTCCTGGAGAGAGGGGCTTGCCCAGGGCTTGAGGCTTCCCTGAGCCCTCTCAAGTCGGGTCCTGGCCCAGTCTGCCCATGAGGCTGGGCCTGAGC
+>17_58088156-58088408
+CAGACAAGGTCTTGAAGCTCCTCATGGGAGGCACTCATTTGAGTGGGGATGTGGCTCCTGGAGAGAGGGGCTTGCCCAGGGCTTGAGGCTTCCCTGAGCCCTCTCAAGTCGGGTCCTGGCCCAGTCTGCCCATGAGGCTGGGCCTGAGCCCCAGCCATTGCCCTGGGATGACCCCTCTTGGGCAGAGGGTTTTGCTTGTGTGTCCTTTGGGGACCCGCCTGAGCCTCCTGTGGGCTGGGAGTGAGCCAGACC
+>17_60344730-60344982
+TCAGACAAGGTCTTGAAGCTCCTCATGGGGGGCACTCATTTGAGTGGGGATGTGGCTCCTGGAGAGAGGGGCTTGCCCAGGGCTTGAGGCTTCCCTGAGCCCTCTCAAGTCGGGTCCTGGCCCAGTCTGCCCATGAGGCTGGGCCTGAGCCCCAGCCATTGCCCTGGGATGACCCCTCTTGGGCAGAGGGTTTTGCTTGTGTGTCCTTTGGGGACCCGCCTGAGCCTCCTGTGGGCTGGGAGTGAGCCAGAC

kevlar/tests/test_call.py

@@ -288,3 +288,23 @@ def test_perfect_match():
     assert calls[0].seqid == 'chr99'
     assert calls[0].position == 2899377
     assert calls[0].info['NC'] == 'perfectmatch'
+
+
+def test_multibest_revcom():
+    contigfile = data_file('multibestrc.contig.fa')
+    contigstream = kevlar.parse_augmented_fastx(kevlar.open(contigfile, 'r'))
+    contigs = list(contigstream)
+
+    gdnafile = data_file('multibestrc.gdna.fa')
+    targets = list(khmer.ReadParser(gdnafile))
+
+    calls = list(call(targets, contigs))
+    assert len(calls) == 4
+
+    coordinates = [c.position + 1 for c in calls]
+    assert coordinates == [34495786, 34583830, 58088279, 60344854]
+    for c in calls:
+        assert c._refr == 'A'
+        assert c._alt == 'G'
+        assert c.window == ('CCTGAGCCCTCTCAAGTCGGGTCCTGGCCCGGTCTGCCCATGAGGCTGG'
+                            'GCCTGAGCCCCA')