Skip to content

Commit

Permalink
handle case with one read and --support 1
Browse files Browse the repository at this point in the history
should solve issue reported in #4
  • Loading branch information
wdecoster committed Jul 2, 2024
1 parent a0709e5 commit 25c43e7
Show file tree
Hide file tree
Showing 3 changed files with 53 additions and 2 deletions.
13 changes: 13 additions & 0 deletions src/genotype.rs
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,19 @@ fn genotype_repeat(
insertions.len().to_string(),
));
}
// handle the special case where there is only 1 read for the repeat (but minimal support is 1 so it passes)
// no need to make a consensus or phase, just report the single read
if insertions.len() == 1 {
return Ok(crate::vcf::VCFRecord::single_read(
&insertions[0],
repeat,
&repeat_ref_seq,
all_insertions,
reads.ps,
flags,
));
}

debug!("{repeat}: Phasing {} insertions", insertions.len(),);
let phased = crate::phase_insertions::split(&insertions, repeat, args.find_outliers);
match phased.hap2 {
Expand Down
4 changes: 2 additions & 2 deletions src/phase_insertions.rs
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ pub fn split(
let mut clusters = vec![];
// create a vector with candidate haplotype-clusters
let mut haplotype_clusters = vec![];
// clusters have to represent at least 20% of the reads
let min_cluster_size = (insertions.len() as f32 / 10.0) as usize;
// clusters have to represent at least 10% of the reads, but never less than 1
let min_cluster_size = max((insertions.len() as f32 / 10.0) as usize, 1);
debug!("{repeat}: Minimum cluster size: {}", min_cluster_size);

for (index, step) in dend.steps().iter().enumerate() {
Expand Down
38 changes: 38 additions & 0 deletions src/vcf.rs
Original file line number Diff line number Diff line change
Expand Up @@ -184,6 +184,44 @@ impl VCFRecord {
allele: (".".to_string(), ".".to_string()),
}
}

pub fn single_read(
seq: &str,
repeat: &crate::repeats::RepeatInterval,
repeat_ref_seq: &str,
all_insertions: Option<Vec<String>>,
ps: Option<u32>,
flag: Vec<String>,
) -> VCFRecord {
let somatic_info_field = match all_insertions {
Some(somatic_insertions) => {
format!(";SEQS={}", somatic_insertions.join(","))
}
None => "".to_string(),
};
VCFRecord {
chrom: repeat.chrom.clone(),
start: repeat.start,
end: repeat.end,
ref_seq: repeat_ref_seq.to_string(),
alt_seq: Some(seq.to_string()),
length: ((seq.len() as u32 - (repeat.end - repeat.start)).to_string(), ".".to_string()),
full_length: ((seq.len() as u32).to_string(), ".".to_string()),
support: ("1".to_string(), ".".to_string()),
std_dev: (".".to_string(), ".".to_string()),
score: (".".to_string(), ".".to_string()),
somatic_info_field,
outliers: "".to_string(),
ps,
flags : if flag.is_empty() {
"".to_string()
} else {
format!("{};", flag.join(";"))
},
allele: (".".to_string(), ".".to_string()),
}

}
}

impl fmt::Display for VCFRecord {
Expand Down

0 comments on commit 25c43e7

Please sign in to comment.