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" diff --git a/README.md b/README.md index 318cdeb..49eb776 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,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. -> Talenti A. \ No newline at end of file +> 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. diff --git a/src/plink/mod.rs b/src/plink/mod.rs index 4e5d5a4..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,27 +16,45 @@ 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, + kept_samples: &[bool], + samples_mat_idx: FxHashMap, + 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 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 log::info!("Load genotypes..."); - - for snp_row in gts.iter_mut() { + let mut kept_snp_idx = 0; + for keep_site in kept_sites.iter() { bed_file .read_exact(&mut snp_bytes) .expect("Cannot read bytes"); - for (sample_idx, slot) in snp_row.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; + if *keep_site { + 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."); @@ -44,26 +63,44 @@ 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, + kept_samples: &[bool], + samples_mat_idx: FxHashMap, + 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 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_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, 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 { + 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_in_mat] = bits; + kept_snp_idx += 1; + } } } log::info!("Genotypes loaded."); @@ -108,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() { @@ -183,15 +249,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, &kept_samples, new_indexes, &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, &kept_samples, new_indexes, &keep_vec) + }; // Filter the dataset let mut monom_gt2 = 0;