Lines
91.51 %
Functions
50 %
Branches
100 %
use std::path::PathBuf;
use color_eyre::{eyre::eyre, Result};
use plotters::prelude::*;
use crate::args::{GraphArgs, Selection, StandardArgs};
use crate::read_vcf::read_vcf_to_matrix;
use crate::structs::{HapVariant, PhasedMatrix, SharedLength};
use crate::utils::{
push_to_output, push_to_output_mbah, select_carrier_haplotypes, select_only_longest_alleles,
shared_lengths_by_majority,
};
use crate::core::{open_csv_writer, parse_snp_coord};
#[doc(hidden)]
#[tracing::instrument]
pub fn run(
args: StandardArgs,
graph_args: GraphArgs,
x_start: Option<i64>,
x_stop: Option<i64>,
) -> Result<()> {
if args.selection == Selection::Unphased {
return Err(eyre!("Running with unphased data is not supported."))
}
let mut decay_csv = args.output.clone();
push_to_output(&args, &mut decay_csv, "shares_of_mbah", "csv");
let mut decay_writer = open_csv_writer(decay_csv)?;
let mut ancestral_haplo_csv = args.output.clone();
push_to_output_mbah(&args, &mut ancestral_haplo_csv, "mbah", "csv");
let mut ancestral_writer = open_csv_writer(ancestral_haplo_csv)?;
let (contig, variant_pos) = parse_snp_coord(&args.coords)?;
let mut vcf = 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)
let shared_lengths = shared_lengths_by_majority(&vcf, vcf.variant_idx());
let only_longest = args.selection == Selection::OnlyLongest;
print_haplotype_lengths_csv(&shared_lengths, &vcf, only_longest, &mut decay_writer)?;
let haplotype = find_ancestral_haplotype(contig, &shared_lengths, &vcf)?;
print_ancestral_haplotype_csv(&haplotype, contig, &mut ancestral_writer)?;
draw_plot(
&args,
graph_args,
args.output.clone(),
shared_lengths,
only_longest,
x_start,
x_stop,
)?;
Ok(())
fn draw_plot(
st_args: &StandardArgs,
args: GraphArgs,
mut output: PathBuf,
sl: Vec<SharedLength>,
only_longest: bool,
push_to_output(st_args, &mut output, "shares_of_mbah", "png");
let root = BitMapBackend::new(&output, (args.width, args.height)).into_drawing_area();
let mut sl = match only_longest {
true => select_only_longest_alleles(&sl),
false => sl.iter().collect(),
sl.sort_by(|a, b| b.end_pos.partial_cmp(&a.end_pos).unwrap());
sl.sort_by(|a, b| a.start_pos.partial_cmp(&b.start_pos).unwrap());
let x_start = match x_start {
Some(x_start) => x_start,
None => {
sl.iter()
.min_by(|a, b| a.start_pos.cmp(&b.start_pos))
.unwrap()
.start_pos
let x_stop = match x_stop {
Some(x_stop) => x_stop,
.max_by(|a, b| a.end_pos.cmp(&b.end_pos))
.end_pos
root.fill(&WHITE)?;
let mut chart = ChartBuilder::on(&root)
.margin(10)
.set_label_area_size(LabelAreaPosition::Left, 60)
.set_label_area_size(LabelAreaPosition::Right, 60)
.set_label_area_size(LabelAreaPosition::Bottom, 40)
.build_cartesian_2d(x_start..x_stop, 0..sl.len())?;
chart
.configure_mesh()
.disable_x_mesh()
.disable_y_mesh()
.x_labels(30)
.max_light_lines(4)
.y_desc("SAMPLE")
.x_desc("POS")
.draw()?;
for (idx, length) in sl.iter().enumerate() {
chart.draw_series(LineSeries::new(
[(length.start_pos, idx), (length.end_pos, idx)],
&BLACK,
))?;
root.present()?;
fn print_ancestral_haplotype_csv(
haplotype: &Vec<HapVariant>,
contig: &str,
writer: &mut csv::Writer<Box<dyn std::io::Write>>,
writer.write_record(vec!["contig", "pos", "ref", "alt", "gt"])?;
for variant in haplotype {
writer.write_record(vec![
contig,
&variant.pos.to_string(),
&variant.reference,
&variant.alt,
&variant.gt.to_string(),
])?;
fn find_ancestral_haplotype(
shared_lengths: &[SharedLength],
vcf: &PhasedMatrix,
) -> Result<Vec<HapVariant>> {
let max = shared_lengths.iter().max_by_key(|l| l.end_pos).unwrap();
let min = shared_lengths.iter().min_by_key(|l| l.start_pos).unwrap();
let left = min.start()..vcf.variant_idx();
let right = vcf.variant_idx()..max.end() + 1;
let mut left = vcf.find_haplotype_for_sample(contig, left, min.index);
let right = vcf.find_haplotype_for_sample(contig, right, max.index);
left.extend(right);
Ok(left)
fn print_haplotype_lengths_csv(
shared_lengths: &Vec<SharedLength>,
writer.write_record(vec!["id", "start", "stop"])?;
match only_longest {
true => {
for length in select_only_longest_alleles(shared_lengths) {
vcf.get_sample_name(length.index),
length.start_pos.to_string(),
length.end_pos.to_string(),
false => {
for length in shared_lengths {
if length.index >= vcf.nsamples() / 2 {
format!("{}_2", vcf.get_sample_name(length.index)),
} else {
format!("{}_1", vcf.get_sample_name(length.index)),