Entering edit mode
17 months ago
raphael.B
▴
520
Hello,
I am just starting Rust and I wanted to use rust_htslib to process a VCF. Idea was to set a filter flag (here AFnfe
) in all records having AF>0.1
.
Yet I have troubles adding this filter: code compiles but fails because the id AFnfe
is not defined in the header. Could you help me finding out a proper way of adding it?
Here is the code I am running:
fn main() {
let path = &"/home/blanchet-r/Bureau/AIC_Genet/tset.vcf";
// open bcf file
let mut bcf = Reader::from_path(path).expect("Error opening file.");
// add filter to header
let mut header_in = Header::from_template(bcf.header());
header_in.push_record(br#"##FILTER=<ID=AFnfe,Description="MAF superior to 0.1">"#);
// open output bcf
let mut writer = bcf::Writer::from_path("/home/blanchet-r/Bureau/AIC_Genet/tset.filt.vcf", &header_in, true ,Format::Vcf).unwrap();
for (i, record_result) in bcf.records().enumerate() {
let mut record = record_result.unwrap();
// Check if the record has at least one AF<0.1
let filt_tag = record.info("AF".as_bytes());
let passes_filt = filt_tag.float().unwrap().iter().any(|filt| filt[0] < 0.1);
// If not, add AFnfe to FILTER
if ! passes_filt{
// TRY: add filter to record header ? (should already have it)
let mut header = Header::from_template(record.header());
header.push_record(br#"##FILTER=<ID=AFnfe,Description="MAF superior to 0.1">"#);
// get filter id and push to FILTER
let af = record.header().name_to_id(b"AFnfe").unwrap();
record.push_filter(&af);
println!("{:?}", record.filters());
}
// write records
writer.write(&record);
}
}
Error:
thread 'main' panicked at 'called `Result::unwrap()` on an `Err` value: BcfUnknownID { id: "AFnfe" }', src/main.rs:185:59
PS: This code runs just fine when changing filter id to "AF", for example, since it is defined in the header from the beginning.
Thanks for the suggestions!
name_to_id
is aHeaderView
method I'm afraid. But even when pushing the filter id to the record header earlier (beforeif
for exemple), the error persists.