misc fixes

compute_features.cpp

@@ -208,6 +208,14 @@ class Igor : Program
             {
                 continue;
             }
+
+            float qual = 0;
+            int32_t n = 0;
+            est->compute_qual(pls, no_samples, ploidy, no_alleles, qual, n);
+            if (n)
+            {
+                bcf_set_qual(v, qual);
+            }
 
             int32_t g[ploidy];
             for (int32_t i=0; i<ploidy; ++i) g[i]=0;
@@ -241,7 +249,7 @@ class Igor : Program
 
             float MLE_HWE_AF[no_alleles];
             float MLE_HWE_GF[no_genotypes];
-            int32_t n = 0;
+            n = 0;
             est->compute_gl_af_hwe(pls, no_samples, ploidy,no_alleles, MLE_HWE_AF, MLE_HWE_GF,  n, 1e-20);
             if (n)
             {

estimator.cpp

@@ -798,3 +798,71 @@ void Estimator::compute_gl_ab(int32_t *pls, int32_t no_samples, int32_t ploidy,
         ab = (0.05+num)/(0.10+denum);
     }
 };
+
+/**
+ * Computes the phred scaled QUAL for a variant.
+ *
+ * @pls        - PHRED genotype likelihoods
+ * @no_samples - number of samples
+ * @ploidy     - ploidy
+ * @no_alleles - number of alleles
+ * @n          - effective sample size
+ * @qual       - PHRED scaled QUAL
+ */
+void Estimator::compute_qual(int32_t *pls, int32_t no_samples, int32_t ploidy,
+            int32_t no_alleles, float &qual, int32_t &n)
+{
+    n = 0;
+    if (ploidy==2)
+    {
+        if (no_alleles==2)
+        {
+            int32_t no_genotypes = 3;
+
+            qual = 0;
+            for (size_t k=0; k<no_samples; ++k)
+            {
+                size_t offset = k*3;
+
+                if (pls[offset]==bcf_int32_missing) continue;
+
+                ++n;
+
+                qual += lt->log10((1-lt->pl2prob(pls[offset])/(lt->pl2prob(pls[offset])+lt->pl2prob(pls[offset+1])+lt->pl2prob(pls[offset+2]))));
+
+            }
+
+            if (!n) return;
+
+            qual = lt->round(-qual*10);
+        }
+        else
+        {
+            //works only for ploidy of 2
+
+            int32_t no_genotypes = bcf_an2gn(no_alleles);
+
+            float gq = 0;
+            for (size_t k=0; k<no_samples; ++k)
+            {
+                size_t offset = k*no_genotypes;
+
+                if (pls[offset]==bcf_int32_missing) continue;
+
+                ++n;
+
+                float denom = 0;
+                for (size_t j=0; j<no_genotypes; ++j)
+                {
+                    denom += lt->pl2prob(pls[offset]);
+                }
+
+                qual += lt->log10((1-lt->pl2prob(pls[offset])/denom));
+            }
+
+            if (!n) return;
+
+            qual = lt->round(-qual*10);
+        }
+    }
+};

estimator.h

@@ -171,6 +171,20 @@ class Estimator
                     float* GF, int32_t no_alleles,
                     float& ab, int32_t& n);
 
+
+    /**
+     * Computes the phred scaled QUAL for a variant.
+     *
+     * @pls        - PHRED genotype likelihoods
+     * @no_samples - number of samples
+     * @ploidy     - ploidy
+     * @no_alleles - number of alleles
+     * @n          - effective sample size
+     * @qual       - PHRED scaled QUAL
+     */
+    void compute_qual(int32_t *pls, int32_t no_samples, int32_t ploidy,
+                int32_t no_alleles, float &qual, int32_t &n);
+
     private:
 };
 

hts_utils.h

@@ -387,4 +387,14 @@ bool bcf_is_passed(bcf_hdr_t *h, bcf1_t *v);
  */
 #define bcf_set_n_sample(v, n) ((v)->n_sample = (n));
 
+/**
+ * Get qual
+ */
+#define bcf_get_qual(v) ((v)->qual)
+
+/**
+ * Set qual
+ */
+#define bcf_set_qual(v, q) ((v)->qual = (q))
+
 #endif