Lines
73.91 %
Functions
12.5 %
Branches
100 %
use std::ops::Range;
use color_eyre::Result;
use ndarray::{s, ArrayBase, Axis, Dim, ViewRepr};
use crate::args::{Selection, StandardArgs};
use crate::core::{open_csv_writer, parse_snp_coord};
use crate::read_vcf::{read_vcf_to_matrix, read_vcf_to_matrix_missing_allowed};
use crate::structs::PhasedMatrix;
use crate::utils::{
select_carrier_haplotypes, select_only_longest_haplotypes, shared_lengths_by_majority, push_to_output,
};
#[doc(hidden)]
#[tracing::instrument]
pub fn run(args: StandardArgs, missing_allowed: bool) -> Result<()> {
let mut output = args.output.clone();
push_to_output(&args, &mut output, "shared_haplotype", "csv");
let writer = open_csv_writer(output)?;
let (contig, variant_pos) = parse_snp_coord(&args.coords)?;
let mut vcf = if missing_allowed {
read_vcf_to_matrix_missing_allowed(&args, contig, variant_pos, None)?
} else {
read_vcf_to_matrix(&args, contig, variant_pos, None)?
if args.selection == Selection::OnlyAlts || args.selection == Selection::OnlyRefs {
vcf = select_carrier_haplotypes(vcf, variant_pos, &args.coords, &args.selection)
}
if args.selection == Selection::OnlyLongest {
let shared_lengths = shared_lengths_by_majority(&vcf, vcf.variant_idx());
vcf = select_only_longest_haplotypes(&shared_lengths, vcf);
find_and_write_shared_haplotype(&vcf, args.selection == Selection::Unphased, writer, contig)?;
Ok(())
#[rustfmt::skip]
pub fn find_and_write_shared_haplotype(
vcf: &PhasedMatrix,
is_rwc: bool,
mut writer: csv::Writer<Box<dyn std::io::Write>>,
contig: &str,
) -> Result<()> {
writer.write_record(vec!["contig", "pos", "ref", "alt", "gt"])?;
if is_rwc {
let (left_gt, right_gt) = find_shared_rwc_run(vcf);
let left_count = left_gt.len();
for (i, gt) in left_gt.into_iter().rev().enumerate() {
let coords = &vcf.coords().get(vcf.variant_idx - left_count + i).unwrap();
write_haplo_record(&mut writer, contig, coords, gt)?;
for (i, gt) in right_gt.into_iter().enumerate() {
let coords = &vcf.coords().get(vcf.variant_idx + i).unwrap();
let range = find_shared_haplotype_phased(vcf);
for (gt, index) in vcf.matrix.slice(s![0, range.clone()]).into_iter().zip(range) {
let coords = &vcf.coords().get(index).unwrap();
write_haplo_record(&mut writer, contig, coords, *gt)?;
fn write_haplo_record(
writer: &mut csv::Writer<Box<dyn std::io::Write + 'static>>,
coords: &crate::structs::Coord,
gt: u8,
writer.write_record(vec![
contig,
&coords.pos.to_string(),
&coords.reference,
&coords.alt,
>.to_string(),
])?;
pub fn find_shared_rwc_run(vcf: &PhasedMatrix) -> (Vec<u8>, Vec<u8>) {
let slice = vcf.matrix.slice(s![.., vcf.variant_idx..]);
let right_gt = find_breakpoint_unphased(slice);
let mut slice = vcf.matrix.slice(s![.., ..vcf.variant_idx]);
slice.invert_axis(Axis(1));
let left_gt = find_breakpoint_unphased(slice);
(left_gt, right_gt)
pub fn find_breakpoint_unphased(matrix: ArrayBase<ViewRepr<&u8>, Dim<[usize; 2]>>) -> Vec<u8> {
let mut gt_vec = vec![];
for i in 0..matrix.ncols() {
let column = matrix.slice(ndarray::s![.., i]);
let ones = column.into_iter().filter(|&n| *n == 1).count();
let alt_af = ones as f32 / matrix.nrows() as f32;
// Check no contrahomozygotes are present
let mut alt_homozygote = false;
let mut ref_homozygote = false;
for i in 0..matrix.nrows() / 2 {
if *column.slice(s![i]).into_scalar() == 0
&& *column.slice(s![i + matrix.nrows() / 2]).into_scalar() == 0
{
ref_homozygote = true
} else if *column.slice(s![i]).into_scalar() == 1
&& *column.slice(s![i + matrix.nrows() / 2]).into_scalar() == 1
alt_homozygote = true
if alt_homozygote && ref_homozygote {
break;
if alt_af >= 0.5 {
gt_vec.push(1_u8);
gt_vec.push(0_u8);
gt_vec
pub fn find_shared_haplotype_phased(vcf: &PhasedMatrix) -> Range<usize> {
let variant_slice = vcf.matrix.slice(s![.., vcf.variant_idx]);
let zeroes = variant_slice.into_iter().filter(|&n| *n == 0).count();
let ones = variant_slice.into_iter().filter(|&n| *n == 1).count();
if zeroes != 0 && ones != 0 {
return vcf.variant_idx..vcf.variant_idx;
let right_count = find_breakpoint(slice);
let left_count = find_breakpoint(slice);
let start = vcf.variant_idx - left_count;
let stop = vcf.variant_idx + right_count;
start..stop
//
pub fn find_breakpoint(matrix: ArrayBase<ViewRepr<&u8>, Dim<[usize; 2]>>) -> usize {
let mut counter = 0;
let zeroes = column.into_iter().filter(|&n| *n == 0).count();
counter += 1;
counter