From 152d8f1e96db336a80fcbfabb475107f7e6937af Mon Sep 17 00:00:00 2001 From: "Kevin R. Thornton" Date: Fri, 2 Oct 2026 12:39:21 -0700 Subject: [PATCH] feat: popgen-htslib crate --- .github/workflows/test.yml | 6 + Cargo.lock | 24 ++- Cargo.toml | 1 + popgen-htslib/Cargo.toml | 18 ++ .../examples/from_vcf_htslib.rs | 0 popgen-htslib/src/lib.rs | 188 ++++++++++++++++++ popgen-htslib/tests/htslib_vcf.rs | 54 +++++ popgen/Cargo.toml | 4 - 8 files changed, 283 insertions(+), 12 deletions(-) create mode 100644 popgen-htslib/Cargo.toml rename {popgen => popgen-htslib}/examples/from_vcf_htslib.rs (100%) create mode 100644 popgen-htslib/src/lib.rs create mode 100644 popgen-htslib/tests/htslib_vcf.rs diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 1b8bf06f..9ecfeed6 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -61,6 +61,12 @@ jobs: cargo hack test --all-targets --target=i686-unknown-linux-gnu --feature-powerset --workspace --exclude python-integration-tests cargo hack test --doc --target=i686-unknown-linux-gnu --feature-powerset --workspace --exclude python-integration-tests if: matrix.os == 'ubuntu-24.04' + - name: valgrind (crates using unsafe) + run: | + sudo apt-get install -y valgrind + cargo install cargo-valgrind + cargo valgrind test --manifest-path popgen-htslib/Cargo.toml + if: matrix.os == 'ubuntu-24.04' fmt: name: rustfmt diff --git a/Cargo.lock b/Cargo.lock index 5a352c3e..37b2b955 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -568,9 +568,9 @@ checksum = "d26c52dbd32dccf2d10cac7725f8eae5296885fb5703b261f7d0a0739ec807ab" [[package]] name = "linux-raw-sys" -version = "0.9.4" +version = "0.12.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cd945864f07fe9f5371a27ad7b52a172b4b499999f1d97574c9fa68373937e12" +checksum = "32a66949e030da00e8c7d4434b251670a91556f4144941d37452769c25d58a53" [[package]] name = "litemap" @@ -733,7 +733,15 @@ dependencies = [ "proptest", "rand", "rand_distr", +] + +[[package]] +name = "popgen-htslib" +version = "0.10.0-alpha.0" +dependencies = [ + "popgen", "rust-htslib", + "tempfile", ] [[package]] @@ -1051,14 +1059,14 @@ dependencies = [ [[package]] name = "rustix" -version = "1.0.8" +version = "1.1.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "11181fbabf243db407ef8df94a6ce0b2f9a733bd8be4ad02b4eda9602296cac8" +checksum = "891efababe418670775f199f0d233d84843c227a0949a883ce15b37c78d6629d" dependencies = [ "bitflags", "errno", "libc", - "linux-raw-sys 0.9.4", + "linux-raw-sys 0.12.1", "windows-sys", ] @@ -1173,14 +1181,14 @@ checksum = "adb6935a6f5c20170eeceb1a3835a49e12e19d792f6dd344ccc76a985ca5a6ca" [[package]] name = "tempfile" -version = "3.23.0" +version = "3.27.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2d31c77bdf42a745371d260a26ca7163f1e0924b64afa0b688e61b5a9fa02f16" +checksum = "32497e9a4c7b38532efcdebeef879707aa9f794296a4f0244f6f69e9bc8574bd" dependencies = [ "fastrand", "getrandom", "once_cell", - "rustix 1.0.8", + "rustix 1.1.5", "windows-sys", ] diff --git a/Cargo.toml b/Cargo.toml index 00b1c601..1805aa91 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,6 +3,7 @@ resolver = "3" members = [ "popgen", + "popgen-htslib", "popgen-noodles", "popgen-tskit", "python-integration-tests" diff --git a/popgen-htslib/Cargo.toml b/popgen-htslib/Cargo.toml new file mode 100644 index 00000000..484983e4 --- /dev/null +++ b/popgen-htslib/Cargo.toml @@ -0,0 +1,18 @@ +[package] +name = "popgen-htslib" +version = "0.10.0-alpha.0" +edition.workspace = true +rust-version.workspace = true + +[lints] +workspace = true + +[dependencies] +popgen = { version = "0.10.0-alpha.0", path = "../popgen" } +rust-htslib = { version = "1.0.0", default-features = false } + +[dev-dependencies] +tempfile = "3.27.0" + +[[example]] +name = "from_vcf_htslib" diff --git a/popgen/examples/from_vcf_htslib.rs b/popgen-htslib/examples/from_vcf_htslib.rs similarity index 100% rename from popgen/examples/from_vcf_htslib.rs rename to popgen-htslib/examples/from_vcf_htslib.rs diff --git a/popgen-htslib/src/lib.rs b/popgen-htslib/src/lib.rs new file mode 100644 index 00000000..97e2c721 --- /dev/null +++ b/popgen-htslib/src/lib.rs @@ -0,0 +1,188 @@ +//! Adapter types for [`rust_htslib`] + +use popgen::AlleleID; +use std::ffi::c_void; + +/// Error type +#[non_exhaustive] +#[derive(Debug)] +pub enum Error { + // NOTE: this is a bad name... + /// Encapsulation of errors from [`rust_htslib`] + RecordError(rust_htslib::errors::Error), + /// Integer error codes from the htslib C API + ErrorCode(i32), +} + +impl std::fmt::Display for Error { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::RecordError(e) => write!(f, "{e}"), + Self::ErrorCode(code) => write!(f, "error code: {code}"), + } + } +} + +impl std::error::Error for Error {} + +impl From for Error { + fn from(value: rust_htslib::errors::Error) -> Self { + Self::RecordError(value) + } +} + +struct GenotypesAdapter(T, *mut c_void); + +impl Iterator for GenotypesAdapter +where + T: Iterator, +{ + type Item = I; + + fn next(&mut self) -> Option { + self.0.next() + } +} + +impl DoubleEndedIterator for GenotypesAdapter +where + T: DoubleEndedIterator, +{ + fn next_back(&mut self) -> Option { + self.0.next_back() + } +} + +impl ExactSizeIterator for GenotypesAdapter +where + T: ExactSizeIterator, +{ + fn len(&self) -> usize { + self.0.len() + } +} + +impl Drop for GenotypesAdapter { + fn drop(&mut self) { + // Safety: we control construction of this type, and we only construct it using a buffer allocated by htslib C. + unsafe { + rust_htslib::htslib::free(self.1); + } + } +} + +/// Iterator over genotypes in a [`Record`](rust_htslib::bcf::Record). +/// +/// The iterator emits iterators over [`Option`] of [`popgen::AlleleID`] in each genotype. +/// The [`Option::None`] variant implies missing data. +/// +/// +/// Use [`Iterator::flatten`] to convert the return value into an iterator over all the +/// alleles called in the record. +pub fn bcf_record_to_genotypes_iter_adapter( + record: &rust_htslib::bcf::Record, +) -> Result>> + '_, Error> { + // NOTE: this implementation does not rely on the rust_htslib safe API + // because the iterator types defined there cannot be properly flattened/aggregated. + // Basically, the borrow checker prevents this API from being written. + // Therefore, we work at the level of htslib C types. + // Here, the borrow checker correctly notes that the iterator is tied to the + // lifetime of the borrowed (rust-side) Record. + use rust_htslib::htslib; + + // The dst and ndst parameters to bcf_get_format_values and related functions require a buffer with provided pointer and length. + // This is malloc'd/realloc'd within htslib, but we need to hold onto the pointer ourselves. + + let mut gt_buf = (std::ptr::null_mut(), 0); + // SAFETY: we are not using this pointer after the header is dropped. + // (see docs of the header.as_ptr fn in the unsafe block below) + assert!(!unsafe { record.header().as_ptr() }.is_null()); + + let n_alleles_in_record = { + // See docs at https://github.com/samtools/htslib/blob/7c5e3e7ebcdf90c8f96afd8a06d75ffa5603e417/htslib/vcf.h#L1135 + // We ask htslib to parse the GT field. + // Within htslib, #define bcf_get_format_int32 is provided to omit the last argument. + let format_values = unsafe { + htslib::bcf_get_format_values( + // The header struct bcf_hdr_t. + record.header().as_ptr(), + // The line/record bcf1_t. + record.inner, + // We want the field tagged "GT". + c"GT".as_ptr(), + // We pass the buffer from earlier. + &mut gt_buf.0, + &mut gt_buf.1, + // We want this parsed as a collection of integers, where an integer is one of the VCF datatypes. + htslib::BCF_HT_INT as i32, + ) + }; + + // The return value is either a negative error code or the number of values written. + match format_values { + -1 => { + return Err(rust_htslib::errors::Error::BcfUndefinedTag { + tag: String::from("GT"), + } + .into()) + } + -2 => { + return Err(rust_htslib::errors::Error::BcfUnexpectedType { + tag: String::from("GT"), + } + .into()) + } + -3 => { + return Err(rust_htslib::errors::Error::BcfMissingTag { + tag: String::from("GT"), + record: record.desc(), + } + .into()) + } + -4 => return Err(rust_htslib::errors::Error::BcfAllocationError.into()), + other if other < 0 => return Err(Error::ErrorCode(other)), + ret => ret, + } + }; + + // https://github.com/samtools/htslib/blob/7c5e3e7ebcdf90c8f96afd8a06d75ffa5603e417/htslib/vcf.h#L1049 + // https://github.com/samtools/htslib/blob/7c5e3e7ebcdf90c8f96afd8a06d75ffa5603e417/htslib/vcf.h#L166 + // We need field n on the struct bcf_fmt_t, describing the maximum number of alleles per sample. + let fmt_inner = + unsafe { htslib::bcf_get_fmt(record.header().as_ptr(), record.inner, c"GT".as_ptr()) }; + assert!(!fmt_inner.is_null()); + + assert!(!gt_buf.0.is_null()); + + // Safety: We called for the parsing of i32 from htslib, and if the pointer is not null, it points to malloc'd memory. + // As long as htslib correctly returns n_alleles_in_record, this slice is valid. + let gt_iter_iter = unsafe { + std::slice::from_raw_parts( + gt_buf.0.cast_const().cast::(), + n_alleles_in_record as usize, + ) + } + // See https://github.com/samtools/htslib/blob/7c5e3e7ebcdf90c8f96afd8a06d75ffa5603e417/vcf.c#L6157 + // If fmt_inner is not null, then it points to a newly or previously initialized bcf_fmt_t. + .chunks(unsafe { fmt_inner.as_ref().unwrap() }.n as usize) + .map(|s| { + // The magic number below is + // rust_htslib::bcf::record::VECTOR_END_INTEGER, + // which is private. + // It indicates the position where data ends for this sample. + s.split(|v| *v == i32::MIN + 1) + .next() + .unwrap() + .iter() + .map(|a| { + // This replicates code in rust_htslib::bcf::record::GenotypeAllele. + if a > &0 { + Some(AlleleID::from(((*a >> 1) - 1) as usize)) + } else { + None + } + }) + }); + + Ok(GenotypesAdapter(gt_iter_iter, gt_buf.0)) +} diff --git a/popgen-htslib/tests/htslib_vcf.rs b/popgen-htslib/tests/htslib_vcf.rs new file mode 100644 index 00000000..faf50db9 --- /dev/null +++ b/popgen-htslib/tests/htslib_vcf.rs @@ -0,0 +1,54 @@ +//! Basic integration tests + +use rust_htslib::bcf; +use rust_htslib::bcf::Read; + +#[test] +fn test_basic_vcf_input_iter() { + use std::io::Write; + + static VCF_FILE: &str = r#"##fileformat=VCFv4.6 +##FORMAT= +##contig= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT s0 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 +chr0 1 . A C . . . GT 0/0 0/1 0/1 0/0 0/1 0/0 0/1 0/0 0/0 0/0 0/0 1/1 1/0 0/1 0/0 ./. 0/1 0/0 +chr0 1 . G A . . . GT 0/0 ./. 0/1 0/0 0/1 0/1 0/0 0/0 0/0 0/0 0/0 0/0 0/1 ./. 0/1 0/1 0/0 0/0"#; + let mut tfile = tempfile::NamedTempFile::new().unwrap(); + tfile.write_all(VCF_FILE.as_bytes()).unwrap(); + let (_, tfile_path) = tfile.into_parts(); + let mut bcf = bcf::Reader::from_path(tfile_path.as_os_str()).expect("Error opening file."); + let mut counts = popgen::SampleAlleleCounts::of_empty_sample_sets(1); + for record_result in bcf.records() { + let record = record_result.unwrap(); + let allele_id_iter = popgen_htslib::bcf_record_to_genotypes_iter_adapter(&record) + .unwrap() + .flatten(); + counts.add_site(allele_id_iter).unwrap(); + } + assert_eq!( + counts + .iter_sample_set(0) + .unwrap() + .filter(|c| c.total_alleles() == 36) + .count(), + 2 + ); + assert_eq!(counts.num_sites(), 2); + let num_non_missing = counts + .iter_sample_set(0) + .unwrap() + .take(1) + .map(|a| a.counts().iter().sum::()) + .collect::>()[0]; + assert_eq!(num_non_missing, 2 * 18 - 2); + let num_non_missing = counts + .iter_sample_set(0) + .unwrap() + .skip(1) + .map(|a| a.counts().iter().sum::()) + .collect::>()[0]; + assert_eq!(num_non_missing, 2 * 18 - 4); + for i in counts.iter_sample_set(0).unwrap() { + assert_eq!(i.counts().len(), 2) + } +} diff --git a/popgen/Cargo.toml b/popgen/Cargo.toml index 77851245..40785080 100644 --- a/popgen/Cargo.toml +++ b/popgen/Cargo.toml @@ -14,10 +14,6 @@ repository = "https://github.com/ThorntonLab/popgen-oxide" rand = { workspace = true } rand_distr = { workspace = true } proptest = { workspace = true } -rust-htslib = { version = "1.0.0", default-features = false } [lints] workspace = true - -[[example]] -name = "from_vcf_htslib"