vcfkit 0.1.1

vcfkit: a simple program for vcf/bcf file manipulation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

use rust_htslib::bcf::{Read, Reader};
use rust_htslib::errors::Error;

#[derive(Debug)]
struct ALLELES {
    count: Vec<Vec<usize>>,
}

impl ALLELES {
    fn new(sample_count: usize, allele_count: usize) -> Self {
        ALLELES {
            count: vec![vec![0; sample_count]; allele_count],
        }
    }

    fn row_sum(&self, allele_count_idx: usize) -> usize {
        self.count[allele_count_idx].iter().sum::<usize>()
    }
    // calc site miss rate
    fn miss_rate(&self, allele_count: usize) -> f64 {
        let mut total = 0;
        for r in 0..=allele_count {
            total += self.row_sum(r);
        }
        self.row_sum(allele_count) as f64 / total as f64
    }
    // calc ref allele frequence
    fn ref_af(&self, allele_count: usize) -> f64 {
        let mut total = 0;
        for r in 0..allele_count {
            total += self.row_sum(r);
        }
        self.row_sum(0) as f64 / total as f64
    }
    // calc alt alleles frequence
    fn alts_af(&self, allele_count: usize) -> Vec<f64> {
        let mut alts: Vec<f64> = vec![];
        let mut total = 0;
        for r in 0..allele_count {
            total += self.row_sum(r);
        }
        for alt in 1..allele_count {
            alts.push(self.row_sum(alt) as f64 / total as f64);
        }
        alts
    }
}

pub fn calc_sites(
    vcf: &Option<&str>,
    //out: &Option<&str>,
) -> Result<(), Error> {
    let mut bcf = if let Some(file) = vcf {
        Reader::from_path(file)?
    } else {
        Reader::from_stdin()?
    };
    let header = bcf.header();
    let sample_count = header.sample_count() as usize;

    println!("chrom\tpos\tqual\tallele_count\tmiss_rate\tref:af\talt:af");
    for rec in bcf.records() {
        let mut alleles: Vec<String> = vec![]; // vec![ref, alt1, alt2 ...]
        if let Ok(ref row) = rec {
            for allele in row.alleles() {
                let mut each = String::new();
                for nt in allele {
                    each.push(char::from(*nt))
                }
                alleles.push(each);
            }
            let pos = row.pos() + 1;
            let allele_count = row.allele_count() as usize;

            let mut geno_count: ALLELES = ALLELES::new(sample_count, allele_count + 1);
            let genotypes = row.genotypes()?;
            for idx in 0..sample_count {
                let genotype = genotypes.get(idx);
                for geno in genotype.iter() {
                    if let Some(x) = geno.index() {
                        geno_count.count[x as usize][idx] += 1; // ref alt1 alt2 ...
                    } else {
                        geno_count.count[allele_count][idx] += 1; // calc miss
                    }
                }
            }
            // fmt output
            print!(
                "{}\t{}\t{}\t{}\t{:.4}\t",
                String::from_utf8(row.header().rid2name(row.rid().unwrap()).unwrap().to_vec())
                    .unwrap(),
                pos,
                row.qual(),
                row.allele_count(),
                geno_count.miss_rate(allele_count)
            );
            for n in 0..allele_count {
                if n == 0 {
                    print!("{}:{:.4}\t", alleles[n], geno_count.ref_af(allele_count));
                } else {
                    if allele_count == 2 {
                        print!("{}:{:.4}", alleles[n],geno_count.alts_af(allele_count)[n - 1]);    
                    } else {
                        print!("{}:{:.4},", alleles[n],geno_count.alts_af(allele_count)[n - 1]);
                    }
                }
            }
            println!();
        }
    }
    Ok(())
}