From a22f9794c73a3a25ec763e8c0f038c01fbc60e49 Mon Sep 17 00:00:00 2001 From: uzbit Date: Wed, 23 Jul 2014 17:10:53 -0700 Subject: [PATCH 01/20] Update bam_tview_curses.c Large scroll up and down with cap K, J --- bam_tview_curses.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/bam_tview_curses.c b/bam_tview_curses.c index ce3f0b7fb..e0bd873ac 100644 --- a/bam_tview_curses.c +++ b/bam_tview_curses.c @@ -225,8 +225,10 @@ static int curses_loop(tview_t* tv) case ' ': pos += tv->mcol; break; case KEY_UP: case 'j': --tv->row_shift; break; + case 'J': tv->row_shift -= 20; break; case KEY_DOWN: case 'k': ++tv->row_shift; break; + case 'K': tv->row_shift += 20; break; case KEY_BACKSPACE: case '\177': pos -= tv->mcol; break; case KEY_RESIZE: getmaxyx(stdscr, tv->mrow, tv->mcol); break; From 54a91c14965d8ba9bd82fdb99351567bcf0b5d5b Mon Sep 17 00:00:00 2001 From: George Hall Date: Tue, 12 Apr 2016 15:14:38 +0100 Subject: [PATCH 02/20] Implemented tab completion I have added a script which determines all Samtools subcommands and the options available for these subcommands. It then generates a file which, when sourced in Bash, allows for tab completion of all subcommands, and, when appropriate, the long options for the relevant subcommand. --- misc/generate_completion_file.sh | 115 +++++++++++++++++++++++++++++++ 1 file changed, 115 insertions(+) create mode 100644 misc/generate_completion_file.sh diff --git a/misc/generate_completion_file.sh b/misc/generate_completion_file.sh new file mode 100644 index 000000000..baa3a6b11 --- /dev/null +++ b/misc/generate_completion_file.sh @@ -0,0 +1,115 @@ +################################################################################ +# Copyright (c) 2016 Genome Research Ltd. +# +# Author: George Hall +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject to +# the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +################################################################################ + + +# This script uses usage messages to compute all Samtools subcommands and their +# corresponding options, and then prints a file allowing for tab completion of +# these subcommands and their options to stdout. The user must then source this +# file from their .bashrc. + +# Usage: +# Step 1: Check SAMTOOLS_BIN is correctly set +# Step 2: $ sh generate_completion_file.sh > /where/to/save/completion/file +# Step 3: Source completion file (probably want to add it to be sourced from .bashrc) + + +# The completion file generated will only work for this binary, change if necessary: +SAMTOOLS_BIN="samtools" + +# Get subcommands: + +SUBCOMMANDS=$($SAMTOOLS_BIN 2>&1 | \ + awk '{if ($1 != "" && \ + $1 != "Program:" && \ + $1 != "Usage:" && \ + $1 != "Version:" && \ + $1 != "Commands:" && \ + $1 != "--") print $1}' | \ + xargs -n 1 -I @ printf @" ") + +SAMTOOLS_VERSION=$(samtools 2>&1 | grep "Version: [[:graph:]]*" | awk '{print $2}') +printf "# Generated by $0 using Samtools version $SAMTOOLS_VERSION\n\n" + +printf "# All subcommands can be tab completed, and those with long options can be tab completed\n" +printf "# once the first '-' has been provided. Otherwise, it is assumed that the user is trying\n" +printf "# to tab complete filenames.\n\n" + +printf "# This file must be sourced for the tab completions to work\n\n" + +printf "_samtools_options()\n" +printf "{\n" +printf "\tlocal cur current_sub opts\n" +printf "\tCOMPREPLY=()\n" +printf "\tcur=\"\${COMP_WORDS[COMP_CWORD]}\"\n\n" + +printf "\tif [[ \$COMP_CWORD == 1 ]] ; then\n" +printf "\t\t# Complete main subcommand\n" +printf "\t\topts=\"$SUBCOMMANDS\"\n\n" + +printf "\t\tCOMPREPLY=(\$(compgen -W \"\${opts}\" -- \${cur}))\n" +printf "\t\treturn 0\n\n" + +printf "\telse\n" +printf "\t\t# Complete options for subcommands\n\n" +printf "\t\tcurrent_sub=\"\${COMP_WORDS[1]}\"\n\n" + +printf "\t\tcase \$current_sub in\n\n" + +for SUB in $SUBCOMMANDS; do + + # This matches with two double dashed (i.e. a long option) followed by alphanums + # or dashes. The reason I couldn't use [[:graph:]] is that it could include commas, + # which are used to separate the short version of long options, which start with + # double dashes when they are >= 2 chars long. + SUB_OPTS=$($SAMTOOLS_BIN $SUB 2>&1 | \ + grep -oh "\(\-\-\)\([[:alnum:]]\|\-\)* " | \ + xargs -I @ printf @) + + # Only print options when option list is non-empty + OPTS_LEN=${#SUB_OPTS} + if [ $OPTS_LEN != 0 ]; then + printf "\t\t\t\"$SUB\")\n" + printf "\t\t\t\t# Completion for $SUB - only triggered after first '-'\n" + printf "\t\t\t\tif [[ \${cur} == -* ]] ; then\n" + printf "\t\t\t\t\topts=\"$SUB_OPTS\"\n\n" + printf "\t\t\t\t\tCOMPREPLY=(\$(compgen -W \"\${opts}\" -- \${cur}))\n" + printf "\t\t\t\t\treturn 0\n" + printf "\t\t\t\tfi\n" + printf "\t\t\t\t;;\n\n" + fi +done + +printf "\t\t\t*)\n" +printf "\t\t\t\t;;\n\n" +printf "\t\tesac\n\n" + +printf "\tfi\n\n" + +printf "\t# If we have not returned by now, assume the user is looking for files\n" +printf "\tcompopt -o filenames # Don't append a space\n" +printf "\tCOMPREPLY=(\$(compgen -f \"\${cur}\"))\n" +printf "\treturn 0\n" +printf "}\n" +printf "complete -F _samtools_options $SAMTOOLS_BIN\n" From fd02fc5614bd7ef6214c4bddcf4d941000670f95 Mon Sep 17 00:00:00 2001 From: George Hall Date: Wed, 13 Apr 2016 12:36:56 +0100 Subject: [PATCH 03/20] Dynamically generate options As was apparently mentioned at the meeting, it is probably better to dynamically generate tab completions rather than to generate a one-off script which the user then stores. This is what this new version does - possible completions for subcommands and long-options are generated on the fly, meaning that it is no longer necessary to try and keep the completion file up-to-date. Additionally, this allows for support of mulitple versions etc, by append lines to the bottom of the file, as noted in the comments. I have removed the previous version. --- misc/generate_completion_file.sh | 115 ------------------------------- misc/samtools_tab_completion | 81 ++++++++++++++++++++++ 2 files changed, 81 insertions(+), 115 deletions(-) delete mode 100644 misc/generate_completion_file.sh create mode 100644 misc/samtools_tab_completion diff --git a/misc/generate_completion_file.sh b/misc/generate_completion_file.sh deleted file mode 100644 index baa3a6b11..000000000 --- a/misc/generate_completion_file.sh +++ /dev/null @@ -1,115 +0,0 @@ -################################################################################ -# Copyright (c) 2016 Genome Research Ltd. -# -# Author: George Hall -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject to -# the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF -# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. -# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY -# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, -# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE -# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -################################################################################ - - -# This script uses usage messages to compute all Samtools subcommands and their -# corresponding options, and then prints a file allowing for tab completion of -# these subcommands and their options to stdout. The user must then source this -# file from their .bashrc. - -# Usage: -# Step 1: Check SAMTOOLS_BIN is correctly set -# Step 2: $ sh generate_completion_file.sh > /where/to/save/completion/file -# Step 3: Source completion file (probably want to add it to be sourced from .bashrc) - - -# The completion file generated will only work for this binary, change if necessary: -SAMTOOLS_BIN="samtools" - -# Get subcommands: - -SUBCOMMANDS=$($SAMTOOLS_BIN 2>&1 | \ - awk '{if ($1 != "" && \ - $1 != "Program:" && \ - $1 != "Usage:" && \ - $1 != "Version:" && \ - $1 != "Commands:" && \ - $1 != "--") print $1}' | \ - xargs -n 1 -I @ printf @" ") - -SAMTOOLS_VERSION=$(samtools 2>&1 | grep "Version: [[:graph:]]*" | awk '{print $2}') -printf "# Generated by $0 using Samtools version $SAMTOOLS_VERSION\n\n" - -printf "# All subcommands can be tab completed, and those with long options can be tab completed\n" -printf "# once the first '-' has been provided. Otherwise, it is assumed that the user is trying\n" -printf "# to tab complete filenames.\n\n" - -printf "# This file must be sourced for the tab completions to work\n\n" - -printf "_samtools_options()\n" -printf "{\n" -printf "\tlocal cur current_sub opts\n" -printf "\tCOMPREPLY=()\n" -printf "\tcur=\"\${COMP_WORDS[COMP_CWORD]}\"\n\n" - -printf "\tif [[ \$COMP_CWORD == 1 ]] ; then\n" -printf "\t\t# Complete main subcommand\n" -printf "\t\topts=\"$SUBCOMMANDS\"\n\n" - -printf "\t\tCOMPREPLY=(\$(compgen -W \"\${opts}\" -- \${cur}))\n" -printf "\t\treturn 0\n\n" - -printf "\telse\n" -printf "\t\t# Complete options for subcommands\n\n" -printf "\t\tcurrent_sub=\"\${COMP_WORDS[1]}\"\n\n" - -printf "\t\tcase \$current_sub in\n\n" - -for SUB in $SUBCOMMANDS; do - - # This matches with two double dashed (i.e. a long option) followed by alphanums - # or dashes. The reason I couldn't use [[:graph:]] is that it could include commas, - # which are used to separate the short version of long options, which start with - # double dashes when they are >= 2 chars long. - SUB_OPTS=$($SAMTOOLS_BIN $SUB 2>&1 | \ - grep -oh "\(\-\-\)\([[:alnum:]]\|\-\)* " | \ - xargs -I @ printf @) - - # Only print options when option list is non-empty - OPTS_LEN=${#SUB_OPTS} - if [ $OPTS_LEN != 0 ]; then - printf "\t\t\t\"$SUB\")\n" - printf "\t\t\t\t# Completion for $SUB - only triggered after first '-'\n" - printf "\t\t\t\tif [[ \${cur} == -* ]] ; then\n" - printf "\t\t\t\t\topts=\"$SUB_OPTS\"\n\n" - printf "\t\t\t\t\tCOMPREPLY=(\$(compgen -W \"\${opts}\" -- \${cur}))\n" - printf "\t\t\t\t\treturn 0\n" - printf "\t\t\t\tfi\n" - printf "\t\t\t\t;;\n\n" - fi -done - -printf "\t\t\t*)\n" -printf "\t\t\t\t;;\n\n" -printf "\t\tesac\n\n" - -printf "\tfi\n\n" - -printf "\t# If we have not returned by now, assume the user is looking for files\n" -printf "\tcompopt -o filenames # Don't append a space\n" -printf "\tCOMPREPLY=(\$(compgen -f \"\${cur}\"))\n" -printf "\treturn 0\n" -printf "}\n" -printf "complete -F _samtools_options $SAMTOOLS_BIN\n" diff --git a/misc/samtools_tab_completion b/misc/samtools_tab_completion new file mode 100644 index 000000000..9a778559d --- /dev/null +++ b/misc/samtools_tab_completion @@ -0,0 +1,81 @@ +################################################################################ +# Copyright (c) 2016 Genome Research Ltd. +# +# Author: George Hall +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject to +# the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +################################################################################ + +# This script determines whether the user is trying to complete a subcommand or +# a long-option, then determines and displays possible completions. For commands +# with long-options, the initial '-' must be first be typed in order to trigger +# the completion mechanism. + +# This file must be sourced for the tab completions to work - it is advisable to +# source it from your .bashrc. By default, it will only perform tab completion for +# a binary called 'samtools'. To enable tab completion for other versions of +# Samtools, or for Samtools binaries with different names, simply append +# 'complete -F _samtools_options ' to the end of this file. This +# file will then need to be sourced again. + +_samtools_options() +{ + + SAMTOOLS_BIN="${COMP_WORDS[0]}" + local cur current_sub opts + COMPREPLY=() + cur="${COMP_WORDS[COMP_CWORD]}" + + if [[ $COMP_CWORD == 1 ]] ; then + + # If on the first word, generate possible subcommands, and tab complete those + + opts=$($SAMTOOLS_BIN 2>&1 | awk '{if ($1 != "" && $1 != "Program:" && \ + $1 != "Usage:" && $1 != "Version:" && $1 != "Commands:" && \ + $1 != "--") print $1}' | xargs -n 1 -I @ printf @" ") + + COMPREPLY=($(compgen -W "${opts}" -- ${cur})) + return 0 + + else + + # Complete long-options (if available) for current subcommand + + current_sub="${COMP_WORDS[1]}" + + opts=$($SAMTOOLS_BIN $current_sub 2>&1 | grep -oh "\(\-\-\)\([[:alnum:]]\|\-\)* " | \ + xargs -I @ printf @) + + if [[ ${cur} == -* ]] ; then + COMPREPLY=($(compgen -W "${opts}" -- ${cur})) + return 0 + fi + + fi + + # If we have not returned by now, assume that the user is looking for files + compopt -o filenames # Don't append a space + COMPREPLY=($(compgen -f "${cur}")) + return 0 + +} + +complete -F _samtools_options samtools + From 248534509594ea30f95eec30cf7f30129d332a6a Mon Sep 17 00:00:00 2001 From: George Hall Date: Wed, 13 Apr 2016 14:21:25 +0100 Subject: [PATCH 04/20] Cleaned up subcommand name retrieval --- misc/samtools_tab_completion | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/misc/samtools_tab_completion b/misc/samtools_tab_completion index 9a778559d..e059aaa67 100644 --- a/misc/samtools_tab_completion +++ b/misc/samtools_tab_completion @@ -47,12 +47,10 @@ _samtools_options() # If on the first word, generate possible subcommands, and tab complete those - opts=$($SAMTOOLS_BIN 2>&1 | awk '{if ($1 != "" && $1 != "Program:" && \ - $1 != "Usage:" && $1 != "Version:" && $1 != "Commands:" && \ - $1 != "--") print $1}' | xargs -n 1 -I @ printf @" ") + opts=$($SAMTOOLS_BIN 2>&1 | grep '^ ' | awk '{print $1}') - COMPREPLY=($(compgen -W "${opts}" -- ${cur})) - return 0 + COMPREPLY=($(compgen -W "${opts}" -- ${cur})) + return 0 else From 193eecb56d255aaa703644a9d4d39fbdd73b3f7e Mon Sep 17 00:00:00 2001 From: George Hall Date: Wed, 13 Apr 2016 14:28:58 +0100 Subject: [PATCH 05/20] Changed '==' to '-eq' --- misc/samtools_tab_completion | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/misc/samtools_tab_completion b/misc/samtools_tab_completion index e059aaa67..a2bb44683 100644 --- a/misc/samtools_tab_completion +++ b/misc/samtools_tab_completion @@ -43,7 +43,7 @@ _samtools_options() COMPREPLY=() cur="${COMP_WORDS[COMP_CWORD]}" - if [[ $COMP_CWORD == 1 ]] ; then + if [[ $COMP_CWORD -eq 1 ]] ; then # If on the first word, generate possible subcommands, and tab complete those From 0bfa8d4249fc4035043e1aa3bb63bca0cdf6f7e3 Mon Sep 17 00:00:00 2001 From: George Hall Date: Thu, 14 Apr 2016 16:55:32 +0100 Subject: [PATCH 06/20] Captialised variables for consistency --- misc/samtools_tab_completion | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/misc/samtools_tab_completion b/misc/samtools_tab_completion index a2bb44683..7c60b8fe4 100644 --- a/misc/samtools_tab_completion +++ b/misc/samtools_tab_completion @@ -39,30 +39,30 @@ _samtools_options() { SAMTOOLS_BIN="${COMP_WORDS[0]}" - local cur current_sub opts + local CUR CURRENT_SUB OPTS COMPREPLY=() - cur="${COMP_WORDS[COMP_CWORD]}" + CUR="${COMP_WORDS[COMP_CWORD]}" if [[ $COMP_CWORD -eq 1 ]] ; then # If on the first word, generate possible subcommands, and tab complete those - opts=$($SAMTOOLS_BIN 2>&1 | grep '^ ' | awk '{print $1}') + OPTS=$($SAMTOOLS_BIN 2>&1 | grep '^ ' | awk '{print $1}') - COMPREPLY=($(compgen -W "${opts}" -- ${cur})) + COMPREPLY=($(compgen -W "${OPTS}" -- ${CUR})) return 0 else # Complete long-options (if available) for current subcommand - current_sub="${COMP_WORDS[1]}" + CURRENT_SUB="${COMP_WORDS[1]}" - opts=$($SAMTOOLS_BIN $current_sub 2>&1 | grep -oh "\(\-\-\)\([[:alnum:]]\|\-\)* " | \ + OPTS=$($SAMTOOLS_BIN $CURRENT_SUB 2>&1 | grep -oh "\(\-\-\)\([[:alnum:]]\|\-\)* " | \ xargs -I @ printf @) - if [[ ${cur} == -* ]] ; then - COMPREPLY=($(compgen -W "${opts}" -- ${cur})) + if [[ ${CUR} == -* ]] ; then + COMPREPLY=($(compgen -W "${OPTS}" -- ${CUR})) return 0 fi @@ -70,7 +70,7 @@ _samtools_options() # If we have not returned by now, assume that the user is looking for files compopt -o filenames # Don't append a space - COMPREPLY=($(compgen -f "${cur}")) + COMPREPLY=($(compgen -f "${CUR}")) return 0 } From 2389a31d83b6baf8dd8cfc6005a32ab5e5798d6a Mon Sep 17 00:00:00 2001 From: George Hall Date: Fri, 15 Apr 2016 12:51:09 +0100 Subject: [PATCH 07/20] Fixed slight bug in default behaviour Now handles files with spaces in their names correctly. In doing so, I have basically made the change mentioned here: https://github.com/samtools/samtools/pull/560#discussion-diff-59540631 I have made the change in a slightly different location, however, as this has the advantage of not allowing the user to expand '/' or similar when they should be entering a subcommand (i.e. when completing the first argument). --- misc/samtools_tab_completion | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/misc/samtools_tab_completion b/misc/samtools_tab_completion index 7c60b8fe4..43de19904 100644 --- a/misc/samtools_tab_completion +++ b/misc/samtools_tab_completion @@ -64,14 +64,12 @@ _samtools_options() if [[ ${CUR} == -* ]] ; then COMPREPLY=($(compgen -W "${OPTS}" -- ${CUR})) return 0 + else + # Assume the user wants normal file name completion + compopt -o default fi - - fi - # If we have not returned by now, assume that the user is looking for files - compopt -o filenames # Don't append a space - COMPREPLY=($(compgen -f "${CUR}")) - return 0 + fi } From f5120befe07a9bcee844b66fbdc27b59a504b234 Mon Sep 17 00:00:00 2001 From: awenger Date: Mon, 8 Dec 2014 23:09:29 -0800 Subject: [PATCH 08/20] Add "-b " flag to samtools cat. Add a "-b " like that in "samtools merge" to "samtools cat". This alleviates any problems with the list of input file name exceeding the number of characters permitted in a single command. [jkb during merge: adjusted memory freeing and usage format.] --- bam_cat.c | 57 ++++++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 50 insertions(+), 7 deletions(-) diff --git a/bam_cat.c b/bam_cat.c index 5c303d1ed..0d8940b0d 100644 --- a/bam_cat.c +++ b/bam_cat.c @@ -531,10 +531,12 @@ int main_cat(int argc, char *argv[]) { bam_hdr_t *h = 0; char *outfn = 0; + char **infns = NULL; // files to concatenate + int infns_size = 0; int c, ret = 0; samFile *in; - while ((c = getopt(argc, argv, "h:o:")) >= 0) { + while ((c = getopt(argc, argv, "h:o:b:")) >= 0) { switch (c) { case 'h': { samFile *fph = sam_open(optarg, "r"); @@ -553,29 +555,61 @@ int main_cat(int argc, char *argv[]) break; } case 'o': outfn = strdup(optarg); break; + case 'b': { + // add file names in "optarg" to the list + // of files to concatenate + int nfns; + char **fns_read = hts_readlines(optarg, &nfns); + if (fns_read) { + infns = realloc(infns, (infns_size + nfns) * sizeof(char*)); + if (infns == NULL) { ret = 1; goto end; } + memcpy(infns+infns_size, fns_read, nfns * sizeof(char*)); + infns_size += nfns; + free(fns_read); + } else { + print_error("cat", "Invalid file list \"%s\"", optarg); + ret = 1; + } + break; + } } } - if (argc - optind < 1) { - fprintf(stderr, "Usage: samtools cat [-h header.sam] [-o out.bam] [...]\n"); + + // Append files specified in argv to the list. + int nargv_fns = argc - optind; + if (nargv_fns > 0) { + infns = realloc(infns, (infns_size + nargv_fns) * sizeof(char*)); + if (infns == NULL) { ret = 1; goto end; } + memcpy(infns + infns_size, argv + optind, nargv_fns * sizeof(char*)); + } + + // Require at least one input file + if (infns_size + nargv_fns == 0) { + fprintf(stderr, "Usage: samtools cat [options] [... ]\n"); + fprintf(stderr, " samtools cat [options] [... ]\n\n"); + fprintf(stderr, "Concatenate BAM or CRAM files, first those in , then those\non the command line.\n\n"); + fprintf(stderr, "Options: -b FILE list of input BAM/CRAM file names, one per line\n"); + fprintf(stderr, " -h FILE copy the header from FILE [default is 1st input file]\n"); + fprintf(stderr, " -o FILE output BAM/CRAM\n"); return 1; } - in = sam_open(argv[optind], "r"); + in = sam_open(infns[0], "r"); if (!in) { - print_error_errno("cat", "failed to open file '%s'", argv[optind]); + print_error_errno("cat", "failed to open file '%s'", infns[0]); return 1; } switch (hts_get_format(in)->format) { case bam: sam_close(in); - if (bam_cat(argc - optind, argv + optind, h, outfn? outfn : "-") < 0) + if (bam_cat(infns_size+nargv_fns, infns, h, outfn? outfn : "-") < 0) ret = 1; break; case cram: sam_close(in); - if (cram_cat(argc - optind, argv + optind, h, outfn? outfn : "-") < 0) + if (cram_cat(infns_size+nargv_fns, infns, h, outfn? outfn : "-") < 0) ret = 1; break; @@ -584,7 +618,16 @@ int main_cat(int argc, char *argv[]) fprintf(stderr, "[%s] ERROR: input is not BAM or CRAM\n", __func__); return 1; } + + end: + if (infns_size > 0) { + int i; + for (i=0; i Date: Wed, 22 Mar 2017 14:38:42 +0000 Subject: [PATCH 09/20] Document the "-b bam-list" option to samtools cat. --- samtools.1 | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/samtools.1 b/samtools.1 index d87a88e4e..a2bce4248 100644 --- a/samtools.1 +++ b/samtools.1 @@ -1482,14 +1482,33 @@ The amount of space available will differ for each CRAM file. .TP \"-------- cat .B cat -samtools cat [-h header.sam] [-o out.bam] [ ... ] +samtools cat [-b list] [-h header.sam] [-o out.bam] [ ... ] -Concatenate BAMs. The sequence dictionary of each input BAM must be identical, -although this command does not check this. This command uses a similar trick -to +Concatenate BAMs or CRAMs. Although this works on either BAM or CRAM, +all input files must be the same format as each other. The sequence +dictionary of each input file must be identical, although this command +does not check this. This command uses a similar trick to .B reheader which enables fast BAM concatenation. +.B OPTIONS: +.RS +.TP 8 +.BI "-b " FOFN +Read the list of input BAM or CRAM files from \fIFOFN\fR. These are +concatenated prior to any files specified on the command line. +Multiple \fB-b\fR \fIFOFN\fR options may be specified to concatenate +multiple lists of BAM/CRAM files. +.TP 8 +.BI "-h " FILE +Uses the SAM header from \fIFILE\fR. By default the header is taken +from the first file to be concatenated. +.TP 8 +.BI "-o " FILE +Write the concatenated output to \fIFILE\fR. By default this is sent +to stdout. +.RE + .TP \"-------- rmdup .B rmdup samtools rmdup [-sS] From ea171c4db9969e7a46f4fb58a986541112b6b1a2 Mon Sep 17 00:00:00 2001 From: James Bonfield Date: Wed, 22 Mar 2017 14:53:43 +0000 Subject: [PATCH 10/20] Remove two tiny memory leaks in samtools merge. --- bam_sort.c | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/bam_sort.c b/bam_sort.c index 0ff2b9305..5dab20e5e 100644 --- a/bam_sort.c +++ b/bam_sort.c @@ -1147,14 +1147,6 @@ int bam_merge_core2(int by_qname, const char *out, const char *mode, print_error(cmd, "couldn't read headers from \"%s\"", headers); goto mem_fail; } - } else { - hout = bam_hdr_init(); - if (!hout) { - print_error(cmd, "couldn't allocate bam header"); - goto mem_fail; - } - hout->text = strdup(""); - if (!hout->text) goto mem_fail; } g_is_by_qname = by_qname; @@ -1504,6 +1496,7 @@ int bam_merge(int argc, char *argv[]) if (fn == NULL) { ret = 1; goto end; } memcpy(fn+fn_size, fn_read, nfiles * sizeof(char*)); fn_size += nfiles; + free(fn_read); } else { print_error("merge", "Invalid file list \"%s\"", optarg); From 692e058392cf9337deac8cbcda74761e73169972 Mon Sep 17 00:00:00 2001 From: James Bonfield Date: Thu, 23 Mar 2017 11:43:35 +0000 Subject: [PATCH 11/20] Cram lossy names (#543) * Added missing input/output options to the man page. * Rewording of input/output options. Given the command line options are --input-fmt-options and --output-fmt-options, it seemed wiser to use the phrases "CRAM input only" and "CRAM output only" when documenting the input/output option strings. It still refers to encoding and decoding elsewhere, but I think this is sufficiently clear. --- samtools.1 | 36 +++++++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 9 deletions(-) diff --git a/samtools.1 b/samtools.1 index a2bce4248..57182deae 100644 --- a/samtools.1 +++ b/samtools.1 @@ -1751,17 +1751,17 @@ It usually is not required for decoding except in the situation of the MD5 not being obtainable via the REF_PATH or REF_CACHE environment variables. .TP .BI decode_md= 0|1 -CRAM decode only; defaults to 1 (on). CRAM does not typically store +CRAM input only; defaults to 1 (on). CRAM does not typically store MD and NM tags, preferring to generate them on the fly. This option controls this behaviour. .TP .BI ignore_md5= 0|1 -CRAM decode only; defaults to 0 (off). When enabled, md5 checksum +CRAM input only; defaults to 0 (off). When enabled, md5 checksum errors on the reference sequence and block checksum errors within CRAM are ignored. Use of this option is strongly discouraged. .TP .BI required_fields= bit-field -CRAM decode only; specifies which SAM columns need to be populated. +CRAM input only; specifies which SAM columns need to be populated. By default all fields are used. Limiting the decode to specific columns can have significant performance gains. The bit-field is a numerical value constructed from the following table. @@ -1783,34 +1783,52 @@ rb l . 0x1000 SAM_RGAUX .TE .TP +.BI name_prefix= string +CRAM input only; defaults to output filename. Any sequences with +auto-generated read names will use \fIstring\fR as the name prefix. +.TP .BI multi_seq_per_slice= 0|1 -CRAM encode only; defaults to 0 (off). By default CRAM generates one +CRAM output only; defaults to 0 (off). By default CRAM generates one container per reference sequence, except in the case of many small references (such as a fragmented assembly). .TP .BI version= major.minor -CRAM encode only. Specifies the CRAM version number. Acceptable +CRAM output only. Specifies the CRAM version number. Acceptable values are "2.1" and "3.0". .TP .BI seqs_per_slice= INT -CRAM encode only; defaults to 10000. +CRAM output only; defaults to 10000. .TP .BI slices_per_container= INT -CRAM encode only; defaults to 1. The effect of having multiple slices +CRAM output only; defaults to 1. The effect of having multiple slices per container is to share the compression header block between multiple slices. This is unlikely to have any significant impact unless the number of sequences per slice is reduced. (Together these two options control the granularity of random access.) .TP .BI embed_ref= 0|1 -CRAM encode only; defaults to 0 (off). If 1, this will store portions +CRAM output only; defaults to 0 (off). If 1, this will store portions of the reference sequence in each slice, permitting decode without having requiring an external copy of the reference sequence. .TP .BI no_ref= 0|1 -CRAM encode only; defaults to 0 (off). If 1, sequences will be stored +CRAM output only; defaults to 0 (off). If 1, sequences will be stored verbatim with no reference encoding. This can be useful if no reference is available for the file. +.TP +.BI use_bzip2= 0|1 +CRAM output only; defaults to 0 (off). Permits use of bzip2 in CRAM +block compression. +.TP +.BI use_lzma= 0|1 +CRAM output only; defaults to 0 (off). Permits use of lzma in CRAM +block compression. +.TP +.BI lossy_names= 0|1 +CRAM output only; defaults to 0 (off). If 1, templates with all +members within the same CRAM slice will have their read names +removed. New names will be automatically generated during decoding. +Also see the \fBname_prefix\fR option. .RE .PP For example: From 8c20354dc78540b3a29217945321752a9db4fc52 Mon Sep 17 00:00:00 2001 From: Blaise Li Date: Fri, 18 Dec 2015 17:26:10 +0100 Subject: [PATCH 12/20] Option -G to exclude based on a set of bits. I wanted to be able to exclude reads having a given combination of bits set in the SAM flag. For instance, to exclude reads where both mates are unmapped, use `samtools view -G 12` Did I miss how to do it with -f and -F ? (Committer: jkb - fixed bam2fq bug to handle -G0 correctly.) --- sam_view.c | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/sam_view.c b/sam_view.c index b92f9b53e..3e3f099e7 100644 --- a/sam_view.c +++ b/sam_view.c @@ -51,6 +51,7 @@ typedef struct samview_settings { int min_mapQ; int flag_on; int flag_off; + int flag_alloff; int min_qlen; int remove_B; uint32_t subsam_seed; @@ -84,6 +85,8 @@ static int process_aln(const bam_hdr_t *h, bam1_t *b, samview_settings_t* settin } if (b->core.qual < settings->min_mapQ || ((b->core.flag & settings->flag_on) != settings->flag_on) || (b->core.flag & settings->flag_off)) return 1; + if (settings->flag_alloff && ((b->core.flag & settings->flag_alloff) == settings->flag_alloff)) + return 1; if (settings->bed && (b->core.tid < 0 || !bed_overlap(settings->bed, h->target_name[b->core.tid], b->core.pos, bam_endpos(b)))) return 1; if (settings->subsam_frac > 0.) { @@ -247,6 +250,7 @@ int main_samview(int argc, char *argv[]) .min_mapQ = 0, .flag_on = 0, .flag_off = 0, + .flag_alloff = 0, .min_qlen = 0, .remove_B = 0, .subsam_seed = 0, @@ -264,7 +268,7 @@ int main_samview(int argc, char *argv[]) strcpy(out_mode, "w"); strcpy(out_un_mode, "w"); while ((c = getopt_long(argc, argv, - "SbBcCt:h1Ho:O:q:f:F:ul:r:?T:R:L:s:@:m:x:U:", + "SbBcCt:h1Ho:O:q:f:F:G:ul:r:?T:R:L:s:@:m:x:U:", lopts, NULL)) >= 0) { switch (c) { case 's': @@ -288,6 +292,7 @@ int main_samview(int argc, char *argv[]) case 'U': fn_un_out = strdup(optarg); break; case 'f': settings.flag_on |= strtol(optarg, 0, 0); break; case 'F': settings.flag_off |= strtol(optarg, 0, 0); break; + case 'G': settings.flag_alloff |= strtol(optarg, 0, 0); break; case 'q': settings.min_mapQ = atoi(optarg); break; case 'u': compress_level = 0; break; case '1': compress_level = 1; break; @@ -571,6 +576,7 @@ static int usage(FILE *fp, int exit_status, int is_long_help) " query sequence >= INT [0]\n" " -f INT only include reads with all bits set in INT set in FLAG [0]\n" " -F INT only include reads with none of the bits set in INT set in FLAG [0]\n" +" -G INT only include reads with not all the bits set in INT set in FLAG [0]\n" " -s FLOAT subsample reads (given INT.FRAC option value, 0.FRAC is the\n" " fraction of templates/read pairs to keep; INT part sets seed)\n" // read processing @@ -654,6 +660,7 @@ static void bam2fq_usage(FILE *to, const char *command) " -2 FILE write paired reads flagged READ2 to FILE\n" " -f INT only include reads with all bits set in INT set in FLAG [0]\n" " -F INT only include reads with none of the bits set in INT set in FLAG [0]\n" +" -G INT only include reads with not all the bits set in INT set in FLAG [0]\n" " -n don't append /1 and /2 to the read name\n"); if (fq) fprintf(to, " -O output quality in the OQ tag if present\n"); @@ -673,7 +680,7 @@ typedef struct bam2fq_opts { char *fnr[3]; char *fn_input; // pointer to input filename in argv do not free bool has12, use_oq, copy_tags; - int flag_on, flag_off; + int flag_on, flag_off, flag_alloff; sam_global_args ga; fastfile filetype; int def_qual; @@ -685,7 +692,7 @@ typedef struct bam2fq_state { FILE *fpr[3]; bam_hdr_t *h; bool has12, use_oq, copy_tags; - int flag_on, flag_off; + int flag_on, flag_off, flag_alloff; fastfile filetype; int def_qual; } bam2fq_state_t; @@ -803,13 +810,14 @@ static bool parse_opts(int argc, char *argv[], bam2fq_opts_t** opts_out) SAM_OPT_GLOBAL_OPTIONS('-', 0, '-', '-', 0, '@'), { NULL, 0, NULL, 0 } }; - while ((c = getopt_long(argc, argv, "0:1:2:f:F:nOs:tv:@:", lopts, NULL)) > 0) { + while ((c = getopt_long(argc, argv, "0:1:2:f:F:G:nOs:tv:@:", lopts, NULL)) > 0) { switch (c) { case '0': opts->fnr[0] = optarg; break; case '1': opts->fnr[1] = optarg; break; case '2': opts->fnr[2] = optarg; break; case 'f': opts->flag_on |= strtol(optarg, 0, 0); break; case 'F': opts->flag_off |= strtol(optarg, 0, 0); break; + case 'G': opts->flag_alloff |= strtol(optarg, 0, 0); break; case 'n': opts->has12 = false; break; case 'O': opts->use_oq = true; break; case 's': opts->fnse = optarg; break; @@ -867,6 +875,7 @@ static bool init_state(const bam2fq_opts_t* opts, bam2fq_state_t** state_out) bam2fq_state_t* state = calloc(1, sizeof(bam2fq_state_t)); state->flag_on = opts->flag_on; state->flag_off = opts->flag_off; + state->flag_alloff = opts->flag_alloff; state->has12 = opts->has12; state->use_oq = opts->use_oq; state->copy_tags = opts->copy_tags; @@ -945,7 +954,8 @@ static inline bool filter_it_out(const bam1_t *b, const bam2fq_state_t *state) { return (b->core.flag&(BAM_FSECONDARY|BAM_FSUPPLEMENTARY) // skip secondary and supplementary alignments || (b->core.flag&(state->flag_on)) != state->flag_on // or reads indicated by filter flags - || (b->core.flag&(state->flag_off)) != 0); + || (b->core.flag&(state->flag_off)) != 0 + || (b->core.flag&(state->flag_alloff) && (b->core.flag&(state->flag_alloff)) == state->flag_alloff)); } From b1b6717e7017e5808560b9f808726049f7485f41 Mon Sep 17 00:00:00 2001 From: James Bonfield Date: Thu, 30 Mar 2017 17:33:24 +0100 Subject: [PATCH 13/20] Improved wording of usage and updated man page. --- sam_view.c | 12 ++++++------ samtools.1 | 12 +++++++++++- 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/sam_view.c b/sam_view.c index 3e3f099e7..2cf7ae894 100644 --- a/sam_view.c +++ b/sam_view.c @@ -574,9 +574,9 @@ static int usage(FILE *fp, int exit_status, int is_long_help) " -l STR only include reads in library STR [null]\n" " -m INT only include reads with number of CIGAR operations consuming\n" " query sequence >= INT [0]\n" -" -f INT only include reads with all bits set in INT set in FLAG [0]\n" -" -F INT only include reads with none of the bits set in INT set in FLAG [0]\n" -" -G INT only include reads with not all the bits set in INT set in FLAG [0]\n" +" -f INT only include reads with all of the bits set in INT set in FLAG [0]\n" +" -F INT only include reads with all of the bits unset in INT set in FLAG [0]\n" +" -G INT only include reads with any of the bits unset in INT set in FLAG [0]\n" " -s FLOAT subsample reads (given INT.FRAC option value, 0.FRAC is the\n" " fraction of templates/read pairs to keep; INT part sets seed)\n" // read processing @@ -658,9 +658,9 @@ static void bam2fq_usage(FILE *to, const char *command) " -0 FILE write paired reads flagged both or neither READ1 and READ2 to FILE\n" " -1 FILE write paired reads flagged READ1 to FILE\n" " -2 FILE write paired reads flagged READ2 to FILE\n" -" -f INT only include reads with all bits set in INT set in FLAG [0]\n" -" -F INT only include reads with none of the bits set in INT set in FLAG [0]\n" -" -G INT only include reads with not all the bits set in INT set in FLAG [0]\n" +" -f INT only include reads with all of the bits set in INT set in FLAG [0]\n" +" -F INT only include reads with all of the bits unset in INT set in FLAG [0]\n" +" -G INT only include reads with any of the bits unset in INT set in FLAG [0]\n" " -n don't append /1 and /2 to the read name\n"); if (fq) fprintf(to, " -O output quality in the OQ tag if present\n"); diff --git a/samtools.1 b/samtools.1 index 57182deae..8b8e49dfb 100644 --- a/samtools.1 +++ b/samtools.1 @@ -159,8 +159,9 @@ The .BR -l , .BR -m , .BR -f , +.BR -F , and -.B -F +.B -G options filter the alignments that will be included in the output to only those alignments that match certain criteria. @@ -325,6 +326,15 @@ present in the FLAG field. can be specified in hex by beginning with `0x' (i.e. /^0x[0-9A-F]+/) or in octal by beginning with `0' (i.e. /^0[0-7]+/) [0]. .TP +.BI "-G " INT +Do not output alignments with all bits set in +.I INT +present in the FLAG field. This is the opposite of \fI-f\fR such +that \fI-f12 -G12\fR is the same as no filtering at all. +.I INT +can be specified in hex by beginning with `0x' (i.e. /^0x[0-9A-F]+/) +or in octal by beginning with `0' (i.e. /^0[0-7]+/) [0]. +.TP .BI "-x " STR Read tag to exclude from output (repeatable) [null] .TP From 8b779df31fbcea305923e293070d4e41233070eb Mon Sep 17 00:00:00 2001 From: James Bonfield Date: Fri, 31 Mar 2017 09:49:09 +0100 Subject: [PATCH 14/20] Another change to view/bam2fq filtering wording. The previous attempt had "unset" in the wrong place possibly, but it's now been reworded in terms of flags instead of bits (potentially permitting the symbolic flag characters again in the future). --- sam_view.c | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/sam_view.c b/sam_view.c index 2cf7ae894..d30b5a18b 100644 --- a/sam_view.c +++ b/sam_view.c @@ -574,9 +574,9 @@ static int usage(FILE *fp, int exit_status, int is_long_help) " -l STR only include reads in library STR [null]\n" " -m INT only include reads with number of CIGAR operations consuming\n" " query sequence >= INT [0]\n" -" -f INT only include reads with all of the bits set in INT set in FLAG [0]\n" -" -F INT only include reads with all of the bits unset in INT set in FLAG [0]\n" -" -G INT only include reads with any of the bits unset in INT set in FLAG [0]\n" +" -f INT only include reads with all of the FLAGs in INT present [0]\n" // F&x == x +" -F INT only include reads with none of the FLAGS in INT present [0]\n" // F&x == 0 +" -G INT only EXCLUDE reads with all of the FLAGs in INT present [0]\n" // !(F&x == x) " -s FLOAT subsample reads (given INT.FRAC option value, 0.FRAC is the\n" " fraction of templates/read pairs to keep; INT part sets seed)\n" // read processing @@ -658,9 +658,9 @@ static void bam2fq_usage(FILE *to, const char *command) " -0 FILE write paired reads flagged both or neither READ1 and READ2 to FILE\n" " -1 FILE write paired reads flagged READ1 to FILE\n" " -2 FILE write paired reads flagged READ2 to FILE\n" -" -f INT only include reads with all of the bits set in INT set in FLAG [0]\n" -" -F INT only include reads with all of the bits unset in INT set in FLAG [0]\n" -" -G INT only include reads with any of the bits unset in INT set in FLAG [0]\n" +" -f INT only include reads with all of the FLAGs in INT present [0]\n" // F&x == x +" -F INT only include reads with none of the FLAGS in INT present [0]\n" // F&x == 0 +" -G INT only EXCLUDE reads with all of the FLAGs in INT present [0]\n" // !(F&x == x) " -n don't append /1 and /2 to the read name\n"); if (fq) fprintf(to, " -O output quality in the OQ tag if present\n"); From 112952092efaa09678470503de751832ba28651f Mon Sep 17 00:00:00 2001 From: James Bonfield Date: Fri, 31 Mar 2017 12:30:13 +0100 Subject: [PATCH 15/20] Improvements for C99/XOPEN compliance. Replaced alloca with malloc. Avoid arithmetic on void*. Ensure unsigned types for bit fields. Use of strings.h for strcasecmp. --- Makefile | 1 + bam2bcf_indel.c | 14 ++++++++++---- bam_cat.c | 3 ++- bam_plcmd.c | 1 + bam_reheader.c | 4 ++-- sam.h | 2 +- sam_view.c | 1 + stats.c | 2 +- 8 files changed, 19 insertions(+), 9 deletions(-) diff --git a/Makefile b/Makefile index a38e0c056..e7629a061 100644 --- a/Makefile +++ b/Makefile @@ -24,6 +24,7 @@ CC = gcc AR = ar CPPFLAGS = +#CFLAGS = -g -Wall -O2 -pedantic -std=c99 -D_XOPEN_SOURCE=600 CFLAGS = -g -Wall -O2 LDFLAGS = LIBS = diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c index 2b42b09d7..9749d5bb7 100644 --- a/bam2bcf_indel.c +++ b/bam2bcf_indel.c @@ -439,10 +439,13 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla } free(ref2); free(query); { // compute indelQ - int *sc, tmp, *sumq; - sc = alloca(n_types * sizeof(int)); - sumq = alloca(n_types * sizeof(int)); - memset(sumq, 0, sizeof(int) * n_types); + int sc_a[16], sumq_a[16]; + int tmp, *sc = sc_a, *sumq = sumq_a; + if (n_types > 16) { + sc = (int *)malloc(n_types * sizeof(int)); + sumq = (int *)malloc(n_types * sizeof(int)); + } + memset(sumq, 0, n_types * sizeof(int)); for (s = K = 0; s < n; ++s) { for (i = 0; i < n_plp[s]; ++i, ++K) { bam_pileup1_t *p = plp[s] + i; @@ -523,6 +526,9 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla //fprintf(stderr, "X pos=%d read=%d:%d name=%s call=%d type=%d seqQ=%d indelQ=%d\n", pos, s, i, bam1_qname(p->b), (p->aux>>16)&0x3f, bca->indel_types[(p->aux>>16)&0x3f], (p->aux>>8)&0xff, p->aux&0xff); } } + + if (sc != sc_a) free(sc); + if (sumq != sumq_a) free(sumq); } free(score1); free(score2); // free diff --git a/bam_cat.c b/bam_cat.c index 0d8940b0d..95498ec36 100644 --- a/bam_cat.c +++ b/bam_cat.c @@ -40,6 +40,7 @@ Illumina. #include #include #include +#include #include "htslib/bgzf.h" #include "htslib/sam.h" @@ -468,7 +469,7 @@ int bam_cat(int nfn, char * const *fn, const bam_hdr_t *h, const char* outbam) } if (in->block_offset < in->block_length) { - if (bgzf_write(fp, in->uncompressed_block + in->block_offset, in->block_length - in->block_offset) < 0) goto write_fail; + if (bgzf_write(fp, (char *)in->uncompressed_block + in->block_offset, in->block_length - in->block_offset) < 0) goto write_fail; if (bgzf_flush(fp) != 0) goto write_fail; } diff --git a/bam_plcmd.c b/bam_plcmd.c index 42fc096ee..d17e9d672 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -31,6 +31,7 @@ DEALINGS IN THE SOFTWARE. */ #include #include #include +#include #include #include #include diff --git a/bam_reheader.c b/bam_reheader.c index 0469c0671..acaebd409 100644 --- a/bam_reheader.c +++ b/bam_reheader.c @@ -91,7 +91,7 @@ int bam_reheader(BGZF *in, bam_hdr_t *h, int fd, goto fail; } if (in->block_offset < in->block_length) { - if (bgzf_write(fp, in->uncompressed_block + in->block_offset, in->block_length - in->block_offset) < 0) goto write_fail; + if (bgzf_write(fp, (char *)in->uncompressed_block + in->block_offset, in->block_length - in->block_offset) < 0) goto write_fail; if (bgzf_flush(fp) < 0) goto write_fail; } while ((len = bgzf_raw_read(in, buf, BUF_SIZE)) > 0) { @@ -246,7 +246,7 @@ int cram_reheader_inplace2(cram_fd *fd, const bam_hdr_t *h, const char *arg_list int32_put_blk(b, header_len); cram_block_append(b, sam_hdr_str(hdr), header_len); // Zero the remaining block - memset(cram_block_get_data(b)+cram_block_get_offset(b), 0, + memset((char *)cram_block_get_data(b)+cram_block_get_offset(b), 0, cram_block_get_uncomp_size(b) - cram_block_get_offset(b)); // Make sure all sizes and byte-offsets are consistent after memset cram_block_set_offset(b, cram_block_get_uncomp_size(b)); diff --git a/sam.h b/sam.h index 513010512..6545e640b 100644 --- a/sam.h +++ b/sam.h @@ -50,7 +50,7 @@ typedef struct { samFile *file; struct { BGZF *bam; } x; // Hack so that fp->x.bam still works bam_hdr_t *header; - short is_write:1; + unsigned short is_write:1; } samfile_t; #ifdef __cplusplus diff --git a/sam_view.c b/sam_view.c index d30b5a18b..62bb8fcdc 100644 --- a/sam_view.c +++ b/sam_view.c @@ -27,6 +27,7 @@ DEALINGS IN THE SOFTWARE. */ #include #include +#include #include #include #include diff --git a/stats.c b/stats.c index e377fba7e..35574ed83 100644 --- a/stats.c +++ b/stats.c @@ -1263,7 +1263,7 @@ void init_regions(stats_t *stats, const char *file) stats->regions[tid].pos = realloc(stats->regions[tid].pos,sizeof(pos_t)*stats->regions[tid].mpos); } - if ( (sscanf(&line.s[i+1],"%d %d",&stats->regions[tid].pos[npos].from,&stats->regions[tid].pos[npos].to))!=2 ) error("Could not parse the region [%s]\n", &line.s[i+1]); + if ( (sscanf(&line.s[i+1],"%u %u",&stats->regions[tid].pos[npos].from,&stats->regions[tid].pos[npos].to))!=2 ) error("Could not parse the region [%s]\n", &line.s[i+1]); if ( prev_tid==-1 || prev_tid!=tid ) { prev_tid = tid; From 64852009960cb863b598ade1eb74f6bf3c3fa470 Mon Sep 17 00:00:00 2001 From: James Bonfield Date: Fri, 7 Apr 2017 09:15:23 +0100 Subject: [PATCH 16/20] Switch to -lpthread instead of -pthread for linking. See https://github.com/samtools/htslib/pull/255 for analogous htslib change. --- Makefile | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Makefile b/Makefile index e7629a061..5bd1ddf8e 100644 --- a/Makefile +++ b/Makefile @@ -143,7 +143,7 @@ libbam.a:$(LOBJS) $(AR) -csru $@ $(LOBJS) samtools: $(AOBJS) libbam.a libst.a $(HTSLIB) - $(CC) -pthread $(ALL_LDFLAGS) -o $@ $(AOBJS) libbam.a libst.a $(HTSLIB_LIB) $(CURSES_LIB) -lm $(ALL_LIBS) + $(CC) $(ALL_LDFLAGS) -o $@ $(AOBJS) libbam.a libst.a $(HTSLIB_LIB) $(CURSES_LIB) -lm $(ALL_LIBS) -lpthread # For building samtools and its test suite only: NOT to be installed. libst.a: $(LIBST_OBJS) @@ -224,28 +224,28 @@ check test: samtools $(BGZIP) $(TEST_PROGRAMS) test/merge/test_bam_translate: test/merge/test_bam_translate.o test/test.o libst.a $(HTSLIB) - $(CC) -pthread $(ALL_LDFLAGS) -o $@ test/merge/test_bam_translate.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) + $(CC) $(ALL_LDFLAGS) -o $@ test/merge/test_bam_translate.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread test/merge/test_rtrans_build: test/merge/test_rtrans_build.o libst.a $(HTSLIB) - $(CC) -pthread $(ALL_LDFLAGS) -o $@ test/merge/test_rtrans_build.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) + $(CC) $(ALL_LDFLAGS) -o $@ test/merge/test_rtrans_build.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread test/merge/test_trans_tbl_init: test/merge/test_trans_tbl_init.o libst.a $(HTSLIB) - $(CC) -pthread $(ALL_LDFLAGS) -o $@ test/merge/test_trans_tbl_init.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) + $(CC) $(ALL_LDFLAGS) -o $@ test/merge/test_trans_tbl_init.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread test/split/test_count_rg: test/split/test_count_rg.o test/test.o libst.a $(HTSLIB) - $(CC) -pthread $(ALL_LDFLAGS) -o $@ test/split/test_count_rg.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) + $(CC) $(ALL_LDFLAGS) -o $@ test/split/test_count_rg.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread test/split/test_expand_format_string: test/split/test_expand_format_string.o test/test.o libst.a $(HTSLIB) - $(CC) -pthread $(ALL_LDFLAGS) -o $@ test/split/test_expand_format_string.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) + $(CC) $(ALL_LDFLAGS) -o $@ test/split/test_expand_format_string.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread test/split/test_filter_header_rg: test/split/test_filter_header_rg.o test/test.o libst.a $(HTSLIB) - $(CC) -pthread $(ALL_LDFLAGS) -o $@ test/split/test_filter_header_rg.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) + $(CC) $(ALL_LDFLAGS) -o $@ test/split/test_filter_header_rg.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread test/split/test_parse_args: test/split/test_parse_args.o test/test.o libst.a $(HTSLIB) - $(CC) -pthread $(ALL_LDFLAGS) -o $@ test/split/test_parse_args.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) + $(CC) $(ALL_LDFLAGS) -o $@ test/split/test_parse_args.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread test/vcf-miniview: test/vcf-miniview.o $(HTSLIB) - $(CC) -pthread $(ALL_LDFLAGS) -o $@ test/vcf-miniview.o $(HTSLIB_LIB) $(ALL_LIBS) + $(CC) $(ALL_LDFLAGS) -o $@ test/vcf-miniview.o $(HTSLIB_LIB) $(ALL_LIBS) -lpthread test_test_h = test/test.h $(htslib_sam_h) From 6d612fea9cb324b15525e172a63270c32420494a Mon Sep 17 00:00:00 2001 From: jenniferliddle Date: Wed, 5 Apr 2017 11:50:09 +0100 Subject: [PATCH 17/20] Added option to produce fastq files from BC or other tags jkb: Removed xstr() and ifmt strdup. --- NEWS | 6 + sam_view.c | 409 +++++++++++++++++++++++++------ samtools.1 | 42 ++++ test/bam2fq/5.1.fq.expected | 68 +++++ test/bam2fq/5.2.fq.expected | 68 +++++ test/bam2fq/5.s.fq.expected | 8 + test/bam2fq/bc.fq.expected | 12 + test/bam2fq/bc_split.fq.expected | 12 + test/dat/bam2fq.004.sam | 65 +++++ test/dat/bam2fq.005.sam | 65 +++++ test/test.pl | 5 +- 11 files changed, 679 insertions(+), 81 deletions(-) create mode 100644 test/bam2fq/5.1.fq.expected create mode 100644 test/bam2fq/5.2.fq.expected create mode 100644 test/bam2fq/5.s.fq.expected create mode 100644 test/bam2fq/bc.fq.expected create mode 100644 test/bam2fq/bc_split.fq.expected create mode 100644 test/dat/bam2fq.004.sam create mode 100644 test/dat/bam2fq.005.sam diff --git a/NEWS b/NEWS index 6b4ee27e5..fc3010772 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,9 @@ +Release x.x (the future) +~~~~~~~~~~~~~~~~~~~~~~~~~ + +* Added options to fastq to create fastq files from BC (or other) + tags. + Release 1.4 (13 March 2017) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/sam_view.c b/sam_view.c index 62bb8fcdc..9f3e998c5 100644 --- a/sam_view.c +++ b/sam_view.c @@ -35,6 +35,7 @@ DEALINGS IN THE SOFTWARE. */ #include #include #include +#include #include "htslib/sam.h" #include "htslib/faidx.h" #include "htslib/kstring.h" @@ -42,6 +43,10 @@ DEALINGS IN THE SOFTWARE. */ #include "htslib/thread_pool.h" #include "samtools.h" #include "sam_opts.h" + +#define DEFAULT_BARCODE_TAG "BC" +#define DEFAULT_QUALITY_TAG "QT" + KHASH_SET_INIT_STR(rg) typedef khash_t(rg) *rghash_t; @@ -656,22 +661,37 @@ static void bam2fq_usage(FILE *to, const char *command) "Usage: samtools %s [options...] \n", command); fprintf(to, "Options:\n" -" -0 FILE write paired reads flagged both or neither READ1 and READ2 to FILE\n" -" -1 FILE write paired reads flagged READ1 to FILE\n" -" -2 FILE write paired reads flagged READ2 to FILE\n" -" -f INT only include reads with all of the FLAGs in INT present [0]\n" // F&x == x -" -F INT only include reads with none of the FLAGS in INT present [0]\n" // F&x == 0 -" -G INT only EXCLUDE reads with all of the FLAGs in INT present [0]\n" // !(F&x == x) -" -n don't append /1 and /2 to the read name\n"); +" -0 FILE write paired reads flagged both or neither READ1 and READ2 to FILE\n" +" -1 FILE write paired reads flagged READ1 to FILE\n" +" -2 FILE write paired reads flagged READ2 to FILE\n" +" -f INT only include reads with all of the FLAGs in INT present [0]\n" // F&x == x +" -F INT only include reads with none of the FLAGS in INT present [0]\n" // F&x == 0 +" -G INT only EXCLUDE reads with all of the FLAGs in INT present [0]\n" // !(F&x == x) +" -n don't append /1 and /2 to the read name\n" +" -N always append /1 and /2 to the read name\n"); if (fq) fprintf(to, -" -O output quality in the OQ tag if present\n"); +" -O output quality in the OQ tag if present\n"); fprintf(to, -" -s FILE write singleton reads to FILE [assume single-end]\n" -" -t copy RG, BC and QT tags to the %s header line\n", +" -s FILE write singleton reads to FILE [assume single-end]\n" +" -t copy RG, BC and QT tags to the %s header line\n", fq ? "FASTQ" : "FASTA"); if (fq) fprintf(to, -" -v INT default quality score if not given in file [1]\n"); +" -v INT default quality score if not given in file [1]\n" +" --i1 FILE write first index reads to FILE\n" +" --i2 FILE write second index reads to FILE\n" +" --barcode-tag TAG Barcode tag [default: " DEFAULT_BARCODE_TAG "]\n" +" --quality-tag TAG Quality tag [default: " DEFAULT_QUALITY_TAG "]\n" +" --index-format STR How to parse barcode and quality tags\n\n"); sam_global_opt_help(to, "-.--.@"); + fprintf(to, +" \n" +" The index-format string describes how to parse the barcode and quality tags, for example:\n" +" i14i8 the first 14 characters are index 1, the next 8 characters are index 2\n" +" n8i14 ignore the first 8 characters, and use the next 14 characters for index 1\n" +" If the tag contains a separator, then the numeric part can be replaced with '*' to mean\n" +" 'read until the separator or end of tag', for example:\n" +" n*i* ignore the left part of the tag until the separator, then use the second part\n" +" of the tag as index 1\n"); } typedef enum { READ_UNKNOWN = 0, READ_1 = 1, READ_2 = 2 } readpart; @@ -680,17 +700,22 @@ typedef struct bam2fq_opts { char *fnse; char *fnr[3]; char *fn_input; // pointer to input filename in argv do not free - bool has12, use_oq, copy_tags; + bool has12, has12always, use_oq, copy_tags; int flag_on, flag_off, flag_alloff; sam_global_args ga; fastfile filetype; int def_qual; + char *barcode_tag; + char *quality_tag; + char *index_file[2]; + char *index_format; } bam2fq_opts_t; typedef struct bam2fq_state { samFile *fp; FILE *fpse; FILE *fpr[3]; + FILE *fpi[2]; bam_hdr_t *h; bool has12, use_oq, copy_tags; int flag_on, flag_off, flag_alloff; @@ -698,6 +723,71 @@ typedef struct bam2fq_state { int def_qual; } bam2fq_state_t; +/* + * get and decode the read from a BAM record + * Why on earth isn't this (and the next few functions) part of htslib? + * The whole point of an API is that we should not have to know the internal structure... + */ + +/* + * reverse a string in place + (stolen from http://stackoverflow.com/questions/8534274/is-the-strrev-function-not-available-in-linux) + */ +static char *reverse(char *str) +{ + int i = strlen(str)-1,j=0; + char ch; + while (i>j) { + ch = str[i]; + str[i]= str[j]; + str[j] = ch; + i--; + j++; + } + return str; +} + +/* return the read, reverse complemented if necessary */ +static char *get_read(const bam1_t *rec) +{ + int len = rec->core.l_qseq + 1; + char *read = calloc(1, len); + char *seq = (char *)bam_get_seq(rec); + int n; + + if (!read) return NULL; + + for (n=0; n < rec->core.l_qseq; n++) { + if (rec->core.flag & BAM_FREVERSE) read[n] = seq_nt16_str[seq_comp_table[bam_seqi(seq,n)]]; + else read[n] = seq_nt16_str[bam_seqi(seq,n)]; + } + if (rec->core.flag & BAM_FREVERSE) reverse(read); + return read; +} + +/* + * get and decode the quality from a BAM record + */ +static char *get_quality(const bam1_t *rec) +{ + char *quality = calloc(1, rec->core.l_qseq + 1); + char *q = (char *)bam_get_qual(rec); + int n; + + if (*q == '\xff') { free(quality); return NULL; } + + for (n=0; n < rec->core.l_qseq; n++) { + quality[n] = q[n]+33; + } + if (rec->core.flag & BAM_FREVERSE) reverse(quality); + return quality; +} + +// +// End of htslib complaints +// + + static readpart which_readpart(const bam1_t *b) { if ((b->core.flag & BAM_FREAD1) && !(b->core.flag & BAM_FREAD2)) { @@ -709,85 +799,60 @@ static readpart which_readpart(const bam1_t *b) } } -// Transform a bam1_t record into a string with the FASTQ representation of it -// @returns false for error, true for success -static bool bam1_to_fq(const bam1_t *b, kstring_t *linebuf, const bam2fq_state_t *state) +/* + * parse the length part from the index-format string + */ +static int getLength(char **s) { - int i; - int32_t qlen = b->core.l_qseq; - assert(qlen >= 0); - uint8_t *seq; - uint8_t *qual = bam_get_qual(b); - const uint8_t *oq = NULL; - if (state->use_oq) { - oq = bam_aux_get(b, "OQ"); - if (oq) oq++; // skip tag type + int n = 0; + while (**s) { + if (**s == '*') { n=-1; (*s)++; break; } + if ( !isdigit(**s)) break; + n = n*10 + ((**s)-'0'); + (*s)++; } - bool has_qual = (qual[0] != 0xff || (state->use_oq && oq)); // test if there is quality + return n; +} + +static bool make_fq_line(const bam1_t *rec, char *seq, char *qual, kstring_t *linebuf, const bam2fq_state_t *state) +{ + int i; linebuf->l = 0; // Write read name - readpart readpart = which_readpart(b); kputc(state->filetype == FASTA? '>' : '@', linebuf); - kputs(bam_get_qname(b), linebuf); + kputs(bam_get_qname(rec), linebuf); // Add the /1 /2 if requested if (state->has12) { + readpart readpart = which_readpart(rec); if (readpart == READ_1) kputs("/1", linebuf); else if (readpart == READ_2) kputs("/2", linebuf); } if (state->copy_tags) { for (i = 0; copied_tags[i]; ++i) { uint8_t *s; - if ((s = bam_aux_get(b, copied_tags[i])) != 0) { - kputc('\t', linebuf); - kputsn(copied_tags[i], 2, linebuf); - kputsn(":Z:", 3, linebuf); - kputs(bam_aux2Z(s), linebuf); + if ((s = bam_aux_get(rec, copied_tags[i])) != 0) { + if (*s == 'Z') { + kputc('\t', linebuf); + kputsn(copied_tags[i], 2, linebuf); + kputsn(":Z:", 3, linebuf); + kputs(bam_aux2Z(s), linebuf); + } } } } kputc('\n', linebuf); - - seq = bam_get_seq(b); - - if (b->core.flag & BAM_FREVERSE) { // read is reverse complemented - for (i = qlen-1; i > -1; --i) { - char c = seq_nt16_str[seq_comp_table[bam_seqi(seq,i)]]; - kputc(c, linebuf); - } - } else { - for (i = 0; i < qlen; ++i) { - char c = seq_nt16_str[bam_seqi(seq,i)]; - kputc(c, linebuf); - } - } + kputs(seq, linebuf); kputc('\n', linebuf); if (state->filetype == FASTQ) { // Write quality kputs("+\n", linebuf); - if (has_qual) { - if (state->use_oq && oq) { - if (b->core.flag & BAM_FREVERSE) { // read is reverse complemented - for (i = qlen-1; i > -1; --i) { - kputc(oq[i], linebuf); - } - } else { - kputs((char*)oq, linebuf); - } - } else { - if (b->core.flag & BAM_FREVERSE) { // read is reverse complemented - for (i = qlen-1; i > -1; --i) { - kputc(33 + qual[i], linebuf); - } - } else { - for (i = 0; i < qlen; ++i) { - kputc(33 + qual[i], linebuf); - } - } - } + if (qual && *qual) { + kputs(qual, linebuf); } else { - for (i = 0; i < qlen; ++i) { + int len = strlen(seq); + for (i = 0; i < len; ++i) { kputc(33 + state->def_qual, linebuf); } } @@ -796,23 +861,147 @@ static bool bam1_to_fq(const bam1_t *b, kstring_t *linebuf, const bam2fq_state_t return true; } +/* + * Create FASTQ lines from the barcode tag using the index-format + */ +static bool tags2fq(bam1_t *rec, const bam2fq_state_t *state, const bam2fq_opts_t* opts) +{ + uint8_t *p; + char *ifmt = opts->index_format; + char *tag = NULL; + char *qual = NULL; + int file_number = 0; + kstring_t linebuf = { 0, 0, NULL }; // Buffer + + // read barcode tag + p = bam_aux_get(rec,opts->barcode_tag); + if (p) tag = bam_aux2Z(p); + + if (!tag) return true; // there is no tag + + // read quality tag + p = bam_aux_get(rec, opts->quality_tag); + if (p) qual = bam_aux2Z(p); + + // Parse the index-format string + while (*ifmt) { + if (file_number > 1) break; // shouldn't happen if we've validated paramaters correctly + char action = *ifmt; // should be 'i' or 'n' + ifmt++; // skip over action + int index_len = getLength(&ifmt); + + char *sub_tag = calloc(1, strlen(tag)+1); + char *sub_qual = calloc(1, strlen(tag)+1); + int n = 0; + + if (index_len < 0) { + // read until separator + while (isalpha(*tag)) { + sub_tag[n] = *tag++; + if (qual) sub_qual[n] = *qual++; + n++; + } + if (*tag) { // skip separator + tag++; + if (qual) qual++; + } + } else { + // read index_len characters + while (index_len-- && *tag) { + sub_tag[n] = *tag++; + if (qual) sub_qual[n] = *qual++; + n++; + } + } + + if (action=='i' && *sub_tag && state->fpi[file_number]) { + make_fq_line(rec, sub_tag, sub_qual, &linebuf, state); + fputs(linebuf.s, state->fpi[file_number++]); + } + free(sub_qual); free(sub_tag); + + } + + free(linebuf.s); + return true; +} + +// Transform a bam1_t record into a string with the FASTQ representation of it +// @returns false for error, true for success +static bool bam1_to_fq(const bam1_t *b, kstring_t *linebuf, const bam2fq_state_t *state) +{ + int32_t qlen = b->core.l_qseq; + assert(qlen >= 0); + const uint8_t *oq = NULL; + char *qual = NULL; + + char *seq = get_read(b); + + if (state->use_oq) { + oq = bam_aux_get(b, "OQ"); + if (oq) { + oq++; + qual = strdup(bam_aux2Z(oq)); + if (b->core.flag & BAM_FREVERSE) { // read is reverse complemented + reverse(qual); + } + } + } else { + qual = get_quality(b); + } + + make_fq_line(b, seq, qual, linebuf, state); + + free(qual); + free(seq); + return true; +} + +static void free_opts(bam2fq_opts_t *opts) +{ + free(opts->barcode_tag); + free(opts->quality_tag); + free(opts->index_format); + free(opts); +} + // return true if valid static bool parse_opts(int argc, char *argv[], bam2fq_opts_t** opts_out) { // Parse args bam2fq_opts_t* opts = calloc(1, sizeof(bam2fq_opts_t)); opts->has12 = true; + opts->has12always = false; opts->filetype = FASTQ; opts->def_qual = 1; + opts->barcode_tag = NULL; + opts->quality_tag = NULL; + opts->index_format = NULL; + opts->index_file[0] = NULL; + opts->index_file[1] = NULL; int c; sam_global_args_init(&opts->ga); static const struct option lopts[] = { SAM_OPT_GLOBAL_OPTIONS('-', 0, '-', '-', 0, '@'), + {"i1", required_argument, NULL, 1}, + {"I1", required_argument, NULL, 1}, + {"i2", required_argument, NULL, 2}, + {"I2", required_argument, NULL, 2}, + {"if", required_argument, NULL, 3}, + {"IF", required_argument, NULL, 3}, + {"index-format", required_argument, NULL, 3}, + {"barcode-tag", required_argument, NULL, 'b'}, + {"quality-tag", required_argument, NULL, 'q'}, { NULL, 0, NULL, 0 } }; - while ((c = getopt_long(argc, argv, "0:1:2:f:F:G:nOs:tv:@:", lopts, NULL)) > 0) { + while ((c = getopt_long(argc, argv, "0:1:2:f:F:G:nNOs:tv:@:", lopts, NULL)) > 0) { switch (c) { + case 'b': opts->barcode_tag = strdup(optarg); break; + case 'q': opts->quality_tag = strdup(optarg); break; + case 1 : opts->index_file[0] = optarg; break; + case 2 : opts->index_file[1] = optarg; break; + case 3 : opts->index_format = strdup(optarg); break; case '0': opts->fnr[0] = optarg; break; case '1': opts->fnr[1] = optarg; break; case '2': opts->fnr[2] = optarg; break; @@ -820,26 +1009,66 @@ static bool parse_opts(int argc, char *argv[], bam2fq_opts_t** opts_out) case 'F': opts->flag_off |= strtol(optarg, 0, 0); break; case 'G': opts->flag_alloff |= strtol(optarg, 0, 0); break; case 'n': opts->has12 = false; break; + case 'N': opts->has12always = true; break; case 'O': opts->use_oq = true; break; case 's': opts->fnse = optarg; break; case 't': opts->copy_tags = true; break; case 'v': opts->def_qual = atoi(optarg); break; - case '?': bam2fq_usage(stderr, argv[0]); free(opts); return false; + case '?': bam2fq_usage(stderr, argv[0]); free_opts(opts); return false; default: if (parse_sam_global_opt(c, optarg, lopts, &opts->ga) != 0) { - bam2fq_usage(stderr, argv[0]); free(opts); return false; + bam2fq_usage(stderr, argv[0]); free_opts(opts); return false; } break; } } if (opts->fnr[1] || opts->fnr[2]) opts->has12 = false; + if (opts->has12always) opts->has12 = true; + + if (!opts->barcode_tag) opts->barcode_tag = strdup(DEFAULT_BARCODE_TAG); + if (!opts->quality_tag) opts->quality_tag = strdup(DEFAULT_QUALITY_TAG); + + int nIndex = 0; + if (opts->index_format) { + char *s; + for (s = opts->index_format; *s; s++) { + if (*s == 'i') nIndex++; + } + } + if (nIndex>2) { + fprintf(stderr,"Invalid index format: more than 2 indexes\n"); + bam2fq_usage(stderr, argv[0]); + free_opts(opts); + return false; + } + + if (opts->index_file[1] && !opts->index_file[0]) { + fprintf(stderr, "Index one specified, but index two not given\n"); + bam2fq_usage(stderr, argv[0]); + free_opts(opts); + return false; + } + + if (nIndex==2 && !opts->index_file[1]) { + fprintf(stderr, "index_format specifies two indexes, but only one index file given\n"); + bam2fq_usage(stderr, argv[0]); + free_opts(opts); + return false; + } + + if (nIndex==1 && !opts->index_file[0]) { + fprintf(stderr, "index_format specifies an index, but no index file given\n"); + bam2fq_usage(stderr, argv[0]); + free_opts(opts); + return false; + } if (opts->def_qual < 0 || 93 < opts->def_qual) { fprintf(stderr, "Invalid -v default quality %i, allowed range 0 to 93\n", opts->def_qual); bam2fq_usage(stderr, argv[0]); - free(opts); - return true; + free_opts(opts); + return false; } const char* type_str = argv[0]; @@ -850,20 +1079,21 @@ static bool parse_opts(int argc, char *argv[], bam2fq_opts_t** opts_out) } else { print_error("bam2fq", "Unrecognised type call \"%s\", this should be impossible... but you managed it!", type_str); bam2fq_usage(stderr, argv[0]); - free(opts); + free_opts(opts); return false; } if ((argc - (optind)) == 0) { + fprintf(stderr, "No input file specified.\n"); bam2fq_usage(stdout, argv[0]); - free(opts); + free_opts(opts); return false; } if ((argc - (optind)) != 1) { fprintf(stderr, "Too many arguments.\n"); bam2fq_usage(stderr, argv[0]); - free(opts); + free_opts(opts); return false; } opts->fn_input = argv[optind]; @@ -925,6 +1155,17 @@ static bool init_state(const bam2fq_opts_t* opts, bam2fq_state_t** state_out) state->fpr[i] = stdout; } } + for (i = 0; i < 2; i++) { + state->fpi[i] = NULL; + if (opts->index_file[i]) { + state->fpi[i] = fopen(opts->index_file[i], "w"); + if (state->fpi[i] == NULL) { + print_error_errno("bam2fq", "Cannot write to i%d file \"%s\"", i+1, opts->index_file[i]); + free(state); + return false; + } + } + } state->h = sam_hdr_read(state->fp); if (state->h == NULL) { @@ -947,6 +1188,12 @@ static bool destroy_state(const bam2fq_opts_t *opts, bam2fq_state_t *state, int* for (i = 0; i < 3; ++i) { if (state->fpr[i] != stdout && fclose(state->fpr[i])) { print_error_errno("bam2fq", "Error closing r%d file \"%s\"", i, opts->fnr[i]); valid = false; } } + for (i = 0; i < 2; i++) { + if (state->fpi[i] && fclose(state->fpi[i])) { + print_error_errno("bam2fq", "Error closing i%d file \"%s\"", i+1, opts->index_file[i]); + valid = false; + } + } free(state); return valid; } @@ -960,7 +1207,7 @@ static inline bool filter_it_out(const bam1_t *b, const bam2fq_state_t *state) } -static bool bam2fq_mainloop_singletontrack(bam2fq_state_t *state) +static bool bam2fq_mainloop_singletontrack(bam2fq_state_t *state, bam2fq_opts_t* opts) { bam1_t* b = bam_init1(); char *current_qname = NULL; @@ -1016,6 +1263,7 @@ static bool bam2fq_mainloop_singletontrack(bam2fq_state_t *state) return false; } score[which_readpart(b)] = b_score; + if (state->fpi[0]) tags2fq(b, state, opts); } } if (!valid) @@ -1033,7 +1281,7 @@ static bool bam2fq_mainloop_singletontrack(bam2fq_state_t *state) return valid; } -static bool bam2fq_mainloop(bam2fq_state_t *state) +static bool bam2fq_mainloop(bam2fq_state_t *state, bam2fq_opts_t* opts) { // process a name collated BAM into fastq bam1_t* b = bam_init1(); @@ -1049,6 +1297,7 @@ static bool bam2fq_mainloop(bam2fq_state_t *state) if (!bam1_to_fq(b, &linebuf, state)) return false; fputs(linebuf.s, state->fpr[which_readpart(b)]); + if (state->fpi[0]) tags2fq(b, state, opts); } free(linebuf.s); bam_destroy1(b); @@ -1069,14 +1318,14 @@ int main_bam2fq(int argc, char *argv[]) if (!init_state(opts, &state)) return EXIT_FAILURE; if (state->fpse) { - if (!bam2fq_mainloop_singletontrack(state)) status = EXIT_FAILURE; + if (!bam2fq_mainloop_singletontrack(state,opts)) status = EXIT_FAILURE; } else { - if (!bam2fq_mainloop(state)) status = EXIT_FAILURE; + if (!bam2fq_mainloop(state,opts)) status = EXIT_FAILURE; } if (!destroy_state(opts, state, &status)) return EXIT_FAILURE; sam_global_args_free(&opts->ga); - free(opts); + free_opts(opts); return status; } diff --git a/samtools.1 b/samtools.1 index 8b8e49dfb..594340b04 100644 --- a/samtools.1 +++ b/samtools.1 @@ -1389,6 +1389,10 @@ Using .B -n causes read names to be left as they are. .TP 8 +.B -N +Always add either '/1' or '/2' to the end of read names +even when put into different files. +.TP 8 .B -O Use quality values from OQ tags in preference to standard quality string if available. @@ -1424,6 +1428,44 @@ present in the FLAG field. .I INT can be specified in hex by beginning with `0x' (i.e. /^0x[0-9A-F]+/) or in octal by beginning with `0' (i.e. /^0[0-7]+/) [0]. +.TP 8 +.BI "-G " INT +Only EXCLUDE reads with all of the bits set in +.I INT +present in the FLAG field. +.I INT +can be specified in hex by beginning with `0x' (i.e. /^0x[0-9A-F]+/) +or in octal by beginning with `0' (i.e. /^0[0-7]+/) [0]. +.TP 8 +.B --i1 FILE +write first index reads to FILE +.TP 8 +.B --i2 FILE +write second index reads to FILE +.TP 8 +.B --barcode-tag TAG +aux tag to find index reads in [default: BC] +.TP 8 +.B --quality-tag TAG +aux tag to find index quality in [default: QT] +.TP 8 +.B --index-format STR +string to describe how to parse the barcode and quality tags. For example: + +.RS +.TP 8 +.B i14i8 +the first 14 characters are index 1, the next 8 characters are index 2 +.TP 8 +.B n8i14 +ignore the first 8 characters, and use the next 14 characters for index 1 + +If the tag contains a separator, then the numeric part can be replaced with '*' to +mean 'read until the separator or end of tag', for example: +.TP 8 +.B n*i* +ignore the left part of the tag until the separator, then use the second part +.RE .RE .TP \"-------- collate diff --git a/test/bam2fq/5.1.fq.expected b/test/bam2fq/5.1.fq.expected new file mode 100644 index 000000000..97fa35b87 --- /dev/null +++ b/test/bam2fq/5.1.fq.expected @@ -0,0 +1,68 @@ +@ref1_grp1_p001/1 +CGAGCTCGGT ++ +!!!!!!!!!! +@ref1_grp1_p002/1 +CTCGGTACCC ++ +########## +@ref1_grp1_p003/1 +GTACCCGGGG ++ +%%%%%%%%%% +@ref1_grp1_p004/1 +CCGGGGATCC ++ +'''''''''' +@ref1_grp1_p005/1 +GGATCCTCTA ++ +)))))))))) +@ref1_grp1_p006/1 +CCTCTAGAGT ++ +++++++++++ +@ref1_grp2_p001/1 +AGCTCGGTAC ++ +"""""""""" +@ref1_grp2_p002/1 +CGGTACCCGG ++ +$$$$$$$$$$ +@ref1_grp2_p003/1 +ACCCGGGGAT ++ +&&&&&&&&&& +@ref1_grp2_p004/1 +GGGGATCCTC ++ +(((((((((( +@ref1_grp2_p005/1 +ATCCTCTAGA ++ +********** +@ref1_grp2_p006/1 +TCTAGAGTCG ++ +,,,,,,,,,, +@ref2_grp3_p001/1 +GTGACACTATAGAAT ++ +~~~~~~~~~~~~~~~ +@ref2_grp3_p002/1 +CTGTTTCCTGTGTGA ++ +{{{{{{{{{{{{{{{ +@ref2_grp3_p003/1 +ACGTMRWSYKVHDBN ++ +0123456789abcd! +@ref12_grp1_p001/1 +TGCAGGCATG ++ +AAAAAAAAAA +@ref12_grp2_p001/1 +CAAGCTTGAG ++ +AAAAAAAAAA diff --git a/test/bam2fq/5.2.fq.expected b/test/bam2fq/5.2.fq.expected new file mode 100644 index 000000000..edbf2e9eb --- /dev/null +++ b/test/bam2fq/5.2.fq.expected @@ -0,0 +1,68 @@ +@ref1_grp1_p001/2 +GTCGACTCTA ++ +---------- +@ref1_grp1_p002/2 +GCAGGTCGAC ++ +////////// +@ref1_grp1_p003/2 +GCCTGCAGGT ++ +1111111111 +@ref1_grp1_p004/2 +GCATGCCTGC ++ +3333333333 +@ref1_grp1_p005/2 +GCTTGCATGC ++ +5555555555 +@ref1_grp1_p006/2 +TCAAGCTTGC ++ +7777777777 +@ref1_grp2_p001/2 +AGGTCGACTC ++ +.......... +@ref1_grp2_p002/2 +CTGCAGGTCG ++ +0000000000 +@ref1_grp2_p003/2 +ATGCCTGCAG ++ +2222222222 +@ref1_grp2_p004/2 +TTGCATGCCT ++ +4444444444 +@ref1_grp2_p005/2 +AAGCTTGCAT ++ +6666666666 +@ref1_grp2_p006/2 +ACTCAAGCTT ++ +8888888888 +@ref2_grp3_p001/2 +CTGTTTCCTGTGTGA ++ +||||||||||||||| +@ref2_grp3_p002/2 +CGCCAAGCTATTTAG ++ +}}}}}}}}}}}}}}} +@ref2_grp3_p003/2 +ACGTMRWSYKVHDBN ++ +0123456789abcd! +@ref12_grp1_p001/2 +CACTATAGAA ++ +BBBBBBBBBB +@ref12_grp2_p001/2 +ATTTAGGTGA ++ +BBBBBBBBBB diff --git a/test/bam2fq/5.s.fq.expected b/test/bam2fq/5.s.fq.expected new file mode 100644 index 000000000..58fb091bc --- /dev/null +++ b/test/bam2fq/5.s.fq.expected @@ -0,0 +1,8 @@ +@ref1_grp2_p002a/1 +CGGTACCCGG ++ +$$$$$$$$$$ +@unaligned_grp3_p001/1 +CACTCGTTCATGACG ++ +0123456789abcde diff --git a/test/bam2fq/bc.fq.expected b/test/bam2fq/bc.fq.expected new file mode 100644 index 000000000..e2efa367a --- /dev/null +++ b/test/bam2fq/bc.fq.expected @@ -0,0 +1,12 @@ +@ref1_grp1_p001 +GT ++ +"" +@ref1_grp1_p002 +TT ++ +"" +@ref1_grp2_p001 +CA ++ +cd diff --git a/test/bam2fq/bc_split.fq.expected b/test/bam2fq/bc_split.fq.expected new file mode 100644 index 000000000..fc9935905 --- /dev/null +++ b/test/bam2fq/bc_split.fq.expected @@ -0,0 +1,12 @@ +@ref1_grp1_p001/1 +GT ++ +"" +@ref1_grp1_p002/1 +CCGG ++ +"""" +@ref1_grp2_p001/1 +CA ++ +cd diff --git a/test/dat/bam2fq.004.sam b/test/dat/bam2fq.004.sam new file mode 100644 index 000000000..aafbd1ba8 --- /dev/null +++ b/test/dat/bam2fq.004.sam @@ -0,0 +1,65 @@ +@HD VN:1.4 SO:queryname +@RG ID:grp1 DS:Group 1 LB:Library 1 SM:Sample +@RG ID:grp2 DS:Group 2 LB:Library 2 SM:Sample +@RG ID:grp3 DS:Group 3 LB:Library 3 SM:Sample +@PG ID:prog1 PN:emacs CL:emacs VN:23.1.1 +@CO The MIT License +@CO +@CO Copyright (c) 2014, 2015 Genome Research Ltd. +@CO +@CO Permission is hereby granted, free of charge, to any person obtaining a copy +@CO of this software and associated documentation files (the "Software"), to deal +@CO in the Software without restriction, including without limitation the rights +@CO to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +@CO copies of the Software, and to permit persons to whom the Software is +@CO furnished to do so, subject to the following conditions: +@CO +@CO The above copyright notice and this permission notice shall be included in +@CO all copies or substantial portions of the Software. +@CO +@CO THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +@CO IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +@CO FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +@CO AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +@CO LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +@CO OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +@CO THE SOFTWARE. +@SQ SN:ref1 LN:56 M5:08c04d512d4797d9ba2a156c1daba468 +@SQ SN:ref2 LN:60 M5:7c35feac7036c1cdef3bee0cc4b21437 +ref1_grp1_p001 99 ref1 1 0 10M = 25 34 CGAGCTCGGT !!!!!!!!!! MD:Z:10 NM:i:0 RG:Z:grp1 BC:Z:ACGT H0:i:1 aa:A:! ab:A:~ fa:f:3.14159 za:Z:Hello world! ha:H:DEADBEEF ba:B:c,-128,0,127 bb:B:C,0,127,255 bc:B:s,-32768,0,32767 bd:B:S,0,32768,65535 be:B:i,-2147483648,0,2147483647 bf:B:I,0,2147483648,4294967295 bg:B:f,2.71828,6.626e-34,2.9979e+09 +ref1_grp1_p001 147 ref1 25 12 10M = 1 -34 TAGAGTCGAC ---------- MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp1_p002 99 ref1 5 2 10M = 29 34 CTCGGTACCC ########## MD:Z:10 NM:i:0 RG:Z:grp1 BC:Z:AATTCCGG H0:i:1 aa:A:a ab:A:z fa:f:4.3597e-18 za:Z:Another string ha:H:2000AD +ref1_grp1_p002 147 ref1 29 14 10M = 5 -34 GTCGACCTGC ////////// MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp1_p003 99 ref1 9 4 10M = 33 34 GTACCCGGGG %%%%%%%%%% MD:Z:10 NM:i:0 fa:f:1.66e-27 RG:Z:grp1 +ref1_grp1_p003 147 ref1 33 16 10M = 9 -34 ACCTGCAGGC 1111111111 MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp1_p004 99 ref1 13 6 10M = 37 34 CCGGGGATCC '''''''''' MD:Z:10 NM:i:0 fa:f:1.38e-23 za:Z:xRG:Z:grp2 RG:Z:grp1 +ref1_grp1_p004 147 ref1 37 18 10M = 13 -34 GCAGGCATGC 3333333333 MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp1_p005 99 ref1 17 8 10M = 41 34 GGATCCTCTA )))))))))) MD:Z:10 NM:i:0 RG:Z:grp1 ia:i:40000 +ref1_grp1_p005 147 ref1 41 20 10M = 17 -34 GCATGCAAGC 5555555555 MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp1_p006 99 ref1 21 10 10M = 45 34 CCTCTAGAGT ++++++++++ MD:Z:10 NM:i:0 RG:Z:grp1 ia:i:255 +ref1_grp1_p006 147 ref1 45 22 10M = 21 -34 GCAAGCTTGA 7777777777 MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp2_p001 99 ref1 3 1 8M2S = 27 34 AGCTCGGTAC """""""""" MD:Z:8 NM:i:0 RG:Z:grp2 BC:Z:TGCA QT:Z:abcd H0:i:1 aa:A:A ab:A:Z fa:f:6.67e-11 za:Z:!"$%^&*() ha:H:CAFE +ref1_grp2_p001 2147 ref1 11 1 8H2M = 27 34 AC "" MD:Z:2 NM:i:0 RG:Z:grp2 BC:Z:TGCA H0:i:1 aa:A:A ab:A:Z fa:f:6.67e-11 za:Z:!"$%^&*() ha:H:CAFE +ref1_grp2_p001 147 ref1 27 13 10M = 3 -34 GAGTCGACCT .......... MD:Z:10 NM:i:0 RG:Z:grp2 +ref1_grp2_p002a 99 ref1 7 3 10M = 31 34 CGGTACCCGG $$$$$$$$$$ MD:Z:10 NM:i:0 fa:f:6.022e+23 RG:Z:grp2 +ref1_grp2_p002 99 ref1 7 3 10M = 31 34 CGGTACCCGG $$$$$$$$$$ MD:Z:10 NM:i:0 fa:f:6.022e+23 RG:Z:grp2 +ref1_grp2_p002 147 ref1 31 15 10M = 7 -34 CGACCTGCAG 0000000000 MD:Z:10 NM:i:0 RG:Z:grp2 +ref1_grp2_p003 99 ref1 11 5 10M = 35 34 ACCCGGGGAT &&&&&&&&&& MD:Z:10 NM:i:0 RG:Z:grp2 ia:i:4294967295 +ref1_grp2_p003 147 ref1 35 17 10M = 11 -34 CTGCAGGCAT 2222222222 MD:Z:10 NM:i:0 RG:Z:grp2 +ref1_grp2_p004 99 ref1 15 7 10M = 39 34 GGGGATCCTC (((((((((( MD:Z:10 NM:i:0 RG:Z:grp2 ia:i:-2147483648 +ref1_grp2_p004 147 ref1 39 19 10M = 15 -34 AGGCATGCAA 4444444444 MD:Z:10 NM:i:0 RG:Z:grp2 +ref1_grp2_p005 99 ref1 19 9 10M = 43 34 ATCCTCTAGA ********** MD:Z:10 NM:i:0 RG:Z:grp2 ia:i:-1000 +ref1_grp2_p005 147 ref1 43 21 10M = 19 -34 ATGCAAGCTT 6666666666 MD:Z:10 NM:i:0 RG:Z:grp2 +ref1_grp2_p006 99 ref1 23 11 10M = 47 34 TCTAGAGTCG ,,,,,,,,,, MD:Z:10 NM:i:0 RG:Z:grp2 ia:i:-1 +ref1_grp2_p006 147 ref1 47 23 10M = 23 -34 AAGCTTGAGT 8888888888 MD:Z:10 NM:i:0 RG:Z:grp2 +ref2_grp3_p001 83 ref2 1 99 15M = 31 45 ATTCTATAGTGTCAC ~~~~~~~~~~~~~~~ MD:Z:15 NM:i:0 RG:Z:grp3 +ref2_grp3_p001 163 ref2 31 99 15M = 1 -45 CTGTTTCCTGTGTGA ||||||||||||||| MD:Z:0T0A0A1C0A0T0G0G0T0C0A1A0G0 NM:i:13 RG:Z:grp3 +ref2_grp3_p002 99 ref2 46 99 15M = 16 -45 CTGTTTCCTGTGTGA {{{{{{{{{{{{{{{ MD:Z:15 NM:i:0 RG:Z:grp3 +ref2_grp3_p002 147 ref2 16 99 15M = 46 45 CTAAATAGCTTGGCG }}}}}}}}}}}}}}} MD:Z:15 NM:i:0 RG:Z:grp3 +ref2_grp3_p003 99 ref2 1 99 15M = 1 15 ACGTMRWSYKVHDBN 0123456789abcd! MD:Z:1T0T0C0T0A0T0A0G0T0G0T0C0A0C0 NM:i:14 RG:Z:grp3 +ref2_grp3_p003 147 ref2 1 99 15M = 1 -15 NVHDBMRSWYKACGT !dcba9876543210 MD:Z:0A0T0T0C0T0A0T0A0G0T0G0T1A0C0 NM:i:14 RG:Z:grp3 +ref12_grp1_p001 97 ref1 36 50 10M ref2 2 0 TGCAGGCATG AAAAAAAAAA MD:Z:10 NM:i:0 RG:Z:grp1 +ref12_grp1_p001 145 ref2 2 50 10M ref1 36 0 TTCTATAGTG BBBBBBBBBB MD:Z:10 NM:i:0 RG:Z:grp1 +ref12_grp2_p001 97 ref1 46 50 10M ref2 12 0 CAAGCTTGAG AAAAAAAAAA MD:Z:10 NM:i:0 RG:Z:grp2 +ref12_grp2_p001 145 ref2 12 50 10M ref1 46 0 TCACCTAAAT BBBBBBBBBB MD:Z:10 NM:i:0 RG:Z:grp2 +unaligned_grp3_p001 77 * 0 0 * * 0 0 CACTCGTTCATGACG 0123456789abcde RG:Z:grp3 diff --git a/test/dat/bam2fq.005.sam b/test/dat/bam2fq.005.sam new file mode 100644 index 000000000..0e0b05b64 --- /dev/null +++ b/test/dat/bam2fq.005.sam @@ -0,0 +1,65 @@ +@HD VN:1.4 SO:queryname +@RG ID:grp1 DS:Group 1 LB:Library 1 SM:Sample +@RG ID:grp2 DS:Group 2 LB:Library 2 SM:Sample +@RG ID:grp3 DS:Group 3 LB:Library 3 SM:Sample +@PG ID:prog1 PN:emacs CL:emacs VN:23.1.1 +@CO The MIT License +@CO +@CO Copyright (c) 2014, 2015 Genome Research Ltd. +@CO +@CO Permission is hereby granted, free of charge, to any person obtaining a copy +@CO of this software and associated documentation files (the "Software"), to deal +@CO in the Software without restriction, including without limitation the rights +@CO to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +@CO copies of the Software, and to permit persons to whom the Software is +@CO furnished to do so, subject to the following conditions: +@CO +@CO The above copyright notice and this permission notice shall be included in +@CO all copies or substantial portions of the Software. +@CO +@CO THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +@CO IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +@CO FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +@CO AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +@CO LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +@CO OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +@CO THE SOFTWARE. +@SQ SN:ref1 LN:56 M5:08c04d512d4797d9ba2a156c1daba468 +@SQ SN:ref2 LN:60 M5:7c35feac7036c1cdef3bee0cc4b21437 +ref1_grp1_p001 99 ref1 1 0 10M = 25 34 CGAGCTCGGT !!!!!!!!!! MD:Z:10 NM:i:0 RG:Z:grp1 BC:Z:AC-GT H0:i:1 aa:A:! ab:A:~ fa:f:3.14159 za:Z:Hello world! ha:H:DEADBEEF ba:B:c,-128,0,127 bb:B:C,0,127,255 bc:B:s,-32768,0,32767 bd:B:S,0,32768,65535 be:B:i,-2147483648,0,2147483647 bf:B:I,0,2147483648,4294967295 bg:B:f,2.71828,6.626e-34,2.9979e+09 +ref1_grp1_p001 147 ref1 25 12 10M = 1 -34 TAGAGTCGAC ---------- MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp1_p002 99 ref1 5 2 10M = 29 34 CTCGGTACCC ########## MD:Z:10 NM:i:0 RG:Z:grp1 BC:Z:AATT+CCGG H0:i:1 aa:A:a ab:A:z fa:f:4.3597e-18 za:Z:Another string ha:H:2000AD +ref1_grp1_p002 147 ref1 29 14 10M = 5 -34 GTCGACCTGC ////////// MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp1_p003 99 ref1 9 4 10M = 33 34 GTACCCGGGG %%%%%%%%%% MD:Z:10 NM:i:0 fa:f:1.66e-27 RG:Z:grp1 +ref1_grp1_p003 147 ref1 33 16 10M = 9 -34 ACCTGCAGGC 1111111111 MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp1_p004 99 ref1 13 6 10M = 37 34 CCGGGGATCC '''''''''' MD:Z:10 NM:i:0 fa:f:1.38e-23 za:Z:xRG:Z:grp2 RG:Z:grp1 +ref1_grp1_p004 147 ref1 37 18 10M = 13 -34 GCAGGCATGC 3333333333 MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp1_p005 99 ref1 17 8 10M = 41 34 GGATCCTCTA )))))))))) MD:Z:10 NM:i:0 RG:Z:grp1 ia:i:40000 +ref1_grp1_p005 147 ref1 41 20 10M = 17 -34 GCATGCAAGC 5555555555 MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp1_p006 99 ref1 21 10 10M = 45 34 CCTCTAGAGT ++++++++++ MD:Z:10 NM:i:0 RG:Z:grp1 ia:i:255 +ref1_grp1_p006 147 ref1 45 22 10M = 21 -34 GCAAGCTTGA 7777777777 MD:Z:10 NM:i:0 RG:Z:grp1 +ref1_grp2_p001 99 ref1 3 1 8M2S = 27 34 AGCTCGGTAC """""""""" MD:Z:8 NM:i:0 RG:Z:grp2 BC:Z:TG+CA QT:Z:ab+cd H0:i:1 aa:A:A ab:A:Z fa:f:6.67e-11 za:Z:!"$%^&*() ha:H:CAFE +ref1_grp2_p001 2147 ref1 11 1 8H2M = 27 34 AC "" MD:Z:2 NM:i:0 RG:Z:grp2 BC:Z:TG-CA H0:i:1 aa:A:A ab:A:Z fa:f:6.67e-11 za:Z:!"$%^&*() ha:H:CAFE +ref1_grp2_p001 147 ref1 27 13 10M = 3 -34 GAGTCGACCT .......... MD:Z:10 NM:i:0 RG:Z:grp2 +ref1_grp2_p002a 99 ref1 7 3 10M = 31 34 CGGTACCCGG $$$$$$$$$$ MD:Z:10 NM:i:0 fa:f:6.022e+23 RG:Z:grp2 +ref1_grp2_p002 99 ref1 7 3 10M = 31 34 CGGTACCCGG $$$$$$$$$$ MD:Z:10 NM:i:0 fa:f:6.022e+23 RG:Z:grp2 +ref1_grp2_p002 147 ref1 31 15 10M = 7 -34 CGACCTGCAG 0000000000 MD:Z:10 NM:i:0 RG:Z:grp2 +ref1_grp2_p003 99 ref1 11 5 10M = 35 34 ACCCGGGGAT &&&&&&&&&& MD:Z:10 NM:i:0 RG:Z:grp2 ia:i:4294967295 +ref1_grp2_p003 147 ref1 35 17 10M = 11 -34 CTGCAGGCAT 2222222222 MD:Z:10 NM:i:0 RG:Z:grp2 +ref1_grp2_p004 99 ref1 15 7 10M = 39 34 GGGGATCCTC (((((((((( MD:Z:10 NM:i:0 RG:Z:grp2 ia:i:-2147483648 +ref1_grp2_p004 147 ref1 39 19 10M = 15 -34 AGGCATGCAA 4444444444 MD:Z:10 NM:i:0 RG:Z:grp2 +ref1_grp2_p005 99 ref1 19 9 10M = 43 34 ATCCTCTAGA ********** MD:Z:10 NM:i:0 RG:Z:grp2 ia:i:-1000 +ref1_grp2_p005 147 ref1 43 21 10M = 19 -34 ATGCAAGCTT 6666666666 MD:Z:10 NM:i:0 RG:Z:grp2 +ref1_grp2_p006 99 ref1 23 11 10M = 47 34 TCTAGAGTCG ,,,,,,,,,, MD:Z:10 NM:i:0 RG:Z:grp2 ia:i:-1 +ref1_grp2_p006 147 ref1 47 23 10M = 23 -34 AAGCTTGAGT 8888888888 MD:Z:10 NM:i:0 RG:Z:grp2 +ref2_grp3_p001 83 ref2 1 99 15M = 31 45 ATTCTATAGTGTCAC ~~~~~~~~~~~~~~~ MD:Z:15 NM:i:0 RG:Z:grp3 +ref2_grp3_p001 163 ref2 31 99 15M = 1 -45 CTGTTTCCTGTGTGA ||||||||||||||| MD:Z:0T0A0A1C0A0T0G0G0T0C0A1A0G0 NM:i:13 RG:Z:grp3 +ref2_grp3_p002 99 ref2 46 99 15M = 16 -45 CTGTTTCCTGTGTGA {{{{{{{{{{{{{{{ MD:Z:15 NM:i:0 RG:Z:grp3 +ref2_grp3_p002 147 ref2 16 99 15M = 46 45 CTAAATAGCTTGGCG }}}}}}}}}}}}}}} MD:Z:15 NM:i:0 RG:Z:grp3 +ref2_grp3_p003 99 ref2 1 99 15M = 1 15 ACGTMRWSYKVHDBN 0123456789abcd! MD:Z:1T0T0C0T0A0T0A0G0T0G0T0C0A0C0 NM:i:14 RG:Z:grp3 +ref2_grp3_p003 147 ref2 1 99 15M = 1 -15 NVHDBMRSWYKACGT !dcba9876543210 MD:Z:0A0T0T0C0T0A0T0A0G0T0G0T1A0C0 NM:i:14 RG:Z:grp3 +ref12_grp1_p001 97 ref1 36 50 10M ref2 2 0 TGCAGGCATG AAAAAAAAAA MD:Z:10 NM:i:0 RG:Z:grp1 +ref12_grp1_p001 145 ref2 2 50 10M ref1 36 0 TTCTATAGTG BBBBBBBBBB MD:Z:10 NM:i:0 RG:Z:grp1 +ref12_grp2_p001 97 ref1 46 50 10M ref2 12 0 CAAGCTTGAG AAAAAAAAAA MD:Z:10 NM:i:0 RG:Z:grp2 +ref12_grp2_p001 145 ref2 12 50 10M ref1 46 0 TCACCTAAAT BBBBBBBBBB MD:Z:10 NM:i:0 RG:Z:grp2 +unaligned_grp3_p001 77 * 0 0 * * 0 0 CACTCGTTCATGACG 0123456789abcde RG:Z:grp3 diff --git a/test/test.pl b/test/test.pl index 0c31877c1..f95a0bcc5 100755 --- a/test/test.pl +++ b/test/test.pl @@ -61,7 +61,6 @@ test_addrprg($opts); test_addrprg($opts, threads=>2); - print "\nNumber of tests:\n"; printf " total .. %d\n", $$opts{nok}+$$opts{nfailed}+$$opts{nxfail}+$$opts{nxpass}; printf " passed .. %d\n", $$opts{nok}; @@ -2240,6 +2239,10 @@ sub test_bam2fq test_cmd($opts, out=>'bam2fq/2.stdout.expected', out_map=>{'1.fq' => 'bam2fq/3.1.fq.expected', '2.fq' => 'bam2fq/3.2.fq.expected', 's.fq' => 'bam2fq/3.s.fq.expected'}, cmd=>"$$opts{bin}/samtools fastq @$threads -s $$opts{path}/s.fq -1 $$opts{path}/1.fq -2 $$opts{path}/2.fq $$opts{path}/dat/bam2fq.002.sam"); # basic 2 output test with singleton tracking with a singleton as last read test_cmd($opts, out=>'bam2fq/2.stdout.expected', out_map=>{'1.fq' => 'bam2fq/4.1.fq.expected', '2.fq' => 'bam2fq/4.2.fq.expected', 's.fq' => 'bam2fq/4.s.fq.expected'}, cmd=>"$$opts{bin}/samtools fastq @$threads -s $$opts{path}/s.fq -1 $$opts{path}/1.fq -2 $$opts{path}/2.fq $$opts{path}/dat/bam2fq.003.sam"); + # tag output test with singleton tracking with a singleton as last read + test_cmd($opts, out=>'bam2fq/2.stdout.expected', out_map=>{'1.fq' => 'bam2fq/4.1.fq.expected', '2.fq' => 'bam2fq/4.2.fq.expected', 's.fq' => 'bam2fq/4.s.fq.expected', 'bc.fq' => 'bam2fq/bc.fq.expected'}, cmd=>"$$opts{bin}/samtools fastq @$threads --barcode-tag BC --index-format 'n2i2' --i1 $$opts{path}/bc.fq -s $$opts{path}/s.fq -1 $$opts{path}/1.fq -2 $$opts{path}/2.fq $$opts{path}/dat/bam2fq.004.sam"); + # tag output test with separators and -N flag + test_cmd($opts, out=>'bam2fq/2.stdout.expected', out_map=>{'1.fq' => 'bam2fq/5.1.fq.expected', '2.fq' => 'bam2fq/5.2.fq.expected', 's.fq' => 'bam2fq/5.s.fq.expected', 'bc_split.fq' => 'bam2fq/bc_split.fq.expected'}, cmd=>"$$opts{bin}/samtools fastq @$threads --barcode-tag BC -N --index-format 'n*i*' --i1 $$opts{path}/bc_split.fq -s $$opts{path}/s.fq -1 $$opts{path}/1.fq -2 $$opts{path}/2.fq $$opts{path}/dat/bam2fq.005.sam"); } sub test_depad From d7b0234d6cd27b790b8bb254d950ba35e277c0fb Mon Sep 17 00:00:00 2001 From: James Bonfield Date: Thu, 20 Apr 2017 16:43:36 +0100 Subject: [PATCH 18/20] Output CRAM as temporary sort format iff there are more than 64k cigar ops. Fixes #667. In theory we could use CRAM for all temporary output files, but currently it is too high in memory usage. See the comments included in this commit for a discussion on how this could be avoided, but for now we are making the smallest change that fixes the issue as there are larger issues yet to resolve in sort memory usage. --- bam_sort.c | 46 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 34 insertions(+), 12 deletions(-) diff --git a/bam_sort.c b/bam_sort.c index 5dab20e5e..be9789c3a 100644 --- a/bam_sort.c +++ b/bam_sort.c @@ -1219,6 +1219,16 @@ int bam_merge_core2(int by_qname, const char *out, const char *mode, if ((translation_tbl+i)->lost_coord_sort && !by_qname) { fprintf(stderr, "[bam_merge_core] Order of targets in file %s caused coordinate sort to be lost\n", fn[i]); } + + // Potential future improvement is to share headers between CRAM files for + // samtools sort (where all headers are identical. + // Eg: + // + // if (i > 1) { + // sam_hdr_free(cram_fd_get_header(fp[i]->fp.cram)); + // cram_fd_set_header(fp[i]->fp.cram, cram_fd_get_header(fp[0]->fp.cram)); + // sam_hdr_incr_ref(cram_fd_get_header(fp[0]->fp.cram)); + // } } // Did we get an @HD line? @@ -1648,18 +1658,30 @@ static void *worker(void *data) name = (char*)calloc(strlen(w->prefix) + 20, 1); if (!name) { w->error = errno; return 0; } sprintf(name, "%s.%.4d.bam", w->prefix, w->index); - if (write_buffer(name, "wbx1", w->buf_len, w->buf, w->h, 0, NULL) < 0) - w->error = errno; - -// Consider using CRAM temporary files if the final output is CRAM. -// Typically it is comparable speed while being smaller. -// hts_opt opt[2] = { -// {"version=3.0", CRAM_OPT_VERSION, {"3.0"}, NULL}, -// {"no_ref", CRAM_OPT_NO_REF, {1}, NULL} -// }; -// opt[0].next = &opt[1]; -// if (write_buffer(name, "wc1", w->buf_len, w->buf, w->h, 0, opt) < 0) -// w->error = errno; + + uint32_t max_ncigar = 0; + int i; + for (i = 0; i < w->buf_len; i++) { + uint32_t nc = w->buf[i]->core.n_cigar; + if (max_ncigar < nc) + max_ncigar = nc; + } + + if (max_ncigar > 65535) { + htsFormat fmt; + memset(&fmt, 0, sizeof(fmt)); + if (hts_parse_format(&fmt, "cram,version=3.0,no_ref,seqs_per_slice=1000") < 0) { + w->error = errno; + free(name); + return 0; + } + + if (write_buffer(name, "wcx1", w->buf_len, w->buf, w->h, 0, &fmt) < 0) + w->error = errno; + } else { + if (write_buffer(name, "wbx1", w->buf_len, w->buf, w->h, 0, NULL) < 0) + w->error = errno; + } free(name); return 0; From 0aa9fab90bddca2ce277c831ed6407f53cc42496 Mon Sep 17 00:00:00 2001 From: James Bonfield Date: Mon, 24 Apr 2017 09:26:58 +0100 Subject: [PATCH 19/20] Bam2fq comment improvements (no code changed). --- sam_view.c | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/sam_view.c b/sam_view.c index 9f3e998c5..9c2d15b08 100644 --- a/sam_view.c +++ b/sam_view.c @@ -724,14 +724,17 @@ typedef struct bam2fq_state { } bam2fq_state_t; /* - * get and decode the read from a BAM record - * Why on earth isn't this (and the next few functions) part of htslib? - * The whole point of an API is that we should not have to know the internal structure... + * Get and decode the read from a BAM record. + * + * TODO: htslib really needs an interface for this. Consider this or perhaps + * bam_get_seq_str (current vs original orientation) and bam_get_qual_str + * functions as string formatted equivalents to bam_get_{seq,qual}? */ /* - * reverse a string in place - (stolen from http://stackoverflow.com/questions/8534274/is-the-strrev-function-not-available-in-linux) + * Reverse a string in place. + * From http://stackoverflow.com/questions/8534274/is-the-strrev-function-not-available-in-linux. + * Author Sumit-naik: http://stackoverflow.com/users/4590926/sumit-naik */ static char *reverse(char *str) { From f13cf4afff4b04fd015bc0512099448d4fd54089 Mon Sep 17 00:00:00 2001 From: James Bonfield Date: Fri, 5 May 2017 12:24:34 +0100 Subject: [PATCH 20/20] NEWS update. (Ideally merge on day on release, editing date if required.) --- NEWS | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/NEWS b/NEWS index fc3010772..58dcd0f88 100644 --- a/NEWS +++ b/NEWS @@ -1,9 +1,28 @@ -Release x.x (the future) -~~~~~~~~~~~~~~~~~~~~~~~~~ +Release 1.4.1 (8th May 2017) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * Added options to fastq to create fastq files from BC (or other) tags. +* Samtools view has gained a -G option to exclude on all bits + set. For example to discard reads where neither end has been + mapped use "-G 12". + +* Samtools cat has a -b option to ease concatenation of many + files. + +* Added misc/samtools_tab_completion for bash auto-completion of + samtools sub-commands. (#560) + +* Samtools tview now has J and K keys for verticale movement by 20 + lines. (#257) + +* Various compilation / portability improvements. + +* Fixed issue with more than 65536 CIGAR operations and SAM/CRAM files. + (#667) + + Release 1.4 (13 March 2017) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~