Skip to content

Commit

Permalink
feat: warn on missing ANN records in mutational burden estimation (#376)
Browse files Browse the repository at this point in the history
### Description

<!--Add a description of your PR here-->

### QC
<!-- Make sure that you can tick the boxes below. -->

* [ ] The PR contains a test case for the changes or the changes are
already covered by an existing test case.
* [x] The documentation at
https://github.com/varlociraptor/varlociraptor.github.io is updated in a
separate PR to reflect the changes or this is not necessary (e.g. if the
change does neither modify the calling grammar nor the behavior or
functionalities of Varlociraptor).
  • Loading branch information
johanneskoester committed Jun 12, 2023
1 parent 508333e commit 32006fb
Showing 1 changed file with 22 additions and 18 deletions.
40 changes: 22 additions & 18 deletions src/estimation/mutational_burden.rs
Original file line number Diff line number Diff line change
Expand Up @@ -15,27 +15,31 @@ use crate::{Event, SimpleEvent};

/// Consider only variants in coding regions.
/// We rely on the ANN field for this.
fn is_valid_variant(rec: &mut bcf::Record) -> Result<bool> {
for ann in rec
.info(b"ANN")
.string()?
.expect("ANN field not found. Annotate VCF with e.g. snpEff.")
.iter()
{
let mut coding = false;
for (i, entry) in ann.split(|c| *c == b'|').enumerate() {
if i == 7 {
coding = entry == b"protein_coding";
fn is_valid_variant(rec: &mut bcf::Record, header: &bcf::header::HeaderView) -> Result<bool> {
if let Some(ann) = rec.info(b"ANN").string()? {
for ann in ann.iter() {
let mut coding = false;
for (i, entry) in ann.split(|c| *c == b'|').enumerate() {
if i == 7 {
coding = entry == b"protein_coding";
}
if i == 13 {
coding &= entry != b"";
}
}
if i == 13 {
coding &= entry != b"";
if coding {
return Ok(true);
}
}
if coding {
return Ok(true);
}
Ok(false)
} else {
warn!(
"No ANN field found in record at {}:{}.",
std::str::from_utf8(header.rid2name(rec.rid().unwrap())?)?,
rec.pos() + 1
);
Ok(false)
}
Ok(false)
}

#[derive(Debug, Clone, Serialize)]
Expand Down Expand Up @@ -127,7 +131,7 @@ pub(crate) fn collect_estimates(
let vafs = rec.format(b"AF").float()?[*id].to_owned();
vafmap.insert(name, vafs);
}
if !is_valid_variant(&mut rec)? {
if !is_valid_variant(&mut rec, &header)? {
info!(
"Skipping variant {}:{} because it is not coding.",
contig, vcfpos
Expand Down

0 comments on commit 32006fb

Please sign in to comment.