Remove obsolete distinction based on "informativeness"

- Remove obsolete distinction based on "informativeness" of a read pair prior to likelihood calculations. - Fix bug in VCF position calculation if the position of a variant was 1.
kehrlab · Aug 28, 2020 · 01d1bb6 · 01d1bb6
1 parent a7a5615
commit 01d1bb6
Show file tree

Hide file tree

Showing 4 changed files with 10 additions and 20 deletions.
diff --git a/Makefile b/Makefile
@@ -14,8 +14,8 @@ SEQAN_LIB=.
 CXXFLAGS+=-I$(SEQAN_LIB) -DSEQAN_HAS_ZLIB=1 -std=c++14 -DSEQAN_DISABLE_VERSION_CHECK
 LDLIBS=-lz -lpthread
 
-DATE=on 2020-07-30
-VERSION=1.2.1
+DATE=on 2020-08-04
+VERSION=1.2.2
 CXXFLAGS+=-DDATE=\""$(DATE)"\" -DVERSION=\""$(VERSION)"\"
 
 # Enable warnings

diff --git a/README.md b/README.md
@@ -55,8 +55,8 @@ Preprint available at bioRxiv 740225; doi: https://doi.org/10.1101/740225
 
 ## Version and License
 ```
-    Last update: 2020-07-30
-    PopDel version: 1.2.1
+    Last update: 2020-08-04
+    PopDel version: 1.2.2
     SeqAn version: 2.3.1 (modified)
     Author: Sebastian Niehus (Sebastian.Niehus[at]bihealth.de)
 ```

diff --git a/popdel_call/genotype_deletion_popdel_call.h b/popdel_call/genotype_deletion_popdel_call.h
@@ -203,18 +203,13 @@ inline Triple<long double> compute_data_likelihoods(String<Triple<long double> >
             int currentDeviation = chromosomeProfiles.getSingleDeviation(rg, it);
             long double refLikelihood = I(hist, currentDeviation - refShift);
             long double delLikelihood = I(hist, currentDeviation - deletion_length);
-            if (refLikelihood < 2 * delLikelihood && delLikelihood < 2 * refLikelihood)
-            {
-                ++it;
-                continue;               // unclear support.
-            }
 //               std::cout << "====================================================" << std::endl;
 //               std::cout << "Pos:\t" << chromosomeProfiles.currentPos << std::endl;
 //               std::cout << "Dev:\t" << currentDeviation << std::endl;
 //               std::cout << "L(REF):\t" << refLikelihood << std::endl;
 //               std::cout << "L(DEL):\t" << delLikelihood << std::endl;
             long double g0 = log(refLikelihood);
-            long double g1 = log(refLikelihood +delLikelihood) - log(2.0);    
+            long double g1 = log(refLikelihood +delLikelihood) - log(2.0);
             long double g2 = log(delLikelihood);
             currentRgWiseDataLikelihoods.i1 += g0;
             currentRgWiseDataLikelihoods.i2 += g1;
@@ -291,20 +286,12 @@ inline Triple<long double> compute_data_likelihoods(Triple<long double> & gtLogs
             long double refLikelihood = I(hist, currentDeviation - refShift);
             long double delLikelihood = I(hist, currentDeviation - deletion_length);
             if (refLikelihood >= 2 * delLikelihood)
-            {
                 ++lad.i1;       // Supporting the reference
-
-            }
             else if (delLikelihood >= 2 * refLikelihood)
-            {
                 ++lad.i3;       // Supporting a deletion
-            }
             else
-            {
                 ++lad.i2;       // Unclear support.
-                ++it;
-                continue;
-            }
+
             logLikelihoods.i1 += log(refLikelihood);
             gtLogs.i1 += log10(refLikelihood);
             logLikelihoods.i2 += log(refLikelihood +delLikelihood) - log(2.0);

diff --git a/popdel_call/vcfout_popdel_call.h b/popdel_call/vcfout_popdel_call.h
@@ -177,7 +177,10 @@ inline VcfRecord buildRecord(const QuantileMap& quantileMap, const Call& call, i
 {
     VcfRecord record;
     record.rID = rID;
-    record.beginPos = call.position;
+    if (call.position > 1)
+        record.beginPos = call.position - 1;
+    else
+        record.beginPos = call.position;
     record.id = ".";
     record.ref = "N";
     record.alt = "<DEL>";