Merge pull request #6 from yhoogstrate/2_1_0

2 1 0

Author: yhoogstrate

Changelog

@@ -1,3 +1,6 @@
+11-sept-2018:  Youri Hoogstrate
+	* v2.1.0 - Adds unit tests and can calculate ROC for Lorenz curves
+
 05-sept-2018:  Youri Hoogstrate
 	* v2.0.0 - does not require to use samtools, does not take all
 	  pysam.samtools.depth output into memory, does not export

README renamed to README.md

@@ -1,10 +1,10 @@
 bam-lorenz-coverage
 ===================
 
-Generate Lorenz and Cumulative Coverage plots on a bam file.
+Generate Lorenz plots and Coverage plots directly from BAM files
 
 Implemented in:
-    Python3 + Matplotlib + Pysam
+ * Python3 + Matplotlib + Pysam
 
 
 Installation:
@@ -27,10 +27,13 @@ Usage: bam-lorenz-coverage [OPTIONS] INPUT_ALIGNMENT_FILE
 Options:
   --version                  Show the version and exit.
   -l, --lorenz-table TEXT    Output table Lorenz-curve (for stdout use: -)
+  -x, --roc                  Output Lorenz-curve ROC to $lorenz_table.roc.txt
+                             [ requires --lorenz-table to be set to file ]
   -c, --coverage-table TEXT  Output table Coverage-graph (for stdout use: -)
   -L, --lorenz-svg TEXT      Output figure Lorenz-curve (SVG).
   -C, --coverage-svg TEXT    Output figure Coverage-graph (SVG).
   --help                     Show this message and exit.
+
 ```
 
 The lowercase arguments (-l, -c) allow extraction of the raw data tables for custom plotting. The uppercase arguments (-L, -C) directly generate a plot. The implemented plot only contains one sample per plot. For multi-sample plots, use the column tables and your imagination.

bin/bam-lorenz-coverage

@@ -12,10 +12,11 @@ from blc.blc import bamlorenzcoverage
 @click.version_option(blc.__version__ + "\n\n" + blc.__license_notice__ + "\n\nCopyright (C) 2018  " + blc.__author__ + ".\n\nFor more info please visit:\n" + blc.__homepage__)
 @click.argument('input_alignment_file', type=click.Path(exists=True))
 @click.option('-l', '--lorenz-table', nargs=1, help='Output table Lorenz-curve (for stdout use: -)')
+@click.option('-x', '--roc', is_flag=True, help='Output Lorenz-curve ROC to $lorenz_table.roc.txt\n[ requires --lorenz-table to be set to file ]')
 @click.option('-c', '--coverage-table', nargs=1, help='Output table Coverage-graph (for stdout use: -)')
 @click.option('-L', '--lorenz-svg', help='Output figure Lorenz-curve (SVG).')
 @click.option('-C', '--coverage-svg', help='Output figure Coverage-graph (SVG).')
-def CLI(lorenz_table, coverage_table, lorenz_svg, coverage_svg, input_alignment_file):
+def CLI(lorenz_table, roc, coverage_table, lorenz_svg, coverage_svg, input_alignment_file):
     b = bamlorenzcoverage()
     idx_observed = b.bam_file_to_idx(input_alignment_file)
 
@@ -41,6 +42,9 @@ def CLI(lorenz_table, coverage_table, lorenz_svg, coverage_svg, input_alignment_
             else:
                 with open(lorenz_table, 'w') as fh:
                     b.export_lorenz_curves(lorenz_curves, fh)
+                if roc:
+                    with open(lorenz_table + '.roc.txt', 'w') as fh:
+                        fh.write(str(lorenz_curves['roc']) + "\n")
 
         if lorenz_svg:
             b.export_lorenz_plot(lorenz_curves, lorenz_svg)

blc/__init__.py

@@ -5,7 +5,7 @@
 """[License: GNU General Public License v3 (GPLv3)]
 """
 
-__version_info__ = ('2', '0', '0')
+__version_info__ = ('2', '1', '0')
 __version__ = '.'.join(__version_info__) if (len(__version_info__) == 3) else '.'.join(__version_info__[0:3]) + "-" + __version_info__[3]
 __author__ = 'Youri Hoogstrate'
 __homepage__ = 'https://github.com/yhoogstrate/bam-lorenz-coverage'

blc/blc.py

@@ -248,22 +248,37 @@ def estimate_lorenz_curves(self, idx_observed):
         total_covered_positions_of_genome = cumu_frac_of_genome
         total_sequenced_bases = cumulative_sequenced_bases
 
+        previous = [0, 0]
+        top = 0
+        denom = total_sequenced_bases * total_covered_positions_of_genome * 2
+
         lorenz_curves = {'fraction_reads': [0.0], 'fraction_genome': [0.0]}
         for depth in sorted(lorenz_idx, reverse=True):
             if depth != 0:
                 lorenz_curves['fraction_reads'].append(round(1.0 * lorenz_idx[depth][1] / total_sequenced_bases, 4))
                 lorenz_curves['fraction_genome'].append(round(1.0 * lorenz_idx[depth][0] / total_covered_positions_of_genome, 4))
 
+                top += (lorenz_idx[depth][1] - previous[1]) * (lorenz_idx[depth][0] - previous[0])
+                top += (lorenz_idx[depth][1] - previous[1]) * (previous[0]) * 2
+
+                previous = lorenz_idx[depth]
+
+        roc = 1.0 * top / denom
+
+        lorenz_curves['roc'] = roc
         return lorenz_curves
 
     def export_lorenz_curves(self, lorenz_curves, output_stream):
         output_stream.write("X-fraction-sequenced-bases\tY-fraction-genome-covered\n")
         for n in range(len(lorenz_curves['fraction_reads'])):
             output_stream.write(str(lorenz_curves['fraction_reads'][n]) + "\t" + str(lorenz_curves['fraction_genome'][n]) + "\n")
 
-    def export_lorenz_plot(self, lorenz_curves, output_file):
+    def export_lorenz_plot(self, lorenz_curves, output_file, sign_digits=3):
         plt.plot([0.0, 1.0], [0.0, 1.0], 'k--')
         plt.plot(lorenz_curves['fraction_reads'], lorenz_curves['fraction_genome'], '-bo')
+        str_roc = round(lorenz_curves['roc'], sign_digits)
+        a = '{0:.' + str(sign_digits) + 'f}'
+        plt.text(0.0, 0.95, 'ROC=' + a.format(str_roc), fontsize=14)
         plt.xlabel('Fraction sequenced bases')
         plt.ylabel('Fraction covered genome')
         plt.savefig(output_file)

requirements.txt

@@ -1,5 +1,5 @@
-click
 numpy
+click
 pysam
 tqdm
 flake8

tests/blc/test_blc_011.sam

@@ -0,0 +1,5 @@
+@HD	VN:1.4
+@SQ	SN:chr1	LN:14
+@CO	user command line: test-data
+read_01	0	chr1	4	100	5M	*	0	0	NNNNN	BCCBC	NH:i:1	HI:i:1	AS:i:232	nM:i:0
+read_02	0	chr1	7	100	5M	*	0	0	NNNNN	BCCBC	NH:i:1	HI:i:1	AS:i:232	nM:i:0

tests/test_class_blc.py

@@ -146,7 +146,7 @@ def test_008_lorenz_01(self):
 
         # print(idx, file=sys.stderr)
         # denote that it only considers the (size of the) sequences described in the SAM header
-        self.assertDictEqual(lc, {'fraction_genome': [0.0, 1.0], 'fraction_reads': [0.0, 1.0]})
+        self.assertDictEqual(lc, {'fraction_genome': [0.0, 1.0], 'fraction_reads': [0.0, 1.0], 'roc': 0.5})
 
     def test_009_lorenz_02(self):
         # everything covered, is at least covered even densely
@@ -161,7 +161,7 @@ def test_009_lorenz_02(self):
 
         # print(idx, file=sys.stderr)
         # denote that it only considers the (size of the) sequences described in the SAM header
-        self.assertDictEqual(lc, {'fraction_genome': [0.0, 1.0], 'fraction_reads': [0.0, 1.0]})
+        self.assertDictEqual(lc, {'fraction_genome': [0.0, 1.0], 'fraction_reads': [0.0, 1.0], 'roc': 0.5})
 
     def test_010_lorenz_03(self):
         # everything covered, is at least covered even densely
@@ -176,7 +176,33 @@ def test_010_lorenz_03(self):
 
         # print(idx, file=sys.stderr)
         # denote that it only considers the (size of the) sequences described in the SAM header
-        self.assertDictEqual(lc, {'fraction_genome': [0.0, 1.0 * 2 / 8, 1.0], 'fraction_reads': [0.0, 1.0 * 4 / 10, 1.0]})
+        self.assertDictEqual(lc, {'fraction_genome': [0.0, 1.0 * 2 / 8, 1.0], 'fraction_reads': [0.0, 1.0 * 4 / 10, 1.0], 'roc': 0.425})
+
+    def test_011_lorenz_03(self):
+        # everything covered, is at least covered even densely
+        #             x x x x x
+        #       x x x x x
+        # - - - - - - - - - - - - - -
+
+        test_id = 'blc_011'
+
+        input_file_sam = TEST_DIR + "test_" + test_id + ".sam"
+        input_file_bam = T_TEST_DIR + "test_" + test_id + ".bam"
+
+        sam_to_sorted_bam(input_file_sam, input_file_bam)
+
+        b = bamlorenzcoverage()
+        idx = b.bam_file_to_idx(input_file_bam)
+
+        # print(idx, file=sys.stderr)
+        # denote that it only considers the (size of the) sequences described in the SAM header
+        self.assertDictEqual(idx, {0: 6, 1: 6, 2: 2})
+
+        lc = b.estimate_lorenz_curves(idx)
+
+        # print(idx, file=sys.stderr)
+        # denote that it only considers the (size of the) sequences described in the SAM header
+        self.assertDictEqual(lc, {'fraction_genome': [0.0, 1.0 * 2 / 8, 1.0], 'fraction_reads': [0.0, 1.0 * 4 / 10, 1.0], 'roc': 0.425})
 
 
 if __name__ == '__main__':