update tools

.github/workflows/linux.yml

@@ -27,14 +27,20 @@ jobs:
       run: |
         git clone --recursive https://github.com/samtools/htslib.git
         cd htslib && autoreconf -i && ./configure --enable-libcurl
-        sudo make -j 4 install
+        sudo make -j6 install
         sudo ldconfig
         cd -
 
-    - name: test
+    - name: unit test
       run: |
         cd test
         # export ASAN_OPTIONS=alloc_dealloc_mismatch=0
         make
         make test
         # find . -name "*.bin" -exec '{}' --success ';'
+
+    - name: tools test
+      run: |
+        cd tools
+        make -j6
+        make test

tools/Makefile

@@ -9,7 +9,7 @@ LIBS     = -lhts -lz -lm -lbz2 -llzma -lcurl -lpthread
 # OBJS     = $(patsubst %.cpp, %.o, $(wildcard *.cpp))
 BINS     = $(patsubst %.cpp, %, $(wildcard *.cpp))
 
-.PHONY: all clean
+.PHONY: all clean test
 
 all: $(BINS)
 
@@ -18,3 +18,7 @@ all: $(BINS)
 
 clean:
 	rm -f $(BINS)
+
+test:
+	./vcf_addDS -i test.vcf.gz -o out.vcf.gz
+	./vcf_addINFO -i test.vcf.gz -o out.vcf.gz

tools/vcf_addDS.cpp

@@ -1,10 +1,10 @@
-// -*- compile-command: "g++ vcf_addDS.cpp -std=c++11 -g -O3 -Wall -lhts" -*-
+// -*- compile-command: "g++ vcf_addDS.cpp -std=c++11 -g -O3 -Wall -I.. -lhts" -*-
 #include "vcfpp.h"
 
 using namespace std;
 using namespace vcfpp;
 
-int main(int argc, char* argv[])
+int main(int argc, char * argv[])
 {
     // ========= helper message and parameters parsing ============================
     std::vector<std::string> args(argv + 1, argv + argc);
@@ -26,40 +26,36 @@ int main(int argc, char* argv[])
         return 1;
     }
     std::string invcf, outvcf = "-", samples = "-", region = "";
-    for (int i = 0; i < args.size(); i++)
+    for(size_t i = 0; i < args.size(); i++)
     {
-        if (args[i] == "-i")
-            invcf = args[++i];
-        if (args[i] == "-o")
-            outvcf = args[++i];
-        if (args[i] == "-s")
-            samples = args[++i];
-        if (args[i] == "-r")
-            region = args[++i];
+        if(args[i] == "-i") invcf = args[++i];
+        if(args[i] == "-o") outvcf = args[++i];
+        if(args[i] == "-s") samples = args[++i];
+        if(args[i] == "-r") region = args[++i];
     }
     // ========= core calculation part ===========================================
     BcfReader vcf(invcf, region, samples);
-    BcfRecord var(vcf.header);
     BcfWriter bw(outvcf, vcf.header);
+    bw.copyHeader(invcf);
     bw.header.addFORMAT("DS", "1", "Float", "Diploid Genotype Dosage"); // add DS tag into the header
     int nsamples = vcf.nsamples;
     vector<float> gps, ds(nsamples);
     vector<char> gts;
-    while (vcf.getNextVariant(var))
+    BcfRecord var(bw.header);
+    while(vcf.getNextVariant(var))
     {
         try
         {
             var.getFORMAT("GP", gps); // try get GP values
-            for (int i = 0; i < nsamples; i++)
+            for(int i = 0; i < nsamples; i++)
             {
                 ds[i] = gps[i * 3 + 1] + gps[i * 3 + 2] * 2;
             }
         }
-        catch (const std::runtime_error& e)
+        catch(const std::runtime_error & e)
         {
             var.getGenotypes(gts);
-            for (int i = 0; i < nsamples; i++)
-                ds[i] = gts[i * 2] + gts[i * 2 + 1];
+            for(int i = 0; i < nsamples; i++) ds[i] = gts[i * 2] + gts[i * 2 + 1];
         }
         var.setFORMAT("DS", ds);
         bw.writeRecord(var);

tools/vcf_addINFO.cpp

@@ -0,0 +1,73 @@
+// -*- compile-command: "g++ vcf_addINFO.cpp -std=c++11 -g -O3 -Wall -I.. -lhts" -*-
+#include "vcfpp.h"
+#include <cmath>
+
+using namespace std;
+using namespace vcfpp;
+
+int main(int argc, char * argv[])
+{
+    // ========= helper message and parameters parsing ============================
+    std::vector<std::string> args(argv + 1, argv + argc);
+    if(argc <= 1 || args[0] == "-h" || args[0] == "-help" || args[0] == "--help")
+    {
+        std::cout << "Author: Zilong-Li ([email protected])\n"
+                  << "Description:\n"
+                  << "     calculate INFO score per site given GP from input vcf file\n\n"
+                  << "Usage example:\n"
+                  << "     " + (std::string)argv[0] + " -i in.bcf \n"
+                  << "     " + (std::string)argv[0] + " -i in.bcf -o out.bcf -s ^S1,S2 -r chr1:1-1000 \n"
+                  << "     bcftools view in.bcf | " + (std::string)argv[0] + " -i - -o out.bcf \n"
+                  << "\nOptions:\n"
+                  << "     -i    input vcf/bcf file\n"
+                  << "     -o    ouput vcf/bcf file [stdout]\n"
+                  << "     -s    list of samples to be included or excluded\n"
+                  << "     -r    specific region to be included\n"
+                  << std::endl;
+        return 1;
+    }
+    std::string invcf, outvcf = "-", samples = "-", region = "";
+    for(size_t i = 0; i < args.size(); i++)
+    {
+        if(args[i] == "-i") invcf = args[++i];
+        if(args[i] == "-o") outvcf = args[++i];
+        if(args[i] == "-s") samples = args[++i];
+        if(args[i] == "-r") region = args[++i];
+    }
+    // ========= core calculation part ===========================================
+    BcfReader vcf(invcf, region, samples);
+    BcfWriter bw(outvcf);
+    bw.copyHeader(invcf);
+    bw.header.addINFO("INFO", "1", "Float", "INFO score given genotype probability");
+    bw.header.addINFO("EAF", "1", "Float", "Estimated Allele Frequency");
+    int N = vcf.nsamples, i{0};
+    vector<float> gps;
+    float info, eaf, eij, fij, a0, a1, thetaHat;
+    BcfRecord var(bw.header);
+    while(vcf.getNextVariant(var))
+    {
+        var.getFORMAT("GP", gps); // try get GP values
+        for(eij = 0.0, fij = 0.0, i = 0; i < N; i++)
+        {
+            a0 = gps[i * 3 + 1] + gps[i * 3 + 2] * 2;
+            a1 = gps[i * 3 + 1] + gps[i * 3 + 2] * 4;
+            eij += a0;
+            fij += a1 - a0 * a0;
+        }
+        eaf = eij / 2 / N;
+        info = 1.0 - fij / (eij * (1 - eaf));
+        thetaHat = std::round(1e2 * eaf) / 1e2;
+        if(thetaHat == 0 || thetaHat == 1)
+            info = 1.0;
+        else if(info < 0)
+            info = 0.0;
+        else
+            info = std::round(1e3 * info) / 1e3;
+        eaf = std::round(1e6 * eaf) / 1e6;
+        var.setINFO("INFO", info);
+        var.setINFO("EAF", eaf);
+        bw.writeRecord(var);
+    }
+
+    return 0;
+}

tools/vcf_multiallelics.cpp

@@ -27,7 +27,7 @@ int main(int argc, char* argv[])
         return 1;
     }
     std::string invcf, outvcf="-", samples = "-", region = "";
-    for (int i = 0; i < args.size(); i++)
+    for(size_t i = 0; i < args.size(); i++)
     {
         if (args[i] == "-i")
             invcf = args[++i];