From aaa192833cb14fb7f637904b4540527f50ddb627 Mon Sep 17 00:00:00 2001 From: Andrea Talenti <23279528+RenzoTale88@users.noreply.github.com> Date: Mon, 2 Feb 2026 15:25:56 +0000 Subject: [PATCH 1/5] Update README --- README.md | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/README.md b/README.md index 318cdeb..c89c9ff 100644 --- a/README.md +++ b/README.md @@ -3,10 +3,6 @@ A rust implementation of the XP-CLR method. This implementation achieves near identical results in a fraction of the run time. The software analyses chromosome 24 of the VarGoat dataset (777,865 total variants, 236,145 used for the analysis, with two groups of 32 and 22 individuals, respectively) in 00m:26s and using 77Mb of memory, versus 55m:20s and 321Mb of the original implementation. -> [!WARNING] -> The software is still under active development, and is not entirely bug free. There are small differences in the results obtained with the original tool. - - ## Installation The compilation of the software requires the following packages to be installed: 1. [openblas](https://github.com/OpenMathLib/OpenBLAS) @@ -14,12 +10,9 @@ The compilation of the software requires the following packages to be installed: 3. [curl](https://curl.se/) 4. [rust](https://rust-lang.org/) -Then, clone and compile the code: +Then, install with cargo: ``` -git clone https://www.github.com/RenzoTale88/xpclrs -pushd xpclrs -cargo install --path . -popd +cargo install xpclrs ``` Check that the package is successfully installed with: @@ -27,6 +20,16 @@ Check that the package is successfully installed with: xpclrs --help ``` +### Install with Docker/Singularity +The software is also available as a docker container in dockerhub. You can install it by pulling the image with docker: +``` +docker pull tale88/xpclrs:latest +``` +Or with singularity: +``` +singularity build xpclrs.sif tale88/xpclrs:latest +``` + ## Input The software requires the following mandatory options: 1. Input genotypes in VCF(.GZ)/BCF format with `-I-`/`--input`. @@ -75,4 +78,3 @@ Can test with the demo data in the original xpclr repository [here](https://gith ## Citation If you use the tool, please cite: > Chen H, Patterson N, Reich D. Population differentiation as a test for selective sweeps. Genome Res. 2010 Mar;20(3):393-402. doi: 10.1101/gr.100545.109. Epub 2010 Jan 19. PMID: 20086244; PMCID: PMC2840981. -> Talenti A. \ No newline at end of file From 2a05a79d897ef72e21d17441ee91d3794fc654fd Mon Sep 17 00:00:00 2001 From: Andrea Talenti <23279528+RenzoTale88@users.noreply.github.com> Date: Mon, 2 Feb 2026 15:29:06 +0000 Subject: [PATCH 2/5] Bop version --- Cargo.lock | 2 +- Cargo.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 69656a3..f80ad6b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2233,7 +2233,7 @@ checksum = "9edde0db4769d2dc68579893f2306b26c6ecfbe0ef499b013d731b7b9247e0b9" [[package]] name = "xpclrs" -version = "1.0.0" +version = "1.0.1" dependencies = [ "anyhow", "clap", diff --git a/Cargo.toml b/Cargo.toml index af99069..5ce69be 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "xpclrs" -version = "1.0.0" +version = "1.0.1" authors = ["Andrea Talenti"] edition = "2021" license = "MIT" From 246b3358216cb35635c1405b43a17b45b6c9e5b5 Mon Sep 17 00:00:00 2001 From: Andrea Talenti <23279528+RenzoTale88@users.noreply.github.com> Date: Mon, 2 Feb 2026 16:20:30 +0000 Subject: [PATCH 3/5] Load target markers in place, rather than all of them --- src/plink/mod.rs | 65 +++++++++++++++++++++++++++++------------------- 1 file changed, 39 insertions(+), 26 deletions(-) diff --git a/src/plink/mod.rs b/src/plink/mod.rs index 4e5d5a4..0100f91 100644 --- a/src/plink/mod.rs +++ b/src/plink/mod.rs @@ -15,9 +15,14 @@ pub struct PlinkData { } // Primary BED mode ==> SNP major, i.e. codify by individual for SNP -fn bed_snp_major(bed_file: &mut File, n_snps: usize, n_samples: usize) -> Vec> { +fn bed_snp_major(bed_file: &mut File, n_samples: usize, kept_sites: &[bool]) -> Vec> { // Initialize genotype matrix S X I - let mut gts = vec![vec![0u8; n_samples]; n_snps]; + let n_kept_snps = kept_sites + .iter() + .filter(|&&v| v) + .collect::>() + .len(); + let mut gts = vec![vec![0u8; n_samples]; n_kept_snps]; // Read genotypes let bytes_per_snp = n_samples.div_ceil(4); // Ensure we have a full byte every time @@ -25,18 +30,22 @@ fn bed_snp_major(bed_file: &mut File, n_snps: usize, n_samples: usize) -> Vec> shift) & 0b11; + if *keep_site { + for (sample_idx, slot) in gts[kept_snp_idx].iter_mut().enumerate() { + let byte_index = sample_idx / 4; // Which byte (0, 1, 2, ...) + let shift = 2 * (sample_idx % 4); // Shift to get the right pair + let bits = (snp_bytes[byte_index] >> shift) & 0b11; - *slot = bits; + *slot = bits; + } + kept_snp_idx += 1; } + } log::info!("Genotypes loaded."); log::info!("Matrix shape: {} x {}.", gts.len(), gts[0].len()); @@ -44,26 +53,34 @@ fn bed_snp_major(bed_file: &mut File, n_snps: usize, n_samples: usize) -> Vec individual major, i.e. codify by SNPs for individual -fn bed_ind_major(bed_file: &mut File, n_snps: usize, n_samples: usize) -> Vec> { +fn bed_ind_major(bed_file: &mut File, n_samples: usize, kept_sites: &[bool]) -> Vec> { // Initialize genotype matrix I X S - let mut gts = vec![vec![0u8; n_samples]; n_snps]; + let n_kept_snps = kept_sites + .iter() + .filter(|&&v| v) + .collect::>() + .len(); + let mut gts = vec![vec![0u8; n_samples]; n_kept_snps]; // Read genotypes - let bytes_per_ind = n_snps.div_ceil(4); + let bytes_per_ind = n_kept_snps.div_ceil(4); let mut ind_bytes = vec![0u8; bytes_per_ind]; // Load n_snps at the same time, representing one individual - for bit_row in gts.iter_mut().take(n_snps) { + for sample_idx in 0..gts[0].len() { bed_file .read_exact(&mut ind_bytes) .expect("Cannot read bytes"); + let mut kept_snp_idx = 0; + for (snp_idx, keep_site) in kept_sites.iter().enumerate() { + if *keep_site { + let byte_index = snp_idx / 4; // Which byte (0, 1, 2, ...) + let shift = 2 * (snp_idx % 4); // Shift to get the right pair + let bits = (ind_bytes[byte_index] >> shift) & 0b11; - for (snp_idx, slot) in bit_row.iter_mut().enumerate() { - let byte_index = snp_idx / 4; // Which byte (0, 1, 2, ...) - let shift = 2 * (snp_idx % 4); // Shift to get the right pair - let bits = (ind_bytes[byte_index] >> shift) & 0b11; - - *slot = bits; + gts[kept_snp_idx][sample_idx] = bits; + kept_snp_idx += 1; + } } } log::info!("Genotypes loaded."); @@ -183,15 +200,11 @@ pub fn read_plink_files( // Define if it is SNP major or individual major let genotypes: Vec> = if magic[2] == 0x01 { log::info!("Reading in SNP-major mode."); - bed_snp_major(&mut bed_file, n_snps, n_samples) + bed_snp_major(&mut bed_file, n_samples, &keep_vec) } else { log::info!("Reading in Individual-major mode."); - bed_ind_major(&mut bed_file, n_snps, n_samples) - } - .into_iter() - .zip(keep_vec) - .filter_map(|(gt_row, keep)| if keep { Some(gt_row) } else { None }) - .collect(); + bed_ind_major(&mut bed_file, n_samples, &keep_vec) + }; // Filter the dataset let mut monom_gt2 = 0; From 6154969f16626a154bf951135b84e660e6bcf22c Mon Sep 17 00:00:00 2001 From: Andrea Talenti <23279528+RenzoTale88@users.noreply.github.com> Date: Tue, 3 Feb 2026 08:32:55 +0000 Subject: [PATCH 4/5] Focused importing of samples --- src/plink/mod.rs | 81 ++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 65 insertions(+), 16 deletions(-) diff --git a/src/plink/mod.rs b/src/plink/mod.rs index 0100f91..34f99ac 100644 --- a/src/plink/mod.rs +++ b/src/plink/mod.rs @@ -4,6 +4,7 @@ This module provides the binary plink-specific I/O functions. use crate::io::{consolidate_list, get_gt_index, GenoData}; use itertools::izip; use itertools::MultiUnzip; +use rustc_hash::FxHashMap; use std::fs::File; use std::io::{self, BufRead, BufReader, Read}; @@ -15,17 +16,23 @@ pub struct PlinkData { } // Primary BED mode ==> SNP major, i.e. codify by individual for SNP -fn bed_snp_major(bed_file: &mut File, n_samples: usize, kept_sites: &[bool]) -> Vec> { +fn bed_snp_major( + bed_file: &mut File, + kept_samples: &[bool], + samples_mat_idx: FxHashMap, + kept_sites: &[bool], +) -> Vec> { // Initialize genotype matrix S X I let n_kept_snps = kept_sites .iter() .filter(|&&v| v) .collect::>() .len(); - let mut gts = vec![vec![0u8; n_samples]; n_kept_snps]; + let n_kept_samples = samples_mat_idx.len(); + let mut gts = vec![vec![0u8; n_kept_samples]; n_kept_snps]; // Read genotypes - let bytes_per_snp = n_samples.div_ceil(4); // Ensure we have a full byte every time + let bytes_per_snp = kept_samples.len().div_ceil(4); // Ensure we have a full byte every time let mut snp_bytes = vec![0u8; bytes_per_snp]; // Initiate progress bar @@ -36,16 +43,19 @@ fn bed_snp_major(bed_file: &mut File, n_samples: usize, kept_sites: &[bool]) -> .read_exact(&mut snp_bytes) .expect("Cannot read bytes"); if *keep_site { - for (sample_idx, slot) in gts[kept_snp_idx].iter_mut().enumerate() { + for (sample_idx, kept_sample) in kept_samples.iter().enumerate() { + if !*kept_sample { + continue; + } + let sample_idx_in_mat = samples_mat_idx.get(&sample_idx).expect("Missing index"); let byte_index = sample_idx / 4; // Which byte (0, 1, 2, ...) let shift = 2 * (sample_idx % 4); // Shift to get the right pair let bits = (snp_bytes[byte_index] >> shift) & 0b11; - *slot = bits; + gts[kept_snp_idx][*sample_idx_in_mat] = bits; } kept_snp_idx += 1; } - } log::info!("Genotypes loaded."); log::info!("Matrix shape: {} x {}.", gts.len(), gts[0].len()); @@ -53,24 +63,34 @@ fn bed_snp_major(bed_file: &mut File, n_samples: usize, kept_sites: &[bool]) -> } // Other BED mode ==> individual major, i.e. codify by SNPs for individual -fn bed_ind_major(bed_file: &mut File, n_samples: usize, kept_sites: &[bool]) -> Vec> { +fn bed_ind_major( + bed_file: &mut File, + kept_samples: &[bool], + samples_mat_idx: FxHashMap, + kept_sites: &[bool], +) -> Vec> { // Initialize genotype matrix I X S let n_kept_snps = kept_sites .iter() .filter(|&&v| v) .collect::>() .len(); - let mut gts = vec![vec![0u8; n_samples]; n_kept_snps]; + let n_kept_samples = samples_mat_idx.len(); + let mut gts = vec![vec![0u8; n_kept_samples]; n_kept_snps]; // Read genotypes let bytes_per_ind = n_kept_snps.div_ceil(4); let mut ind_bytes = vec![0u8; bytes_per_ind]; // Load n_snps at the same time, representing one individual - for sample_idx in 0..gts[0].len() { + for (sample_idx, kept_sample) in kept_samples.iter().enumerate() { bed_file .read_exact(&mut ind_bytes) .expect("Cannot read bytes"); + if !*kept_sample { + continue; + } + let sample_idx_in_mat = samples_mat_idx.get(&sample_idx).expect("Missing index"); let mut kept_snp_idx = 0; for (snp_idx, keep_site) in kept_sites.iter().enumerate() { if *keep_site { @@ -78,7 +98,7 @@ fn bed_ind_major(bed_file: &mut File, n_samples: usize, kept_sites: &[bool]) -> let shift = 2 * (snp_idx % 4); // Shift to get the right pair let bits = (ind_bytes[byte_index] >> shift) & 0b11; - gts[kept_snp_idx][sample_idx] = bits; + gts[kept_snp_idx][*sample_idx_in_mat] = bits; kept_snp_idx += 1; } } @@ -125,18 +145,47 @@ pub fn read_plink_files( }) .collect(); let sample_ids: Vec<&[u8]> = sample_ids_owned.iter().map(|id| id.as_slice()).collect(); - // Number of samples for later use. - let n_samples = sample_ids.len(); // Load sample lists as an u8 array let s1 = consolidate_list(&sample_ids, s1).expect("Failed to subset sampleA"); let s2 = consolidate_list(&sample_ids, s2).expect("Failed to subset sampleB"); // Fetch the indices of each sample in each list - let i1 = get_gt_index(&sample_ids, &s1).expect("Failed to get indeces of sampleA"); - let i2 = get_gt_index(&sample_ids, &s2).expect("Failed to get indeces of sampleB"); + let i1_o = get_gt_index(&sample_ids, &s1).expect("Failed to get indeces of sampleA"); + let i2_o = get_gt_index(&sample_ids, &s2).expect("Failed to get indeces of sampleB"); + // Create a vector of indexes of samples to retain, and a + // Second vector of unique samples to keep + let mut individuals_ixs: Vec = [&i1_o[..], &i2_o[..]].concat(); + individuals_ixs.sort(); + individuals_ixs.dedup(); + let new_indexes = individuals_ixs + .iter() + .enumerate() + .map(|(e, v)| (*v, e)) + .collect::>(); + // Create new group indexes + let i1 = i1_o + .iter() + .map(|&idx| new_indexes.get(&idx).expect("Missing index")) + .copied() + .collect::>(); + let i2 = i2_o + .iter() + .map(|&idx| new_indexes.get(&idx).expect("Missing index")) + .copied() + .collect::>(); + // Generate a vector of retained samples + let kept_samples: Vec = sample_ids + .iter() + .enumerate() + .map(|(idx, _)| individuals_ixs.contains(&idx)) + .collect(); // Print number of samples log::info!("Samples A: {}", i1.len()); log::info!("Samples B: {}", i2.len()); + log::debug!("Old samples indexes A: {:?}", i1_o); + log::debug!("New samples indexes A: {:?}", i1); + log::debug!("Old samples indexes B: {:?}", i2_o); + log::debug!("New samples indexes B: {:?}", i2); // Dies if no samples are retained if s1.is_empty() || s2.is_empty() { @@ -200,10 +249,10 @@ pub fn read_plink_files( // Define if it is SNP major or individual major let genotypes: Vec> = if magic[2] == 0x01 { log::info!("Reading in SNP-major mode."); - bed_snp_major(&mut bed_file, n_samples, &keep_vec) + bed_snp_major(&mut bed_file, &kept_samples, new_indexes, &keep_vec) } else { log::info!("Reading in Individual-major mode."); - bed_ind_major(&mut bed_file, n_samples, &keep_vec) + bed_ind_major(&mut bed_file, &kept_samples, new_indexes, &keep_vec) }; // Filter the dataset From b5198bed17ed8ffde2d1e7ec88f7dca1bd28e668 Mon Sep 17 00:00:00 2001 From: Andrea Talenti <23279528+RenzoTale88@users.noreply.github.com> Date: Mon, 2 Mar 2026 10:31:47 +0000 Subject: [PATCH 5/5] Update citations --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index c89c9ff..49eb776 100644 --- a/README.md +++ b/README.md @@ -78,3 +78,5 @@ Can test with the demo data in the original xpclr repository [here](https://gith ## Citation If you use the tool, please cite: > Chen H, Patterson N, Reich D. Population differentiation as a test for selective sweeps. Genome Res. 2010 Mar;20(3):393-402. doi: 10.1101/gr.100545.109. Epub 2010 Jan 19. PMID: 20086244; PMCID: PMC2840981. +> The original xpclr tool [here](https://github.com/hardingnj/xpclr) +> Talenti A. XPCLRS: Fast Selection Signature Detection Using Cross-Population Composite Likelihood Ratio. bioRxiv 2026.02.27.708459. doi: 10.64898/2026.02.27.708459.