Merge pull request #270 from PolinaBevad/fix_issue_267_ref_MR

Fixed issue 267: IUPAC ambiguity codes in reference

Author: pcingola

VarDict

@@ -15,6 +15,8 @@
 import java.util.regex.Pattern;
 import java.util.regex.Matcher;
 import java.util.*;
+import java.util.stream.Collectors;
+import java.util.stream.Stream;
 
 import static com.astrazeneca.vardict.data.scopedata.GlobalReadOnlyScope.instance;
 import static com.astrazeneca.vardict.Utils.*;
@@ -35,6 +37,21 @@ public class ToVarsBuilder implements Module<RealignedVariationData, AlignedVars
     private Map<Integer, Character> ref;
     private Double duprate;
 
+    // Map of IUPAC ambiguity codes that we can observe in reference.
+    // By VCF 4.3 specification they aren't allowed and must be reduced to a first alphabetically base.
+    private Map<String, String> IUPAC_AMBIGUITY_CODES = Stream.of(new String[][] {
+            {"M","A"},
+            {"R","A"},
+            {"W","A"},
+            {"S","C"},
+            {"Y","C"},
+            {"K","G"},
+            {"V","A"},
+            {"H","A"},
+            {"D","A"},
+            {"B","C"},
+    }).collect(Collectors.toMap(key -> key[0], key -> key[1]));
+
     public Map<Integer, VariationMap<String, Variation>> getInsertionVariants() {
         return insertionVariants;
     }
@@ -903,7 +920,7 @@ else if (deletionLength < instance().conf.SVMINLEN) {
                 vref.shift3 = shift3;
                 vref.startPosition = startPosition;
                 vref.endPosition = endPosition;
-                vref.refallele = refallele;
+                vref.refallele = validateRefallele(refallele);
                 vref.varallele = varallele;
                 vref.genotype = genotype;
                 vref.totalPosCoverage = totalPosCoverage;
@@ -954,8 +971,8 @@ else if (deletionLength < instance().conf.SVMINLEN) {
             vref.highQualityReadsFrequency = roundHalfEven("0.0000", vref.highQualityReadsFrequency);
             String referenceBase = ref.containsKey(position) ? ref.get(position).toString() : ""; // $r
             //both refallele and varallele are 1 base from reference string
-            vref.refallele = referenceBase;
-            vref.varallele = referenceBase;
+            vref.refallele = validateRefallele(referenceBase);
+            vref.varallele = validateRefallele(referenceBase);
             vref.genotype = referenceBase + "/" + referenceBase;
             vref.leftseq = "";
             vref.rightseq = "";
@@ -976,6 +993,23 @@ else if (deletionLength < instance().conf.SVMINLEN) {
         }
     }
 
+    /**
+     * Validate reference allele according to VCF 4.3 specification in case if IUPAC ambiguity codes are present
+     * in reference.
+     * @param refallele sequence of reference bases that covers variant
+     * @return reference allele sequence where IUPAC ambuguity bases are changed to the one that is
+     * first alphabetically.
+     */
+    String validateRefallele(String refallele) {
+        for (int i = 0; i < refallele.length(); i++) {
+            String refBase = substr(refallele, i, 1);
+            if (IUPAC_AMBIGUITY_CODES.containsKey(refBase)) {
+                refallele = refallele.replaceFirst(refBase, IUPAC_AMBIGUITY_CODES.get(refBase));
+            }
+        }
+        return refallele;
+    }
+
     /**
      * Microsatellite instability
      * Tandemly repeated short sequence motifs ranging from 1– 6(8 in our case) base pairs are called microsatellites.

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

@@ -65,6 +65,11 @@ public void accept(Scope<AlignedVarsData> mapScope) {
                         if (vref.refallele.contains("N")) {
                             continue;
                         }
+                        if (vref.refallele.equals(vref.varallele)) {
+                            if (!conf.doPileup) {
+                                continue;
+                            }
+                        }
                         vref.vartype = vref.varType();
                         if (!vref.isGoodVar(variantsOnPosition.referenceVariant, vref.vartype, mapScope.splice)) {
                             if (!conf.doPileup) {

src/main/java/com/astrazeneca/vardict/postprocessmodules/SimplePostProcessModule.java

@@ -11,10 +11,7 @@
 import org.testng.annotations.BeforeMethod;
 import org.testng.annotations.Test;
 
-import java.util.ArrayList;
-import java.util.HashMap;
-import java.util.List;
-import java.util.Map;
+import java.util.*;
 
 import static org.testng.Assert.assertEquals;
 
@@ -132,4 +129,24 @@ public void createVariant() {
 
         assertEquals(var.get(0).toString(), expectedVariant.toString());
     }
+
+    @Test
+    public void testValidateRefAllele(){
+        List<String> alleles = Arrays.asList("A", "C", "G", "T", "N", "M", "R", "W", "S", "Y", "K", "V", "H", "D", "B");
+        List<String> expected = Arrays.asList("A", "C", "G", "T", "N", "A", "A", "A", "C", "C", "G", "A", "A", "A", "C");
+        List<String> validatedAlleles = new ArrayList<>();
+        for (String allele: alleles) {
+            validatedAlleles.add(toVarsBuilder.validateRefallele(allele));
+        }
+        assertEquals(validatedAlleles, expected);
+
+        List<String> alleles_complex = Arrays.asList("ANYCGT", "MRACT", "CCGKBG");
+        List<String> expected_complex = Arrays.asList("ANCCGT", "AAACT", "CCGGCG");
+
+        validatedAlleles = new ArrayList<>();
+        for (String allele: alleles_complex) {
+            validatedAlleles.add(toVarsBuilder.validateRefallele(allele));
+        }
+        assertEquals(validatedAlleles, expected_complex);
+    }
 }