samtools-help — bug report, feature request, troubleshooting and anything else

[Samtools-help] SAMtools and HTSlib release 1.22.1

[Robert D. <rm...@sa...>]

SAMtools and HTSlib version 1.22.1 are now available from GitHub and
SourceForge.  This fixes crashes that could happen when reading malformed
CRAM files, or when writing files with very long secondary alignments as
CRAM 3.1.  It also fixes some issues that caused builds to fail on
some older platforms, or where libcurl was unavailable.

https://github.com/samtools/htslib/releases/tag/1.22.1
https://github.com/samtools/samtools/releases/tag/1.22.1
https://sourceforge.net/projects/samtools/

The main changes are listed below:

-------------------------------------------------------------------------
htslib - changes v1.22.1
-------------------------------------------------------------------------

Bug Fixes
---------

* SECURITY fix: Prevent CRAM byte_array decoder from overflowing its output
   buffer.  This could be triggered by certain malformed CRAM inputs.
   (PR #1934)

* Two fixes for crashes reported when trying to save data with very long
   alignment records with sequence '*' as CRAM 3.1:

   - The htscodecs submodule is updated to v1.6.4. This includes a fix to the
     rANS encoder to prevent it from failing on these inputs.
     (PR #1935.  Reported by Martin Pollard)

   - Improved error handling in cram_compress_block2().  If the
     previously-chosen CRAM compression method starts to fail, it will now
     try other methods instead of giving up immediately.
     (PR #1931.  Reported by Martin Pollard)

* Fix warnings due to the wrong datatype being passed to curl_easy_setopt()
   (PR #1925.  Thanks to John Marshall)

* Prevent instances of `memcpy(out, NULL, 0)`, which is strictly undefined
   behaviour.
   (PR #1930.  Thanks to Ben Lawrence).

Build Changes
-------------

* Fixed compilation against older glibc / macOS SDKs that incorrectly
   suppressed some symbols if _XOPEN_SOURCE was defined.
   (PR #1928.  Reported by John Marshall)

* Fixed ref-cache configure check for libcurl, so that if libcurl is not
   available, or turned off by `./configure --disable-libcurl`, the
   `ref-cache` build will be automatically disabled as well.
   (PR #1929, fixes #1926.  Reported by biounix)

-------------------------------------------------------------------------
samtools - changes v1.22.1
-------------------------------------------------------------------------

Bug fixes:

* SECURITY: Fix use-after-free bug in `samtools mpileup -a` due to reference
   sequences being discarded too early.  This could happen when the `-a`
   option was used, and all the alignments for one of the references started
   at the same position.  This caused mpileup to try to load the next
   reference before it had finished writing data for the previous one out.
   (PR #2229, fixes #2227.  Reported by Pouya Kheradpour)

* This release bundles htslib-1.22.1 and htscodecs v1.6.4.  Fixes a possible
   buffer overflow on some invalid CRAM inputs; and a failure when trying to
   save data with very long alignment records with sequence '*' as CRAM
   version 3.1.

   See the HTSlib and htscodecs NEWS files for details.

Documentation:

* Clarify that `-@` starts additional threads in `samtools index` help, and
   fix formatting issues in the `samtools reset` manual page.
   (PR #2225)

Build changes:

* Added settings to limit memory used by `samtools sort` when running tests.
   (PR #2226)


----------------------------------------------------------------------
The Wellcome Sanger Institute is operated by Genome Research Limited, a charity registered in England with number 1021457 and a company registered in England with number 2742969, whose registered office is Wellcome Sanger Institute, Wellcome Genome Campus, Hinxton, CB10 1SA.

[Samtools-help] Release 1.22

[Robert D. <rm...@sa...>]

Samtools (and HTSlib and BCFtools) version 1.22 is now available from
GitHub and SourceForge.

https://github.com/samtools/htslib/releases/tag/1.22
https://github.com/samtools/samtools/releases/tag/1.22
https://github.com/samtools/bcftools/releases/tag/1.22
https://sourceforge.net/projects/samtools/

The main changes are listed below:

-------------------------------------------------------------------------
htslib - changes v1.22
-------------------------------------------------------------------------

Note this release changes the default output CRAM version from 3.0 to
3.1.HTSlib and SAMtools have been able to read CRAM 3.1 since version
1.12, however other tools may not yet be able to cope.  We know Noodles
reads CRAM3.1 and htsjdk has a draft implementation that has not yet been
released.

HTSlib has options for modifying the output formats, which are exposed
in SAMtools.  When specifying an output format you can explicitly change
the version via e.g. `samtools view -O cram,version=3.0 ...`.

Further documentation on this change can be found at
https://www.htslib.org/benchmarks/CRAM.html

HTSlib no longer fetches CRAM reference data from EBI's server by default.
Your organisation may wish to set up local infrastructure to supply
reference sequences, e.g., using the new ref-cache tool included in this
HTSlib release. See the REF_CACHE and REF_PATH environment variables
documented in https://www.htslib.org/doc/reference_seqs.html and the
SAMtools manpage for details.

Updates
-------

* NEW. Add ref-cache, a caching proxy for reference sequences.  This is a
   local server of reference sequences, for use when encoding or decoding CRAM
   files that use reference-based compression. (PR #1911, PR #1921, PR #1922)

* Add support for matching VCF lines by ID. (PR #1844, addresses issue
   bcftools#1739 reported by Han Cao)

* Make it possible to test for VCF_REF as declared in the documentation.
   (PR #1879)

* Updated VCF code to work with VCF 4.4 prefixed phasing info. (PR#1861,
   fixes #1847.  Reported by John Marshall)

* Use the highest VCF version when merging headers. (PR#1912, see
   bcftools#2395 and bcftools#2404)

* Update RLEN calculation for VCF 4.4 and 4.5. (PR#1897, fixes #1820.
   Reported by Dave Lawrence)

* Convert U to T instead of U to N when sam_parsing.  Though SAM format
   itself can contain U the BAM format cannot. (PR #1854, fixes samtools#2131
   reported by James Ferguson)

* Add an hts_crc32 function to use zlib or libdeflate.  The libdeflate crc32
   function is faster than native zlib and should be used when available. (PR
   #1850)

* Increase the input block size for bgzip. This deals with a slow down
   introduced in PR #1493 when reading from a pipe. (PR #1768, fixes
   #1767.  Reported by Konstantin Riege)

* Allow BYTE_ARRAY_STOP to work on non-zero STOP code with TOK3.  Although
   the htscodecs name tokeniser uses a NUL between names there is no reason
   why another value could not be used.  This change lets CRAM recognise other
   separator values. (PR #1871)

* Remove cram seek ability to do range queries via SEEK_CUR.  A probable
   misfeature from the original implementation. (PR #1878, fixes #1877.
   Reported by Rick Wertenbroek)

* Add hts_tpool_worker_id() API.  This may be used to associate data with a
   thread rather than to a job. (PR #1875)

* Update bcf_synced_reader to use htsFile. (PR #1868, implements #1862.
   Requested by Brent Pedersen)

* Exit with return value 1 on tabix parse error.  This previously returned 0.
   (PR #1887, fixes #1885.  Reported by Fan-iX)

* Automatically recognise BED vs TSV files and add the option -C, --coords to
   set index positions (1 or 0 based coordinates) in annot-tsv. (PR #1894)

* Reading SQ lines with multiple differing LN will now fail.  Such lines are
   invalid (by the spec) and previous handling was inconsistent. (PR #1882,
   fixes #1866)

* Return errors instead of EOF after all I/O errors etc in
   hts_itr_multi_next/sam_itr_next/sam_read1/vcf_parse/bcf_read. (PR#1899.
    Thanks to John Marshall)

* Remove UR:file:// and UR:ftp:// from ref search path, plus REF_PATH to EBI.
   Removing EBI as the default fallback when REF_PATH not set prevents the
   unintended DDOS on EBI's servers. (PR#1881. PR#1915, fixes oss-fuzz
   issue 418125747)

Build Changes
-------------

* Detect the presence of getauxval() and elf_aux_info() for *BSD variants.
   (PR #1835, thanks to Brad Smith)

* Make HAVE_ATTRIBUTE_TARGET check also check that SSSE3 intrinsics work.
   Mainly for use with old compilers. (PR #1886, fixes #1838 and
   pysam-developers/pysam#1327.  Thanks to John Marshall)

* Fix broken tests due to MSYS2 changes. Due to changes in how MSYS2 perl
   reported the identity of the OS it was built for, our tests were failing to
   adapt to the Windows style file locations. (PR #1892)

* Updated htscodecs submodule to version 1.6.3 (PR #1917)

* Fix the script used to build the symbol version file. (PR #1918)

Bug fixes
---------

* Fix possible 1 byte underflow in find_file_extension(). Fixes an issue
   reported by OSS-Fuzz. (PR #1840, fixes oss-fuzz id 71740)

* Replace home-brew string end searching with memchr() to speed up looking at
   long aux tags. (PR #1842)

* Prevent segfault on empty tbi index.  This could happen when a VCF file has
   a header but no data lines. (PR #1845, fixes bcftools#2286.  Reported by
   Devon Ryan)

* Fix CRAM embed_ref=2 with seqs overlapping ref end. (PR #1848 and PR #1849
   which fixed oss-fuzz issue 372547397)

* Fix sam_hdr_remove_line_pos() not dealing with the 0 index position
   properly. (PR #1853.  Thanks to Julian Regalado Perez)

* Fix threaded sam_read1() after EOF.  Prevents sam_read1() getting stuck
   when trying to read after EOF and waiting forever for data that is never
   going to arrive. (PR #1856, fixes #1855.  Reported by Yan Gao)

* Fix a bug in breakend detection. It was incorrectly assuming that the ALT
   allele is of equal length to REF allele, but the VCF specification allows
   breakend insertions. (PR #1858, fixes bcftools#2317.  Reported by
   Nicolai von Kügelgen).

* Fix cram_encode fuzzer issue caused by negative reference lengths.
   Reported by OSS-Fuzz. (PR #1863 fixes oss-fuzz issue 382922241)

* Fixed a typo in vcf.h. (PR #1870, thanks to Yu Wang)

* Reset variant types after updating alleles with bcf_update_alleles() or
   bcf_update_alleles_str().  Prevents an out-of-bounds access by bcftools
   consensus. (PR #1883)

* Recognize T > A[chr15:12345[ breakend type in VCF. (PR#1903, fixes
   bcftools#2389.  Reported by Dennis Hendriksen)

* Fix possible buffer overruns in expand_path(). (PR#1907)

Documentation updates
---------------------

* Add instructions to INSTALL for FreeBSD, NetBSD and OpenBSD. (PR #1843)

* Clarify bam_set1() parameter documentation to note that quality values do
   not have the ASCII 33 offset. (PR #1891.  Thanks to Chris Wright)

* Fixed incorrectly named table in bam1_t structure documentation. (PR #1923.
   Thanks to Julian Hess)

-------------------------------------------------------------------------
samtools - changes v1.22
-------------------------------------------------------------------------

Note this release changes the default output CRAM version from 3.0 to
3.1.HTSlib and SAMtools have been able to read CRAM 3.1 since version
1.12, however other tools may not yet be able to cope.  We know Noodles
reads CRAM3.1 and htsjdk has a draft implementation, but not yet released.

HTSlib has options for modifying the output formats, which are exposed
in SAMtools.  When specifying an output format you can explicitly change
the version via e.g. `samtools view -O cram,version=3.0 ...`.

Further documentation on this change can be found at
https://www.htslib.org/benchmarks/CRAM.html

HTSlib no longer fetches CRAM reference data from EBI's server by default.
Your organisation may wish to set up local infrastructure to supply
reference sequences, e.g., using the new ref-cache tool included in this
HTSlib release. See the REF_CACHE and REF_PATH environment variables
documented in https://www.htslib.org/doc/reference_seqs.html and the
SAMtools manpage for details.

New work and changes:

* New `samtools checksum` command.  This checksums sequence, name, quality
   and barcode tags in an order and orientation agnostic manner, to facilitate
   validation of no data loss between raw fastq (or unmapped crams) through
   alignment, duplication marking, sorting, and other processing operations to
   get to the final aligned bam/cram. (PR#2122)

* Extend `samtools sort -M` to distinguish between mapped and unmapped files.
   (PR#2110, fixes #2105.  Reported by Armin Töpfer)

* Allow the `samtools sort` "merging from..." message to be silenced. Setting
   the verbosity to 0 or 1 will now silence this message. (PR#2197, resolves
   #2185.  Requested by Alex Predeus)

* Add `--save-counts` option to `samtools view`.  Adds an option to store
   counts of records processed, accepted and rejected by filtering to a file.
   (PR#2120, resolves #2038.  Requested by Chang Y)

* `samtools fasta` and `fastq` can now make faidx/fqidx indexes while writing
   using the `--write-index` option. (PR#2125, resolves #2118.  Requested by
   Filipe G. Vieira)

* Add a warning for `samtools fastq` on coordinate sorted data. (PR#2176,
   fixes #2169 and #2161.  Reported by wook2014)

* `samtools tview` add `-i` to hide inserts. (PR#2123.  Thanks to
   Benjamin Bræstrup Sayoc)

* Show optional headers with `samtools bedcov -H`. (PR#2140, fixes #2126.
   Reported by biounix)

* `samtools consensus` now supports proper multi-threading.  Previously
   this was restricted to decompression only, but it should now scale
   better. (PR#2174, supersedes PR#2141)

* Add `samtools consensus -T ref.fa` functionality.  This reports the
   reference value if a consensus value cannot be calculated. (PR#2153,
   fixes an additional request in #1915)

* In `samtools consensus`, do not use consensus N for "*" (absent) calls that
   are masked due to insufficient depth. (PR#2204, fixes #2167.  Reported by
   sanschiffre)

* Improve `plot-bamstats` quality plots. (PR#2143 combined with PR#2116
   (thanks to James Gilbert))

* Make `reheader -h` use /tmp and honour TMPDIR. (PR#2168, related to #2165.
   Reported by Zhang Yuanfeng)

* Set sort order header tag to unsorted when ordering is lost during
   `samtools merge`. (PR#2173, fixes #2159.  Reported by Filipe G. Vieira)

* Protect against merging CRAM files with different headers. (PR#2220, fixes
   #2218.  Reported by Kevin Lewis)

* `samtools stats` bug-fix to checksum calculation for quality values.  This
   corrects the checksums but in turn makes the calculated value different to
   that reported by previous samtools versions. (PR#2193, fixes #2187)

* Clarification for `samtools stats` when used on files with different sort
   orders. (PR#2198, fixes #2177.  Reported by Filipe G. Vieira)

* In `samtools stats`, dovetailed (completely overlapping) read pairs are
   now always counted as inward-oriented.  Previously they could have been
   inwards or outwards depending on read ordering. (PR#2216, resolves #2210.
    Requested by Pontus Hüer)

Documentation:

* Correct the example for 1:1 `samtools consensus` coords. (PR#2113, fixes
   #2111.  Reported by schorlton-bugseq)

* Documents the fastq format options used in SAMtools and HTSlib. (PR#2123,
   fixes #2121)

* Remove mention of threads from `samtools cat` man page. (PR#2162, fixes
   #2160.  Reported by Brandon Pickett)

* Update `samtools merge` man page to include `--template-coordinate`.
   (PR#2164.  Thanks to Nils Homer)

* Revised CRAM reference sequence documentation in the samtools man page.
   (PR#2178)

* Added fish shell completion and renamed completion for bash shell.  These
   files can be copied to appropriate directories by the user.  For full
   functionality it requires Python3.5+ and installed samtools manpages.
   (PR#2203.  Thanks to LunarEclipse363)

* Fix URL printed by the `seq_cache_populate.pl` script. (PR#2222.  Thanks to
   Charles Plessy)

Bug fixes:

* `samtools consensus` previously could give different results for BAM and
   CRAM files with the same content.  This was because MD/NM tag generation
   was disabled in CRAM, but the `decode_md=0` option did nothing with BAM.
   Note with `--no-adj-MQ` both BAM and CRAM gave identical results. Now
   use `--input-fmt-option decode_md=0` to get the old CRAM behaviour.
   Otherwise, both BAM and CRAM will be utilising MD/NM to locally modify
   mapping quality. (PR#2156)

* `samtools consensus` without `-a` previously still padded with leading Ns
   in some cases.  It now consistently removes both leading and trailing Ns.
   Use "-a" if you want all reference bases displayed. (Part of PR#2174 above)

* Change how `markdup` looks for single reads.  Due to changes to `fixmate`
   in 1.21 `markdup` no longer recognised single reads that would have
   normally have been part of a pair. (PR#2117, fixes #2117.  Reported by
   Kristy Horan)

* Fix `samtools merge` crash on BAM files with malformed headers. (PR#2128,
   fixes #2127.  Reported by Frostb1te)

* Fix `faidx --write-index` invalid free. (PR#2147, fixes 2142.  Reported by
   Alex Leonard)

* Fix `samtools fastq -i` to force CRAM aux tag decoding. (PR#2155, fixes
   #2155.  Reported by Alex Leonard)

Non user-visible changes and build improvements:

* Improve htslib#PRnum support for Cirrus-CI and GitHub Actions. (PR#2115)

* Fix broken tests due to MSYS2 changes. Due to changes in how MSYS2 perl
   reported the identity of the OS it was built for, our tests were failing to
   adapt to the Windows style file locations. (PR #2196)

* Upgrade to `_XOPEN_SOURCE=700`, to match HTSlib.  Also replace `usleep()`
   with `nanosleep()`. (PR#2221)

-------------------------------------------------------------------------
bcftools - changes v1.22
-------------------------------------------------------------------------

Changes affecting the whole of bcftools, or multiple commands:

* Add support for matching lines by ID via the --pair-logic and --collapse
   options (#1739)

* The -i/-e filtering expressions

     - The expressions now properly match the regex negation of missing
       values, e.g. -i 'TAG!~"\."' (#2355)

     - Added support for Fisher's exact test

* Add the option `-v, --verbosity INT` to all bcftools commands and
   plugins. Verbosity values bigger than 3 are passed to the underlying
   HTSlib library so that the user can investigate network issues and
   other problems occurring at the library level.

Changes affecting specific commands:

* bcftools annotate

     - Fix Number in the header definition of transferred FILTER and ID tags
       (#2335)

* bcftools call

     - The `-s, --samples` option was not working properly, now also
       supporting sample negation as advertised in the manual page, e.g.
       `-s ^sample1,sample2` to include all samples but sample1 and
       sample2 (#2380)

* bcftools consensus

     - Preserve entire missing gVCF blocks with --missing (#2350)

     - Fixed a bug, the `-S, --samples-file` option is no longer ignored
       (#2398)

* bcftools convert

     - The command `convert --gvcf2vcf` was not filling the REF allele when
       BCF was output (#243)

* bcftools csq

     - Check the input GFF for features outside transcript boundaries and
       extend the transcript to contain the feature fully (#2323)

     - Add experimental support for alternative genetic code tables,
       accessible via a new option `-C, --genetic-code` (#2368)

     - Change in the `--unify-chr-names` option, no automatic sequence
       name modification is attempted anymore, the prefixes to trim must
       be given explictly. For example, if run with
       `--unify-chr-names chr,Chromosome,-`, the program will trim the
       "chr" prefix in the VCF, "Chromosome" in the GFF, leaving the
       fasta unchanged (#2378)

* bcftools +fill-tags

     - Thanks to the extension of filtering expressions with Fisher's exact
       test, the plugin can now be used to add FT annotation (#1582)

* bcftools merge

     - Preserve phasing in half-missing genotypes (#2331)

     - The option `--merge none` is expected to create no new multiallelic
       sites, but it should allow to merge, say, A>C with A>C,AT (#2333)

     - Make `--merge both` work with indel-only records; for example, the
       multiallelic site G>GT,T should be merged with G>GT (#2339)

     - Do not merge symbolic alleles unless they have not just the same type,
       eg. <DEL>, but also length, i.e the INFO/END coordinate (#2362)

     - Fix a bug where an incorrectly formatted gVCF file with overlapping
       blocks would trigger an infinite loop in the program (#2410)

* bcftools mpileup

     - The -r/-R option newly merge overlapping regions, preventing the output
       of duplicate sites

* bcftools norm

     - Print the number of removed duplicate sites in the final statistics
       (#2346)

     - Preserve the original alleles in `--old-rec-tag` when `--check-ref s`
       requested (#2357)

     - Print a warning when INFO/SVLEN is not defined as Number=A (#2371)

* plot-vcfstats

     - Make the option `-s, --sample-names` functional again (#2353)

* bcftools +prune

     - New option to remove or annotate clusters of sites within a window

* bcftools query

     - The functions used in -i/-e filtering expressions (such as SUM,
       MEDIAN, etc) can be now used in formatting expressions (#2271). If
       the VCF contains INFO/AD and FORMAT/AD, try:

  bcftools query test.vcf -f '%CHROM:%POS \t [ %AD]  \t [ %sSUM(FMT/AD)]'
  bcftools query test.vcf -f '%CHROM:%POS \t [ %AD]  \t [ %SUM(FMT/AD)]'
  bcftools query test.vcf -f '%CHROM:%POS \t [ %AD]  \t   %SUM(FMT/AD)'
  bcftools query test.vcf -f '%CHROM:%POS \t [ %AD]  \t   %SUM(INFO/AD)'

     - Make it possible to refer to the ID column from the FORMAT expression
       (#2337)

  bcftools query test.vcf -f 'ID=%ID  ID=[ %/ID]  vs  FMT_ID=[ %ID]'

* bcftools roh

     - New visualization tool misc/roh-viz, see below

* bcftools +setGT

     - Support for setting missing genotypes with arbitrary ploidy
       via `-n c:./.` (#2303)

* bcftools +split-vep

     - The `-s, --select` option was extended to print only one consequence.
       Previously it was possible to select a single transcript (e.g., the
       one with the worst consequence), and it was possible to filter by
       consequence severity (e.g., missing or worse), but in some cases
       multiple consequences are reported within a single transcript (e.g.,
       start_lost&splice_region). The extended option allows to print the
       worst part, for example as
          --select primary:missense+:worst

* bcftools +trio-dnm2

     - Fix a problem with --strictly-novel option which would neglect the
       presence of the apparent de novo allele in the father for male
       offspring

     - Fix a problem with uncalled mosaic chrX variants in males

* roh-viz

     - HTML/JavaScript visualization of bcftools/roh output and homozygosity
       rate.

* bcftools +vrfs

     - New experimental plugin for scoring variants and assess site noisiness
       (variant read frequency profiles) from a large number of unaffected
       parental samples

----------------------------------------------------------------------
The Wellcome Sanger Institute is operated by Genome Research Limited, a charity registered in England with number 1021457 and a company registered in England with number 2742969, whose registered office is Wellcome Sanger Institute, Wellcome Genome Campus, Hinxton, CB10 1SA.

Re: [Samtools-help] mpileup differences between samtools and bcftools [EXT]

[Robert D. <rm...@sa...>]

On Tue, 22 Oct 2024, Nate Brown wrote:

> I'm running a pipeline that uses a very old version of samtools mpileup
> version 1.6.  When we switched to the more modern bcftools mpileup version
> 1.16 we noticed that the data changed slightly, for example some of the
> lower frequency calls that used to have 2 supporting reads now only have 1
> or vice versa.  Although the shifts in data are very subtle, we would like
> to set the bcftools mpileup commands so that it produces an outputs exactly
> like the old samtools mpileup.  I don't expect anyone here to know what
> these exact parameters are, but I was wondering if anyone had any general
> guidance on how the mpileup command differs between samtools and bcftools.
> Searching the readme.txt files between releases it seems that mpileup had a
> number of updates on how indels were detected and I see some of the
> defaults changed (i.e. min-iread default is 2 instead of 1, gap-frac is
> 0.05 instead of 0.002, etc. etc.)  However, even when I explicitly set the
> values for the arguments in which the defaults changed, I can never exactly
> mimic the old output.  Would it even be possible to do this?  Any advice
> would be appreciated.  Thanks!

I'm afraid it's unlikely that this is possible, due to the number of 
changes that have been made to bcftools mpileup since it was forked from 
the original in samtools.  In fact, this was one of the reasons for 
removing the VCF/BCF output from samtools mpileup - it was not generating 
some of the tags that were expected by bcftools call.

Looking at the change-log, you could try options like --config 1.12, -U 
and --full-BAQ (if using BAQ).  I wouldn't be surprised if you still see 
differences though.

Rob Davies              rm...@sa...
The Sanger Institute    http://www.sanger.ac.uk/
Hinxton, Cambs.,        Tel. +44 (1223) 834244
CB10 1SA, U.K.          Fax. +44 (1223) 494919

----------------------------------------------------------------------
The Wellcome Sanger Institute is operated by Genome Research Limited, a charity registered in England with number 1021457 and a company registered in England with number 2742969, whose registered office is Wellcome Sanger Institute, Wellcome Genome Campus, Hinxton, CB10 1SA.

[Samtools-help] mpileup differences between samtools and bcftools

[Nate B. <nat...@gm...>]

I'm running a pipeline that uses a very old version of samtools mpileup
version 1.6.  When we switched to the more modern bcftools mpileup version
1.16 we noticed that the data changed slightly, for example some of the
lower frequency calls that used to have 2 supporting reads now only have 1
or vice versa.  Although the shifts in data are very subtle, we would like
to set the bcftools mpileup commands so that it produces an outputs exactly
like the old samtools mpileup.  I don't expect anyone here to know what
these exact parameters are, but I was wondering if anyone had any general
guidance on how the mpileup command differs between samtools and bcftools.
Searching the readme.txt files between releases it seems that mpileup had a
number of updates on how indels were detected and I see some of the
defaults changed (i.e. min-iread default is 2 instead of 1, gap-frac is
0.05 instead of 0.002, etc. etc.)  However, even when I explicitly set the
values for the arguments in which the defaults changed, I can never exactly
mimic the old output.  Would it even be possible to do this?  Any advice
would be appreciated.  Thanks!

Nate

[Samtools-help] Release 1.21

[Robert D. <rm...@sa...>]

Samtools (and HTSlib and BCFtools) version 1.21 is now available from
GitHub and SourceForge.

https://github.com/samtools/htslib/releases/tag/1.21
https://github.com/samtools/samtools/releases/tag/1.21
https://github.com/samtools/bcftools/releases/tag/1.21
https://sourceforge.net/projects/samtools/  The main changes are listed below:


------------------------------------------------------------------------------
htslib - changes v1.21
------------------------------------------------------------------------------

The primary user-visible changes in this release are updates to the annot-tsv
tool and some speed improvements.  Full details of other changes and bugs
fixed are below.

Notice: this is the last SAMtools / HTSlib release where CRAM 3.0 will
be the default CRAM version.  From the next we will change to CRAM 3.1
unless the version is explicitly specified, for example using
"samtools view -O cram,version=3.0".

Updates
-------

* Extend annot-tsv with several new command line options.
     --delim permits use of other delimiters.
     --headers for selection of other header formats.
     --no-header-idx to suppress column index numbers in header.
   Also removed -h as it is now short for --headers.  Note --help still works.
   (PR #1779)

* Allow annot-tsv -a to rename annotations. (PR #1709)

* Extend annot-tsv --overlap to be able to specify the overlap fraction
   separately for source and target. (PR #1811)

* Added new APIs to facilitate low-level CRAM container manipulations,
   used by the new "samtools cat" region filtering code. Functions are:
     cram_container_get_coords()
     cram_filter_container()
     cram_index_extents()
     cram_container_num2offset()
     cram_container_offset2num()
     cram_num_containers()
     cram_num_containers_between()
   Also improved cram_index_query() to cope with HTS_IDX_NOCOOR regions.
   (PR #1771)

* Bgzip now retains file modification and access times when compressing and
   decompressing. (PR #1727, fixes #1718.  Requested by Gert Hulselmans.)

* Use FNV1a for string hashing in khash.  The old algorithm was particularly
   weak with base-64 style strings and lead to a large number of collisions.
   (PR #1806.  Fixes samtools/samtools#2066, reported by
   Hans-Joachim Ruscheweyh)

* Improve the speed of the nibble2base() function on Intel (PR #1667, PR
   #1764, PR #1786, PR #1802, thanks to Ruben Vorderman) and ARM (PR #1795,
   thanks to John Marshall).

* bgzf_getline() will now warn if it encounters UTF-16 data. (PR #1487,
   thanks to John Marshall)

* Speed up bgzf_read().  While this does not reduce CPU significantly, it
   does increase the maximum parallelism available permitting 10-15% faster
   decoding. (PR #1772, PR #1800, Issue #1798)

* Speed up faidx by use of better isgraph methods (PR #1797) and whole-line
   reading (PR #1799, thanks to John Marshall).

* Speed up kputll() function, speeding up BAM -> SAM conversion by about 5%
   and also samtools depth. (PR #1805)

* Added more example code, covering fasta/fastq indexing, tabix indexing and
   use of the thread pool. (PR #1666)

Build Changes
-------------

* Code warning fixes for pedantic compilers (PR #1777) and avoid some
   undefined behaviour (PR #1810, PR #1816, PR #1828).

* Windows based CI has been migrated from AppVeyor to GitHub Actions. (PR
   #1796, PR #1803, PR #1808)

* Miscellaneous minor build infrastructure and code fixes. (PR #1807, PR
   #1829, both thanks to John Marshall)

* Updated htscodecs submodule to version 1.6.1 (PR #1828)

* Fixed an awk script in the Makefile that only worked with gawk. (PR #1831)

Bug fixes
---------

* Fix small OSS-Fuzz reported issues with CRAM encoding and long CIGARS
   and/or illegal positions. (PR #1775, PR #1801, PR #1817)

* Fix issues with on-the-fly indexing of VCF/BCF (bcftools --write-index)
   when not using multiple threads. (PR #1837. Fixes samtools/bcftools#2267,
   reported by Giulio Genovese)

* Stricter limits on POS / MPOS / TLEN in sam_parse1().  This fixes a
   signed overflow reported by OSS-Fuzz and should help prevent other
   as-yet undetected bugs. (PR #1812)

* Check that the underlying file open worked for preload: URLs.  Fixes a NULL
   pointer dereference reported by OSS-Fuzz. (PR #1821)

* Fix an infinite loop in hts_itr_query() when given extremely large
   positions which cause integer overflow.  Also adds hts_bin_maxpos() and
   hts_idx_maxpos() functions. (PR #1774, thanks to John Marshall and
   reported by Jesus Alberto Munoz Mesa)

* Fix an out of bounds read in hts_itr_multi_next() when switching
   chromosomes.  This bug is present in releases 1.11 to 1.20. (PR #1788.
   Fixes samtools/samtools#2063, reported by acorvelo)

* Work around parsing problems with colons in CHROM names. Fixes
   samtools/bcftools#2139.  (PR #1781, John Marshall / James Bonfield)

* Correct the CPU detection for Mac OS X 10.7.  cpuid is used by htscodecs
   (see samtools/htscodecs#116), and the corresponding changes in htslib are
   PR #1785.  Reported by Ryan Carsten Schmidt.

* Make BAM zero-length intervals work the same as CRAM; permitted and
   returning overlapping records. (PR #1787.  Fixes samtools/samtools#2060,
   reported by acorvelo)

* Replace assert() with abort() in BCF synced reader.  This is not an ideal
   solution, but it gives consistent behaviour when compiling with or without
   NDEBUG. (PR #1791, thanks to Martin Pollard)

* Fixed failure to change the write block size on compressed SAM or VCF files
   due to an internal type confusion. (PR #1826)

* Fixed an out-of-bounds read in cram_codec_iter_next() (PR #1832)

------------------------------------------------------------------------------
samtools - changes v1.21
------------------------------------------------------------------------------

Notice:

* This is the last SAMtools / HTSlib release where CRAM 3.0 will be the
   default CRAM version.  From the next we will change to CRAM 3.1 unless
   the version is explicitly specified, for example using
   "samtools view -O cram,version=3.0".

New work and changes:

* `samtools reset` now removes a set of predefined auxtags, as these tags
   are no longer valid after the reset operation.  This behaviour can be
   overridden if desired. (PR #2034, fixes #2011.  Reported by Felix Lenner)

* `samtools reset` now also removes duplicate flags. (PR #2047.  Reported by
   Kevin Lewis)

* Region and section/part filtering added to CRAM `samtools cat`.  Region
   filtering permits `samtools cat` to produce new CRAMs that only cover a
   specified region. (PR #2035)

* Added a report of the number of alignments for each primer to `samtools
   ampliconclip`. (PR #2039, PR #2101, feature request #2033.  Thanks to
   Brad Langhorst)

* Make `ampliconclip` primer counts output deterministic. (PR #2081)

* `samtools fixmate` no longer removes the PAIRED flag from reads that have
   no mate.  This is done on the understanding that the PAIRED flag is a
   sequencing technology indicator not a feature of alignment.  This is a
   change to previous `fixmate` behaviour. (PR #2056, fixes #2052.  Reported
   by John Wiedenhoeft)

* Added bgzf compressed FASTA output to `samtools faidx`. (PR #2067, fixes
   #2055. Requested by Filipe G Vieira)

* Optimise `samtools depth` histogram incrementing code. (PR #2078)

* In `samtools merge` zero pad unique suffix IDs. (PR #2087, fixes #2086.
   Thanks to Chris Wright)

* `samtools idxstats` now accepts the `-X` option, making it easier to
   specify the location of the index file. (PR #2093, feature request
   #2071.  Requested by Samuel Chen)

* Improved documentation for the mpileup `--adjust-MQ` option. (PR #2098.
   Requested by Georg Langebrake)

Bug fixes:

* Avoid `tview` buffer overflow for positions with >= 14 digits. (PR #2032.
   Thanks to John Marshall. Reported on bioconda/bioconda-recipes#47137 by
   jmunoz94)

* Added file name and error message to 'error closing output file' error in
   `samtools sort`. (PR #2050, fixes #2049.  Thanks to Joshua C Randall).

* Fixed hard clip trimming issue in `ampliconclip` where right-hand side
   qualities were being removed from left-hand side trims. (PR #2053, fixes
   #2048.  Reported by Duda5)

* Fixed a bug in `samtools merge --template-coordinate` where the wrong
   heap was being tested. (PR #2062.  Thanks to Nils Homer.  Reported on
   ng-core/fastquorum#52 by David Mas-Ponte)

* Do not look at chr "*" for unmapped-placed reads with `samtools view
   --fetch-pairs`.  This was causing a significant slowdown when
   `--fetch-pairs` was being used. (PR #2070, fixes #2059.  Reported by
   acorvelo)

* Fixed bug which could cause `samtools view -L` to give incomplete output
   when the BED file contained nested target locations. (PR #2107, fixes
   #2104.  Reported by geertvandeweyer)

* Enable `samtools coverage` to handle alignments that do not have quality
   score data.  This was causing memory access problems. (PR #2083, fixes
   #2076.  Reported by Matthew Colpus)

* Fix undefined behaviour in `samtools fastq` with empty QUAL. (PR #2084)

* In `plot-bamstats` fixed read-length plot for data with limited variations
   in length. Lack of data was causing gnuplot problems. (PR #2085, fixes
   #2068.  Reported by mariyeta)

* Fixed an accidental fall-through that caused `samtools split -p` to also
   enable `--no-PG`. (PR #2101)

* Fixed an overflow that caused `samtools consensus -m simple` to give
   incorrect output when the input coverage exceeded several million reads
   deep. (PR #2099, fixes #2095.  Reported by Dylan Lawrence)

Non user-visible changes and build improvements:

* Work around address sanitizer going missing from the Cirrus CI ubuntu clang
   compiler by moving the address sanitizer build to gcc. Fix warnings from
   the new clang compiler. (PR #2043)

* Windows based CI has been migrated from AppVeyor to GitHub Actions. (PR
   #2072, PR #2108)

* Turn on more warning options in Cirrus-CI builds, ensure everything builds
   with `-Werror`, and add undefined behaviour checks to the address sanitizer
   test. (PR #2101, PR #2103, PR #2109)

* Tidy up Makefile dependencies and untracked test files. (PR #2106.  Thanks
   to John Marshall)

------------------------------------------------------------------------------
bcftools - changes v1.21
------------------------------------------------------------------------------

Changes affecting the whole of bcftools, or multiple commands:

* Support multiple semicolon-separated strings when filtering by ID
   using -i/-e (#2190).
   For example, `-i 'ID="rs123"'` now correctly matches `rs123;rs456`

* The filtering expression ILEN can be positive (insertion), negative
   (deletion), zero (balanced substitutions), or set to missing value
   (symbolic alleles).

* bcftools query
* bcftools +split-vep

     - The columns indices printed by default with `-H` (e.g., "#[1]CHROM")
       can be now suppressed by giving the option twice `-HH` (#2152)

Changes affecting specific commands:

* bcftools annotate

     - Support dynamic variables read from a tab-delimited annotation file
       (#2151). For example, in the two cases below the field 'STR' from the -a
       file is required to match the INFO/TAG in VCF. In the first example the
       alleles REF,ALT must match, in the second example they are ignored. The
       option -k is required to output also records that were not annotated:

        bcftools annotate -a ann.tsv.gz \
           -c CHROM,POS,REF,ALT,SCORE,~STR -i'TAG={STR}' -k in.vcf

        bcftools annotate -a ann.tsv.gz \
           -c CHROM,POS,-,-,SCORE,~STR     -i'TAG={STR}' -k in.vcf

     - When adding Type=String annotations from a tab-delimited file, encode
       characters with special meaning using percent encoding (';', '=' in
       INFO and ':' in FORMAT) (#2202)

* bcftools consensus

     - Allow to apply a reference allele which overlaps a previous deletion,
       there is no need to complain about overlapping alleles in such case

     - Fix a bug which required `-s -` to be present even when there were no
       samples in the VCF (#2260)

* bcftools csq

     - Fix a rare bug where indel combined with a substitution ending at
       exon boundary is incorrectly predicted to have 'inframe' rather
       than 'frameshift' consequence (#2212)

* bcftools gtcheck

     - Fix a segfault with --no-HWE-prob. The bug was introduced with the
       output format change in 1.19 which replaced the DC section with DCv2
       (#2180)

     - The number of matching genotypes in the DCv2 output was not calculated
       correctly with non-zero `-E, --error-probability`. Consequently, also
       the average HWE score was incorrect. The main output, the discordance
       score, was not affected by the bug

* bcftools +mendelian2

     - Include the number of good cases where at least one of the trio
       genotypes has an alternate allele (#2204)

     - Fix the error message which would report the wrong sample when
       non-existent sample is given. Note that bug only affected the error
       message, the program otherwise assigns the family members correctly
       (#2242)

* bcftools merge

     - Fix a severe bug in merging of FORMAT fields with Number=R and Number=A
       values. For example, rows with high-coverage FORMAT/AD values (bigger
       or equal to 128) could have been assigned to incorrect samples. The bug
       was introduced in version 1.19. For details see #2244.

* bcftools mpileup

     - Return non-zero error code when the input BAM/CRAM file is truncated
       (#2177)

     - Add FORMAT/AD annotation by default, disable with `-a -AD`

* bcftools norm

     - Support realignment of symbolic <DUP.*> alleles, similarly to <DEL.*>
       added previously (#1919,#2145)

     - Fix in reporting reference allele genotypes with `--multi-overlaps .`
       (#2160)

     - Support of duplicate removal of symbolic alleles of the same type but
       different SVLEN (#2182)

     - New `-S, --sort` switch to optionally sort output records by allele
       (#1484)

     - Add the `-i/-e` filtering options to select records for normalization.
       Note duplicate removal ignores this option.

     - Fix a bug where `--atomize` would not fill GT alleles for atomized SNVs
       followed by an indel (#2239)

* bcftools +remove-overlaps

     - Revamp the program to allow greater flexibility, with the following new
       options:

       -M, --mark-tag TAG   Mark -m sites with INFO/TAG

       -m, --mark EXPR      Mark (if also -M is present)
                            or remove sites [overlap]

             dup       .. all overlapping sites
             overlap   .. overlapping sites
             min(QUAL) .. mark sites with lowest QUAL until overlaps
                          are resolved

       --missing EXPR       Value to use for missing tags
                            with -m 'min(QUAL)'

             0   .. the default
             DP  .. heuristics, scale maximum QUAL value proportionally
                    to INFO/DP

       --reverse            Apply the reverse logic, for example preserve
                            duplicates instead of removing

       -O, --output-type t  t: plain list of sites (chr,pos),
                            tz: compressed list

* bcftools +tag2tag

     - The conversions --LXX-to-XX, --XX-to-LXX were working but specific
       cases such as --LAD-to-AD were not.

     - Print more informative error message when source tag type violiates VCF
       specification

* bcftools +trio-dnm2

     - Better handling of the --strictly-novel functionality, especically with
       respect to chrX inheritance


----------------------------------------------------------------------
The Wellcome Sanger Institute is operated by Genome Research Limited, a charity registered in England with number 1021457 and a company registered in England with number 2742969, whose registered office is Wellcome Sanger Institute, Wellcome Genome Campus, Hinxton, CB10 1SA.

[Samtools-help] Release 1.20

[Robert D. <rm...@sa...>]

Samtools (and HTSlib and BCFtools) version 1.20 is now available from
GitHub and SourceForge.

https://github.com/samtools/htslib/releases/tag/1.20
https://github.com/samtools/samtools/releases/tag/1.20
https://github.com/samtools/bcftools/releases/tag/1.20
https://sourceforge.net/projects/samtools/ 
The main changes are listed below:


------------------------------------------------------------------------------
htslib - changes v1.20
------------------------------------------------------------------------------

Updates
-------

* When working on named files, bgzip now sets the modified and access times
   of the output files it makes to match those of the corresponding input.
   (PR #1727, feature request #1718.  Requested by Gert Hulselmans)

* It's now possible to use a -o option to specify the output file name in
   bgzip. (PR #1747, feature request #1726.  Requested by Gert Hulselmans)

* Improved error faidx error messages. (PR #1743, thanks to Nick Moore)

* Faster reading of SAM array (type "B") tags.  These often turn up in ONT
   and PacBio data. (PR #1741)

* Improved validity checking of base modification tags. (PR #1749)

* mpileup overlap removal now works where one read has a deletion. (PR #1751,
   fixes samtools/samtools#1992.  Reported by Long Tian)

* The S3 plugin can now find buckets via S3 access point aliases. (PR #1756,
   thanks to Matt Pawelczyk; fixes samtools/samtools#1984.  Reported by
   Albert Li)

* Added a --threads option (and -@ short option) to tabix. (PR #1755, feature
   request #1735.  Requested by Dan Bolser)

* tabix can now index Graph Alignment Format (GAF) files. (See
   https://github.com/lh3/gfatools/blob/master/doc/rGFA.md) (PR #1763,
   thanks to Adam Novak)

Bug fixes
---------

* Security fix: Prevent possible heap overflow in cram_encode_aux() on bad
   RG:Z tags. (PR #1737)

* Security fix: Prevent attempts to call a NULL pointer if certain URL
   schemes are used in CRAM @SQ UR: tags. (PR #1757)

* Security fix: Fixed a bug where following certain AWS S3 redirects could
   downgrade the connection from TLS (i.e. https://) to unencrypted http://.
   This could happen when using path-based URLs and AWS_DEFAULT_REGION was set
   to a region other that the one where the data was stored. (PR #1762, fixes
   #1760. Reported by andaca)

* Fixed arithmetic overflow when loading very long references for CRAM. (PR
   #1738, fixes #1738.  Reported by Shane McCarthy)

* Fixed faidx and CRAM reference look-ups on compressed fasta where the .fai
   index file was present, but the .gzi index of compressed offsets was not.
   (PR #1745, fixes #1744.  Reported by Theodore Li)

* Fixed BCF indexing on-the-fly bug which produced invalid indexes when
   using multiple compression threads. (PR #1742, fixes #1740.  Reported
   by graphenn)

* Ensure that pileup destructors are called by bam_plp_destroy(), to prevent
   memory leaks. (PR #1749, PR #1754)

* Ensure on-the-fly index timestamps are always older than the data file.
   Previously the files could be closed out of order, leading to warnings
   being printed when using the index. (PR #1753, fixes #1732.  Reported by
   Gert Hulselmans)

* To prevent data corruption when reading (strictly invalid) VCF files with
   duplicated FORMAT tags, all but the first copy of the data associated with
   the tag are now dropped with a warning. (PR #1752, PR #1761, fixes #1733.
   Reported by anthakki)

* Fixed a bug introduced in release 1.19 (PR #1689) which broke variant
   record data if it tried to remove an over-long tag. (PR #1752, PR #1761)

* Changed error to warning when complaining about use of the CG tag in SAM or
   CRAM files. (PR #1758, fixes samtools/samtools#2002)

------------------------------------------------------------------------------
samtools - changes v1.20
------------------------------------------------------------------------------

* Added a `--max-depth` option to `bedcov`, for more control over the depth
   limit used when calculating the pileup.  Previously this limit was set at
   64000; now it is set to over 2 billion, so effectively all bases will be
   counted. (PR #1970, fixes #1950.  Reported by ellisjj)

* Added `mpileup --output-extra RLEN` to display the unclipped read length.
   (PR #1971, feature request #1959.  Requested by Feng Tian)

* Improved checking of symbolic flag names (e.g. UNMAP) passed to samtools.
   (PR #1981, fixes #1977.  Reported by Ilya Shlyakhter)

* The `samtools consensus --min-depth` option now works for the Bayesian mode
   as well as the simple one. (PR #1989, feature request #1982.  Requested by
   Gautier Richard)

* It's now possible to use the `samtools fastq` `-d tag:val` option multiple
   times, allowing matches on more than one tag/value.  It also gets a `-D`
   option which allows the values to be listed in a file. (PR #1993, feature
   request #1958.  Requested by Tristan Lefebure)

* Added `samtools fixmate` `-M` option to sanity check base modification
   (`ML`, `MM`, `MN`) tags, and where necessary adjust modification data on
   hard-clipped records. (PR #1990)

* Made `mpileup` run faster. (PR #1995)

* `samtools import` now adds a `@PG` header to the files it makes. As with
   other sub-commands, this can be disabled by using `--no-PG`. (PR #2008.
   Requested by Steven Leonard)

* The `samtools split` `-d` option to split by tag value now works on tags
   with integer values. (PR #2005, feature request #1956.  Requested by
   Alex Leonard)

* Adjusted `samtools sort -n` (by name) so that primary reads are always
   sorted before secondary / supplementary. (PR #2012, feature request
   #2010.  Requested by Stijn van Dongen)

* Added `samtools bedcov` `-H` option to print column headers in the output.
   (PR #2025.  Thanks to Dr. K. D. Murray)

Documentation:

* Added a note that BAQ is applied before filtering and overlap removal
   during mpileup processing. (PR #1988, fixes #1985.  Reported by
   Joseph Galasso)

* Added 3.1 to the list of supported CRAM versions in the samtools manual
   page. (PR #2009.  Thanks to Andrew Thrasher)

* Made assorted improvements to ampliconclip, flagstat and markdup manual
   pages. (PR #2014)

Bug Fixes:

* Security fix: Fixed double free that could occur if bed file indexing
   failed due to running out of memory.  This bug first appeared in version
   1.19.1. (PR #2026)

* Corrected error message printed when faidx fails to load the fai index. (PR
   #1987.  Thanks to Nick Moore)

* Fixed bug introduced in release 1.4 that caused incorrect reference bases
   to be printed by `samtools mpileup -a -f ref.fa` in the zero-depth regions
   at the end of each reference. (PR #2019, fixes #2018.  Reported by
   Joe Georgeson)

* Fixed a samtools view usage crash on MinGW when given invalid options. (PR
   #2030, fixes #2029.  Reported by Divon Lan)

Non user-visible changes and build improvements:

* Added tests to ensure that CRAM compression is working properly. (PR #1969,
   part of fix for #1968.  Reported by Clockris)

------------------------------------------------------------------------------
bcftools - changes v1.20
------------------------------------------------------------------------------

Changes affecting the whole of bcftools, or multiple commands:

* Add short option -W for --write-index. The option now accepts an optional
   parameter which allows to choose between TBI and CSI index format.

Changes affecting specific commands:

* bcftools consensus

     - Add new --regions-overlap option which allows to take into account
       overlapping deletions that start out of the fasta file target region.

* bcftools isec

     - Add new option `-l, --file-list` to read the list of file names from a
       file

* bcftools merge

     - Add new option `--force-single` to support single-file edge case
       (#2100)

* bcftools mpileup

     - Add new option --indels-cns for an alternative indel calling model,
       which should increase the speed on long read data (thanks to using
       edlib) and the precision (thanks to a number of heuristics).

* bcftools norm

     - Change the order of atomization and multiallelic splitting (when both
       -a,-m are given) from "atomize first, then split" to "split first, then
       atomize". This usually results in a simpler VCF representation. The
       previous behaviour can be achieved by explicitly streaming the output
       of the --atomize command into the --multiallelics splitting command.

     - Fix Type=String multiallelic splitting for Number=A,R,G tags with
       incorrect number of values.

     - Merging into multiallelic sites with `bcftools norm -m +indels` did
       not work. This is now fixed and the merging is now more strict about
       variant types, for example complex events, such as AC>TGA, are not
       considered as indels anymore (#2084)

* bcftools reheader

     - Allow reading the input file from a stream with --fai (#2088)

* bcftools +setGT

     - Support for custom genotypes based on the allele with higher depth,
       such as `--new-gt c:0/X` custom genotypes (#2065)

* bcftools +split-vep

     - When only one of the tags is present, automatically choose INFO/BCSQ
       (the default tag name produced by `bcftools csq`) or INFO/CSQ (produced
       by VEP). When both tags are present, use the default INFO/CSQ.

     - Transcript selection by MANE, PICK, and user-defined transcripts, for
       example

        --select CANONICAL=YES
        --select MANE_SELECT!=""
        --select PolyPhen~probably_damaging

     - Select all matching transcripts via --select, not just one

     - Change automatic type parsing of VEP fields DNA_position, CDS_position,
       and Protein_position from Integer to String, as it can be of the form
       "8586-8599/9231". The type Integer can be still enforced with
       `-c cDNA_position:int,CDS_position:int,Protein_position:int`.

     - Recognize `-c field:str`, not just `-c field:string`, as advertised in
       the usage page

     - Fix a bug which made filtering expression containing missing values
       crash (#2098)

* bcftools stats

     - When GT is missing but AD is present, the program determines the
       alternate allele from AD. However, if the AD tag has incorrect number
       of values, the program would exit with an error printing "Requested
       allele outside valid range". This is now fixed by taking into account
       the actual number of ALT alleles.

* bcftools +tag2tag

     - Support for conversion from tags using localized alleles (e.g. LPL,
       LAD) to the family of standard tags (PL, AD)

* bcftools +trio-dnm2

     - Extend --strictly-novel to exclude cases where the non-Mendelian
       allele is the reference allele. The change is motivated by the
       observation that this class of variants is enriched for errors
       (especially for indels), and better corresponds with the option name.



-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is Wellcome Sanger Institute, Wellcome Genome Campus, 
 Hinxton, CB10 1SA.

[Samtools-help] SAMtools release 1.19.2

[Robert D. <rm...@sa...>]

Samtools release 1.19.2 is now available from GitHub and SourceForge.  It 
fixes an error in 1.19.1 that broke filtering on unordered BED files.

https://github.com/samtools/samtools/releases/tag/1.19.2
https://sourceforge.net/projects/samtools/

------------------------------------------------------------------------------
samtools - changes v1.19.2
------------------------------------------------------------------------------

Bug Fixes:

* Fixed a regression in 1.19.1 that broke BED filtering for inputs where the
   region start positions for the same reference were not sorted in ascending
   order. (PR #1975, fixes #1974.  Reported by Anže Starič)



-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is Wellcome Sanger Institute, Wellcome Genome Campus, 
 Hinxton, CB10 1SA.

[Samtools-help] SAMtools and HTSlib release 1.19.1

[Robert D. <rm...@sa...>]

Samtools and HTSlib release 1.19.1 are now available from
GitHub and SourceForge.  This release fixes a regression
in version 1.19 that caused written CRAM files to be much
bigger than necessary.  It also fixes a number of other
bugs listed below, including one that could cause crashes
or incorrect results when filtering by regions in a BED
file.

https://github.com/samtools/htslib/releases/tag/1.19.1
https://github.com/samtools/samtools/releases/tag/1.19.1
https://sourceforge.net/projects/samtools/

------------------------------------------------------------------------------
htslib - changes v1.19.1
------------------------------------------------------------------------------

* Fixed a regression in release 1.19 that caused all aux records to be stored
   uncompressed in CRAM files.  The resulting files were correctly formatted,
   but bigger than they needed to be. (PR#1729, fixes samtools#1968.  Reported
   by Clockris)

* Fixed possible out-of-bounds reads due to an incorrect check on B tag
   lengths in cram_encode_aux().  (PR#1725)

* Fixed an incorrect check on tag length which could fail to catch a two byte
   out-of-bounds read in bam_get_aux(). (PR#1728)

* Made errors reported by hts_open_format() less confusing when it can't open
   the reference file.  (PR#1724, fixes #1723.  Reported by Alex Leonard)

* Made hts_close() fail more gracefully if it's passed a NULL pointer
   (PR#1724)

------------------------------------------------------------------------------
samtools - changes v1.19.1
------------------------------------------------------------------------------

Bug Fixes:

* Fixed a possible array bounds violation when looking up regions in a BED
   file (e.g. using `samtools view -L`).  This could lead to crashes or the
   return of incomplete results if the BED file contained a large number of
   entries all referencing low positions on a chromosome. (PR #1962, fixes
   #1961.  Reported by geertvandeweyer)

* Fixed a crash in samtools stats that occurred when trying to clean up
   after it was unable to open a CRAM reference file. (PR #1957, fixes
   crash reported in samtools/htslib#1723.  Reported by Alex Leonard)

Documentation:

* Fixed inverted logic in the `samtools consensus --show-del` manual page
   description. (PR #1955, fixes #1951.  Reported by Mikhail Schelkunov)

* Added a description of the MPC section to the `samtools stats` manual page.
   (PR #1963, fixes #1954.  Reported by litun-fkby)



-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is Wellcome Sanger Institute, Wellcome Genome Campus, 
 Hinxton, CB10 1SA.

[Samtools-help] Release 1.19

[Robert D. <rm...@sa...>]

Samtools (and HTSlib and BCFtools) version 1.19 is now available from
GitHub and SourceForge.

https://github.com/samtools/htslib/releases/tag/1.19
https://github.com/samtools/samtools/releases/tag/1.19
https://github.com/samtools/bcftools/releases/tag/1.19
https://sourceforge.net/projects/samtools/

The main changes are listed below:

------------------------------------------------------------------------------
htslib - changes v1.19
------------------------------------------------------------------------------

Updates
-------

* A temporary work-around has been put in the VCF parser so that it is less
   likely to fail on rows with a large number of ALT alleles, where Number=G
   tags like PL can expand beyond the 2Gb limit enforced by HTSlib.  For now,
   where this happens the offending tag will be dropped so the data can be
   processed, albeit without the likelihood data.

   In future work, the library will instead convert such tags into their
   local alternatives (see https://github.com/samtools/hts-specs/pull/434).

* New program. Adds annot-tsv which annotates regions in a destination file
   with texts from overlapping regions in a source file. (PR#1619)

* Change bam_parse_cigar() so that it can modify existing BAM records.  This
   makes more useful as public API.  Previously it could only handle partially
   formed BAM records. (PR#1651, fixes #1650. Reported by Oleksii Nikolaienko)

* Add "uncompressed" to hts_format_description() where appropriate.  This
   adds an "uncompressed" description to uncompressed files that would
   normally be compressed, such as BAM and BCF. (PR#1656, in relation to
   samtools#1884.  Thanks to John Marshall)

* Speed up to the VCF parser and writer. (PR#1644 and PR#1663)

* Add an hclen (hard clip length) SAM filter function. (PR#1660, with
   reference to samtools#813)

* Avoid really closing stdin/stdout in hclose()/hts_close()/et al. See
   discussion in PR for details. (PR#1665.  Thanks to John Marshall)

* Add support to handle multiple files in bgzip. (PR#1658, fixes #1642.
   Requested by bw2)

* Enable auto-vectorisation in CRAM 3.1 codecs.  Speeds decoding on some
   sequencing platform data. (PR#1669)

* Speed up removal of lines in large headers. (PR#1662, fixes #1460.
   Reported by Anže Starič)

* Apply seqtk PR to improve kseq.h parsing performance.  Port of
   Fabian Klötzl's (kloetzl) lh3/seqtk#123 and attractivechaos/klib#173
   to HTSlib. (PR#1674.  Thanks to John Marshall)

Build changes
-------------

* Updated htscodecs submodule to 1.6.0. (PR#1685, PR#1717, PR#1719)

* Apply the packed attribute to uint*_u types for Clang to prevent
   -fsanitize=alignment failures. (PR#1667.  Thanks to Fangrui Song)

* Fuzz testing improvements. (PR#1664)

* Add C++ casts for external headers in klist.h and kseq.h. (PR#1683.  See
   also PR#1674 and PR#1682)

* Add test case compiling the public headers as C++. (PR#1682.  Thanks to
   John Marshall)

* Enable optimisation level -O3 for SAM QUAL+33 formatting. (PR#1679)

* Make compiler flag detection work with zig cc. (PR#1687)

* Fix unused value warnings when built with NDEBUG. (PR#1688)

* Remove some disused Makefile variables, fix typos and a warning.  Improve
   bam_parse_basemod() documentation. (PR#1705, Thanks to John Marshall)

Bug fixes
---------

* Fail bgzf_useek() when offset is above block limits. (PR#1668)

* Fix multi-threaded on-the-fly indexing problems. (PR#1672, fixes
   samtools#1861 and bcftools#1985.  Reported by Mark Ebbert and lacek)

* Fix hfile_libcurl small seek bug. (PR#1676, fixes samtools#1918.  Also may
   fix #1037, #1625 and samtools#1622. Reported by Alex Reynolds, Mark Walker,
   Arthur Gilly and skatragadda-nygc. Thanks to John Marshall)

* Fix a minor memory leak in malformed CRAM EXTERNAL blocks. [fuzz] (PR#1671)

* Fix a cram decode hang from block_resize(). (PR#1680. Reported by
   Sebastian Deorowicz)

* Cram fuzzing improvements.  Fixes a number of cram errors. (PR#1701,
   fixes #1691, #1692, #1693, #1696, #1697, #1698, #1699 and #1700. Thanks
   to Octavio Galland for finding and reporting all these)

* Fix crypt4gh redirection. (PR#1675, fixes  grbot/crypt4gh-tutorial#2.
   Reported by hth4)

* Fix PG header linking when records make a loop. (PR#1702, fixes #1694.
   Reported by Octavio Galland)

* Prevent issues with no-stored-sequence records in CRAM files, by ensuring
   they are accounted for properly in block size calculations, and by limiting
   the maximum query length in the CIGAR data.  Originally seen as an overflow
   by OSS-Fuzz / UBSAN, it turned out this could lead to excessive time and
   memory use by HTSlib, and could result in it writing out unreadable CRAM
   files. (PR#1710)

* Fix some illegal shifts and integer overflows found by OSS-Fuzz / UBSAN.
   (PR#1707, PR#1712, PR#1713)

------------------------------------------------------------------------------
samtools - changes v1.19
------------------------------------------------------------------------------

New work and changes:

* Samtools coverage: add a new --plot-depth option to draw depth (of
   coverage) rather than the percentage of bases covered. (PR #1910.
   Thanks to Pierre Lindenbaum)

* Samtools merge / sort: add a lexicographical name-sort option via the -N
   option.  The "natural" alpha-numeric sort is still available via -n.
   (PR #1900, fixes #1500.  Reported by Steve Huang)

* Samtools view: add -N ^NAME_FILE and -R ^RG_FILE options.  The standard -N
   and -R options only output reads matching a specified list of names or
   read-groups.  With a caret (^) prefix these may be negated to only output
   read not matching the specified files. (PR #1896, fixes #1895.  Suggested
   by Feng Tian)

* Cope with Htslib's change to no longer close stdout on hts_close. Htslib
   companion PR is samtools/htslib#1665. (PR #1909.  Thanks to John Marshall)

* Plot-bamstats: add a new plot of the read lengths ("RL") from samtools
   stats output. (PR #1922, fixes #1824.  Thanks to @erboone, suggested by
   Martin Pollard)

* Samtools split: support splitting files based on the contents of auxiliary
   tags.  Also adds a -M option to limit the number of files split can make,
   to avoid accidental resource over-allocation, and fixes some issues with
   --write-index. (PR #1222, PR #1933, fixes #1758.  Thanks to Valeriu Ohan,
   suggested by Scott Norton)

Bug Fixes:

* Samtools stats: empty barcode tags are now treated as having no barcode.
   (PR #1929, fixes #1926.  Reported by Jukka Matilainen)

* Samtools cat: add support for non-seekable streams.  The file format
   detection blocked pipes from working before, but now files may be
   non-seekable such as stdin or a pipe. (PR #1930, fixes #1731.  Reported
   by Julian Hess)

* Samtools mpileup -aa (absolutely all positions) now produces an output even
   when given an empty input file. (PR #1939.  Reported by Chang Y)

* Samtools markdup: speed up optical duplicate tagging on regions with very
   deep data. (PR #1952)

Documentation:

* Samtools mpileup: add more usage examples to the man page. (PR #1913,
   fixes #1801)

* Samtools fastq: explicitly document the order that filters apply.
   (PR #1907)

* Samtools merge: fix example output to use an uppercase RG PL field.
   (PR #1917.  Thanks to John Marshall.  Reported by Michael Macias)

* Add hclen SAM filter documentation. (PR #1902.  See also
   samtools/htslib#1660)

* Samtools consensus: remove the -5 option from documentation.  This option
   was renamed before the consensus subcommand was merged, but accidentally
   left in the man page. (PR #1924)

------------------------------------------------------------------------------
bcftools - changes v1.19
------------------------------------------------------------------------------

Changes affecting the whole of bcftools, or multiple commands:

* Filtering expressions can be given a file with list of strings to match,
   this was previously possible only for the ID column. For example

   ID=@file            .. selects lines with ID present in the file
   INFO/TA...@fi...  .. selects lines where TAG has a string value listed
                          in the file
   INFO/TAG!=@file.txt .. TAG must not have a string value
                          listed in the file

* Allow to query REF,ALT columns directly, for example

   -e 'REF="N"'

Changes affecting specific commands:

* bcftools annotate

     - Fix `bcftools annotate --mark-sites`, VCF sites overlapping regions in
       a BED file were not annotated (#1989)

     - Add flexibility to FILTER column transfers and allow transfers within
       the same file, across files, and in combination. For examples see
       http://samtools.github.io/bcftools/howtos/annotate.html#transfer_filter_to_info

* bcftools call

     - Output MIN_DP rather than MinDP in gVCF mode

     - New `-*, --keep-unseen-allele` option to output the unobserved
       allele <*>, intended for gVCF.

* bcftools head

     - New `-s, --samples` option to include the #CHROM header line with
       samples.

* bcftools gtcheck

     - Add output options `-o, --output` and `-O, --output-type`

     - Add filtering options `-i, --include` and `-e, --exclude`

     - Rename the short option `-e, --error-probability` from lower case to
       upper case `-E, --error-probability`

     - Changes to the output format, replace the DC section with DCv2:

         - adds a new column for the number of matching genotypes

         - The --error-probability is newly interpreted as the probability of
           erroneous allele rather than genotype. In other words, the
           calculation of the discordance score now considers the probability
           of genotyping error to be different for HOM and HET genotypes,
           i.e. P(0/1|dsg=0) > P(1/1|dsg=0).

         - fixes in HWE score calculation plus output average HWE score rather
           than absolute HWE score

         - better description of fields

* bcftools merge

     - Add `-m` modifiers to suppress the output of the unseen allele <*> or
       <NON_REF> at variant sites (e.g. `-m both,*`) or all sites
       (e.g. `-m both,**`)

* bcftools mpileup

     - Output MIN_DP rather than MinDP in gVCF mode

* bcftools norm

     - Add the number of joined lines to the summary output, for example

       Lines   total/split/joined/realigned/skipped:  6/0/3/0/0

     - Allow combining -m and -a with --old-rec-tag (#2020)

     - Symbolic <DEL> alleles caused norm to expand REF to the full length of
       the deletion. This was not intended and problematic for long deletions,
       the REF allele should list one base only (#2029)

* bcftools query

     - Add new `-N, --disable-automatic-newline` option for pre-1.18 query
       formatting behavior when newline would not be added when missing

     - Make the automatic addition of the newline character in a more
       predictable way and, when missing, always put it at the end of the
       expression. In version 1.18 it could be added at the end of the
       expression (for per-site expressions) or inside the square brackets
       (for per-sample expressions). The new behavior is:

         - if the formatting expression contains a newline character, do
           nothing

         - if there is no newline character and -N,
           --disable-automatic-newline is given, do nothing

         - if there is no newline character and -N is not given, insert
           newline at the end of the expression

  See #1969 for details

     - Add new `-F, --print-filtered` option to output a default string for
       samples that would otherwise be filtered by `-i/-e` expressions.

     - Include sample name in the output header with `-H` whenever it makes
       sense (#1992)

* bcftools +spit-vep

     - Fix on the fly filtering involving numeric subfields, e.g. `-i
       'MAX_AF<0.001'` (#2039)

     - Interpret default column type names (--columns-types) as entire
       strings, rather than substrings to avoid unexpected spurious
       matches (i.e. internally add ^ and $ to all field names)

* bcftools +trio-dnm2

     - Do not flag paternal genotyping errors as de novo mutations.
       Specifically, when father's chrX genotype is 0/1 and mother's 0/0,
       0/1 in the child will not be marked as DNM.

* bcftools view

     - Add new `-A, --trim-unseen-allele` option to remove the unseen allele
       <*> or <NON_REF> at variant sites (`-A`) or all sites (`-AA`)



-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is Wellcome Sanger Institute, Wellcome Genome Campus, 
 Hinxton, CB10 1SA.

Re: [Samtools-help] new user

[Tom B. <tb...@um...>]

Chris  -

There must be oodles of journal articles over the years since about 2010 in
the biology literatrure describing this analysis for non-human species.  I
have no experience myself with non-human species.  A very vanilla workflow
using bwa-mem and samtools is briefly described at
https://www.htslib.org/workflow/wgs-call.html.  The Broad Institute GATK
suite of tools has very extensive and detailed documentation.
https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels-
is one example.  But as an introduction, the right journal article will
probably be more useful than either of these technical web pages.

-  thanks  -  tom blackwell  -



On Tue, Oct 24, 2023 at 6:10 AM Christopher Blair <cb...@gc...>
wrote:

> Hi all,
>
> I was hoping that someone could point me in the right direction for
> processing WGS data for phylogenetics and population genetics. I will soon
> have 10x whole genome data for about 27 individuals of a lizard species. I
> will need to go from raw reads to multiple sequence alignments and
> SNPs/genotypes. There is a reference genome available from a closely
> related species. Any guidance would be appreciated.
>
>
> Chris
>
>
>
>
> <{===================~~~~~~~~~~~
>           /\                      /\
> Christopher Blair, Ph.D.
> Executive Officer, Biology PhD Program
>
> Professor
>
> Graduate Center and
>
> New York City College of Technology
> The City University of New York
> cb...@ci... <cb...@gc...>; cb...@gc...
> Website: https://sites.google.com/site/christopherblairphd/home
> Office: 4315 (GC); Academic Complex 501D (City Tech)
> Lab: Academic Complex 401 (City Tech)
> _______________________________________________
> Samtools-help mailing list
> Sam...@li...
> https://lists.sourceforge.net/lists/listinfo/samtools-help
>

[Samtools-help] new user

[Christopher B. <cb...@gc...>]

Hi all,

I was hoping that someone could point me in the right direction for processing WGS data for phylogenetics and population genetics. I will soon have 10x whole genome data for about 27 individuals of a lizard species. I will need to go from raw reads to multiple sequence alignments and SNPs/genotypes. There is a reference genome available from a closely related species. Any guidance would be appreciated.


Chris





<{===================~~~~~~~~~~~
          /\                      /\
Christopher Blair, Ph.D.
Executive Officer, Biology PhD Program

Professor

Graduate Center and

New York City College of Technology
The City University of New York
cb...@ci...<mailto:cb...@gc...>; cb...@gc...<mailto:cb...@gc...>
Website: https://sites.google.com/site/christopherblairphd/home
Office: 4315 (GC); Academic Complex 501D (City Tech)
Lab: Academic Complex 401 (City Tech)

[Samtools-help] Release 1.18

[Robert D. <rm...@sa...>]

Samtools (and HTSlib and BCFtools) version 1.18 is now available from
GitHub and SourceForge.

https://github.com/samtools/htslib/releases/tag/1.18
https://github.com/samtools/samtools/releases/tag/1.18
https://github.com/samtools/bcftools/releases/tag/1.18
https://sourceforge.net/projects/samtools/ 
The main changes are listed below:

------------------------------------------------------------------------------
htslib - changes v1.18
------------------------------------------------------------------------------

Updates
-------

* Using CRAM 3.1 no longer gives a warning about the specification being
   draft.  Note CRAM 3.0 is still the default output format. (PR#1583)

* Replaced use of sprintf with snprintf, to silence potential warnings from
   Apple's compilers and those who implement similar checks. (PR#1594,
   fixes #1586. Reported by Oleksii Nikolaienko)

* Fastq output will now generate empty records for reads with no sequence
   data (i.e. sequence is "*" in SAM format). (PR#1576, fixes
   samtools/samtools#1576.  Reported by Nils Homer)

* CRAM decoding speed-ups. (PR#1580)

* A new MN aux tag can now be used to verify that MM/ML base modification
   data has not been broken by hard clipping. (PR#1590, PR#1612. See also PR
   samtools/hts-specs#714 and issue samtools/hts-specs#646. Reported by
   Jared Simpson)

* The base modification API has been improved to make it easier for
   callers to tell unchecked bases from unmodified ones. (PR#1636, fixes
   #1550. Requested by Chris Wright)

* A new bam_mods_queryi() API has been added to return additional data about
   the i-th base modification returned by bam_mods_recorded(). (PR#1636, fixes
   #1550 and #1635.  Requested by Jared Simpson)

* Speed up index look-ups for whole-chromosome queries. (PR#1596)

* Mpileup now merges adjacent (mis)match CIGAR operations, so CIGARs using
   the X/= operators give the same results as if the M operator was used.
   (PR#1607, fixes #1597.  Reported by Marcel Martin)

* It's now possible to call bcf_sr_set_regions() after adding readers using
   bcf_sr_add_reader() (previously this returned an error).  Doing so will
   discard any unread data, and reset the readers so they iterate over the
   new regions.  (PR#1624, fixes samtools/bcftools#1918.  Reported by
   Gregg Thomas)

* The synced BCF reader can now accept regions with reference names
   including colons and hyphens, by enclosing them in curly braces.  For
   example, {chr_part:1-1001}:10-20 will return bases 10 to 20 from
   reference "chr_part:1-1001".  (PR#1630, fixes #1620.  Reported by Bren)

* Add a "samples" directory with code demonstrating usage of HTSlib plus a
   tutorial document. (PR#1589)

Build changes
-------------

* Htscodecs has been updated to 1.5.1 (PR#1654)

* Htscodecs SIMD code now works with Apple multiarch binaries. (PR#1587,
   HTSlib fix for samtools/htscodecs#76. Reported by John Marshall)

* Improve portability of "expr" usage in version.sh. (PR#1593, fixes #1592.
   Reported by John Marshall)

* Improve portability to *BSD targets by ensuring _XOPEN_SOURCE is
   defined correctly and that source files properly include "config.h".
   Perl scripts also now all use #!/usr/bin/env instead of assuming that
   it's in /usr/bin/perl. (PR#1628, fixes #1606. Reported by
   Robert Clausecker)

* Fixed NAME entry in htslib-s3-plugin man page so the whatis and apropos
   commands find it.  (PR#1634, thanks to Étienne Mollier)

* Assorted dependency tracking fixes.  (PR#1653, thanks to John Marshall)

Documentation updates
---------------------

* Changed Alpine build instructions as they've switched back to using
   openssl. (PR#1609)

* Recommend using -rdynamic when statically linking a libhts.a with plugins
   enabled.  (PR#1611, thanks to John Marshall.  Fixes #1600, reported by
   Jack Wimberley)

* Fixed example in docs for sam_hdr_add_line().  (PR#1618, thanks to kojix2)

* Improved test harness for base modifications API. (PR#1648)

Bug fixes
---------

* Fix a major bug when searching against a CRAM index where one container
   has start and end coordinates entirely contained within the previous
   container. This would occasionally miss data, and sometimes return much
   more than required.  The bug affected versions 1.11 to 1.17, although the
   change in 1.11 was bug-fixing multi-threaded index queries.  This bug did
   not affect index building.  There is no need to reindex your CRAM files.
   (PR#1574, PR#1640. Fixes #1569, #1639, samtools/samtools#1808,
   samtools/samtools#1819.  Reported by xuxif, Jens Reeder and Jared Simpson)

* Prevent CRAM blocks from becoming too big in files with short sequences but
   very long aux tags.  (PR #1613)

* Fix bug where the CRAM decoder for CONST_INT and CONST_BYTE codecs may
   incorrectly look for extra data in the CORE block. Note that this bug only
   affected the experimental CRAM v4.0 decoder. (PR#1614)

* Fix crypt4gh redirection so it works in conjunction with non-file IO, such
   as using htsget. (PR#1577)

* Improve error checking for the VCF POS column, when facing invalid data.
   (PR#1575, replaces #1570 originally reported and fixed by Colin Nolan.)

* Improved error checking on VCF indexing to validate the data is BGZF
   compressed. (PR#1581)

* Fix bug where bin number calculation could overflow when making iterators
   over regions that go to the end of a chromosome. (PR#1595)

* Backport attractivechaos/klib#78 (by Pall Melsted) to HTSlib. Prevents
   infinite loops in kseq_read() when reading broken gzip files. (PR#1582,
   fixes #1579.  Reported by Goran Vinterhalter)

* Backport attractivechaos/klib@384277a (by innoink) to HTSlib. Fixes the
   kh_int_hash_func2() macro definition. (PR#1599, fixes #1598.  Reported by
   fanxinping)

* Remove a compilation warning on systems with newer libcurl releases.
   (PR#1572)

* Windows: Fixed BGZF EOF check for recent MinGW releases. (PR#1601, fixes
   samtools/bcftools#1901)

* Fixed bug where tabix would not return the correct regions for files
   where the column ordering is end, ..., begin instead of begin, ...,
   end. (PR#1626, fixes #1622.  Reported by Hiruna Samarakoon)

* sam_format_aux1() now always NUL-terminates Z/H tags. (PR#1631)

* Ensure base modification iterator is reset when no MM tag is present.
   (PR#1631, PR#1647)

* Fix segfault when attempting to write an uncompressed BAM file opened
   using hts_open(name, "wbu").  This was attempting to write BAM data
   without wrapping it in BGZF blocks, which is invalid according to the BAM
   specification.  "wbu" is now internally converted to "wb0" to output
   uncompressed data wrapped in BGZF blocks. (PR#1632, fixes #1617. Reported
   by Joyjit Daw)

* Fixed over-strict bounds check in probaln_glocal() which caused it to make
   sub-optimal alignments when the requested band width was greater than the
   query length.  (PR#1616, fixes #1605.  Reported by Jared Simpson)

* Fixed possible double frees when handling errors in bcf_hdr_add_hrec(), if
   particular memory allocations fail. (PR#1637)

* Ensure that bcf_hdr_remove() clears up all pointers to the items removed
   from dictionaries.  Failing to do this could have resulted in a call
   requesting a deleted item via bcf_hdr_get_hrec() returning a stale
   pointer. (PR#1637)

* Stop the gzip decompresser from finishing prematurely when an empty gzip
   block is followed by more data. (PR#1643, PR#1646)

------------------------------------------------------------------------------
samtools - changes v1.18
------------------------------------------------------------------------------

New work and changes:

* Add minimiser sort option to collate by an indexed fasta.  Expand the
   minimiser sort to arrange the minimiser values in the same order as
   they occur in the reference genome. This is acts as an extremely crude
   and simplistic read aligner that can be used to boost read compression.
   (PR#1818)

* Add a --duplicate-count option to markdup.  Adds the number of duplicates
   (including itself) to the original read in a 'dc' tag. (PR#1816. Thanks to
   wulj2)

* Make calmd handle unaligned data or empty files without throwing an error.
   This is to make pipelines work more smoothly.  A warning will still be
   issued. (PR#1841, fixes #1839.  Reported by Filipe G. Vieira)

* Consistent, more comprehensive flag filtering for fasta/fastq.  Added
   --rf/--incl[ude]-flags and long options for -F (--excl[ude]-flags and -f
   (--require-flags). (PR#1842.  Thanks to Devang Thakkar)

* Apply fastq --input-fmt-option settings.  Previously any options
   specified were not being applied to the input file. (PR#1855.  Thanks
   to John Marshall)

* Add fastq -d TAG[:VAL] check.  This mirrors view -d and will only output
   alignments that match TAG (and VAL if specified). (PR#1863, fixes #1854.
   Requested by Rasmus Kirkegaard)

* Extend import --order TAG to --order TAG:length.  If length is specified,
   the tag format goes from integer to a 0-padded string format.  This is a
   workaround for BAM and CRAM that cannot encode an order tag of over 4
   billion records. (PR#1850, fixes #1847.  Reported by Feng Tian)

* New -aa mode for consensus.  This works like the -aa option in depth
   and mpileup. The single 'a' reports all bases in contigs covered by
   alignments. Double 'aa' (or '-a -a') reports Ns even for the references
   with no alignments against them. (PR#1851, fixes #1849.  Requested by
   Tim Fennell)

* Add long option support to samtools index. (PR#1872, fixes #1869.  Reported
   by Jason Bacon)

* Be consistent with rounding of "average length" in samtools stats.
   (PR#1876, fixes #1867.  Reported by Jelinek-J)

* Add option to ampliconclip that marks reads as unmapped when they do not
   have enough aligned bases left after clipping.  Default is to unmap reads
   with zero aligned bases. (PR#1865, fixes #1856.  Requested by ces)

Bug Fixes:

* [From HTSLib] Fix a major bug when searching against a CRAM index where
   one container has start and end coordinates entirely contained within the
   previous container. This would occasionally miss data, and sometimes
   return much more than required.  The bug affected versions 1.11 to 1.17,
   although the change in 1.11 was bug-fixing multi-threaded index queries.
   This bug did not affect index building.  There is no need to reindex your
   CRAM files. (PR#samtools/htslib#1574, PR#samtools/htslib#1640. Fixes
   #samtools/htslib#1569, #samtools/htslib#1639, #1808, #1819.  Reported
   by xuxif, Jens Reeder and Jared Simpson)

* Fix a sort -M bug (regression) when merging sub-blocks.  Data was valid but
   in a poor order for compression. (PR#1812)

* Fix bug in split output format.  Now SAM and CRAM format can chosen as well
   as BAM.  Also a documentation change, see below. (PR#1821)

* Add error checking to view -e filter expression code.  Invalid expressions
   were not returning an error code. (PR#1833, fixes #1829.  Reported by
   Steve Huang)

* Fix reheader CRAM output version.  Sets the correct CRAM output version for
   non-3.0 CRAMs. (PR#1868, fixes #1866.  Reported by John Marshall)

Documentation:

* Expand the default filtering information on the mpileup man page. (PR#1802,
   fixes #1801.  Reported by gevro)

* Add an explanation of the default behaviour of split files on generating a
   file for reads with missing or unrecognised RG tags.  Also a small bug fix,
   see above. (PR#1821, fixes #1817.  Reported by Steve Huang)

* In the INSTALL instructions, switched back to openssl for Alpine.  This
   matches the current Alpine Linux practice. (PR#1837, see htslib#1591.
   Reported by John Marshall)

* Fix various typos caught by lintian parsers. (PR#1877.  Thanks to
   Étienne Mollier)

* Document consensus --qual-calibration option. (PR#1880, fixes #1879.
   Reported by John Marshall)

* Updated the page about samtools duplicate marking with more detail at
   www.htslib.org/algorithms/duplicate.html

Non user-visible changes and build improvements:

* Removed a redundant line that caused a warning in gcc-13. (PR#1838)

------------------------------------------------------------------------------
bcftools - changes v1.18
------------------------------------------------------------------------------

Changes affecting the whole of bcftools, or multiple commands:

* Support auto indexing during writing BCF and VCF.gz via new `--write-index`
   option

Changes affecting specific commands:

* bcftools annotate

     - The `-m, --mark-sites` option can be now used to mark all sites without
       the need to provide the `-a` file (#1861)

     - Fix a bug where the `-m` function did not respect the `--min-overlap`
       option (#1869)

     - Fix a bug when update of INFO/END results in assertion error (#1957)

* bcftools concat

     - New option `--drop-genotypes`

* bcftools consensus

     - Support higher-ploidy genotypes with `-H, --haplotype` (#1892)

     - Allow `--mark-ins` and `--mark-snv` with a character, similarly to
       `--mark-del`

* bcftools convert

     - Support for conversion from tab-delimited files (CHROM,POS,REF,ALT) to
       sites-only VCFs

* bcftools csq

     - New `--unify-chr-names` option to automatically unify different
       chromosome naming conventions in the input GFF, fasta and VCF files
       (e.g. "chrX" vs "X")

     - More versatility in parsing various flavors of GFF

     - A new `--dump-gff` option to help with debugging and investigating the
       internals of hGFF parsing

     - When printing consequences in nonsense mediated decay transcripts,
       include 'NMD_transcript' in the consequence part of the annotation.
       This is to make filtering easier and analogous to VEP annotations.
       For example the consequence annotation
       3_prime_utr|PCGF3|ENST00000430644|NMD is newly printed as
       3_prime_utr&NMD_transcript|PCGF3|ENST00000430644|NMD

* bcftools gtcheck

     - Add stats for the number of sites matched in the GT-vs-GT, GT-vs-PL,
       etc modes. This information is important for interpretation of the
       discordance score, as only the GT-vs-GT matching can be interpreted as
       the number of mismatching genotypes.

* bcftools +mendelian2

     - Fix in command line argument parsing, the `-p` and `-P` options were
       not functioning (#1906)

* bcftools merge

     - New `-M, --missing-rules` option to control the behavior of merging of
       vector tags to prevent mixtures of known and missing values in tags
       when desired

     - Use values pertaining to the unknown allele (<*> or <NON_REF>) when
       available to prevent mixtures of known and missing values (#1888)

     - Revamped line matching code to fix problems in gVCF merging where split
       gVCF blocks would not update genotypes (#1891, #1164).

* bcftool mpileup

     - Fix a bug in --indels-v2.0 which caused an endless loop when CIGAR
       operator 'H' or 'P' was encountered

* bcftools norm

     - The `-m, --multiallelics +` mode now preserves phasing (#1893)

     - Symbolic <DEL.*> alleles are now normalized too (#1919)

     - New `-g, --gff-annot` option to right-align indels in forward
       transcripts to follow HGVS 3'rule (#1929)

* bcftools query

     - Force newline character in formatting expression when not given
       explicitly

     - Fix `-H` header output in formatting expressions containing newlines

* bcftools reheader

     - Make `-f, --fai` aware of long contigs not representable by 32-bit
       integer (#1959)

* bcftools +split-vep

     - Prevent a segfault when `-i/-e` use a VEP subfield not included in `-f`
       or `-c` (#1877)

     - New `-X, --keep-sites` option complementing the existing `-x,
       --drop-sites` options

     - Force newline character in formatting expression when not given
       explicitly

     - Fix a subtle ambiguity: identical rows must be returned when `-s` is
       applied regardless of `-f` containing the `-a` VEP tag itself or not.

* bcftools stats

     - Collect new VAF (variant allele frequency) statistics from FORMAT/AD
       field

     - When counting transitions/transversions, consider also alternate het
       genotypes

* plot-vcfstats

     - Add three new VAF plots



-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is Wellcome Sanger Institute, Wellcome Genome Campus, 
 Hinxton, CB10 1SA.

Re: [Samtools-help] support to CRAM gzip [EXT]

[James B. <jk...@sa...>]

On Wed, Jul 05, 2023 at 03:11:50PM +0200, Adrian wrote:
> I noticed samtools works perfectly fine with gzip SAM and BAM files, but it
> does not work (to my knowledge) with gzip CRAM. Is there any reason for
> this? Is there any plan to implement support to gzip CRAM directly?

CRAM is already superior compression to gzip, so running gzip on the
CRAM will give minimal if any benefit.

James

-- 
James Bonfield (jk...@sa...)
The Sanger Institute, Hinxton, Cambs, CB10 1SA


-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is Wellcome Sanger Institute, Wellcome Genome Campus, 
 Hinxton, CB10 1SA.

[Samtools-help] support to CRAM gzip

[Adrian <ace...@gm...>]

Hi,

I noticed samtools works perfectly fine with gzip SAM and BAM files, but it
does not work (to my knowledge) with gzip CRAM. Is there any reason for
this? Is there any plan to implement support to gzip CRAM directly?

Regards

Adrian

[Samtools-help] Release 1.17

[Robert D. <rm...@sa...>]

Samtools (and HTSlib and BCFtools) version 1.17 is now available from
GitHub and SourceForge.

https://github.com/samtools/htslib/releases/tag/1.17
https://github.com/samtools/samtools/releases/tag/1.17
https://github.com/samtools/bcftools/releases/tag/1.17 
https://sourceforge.net/projects/samtools/

The main changes are listed below:

------------------------------------------------------------------------------
htslib - changes v1.17
------------------------------------------------------------------------------

* A new API for iterating through a BAM record's aux field. (PR#1354,
   addresses #1319.  Thanks to John Marshall)

* Text mode for bgzip. Allows bgzip to compress lines of text with block
   breaks at newlines. (PR#1493, thanks to Mike Lin for the initial version
   PR#1369)

* Make tabix support CSI indices with large positions.  Unlike SAM and VCF
   files, BED files do not set a maximum reference length which hindered CSI
   support.  This change sets an arbitrary large size of 100G to enable it to
   work. (PR#1506)

* Add a fai_line_length function.  Exposes the internal line-wrap length.
   (PR#1516)

* Check for invalid barcode tags in fastq output. (PR#1518, fixes
   samtools#1728.  Reported by Poshi)

* Warn if reference found in a CRAM file is not contained in the specified
   reference file. (PR#1517 and PR#1521, adds diagnostics for #1515. Reported
   by Wei WeiDeng)

* Add a faidx_seq_len64 function that can return sequence lengths longer than
   INT_MAX.  At the same time limit faidx_seq_len to INT_MAX output.  Also add
   a fai_adjust_region to ensure given ranges do not go beyond the end of the
   requested sequence. (PR#1519)

* Add a bcf_strerror function to give text descriptions of BCF errors.
   (PR#1510)

* Add CRAM SQ/M5 header checking when specifying a fasta file.  This is
   to prevent creating a CRAM that cannot be decoded again. (PR#1522.  In
   response to samtools#1748 though not a direct fix)

* Improve support for very long input lines (> 2Gbyte).  This is mostly
   useful for tabix which does not do much interpretation of its input.
   (PR#1542, a partial fix for #1539)

* Speed up load_ref_portion.  This function has been sped up by about 7x,
   which speeds up low-depth CRAM decoding by about 10%. (PR#1551)

* Expand CRAM API to cope with new samtools cram_size command. (PR#1546)

* Merges neighbouring I and D ops into one op within pileup. This means
   4M1D1D1D3M is reported as 4M3D3M.   Fixing this in sam.c means not only is
   samtools mpileup now looking better, but any tool using the mpileup API
   will be getting consistent results. (PR#1552, fixes the last remaining
   part of samtools#139)

* Update the API documentation for bgzf_mt as it refered to a previous
   iteration. (PR#1556, fixes #1553.  Reported by Raghavendra Padmanabhan)

Build changes
-------------

* Use POSIX grep in testing as egrep and fgrep are considered obsolete.
   (PR#1509, thanks to David Seifert)

* Switch to building libdefalte with cmake for Cirris CI. (PR#1511)

* Ensure strings in config_vars.h are escaped correctly. (PR#1530, fixes
   #1527. Reported by Lucas Czech)

* Easier modification of shared library permissions during install. (PR#1532,
   fixes #1525. Reported by StephDC)

* Fix build on ancient compilers.  Added -std=gnu90 to build tests so older C
   compilers will still be happy. (PR#1524, fixes #1523.  Reported by
   Martin Jakt)

* Switch MacOS CI tests to an ARM-based image. (PR#1536)

* Cut down the number of embed_ref=2 tests that get run. (PR#1537)

* Add symbol versions to libhts.so.  This is to aid package developers.
   (PR#1560 addresses #1505, thanks to John Marshall. Reported by
   Stefan Bruens)

* htscodecs now updated to v1.4.0. (PR#1563)

* Cleaned up misleading system error reports in test_bgzf. (PR#1565)

Bug fixes
---------

* VCF. Fix n-squared complexity in sample line with many adjacent tabs
   [fuzz]. (PR#1503)

* Improved bcftools detection and reporting of bgzf decode errors. (PR#1504,
   thanks to Lilian Janin. PR#1529 thanks to Bergur Ragnarsson, fixes #1528.
   PR#1554)

* Prevent crash when the only FASTA entry has no sequence [fuzz]. (PR#1507)

* Fixed typo in sam.h documentation. (PR#1512, thanks to kojix2)

* Fix buffer read-overrun in bam_plp_insertion_mod. (PR#1520)

* Fix hash keys being left behind by bcf_hdr_remove. (PR#1535, fixes #1533.
   Reported by Giulio Genovese in #842)

* Make bcf_hdr_idinfo_exists more robust by checking id value exists.
   (PR#1544, fixes #1538.  Reported by Giulio Genovese)

* CRAM improvements. Fixed crash with multi-threaded CRAM.  Fixed a bug
   in the codec parameter learning for CRAM 3.1 name tokeniser. Fixed Cram
   compression container substitution matrix generation, (PR#1558, PR#1559
   and PR#1562)

------------------------------------------------------------------------------
samtools - changes v1.17
------------------------------------------------------------------------------

New work and changes:

* New samtools reset subcommand.  Removes alignment information.  Alignment
   location, CIGAR, mate mapping and flags are updated. If the alignment was
   in reverse direction, sequence and its quality values are reversed and
   complemented and the reverse flag is reset.  Supplementary and secondary
   alignment data are discarded. (PR#1767, implements #1682. Requested by dkj)

* New samtools cram-size subcommand.  It writes out metrics about a CRAM file
   reporting aggregate sizes per block "Content ID" fields, the data-series
   contained within them, and the compression methods used. (PR#1777)

* Added a --sanitize option to fixmate and view.  This performs some sanity
   checks on the state of SAM record fields, fixing up common mistakes made by
   aligners. (PR#1698)

* Permit 1 thread with samtools view.  All other subcommands already allow
   this and it does provide a modest speed increase. (PR#1755, fixes #1743.
   Reported by Goran Vinterhalter)

* Add CRAM_OPT_REQUIRED_FIELDS option for view -c.  This is a big speed up
   for CRAM (maybe 5-fold), but it depends on which filtering options are
   being used. (PR#1776, fixes #1775. Reported by Chang Y)

* New filtering options in samtools depth.  The new --excl-flags option is a
   synonym for -G, with --incl-flags and --require-flags added to match view
   logic. (PR#1718, fixes #1702. Reported by Dario Beraldi)

* Speed up calmd's slow handling of non-position-sorted data by adding
   caching. This uses more memory but is only activated when needed.
   (PR#1723, fixes #1595. Reported by lxwgcool)

* Improve samtools consensus for platforms with instrument specific
   profiles, considerably helping for data with very different indel error
   models and providing base quality recalibration tables. On PacBio HiFi,
   ONT and  Ultima Genomics consensus qualities are also redistributed
   within homopolymers and the likelihood of nearby indel errors is raised.
   (PR#1721, PR#1733)

* Consensus --mark-ins option.  This permits he consensus output to include a
   markup indicating the next base is an insertion. This is necessary as we
   need a way of outputting both consensus and also how that consensus marries
   up with the reference coordinates. (PR#1746)

* Make faidx/fqidx output line length default to the input line length.
   (PR#1738, fixes #1734. Reported by John Marshall)

* Speed up optical duplicate checking where data has a lot of duplicates
   compared to non-duplicates. (PR#1779, fixes #1771. Reported by Poshi)

* For collate use TMPDIR environment variable, when looking for a temporary
   folder. (PR#1782, based on PR#1178 and fixes #1172.  Reported by
   Martin Pollard)

Bug Fixes:

* Fix stats breakage on long deletions when given a reference. (PR#1712,
   fixes #1707. Reported by John Didion)

* In ampliconclip, stop hard clipping from wrongly removing entire reads.
   (PR#1722, fixes #1717. Reported by Kevin Xu)

* Fix bug in ampliconstats where references mentioned in the input file
   headers but not in the bed file would cause it to complain that the SAM
   headers were inconsistent. (PR#1727, fixes #1650. Reported by jPontix)

* Fixed SEGV in samtools collate when no filename given. (PR#1724)

* Changed the default UMI barcode regex in markdup.  The old regex was too
   restrictive.  This version will at least allow the default read name UMI
   as given in the Illumina example documentation. (PR#1737, fixes #1730.
   Reported by yloemie)

* Fix samtools consensus buffer overrun with MD:Z handling. (PR#1745, fixes
   #1744. Reported by trilisser)

* Fix a buffer read-overflow in mpileup and tview on sequences with seq "*".
   (PR#1747)

* Fix view -X command line parsing that was broken in 1.15. (PR#1772, fixes
   #1720.  Reported by Francisco Rodríguez-Algarra and Miguel Machado)

* Stop samtools view -d from reporting meaningless system errors when tag
   validation fails. (PR#1796)

Documentation:

* Add a description of the samtools tview display layout to the man page.
   Documents . vs , and upper vs lowercase. Adds a -s sample example, and
   documents the -w option. (PR#1765, fixes #1759. Reported by
   Lucas Ferreira da Silva)

* Clarify intention of samtools fasta/q in man page and soft vs hard
   clipping. (PR#1794, fixes #1792. Reported by Ryan Lorig-Roach)

* Minor fix to wording of mpileup --rf usage and man page. (PR#1795, fixes
   #1791. Reported by Luka Pavageau)

Non user-visible changes and build improvements:

* Use POSIX grep in testing as egrep and fgrep are considered obsolete.
   (PR#1726, thanks to David Seifert)

* Switch MacOS CI tests to an ARM-based image. (PR#1770)

------------------------------------------------------------------------------
bcftools - changes v1.17
------------------------------------------------------------------------------

Changes affecting the whole of bcftools, or multiple commands:

* The -i/-e filtering expressions

     - Error checks were added to prevent incorrect use of vector arithmetics.
       For example, when evaluating the sum of two vectors A and B, the
       resulting vector could contain nonsense values when the input vectors
       were not of the same length. The fix introduces the following logic:

         - evaluate to C_i = A_i + B_i when length(A)==B(A) and set
           length(C)=length(A)

         - evaluate to C_i = A_i + B_0 when length(B)=1 and set
           length(C)=length(A)

         - evaluate to C_i = A_0 + B_i when length(A)=1 and set
           length(C)=length(B)

         - throw an error when length(A)!=length(B) AND length(A)!=1 AND
           length(B)!=1

     - Arrays in Number=R tags can be now subscripted by alleles found in
       FORMAT/GT. For example,

  FORMAT/AD[GT] > 10        .. require support of more than 10 reads for
                               each allele
  FORMAT/AD[0:GT] > 10      .. same as above, but in the first sample
  sSUM(FORMAT/AD[GT]) > 20  .. require total sample depth bigger than 20

* The commands `consensus -H` and `+split-vep -H`

     - Drop unnecessary leading space in the first header column and newly
       print `#[1]columnName` instead of the previous `# [1]columnName`
       (#1856)

Changes affecting specific commands:

* bcftools +allele-length

     - Fix overflow for indels longer than 512bp and aggregate alleles equal
       or larger than that in the same bin (#1837)

* bcftools annotate

     - Support sample reordering of annotation file (#1785)

     - Restore lost functionality of the --pair-logic option (#1808)

* bcftools call

     - Fix a bug where too many alleles passed to `-C alleles` via `-T` caused
       memory corruption (#1790)

     - Fix a bug where indels constrained with `-C alleles -T` would sometimes
       be missed (#1706)

* bcftools consensus

     - BREAKING CHANGE: the option `-I, --iupac-codes` newly outputs IUPAC
       codes based on FORMAT/GT of all samples. The `-s, --samples` and `-S,
       --samples-file` options can be used to subset samples. In order to
       ignore samples and consider only the REF and ALT columns (the original
       behavior prior to 1.17), run with `-s -` (#1828)

* bcftools convert

     - Make variantkey conversion work for sites without an ALT allele (#1806)

* bcftool csq

     - Fix a bug where a MNV with multiple consequences (e.g. missense +
       stop_gained) would report only the less severe one (#1810)

     - GFF file parsing was made slightly more flexible, newly ids can be just
       'XXX' rather than, for example, 'gene:XXX'

     - New gff2gff perl script to fix GFF formatting differences

* bcftools +fill-tags

     - More of the available annotations are now added by the `-t all` option

* bcftools +fixref

     - New INFO/FIXREF annotation

     - New -m swap mode

* bcftools +mendelian

     - The +mendelian plugin has been deprecated and replaced with
       +mendelian2. The function of the plugin is the same but the command
       line options and the output format has changed, and for this was
       introduced as a new plugin.

* bcftools mpileup

     - Most of the annotations generated by mpileup are now optional via the
       `-a, --annotate` option and add several new (mostly experimental)
       annotations.

     - New option `--indels-2.0` for an EXPERIMENTAL indel calling model.
       This model aims to address some known deficiencies of the current
       indel calling algorithm, specifically, it uses diploid reference
       consensus sequence. Note that in the current version it has the
       potential to increase sensitivity but at the cost of decreased
       specificity.

     - Make the FS annotation (Fisher exact test strand bias) functional and
       remove it from the default annotations

* bcftools norm

     - New --multi-overlaps option allows to set overlapping alleles either to
       the ref allele (the current default) or to a missing allele (#1764 and
       #1802)

     - Fixed a bug in `-m -` which does not split missing FORMAT values
       correctly and could lead to empty FORMAT fields such as `::` instead
       of the correct `:.:` (#1818)

     - The `--atomize` option previously would not split complex indels such
       as C>GGG. Newly these will be split into two records C>G and C>CGG
       (#1832)

* bcftools query

     - Fix a rare bug where the printing of SAMPLE field with `query` was
       incorrectly suppressed when the `-e` option contained a sample
       expression while the formatting query did not. See #1783 for details.

* bcftools +setGT

     - Add new `--new-gt X` option (#1800)

     - Add new `--target-gt r:FLOAT` option to randomly select a proportion of
       genotypes (#1850)

     - Fix a bug where `-t ./x` mode was advertised as selecting both phased
       and unphased half-missing genotypes, but was in fact selecting only
       unphased genotypes (#1844)

* bcftools +split-vep

     - New options `-g, --gene-list` and `--gene-list-fields` which allow to
       prioritize consequences from a list of genes, or restrict output to the
       listed genes

     - New `-H, --print-header` option to print the header with `-f`

     - Work around a bug in the LOFTEE VEP plugin used to annotate gnomAD
       VCFs. There the LoF_info subfield contains commas which, in
       general, makes it impossible to parse the VEP subfields. The
       +split-vep plugin can now work with such files, replacing the
       offending commas with slash (/) characters. See also
       https://github.com/Ensembl/ensembl-vep/issues/1351

     - Newly the `-c, --columns` option can be omitted when a subfield is used
       in `-i/-e` filtering expression. Note that `-c` may still have to be
       given when it is not possible to infer the type of the subfield. Note
       that this is an experimental feature.

* bcftools stats

     - The per-sample stats (PSC) would not be computed when `-i/-e` filtering
       options and the `-s -` option were given but the expression did not
       include sample columns (1835)

* bcftools +tag2tag

     - Revamp of the plugin to allow wider range of tag conversions,
       specifically all combinations from FORMAT/GL,PL,GP to
       FORMAT/GL,PL,GP,GT

* bcftools +trio-dnm2

     - New `-n, --strictly-novel` option to downplay alleles which violate
       Mendelian inheritance but are not novel

     - Allow to set the `--pn` and `--pns` options separately for SNVs and
       indels and make the indel settings more strict by default

     - Output missing FORMAT/VAF values in non-trio samples, rather than
       random nonsense values

* bcftools +variant-distance

     - New option `-d, --direction` to choose the directionality: forward,
       reverse, nearest (the default) or both (#1829)



-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is 215 Euston Road, London, NW1 2BE. 

[Samtools-help] Problem with RNAME field in bam

[gilsanan <gil...@i3...>]

Hi, 

I'm new using samtools and I'm sorry in advance if my question has an
obvious answer... I did not manage to find any solution online or in the
samtools documentation. 

I've been trying to filter the reads in a BAM file by the reference
name. This field appears in the 3rd column of the header in the bam
file. I checked my BAM file and it seems correct. 

This is an example of one read in my file:

> f5f22aad-bb4a-4ec6-805c-4f290cf5d116 16 MA1#0001 51 60 35S188M1D11M76S * 0 0 ATTGCTAATGATACGGCGACCACCGAGATCTACACGTTCAGAGTTCTACAGTCCGACGATCCATTGTTCACCTCACGCAGACGAGCTCTTAGACAAGTGACAGTGTCCACTCGAGACTGGATATTCTCCAGTTACATACCAATAATAACAAATACTTAACATTGTACAAGAATACCCGATTCTGCGTCCAACCTCATTGAAACCCTCCTTTGGGTGGAATTCTGGGTGCCAAGGAATCCAGTCACCACGATATCTCGTATGCCGTCTTCTGCTTGAGCAATACGTAACTGAACGAAGTGCCGACAAAAGA 4>?ABEEEDC@@@?@>>>>>A@?ABBCAAAABCBCCDBBBDBDCBABA@@ABCCCC?>>>>HDDCDAAAA@@@@@@BCBBAABCCD?==<<?@ABCBC@@??>?@BBICBBABFHFDDDIQHHEEDCCDCBBCBCDDDEGFHHGEDDEGIEEEEEECCC@@AACLEFCBBBBAAABCGIKB@?>>?B@>>>>ADBCCCA?;;;<>ACC?>>?A?>??<=,++(78;:2-,+(''((+-2CBAB?>???FCB?>?@?DEC::::D///,-...-?>GEFDAAAAH6555;50,*'%$#"""##$$####"" NM:i:1 ms:i:392 AS:i:392 nn:i:0 tp:A:P cm:i:31 s1:i:168 s2:i:51 de:f:0.005 rl:i:18

where MA1#0001 is the name of one of my references. 

I was trying to filter the file by reference name with no luck (always
empty outputs). I checked the statistics then using the following line:

> samtools idxstats my_file.bam

The output shows 0 reads mapped to every reference, here a fragment of
the output:

> MA1#0001 300 1099 0
> MA1#0002 300 2325 0
> MA1#0003 300 3757 0
> MA1#0004 300 1686 0
> MA1#0005 300 1856 0
> MA1#0006 300 3051 0
> MA1#0007 300 2965 0

Any idea about why this is happening an how to solve it?  Many thanks!! 

Cheers, 

Eva

[Samtools-help] Samtools installation problem

[Trine R. <tri...@uk...>]

Hi,
I have without much luck tried to install samtools on my pc.
I am using WSL2.

when I type the command "make" during installation I get the following error:

/usr/bin/ld: cannot open output file samtools: Is a directory
collect2: error: ld returned 1 exit status
make: *** [Makefile:145: samtools] Error 1

Does someone have a solution as to what my problem could be?,

Best,
Trine

[Samtools-help] Fwd: Humble request and help_replacing variants that failed quality metrics with "N" with bcftools [EXT]

[Jennifer L. <js...@sa...>]

Good morning Krishna,

I’m afraid I’m not a bcftools developer, so I have forwarded your problem to the samtools-help mailing list, in the hope that somebody there will be able to help you.

Jennifer


Begin forwarded message:

From: Krishna Prasad Acharya <kp...@ca...<mailto:kp...@ca...>>
Subject: Humble request and help_replacing variants that failed quality metrics with "N" with bcftools [EXT]
Date: 2 November 2022 at 04:33:13 GMT
To: "jen...@sa...<mailto:jen...@sa...>" <jen...@sa...<mailto:jen...@sa...>>

Dear Dr. Jennifer,
Good Morning!



I am Krishna Prasad Acharya, an MPhil student at Cambridge Vet School, Cambridge, UK. I am trying to use variant calling using the metrics of DP>=10, %QUAL>=30 and mapping MQ>30 and prepare a whole genome alignment. Everything was working well. However, I had an issue replacing "variants that did not pass filtering" with "N". Could you please help me fix this issue? I am requesting you because you are the developer of bcftools, so I thought you can help me to fix this issue.  Attached is the file containing the code.



I hope you will help me fix this issue.



Thank you.


Regards,
Krishna Prasad Acharya
MPhil Student
Cambridge Vet School
Google Scholar:     https://scholar.google.com/citations?user=N7lg6bAAAAAJ&hl=en [scholar.google.com]<https://urldefense.proofpoint.com/v2/url?u=https-3A__scholar.google.com_citations-3Fuser-3DN7lg6bAAAAAJ-26hl-3Den&d=DwMFAw&c=D7ByGjS34AllFgecYw0iC6Zq7qlm8uclZFI0SqQnqBo&r=KoMRxlOPpFl6Xyyr0RoXf3uy3H2lQZNZP8OyP67Mr_4&m=3SiG0Zle08hFl3VmM7TerK5PaHKntciWlWeZlj5SpJYPpXxaUbMNa5vnHpVWVbnX&s=9bSdTvDBps8yxIPQyBfmsWUMvsMFeg_BUmhQX0gsrjY&e=>
Researchgate:        https://www.researchgate.net/profile/Krishna-Prasad-Acharya-2 [researchgate.net]<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.researchgate.net_profile_Krishna-2DPrasad-2DAcharya-2D2&d=DwMFAw&c=D7ByGjS34AllFgecYw0iC6Zq7qlm8uclZFI0SqQnqBo&r=KoMRxlOPpFl6Xyyr0RoXf3uy3H2lQZNZP8OyP67Mr_4&m=3SiG0Zle08hFl3VmM7TerK5PaHKntciWlWeZlj5SpJYPpXxaUbMNa5vnHpVWVbnX&s=EZvrKzjn09xsHppj_UhftYTnFm4ROC4r-_2kqg-uOvE&e=>
Scopus                    https://www.scopus.com/authid/detail.uri?authorId=57226349259 [scopus.com]<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.scopus.com_authid_detail.uri-3FauthorId-3D57226349259&d=DwMFAw&c=D7ByGjS34AllFgecYw0iC6Zq7qlm8uclZFI0SqQnqBo&r=KoMRxlOPpFl6Xyyr0RoXf3uy3H2lQZNZP8OyP67Mr_4&m=3SiG0Zle08hFl3VmM7TerK5PaHKntciWlWeZlj5SpJYPpXxaUbMNa5vnHpVWVbnX&s=GXpodFwdiY4rLWbHLsHW6K2YpFbLOWk32TybD106t4M&e=>
ORCID:                 https://orcid.org/0000-0002-4208-7592 [orcid.org]<https://urldefense.proofpoint.com/v2/url?u=https-3A__orcid.org_0000-2D0002-2D4208-2D7592&d=DwMFAw&c=D7ByGjS34AllFgecYw0iC6Zq7qlm8uclZFI0SqQnqBo&r=KoMRxlOPpFl6Xyyr0RoXf3uy3H2lQZNZP8OyP67Mr_4&m=3SiG0Zle08hFl3VmM7TerK5PaHKntciWlWeZlj5SpJYPpXxaUbMNa5vnHpVWVbnX&s=m0jeaYXSPt9j7u6RTOhfl3TS1y8Vx1P62If4OG34XJI&e=>
Publons:                https://publons.com/researcher/3703799/krishna-prasad-acharya/publications/ [publons.com]<https://urldefense.proofpoint.com/v2/url?u=https-3A__publons.com_researcher_3703799_krishna-2Dprasad-2Dacharya_publications_&d=DwMFAw&c=D7ByGjS34AllFgecYw0iC6Zq7qlm8uclZFI0SqQnqBo&r=KoMRxlOPpFl6Xyyr0RoXf3uy3H2lQZNZP8OyP67Mr_4&m=3SiG0Zle08hFl3VmM7TerK5PaHKntciWlWeZlj5SpJYPpXxaUbMNa5vnHpVWVbnX&s=UmLDrWTCkujcZ-o2jT3PkZUxWGQoa9jsh2LYVGgSeCo&e=>



"Imagination is more powerful than knowledge"
                                      Albert Einstein







-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is 215 Euston Road, London, NW1 2BE.

Re: [Samtools-help] Question about samtools mpileup and overlapping reads [EXT]

[James B. <jk...@sa...>]

On Tue, Oct 25, 2022 at 04:25:32PM +0200, Carla Mariner Llicer wrote:
> I'm using Samtools version 1.3.1. I was analysing some illumina paired-end

Please use a modern version.  1.3 is 6 year old, with many bug fixes
and improvements made since then.  This particular issue was fixed in
1.13.

However, also see wish-list item
https://github.com/samtools/samtools/issues/1354 for potential issues
still associated with random selection vs by-quality.

James

-- 
James Bonfield (jk...@sa...)
The Sanger Institute, Hinxton, Cambs, CB10 1SA


-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is 215 Euston Road, London, NW1 2BE. 

[Samtools-help] Question about samtools mpileup and overlapping reads

[Carla M. L. <cma...@gm...>]

Hi!
I'm using Samtools version 1.3.1. I was analysing some illumina paired-end
sequenced samples and I realised that when running mpileup with this
version, by default, samtools discards all the Minus Strand reads when
there are overlapping reads. Is there any reason why always discards Minus
Strand reads? Does it make sense?  What happens is that I have some real
SNPs in regions where all the reads overlap so, when it discards all the
Minus Strand reads, the SNPs appear only in one strand and then it doesn't
pass the downstream filters that we apply. I know that with the
option --ignore-overlaps I can keep all the overlapping reads so my
question is Why samtools does not discard Positive strands or Minus Strands
in a random way by default?

Thanks for your attention,

Carla

Re: [Samtools-help] consensus not recognized [EXT]

[James B. <jk...@sa...>]

On Wed, Sep 07, 2022 at 09:27:17PM +0000, Leal, Joel via Samtools-help wrote:
> Thanks for your answer. Does it mean I need to install v1.16.1 instead of v1.2?

Yes.

See https://github.com/samtools/samtools/releases for the list of
what's changed with each release.  A quick search for consensus in
that page shows it appeared with 1.15.

Wherever you read that it was in 1.2 was simply wrong.

James

-- 
James Bonfield (jk...@sa...)
The Sanger Institute, Hinxton, Cambs, CB10 1SA


-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is 215 Euston Road, London, NW1 2BE. 

Re: [Samtools-help] consensus not recognized

[Tom B. <tb...@um...>]

Joel  -

I strongly suspect that the problem is "Which is the higher version number,
samtools v1.16.1 or samtools v1.2 ?"  In this case it is v1.16.1.  Our
sysadmins have had trouble with this too and have "updated" us back to
version 1.2.

-  thanks  -  tom blackwell  -

On Wed, Sep 7, 2022 at 4:40 PM Leal, Joel via Samtools-help <
sam...@li...> wrote:

> Hi,
>
> I am trying to generate a consensus sequence from a .BAM file using
> Samtools. I found it was developed for Samtools1.2. However, when calling
> the function consensus, it is not recognized. The function mpileup is
> recognized.
>
> Could you tell me what is the issue, please?
>
> Thanks,
> Joel
>
>
> _______________________________________________
> Samtools-help mailing list
> Sam...@li...
> https://lists.sourceforge.net/lists/listinfo/samtools-help
>

Re: [Samtools-help] consensus not recognized

[Leal, J. <jd...@he...>]

Thanks for your answer. Does it mean I need to install v1.16.1 instead of v1.2?

Joel
________________________________
From: Tom Blackwell <tb...@um...>
Sent: Wednesday, September 7, 2022 2:01 PM
To: Leal, Joel <jd...@he...>
Cc: sam...@li... <sam...@li...>
Subject: Re: [Samtools-help] consensus not recognized

Joel  -

I strongly suspect that the problem is "Which is the higher version number, samtools v1.16.1 or samtools v1.2 ?"  In this case it is v1.16.1.  Our sysadmins have had trouble with this too and have "updated" us back to version 1.2.

-  thanks  -  tom blackwell  -

On Wed, Sep 7, 2022 at 4:40 PM Leal, Joel via Samtools-help <sam...@li...<mailto:sam...@li...>> wrote:
Hi,

I am trying to generate a consensus sequence from a .BAM file using Samtools. I found it was developed for Samtools1.2. However, when calling the function consensus, it is not recognized. The function mpileup is recognized.

Could you tell me what is the issue, please?

Thanks,
Joel

[cid:18319bef83bcb971f161]
_______________________________________________
Samtools-help mailing list
Sam...@li...<mailto:Sam...@li...>
https://lists.sourceforge.net/lists/listinfo/samtools-help<https://urldefense.com/v3/__https://lists.sourceforge.net/lists/listinfo/samtools-help__;!!LLK065n_VXAQ!kj4Uu5XaZVIF-MsGG3-Os_mgJQ7N-eLd0ETcyVrgDMB-rnlkpvf829NXdb_16i7peqfKgOr5GhwbhkHrd_g$>

[Samtools-help] consensus not recognized

[Leal, J. <jd...@he...>]

Hi,

I am trying to generate a consensus sequence from a .BAM file using Samtools. I found it was developed for Samtools1.2. However, when calling the function consensus, it is not recognized. The function mpileup is recognized.

Could you tell me what is the issue, please?

Thanks,
Joel

[cid:ac963de6-540b-40c8-a8c8-0cfd7d336023]

[Samtools-help] SAMtools release 1.16.1

[Robert D. <rm...@sa...>]

Samtools version 1.16.1 is now available from GitHub and SourceForge.
This release fixes some bugs in the new template-coordinate sort feature.

https://sourceforge.net/projects/samtools/
https://github.com/samtools/samtools/releases/tag/1.16.1

The main changes are listed below:

------------------------------------------------------------------------------
samtools - changes v1.16.1
------------------------------------------------------------------------------

Bug fixes:

  * Fixed a bug with the template-coordinate sort which caused incorrect
    ordering when using threads, or processing large files that don't fit
    completely in memory. (PR#1703, thanks to Nils Homer)

  * Fixed a crash that occurred when trying to use `samtools merge` in
    template-coordinate mode. (PR#1705, thanks to Nils Homer)



-- 
 The Wellcome Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is 215 Euston Road, London, NW1 2BE.