| Previous changeset 1:d1fdfaae5dbe (2014-11-21) Next changeset 3:57b3ea22aff3 (2020-08-11) |
|
Commit message:
v0.0.5 - Supports max coverage depth in recent samtools. Expects samtools 1.4.1 via Conda, not via Tool Shed. |
|
modified:
tools/coverage_stats/README.rst tools/coverage_stats/coverage_stats.py tools/coverage_stats/coverage_stats.xml |
|
added:
test-data/ex1.coverage_stats.md50.tabular |
|
removed:
tools/coverage_stats/tool_dependencies.xml |
| b |
| diff -r d1fdfaae5dbe -r 7254ece0c0ff test-data/ex1.coverage_stats.md50.tabular --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/ex1.coverage_stats.md50.tabular Thu May 11 12:16:10 2017 -0400 |
| b |
| @@ -0,0 +1,4 @@ +#identifer length mapped placed min_cov max_cov mean_cov +chr1 1575 1446 18 0 50 32.06 +chr2 1584 1789 17 0 50 38.70 +* 0 0 0 0 0 0.00 |
| b |
| diff -r d1fdfaae5dbe -r 7254ece0c0ff tools/coverage_stats/README.rst --- a/tools/coverage_stats/README.rst Fri Nov 21 09:43:58 2014 -0500 +++ b/tools/coverage_stats/README.rst Thu May 11 12:16:10 2017 -0400 |
| b |
| @@ -1,7 +1,7 @@ BAM coverage statistics using samtools idxstats and depth ========================================================= -This tool is copyright 2014 by Peter Cock, The James Hutton Institute +This tool is copyright 2014-2017 by Peter Cock, The James Hutton Institute (formerly SCRI, Scottish Crop Research Institute), UK. All rights reserved. See the licence text below. @@ -14,14 +14,15 @@ Automated Installation ====================== -This should be straightforward, Galaxy should automatically download and install -samtools 0.1.19 if required. +This should be straightforward, provided your Galaxy is setup to use Conda and +BioConda for dependencies - Galaxy should then automatically download and +install the expected version of samtools. Manual Installation =================== -This expects samtools to be on the ``$PATH``, and was tested using v0.1.19. +This expects samtools to be on the ``$PATH``, and was tested using v1.3.1. To install the wrapper copy or move the following files under the Galaxy tools folder, e.g. in a ``tools/coverage_stats`` folder: @@ -50,6 +51,20 @@ Version Changes ------- ---------------------------------------------------------------------- v0.0.1 - Initial public release +v0.0.2 - Cope with samtools' default depth limit using modified samtools, + see https://github.com/samtools/samtools/pull/322 +v0.0.3 - Cope with no coverage in final contigs. +v0.0.4 - Reorder XML elements (internal change only). + - Planemo for Tool Shed upload (``.shed.yml``, internal change only). +v0.0.5 - Expose new ``samtools depth -d ...`` argument added in samtools v1.3 + - Depends on samtools v1.4.1 via Conda, which IS NOT AVAILABLE as a + legacy Tool Shed package. + - Apply the new maximum depth parameter within the script to ensure + excess coverage is clear by getting the max coverage equal to the + max depth setting (the raw output from samtools is more fuzzy). + - Report total length and overall mean coverage in stdout. + - Use ``<command detect_errors="aggressive">`` (internal change only). + - Single quote command line arguments (internal change only). ======= ====================================================================== @@ -59,22 +74,27 @@ Development is on this GitHub repository: https://github.com/peterjc/pico_galaxy/tree/master/tools/coverage_stats -For making the "Galaxy Tool Shed" http://toolshed.g2.bx.psu.edu/ tarball use -the following command from the Galaxy root folder:: +For pushing a release to the test or main "Galaxy Tool Shed", use the following +Planemo commands (which requires you have set your Tool Shed access details in +``~/.planemo.yml`` and that you have access rights on the Tool Shed):: + + $ planemo shed_update -t testtoolshed --check_diff tools/coverage_stats/ + ... + +or:: - $ tar -czf coverage_stats.tar.gz tools/coverage_stats/README.rst tools/coverage_stats/coverage_stats.xml tools/coverage_stats/coverage_stats.py tools/coverage_stats/tool_dependencies.xml test-data/ex1.bam test-data/ex1.coverage_stats.tabular test-data/coverage_test.bam test-data/coverage_test.coverage_stats.tabular + $ planemo shed_update -t toolshed --check_diff tools/coverage_stats/ + ... -Check this worked:: +To just build and check the tar ball, use:: - $ tar -tzf coverage_stats.tar.gz + $ planemo shed_upload --tar_only tools/coverage_stats/ + ... + $ tar -tzf shed_upload.tar.gz tools/coverage_stats/README.rst tools/coverage_stats/coverage_stats.xml tools/coverage_stats/coverage_stats.py - tools/coverage_stats/tool_dependencies.xml - test-data/ex1.bam - test-data/ex1.coverage_stats.tabular - test-data/coverage_test.bam - test-data/coverage_test.coverage_stats.tabular + ... Licence (MIT) |
| b |
| diff -r d1fdfaae5dbe -r 7254ece0c0ff tools/coverage_stats/coverage_stats.py --- a/tools/coverage_stats/coverage_stats.py Fri Nov 21 09:43:58 2014 -0500 +++ b/tools/coverage_stats/coverage_stats.py Thu May 11 12:16:10 2017 -0400 |
| [ |
| b'@@ -6,39 +6,61 @@\n * Input BAI filename (via Galaxy metadata)\n * Output tabular filename\n \n+Optional fourth argument:\n+ * Max coverage depth (integer)\n+\n This messes about with the filenames to make samtools happy, then\n runs "samtools idxstats" and "samtools depth", captures and combines\n the output to the desired output tabular file.\n+\n+Because "samtools depth" treats the max depth a little fuzzily, this\n+tool tries to account for this and applies a clear max-depth cut off.\n """\n-import sys\n+\n import os\n import subprocess\n+import sys\n import tempfile\n \n if "-v" in sys.argv or "--version" in sys.argv:\n- #Galaxy seems to invert the order of the two lines\n- print("BAM coverage statistics v0.0.1")\n+ # Galaxy seems to invert the order of the two lines\n+ print("BAM coverage statistics v0.0.5")\n cmd = "samtools 2>&1 | grep -i ^Version"\n sys.exit(os.system(cmd))\n \n-def stop_err(msg, error_level=1):\n- """Print error message to stdout and quit with given error level."""\n- sys.stderr.write("%s\\n" % msg)\n- sys.exit(error_level)\n-\n-if len(sys.argv) != 4:\n- stop_err("Require three arguments: BAM, BAI, tabular filenames")\n-\n-bam_filename, bai_filename, tabular_filename = sys.argv[1:]\n+# TODO - Proper command line API\n+if len(sys.argv) == 4:\n+ bam_filename, bai_filename, tabular_filename = sys.argv[1:]\n+ max_depth = "8000"\n+elif len(sys.argv) == 5:\n+ bam_filename, bai_filename, tabular_filename, max_depth = sys.argv[1:]\n+else:\n+ sys.exit("Require 3 or 4 arguments: BAM, BAI, tabular filename, [max depth]")\n \n if not os.path.isfile(bam_filename):\n- stop_err("Input BAM file not found: %s" % bam_filename)\n+ sys.exit("Input BAM file not found: %s" % bam_filename)\n+if bai_filename == "-":\n+ # Make this optional for ease of use at the command line by hand:\n+ if os.path.isfile(bam_filename + ".bai"):\n+ bai_filename = bam_filename + ".bai"\n if not os.path.isfile(bai_filename):\n if bai_filename == "None":\n- stop_err("Error: Galaxy did not index your BAM file")\n- stop_err("Input BAI file not found: %s" % bai_filename)\n+ sys.exit("Error: Galaxy did not index your BAM file")\n+ sys.exit("Input BAI file not found: %s" % bai_filename)\n \n-#Assign sensible names with real extensions, and setup symlinks:\n+try:\n+ max_depth = int(max_depth)\n+except ValueError:\n+ sys.exit("Bad argument for max depth: %r" % max_depth)\n+if max_depth < 0:\n+ sys.exit("Bad argument for max depth: %r" % max_depth)\n+\n+# fuzz factor to ensure can reach max_depth, e.g. region with\n+# many reads having a deletion present could underestimate the\n+# coverage by capping the number of reads considered\n+depth_margin = 100\n+\n+# Assign sensible names with real extensions, and setup symlinks:\n tmp_dir = tempfile.mkdtemp()\n bam_file = os.path.join(tmp_dir, "temp.bam")\n bai_file = os.path.join(tmp_dir, "temp.bam.bai")\n@@ -50,27 +72,58 @@\n assert os.path.isfile(bai_file), bai_file\n assert os.path.isfile(bam_file + ".bai"), bam_file\n \n+\n def clean_up():\n+ """Remove our temporary files and directory."""\n os.remove(idxstats_filename)\n os.remove(depth_filename)\n os.remove(bam_file)\n os.remove(bai_file)\n os.rmdir(tmp_dir)\n \n+\n+def samtools_depth_opt_available():\n+ """Determine if samtools depth supports maximum coverage argument."""\n+ child = subprocess.Popen(["samtools", "depth"],\n+ stdout=subprocess.PIPE, stderr=subprocess.STDOUT)\n+ # Combined stdout/stderr in case samtools is ever inconsistent\n+ output, tmp = child.communicate()\n+ assert tmp is None\n+ # Expect to find this line in the help text, exact wording could change:\n+ # -d/-m <int> maximum coverage depth [8000]\n+ return " -d/-m " in output\n+\n+\n+depth_hack = False\n+if not samtools_depth_opt_available():\n+ if max_depth + depth_margin <= 8000:\n+ sys.stderr.write("WARNING: The version of samtools depth installed does not "\n+ "support the'..b'ithout the -d option, but hard coded default\n+ # of 8000 should be fine even allowing a margin for fuzzy output\n+ cmd = \'samtools depth "%s" > "%s"\' % (bam_file, depth_filename)\n+else:\n+ cmd = \'samtools depth -d %i "%s" > "%s"\' % (max_depth + depth_margin, bam_file, depth_filename)\n return_code = os.system(cmd)\n if return_code:\n clean_up()\n- stop_err("Return code %i from command:\\n%s" % (return_code, cmd))\n+ sys.exit("Return code %i from command:\\n%s" % (return_code, cmd))\n+\n \n def load_total_coverage(depth_handle, identifier, length):\n """Parse some of the \'samtools depth\' output for coverages.\n@@ -86,18 +139,23 @@\n # print("%s coverage calculation, length %i, ..." % (identifier, length))\n \n if depth_ref is None:\n- # Right at start of file!\n+ # Right at start of file / new contig\n line = depth_handle.readline()\n+ # Are we at the end of the file?\n+ if not line:\n+ # Must be at the end of the file.\n+ # This can happen if the file contig(s) had no reads mapped\n+ return 0, 0, 0.0, 0\n depth_ref, depth_pos, depth_reads = line.rstrip("\\n").split()\n depth_pos = int(depth_pos)\n- depth_reads = int(depth_reads)\n+ depth_reads = min(max_depth, int(depth_reads))\n # Can now treat as later references where first line cached\n elif identifier != depth_ref:\n # Infer that identifier had coverage zero,\n # and so was not in the \'samtools depth\'\n # output.\n # print("%s appears to have no coverage at all" % identifier)\n- return 0, 0, 0.0\n+ return 0, 0, 0.0, 0\n \n # Good, at start of expected reference\n bases = depth_reads\n@@ -115,7 +173,7 @@\n for line in depth_handle:\n ref, pos, depth = line.rstrip("\\n").split()\n pos = int(pos)\n- depth = int(depth)\n+ depth = min(max_depth, int(depth))\n if ref != identifier:\n # Reached the end of this identifier\'s coverage\n # so cache this ready for next identifier\n@@ -135,7 +193,8 @@\n # print("%s has no coverage at end" % identifier)\n min_cov = 0\n mean_cov = bases / float(length)\n- return min_cov, max_cov, mean_cov\n+ return min_cov, max_cov, mean_cov, bases\n+\n \n # Parse and combine the output\n out_handle = open(tabular_filename, "w")\n@@ -148,6 +207,8 @@\n depth_pos = 0\n depth_reads = 0\n \n+global_bases = 0\n+global_length = 0\n for line in idxstats_handle:\n identifier, length, mapped, placed = line.rstrip("\\n").split()\n length = int(length)\n@@ -157,8 +218,12 @@\n min_cov = 0\n max_cov = 0\n mean_cov = 0.0\n+ bases = 0\n else:\n- min_cov, max_cov, mean_cov = load_total_coverage(depth_handle, identifier, length)\n+ min_cov, max_cov, mean_cov, bases = load_total_coverage(depth_handle, identifier, length)\n+ if max_cov > max_depth:\n+ sys.exit("Using max depth %i yet saw max coverage %i for %s"\n+ % (max_depth, max_cov, identifier))\n out_handle.write("%s\\t%i\\t%i\\t%i\\t%i\\t%i\\t%0.2f\\n"\n % (identifier, length, mapped, placed,\n min_cov, max_cov, mean_cov))\n@@ -168,15 +233,19 @@\n depth_handle.close()\n out_handle.close()\n clean_up()\n- stop_err("Problem with coverage for %s, expect min_cov <= mean_cov <= max_cov"\n+ sys.exit("Problem with coverage for %s, expect min_cov <= mean_cov <= max_cov"\n % identifier)\n+ global_length += length\n+ global_bases += bases\n \n idxstats_handle.close()\n depth_handle.close()\n out_handle.close()\n \n+print("Total reference length %i with overall mean coverage %0.2f" % (global_length, float(global_bases) / global_length))\n+\n # Remove the temp symlinks and files:\n clean_up()\n \n if depth_ref is not None:\n- stop_err("Left over output from \'samtools depth\'? %r" % depth_ref)\n+ sys.exit("Left over output from \'samtools depth\'? %r" % depth_ref)\n' |
| b |
| diff -r d1fdfaae5dbe -r 7254ece0c0ff tools/coverage_stats/coverage_stats.xml --- a/tools/coverage_stats/coverage_stats.xml Fri Nov 21 09:43:58 2014 -0500 +++ b/tools/coverage_stats/coverage_stats.xml Thu May 11 12:16:10 2017 -0400 |
| b |
| @@ -1,29 +1,35 @@ -<tool id="coverage_stats" name="BAM coverage statistics" version="0.0.1"> +<tool id="coverage_stats" name="BAM coverage statistics" version="0.0.5"> <description>using samtools idxstats and depth</description> <requirements> - <requirement type="binary">samtools</requirement> - <requirement type="package" version="0.1.19">samtools</requirement> + <requirement type="package" version="1.4.1">samtools</requirement> </requirements> - <version_command interpreter="python">coverage_stats.py --version</version_command> - <command interpreter="python">coverage_stats.py "$input_bam" "${input_bam.metadata.bam_index}" "$out_tabular"</command> + <version_command> +python $__tool_directory__/coverage_stats.py --version + </version_command> + <command detect_errors="aggressive"> +python $__tool_directory__/coverage_stats.py '$input_bam' '${input_bam.metadata.bam_index}' '$out_tabular' '$max_depth' + </command> <inputs> <param name="input_bam" type="data" format="bam" label="Input BAM file" /> + <param name="max_depth" type="integer" min="0" max="10000000" label="Max depth" value="8000" /> </inputs> <outputs> <data name="out_tabular" format="tabular" label="$input_bam.name (coverage stats)" /> </outputs> - <stdio> - <!-- Assume anything other than zero is an error --> - <exit_code range="1:" /> - <exit_code range=":-1" /> - </stdio> <tests> <test> <param name="input_bam" value="ex1.bam" ftype="bam" /> + <param name="max_depth" value="123" /> <output name="out_tabular" file="ex1.coverage_stats.tabular" ftype="tabular" /> </test> <test> + <param name="input_bam" value="ex1.bam" ftype="bam" /> + <param name="max_depth" value="50" /> + <output name="out_tabular" file="ex1.coverage_stats.md50.tabular" ftype="tabular" /> + </test> + <test> <param name="input_bam" value="coverage_test.bam" ftype="bam" /> + <param name="max_depth" value="123" /> <output name="out_tabular" file="coverage_test.coverage_stats.tabular" ftype="tabular" /> </test> </tests> @@ -47,8 +53,8 @@ 2 Reference sequence length 3 Number of mapped reads 4 Number of placed but unmapped reads (typically unmapped partners of mapped reads) - 5 Minimum coverage - 6 Maximum coverage + 5 Minimum coverage (per base of reference) + 6 Maximum coverage (per base of reference) 7 Mean coverage (given to 2 dp) ====== ================================================================================= @@ -77,9 +83,15 @@ .. class:: warningmark -**Note**. There is an internal hard limit of 8000 for the pileup routine in -samtools, meaning the reported coverage from ``samtools depth`` will show -maximum coverage depths *around* 8000. +**Note**. If using this on a mapping BAM file, beware that the coverage counting is +done per base of the reference. This means if your reference has any extra bases +compared to the reads being mapped, those bases will be skipped by CIGAR D operators +and these "extra" bases can have an extremely low coverage, giving a potentially +misleading ``min_cov`` values. A sliding window coverage may be more appropriate. + +**Note**. Up until samtools 1.2, there was an internal hard limit of 8000 for the +pileup routine, meaning the reported coverage from ``samtools depth`` would show +maximum coverage depths *around* 8000. This is now a run time option. **Citation** |
| b |
| diff -r d1fdfaae5dbe -r 7254ece0c0ff tools/coverage_stats/tool_dependencies.xml --- a/tools/coverage_stats/tool_dependencies.xml Fri Nov 21 09:43:58 2014 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -1,6 +0,0 @@ -<?xml version="1.0"?> -<tool_dependency> - <package name="samtools" version="0.1.19"> - <repository changeset_revision="923adc89c666" name="package_samtools_0_1_19" owner="iuc" toolshed="https://toolshed.g2.bx.psu.edu" /> - </package> -</tool_dependency> |