Lines
80.52 %
Functions
17.86 %
Branches
100 %
use color_eyre::Result;
use ndarray::{s, Array1, Array2, Axis};
use ndarray_rand::{RandomExt, SamplingStrategy};
use rayon::prelude::*;
use crate::args::{Selection, StandardArgs};
use crate::read_vcf::{read_vcf_to_contra_matrix, read_vcf_to_matrix};
use crate::utils::{
select_carrier_haplotypes, select_only_longest_haplotypes, shared_lengths_by_majority, push_to_output,
};
use crate::core::{open_csv_writer, parse_snp_coord};
#[doc(hidden)]
#[tracing::instrument]
pub fn run(
args: StandardArgs,
niters: i64,
sample_size: usize,
) -> Result<()> {
let mut output = args.output.clone();
push_to_output(&args, &mut output, "random_sampled_shared_haplotype", "csv");
let mut writer = open_csv_writer(output)?;
let (contig, variant_pos) = parse_snp_coord(&args.coords)?;
let (variant_matrix, pos_vec) = match args.selection {
Selection::OnlyAlts | Selection::OnlyRefs => {
let vcf = read_vcf_to_matrix(&args, contig, variant_pos, None)?;
let vcf = select_carrier_haplotypes(vcf, variant_pos, &args.coords, &args.selection);
let pos_vec = vcf.coords().iter().map(|c| c.pos).collect::<Vec<i64>>();
(vcf.matrix, pos_vec)
}
Selection::OnlyLongest => {
let mut vcf = read_vcf_to_matrix(&args, contig, variant_pos, None)?;
let shared_lengths = shared_lengths_by_majority(&vcf, vcf.variant_idx());
vcf = select_only_longest_haplotypes(&shared_lengths, vcf);
Selection::All => {
Selection::Unphased => read_vcf_to_contra_matrix(&args, contig, None)?
tracing::debug!("Finished reading the matrix with selection {:?}.", args.selection);
writer.write_record(vec!["start", "snp_count", "bp_length"])?;
// Multithreaded random sampling of the variant matrix
let lengths = (0..niters)
.into_par_iter()
.flat_map(|_| {
let sampled_matrix = variant_matrix.sample_axis(
Axis(0),
sample_size,
SamplingStrategy::WithoutReplacement,
);
let shared_sequence =
shared_sequences(&sampled_matrix, args.selection != Selection::Unphased);
find_shared_sequence_lengths(shared_sequence, &pos_vec, variant_pos)
})
.collect::<Vec<(i64, i64, i64)>>();
tracing::debug!("Finished processing lengths.");
for (start, snp_count, bp_length) in lengths {
writer.write_record(vec![
start.to_string(),
snp_count.to_string(),
bp_length.to_string(),
])?;
Ok(())
// Create a 1-D array from 2-D a genotype matrix that denotes all the positions where the samples
// have a contrahomozygote in between them. A 0 denotes no contrahomozygotes are present and 1
// denotes the opposite
//
pub fn shared_sequences(matrix: &Array2<u8>, phased: bool) -> Array1<u8> {
let mut shared_sequence = Array1::<u8>::from_elem(matrix.ncols(), 0);
// Flatten a 2-D matrix into an array. If both 0 and 2 exist denote as 1 else remain as 0
for i in 0..matrix.ncols() {
let column = matrix.slice(ndarray::s![.., i]);
if phased {
if column.into_iter().any(|v| *v == 0) && column.into_iter().any(|v| *v == 1) {
shared_sequence.slice_mut(s![i]).fill(1);
} else if column.into_iter().any(|v| *v == 0) && column.into_iter().any(|v| *v == 2) {
shared_sequence
pub fn find_shared_sequence_lengths(
shared_sequence: Array1<u8>,
pos_coords: &[i64],
variant_pos: i64,
) -> Option<(i64, i64, i64)> {
let mut shared_lengths_and_pos = None;
let mut counter = 1;
.iter()
.enumerate()
.fold(1, |last, (index, curr)| {
if last == 0 && *curr == 0 {
counter += 1;
} else if last == 0 && *curr == 1 {
let start = pos_coords.get(index - counter).unwrap();
let stop = pos_coords.get(index).unwrap();
// Find the sequence that overlaps the targeted variant
if *start < variant_pos && *stop > variant_pos {
shared_lengths_and_pos = Some((*start, counter as i64, stop - start));
counter = 1;
*curr
});
counter -= 1;
if counter > 0 {
let index = shared_sequence.len() - 1;
shared_lengths_and_pos