Merge pull request #28 from martinghunt/bug_when_samtools_zero_depth

Bug when samtools zero depth
sanger-pathogens · May 12, 2015 · 2b6417d · 2b6417d
2 parents 9e961d9 + e1472bf
commit 2b6417d
Show file tree

Hide file tree

Showing 3 changed files with 19 additions and 13 deletions.
diff --git a/ariba/cluster.py b/ariba/cluster.py
@@ -136,8 +136,8 @@ def __init__(self,
 
 
     def _get_read_counts(self):
-        count1 = pyfastaq.tasks.count_sequences(self.reads1) 
-        count2 = pyfastaq.tasks.count_sequences(self.reads2) 
+        count1 = pyfastaq.tasks.count_sequences(self.reads1)
+        count2 = pyfastaq.tasks.count_sequences(self.reads2)
         if count1 == count2:
             return count1 + count2
         else:
@@ -160,7 +160,7 @@ def _get_total_alignment_score(self, gene_name):
             bowtie2=self.bowtie2_exe,
             bowtie2_preset=self.bowtie2_preset,
             verbose=self.verbose,
-        ) 
+        )
 
         score = mapping.get_total_alignment_score(tmp_bam)
         os.unlink(tmp_bam)
@@ -795,7 +795,7 @@ def _get_samtools_variant_positions(self):
         positions = [(t[0], int(t[1]) - 1) for t in positions]
         pyfastaq.utils.close(f)
         return positions
-        
+
 
     def _get_samtools_variants(self, positions=None):
         if positions is None:
@@ -870,11 +870,17 @@ def _make_report_lines(self):
                 if t is not None:
                     effect, new_bases = t
                     for v in variants:
-                        depths = self._get_assembly_read_depths(contig, v.qry_start) 
+                        depths = self._get_assembly_read_depths(contig, v.qry_start)
                         if depths is None:
-                            raise Error('Error getting depth info on contig "' + contig + '" at position ' + str(v.qry_start + 1) + ' from file ' + self.final_assembly_read_depths)
-                        else:
-                            ref_base, alt_base, ref_counts, alt_counts = depths
+                            # this happens with low coverage contigs. It can get assembled, but
+                            # there are some bases that do not have reads mapped to them.
+                            # If mummer called a variant at one of these, then we're looking
+                            # for read dpeth where there is none.
+                            print('Warning: could not get read depth info on contig "' + contig + '" at position ', str(v.qry_start + 1), 'from file', self.final_assembly_read_depths, file=sys.stderr)
+                            print(' - a variant was called at this position using nucmer, but there is no read depth (probably a mapping artifact)', file=sys.stderr)
+                            depths = ['.'] * 4
+
+                        ref_base, alt_base, ref_counts, alt_counts = depths
 
                         self.report_lines.append([
                             self.gene.id,
@@ -902,10 +908,10 @@ def _make_report_lines(self):
 
                         if contig in samtools_variants and v.qry_start in samtools_variants[contig]:
                             del samtools_variants[contig][v.qry_start]
-                            if len(samtools_variants[contig]) == 0: 
+                            if len(samtools_variants[contig]) == 0:
                                 del samtools_variants[contig]
 
-            if contig in samtools_variants: 
+            if contig in samtools_variants:
                 for pos in samtools_variants[contig]:
                     ref_base, alt_base, ref_counts, alt_counts = samtools_variants[contig][pos]
                     self.report_lines.append(
@@ -930,7 +936,7 @@ def _make_report_lines(self):
                         alt_counts
                       ]
                     )
-                        
+
         if len(self.report_lines) == 0:
             for contig in self.percent_identities:
                 self.report_lines.append([

diff --git a/ariba/common.py b/ariba/common.py
@@ -1,7 +1,7 @@
 import sys
 import subprocess
 
-version = '0.3.1'
+version = '0.3.2'
 
 def syscall(cmd, allow_fail=False, verbose=False):
     if verbose:

diff --git a/setup.py b/setup.py
@@ -7,7 +7,7 @@
 
 setup(
     name='ariba',
-    version='0.3.1',
+    version='0.3.2',
     description='ARIBA: Antibiotic Resistance Identification By Assembly',
     packages = find_packages(),
     author='Martin Hunt',