rust_htslib unrecognized filter ID
1
2
Entering edit mode
10 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.

htslib rust SNP filtering • 621 views
ADD COMMENT
0
Entering edit mode
10 months ago

You may have already tried this, but I would try:

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 = header().name_to_id(b"AFnfe").unwrap();

or

let af = header_in.name_to_id(b"AFnfe").unwrap();
ADD COMMENT
0
Entering edit mode

Thanks for the suggestions! name_to_idis a HeaderView method I'm afraid. But even when pushing the filter id to the record header earlier (before if for exemple), the error persists.

ADD REPLY

Login before adding your answer.

Traffic: 1690 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6