This is a test version of Biostars. For the public version, visit https://www.biostars.org.
rust_htslib unrecognized filter ID

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

1 answer

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();

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.

Log in to answer this question.