Lines
76.3 %
Functions
12.5 %
Branches
100 %
//! Contrahomozygosity
//!
//! This tool checks for genome wide runs of homozygosity
//! To run the command for a $filename using 8 threads use:
//! ```bash
//! hatk homozygosity --file $filename -t 8
//! ```
//! To see the source code of the program, click on the upper right.
use std::{collections::BTreeMap, path::PathBuf};
use color_eyre::{eyre::eyre, Result};
use crossbeam_channel::{unbounded, Sender};
use rayon::prelude::*;
use rust_htslib::bcf::Read;
use crate::core::{read_tabix, spawn_csv_collector};
#[doc(hidden)]
#[tracing::instrument]
pub fn run(
file: PathBuf,
samples_path: Option<PathBuf>,
mut output: PathBuf,
info_limit: Option<f32>,
sort: bool,
min_marker_len: usize,
prefix: Option<String>,
) -> Result<()> {
let tabix = read_tabix(&file)?;
let contigs: Vec<String> = tabix.keys().cloned().collect();
let (tx, rx) = unbounded();
if let Some(prefix) = &prefix {
output.push(format!("{prefix}_genome_wide_rwc.csv"));
} else {
output.push("genome_wide_rwc.csv");
}
let header = ["chr", "pos", "snp_count", "bp_length"]
.iter()
.map(|&s| s.to_string())
.collect();
let handle = spawn_csv_collector(rx, header, Some(output), sort);
contigs
.par_iter()
.map_with((&file, tx), |(path, tx), contig| {
// Fetch readers
let mut reader = crate::read_vcf::get_reader(path, contig, None)?;
let samples = crate::read_vcf::get_samples(reader.header())?;
let wanted = crate::io::read_sample_ids(&samples_path)?;
let sample_indexes = crate::utils::filter_samples(&samples, wanted);
let records = reader.records();
// Program logic
let (merged_gt_vec, pos_vec) =
read_vcf_to_rwc_vec(records, &sample_indexes, info_limit)?;
let map = rwc_vec_to_btreemap(merged_gt_vec);
write_out_btreemap(map, pos_vec, min_marker_len, contig, tx)?;
Ok(())
})
.collect::<Result<()>>()?;
let join = handle.join().unwrap();
join?;
pub fn read_vcf_to_rwc_vec<R: Read>(
records: rust_htslib::bcf::Records<R>,
sample_indexes: &[usize],
) -> Result<(Vec<u8>, Vec<i64>)> {
let mut haplotypes_vec = vec![];
let mut pos_vec = vec![];
for record in records {
let record = record?;
if record.alleles().len() > 2 {
return Err(eyre!("VCF has multiallelics. Use bcftools norm."));
let passed_check = if let Some(info_limit) = info_limit {
match record.info(b"INFO").float()? {
Some(buffer) => buffer[0] >= info_limit,
None => return Err(eyre!("No INFO score on a variant at: {}", record.pos())),
true
};
if passed_check {
let gts = record.genotypes()?;
let found_contrahomozygote = sample_indexes
.map(|x| {
let gt_vec: Vec<Option<u32>> =
gts.get(*x).iter().map(|gt| gt.index()).collect();
match (gt_vec[0], gt_vec[1]) {
(Some(0), Some(0)) => 0,
(Some(1), Some(1)) => 2,
_ => 1,
.scan((false, false), |acc, gt| {
if gt == 0 {
acc.0 = true
} else if gt == 2 {
acc.1 = true
Some(*acc)
.find(|acc| acc == &(true, true));
match found_contrahomozygote {
Some(_) => haplotypes_vec.push(1),
None => haplotypes_vec.push(0),
pos_vec.push(record.pos() + 1);
Ok((haplotypes_vec, pos_vec))
pub fn rwc_vec_to_btreemap(merged_gt_vec: Vec<u8>) -> BTreeMap<usize, (usize, usize)> {
let mut counter = 1;
let mut map = BTreeMap::<usize, (usize, usize)>::new();
merged_gt_vec
.enumerate()
.fold(1, |last, (index, curr)| {
// The last and current element are 0, no contrahomozygotes are present
if last == 0 && *curr == 0 {
counter += 1;
// If the last element is 0 and the current is not, a shared sequence has ended
} else if last == 0 && *curr == 1 {
// check no duplicate with longer sequence exists
if let Some((old_counter, _)) = map.insert(index - counter, (counter, index)) {
if old_counter > counter {
map.insert(index - old_counter, (old_counter, index));
counter = 1;
*curr
});
counter -= 1;
if counter > 0 {
let index = merged_gt_vec.len() - 1;
map
pub fn write_out_btreemap(
map: BTreeMap<usize, (usize, usize)>,
pos_vec: Vec<i64>,
min_run_length: usize,
contig: &str,
tx: &mut Sender<Vec<String>>,
for (start, (count, stop)) in map {
if count >= min_run_length {
let start = &pos_vec[start];
let stop = &pos_vec[stop];
tx.send(vec![
contig.to_string(),
start.to_string(),
count.to_string(),
(stop - start).to_string(),
])
.unwrap();