new things

- fixes
- sumcov

Author: Dom Laetsch

blobtools

@@ -11,6 +11,7 @@ commands:
 
   comparecov    compare BlobDB cov(s) to additional cov file
   bam2cov       generate cov file from bam file
+  sumcov        sum coverage from multiple COV files
 
 -h --help    show this
 
@@ -40,5 +41,7 @@ if __name__ == '__main__':
         exit(call(['python', main_dir + 'bam2cov.py'] + argv))
     elif args['<command>'] == 'comparecov':
         exit(call(['python', main_dir + 'comparecov.py'] + argv))  
+    elif args['<command>'] == 'sumcov':
+        exit(call(['python', main_dir + 'sumcov.py'] + argv))  
     else:
         exit("%r is not a blobtools command. See 'blobtools -h'." % args['<command>'])

comparecov.py

@@ -179,7 +179,7 @@
         if (plotObj.title):
             plotObj.title = "%s.%s.%s" % (title, taxrule, cov_lib)
 
-        out_f = "%s.%s.%s.p%s.%s" % (title, hist_type, rank, max_group_plot, cov_lib)
+        out_f = "%s.%s.%s.p%s.%s.%s" % (title, hist_type, rank, max_group_plot, scale, cov_lib)
         if out_prefix:
             out_f = "%s.%s" % (out_prefix, out_f)
         #if catcolour_dict:

lib/BtLog.py

@@ -22,7 +22,7 @@ def error(message, *argv):
 def progress(iteration, steps, max_value):
     if int(iteration == max_value):
         sys.stdout.write('\r')
-        print "[PROGRESS]\t: %d%%" % (100)
+        print "[PROGRESS]\t: %d%%" % (100),
     else:
         if int(iteration) % int(steps) == 0:
             sys.stdout.write('\r')

lib/BtPlot.py

@@ -132,13 +132,17 @@ def set_format_scatterplot(axScatter, max_cov):
     axScatter.xaxis.labelpad = 20
     return axScatter
 
-def set_format_covplot(axScatter, x_max, y_max, x_label, y_label):
-    axScatter.set_xlabel(x_label, fontsize=35)
+def set_format_covplot(axScatter, x_max, y_max, x_label, y_label, scale, cov_lib):
+    axScatter.set_xlabel(x_label + " " + cov_lib, fontsize=35)
     axScatter.set_ylabel(y_label, fontsize=35)
     axScatter.grid(True, which="major", lw=2., color=WHITE, linestyle='-') 
     axScatter.set_axisbelow(True)
-    axScatter.set_xlim( (0.01, x_max+100) )
-    axScatter.set_ylim( (0.01, y_max+100) ) # This sets the max-Coverage so that all libraries + sum are at the same scale
+    if scale == 'log':
+        axScatter.set_xlim( (0.01, x_max+100) )
+        axScatter.set_ylim( (0.01, y_max+100) ) 
+    else:
+        axScatter.set_xlim( (-10, x_max + (x_max*0.1) ) )
+        axScatter.set_ylim( (-10, y_max + (y_max*0.1) ) ) 
     axScatter.xaxis.labelpad = 20
     axScatter.xaxis.labelpad = 20
     return axScatter
@@ -363,7 +367,7 @@ def plotScatterCov(self, cov_lib, cov_dict, info_flag, x_label, y_label, scale,
         # Setting up plots and axes
         plt.figure(1, figsize=(35,35), dpi=400)
         axScatter = plt.axes(rect_scatter, axisbg=BGGREY, yscale = scale, xscale = scale)
-        axScatter = set_format_covplot(axScatter, x_max, y_max, x_label, y_label)
+        axScatter = set_format_covplot(axScatter, x_max, y_max, x_label, y_label, scale, cov_lib)
         axHistx = plt.axes(rect_histx, axisbg=BGGREY)
         axHistx = set_format_hist_x_cov(axHistx, axScatter, scale)
         axHisty = plt.axes(rect_histy, axisbg=BGGREY)

sumcov.py

@@ -0,0 +1,63 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""usage: blobtools sumcov          -i FASTA -c COV... 
+                                    [-h|--help] 
+
+    Options:
+        -h --help                   show this
+
+        -i, --infile FASTA          FASTA file of assembly. Headers are split at whitespaces.  
+        -c, --cov COV               COV file(s) to sum up
+"""
+
+from __future__ import division
+from docopt import docopt
+import lib.BtCore as bt
+import lib.BtLog as BtLog
+import lib.BtIO as BtIO
+import lib.BtPlot as BtPlot
+from os.path import dirname, isfile
+
+def readFasta(infile):
+    header = ''
+    with open(infile) as fh: 
+        for l in fh: 
+            if l[0] == '>':
+                yield l[1:-1].split()[0]
+            else:
+                pass
+
+def parseFasta(infile):
+    fasta_order = []
+    fasta_dict = {}
+    for name in readFasta(infile):
+        fasta_order.append(name)
+        fasta_dict[name] = 0.0
+    return fasta_dict, fasta_order
+
+if __name__ == '__main__':
+    main_dir = dirname(__file__)
+    #print data_dir
+    args = docopt(__doc__)
+    assembly_f = args['--infile']
+    cov_fs = args['--cov']
+    
+    fasta_dict = {}
+    fasta_order = []
+    if not isfile(assembly_f):
+        BtLog.error('0', assembly_f)
+    else:
+        fasta_dict, fasta_order = parseFasta(assembly_f)
+    
+    for cov_f in cov_fs:
+        if not isfile(cov_f):
+            BtLog.error('0', cov_f)
+        else:
+            lib_cov_dict = BtPlot.parseCovFile(cov_f)
+            for name in fasta_order:
+                fasta_dict[name] = fasta_dict.get(name, 0.0) + lib_cov_dict[name]
+                    
+    
+    for name in fasta_order:
+        print "%s\t%s" % (name, fasta_dict[name])