Merge pull request #291 from PolinaBevad/issue_285_and_288

Issues 285 and 288 fixes

Author: pcingola

VarDict

@@ -30,6 +30,13 @@
  * and variants in region.
  */
 public class ToVarsBuilder implements Module<RealignedVariationData, AlignedVarsData> {
+    // Left 20 and right 20 bases in reference to output in VCF file
+    public static final int REF_20_BASES = 20;
+    // Number of bases where to check for MSIs
+    public static final int REF_30_BASES = 30;
+    public static final int REF_50_BASES = 50;
+    public static final int REF_70_BASES = 70;
+
     private Region region;
     private Map<Integer, Integer> refCoverage;
     private Map<Integer, VariationMap<String, Variation>> insertionVariants;
@@ -196,10 +203,10 @@ private boolean isTheSameVariationOnRef(int position, VariationMap<String, Varia
      */
     private MSI proceedVrefIsDeletion(int position, int dellen) {
         //left 70 bases in reference sequence
-        String leftseq = joinRef(ref, (position - 70 > 1 ? position - 70 : 1), position - 1); // left 10 nt
+        String leftseq = joinRef(ref, (Math.max(position - REF_70_BASES, 1)), position - 1); // left 10 nt
         int chr0 = getOrElse(instance().chrLengths, region.chr, 0);
         //right 70 + dellen bases in reference sequence
-        String tseq = joinRef(ref, position, position + dellen + 70 > chr0 ? chr0 : position + dellen + 70);
+        String tseq = joinRef(ref, position, Math.min(position + dellen + REF_70_BASES, chr0));
 
         //Try to adjust for microsatellite instability
         MSI msiData = findMSI(substr(tseq, 0, dellen), substr(tseq, dellen), leftseq);
@@ -223,7 +230,7 @@ private MSI proceedVrefIsDeletion(int position, int dellen) {
 
     /**
      * For insertion variant calculates shift3 (number of bases to be shifted to 3' for deletions due to alternative alignment),
-     * msi (which indicates Microsatellite instability) and msint (MicroSattelite unit length in base pairs) fields.
+     * msi (which indicates Microsatellite instability) and msint (MicroSatellite unit length in base pairs) fields.
      * @param position position to seek in reference to get left sequence of variant
      * @param vn variant description string
      * @return tuple of msi, shift3 and msint.
@@ -232,10 +239,10 @@ private MSI proceedVrefIsInsertion(int position, String vn) {
         //variant description string without first symbol '+'
         String tseq1 = vn.substring(1);
         //left 50 bases in reference sequence
-        String leftseq = joinRef(ref, position - 50 > 1 ? position - 50 : 1, position); // left 10 nt
+        String leftseq = joinRef(ref, Math.max(position - REF_50_BASES, 1), position); // left 10 nt
         int x = getOrElse(instance().chrLengths, region.chr, 0);
         //right 70 bases in reference sequence
-        String tseq2 = joinRef(ref, position + 1, (position + 70 > x ? x : position + 70));
+        String tseq2 = joinRef(ref, position + 1, (Math.min(position + REF_70_BASES, x)));
 
         MSI msiData = findMSI(tseq1, tseq2, leftseq);
         double msi = msiData.msi;
@@ -741,9 +748,9 @@ else if (deletionLength < instance().conf.SVMINLEN) {
                     }
                 } else { //Not insertion/deletion variant. SNP or MNP
                     //Find MSI adjustment
-                    String tseq1 = joinRef(ref, position - 30 > 1 ? position - 30 : 1, position + 1);
+                    String tseq1 = joinRef(ref, Math.max(position - REF_30_BASES, 1), position + 1);
                     int chr0 = getOrElse(instance().chrLengths, region.chr, 0);
-                    String tseq2 = joinRef(ref, position + 2, position + 70 > chr0 ? chr0 : position + 70);
+                    String tseq2 = joinRef(ref, position + 2, Math.min(position + REF_70_BASES, chr0));
 
                     MSI msiData = findMSI(tseq1, tseq2, null);
                     msi = msiData.msi;
@@ -856,9 +863,7 @@ else if (deletionLength < instance().conf.SVMINLEN) {
                 if (cutSite != 0 && (refallele.length() != varallele.length() && substr(refallele, 0, 1).equals(substr(varallele, 0, 1)))) {
                     if (!(startPosition == cutSite || endPosition == cutSite)) {
                         int n = 0;
-                        int dis = abs(cutSite - startPosition) < abs(cutSite - endPosition)
-                                ? abs(cutSite - startPosition)
-                                : abs(cutSite - endPosition);
+                        int dis = Math.min(abs(cutSite - startPosition), abs(cutSite - endPosition));
                         if (startPosition < cutSite) {
                             while(startPosition + n < cutSite && n < shift3 && endPosition + n != cutSite) {
                                 n++;
@@ -899,10 +904,10 @@ else if (deletionLength < instance().conf.SVMINLEN) {
                 }
 
                 //preceding reference sequence
-                vref.leftseq = joinRef(ref, startPosition - 20 < 1 ? 1 : startPosition - 20, startPosition - 1); // left 20 nt
+                vref.leftseq = joinRef(ref, Math.max(startPosition - REF_20_BASES, 1), startPosition - 1); // left 20 nt
                 int chr0 = getOrElse(instance().chrLengths, region.chr, 0);
                 //following reference sequence
-                vref.rightseq = joinRef(ref, endPosition + 1, endPosition + 20 > chr0 ? chr0 : endPosition + 20); // right 20 nt
+                vref.rightseq = joinRef(ref, endPosition + 1, Math.min(endPosition + REF_20_BASES, chr0)); // right 20 nt
                 //genotype description string
                 String genotype = genotype1current + "/" + genotype2;
                 //remove '&' and '#' symbols from genotype string
@@ -960,6 +965,29 @@ else if (deletionLength < instance().conf.SVMINLEN) {
             updateRefVariant(position, totalPosCoverage, vref, debugLines,
                     referenceForwardCoverage, referenceReverseCoverage);
         }
+
+        if (instance().conf.doPileup && variationsAtPos.referenceVariant != null) {
+            Variant refVar = variationsAtPos.referenceVariant;
+            fillReferenceVarInPileup(position, refVar);
+        }
+    }
+
+    /**
+     * Fill reference variant fields for pileup: positions, varallele and refallele
+     * @param position current position in reference
+     * @param refVar reference Variant
+     */
+    private void fillReferenceVarInPileup(int position, Variant refVar) {
+        refVar.startPosition = position;
+        refVar.endPosition = position;
+
+        char emptyChar = 0;
+        if (refVar.refallele.equals("")) {
+            refVar.refallele = ref.getOrDefault(position, emptyChar).toString();
+        }
+        if (refVar.varallele.equals("")) {
+            refVar.varallele = ref.getOrDefault(position, emptyChar).toString();
+        }
     }
 
     private void updateRefVariant(int position, int totalPosCoverage, Variant vref, List<String> debugLines,

src/main/java/com/astrazeneca/vardict/modules/ToVarsBuilder.java

@@ -2452,10 +2452,17 @@ void adjRefFactor(Variation ref, double factor_f) {
         ref.varsCount -= (int) (factor_f * ref.varsCount);
         ref.highQualityReadsCount -= (int) (factor_f * ref.highQualityReadsCount);
         ref.lowQualityReadsCount -= (int) (factor_f * ref.lowQualityReadsCount);
+
+        // Adjust mean mapping quality, mean quality and mean position only on number of changed counts
         double factorCnt = oldVarsCount != 0 ? Math.abs((ref.varsCount - oldVarsCount)) / (double) oldVarsCount : 1;
-        ref.meanPosition -= ref.meanPosition * factor_f * factorCnt;
-        ref.meanQuality -= ref.meanQuality * factor_f * factorCnt;
-        ref.meanMappingQuality -= ref.meanMappingQuality * factor_f * factorCnt;
+        // Factors must be the same sign
+        if (factor_f < 0 && factorCnt > 0 || factor_f > 0 && factorCnt < 0){
+            factorCnt = -factorCnt;
+        }
+        ref.meanPosition -= ref.meanPosition * factorCnt;
+        ref.meanQuality -= ref.meanQuality * factorCnt;
+        ref.meanMappingQuality -= ref.meanMappingQuality * factorCnt;
+
         ref.numberOfMismatches -= factor_f * ref.numberOfMismatches;
         ref.varsCountOnForward -= (int) (factor_f * ref.varsCountOnForward);
         ref.varsCountOnReverse -= (int) (factor_f * ref.varsCountOnReverse);