From 25e366d523d3dc51f8da4fc6c48e88b726cf7209 Mon Sep 17 00:00:00 2001 From: =?utf8?q?=C3=89tienne=20Mollier?= Date: Fri, 14 Nov 2025 12:32:40 +0100 Subject: [PATCH] New upstream version 0.23.3+ds --- .cirrus.yml | 82 ------- Containerfile | 13 - INSTALL | 2 +- NEWS | 110 ++++++++- README.rst | 4 +- bcftools/bcftools.pysam.h | 6 + devtools/buildwheels.sh | 73 ------ devtools/conda-recipe/build.sh | 8 - devtools/conda-recipe/meta.yaml | 29 --- devtools/emulate-tools.py | 59 +++++ devtools/install-CGAT-tools.sh | 282 ---------------------- devtools/install-prerequisites.sh | 13 + devtools/run_tests_travis.sh | 126 ---------- doc/installation.rst | 2 +- import/pysam.h | 6 + pyproject.toml | 6 +- pysam/conftest_cstd.c | 11 + pysam/libcalignedsegment.pxd | 21 +- pysam/libcalignedsegment.pyi | 28 ++- pysam/libcalignedsegment.pyx | 389 ++++++++++++++++++------------ pysam/libcalignmentfile.pyx | 51 ++-- pysam/libcbcf.pyx | 56 +++-- pysam/libchtslib.pxd | 6 +- pysam/libctabix.pyx | 18 +- pysam/libcutils.pxd | 2 + pysam/libcutils.pyx | 15 +- pysam/version.py | 2 +- samtools/samtools.pysam.h | 6 + setup.py | 89 ++++--- tests/AlignedSegment_test.py | 130 ++++++++++ tests/compile_test.py | 21 +- tests/refactoring.pl | 60 ----- tests/samtools_test.py | 10 + 33 files changed, 784 insertions(+), 952 deletions(-) delete mode 100644 .cirrus.yml delete mode 100644 Containerfile delete mode 100755 devtools/buildwheels.sh delete mode 100644 devtools/conda-recipe/build.sh delete mode 100644 devtools/conda-recipe/meta.yaml create mode 100755 devtools/emulate-tools.py delete mode 100755 devtools/install-CGAT-tools.sh delete mode 100755 devtools/run_tests_travis.sh create mode 100644 pysam/conftest_cstd.c delete mode 100644 tests/refactoring.pl diff --git a/.cirrus.yml b/.cirrus.yml deleted file mode 100644 index 1ad3445..0000000 --- a/.cirrus.yml +++ /dev/null @@ -1,82 +0,0 @@ -build_wheels_task: - only_if: $CIRRUS_BRANCH =~ "release/.*" || $CIRRUS_TAG =~ "v0\..*" - - matrix: - - compute_engine_instance: - image_project: cirrus-images - image: family/docker-builder-arm64 - architecture: arm64 - platform: linux - matrix: - - name: Build ARM Linux py3.6-9 wheels - env: - CIBW_BUILD: "cp36-* cp37-* cp38-* cp39-*" - - name: Build ARM Linux py3.10-13 wheels - env: - CIBW_BUILD: "cp310-* cp311-* cp312-* cp313-*" - - - name: Build ARM macOS wheels - macos_instance: - image: ghcr.io/cirruslabs/macos-sonoma-base:latest - env: - CIBW_BUILD: "cp39-* cp310-* cp311-* cp312-* cp313-*" - - alias: build_wheels - - env: - CIRRUS_CLONE_DEPTH: 1 - - VENV: $HOME/relenv - PATH: $VENV/bin:$PATH - - CIBW_SKIP: "*-musllinux_*" - CIBW_BUILD_VERBOSITY: 1 - - # Avoid linking with non-system library libdeflate.dylib - CIBW_ENVIRONMENT_MACOS: HTSLIB_CONFIGURE_OPTIONS="--without-libdeflate" - - CIBW_MANYLINUX_AARCH64_IMAGE: manylinux_2_28 - - install_script: | - python3 -m venv $VENV - pip3 install cibuildwheel==2.22.0 - - build_script: | - cibuildwheel - - wheels_artifacts: - path: wheelhouse/*.whl - -upload_pypi_task: - only_if: $CIRRUS_BRANCH =~ "release/.*" || $CIRRUS_TAG =~ "v0\..*" - depends_on: build_wheels - - name: Publish ARM wheels - - container: - image: python:latest - - env: - CIRRUS_CLONE_DEPTH: 1 - API_BASEURL: https://api.cirrus-ci.com/v1 - TWINE_USERNAME: __token__ - - install_script: | - python3 -m pip install twine - - get_artifacts_script: | - curl -sSLO $API_BASEURL/artifact/build/$CIRRUS_BUILD_ID/wheels.zip - unzip -q wheels.zip - - upload_script: | - case "$CIRRUS_TAG" in - v0.*) - export TWINE_REPOSITORY=pypi TWINE_PASSWORD=$PYPI_TOKEN ;; - *) - export TWINE_REPOSITORY=testpypi TWINE_PASSWORD=$TESTPYPI_TOKEN ;; - esac - - echo Uploading wheels to $TWINE_REPOSITORY... - - python3 -m twine check wheelhouse/*.whl - python3 -m twine upload --disable-progress-bar wheelhouse/*.whl diff --git a/Containerfile b/Containerfile deleted file mode 100644 index 633f07e..0000000 --- a/Containerfile +++ /dev/null @@ -1,13 +0,0 @@ -FROM ubi8:latest - -RUN yum update \ - && yum install -y python3-pip python3-devel pigz \ - && cd /usr/local/bin \ - && ln -s /usr/bin/python3 python \ - && pip3 --no-cache-dir install --upgrade pip \ - && yum clean all \ - && echo "system packages installed" - -RUN python -m pip install pysam - -WORKDIR /opt/ diff --git a/INSTALL b/INSTALL index 6dd0163..1f3677a 100644 --- a/INSTALL +++ b/INSTALL @@ -13,8 +13,8 @@ Conda installation To install pysam in your current conda_ environment, type:: - conda config --add channels r conda config --add channels bioconda + conda config --add channels conda-forge conda install pysam This will install pysam from the bioconda_ channel and automatically diff --git a/NEWS b/NEWS index a952ef4..17c916c 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,106 @@ .. An online version of the release history can be found here: .. http://pysam.readthedocs.io/en/latest/release.html +Release 0.23.3 +============== + +.. rubric:: 10 June 2025 + +This is a bugfix release, still wrapping htslib/samtools/bcftools 1.21. + +It has been tested with Python versions 3.8 through 3.13, and wheels are +available via PyPI_ for all of those Python versions. Wheels are built for +macOS and Linux (manylinux_2_28 and musllinux_1_2) on both ARM and x86-64. + +The final pysam release that supported Python 3.6 and 3.7 was v0.23.0. + +Bugs fixed: + +* Worked around Cython 3.1.2 bug whereby a new Cython-internal function + incorrectly triggered the pysam build system's duplicate symbol check. + + +Release 0.23.2 +============== + +.. rubric:: 6 June 2025 + +This is a bugfix release, still wrapping htslib/samtools/bcftools 1.21. + +It has been tested with Python versions 3.8 through 3.13, and wheels are +available via PyPI_ for all of those Python versions. Wheels are built for +macOS and Linux (manylinux_2_28 and musllinux_1_2) on both ARM and x86-64. + +Bugs fixed: + +* Pysam 0.23.1 inadvertently broke binary compatibility for Cython projects + by changing the size of :class:`.AlignedSegment`. That class has been + restored to its previous size, so v0.23.2 restores binary compatibility + with Cython projects previously compiled against v0.23.0 and earlier. + (This bug does not affect pure Python projects using pysam.) + +* Improved I/O exception messages in :class:`.AlignmentFile` and + :class:`.VariantFile`, and ensured that error codes returned by all + HTSlib function invocations result in an exception being raised. + + +Release 0.23.1 +============== + +.. rubric:: 28 May 2025 + +This is a bugfix release, still wrapping htslib/samtools/bcftools 1.21. + +It has been tested with Python versions 3.8 through 3.13, and wheels are +available via PyPI_ for all of those Python versions. Wheels are built for +macOS and Linux (manylinux_2_28 and musllinux_1_2) on both ARM and x86-64. + + This release inadvertently broke binary compatibility for Cython projects. + This should not affect Python projects using pysam. However Cython projects + should update to release 0.23.2. + +Bugs fixed: + +* Worked around Cython 3.1.0 behaviour change so that ``pysam.CMATCH`` etc + remain as synonyms for ``pysam.CIGAR_OPS.CMATCH`` etc. Note that a future + release will remove these synonyms to comply with Cython's new convention. + (#1339, PR #1340, reported by Adam Taranto, Andy Bond, et al) + +* The Python interpreter no longer exits abruptly on x86-64 Linux + when ``pysam.bcftools.*()`` commands terminate with error messages. + (#1333, #1335, reported by Maximilian Blacher and Ben Habermeyer) + +* Fixed :meth:`.AlignedSegment.get_aligned_pairs` type hint, which now + has overloaded hints reflecting the different tuples returned. + (PR #1342, thanks to Gosuke Shibahara) + +* Corrected "invalid type for record" exception message. + (PR #1334, thanks to Gaik Tamazian) + +* Fixed :class:`.TabixIterator` error handling. (#1328, reported by @Fan-iX) + +* Fixed bugs in the caching of :attr:`.AlignedSegment.query_sequence` and + :attr:`.AlignedSegment.query_qualities` values. + +* Pysam once again builds on CentOS 7, as it now works around limitations + in that obsolete platform's obsolete GCC 4.8.x system compiler. + (#1327, reported by Mario Fasold) + +New functionality: + +* New :attr:`.AlignedSegment.query_qualities_str` property enables + accessing the QUAL field as an ASCII-encoded base quality string. + The existing :attr:`.AlignedSegment.query_qualities` property can + now take an array or such a string when setting the value. + (PR #1324, PR #1341, requested by Nils Homer) + +* Improved :attr:`.AlignedSegment.cigarstring` performance. + (PR #1295, thanks to @limwz01) + +* Wheels are now built for musllinux in addition to manylinux and macOS. + (Requested by Nils Homer) + + Release 0.23.0 ============== @@ -696,11 +796,11 @@ Overview The 0.9.0 release upgrades htslib to htslib 1.3 and numerous other enhancements and bugfixes. See below for a detailed list. -`Htslib 1.3 `_ -comes with additional capabilities for remote file access which depend -on the presence of optional system libraries. As a consequence, the -installation script :file:`setup.py` has become more complex. For an -overview, see :ref:`installation`. We have tested installation on +The `Htslib 1.3 `_ +release comes with additional capabilities for remote file access which +depend on the presence of optional system libraries. As a consequence, +the installation script :file:`setup.py` has become more complex. For +an overview, see :ref:`installation`. We have tested installation on linux and OS X, but could not capture all variations. It is possible that a 0.9.1 release might follow soon addressing installation issues. diff --git a/README.rst b/README.rst index 89a51f3..68f2fd5 100644 --- a/README.rst +++ b/README.rst @@ -15,9 +15,9 @@ includes an interface for tabix_. If you are using the conda packaging manager (e.g. miniconda or anaconda), you can install pysam from the `bioconda channel `_:: - conda config --add channels defaults - conda config --add channels conda-forge conda config --add channels bioconda + conda config --add channels conda-forge + conda config --set channel_priority strict conda install pysam Installation through bioconda is the recommended way to install pysam diff --git a/bcftools/bcftools.pysam.h b/bcftools/bcftools.pysam.h index e6717bb..a4e2a42 100644 --- a/bcftools/bcftools.pysam.h +++ b/bcftools/bcftools.pysam.h @@ -69,6 +69,12 @@ extern int bcftools_main(int argc, char *argv[]); #define bam_smpl_destroy bcftools_bam_smpl_destroy #define read_file_list bcftools_read_file_list +/*! A non-static error() function name is used in bcftools, which collides + with glibc's error() function and leads to the wrong function being called + on some platforms. #define this name with a prefix to avoid this collision. + */ +#define error bcftools_error + #endif #endif diff --git a/devtools/buildwheels.sh b/devtools/buildwheels.sh deleted file mode 100755 index 4d12566..0000000 --- a/devtools/buildwheels.sh +++ /dev/null @@ -1,73 +0,0 @@ -#!/bin/bash -# -# Build manylinux1 wheels for pysam. Based on the example at -# -# -# It is best to run this in a fresh clone of the repository! -# -# Before running, make sure to update image: -# -# docker pull quay.io/pypa/manylinux1_x86_64 -# -# Run this within the repository root: -# docker run --rm -v $(pwd):/io quay.io/pypa/manylinux1_x86_64 /io/buildwheels.sh -# -# The wheels will be put into the wheelhouse/ subdirectory. -# -# For interactive tests: -# docker run -it -v $(pwd):/io quay.io/pypa/manylinux1_x86_64 /bin/bash - -set -xeuo pipefail - -# For convenience, if this script is called from outside of a docker container, -# it starts a container and runs itself inside of it. -if ! grep -q docker /proc/1/cgroup; then - # We are not inside a container - exec docker run --rm -v $(pwd):/io quay.io/pypa/manylinux1_x86_64 /io/$0 -fi - -yum install -y zlib-devel bzip2-devel xz-devel - -# Without libcurl support, htslib can open files from HTTP and FTP URLs. -# With libcurl support, it also supports HTTPS and S3 URLs, but libcurl needs a -# current version of OpenSSL, and we do not want to be responsible for -# updating the wheels as soon as there are any security issues. So disable -# libcurl for now. -# See also . -# -export HTSLIB_CONFIGURE_OPTIONS="--disable-libcurl" - -PYBINS="/opt/python/*/bin" -echo $PYBINS - -for PYBIN in ${PYBINS}; do - ${PYBIN}/pip install -r /io/requirements.txt - ${PYBIN}/pip wheel /io/ -w wheelhouse/ -done - -# Bundle external shared libraries into the wheels -# -# The '-L ""' option is a workaround. By default, auditwheel puts all external -# libraries (.so files) into a .libs directory and sets the RUNPATH to $ORIGIN/.libs. -# When HTSLIB_MODE is 'shared' (now the default), then all so libraries part of -# pysam require that RUNPATH is set to $ORIGIN (without the .libs). It seems -# auditwheel overwrites $ORIGIN with $ORIGIN/.libs. This workaround makes -# auditwheel keeps the RUNPATH at "$ORIGIN" and put all the external libraries into -# the pysam directory. -for whl in wheelhouse/*.whl; do - auditwheel repair -L "" $whl -w /io/wheelhouse/ -done - -# mkdir -p /io/wheelhouse -# cp wheelhouse/*.whl /io/wheelhouse - -# Created files are owned by root, so fix permissions. -chown -R --reference=/io/setup.py /io/wheelhouse/ - -# TODO Install packages and test them -for PYBIN in ${PYBINS}; do - ${PYBIN}/pip install pysam --no-index -f /io/wheelhouse - # smoketest - (cd $HOME; ${PYBIN}/python -c 'import pysam') - # todo: add more tests -done diff --git a/devtools/conda-recipe/build.sh b/devtools/conda-recipe/build.sh deleted file mode 100644 index 32b67db..0000000 --- a/devtools/conda-recipe/build.sh +++ /dev/null @@ -1,8 +0,0 @@ -#!/bin/bash - -# Use internal htslib -chmod a+x ./htslib/configure -export CFLAGS="-I${PREFIX}/include/curl/ -I${PREFIX}/include -L${PREFIX}/lib" -export HTSLIB_CONFIGURE_OPTIONS="--disable-libcurl" - -$PYTHON setup.py install diff --git a/devtools/conda-recipe/meta.yaml b/devtools/conda-recipe/meta.yaml deleted file mode 100644 index 4e57895..0000000 --- a/devtools/conda-recipe/meta.yaml +++ /dev/null @@ -1,29 +0,0 @@ -package: - name: pysam - version: 0.8.5 - -source: - path: ../../ - -build: - number: 0 - -requirements: - build: - - python - - setuptools - - zlib - - cython - - run: - - python - - zlib - -test: - imports: - - pysam - -about: - home: https://github.com/pysam-developers/pysam - license: MIT - summary: Pysam is a python module for reading and manipulating Samfiles. It's a lightweight wrapper of the samtools C-API. Pysam also includes an interface for tabix. diff --git a/devtools/emulate-tools.py b/devtools/emulate-tools.py new file mode 100755 index 0000000..30dd75a --- /dev/null +++ b/devtools/emulate-tools.py @@ -0,0 +1,59 @@ +#!/usr/bin/env python3 +""" +This script can be symlinked to samtools, bcftools, bgzip, or tabix. +When invoked under one of those names, it will emulate that tool's +behaviour by using pysam's facilities. +""" + +import argparse +import gzip +import os +import sys +import tempfile + +import pysam + +command = os.path.basename(sys.argv[0]) + +if command in ("samtools", "bcftools"): + if len(sys.argv) > 1: + try: + tool = pysam.utils.PysamDispatcher(command, sys.argv[1]) + tool(*sys.argv[2:], catch_stdout=None) + print(tool.stderr, end="", file=sys.stderr) + except pysam.utils.SamtoolsError as e: + sys.exit(f"emulate-tools.py: {e}") + + else: + version = getattr(pysam.version, f"__{command}_version__") + print(f"Program: {command}\nVersion: {version}", file=sys.stderr) + +else: + parser = argparse.ArgumentParser() + parser.add_argument("-c", "--stdout", action="store_true") + parser.add_argument("-d", "--decompress", action="store_true") + parser.add_argument("-f", "--force", action="store_true") + parser.add_argument("-p", "--preset") + parser.add_argument("input_file", nargs="?") + opt = parser.parse_args() + + if command == "bgzip": + if opt.decompress: + with gzip.open(sys.stdin.buffer, "rb") as f: + sys.stdout.buffer.write(f.read()) + + elif opt.input_file and opt.stdout: + pysam.tabix_compress(opt.input_file, "-", force=True) + + else: + f = tempfile.NamedTemporaryFile(delete=False) + f.write(sys.stdin.buffer.read()) + f.close() + pysam.tabix_compress(f.name, "-", force=True) + os.remove(f.name) + + elif command == "tabix": + pysam.tabix_index(opt.input_file, preset=opt.preset, force=opt.force) + + else: + sys.exit(f"emulate-tools.py: unknown command {command!r}") diff --git a/devtools/install-CGAT-tools.sh b/devtools/install-CGAT-tools.sh deleted file mode 100755 index e45d391..0000000 --- a/devtools/install-CGAT-tools.sh +++ /dev/null @@ -1,282 +0,0 @@ -#!/usr/bin/env bash - -# function to detect the Operating System -detect_os(){ - -if [ -f /etc/os-release ]; then - - OS=$(cat /etc/os-release | awk '/VERSION_ID/ {sub("="," "); print $2;}' | sed 's/\"//g' | awk '{sub("\\."," "); print $1;}') - if [ "$OS" != "12" ] ; then - - echo - echo " Ubuntu version not supported " - echo - echo " Only Ubuntu 12.x has been tested so far " - echo - exit 1; - - fi - - OS="ubuntu" - -elif [ -f /etc/system-release ]; then - - OS=$(cat /etc/system-release | awk ' {print $4;}' | awk '{sub("\\."," "); print $1;}') - if [ "$OS" != "6" ] ; then - echo - echo " Scientific Linux version not supported " - echo - echo " Only 6.x Scientific Linux has been tested so far " - echo - exit 1; - fi - - OS="sl" - -else - - echo - echo " Operating system not supported " - echo - echo " Exiting installation " - echo - exit 1; - -fi -} # detect_os - -# message to display when the OS is not correct -sanity_check_os() { - echo - echo " Unsupported operating system: $OS " - echo " Installation aborted " - echo - exit 1; -} # sanity_check_os - -# function to install operating system dependencies -install_os_packages() { - -if [ "$OS" == "ubuntu" -o "$OS" == "travis" ] ; then - - echo - echo " Installing packages for Ubuntu " - echo - - apt-get install -y gcc g++ - -elif [ "$OS" == "sl" ] ; then - - echo - echo " Installing packages for Scientific Linux " - echo - - yum -y install gcc zlib-devel gcc-c++ - -else - - sanity_check_os - -fi # if-OS -} # install_os_packages - -# function to install Python dependencies -install_python_deps() { - -if [ "$OS" == "ubuntu" -o "$OS" == "sl" ] ; then - - echo - echo " Installing Python dependencies for $1 " - echo - - # Create virtual environment - cd - mkdir CGAT - cd CGAT - wget --no-check-certificate https://pypi.python.org/packages/source/v/virtualenv/virtualenv-1.10.1.tar.gz - tar xvfz virtualenv-1.10.1.tar.gz - rm virtualenv-1.10.1.tar.gz - cd virtualenv-1.10.1 - python virtualenv.py cgat-venv - source cgat-venv/bin/activate - - # Install Python prerequisites - pip install cython - -elif [ "$OS" == "travis" ] ; then - # Travis-CI provides a virtualenv with Python 2.7 - echo - echo " Installing Python dependencies in travis " - echo - - # Install Python prerequisites - pip install cython - pip install nose - -else - - sanity_check_os - -fi # if-OS -} # install_python_deps - -# common set of tasks to prepare external dependencies -nosetests_external_deps() { -echo -echo " Running nosetests for $1 " -echo - -pushd . - -# create a new folder to store external tools -mkdir -p $HOME/CGAT/external-tools - -# install samtools -cd $HOME/CGAT/external-tools -curl -L http://downloads.sourceforge.net/project/samtools/samtools/1.3/samtools-1.3.tar.bz2 > samtools-1.3.tar.bz2 -tar xjf samtools-1.3.tar.bz2 -cd samtools-1.3 -make -PATH=$PATH:$HOME/CGAT/external-tools/samtools-1.3 - -echo "installed samtools" -samtools --version - -if [ $? != 0 ]; then - exit 1 -fi - -# install bcftools -cd $HOME/CGAT/external-tools -curl -L https://github.com/samtools/bcftools/releases/download/1.3/bcftools-1.3.tar.bz2 > bcftools-1.3.tar.bz2 -tar xjf bcftools-1.3.tar.bz2 -cd bcftools-1.3 -make -PATH=$PATH:$HOME/CGAT/external-tools/bcftools-1.3 - -echo "installed bcftools" -bcftools --version - -if [ $? != 0 ]; then - exit 1 -fi - -popd - -} # nosetests_external_deps - - -# function to run nosetests -run_nosetests() { - -echo -echo " Running nosetests for $1 " -echo - -# prepare external dependencies -nosetests_external_deps $OS - -# install code -python setup.py install - -# change into tests directory. Otherwise, -# 'import pysam' will import the repository, -# not the installed version. This causes -# problems in the compilation test. -cd tests - -# create auxiliary data -echo -echo 'building test data' -echo -make -C pysam_data all -make -C cbcf_data all -make -C tabix_data all - -# run nosetests -# -s: do not capture stdout, conflicts with pysam.dispatch -# -v: verbose output -nosetests -s -v - -} # run_nosetests - -# function to display help message -help_message() { -echo -echo " Use this script as follows: " -echo -echo " 1) Become root and install the operating system* packages: " -echo " ./install-CGAT-tools.sh --install-os-packages" -echo -echo " 2) Now, as a normal user (non root), install the Python dependencies**: " -echo " ./install-CGAT-tools.sh --install-python-deps" -echo -echo " At this stage the CGAT Code Collection is ready to go and you do not need further steps. Please type the following for more information:" -echo " source $HOME/CGAT/virtualenv-1.10.1/cgat-venv/bin/activate" -echo " cgat --help " -echo -echo " The CGAT Code Collection tests the software with nosetests. If you are interested in running those, please continue with the following steps:" -echo -echo " 3) Become root to install external tools and set up the environment: " -echo " ./install-CGAT-tools.sh --install-nosetests-deps" -echo -echo " 4) Then, back as a normal user (non root), run nosetests as follows:" -echo " ./install-CGAT-tools.sh --run-nosetests" -echo -echo " NOTES: " -echo " * Supported operating systems: Ubuntu 12.x and Scientific Linux 6.x " -echo " ** An isolated virtual environment will be created to install Python dependencies " -echo -exit 1; -} - - -# the main script starts here - -if [ $# -eq 0 -o $# -gt 1 ] ; then - - help_message - -else - - if [ "$1" == "--help" ] ; then - - help_message - - elif [ "$1" == "--travis" ] ; then - - OS="travis" - install_os_packages - install_python_deps - run_nosetests - - elif [ "$1" == "--install-os-packages" ] ; then - - detect_os - install_os_packages - - elif [ "$1" == "--install-python-deps" ] ; then - - detect_os - install_python_deps - - elif [ "$1" == "--install-nosetests-deps" ] ; then - - detect_os - install_nosetests_deps - - elif [ "$1" == "--run-nosetests" ] ; then - - detect_os - run_nosetests - - else - - echo - echo " Incorrect input parameter: $1 " - help_message - - fi # if argument 1 - -fi # if number of input parameters - diff --git a/devtools/install-prerequisites.sh b/devtools/install-prerequisites.sh index 83be6aa..e78ef59 100755 --- a/devtools/install-prerequisites.sh +++ b/devtools/install-prerequisites.sh @@ -12,6 +12,7 @@ elif test -x /usr/bin/yum; then else echo Installing non-test prerequisites via yum... yum -y install zlib-devel bzip2-devel xz-devel curl-devel openssl-devel + emulate=yes fi elif test -d /etc/dpkg; then @@ -23,6 +24,7 @@ elif test -x /sbin/apk; then echo Installing non-test prerequisites via apk... apk update apk add zlib-dev bzip2-dev xz-dev curl-dev openssl-dev + emulate=yes elif test -x ${HOMEBREW_PREFIX-/usr/local}/bin/brew; then echo Installing prerequisites via brew... @@ -32,3 +34,14 @@ elif test -x ${HOMEBREW_PREFIX-/usr/local}/bin/brew; then else echo No package manager detected fi + +if test -n "$emulate" && test $# -ge 2; then + emulator=$1 + bindir=$2 + echo Creating symlinks to $emulator in $bindir... + mkdir -p $bindir + ln -s $emulator $bindir/samtools + ln -s $emulator $bindir/bcftools + ln -s $emulator $bindir/bgzip + ln -s $emulator $bindir/tabix +fi diff --git a/devtools/run_tests_travis.sh b/devtools/run_tests_travis.sh deleted file mode 100755 index 1f14fc3..0000000 --- a/devtools/run_tests_travis.sh +++ /dev/null @@ -1,126 +0,0 @@ -#!/usr/bin/env bash - -# test script for pysam. -# The script performs the following tasks: -# 1. Setup a conda environment and install dependencies via conda -# 2. Build pysam via the conda recipe -# 3. Build pysam via setup.py from repository -# 4. Run tests on the setup.py version -# 5. Additional build tests -# 5.1 pip install with cython -# 5.2 pip install without cython -# 5.3 pip install without cython and without configure options - -pushd . - -WORKDIR=`pwd` - -#Install miniconda python -if [ $TRAVIS_OS_NAME == "osx" ]; then - wget -q https://repo.continuum.io/miniconda/Miniconda3-latest-MacOSX-x86_64.sh -O Miniconda3.sh -else - wget -q https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O Miniconda3.sh --no-check-certificate # Default OS versions are old and have SSL / CERT issues -fi - -bash Miniconda3.sh -b - -# Create a new conda environment with the target python version -~/miniconda3/bin/conda install conda-build -y -~/miniconda3/bin/conda create -q -y --name testenv python=$CONDA_PY cython numpy pytest psutil pip - -# activate testenv environment -source ~/miniconda3/bin/activate testenv - -conda config --add channels defaults -conda config --add channels bioconda -conda config --add channels conda-forge - -# pin versions, so that tests do not fail when pysam/htslib out of step -# add htslib dependencies -# NB: force conda-forge:blas due to conda/conda#7548 -conda install -y "samtools>=1.11" "bcftools>=1.11" "htslib>=1.11" xz curl bzip2 "conda-forge::blas=*=openblas" - -# As HTSLIB_MODE is (defaulted to) 'shared', ensure we don't pick up -# the external headers from the Conda-installed htslib package. -mv $CONDA_PREFIX/include/htslib $CONDA_PREFIX/include/htslib-disable - -export HTSLIB_CONFIGURE_OPTIONS="--disable-libcurl" - -echo "show samtools, htslib, and bcftools versions" -samtools --version -htsfile --version -bcftools --version - -# Try building conda recipe first -~/miniconda3/bin/conda-build devtools/conda-recipe/ --python=$CONDA_PY - -# install code from the repository via setup.py -echo -echo "============ installing via setup.py from repository ============" -echo -python setup.py install || exit - -# create auxiliary data -echo -echo 'building test data' -echo -make -C tests/pysam_data -make -C tests/cbcf_data -make -C tests/tabix_data - -# echo any limits that are in place -ulimit -a - -# run tests -pytest - -if [ $? != 0 ]; then - exit 1 -fi - -# build source tar-ball. Make sure to run 'build' target so that .pyx -# files are cythonized. -python setup.py build sdist - -if [ $? != 0 ]; then - exit 1 -fi - -# check for presence of config.h files -echo "checking for presence of config.h files in tar-ball" -tar -tvzf dist/pysam-*.tar.gz | grep "config.h$" - -if [ $? != 1 ]; then - echo "ERROR: found config.h in tar-ball" - tar -tvzf dist/pysam-*.tar.gz | grep "config.h%" - exit 1 -fi - -# test pip installation from tar-ball with cython -echo "pip installing with cython" -pip install --verbose --no-deps --no-binary=:all: dist/pysam-*.tar.gz - -if [ $? != 0 ]; then - exit 1 -fi - -# attempt pip installation without cython -echo "pip installing without cython" -~/miniconda3/bin/conda remove -y cython -~/miniconda3/bin/conda list -echo "python is" `which python` -pip install --verbose --no-deps --no-binary=:all: --force-reinstall --upgrade dist/pysam-*.tar.gz - -if [ $? != 0 ]; then - exit 1 -fi - -# attempt pip installation without cython and without -# command line options -echo "pip installing without cython and no configure options" -export HTSLIB_CONFIGURE_OPTIONS="" -pip install --verbose --no-deps --no-binary=:all: --force-reinstall --upgrade dist/pysam-*.tar.gz - -if [ $? != 0 ]; then - exit 1 -fi diff --git a/doc/installation.rst b/doc/installation.rst index 6282fb5..9cce3f0 100644 --- a/doc/installation.rst +++ b/doc/installation.rst @@ -12,8 +12,8 @@ Conda installation To install pysam in your current conda_ environment, type:: - conda config --add channels r conda config --add channels bioconda + conda config --add channels conda-forge conda install pysam This will install pysam from the bioconda_ channel and automatically diff --git a/import/pysam.h b/import/pysam.h index da07281..94cf203 100644 --- a/import/pysam.h +++ b/import/pysam.h @@ -69,6 +69,12 @@ extern int @pysam@_main(int argc, char *argv[]); #define bam_smpl_destroy @pysam@_bam_smpl_destroy #define read_file_list @pysam@_read_file_list +/*! A non-static error() function name is used in bcftools, which collides + with glibc's error() function and leads to the wrong function being called + on some platforms. #define this name with a prefix to avoid this collision. + */ +#define error @pysam@_error + #endif #endif diff --git a/pyproject.toml b/pyproject.toml index 1f89f9b..aa99371 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ license = { text = "MIT License" } authors = [ { name = "Andreas Heger", email = "andreas.heger@gmail.com"} ] -requires-python = ">=3.6" +requires-python = ">=3.8" dynamic = [ "classifiers", @@ -22,7 +22,7 @@ requires = ["setuptools>=59.0", "Cython>=0.29.12,<4"] build-backend = "setuptools.build_meta:__legacy__" [tool.cibuildwheel] -before-all = "{project}/devtools/install-prerequisites.sh" +before-all = "{project}/devtools/install-prerequisites.sh {project}/devtools/emulate-tools.py /usr/local/bin" # Necessary until we build libhts.a out-of-tree from within build_temp before-build = "make -C {project}/htslib distclean" @@ -32,7 +32,7 @@ test-command = "REF_PATH=: pytest {project}/tests" [tool.tox] legacy_tox_ini = """ [tox] - envlist = py36, py311 + envlist = py38, py311 [testenv] deps = pytest diff --git a/pysam/conftest_cstd.c b/pysam/conftest_cstd.c new file mode 100644 index 0000000..5824847 --- /dev/null +++ b/pysam/conftest_cstd.c @@ -0,0 +1,11 @@ +#include +#include + +// C++-style comment; for-decl; optind + +int main(int argc, char **argv) { + int sum = 0; + for (int i = optind; i < argc; i++) + sum += strlen(argv[i]); + return sum; +} diff --git a/pysam/libcalignedsegment.pxd b/pysam/libcalignedsegment.pxd index e14cbb1..1c36319 100644 --- a/pysam/libcalignedsegment.pxd +++ b/pysam/libcalignedsegment.pxd @@ -31,6 +31,18 @@ from pysam.libcalignmentfile cimport AlignmentFile, AlignmentHeader ctypedef AlignmentFile AlignmentFile_t +cdef class _AlignedSegment_Cache: # For internal use only + cdef clear_query_sequences(self) + cdef clear_query_qualities(self) + + cdef object query_sequence + cdef object query_alignment_sequence + cdef object query_qualities + cdef object query_qualities_str + cdef object query_alignment_qualities + cdef object query_alignment_qualities_str + + # Note: need to declare all C fields and methods here cdef class AlignedSegment: @@ -41,10 +53,11 @@ cdef class AlignedSegment: cdef readonly AlignmentHeader header # caching of array properties for quick access - cdef object cache_query_qualities - cdef object cache_query_alignment_qualities - cdef object cache_query_sequence - cdef object cache_query_alignment_sequence + cdef _AlignedSegment_Cache cache + + cdef object unused1 + cdef object unused2 + cdef object unused3 # add an alignment tag with value to the AlignedSegment # an existing tag of the same name will be replaced. diff --git a/pysam/libcalignedsegment.pyi b/pysam/libcalignedsegment.pyi index fdd142a..5665d3e 100644 --- a/pysam/libcalignedsegment.pyi +++ b/pysam/libcalignedsegment.pyi @@ -83,6 +83,7 @@ class AlignedSegment: template_length: int query_sequence: Optional[str] query_qualities: Optional[array] + query_qualities_str: Optional[str] bin: int is_paired: bool is_proper_pair: bool @@ -121,6 +122,8 @@ class AlignedSegment: @property def query_alignment_qualities(self) -> Optional[array]: ... @property + def query_alignment_qualities_str(self) -> Optional[str]: ... + @property def query_alignment_start(self) -> int: ... @property def query_alignment_end(self) -> int: ... @@ -135,9 +138,28 @@ class AlignedSegment: def get_reference_sequence(self) -> str: ... def get_forward_sequence(self) -> Optional[str]: ... def get_forward_qualities(self) -> Optional[array]: ... - def get_aligned_pairs( - self, matches_only: bool = ..., with_seq: bool = ..., with_cigar: bool = ... - ) -> List[Tuple[int, int]]: ... + + @overload + def get_aligned_pairs(self, matches_only: Literal[True], with_seq: Literal[False] = ..., with_cigar: Literal[False] = ...) -> List[Tuple[int, int]]: ... + @overload + def get_aligned_pairs(self, matches_only: Literal[True], with_seq: Literal[False], with_cigar: Literal[True]) -> List[Tuple[int, int, CIGAR_OPS]]: ... + @overload + def get_aligned_pairs(self, matches_only: Literal[True], with_seq: Literal[True], with_cigar: Literal[False] = ...) -> List[Tuple[int, int, str]]: ... + @overload + def get_aligned_pairs(self, matches_only: Literal[True], with_seq: Literal[True], with_cigar: Literal[True]) -> List[Tuple[int, int, str, CIGAR_OPS]]: ... + + @overload + def get_aligned_pairs(self, matches_only: bool = ..., with_seq: Literal[False] = ..., with_cigar: Literal[False] = ...) -> List[Tuple[Optional[int], Optional[int]]]: ... + @overload + def get_aligned_pairs(self, matches_only: bool, with_seq: Literal[False], with_cigar: Literal[True]) -> List[Tuple[Optional[int], Optional[int], CIGAR_OPS]]: ... + @overload + def get_aligned_pairs(self, matches_only: bool, with_seq: Literal[True], with_cigar: Literal[False] = ...) -> List[Tuple[Optional[int], Optional[int], Optional[str]]]: ... + @overload + def get_aligned_pairs(self, matches_only: bool, with_seq: Literal[True], with_cigar: Literal[True]) -> List[Tuple[Optional[int], Optional[int], Optional[str], CIGAR_OPS]]: ... + + @overload + def get_aligned_pairs(self, matches_only: bool = ..., with_seq: bool = ..., with_cigar: bool = ...) -> List[Tuple]: ... + def get_blocks(self) -> List[Tuple[int, int]]: ... def get_overlap(self, start: int, end: int) -> Optional[int]: ... def get_cigar_stats(self) -> Tuple[array, array]: ... diff --git a/pysam/libcalignedsegment.pyx b/pysam/libcalignedsegment.pyx index c9c2773..ceffc3e 100644 --- a/pysam/libcalignedsegment.pyx +++ b/pysam/libcalignedsegment.pyx @@ -65,7 +65,6 @@ cimport cython from cpython cimport array as c_array from cpython cimport PyBytes_FromStringAndSize from libc.string cimport memset, strchr -from cpython cimport array as c_array from libc.stdint cimport INT8_MIN, INT16_MIN, INT32_MIN, \ INT8_MAX, INT16_MAX, INT32_MAX, \ UINT8_MAX, UINT16_MAX, UINT32_MAX @@ -73,8 +72,6 @@ from libc.stdint cimport INT8_MIN, INT16_MIN, INT32_MIN, \ from pysam.libchtslib cimport HTS_IDX_NOCOOR from pysam.libcutils cimport force_bytes, force_str, \ charptr_to_str, charptr_to_bytes -from pysam.libcutils cimport qualities_to_qualitystring, qualitystring_to_array, \ - array_to_qualitystring # Constants for binary tag conversion cdef char * htslib_types = 'cCsSiIf' @@ -556,28 +553,6 @@ cdef inline bytes getSequenceInRange(bam1_t *src, return charptr_to_bytes(seq) -cdef inline object getQualitiesInRange(bam1_t *src, - uint32_t start, - uint32_t end): - """return python array of quality values from a bam1_t object""" - - cdef uint8_t * p - cdef uint32_t k - - p = pysam_bam_get_qual(src) - if p[0] == 0xff: - return None - - # 'B': unsigned char - cdef c_array.array result = array.array('B', [0]) - c_array.resize(result, end - start) - - # copy data - memcpy(result.data.as_voidptr, &p[start], end - start) - - return result - - ##################################################################### ## factory methods for instantiating extension classes cdef class AlignedSegment @@ -588,6 +563,7 @@ cdef AlignedSegment makeAlignedSegment(bam1_t *src, cdef AlignedSegment dest = AlignedSegment.__new__(AlignedSegment) dest._delegate = bam_dup1(src) dest.header = header + dest.cache = _AlignedSegment_Cache() return dest @@ -920,6 +896,22 @@ cdef _alignedpairs_with_seq_cigar(qpos, pos, ref_seq, uint32_t r_idx, int op): return (qpos, pos, ref_base, CIGAR_OPS(op)) +cdef class _AlignedSegment_Cache: + def __cinit__(self): + self.clear_query_sequences() + self.clear_query_qualities() + + cdef clear_query_sequences(self): + self.query_sequence = NotImplemented + self.query_alignment_sequence = NotImplemented + + cdef clear_query_qualities(self): + self.query_qualities = NotImplemented + self.query_qualities_str = NotImplemented + self.query_alignment_qualities = NotImplemented + self.query_alignment_qualities_str = NotImplemented + + cdef class AlignedSegment: '''Class representing an aligned segment. @@ -971,10 +963,7 @@ cdef class AlignedSegment: self._delegate.core.mpos = -1 # caching for selected fields - self.cache_query_qualities = None - self.cache_query_alignment_qualities = None - self.cache_query_sequence = None - self.cache_query_alignment_sequence = None + self.cache = _AlignedSegment_Cache() self.header = header @@ -1122,6 +1111,7 @@ cdef class AlignedSegment: cdef AlignedSegment dest = cls.__new__(cls) dest._delegate = calloc(1, sizeof(bam1_t)) dest.header = header + dest.cache = _AlignedSegment_Cache() cdef kstring_t line line.l = line.m = len(sam) @@ -1313,15 +1303,30 @@ cdef class AlignedSegment: empty string. ''' def __get__(self): - c = self.cigartuples - if c is None: + cdef bam1_t *src = self._delegate + if pysam_get_n_cigar(src) == 0: return None - # reverse order - else: - return "".join([ "%i%c" % (y,CODE2CIGAR[x]) for x,y in c]) + + cdef kstring_t buf + buf.l = buf.m = 0 + buf.s = NULL + + cdef uint32_t *cigar_p = pysam_bam_get_cigar(src) + cdef uint32_t op, l + cdef int k + for k from 0 <= k < pysam_get_n_cigar(src): + op = cigar_p[k] & BAM_CIGAR_MASK + l = cigar_p[k] >> BAM_CIGAR_SHIFT + kputl(l, &buf) + kputc(CODE2CIGAR[op], &buf) + + try: + return buf.s[:buf.l].decode("ascii") + finally: + free(buf.s) def __set__(self, cigar): - if cigar is None or len(cigar) == 0: + if cigar is None or len(cigar) == 0 or cigar == "*": self.cigartuples = [] else: parts = CIGAR_REGEX.findall(cigar) @@ -1420,19 +1425,20 @@ cdef class AlignedSegment: has aligned the read to the reverse strand.) """ def __get__(self): - if self.cache_query_sequence: - return self.cache_query_sequence + if self.cache.query_sequence is not NotImplemented: + return self.cache.query_sequence cdef bam1_t * src cdef char * s src = self._delegate if src.core.l_qseq == 0: + self.cache.query_sequence = None return None - self.cache_query_sequence = force_str(getSequenceInRange( + self.cache.query_sequence = force_str(getSequenceInRange( src, 0, src.core.l_qseq)) - return self.cache_query_sequence + return self.cache.query_sequence def __set__(self, seq): # samtools manages sequence and quality length memory together @@ -1443,7 +1449,7 @@ cdef class AlignedSegment: cdef int l, k cdef Py_ssize_t nbytes_new, nbytes_old - if seq == None: + if seq is None or len(seq) == 0 or seq == "*": l = 0 else: l = len(seq) @@ -1484,11 +1490,8 @@ cdef class AlignedSegment: p = pysam_bam_get_qual(src) memset(p, 0xff, l) - self.cache_query_sequence = force_str(seq) - - # clear cached values for quality values - self.cache_query_qualities = None - self.cache_query_alignment_qualities = None + self.cache.clear_query_sequences() + self.cache.clear_query_qualities() property query_qualities: """read sequence base qualities, including :term:`soft clipped` bases @@ -1501,62 +1504,113 @@ cdef class AlignedSegment: Note that to set quality scores the sequence has to be set beforehand as this will determine the expected length of the - quality score array. - - This method raises a ValueError if the length of the - quality scores and the sequence are not the same. + quality score array. Setting will raise a ValueError if the + length of the new quality scores is not the same as the + length of the existing sequence. + Quality scores to be set may be specified as a Python array + or other iterable of ints, or as a string of ASCII-encooded + FASTQ/SAM-style base quality characters. """ def __get__(self): + if self.cache.query_qualities is not NotImplemented: + return self.cache.query_qualities - if self.cache_query_qualities: - return self.cache_query_qualities + cdef bam1_t *src = self._delegate + cdef int qual_len = src.core.l_qseq + cdef uint8_t *qual = pysam_bam_get_qual(src) - cdef bam1_t * src - cdef char * q + if qual_len == 0 or qual[0] == 0xff: + self.cache.query_qualities = None + return None - src = self._delegate + cdef c_array.array qual_array = array.array('B') + c_array.resize(qual_array, qual_len) + memcpy(qual_array.data.as_uchars, qual, qual_len) + self.cache.query_qualities = qual_array + return qual_array - if src.core.l_qseq == 0: - return None + def __set__(self, new_qual): + if isinstance(new_qual, str): + self.query_qualities_str = new_qual + return - self.cache_query_qualities = getQualitiesInRange(src, 0, src.core.l_qseq) - return self.cache_query_qualities + cdef bam1_t *src = self._delegate + cdef int qual_len = src.core.l_qseq + cdef uint8_t *qual = pysam_bam_get_qual(src) - def __set__(self, qual): + cdef int new_qual_len = len(new_qual) if new_qual is not None else 0 + if new_qual_len == 0: + if qual_len != 0: memset(qual, 0xff, qual_len) + self.cache.clear_query_qualities() + return - # note that memory is already allocated via setting the sequence - # hence length match of sequence and quality needs is checked. - cdef bam1_t * src - cdef uint8_t * p - cdef int l + if new_qual_len != qual_len: + raise ValueError(f"quality ({new_qual_len}) and sequence ({qual_len}) length mismatch") - src = self._delegate - p = pysam_bam_get_qual(src) - if qual is None or len(qual) == 0: - # if absent and there is a sequence: set to 0xff - memset(p, 0xff, src.core.l_qseq) + if isinstance(new_qual, array.array) and new_qual.typecode == 'B': + memcpy(qual, ( new_qual).data.as_uchars, qual_len) + self.cache.clear_query_qualities() return - # check for length match - l = len(qual) - if src.core.l_qseq != l: - raise ValueError( - "quality and sequence mismatch: %i != %i" % - (l, src.core.l_qseq)) + cdef uint8_t *s = qual + cdef uint8_t q + for q in new_qual: + s[0] = q + s += 1 + + self.cache.clear_query_qualities() + + property query_qualities_str: + """read sequence base qualities, including :term:`soft clipped` bases, + returned as an ASCII-encoded string similar to that in FASTQ or SAM files, + or None if base qualities are not present. + + Note that to set quality scores the sequence has to be set beforehand + as this will determine the expected length of the quality score string. + Setting will raise a ValueError if the length of the new quality scores + is not the same as the length of the existing sequence. + """ + def __get__(self): + if self.cache.query_qualities_str is not NotImplemented: + return self.cache.query_qualities_str + + cdef bam1_t *src = self._delegate + cdef int qual_len = src.core.l_qseq + cdef uint8_t *qual = pysam_bam_get_qual(src) - # create a python array object filling it - # with the quality scores + if qual_len == 0 or qual[0] == 0xff: + self.cache.query_qualities_str = None + return None - # NB: should avoid this copying if qual is - # already of the correct type. - cdef c_array.array result = c_array.array('B', qual) + cdef bytes qual_bytes = qual[:qual_len] + cdef char *s = qual_bytes + cdef int i + for i in range(qual_len): s[i] += 33 - # copy data - memcpy(p, result.data.as_voidptr, l) + self.cache.query_qualities_str = qual_bytes.decode('ascii') + return self.cache.query_qualities_str - # save in cache - self.cache_query_qualities = qual + def __set__(self, new_qual): + cdef bam1_t *src = self._delegate + cdef int qual_len = src.core.l_qseq + cdef uint8_t *qual = pysam_bam_get_qual(src) + + cdef int new_qual_len = len(new_qual) if new_qual is not None else 0 + if new_qual_len == 0 or new_qual == "*": + if qual_len != 0: memset(qual, 0xff, qual_len) + self.cache.clear_query_qualities() + return + + if new_qual_len != qual_len: + raise ValueError(f"quality ({new_qual_len}) and sequence ({qual_len}) length mismatch") + + cdef bytes new_qual_bytes = new_qual.encode('ascii') + cdef const char *s = new_qual_bytes + cdef int i + for i in range(qual_len): qual[i] = s[i] - 33 + + self.cache.clear_query_qualities() property bin: """properties bin""" @@ -1735,8 +1789,8 @@ cdef class AlignedSegment: """ def __get__(self): - if self.cache_query_alignment_sequence: - return self.cache_query_alignment_sequence + if self.cache.query_alignment_sequence is not NotImplemented: + return self.cache.query_alignment_sequence cdef bam1_t * src cdef uint32_t start, end @@ -1744,14 +1798,14 @@ cdef class AlignedSegment: src = self._delegate if src.core.l_qseq == 0: + self.cache.query_alignment_sequence = None return None start = getQueryStart(src) end = getQueryEnd(src) - self.cache_query_alignment_sequence = force_str( - getSequenceInRange(src, start, end)) - return self.cache_query_alignment_sequence + self.cache.query_alignment_sequence = force_str(getSequenceInRange(src, start, end)) + return self.cache.query_alignment_sequence property query_alignment_qualities: """aligned query sequence quality values (None if not present). These @@ -1766,26 +1820,44 @@ cdef class AlignedSegment: needs to be subtracted. This property is read-only. - """ def __get__(self): + if self.cache.query_alignment_qualities is not NotImplemented: + return self.cache.query_alignment_qualities - if self.cache_query_alignment_qualities: - return self.cache_query_alignment_qualities + cdef object full_qual = self.query_qualities + if full_qual is None: + self.cache.query_alignment_qualities = None + return None - cdef bam1_t * src - cdef uint32_t start, end + cdef bam1_t *src = self._delegate + cdef uint32_t start = getQueryStart(src) + cdef uint32_t end = getQueryEnd(src) + self.cache.query_alignment_qualities = full_qual[start:end] + return self.cache.query_alignment_qualities - src = self._delegate + property query_alignment_qualities_str: + """aligned query sequence quality values, returned as an ASCII-encoded string + similar to that in FASTQ or SAM files, or None if base qualities are not present. + These are the quality values that correspond to :attr:`query_alignment_sequence`, + i.e., excluding qualities corresponding to soft-clipped bases. - if src.core.l_qseq == 0: + This property is read-only. + """ + def __get__(self): + if self.cache.query_alignment_qualities_str is not NotImplemented: + return self.cache.query_alignment_qualities_str + + cdef object full_qual = self.query_qualities_str + if full_qual is None: + self.cache.query_alignment_qualities_str = None return None - start = getQueryStart(src) - end = getQueryEnd(src) - self.cache_query_alignment_qualities = \ - getQualitiesInRange(src, start, end) - return self.cache_query_alignment_qualities + cdef bam1_t *src = self._delegate + cdef uint32_t start = getQueryStart(src) + cdef uint32_t end = getQueryEnd(src) + self.cache.query_alignment_qualities_str = full_qual[start:end] + return self.cache.query_alignment_qualities_str property query_alignment_start: """start index of the aligned query portion of the sequence (0-based, @@ -2190,29 +2262,29 @@ cdef class AlignedSegment: field will always be 0. (Accessing this field via index -1 avoids changes if more CIGAR operators are added in future.) - +-----+----------------+--------+ - |M |pysam.CMATCH |0 | - +-----+----------------+--------+ - |I |pysam.CINS |1 | - +-----+----------------+--------+ - |D |pysam.CDEL |2 | - +-----+----------------+--------+ - |N |pysam.CREF_SKIP |3 | - +-----+----------------+--------+ - |S |pysam.CSOFT_CLIP|4 | - +-----+----------------+--------+ - |H |pysam.CHARD_CLIP|5 | - +-----+----------------+--------+ - |P |pysam.CPAD |6 | - +-----+----------------+--------+ - |= |pysam.CEQUAL |7 | - +-----+----------------+--------+ - |X |pysam.CDIFF |8 | - +-----+----------------+--------+ - |B |pysam.CBACK |9 | - +-----+----------------+--------+ - |NM |NM tag |10 or -1| - +-----+----------------+--------+ + +-----+--------------------------+--------+ + |M |pysam.CIGAR_OPS.CMATCH |0 | + +-----+--------------------------+--------+ + |I |pysam.CIGAR_OPS.CINS |1 | + +-----+--------------------------+--------+ + |D |pysam.CIGAR_OPS.CDEL |2 | + +-----+--------------------------+--------+ + |N |pysam.CIGAR_OPS.CREF_SKIP |3 | + +-----+--------------------------+--------+ + |S |pysam.CIGAR_OPS.CSOFT_CLIP|4 | + +-----+--------------------------+--------+ + |H |pysam.CIGAR_OPS.CHARD_CLIP|5 | + +-----+--------------------------+--------+ + |P |pysam.CIGAR_OPS.CPAD |6 | + +-----+--------------------------+--------+ + |= |pysam.CIGAR_OPS.CEQUAL |7 | + +-----+--------------------------+--------+ + |X |pysam.CIGAR_OPS.CDIFF |8 | + +-----+--------------------------+--------+ + |B |pysam.CIGAR_OPS.CBACK |9 | + +-----+--------------------------+--------+ + |NM |NM tag |10 or -1| + +-----+--------------------------+--------+ If no cigar string is present, empty arrays will be returned. @@ -2267,27 +2339,27 @@ cdef class AlignedSegment: The operations are: - +-----+----------------+-----+ - |M |pysam.CMATCH |0 | - +-----+----------------+-----+ - |I |pysam.CINS |1 | - +-----+----------------+-----+ - |D |pysam.CDEL |2 | - +-----+----------------+-----+ - |N |pysam.CREF_SKIP |3 | - +-----+----------------+-----+ - |S |pysam.CSOFT_CLIP|4 | - +-----+----------------+-----+ - |H |pysam.CHARD_CLIP|5 | - +-----+----------------+-----+ - |P |pysam.CPAD |6 | - +-----+----------------+-----+ - |= |pysam.CEQUAL |7 | - +-----+----------------+-----+ - |X |pysam.CDIFF |8 | - +-----+----------------+-----+ - |B |pysam.CBACK |9 | - +-----+----------------+-----+ + +-----+--------------------------+-----+ + |M |pysam.CIGAR_OPS.CMATCH |0 | + +-----+--------------------------+-----+ + |I |pysam.CIGAR_OPS.CINS |1 | + +-----+--------------------------+-----+ + |D |pysam.CIGAR_OPS.CDEL |2 | + +-----+--------------------------+-----+ + |N |pysam.CIGAR_OPS.CREF_SKIP |3 | + +-----+--------------------------+-----+ + |S |pysam.CIGAR_OPS.CSOFT_CLIP|4 | + +-----+--------------------------+-----+ + |H |pysam.CIGAR_OPS.CHARD_CLIP|5 | + +-----+--------------------------+-----+ + |P |pysam.CIGAR_OPS.CPAD |6 | + +-----+--------------------------+-----+ + |= |pysam.CIGAR_OPS.CEQUAL |7 | + +-----+--------------------------+-----+ + |X |pysam.CIGAR_OPS.CDIFF |8 | + +-----+--------------------------+-----+ + |B |pysam.CIGAR_OPS.CBACK |9 | + +-----+--------------------------+-----+ .. note:: The output is a list of (operation, length) tuples, such as @@ -2794,11 +2866,11 @@ cdef class AlignedSegment: def __set__(self, v): self.query_sequence = v property qual: - """deprecated, use :attr:`query_qualities` instead.""" + """deprecated, use :attr:`query_qualities` or :attr:`query_qualities_str` instead.""" def __get__(self): - return array_to_qualitystring(self.query_qualities) + return self.query_qualities_str def __set__(self, v): - self.query_qualities = qualitystring_to_array(v) + self.query_qualities_str = v property alen: """deprecated, use :attr:`reference_length` instead.""" def __get__(self): @@ -2822,34 +2894,24 @@ cdef class AlignedSegment: instead.""" def __get__(self): return self.query_alignment_sequence - def __set__(self, v): - self.query_alignment_sequence = v property qqual: - """deprecated, use :attr:`query_alignment_qualities` + """deprecated, use :attr:`query_alignment_qualities` or :attr:`query_alignment_qualities_str` instead.""" def __get__(self): - return array_to_qualitystring(self.query_alignment_qualities) - def __set__(self, v): - self.query_alignment_qualities = qualitystring_to_array(v) + return self.query_alignment_qualities_str property qstart: """deprecated, use :attr:`query_alignment_start` instead.""" def __get__(self): return self.query_alignment_start - def __set__(self, v): - self.query_alignment_start = v property qend: """deprecated, use :attr:`query_alignment_end` instead.""" def __get__(self): return self.query_alignment_end - def __set__(self, v): - self.query_alignment_end = v property qlen: """deprecated, use :attr:`query_alignment_length` instead.""" def __get__(self): return self.query_alignment_length - def __set__(self, v): - self.query_alignment_length = v property mrnm: """deprecated, use :attr:`next_reference_id` instead.""" def __get__(self): @@ -3409,10 +3471,16 @@ cpdef enum SAM_FLAGS: FSUPPLEMENTARY = 2048 +# TODO Remove these and remove the enumerators from __all__ +globals().update(getattr(CIGAR_OPS, "__members__")) +globals().update(getattr(SAM_FLAGS, "__members__")) + + __all__ = [ "AlignedSegment", "PileupColumn", "PileupRead", + "CIGAR_OPS", "CMATCH", "CINS", "CDEL", @@ -3423,6 +3491,7 @@ __all__ = [ "CEQUAL", "CDIFF", "CBACK", + "SAM_FLAGS", "FPAIRED", "FPROPER_PAIR", "FUNMAP", diff --git a/pysam/libcalignmentfile.pyx b/pysam/libcalignmentfile.pyx index 9a92d07..4373a75 100644 --- a/pysam/libcalignmentfile.pyx +++ b/pysam/libcalignmentfile.pyx @@ -65,13 +65,13 @@ import re import warnings import array from libc.errno cimport errno, EPIPE -from libc.string cimport strcmp, strpbrk, strerror +from libc.string cimport strcmp, strpbrk from libc.stdint cimport INT32_MAX from cpython cimport array as c_array from pysam.libcutils cimport force_bytes, force_str, charptr_to_str -from pysam.libcutils cimport encode_filename, from_string_and_size +from pysam.libcutils cimport OSError_from_errno, encode_filename, from_string_and_size from pysam.libcalignedsegment cimport makeAlignedSegment, makePileupColumn from pysam.libchtslib cimport HTSFile, hisremote, sam_index_load2, sam_index_load3, \ HTS_IDX_SAVE_REMOTE, HTS_IDX_SILENT_FAIL @@ -302,7 +302,7 @@ cdef class AlignmentHeader(object): data = header_dict[record] if not isinstance(data, VALID_HEADER_TYPES[record]): raise ValueError( - "invalid type for record %s: %s, expected %s".format( + "invalid type for record {}: {}, expected {}".format( record, type(data), VALID_HEADER_TYPES[record])) if isinstance(data, Mapping): lines.append(build_header_line(data, record)) @@ -818,6 +818,7 @@ cdef class AlignmentFile(HTSFile): cdef char *cindexname = NULL cdef char *cmode = NULL cdef bam_hdr_t * hdr = NULL + cdef int ret if threads > 1 and ignore_truncation: # This won't raise errors if reaching a truncated alignment, @@ -921,9 +922,7 @@ cdef class AlignmentFile(HTSFile): if self.htsfile == NULL: if errno: - raise IOError(errno, "could not open alignment file `{}`: {}".format( - force_str(filename), - force_str(strerror(errno)))) + raise OSError_from_errno("Could not open alignment file", filename) else: raise ValueError("could not open alignment file `{}`".format(force_str(filename))) if format_options and len(format_options): @@ -939,7 +938,9 @@ cdef class AlignmentFile(HTSFile): if "b" in mode or "c" in mode or "h" in mode: hdr = self.header.ptr with nogil: - sam_hdr_write(self.htsfile, hdr) + ret = sam_hdr_write(self.htsfile, hdr) + if ret < 0: + raise OSError_from_errno("Could not write headers", filename) elif mode[0] == "r": # open file for reading @@ -947,8 +948,7 @@ cdef class AlignmentFile(HTSFile): if self.htsfile == NULL: if errno: - raise IOError(errno, "could not open alignment file `{}`: {}".format(force_str(filename), - force_str(strerror(errno)))) + raise OSError_from_errno("Could not open alignment file", filename) else: raise ValueError("could not open alignment file `{}`".format(force_str(filename))) @@ -1014,7 +1014,7 @@ cdef class AlignmentFile(HTSFile): if not self.index and (cindexname or require_index): if errno: - raise IOError(errno, force_str(strerror(errno))) + raise OSError_from_errno("Could not open index file", self.index_filename) else: raise IOError('unable to open index file `%s`' % self.index_filename) @@ -1679,12 +1679,11 @@ cdef class AlignmentFile(HTSFile): self.header = None - if ret < 0: - global errno - if errno == EPIPE: - errno = 0 + if ret < 0 and errno != EPIPE: + if isinstance(self.filename, (str, bytes)): + raise OSError_from_errno("Closing failed", self.filename) else: - raise IOError(errno, force_str(strerror(errno))) + raise OSError_from_errno("Closing failed") def __dealloc__(self): cdef int ret = 0 @@ -1703,12 +1702,11 @@ cdef class AlignmentFile(HTSFile): bam_destroy1(self.b) self.b = NULL - if ret < 0: - global errno - if errno == EPIPE: - errno = 0 + if ret < 0 and errno != EPIPE: + if isinstance(self.filename, (str, bytes)): + raise OSError_from_errno("Closing failed", self.filename) else: - raise IOError(errno, force_str(strerror(errno))) + raise OSError_from_errno("Closing failed") cpdef int write(self, AlignedSegment read) except -1: ''' @@ -1742,8 +1740,7 @@ cdef class AlignmentFile(HTSFile): # when ret == -1 we get a "SystemError: error return without # exception set". if ret < 0: - raise IOError( - "sam_write1 failed with error code {}".format(ret)) + raise IOError("sam_write1 failed with error code {}".format(ret)) return ret @@ -2298,14 +2295,14 @@ cdef class IteratorRowSelection(IteratorRow): # end iteration if out of positions if self.current_pos >= len(self.positions): return -1 + cdef int ret cdef uint64_t pos = self.positions[self.current_pos] with nogil: - bgzf_seek(hts_get_bgzfp(self.htsfile), - pos, - 0) + ret = bgzf_seek(hts_get_bgzfp(self.htsfile), pos, 0) + if ret < 0: + raise OSError_from_errno("Can't seek", self.samfile.filename) self.current_pos += 1 - cdef int ret cdef bam_hdr_t * hdr = self.header.ptr with nogil: ret = sam_read1(self.htsfile, @@ -2911,7 +2908,7 @@ cdef class IndexedReads: with nogil: self.htsfile = hts_open(cfilename, 'r') if self.htsfile == NULL: - raise OSError("unable to reopen htsfile") + raise OSError_from_errno("Unable to reopen file", cfilename) # need to advance in newly opened file to position after header # better: use seek/tell? diff --git a/pysam/libcbcf.pyx b/pysam/libcbcf.pyx index c939887..f1eac7c 100644 --- a/pysam/libcbcf.pyx +++ b/pysam/libcbcf.pyx @@ -143,7 +143,7 @@ cdef tuple METADATA_LENGTHS = ('FIXED', 'VARIABLE', 'A', 'G', 'R') ######################################################################## from pysam.libcutils cimport force_bytes, force_str, charptr_to_str, charptr_to_str_w_len -from pysam.libcutils cimport encode_filename, from_string_and_size, decode_bytes +from pysam.libcutils cimport OSError_from_errno, encode_filename, from_string_and_size, decode_bytes ######################################################################## @@ -2158,10 +2158,12 @@ cdef class VariantHeader(object): quoted = not isinstance(value, unquoted_str) and key not in ("ID", "Number", "Type") key = force_bytes(key) - bcf_hrec_add_key(hrec, key, len(key)) + if bcf_hrec_add_key(hrec, key, len(key)) < 0: + raise MemoryError("Could not allocate VCF header record") value = force_bytes(str(value)) - bcf_hrec_set_val(hrec, hrec.nkeys-1, value, len(value), quoted) + if bcf_hrec_set_val(hrec, hrec.nkeys-1, value, len(value), quoted) < 0: + raise MemoryError("Could not allocate VCF header record") except: bcf_hrec_destroy(hrec) raise @@ -4116,42 +4118,42 @@ cdef class VariantFile(HTSFile): if not self.htsfile or not self.header: return + cdef int ret + # Write header if no records were written if self.htsfile.is_write and not self.header_written: with nogil: - bcf_hdr_write(self.htsfile, self.header.ptr) + ret = bcf_hdr_write(self.htsfile, self.header.ptr) + if ret < 0 and errno != EPIPE: + raise OSError_from_errno("Can't write headers", self.filename) - cdef int ret = hts_close(self.htsfile) + ret = hts_close(self.htsfile) self.htsfile = NULL self.header = self.index = None - if ret < 0: - global errno - if errno == EPIPE: - errno = 0 - else: - raise IOError(errno, force_str(strerror(errno))) + if ret < 0 and errno != EPIPE: + raise OSError_from_errno("Closing failed", self.filename) def close(self): """closes the :class:`pysam.VariantFile`.""" if not self.htsfile: return + cdef int ret + # Write header if no records were written if self.htsfile.is_write and not self.header_written: with nogil: - bcf_hdr_write(self.htsfile, self.header.ptr) + ret = bcf_hdr_write(self.htsfile, self.header.ptr) + if ret < 0 and errno != EPIPE: + raise OSError_from_errno("Can't write headers", self.filename) - cdef int ret = hts_close(self.htsfile) + ret = hts_close(self.htsfile) self.htsfile = NULL self.header = self.index = None - if ret < 0: - global errno - if errno == EPIPE: - errno = 0 - else: - raise IOError(errno, force_str(strerror(errno))) + if ret < 0 and errno != EPIPE: + raise OSError_from_errno("Closing failed", self.filename) def __iter__(self): if not self.is_open: @@ -4188,7 +4190,7 @@ cdef class VariantFile(HTSFile): elif ret == -2: raise IOError('truncated file') elif errno: - raise IOError(errno, strerror(errno)) + raise OSError_from_errno("Unable to fetch next record", self.filename) else: raise IOError('unable to fetch next record') @@ -4335,7 +4337,7 @@ cdef class VariantFile(HTSFile): if not self.htsfile: if errno: - raise IOError(errno, 'could not open variant file `{}`: {}'.format(filename, force_str(strerror(errno)))) + raise OSError_from_errno("Could not open variant file", filename) else: raise ValueError('could not open variant file `{}`'.format(filename)) @@ -4485,15 +4487,19 @@ cdef class VariantFile(HTSFile): raise ValueError('record must not be None') if not self.is_open: - return ValueError('I/O operation on closed file') + raise ValueError('I/O operation on closed file') if not self.htsfile.is_write: raise ValueError('cannot write to a Variantfile opened for reading') + cdef int ret + if not self.header_written: self.header_written = True with nogil: - bcf_hdr_write(self.htsfile, self.header.ptr) + ret = bcf_hdr_write(self.htsfile, self.header.ptr) + if ret < 0: + raise OSError_from_errno("Can't write headers", self.filename) #if record.header is not self.header: # record.translate(self.header) @@ -4506,13 +4512,11 @@ cdef class VariantFile(HTSFile): # Sync END annotation before writing bcf_sync_end(record) - cdef int ret - with nogil: ret = bcf_write1(self.htsfile, self.header.ptr, record.ptr) if ret < 0: - raise IOError(errno, strerror(errno)) + raise OSError_from_errno("Can't write record", self.filename) return ret diff --git a/pysam/libchtslib.pxd b/pysam/libchtslib.pxd index 00075cc..99ea157 100644 --- a/pysam/libchtslib.pxd +++ b/pysam/libchtslib.pxd @@ -1871,10 +1871,10 @@ cdef extern from "htslib/vcf.h" nogil: # bcf_hrec_t *bcf_hdr_get_hrec(const bcf_hdr_t *hdr, int type, const char *key, const char *value, const char *str_class) bcf_hrec_t *bcf_hrec_dup(bcf_hrec_t *hrec) - void bcf_hrec_add_key(bcf_hrec_t *hrec, const char *str, int len) - void bcf_hrec_set_val(bcf_hrec_t *hrec, int i, const char *str, int len, int is_quoted) + int bcf_hrec_add_key(bcf_hrec_t *hrec, const char *str, size_t len) + int bcf_hrec_set_val(bcf_hrec_t *hrec, int i, const char *str, size_t len, int is_quoted) int bcf_hrec_find_key(bcf_hrec_t *hrec, const char *key) - void hrec_add_idx(bcf_hrec_t *hrec, int idx) + int hrec_add_idx(bcf_hrec_t *hrec, int idx) void bcf_hrec_destroy(bcf_hrec_t *hrec) #************************************************************************ diff --git a/pysam/libctabix.pyx b/pysam/libctabix.pyx index b38cf55..deae908 100644 --- a/pysam/libctabix.pyx +++ b/pysam/libctabix.pyx @@ -658,6 +658,12 @@ cdef class TabixIterator: return retval + def _itr_error(self, err: int): + if err == -5: + return IOError("iteration on closed file") + else: + return ValueError(f"iteration failed (error code {err})") + def __next__(self): """python version of next(). @@ -665,10 +671,8 @@ cdef class TabixIterator: """ cdef int retval = self.__cnext__() - if retval == -5: - raise IOError("iteration on closed file") - elif retval < 0: - raise StopIteration + if retval < 0: + raise StopIteration if retval == -1 else self._itr_error(retval) return charptr_to_str(self.buffer.s, self.encoding) @@ -708,10 +712,8 @@ cdef class TabixIteratorParsed(TabixIterator): """ cdef int retval = self.__cnext__() - if retval == -5: - raise IOError("iteration on closed file") - elif retval < 0: - raise StopIteration + if retval < 0: + raise StopIteration if retval == -1 else self._itr_error(retval) return self.parser.parse(self.buffer.s, self.buffer.l) diff --git a/pysam/libcutils.pxd b/pysam/libcutils.pxd index 1bce057..7e2e098 100644 --- a/pysam/libcutils.pxd +++ b/pysam/libcutils.pxd @@ -9,6 +9,8 @@ cpdef parse_region(contig=*, start=*, stop=*, region=*, reference=*, end=*) cdef int libc_whence_from_io(int whence) +cdef OSError_from_errno(message, filename=*) + ######################################################################### # Utility functions for quality string conversions diff --git a/pysam/libcutils.pyx b/pysam/libcutils.pyx index 72eaee5..7a03f3c 100644 --- a/pysam/libcutils.pyx +++ b/pysam/libcutils.pyx @@ -14,12 +14,12 @@ from cpython cimport PyBytes_Check, PyUnicode_Check from cpython cimport array as c_array from libc.errno cimport errno from libc.stdlib cimport calloc, free -from libc.string cimport strncpy +from libc.string cimport strerror, strncpy from libc.stdint cimport INT32_MAX, int32_t from libc.stdio cimport fprintf, stderr, fflush from libc.stdio cimport stdout as c_stdout from posix.fcntl cimport open as c_open, O_WRONLY, O_CREAT, O_TRUNC -from posix.unistd cimport SEEK_SET, SEEK_CUR, SEEK_END +from posix.unistd cimport dup as c_dup, SEEK_SET, SEEK_CUR, SEEK_END, STDOUT_FILENO from pysam.libcsamtools cimport samtools_dispatch, samtools_set_stdout, samtools_set_stderr, \ samtools_close_stdout, samtools_close_stderr, samtools_set_stdout_fn @@ -175,6 +175,12 @@ cdef decode_bytes(bytes s, encoding=None, errors=None): return s.decode(encoding or TEXT_ENCODING, errors or ERROR_HANDLER) +cdef OSError_from_errno(message, filename=None): + cdef int err = errno + if filename is not None: filename = os.fsdecode(filename) + return OSError(err, f"{message}: {strerror(err).decode()}", filename) + + cpdef parse_region(contig=None, start=None, stop=None, @@ -324,7 +330,7 @@ def _pysam_dispatch(collection, stdout_h = c_open(force_bytes(stdout_f), O_WRONLY|O_CREAT|O_TRUNC, 0666) if stdout_h == -1: - raise OSError(errno, "error while opening file for writing", stdout_f) + raise OSError_from_errno("Could not redirect standard output", stdout_f) samtools_set_stdout_fn(force_bytes(stdout_f)) bcftools_set_stdout_fn(force_bytes(stdout_f)) @@ -361,7 +367,8 @@ def _pysam_dispatch(collection, else: samtools_set_stdout_fn("-") bcftools_set_stdout_fn("-") - stdout_h = c_open(b"/dev/null", O_WRONLY) + if catch_stdout is None: stdout_h = c_dup(STDOUT_FILENO) + else: stdout_h = c_open(b"/dev/null", O_WRONLY) # setup the function call to samtools/bcftools main cdef char ** cargs diff --git a/pysam/version.py b/pysam/version.py index bdea9ce..0e26758 100644 --- a/pysam/version.py +++ b/pysam/version.py @@ -1,5 +1,5 @@ # pysam versioning information -__version__ = "0.23.0" +__version__ = "0.23.3" __samtools_version__ = "1.21" __bcftools_version__ = "1.21" diff --git a/samtools/samtools.pysam.h b/samtools/samtools.pysam.h index cb63b60..5fedbea 100644 --- a/samtools/samtools.pysam.h +++ b/samtools/samtools.pysam.h @@ -69,6 +69,12 @@ extern int samtools_main(int argc, char *argv[]); #define bam_smpl_destroy samtools_bam_smpl_destroy #define read_file_list samtools_read_file_list +/*! A non-static error() function name is used in bcftools, which collides + with glibc's error() function and leads to the wrong function being called + on some platforms. #define this name with a prefix to avoid this collision. + */ +#define error samtools_error + #endif #endif diff --git a/setup.py b/setup.py index 64093c2..455b42a 100644 --- a/setup.py +++ b/setup.py @@ -28,9 +28,9 @@ from setuptools.command.sdist import sdist from setuptools.extension import Extension try: - from setuptools.errors import LinkError + from setuptools.errors import CompileError, LinkError except ImportError: - from distutils.errors import LinkError + from distutils.errors import CompileError, LinkError try: from Cython.Distutils import build_ext @@ -99,11 +99,14 @@ def run_nm_defined_symbols(objfile): if symtype not in "UFNWw": if IS_DARWIN: # On macOS, all symbols have a leading underscore - symbols.add(sym.lstrip('_')) + symbols.add(sym[1:] if sym.startswith("_") else sym) else: # Ignore symbols such as _edata (present in all shared objects) if sym[0] not in "_$.@": symbols.add(sym) + # Work around Cython 3.1.2 bug whereby this function is not static + symbols.discard("__pyx_CommonTypesMetaclass_get_module") + return symbols @@ -160,7 +163,7 @@ def write_configvars_header(filename, ext, prefix): log.info("creating %s for '%s' extension", filename, ext.name) with open(filename, "w") as outf: for var, value in config.items(): - outf.write('#define {}_{} "{}"\n'.format(prefix, var, value)) + outf.write(f'#define {prefix}_{var} "{value}"\n') @contextmanager @@ -170,12 +173,12 @@ def set_compiler_envvars(): if var in os.environ: if var == 'CFLAGS' and 'CCSHARED' in sysconfig.get_config_vars(): os.environ[var] += ' ' + sysconfig.get_config_var('CCSHARED') - print("# pysam: (env) {}={}".format(var, os.environ[var])) + print(f"# pysam: (env) {var}={os.environ[var]}") elif var in sysconfig.get_config_vars(): value = sysconfig.get_config_var(var) if var == 'CFLAGS' and 'CCSHARED' in sysconfig.get_config_vars(): value += ' ' + sysconfig.get_config_var('CCSHARED') - print("# pysam: (sysconfig) {}={}".format(var, value)) + print(f"# pysam: (sysconfig) {var}={value}") os.environ[var] = value tmp_vars += [var] @@ -200,8 +203,7 @@ def configure_library(library_dir, env_options=None, options=[]): env_options = "--disable-bz2" if not os.path.exists(configure_script): - raise ValueError( - "configure script {} does not exist".format(configure_script)) + raise ValueError(f"configure script {configure_script!r} does not exist") with changedir(library_dir), set_compiler_envvars(): if env_options is not None: @@ -273,6 +275,22 @@ class cy_build_ext(build_ext): if errors > 0: raise LinkError("symbols defined in multiple extensions") + def c99_compile_args(self): + """Determines whether any compiler flags are needed to ensure C99 compilation.""" + compiler = getattr(self.compiler, "compiler", "C compiler") + if isinstance(compiler, list): compiler = compiler[0] + log.info("checking for %s option to enable C99 features...", compiler) + for flags in [None, ["-std=c99"], ["-std=gnu99"]]: + try: + self.compiler.compile(["pysam/conftest_cstd.c"], output_dir=self.build_temp, extra_preargs=flags) + log.info("%s option to enable C99 features: %s", compiler, " ".join(flags) if flags else "none needed") + return flags + except CompileError: + log.info("(ignoring errors from test probes)") + + log.error("%s cannot compile C99 source code", compiler) + return None + def run(self): if sys.platform == 'darwin': ldshared = os.environ.get('LDSHARED', sysconfig.get_config_var('LDSHARED')) @@ -287,6 +305,19 @@ class cy_build_ext(build_ext): except subprocess.CalledProcessError: log.warning("skipping symbol collision check (invoking nm failed)") + def build_extensions(self): + c99_flags = self.c99_compile_args() + if c99_flags: + executables = {} + for executable in ["compiler", "compiler_so"]: + command = getattr(self.compiler, executable, None) + if command: + if isinstance(command, list): executables[executable] = command + c99_flags + elif isinstance(command, str): executables[executable] = f"{command} {' '.join(c99_flags)}" + self.compiler.set_executables(**executables) + + super().build_extensions() + def build_extension(self, ext): if isinstance(ext, CyExtension) and ext._init_func: @@ -311,7 +342,7 @@ class cy_build_ext(build_ext): ext.extra_link_args += ['-dynamiclib', '-rpath', '@loader_path', '-Wl,-headerpad_max_install_names', - '-Wl,-install_name,%s' % library_path, + f'-Wl,-install_name,{library_path}', '-Wl,-x'] else: if not ext.extra_link_args: @@ -393,24 +424,21 @@ config_headers = ["samtools/config.h", # the .pyx files. If no cython is available, the C-files included in the # distribution will be used. if HAVE_CYTHON: - print("# pysam: cython is available - using cythonize if necessary") + print(f"# pysam: Cython {cython.__version__} is available - using cythonize if necessary") source_pattern = "pysam/libc%s.pyx" else: - print("# pysam: no cython available - using pre-compiled C") + print("# pysam: no Cython available - using pre-compiled C") source_pattern = "pysam/libc%s.c" # Exit if there are no pre-compiled files and no cython available fn = source_pattern % "htslib" if not os.path.exists(fn): raise ValueError( - "no cython installed, but can not find {}." - "Make sure that cython is installed when building " - "from the repository" - .format(fn)) - -print("# pysam: htslib mode is {}".format(HTSLIB_MODE)) -print("# pysam: HTSLIB_CONFIGURE_OPTIONS={}".format( - HTSLIB_CONFIGURE_OPTIONS)) + f"no Cython installed, but cannot find {fn}. " + "Make sure that Cython is installed when building from the repository") + +print(f"# pysam: htslib mode is {HTSLIB_MODE}") +print(f"# pysam: HTSLIB_CONFIGURE_OPTIONS={HTSLIB_CONFIGURE_OPTIONS}") htslib_configure_options = None if HTSLIB_MODE in ['shared', 'separate']: @@ -425,8 +453,7 @@ if HTSLIB_MODE in ['shared', 'separate']: "--disable-libcurl"]) HTSLIB_SOURCE = "builtin" - print("# pysam: htslib configure options: {}".format( - str(htslib_configure_options))) + print(f"# pysam: htslib configure options: {htslib_configure_options}") config_headers += ["htslib/config.h"] if htslib_configure_options is None: @@ -441,7 +468,7 @@ if HTSLIB_MODE in ['shared', 'separate']: htslib_make_options = run_make_print_config() for key, value in htslib_make_options.items(): - print("# pysam: htslib_config {}={}".format(key, value)) + print(f"# pysam: htslib_config {key}={value}") external_htslib_libraries = ['z'] if "LIBS" in htslib_make_options: @@ -478,11 +505,11 @@ elif HTSLIB_MODE == 'shared': htslib_library_dirs = ["."] # when using setup.py develop? htslib_include_dirs = ['htslib'] else: - raise ValueError("unknown HTSLIB value '%s'" % HTSLIB_MODE) + raise ValueError(f"unknown HTSLIB value {HTSLIB_MODE!r}") # build config.py with open(os.path.join("pysam", "config.py"), "w") as outf: - outf.write('HTSLIB = "{}"\n'.format(HTSLIB_SOURCE)) + outf.write(f'HTSLIB = "{HTSLIB_SOURCE}"\n') config_values = collections.defaultdict(int) if HTSLIB_SOURCE == "builtin": @@ -502,8 +529,8 @@ with open(os.path.join("pysam", "config.py"), "w") as outf: "HAVE_LIBDEFLATE", "HAVE_LIBLZMA", "HAVE_MMAP"]: - outf.write("{} = {}\n".format(key, config_values[key])) - print("# pysam: config_option: {}={}".format(key, config_values[key])) + outf.write(f"{key} = {config_values[key]}\n") + print(f"# pysam: config_option: {key}={config_values[key]}") # create empty config.h files if they have not been created automatically # or created by the user: @@ -542,13 +569,15 @@ if os.environ.get("CIBUILDWHEEL", "0") == "1": suffix = sysconfig.get_config_var('EXT_SUFFIX') internal_htslib_libraries = [ - os.path.splitext("chtslib{}".format(suffix))[0]] + os.path.splitext(f"chtslib{suffix}")[0], + ] internal_samtools_libraries = [ - os.path.splitext("csamtools{}".format(suffix))[0], - os.path.splitext("cbcftools{}".format(suffix))[0], + os.path.splitext(f"csamtools{suffix}")[0], + os.path.splitext(f"cbcftools{suffix}")[0], ] internal_pysamutil_libraries = [ - os.path.splitext("cutils{}".format(suffix))[0]] + os.path.splitext(f"cutils{suffix}")[0], + ] libraries_for_pysam_module = external_htslib_libraries + internal_htslib_libraries + internal_pysamutil_libraries diff --git a/tests/AlignedSegment_test.py b/tests/AlignedSegment_test.py index af0e270..0e85060 100644 --- a/tests/AlignedSegment_test.py +++ b/tests/AlignedSegment_test.py @@ -236,6 +236,42 @@ class TestAlignedSegment(ReadTest): self.assertEqual(pysam.qualities_to_qualitystring(a.query_qualities), s[5:10]) + def testClearSequence(self): + a = pysam.AlignedSegment() + a.query_sequence = "ATGC" + self.assertEqual(a.query_sequence, "ATGC") + a.query_sequence = None + self.assertEqual(a.query_length, 0) + + a.query_sequence = "ATGC" + self.assertEqual(a.query_sequence, "ATGC") + a.query_sequence = "" + self.assertEqual(a.query_length, 0) + + a.query_sequence = "ATGC" + self.assertEqual(a.query_sequence, "ATGC") + a.query_sequence = "*" + self.assertEqual(a.query_length, 0) + + def testUpdateSequenceEffects1(self): + a = self.build_read() + a.query_sequence = "ATGCATGC" + a.cigarstring = "1S5M2S" + self.assertEqual(a.query_alignment_sequence, "TGCAT") + + a.query_sequence = "AATTGGCC" + self.assertEqual(a.query_alignment_sequence, "ATTGG") + + def testUpdateSequenceEffects2(self): + a = self.build_read() + a.query_sequence = "ATGCATGC" + a.cigarstring = "1S5M2S" + self.assertEqual(a.query_alignment_sequence, "TGCAT") + + a.query_sequence = "*" + self.assertIsNone(a.query_sequence) + self.assertIsNone(a.query_alignment_sequence) + def testUpdateQual(self): """Ensure SEQ and QUAL updates leading to absent QUAL set all bytes to 0xff""" @@ -273,6 +309,80 @@ class TestAlignedSegment(ReadTest): self.assertEqual(qual, b'\xff' * l_seq) + def testClearQual(self): + a = pysam.AlignedSegment() + a.query_sequence = "ATGC" + a.query_qualities = pysam.qualitystring_to_array("qrst") + a.query_qualities = None + self.assertIsNone(a.query_qualities) + + def testClearQualStr(self): + a = pysam.AlignedSegment() + a.query_sequence = "ATGC" + a.query_qualities_str = "qrst" + self.assertEqual(a.query_qualities, pysam.qualitystring_to_array("qrst")) + self.assertEqual(a.query_qualities_str, "qrst") + + a.query_qualities_str = None + self.assertIsNone(a.query_qualities) + self.assertIsNone(a.query_qualities_str) + + a.query_qualities_str = "qrst" + a.query_qualities_str = "" + self.assertIsNone(a.query_qualities) + self.assertIsNone(a.query_qualities_str) + + a.query_qualities_str = "qrst" + a.query_qualities_str = "*" + self.assertIsNone(a.query_qualities) + self.assertIsNone(a.query_qualities_str) + + def testUpdateQualArrayB(self): + a = pysam.AlignedSegment() + a.query_sequence = "ATGC" + a.query_qualities = array.array('B', [80, 81, 82, 83]) + self.assertEqual(len(a.query_qualities), 4) + self.assertEqual(a.query_qualities_str, "qrst") + + def testUpdateQualArrayI(self): + a = pysam.AlignedSegment() + a.query_sequence = "ATGC" + a.query_qualities = array.array('I', [80, 81, 82, 83]) + self.assertEqual(len(a.query_qualities), 4) + self.assertEqual(a.query_qualities_str, "qrst") + + def testUpdateQualList(self): + a = pysam.AlignedSegment() + a.query_sequence = "ATGC" + qual = [80, 81, 82, 83] + a.query_qualities = qual + qual.pop() + self.assertEqual(len(a.query_qualities), 4) + self.assertEqual(a.query_qualities_str, "qrst") + + def testUpdateQualString(self): + a = pysam.AlignedSegment() + a.query_sequence = "ATGC" + a.query_qualities = "qrst" + self.assertEqual(len(a.query_qualities), 4) + self.assertEqual(a.query_qualities_str, "qrst") + self.assertEqual(a.qual, "qrst") + + def testUpdateQualString2(self): + a = pysam.AlignedSegment() + a.query_sequence = "ATGC" + a.query_qualities_str = "qrst" + self.assertEqual(len(a.query_qualities), 4) + self.assertEqual(a.query_qualities_str, "qrst") + self.assertEqual(a.qual, "qrst") + + def testUpdateQualTuple(self): + a = pysam.AlignedSegment() + a.query_sequence = "ATGC" + a.query_qualities = (80, 81, 82, 83) + self.assertEqual(len(a.query_qualities), 4) + self.assertEqual(a.query_qualities_str, "qrst") + def testLargeRead(self): """build an example read.""" @@ -955,6 +1065,16 @@ class TestCigar(ReadTest): r.cigarstring = None self.assertEqual(r.cigarstring, None) + r.cigarstring = "40M" + self.assertEqual(r.cigartuples, [(0, 40)]) + r.cigarstring = "" + self.assertEqual(r.cigarstring, None) + + r.cigarstring = "40M" + self.assertEqual(r.cigartuples, [(0, 40)]) + r.cigarstring = "*" + self.assertEqual(r.cigarstring, None) + def testCigar(self): r = self.build_read() self.assertEqual(r.cigartuples, [(0, 10), (2, 1), (0, 9), (1, 1), (0, 20)]) @@ -992,6 +1112,16 @@ class TestCigarStats(ReadTest): expected[1][i] = 1 self.assertEqual([list(x) for x in a.get_cigar_stats()], expected) + for i in range(1, 100): + cigarstring = "".join("10{}".format(x) + for x in iter("MIDNSHP=X")) * i + a.cigarstring = cigarstring + self.assertEqual(a.cigarstring, cigarstring) + expected = [[i * 10 for j in range(len("MIDNSHP=X"))] + [0, 0], + [i for j in range(len("MIDNSHP=X"))] + [0, 0]] + obtained = [list(x) for x in a.get_cigar_stats()] + self.assertEqual(obtained, expected) + a.cigarstring = "10M" a.set_tag("NM", 5) self.assertEqual( diff --git a/tests/compile_test.py b/tests/compile_test.py index 2ea3343..2de1b1e 100644 --- a/tests/compile_test.py +++ b/tests/compile_test.py @@ -10,11 +10,12 @@ pysam and tabix works. import os import pytest import pysam -from TestUtils import make_data_files, BAM_DATADIR, TABIX_DATADIR +from TestUtils import make_data_files, BAM_DATADIR, CBCF_DATADIR, TABIX_DATADIR def setUpModule(): make_data_files(BAM_DATADIR) + make_data_files(CBCF_DATADIR) make_data_files(TABIX_DATADIR) @@ -50,3 +51,21 @@ def test_gtf(): nread = _compile_test.testCountGTF( pysam.Tabixfile(input_filename)) assert nread == 237 + + +class TestBinaryCompatibility: + def test_alignments(self): + fp = pysam.AlignmentFile(os.path.join(BAM_DATADIR, "ex1.bam")) + hdr = pysam.AlignmentHeader() + aln = pysam.AlignedSegment() + + assert fp.__sizeof__() == 120 + assert hdr.__sizeof__() == 24 + assert aln.__sizeof__() == 72 + + def test_variants(self): + fp = pysam.VariantFile(os.path.join(CBCF_DATADIR, "example_vcf43.vcf")) + hdr = pysam.VariantHeader() + + assert fp.__sizeof__() == 120 + assert hdr.__sizeof__() == 32 diff --git a/tests/refactoring.pl b/tests/refactoring.pl deleted file mode 100644 index 3300111..0000000 --- a/tests/refactoring.pl +++ /dev/null @@ -1,60 +0,0 @@ -while () { - # Samfile refactoring - s/Samfile/AlignmentFile/g; - - # AlignedRead refactoring - s/AlignedRead/AlignedSegment/g; - - # Do these patterns first as they match - # the new names - s/\.query/\.query_alignment_sequence/g; - s/\.positions/\.getReferencePositions()/g; - - # Tabixfile, etc - s/Tabixfile/TabixFile/g; - s/Fastafile/FastaFile/g; - s/Fastqfile/FastqFile/g; - - # basic attributes - s/\.qname/\.query_name/g; - s/\.tid/\.reference_id/g; - s/\.pos/\.reference_start/g; - s/\.mapq/\.mapping_quality/g; - s/\.rnext/\.next_reference_id/g; - s/\.pnext/\.next_reference_start/g; - s/\.tlen/\.query_length/g; - s/\.seq/\.query_sequence/g; - if (/\.qual =/) { - s/([[\].0-9a-zA-Z]*)\.qual = (\S*)/$1.query_qualities = pysam.fromQualityString($2)/g; - } else { - s/([[\].0-9a-zA-Z]*)\.qual/pysam.toQualityString($1\.query_qualities)/g; - } - s/\.alen/\.reference_length/g; - s/\.aend/\.reference_end/g; - s/\.rlen/\.query_alignment_length/g; - s/([[\].0-9a-zA-Z]*)\.qqual/pysam.toQualityString($1\.query_alignment_qualities)/g; - s/\.qstart/\.query_alignment_start/g; - s/\.qend/\.query_alignment_end/g; - s/\.qlen/\.query_alignment_length/g; - s/\.mrnm/\.next_reference_id/g; - s/\.rnext/\.next_reference_id/g; - s/\.mpos/\.next_reference_start/g; - s/\.rname/\.reference_id/g; - s/\.isize/\.query_length/g; - s/\.cigar/\.cigartuples/g unless (/\.cigarstring/); - - s/\.blocks/\.getBlocks()/g; - s/\.aligned_pairs/\.getAlignedPairs()/g; - s/\.inferred_length/\.getInferredQueryLength()/g; - - s/\.overlap()/\.getOverlap()/g; - - # PileupProxy - s/\.n([^a-zA-Z])/\.nsegments$1/g; - - # if (/\.mrnm/ || /\.rnext/ || /\.mpos/ || /\.rname/) - # { - # warn "Deprecated tag $& at line $.\n"; - # } - print; -} diff --git a/tests/samtools_test.py b/tests/samtools_test.py index 182e919..fc31d8c 100644 --- a/tests/samtools_test.py +++ b/tests/samtools_test.py @@ -291,6 +291,16 @@ class EmptyIndexTest(unittest.TestCase): "exdoesntexist.bam") +class ExerciseSubcommands(unittest.TestCase): + def testFailingSamtools(self): + with self.assertRaises(pysam.SamtoolsError): + pysam.samtools.view("nonexistent.bam") + + def testFailingBCFtools(self): + with self.assertRaises(pysam.SamtoolsError): + pysam.bcftools.view("nonexistent.vcf") + + if sys.platform != "darwin": # fails with segfault with htslib 1.5 on Osx, an issue with flockfile # issue seems to be with repeated calls to interface -- 2.39.5