protected (some) complex variants from left-trimming in normalize

test/normalize/01_IN.vcf

@@ -94,6 +94,7 @@
 ##INFO=<ID=refseq.positionType,Number=1,Type=String,Description="RefSeq genome type position">
 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+20	61240	.	CAT	CG	1.01978e-13	.	AC=0;AF=0;AN=3	GT:GQ:DP:DPR:RO:QR:AO:QA	0/0/0:136.293:639:639,27:612:24483:27:1091
 20	421808	.	A	ACCA	.	PASS	VC=INDEL;AC=24;AF=0.08;AN=316;refseq.name=NM_144628;refseq.positionType=intron
 20	1292033	.	C	CTTGT	.	PASS	VC=INDEL;AC=28;AF=0.1;AN=276;refseq.name=NM_080489;refseq.positionType=intron
 20	1340527	.	T	TGTC	.	PASS	VC=INDEL;AC=56;AF=0.18;AN=316

test/normalize/01_OUT.stderr

@@ -29,7 +29,7 @@ stats: biallelic
        total no. multiallelic normalized        : 0
 
        total no. variants normalized            : 80
-       total no. variants observed              : 194
+       total no. variants observed              : 195
        total no. reference observed             : 0
 
 Time elapsed <stripped>

test/normalize/01_OUT.vcf

@@ -95,6 +95,7 @@
 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
 ##INFO=<ID=OLD_VARIANT,Number=.,Type=String,Description="Original chr:pos:ref:alt encoding">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+20	61240	.	CAT	CG	1.01978e-13	.	AC=0;AF=0;AN=3
 20	421805	.	T	TCCA	.	PASS	VC=INDEL;AC=24;AF=0.08;AN=316;refseq.name=NM_144628;refseq.positionType=intron;OLD_VARIANT=20:421808:A/ACCA
 20	1292033	.	C	CTTGT	.	PASS	VC=INDEL;AC=28;AF=0.1;AN=276;refseq.name=NM_080489;refseq.positionType=intron
 20	1340527	.	T	TGTC	.	PASS	VC=INDEL;AC=56;AF=0.18;AN=316

variant_manip.cpp

@@ -90,15 +90,15 @@ int32_t VariantManip::is_not_ref_consistent(bcf_hdr_t *h, bcf1_t *v)
             }
         }
     }
-    
+
     if (is_not_consistent)
     {
-       fprintf(stderr, "[%s:%d %s] reference bases not consistent: %s:%d-%d  %s(REF) vs %s(FASTA)\n", __FILE__, __LINE__, __FUNCTION__, 
-                                                                                chrom, pos0, pos0+rlen-1, vcf_ref, ref);       
-    }    
+       fprintf(stderr, "[%s:%d %s] reference bases not consistent: %s:%d-%d  %s(REF) vs %s(FASTA)\n", __FILE__, __LINE__, __FUNCTION__,
+                                                                                chrom, pos0, pos0+rlen-1, vcf_ref, ref);
+    }
 
     if (ref_len) free(ref);
-   
+
     return is_not_consistent;
 }
 
@@ -585,33 +585,40 @@ void VariantManip::right_trim_or_left_extend(std::vector<std::string>& alleles,
 /**
  * Left trims a variant.
  */
-void VariantManip::left_trim(std::vector<std::string>& alleles, int32_t& pos1, int32_t& left_trimmed)
-{
-    bool to_left_trim =  true;
+void VariantManip::left_trim(std::vector<std::string>& alleles, int32_t& pos1, int32_t& left_trimmed) {
+  bool to_left_trim = true;
+
+  while (to_left_trim) {
+    // Checks if left trimmable.
+    for (size_t i = 0; i < alleles.size(); ++i) {
+      if (alleles[i].size() == 1 || alleles[i].at(0) != alleles[0].at(0)) {
+        to_left_trim = false;
+        break;
+      }
+    }
 
-    while (to_left_trim)
-    {
-        //checks if left trimmable.
-        for (size_t i=0; i<alleles.size(); ++i)
-        {
-            if (alleles[i].size()==1 || alleles[i].at(0)!=alleles[0].at(0))
-            {
-                to_left_trim = false;
-                break;
-            }
+    // The above rule does not prevent the trimming of complex alleles, e.g., CAT -> CG.
+    // Correct that in the case of two alleles. If allele lengths are different, require that
+    // they share at least one common base.
+    if (alleles.size() == 2) {
+      for (size_t i = 0; i < alleles[0].length() and to_left_trim; ++i) {
+        for (size_t j = 0; j < alleles[1].length() and to_left_trim; ++j) {
+          if (alleles[0].at(i) == alleles[0].at(j)) {
+            to_left_trim = false;
+          }
         }
+      }
+    }
 
-        if (to_left_trim)
-        {
-            for (size_t i=0; i<alleles.size(); ++i)
-            {
-                alleles[i].erase(0, 1);
-            }
+    if (to_left_trim) {
+      for (size_t i = 0; i < alleles.size(); ++i) {
+        alleles[i].erase(0, 1);
+      }
 
-            ++pos1;
-            ++left_trimmed;
-        }
+      ++pos1;
+      ++left_trimmed;
     }
+  }
 };
 
 /**
@@ -769,4 +776,4 @@ void VariantManip::generate_probes(const char* chrom,
         //std::cerr << probes[0] << "\n" << probes[1] << "\n";
         generate_probes(chrom, pos1, flankLength, currentDiff, length, gald, diff, alleles, probes);
     }
-}
+}