Lines
93.42 %
Functions
7.5 %
Branches
100 %
use std::path::PathBuf;
use color_eyre::{eyre::eyre, Result};
use crate::args::{Selection, StandardArgs};
use crate::read_vcf::read_vcf_to_matrix;
use crate::structs::{HapVariant, PhasedMatrix};
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)]
pub fn run(args: StandardArgs, haplotype_path: PathBuf) -> Result<()> {
let (contig, variant_pos) = parse_snp_coord(&args.coords)?;
let mut csv_output = args.output.clone();
push_to_output(&args, &mut csv_output, "haplotype_check", "csv");
let mut writer = open_csv_writer(csv_output)?;
let ht = crate::io::read_haplotype_file(haplotype_path)?;
let start = ht.first().unwrap();
let end = ht.last().unwrap();
tracing::debug!("Reading into a sample-variant matrix");
let vcf = match args.selection {
Selection::OnlyAlts | Selection::OnlyRefs => {
let vcf = read_vcf_to_matrix(
&args,
contig,
variant_pos,
Some((start.pos as u64, end.pos as u64)),
)?;
select_carrier_haplotypes(vcf, variant_pos, &args.coords, &args.selection)
}
Selection::OnlyLongest => {
let vcf = read_vcf_to_matrix(&args, contig, variant_pos, None)?;
let shared_lengths = shared_lengths_by_majority(&vcf, vcf.variant_idx());
select_only_longest_haplotypes(&shared_lengths, vcf)
Selection::All => read_vcf_to_matrix(
)?,
Selection::Unphased => return Err(eyre!("Running with unphased data is not supported."))
tracing::debug!("Finished reading the matrix");
let matching_indexes = identical_haplotype_count(&vcf, &ht);
write_matches_to_csv(&matching_indexes, &mut writer, &vcf)?;
tracing::debug!(
"Haplotype found in {}/{} of alleles.",
matching_indexes.len(),
vcf.samples().len()
);
Ok(())
fn write_matches_to_csv(
matching_indexes: &[usize],
writer: &mut csv::Writer<Box<dyn std::io::Write>>,
vcf: &PhasedMatrix,
) -> Result<()> {
writer.write_record(vec!["id", "match"])?;
for (idx, sample) in vcf.samples().iter().enumerate() {
let allele = if idx >= vcf.nsamples() / 2 { 2 } else { 1 };
writer.write_record(vec![
format!("{sample}_{allele}"),
format!("{}", matching_indexes.contains(&idx)),
])?;
pub fn identical_haplotype_count(vcf: &PhasedMatrix, ht: &[HapVariant]) -> Vec<usize> {
let indexes: Vec<(usize, usize)> = ht
.iter()
.enumerate()
.filter_map(|(coords_i, h)| {
// Check the PartialEq implementation between HapVariant and Coord
if let Some(vcf_i) = vcf.coords().iter().position(|c| c == h) {
Some((vcf_i, coords_i))
} else {
tracing::warn!("Haplotype variant {:?} not found in the vcf", h);
None
})
.collect();
if indexes.is_empty() {
return vec![];
let mut matching = vec![true; vcf.matrix.nrows()];
for (n, value) in matching.iter_mut().enumerate() {
'inner: for (vcf_i, coords_i) in &indexes {
if vcf.matrix[[n, *vcf_i]] != ht[*coords_i].gt {
*value = false;
break 'inner;
matching
.filter_map(|(i, m)| if *m { Some(i) } else { None })
.collect()