Merge pull request #271 from PolinaBevad/fix_issue_265_pileup_missed_position

pcingola · web-flow · commit db17547f6e50 · 2019-11-07T10:46:18.000-05:00
Fixed issue 265: missed position in pileup when insertion/complex occurs
diff --git a/src/main/java/com/astrazeneca/vardict/modules/ToVarsBuilder.java b/src/main/java/com/astrazeneca/vardict/modules/ToVarsBuilder.java
@@ -625,7 +625,7 @@ && getNonInsertionVariants().get(position + 1).containsKey(ref.get(position + 1)
             referenceForwardCoverage = tpref.varsCountOnForward;
             referenceReverseCoverage = tpref.varsCountOnReverse;
         }
-
+        List<Integer> positionsForChangedRefVariant = new ArrayList<>();
         // only reference reads are observed.
         if (variationsAtPos.variants.size() > 0) { //Condition: non-reference variants are found
             //Loop over non-reference variants
@@ -936,61 +936,75 @@ else if (deletionLength < instance().conf.SVMINLEN) {
                 }
 
                 adjustVariantCounts(position, vref);
-
-                //Construct debug lines
-                if (instance().conf.debug) {
-                    StringBuilder sb = new StringBuilder();
-                    for (String str : debugLines) {
-                        if (sb.length() > 0) {
-                            sb.append(" & ");
-                        }
-                        sb.append(str);
-                    }
-                    vref.DEBUG = sb.toString();
+                if (startPosition != position && instance().conf.doPileup) {
+                    positionsForChangedRefVariant.add(position);
                 }
+                constructDebugLines(debugLines, vref);
             }
             //TODO: It is a "lazy" solution because current logic in realignment methods can't be changed simply for --nosv option
             if (instance().conf.disableSV) {
                 variationsAtPos.variants.removeIf(vref -> ANY_SV.matcher(vref.varallele).find());
             }
-        } else if (variationsAtPos.referenceVariant != null) {
-            Variant vref = variationsAtPos.referenceVariant; //no variant reads are detected.
-            vref.totalPosCoverage = totalPosCoverage;
-            vref.positionCoverage = 0;
-            vref.frequency = 0;
-            vref.refForwardCoverage = referenceForwardCoverage;
-            vref.refReverseCoverage = referenceReverseCoverage;
-            vref.varsCountOnForward = 0;
-            vref.varsCountOnReverse = 0;
-            vref.msi = 0;
-            vref.msint = 0;
-            vref.strandBiasFlag += ";0";
-            vref.shift3 = 0;
-            vref.startPosition = position;
-            vref.endPosition = position;
-            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 = validateRefallele(referenceBase);
-            vref.varallele = validateRefallele(referenceBase);
-            vref.genotype = referenceBase + "/" + referenceBase;
-            vref.leftseq = "";
-            vref.rightseq = "";
-            vref.duprate = duprate;
-            //Construct debug lines
-            if (instance().conf.debug) {
-                StringBuilder sb = new StringBuilder();
-                for (String str : debugLines) {
-                    if (sb.length() > 0) {
-                        sb.append(" & ");
-                    }
-                    sb.append(str);
-                }
-                vref.DEBUG = sb.toString();
-            }
+        } else if (variationsAtPos.referenceVariant != null ) {  //no variant reads are detected
+            Variant vref = variationsAtPos.referenceVariant;
+            updateRefVariant(position, totalPosCoverage, vref, debugLines,
+                    referenceForwardCoverage, referenceReverseCoverage);
         } else {
             variationsAtPos.referenceVariant = new Variant();
         }
+
+        // Update reference variants if there were indels and start position were changed, so after update
+        // ref variants can be output in pileup
+        if (positionsForChangedRefVariant.contains(position) && variationsAtPos.referenceVariant != null) {
+            Variant vref = variationsAtPos.referenceVariant;
+            updateRefVariant(position, totalPosCoverage, vref, debugLines,
+                    referenceForwardCoverage, referenceReverseCoverage);
+        }
+    }
+
+    private void updateRefVariant(int position, int totalPosCoverage, Variant vref, List<String> debugLines,
+                                  int referenceForwardCoverage, int referenceReverseCoverage) {
+        vref.totalPosCoverage = totalPosCoverage;
+        vref.positionCoverage = 0;
+        vref.frequency = 0;
+        vref.refForwardCoverage = referenceForwardCoverage;
+        vref.refReverseCoverage = referenceReverseCoverage;
+        vref.varsCountOnForward = 0;
+        vref.varsCountOnReverse = 0;
+        vref.msi = 0;
+        vref.msint = 0;
+        vref.strandBiasFlag += ";0";
+        vref.shift3 = 0;
+        vref.startPosition = position;
+        vref.endPosition = position;
+        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 = validateRefallele(referenceBase);
+        vref.varallele = validateRefallele(referenceBase);
+        vref.genotype = referenceBase + "/" + referenceBase;
+        vref.leftseq = "";
+        vref.rightseq = "";
+        vref.duprate = duprate;
+
+        constructDebugLines(debugLines, vref);
+    }
+
+    /**
+     * Construct DEBUG lines for the variant
+     */
+    private void constructDebugLines(List<String> debugLines, Variant vref) {
+
+        if (instance().conf.debug) {
+            StringBuilder sb = new StringBuilder();
+            for (String str : debugLines) {
+                if (sb.length() > 0) {
+                    sb.append(" & ");
+                }
+                sb.append(str);
+            }
+            vref.DEBUG = sb.toString();
+        }
     }
 
     /**
diff --git a/src/main/java/com/astrazeneca/vardict/postprocessmodules/SimplePostProcessModule.java b/src/main/java/com/astrazeneca/vardict/postprocessmodules/SimplePostProcessModule.java
@@ -47,7 +47,7 @@ public void accept(Scope<AlignedVarsData> mapScope) {
                     }
                 }
 
-                if (variantsOnPosition != null && variantsOnPosition.variants.isEmpty()) {
+                if (variantsOnPosition.variants.isEmpty()) {
                     if (!conf.doPileup) {
                         continue;
                     }
@@ -70,6 +70,16 @@ public void accept(Scope<AlignedVarsData> mapScope) {
                                 continue;
                             }
                         }
+                        if (vref.startPosition != position && conf.doPileup && vvar.size() == 1) {
+                            Variant refVar = variantsOnPosition.referenceVariant;
+                            if (refVar == null) {
+                                SimpleOutputVariant outputVariant = new SimpleOutputVariant(refVar, mapScope.region, variantsOnPosition.sv, position);
+                                variantPrinter.print(outputVariant);
+                                refVar = new Variant();
+                            }
+                            refVar.vartype = "";
+                            vrefs.add(refVar);
+                        }
                         vref.vartype = vref.varType();
                         if (!vref.isGoodVar(variantsOnPosition.referenceVariant, vref.vartype, mapScope.splice)) {
                             if (!conf.doPileup) {
diff --git a/testdata/integrationtestcases/Simple;hard_clip_case.fa;hard_clip_next_to_del_test1.bam;test;6674-6824;-f 0.0 -p -r 1.txt b/testdata/integrationtestcases/Simple;hard_clip_case.fa;hard_clip_next_to_del_test1.bam;test;6674-6824;-f 0.0 -p -r 1.txt
@@ -3,6 +3,7 @@ hard_clip_next_to_del_test1	testbed	test	6768	6768	G	G	1	0	1	0	0	0	G/G	0	0;0	20.
 hard_clip_next_to_del_test1	testbed	test	6769	6769	G	G	1	0	1	0	0	0	G/G	0	0;0	19.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	0	0	test:6675-6824		0	0
 hard_clip_next_to_del_test1	testbed	test	6770	6770	G	G	1	0	1	0	0	0	G/G	0	0;0	18.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	0	0	test:6675-6824		0	0
 hard_clip_next_to_del_test1	testbed	test	6771	6787	GCCTCTTGCCTGGCAGC	CCTCTTGCCTGGCAGCA	1	1	0	0	1	0	CCTCTTGCCTGGCAGCACCTCTTGCCTGGCAGC/CCTCTTGCCTGGCAGCA	1.0000	0;0	1.0	0	90.0	0	60.0	2.000	1.0000	0	0	2.000	1	0	1	1	AGCTGGATATTGTGGCTGGG	ACAAAGCGCAGAGGTCCCCT	test:6675-6824	Complex	0	0
+hard_clip_next_to_del_test1	testbed	test	6692	6692	C	C	1	0	1	0	0	0	C/C	0	0;0	1.0	0	90.0	0	60.0	0	0	0	0	0	0	0	0	0	0	0	test:6675-6824		0	0
 hard_clip_next_to_del_test1	testbed	test	6693	6695	CTC	TTCT	1	1	1	0	1	0	CCTC/+5	1.0000	0;0	1.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	TGAGCAGAGGCAAGGCACCC	TTACGGGAAGCCCTTCTGGC	test:6675-6824	Complex	0	0
 hard_clip_next_to_del_test1	testbed	test	6696	6767	TTACGGGAAGCCCTTCTGGCGCTCACTCAGGGCAGCAGCTTCCCAAGCCAAGCCAAGCTGGATATTGTGGCT	GGCGCTTCTGGCGCTCACTCAGGGCAGCAGCTTCCCAAGCCAAGCCAAGCTGGATATTGTGGCTG	1	1	0	0	1	0	-7GGCGCTTCTGGCGCTCACTCAGGGCAGCAGCTTCCCAAGCCAAGCCAAGCTGGATATTGTGGCTGAAGCCCTTCTGGCGCTCACTCAGGGCAGCAGCTTCCCAAGCCAAGCCAAGCTGGATATTGTGGCT/-7GGCGCTTCTGGCGCTCACTCAGGGCAGCAGCTTCCCAAGCCAAGCCAAGCTGGATATTGTGGCTG	1.0000	0;0	5.0	0	90.0	0	60.0	2.000	1.0000	0	0	3.000	1	0	1	1	GCAGAGGCAAGGCACCCCTC	GGGGCCTCTTGCCTGGCAGC	test:6675-6824	Complex	0	0
 
diff --git a/testdata/integrationtestcases/Simple;hard_clip_case.fa;hard_clip_next_to_del_test2.bam;test;6434-6564;-f 0.0 -p -r 1.txt b/testdata/integrationtestcases/Simple;hard_clip_case.fa;hard_clip_next_to_del_test2.bam;test;6434-6564;-f 0.0 -p -r 1.txt
@@ -18,6 +18,7 @@ hard_clip_next_to_del_test2	testbed	test	6497	6497	C	C	1	0	1	0	0	0	C/C	0	0;0	49.
 hard_clip_next_to_del_test2	testbed	test	6498	6498	C	C	1	0	1	0	0	0	C/C	0	0;0	48.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	0	0	test:6434-6564		0	0
 hard_clip_next_to_del_test2	testbed	test	6499	6499	C	C	1	0	1	0	0	0	C/C	0	0;0	47.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	0	0	test:6434-6564		0	0
 hard_clip_next_to_del_test2	testbed	test	6500	6501	CT	T	1	1	0	0	1	0	-1TCCACACTCCACAC/-1TCCACAC	1.0000	0;0	40.0	0	90.0	0	60.0	2.000	1.0000	0	0	5.000	1	0	1	1	CAACGCGTCTCTTTTTCCCC	CCACACAAAATAAATGAAAA	test:6434-6564	Complex	0	0
+hard_clip_next_to_del_test2	testbed	test	6442	6442	A	A	1	0	1	0	0	0	A/A	0	0;0	1.0	0	90.0	0	60.0	0	0	0	0	0	0	0	0	0	0	0	test:6434-6564		0	0
 hard_clip_next_to_del_test2	testbed	test	6443	6447	ATGTG	TTGTGA	1	1	1	0	1	0	AATGTG/+7	1.0000	0;0	1.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	CTTTAGTACACAGTTAATAA	AAGTTGTCCTGTGCCCTCCC	test:6434-6564	Complex	0	0
 hard_clip_next_to_del_test2	testbed	test	6508	6508	A	A	1	0	1	0	0	0	A/A	0	0;0	39.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	0	0	test:6434-6564		0	0
 hard_clip_next_to_del_test2	testbed	test	6509	6509	A	A	1	0	1	0	0	0	A/A	0	0;0	38.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	0	0	test:6434-6564		0	0
@@ -29,5 +30,6 @@ hard_clip_next_to_del_test2	testbed	test	6520	6520	A	A	1	0	1	0	0	0	A/A	0	0;0	27.
 hard_clip_next_to_del_test2	testbed	test	6521	6524	ATTT	TTTA	1	1	0	0	1	0	TTTATTT/TTTA	1.0000	0;0	23.0	0	90.0	0	60.0	2.000	1.0000	0	2	3.000	1	0	1	1	TCCACACAAAATAAATGAAA	AAAAGACAGATTTTTTTTTT	test:6434-6564	Complex	0	0
 hard_clip_next_to_del_test2	testbed	test	6525	6525	A	A	1	0	1	0	0	0	A/A	0	0;0	22.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	0	0	test:6434-6564		0	0
 hard_clip_next_to_del_test2	testbed	test	6526	6526	A	A	1	0	1	0	0	0	A/A	0	0;0	21.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	0	0	test:6434-6564		0	0
+hard_clip_next_to_del_test2	testbed	test	6527	6527	A	A	1	0	1	0	0	0	A/A	0	0;0	20.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	0	0	test:6434-6564		0	0
 hard_clip_next_to_del_test2	testbed	test	6528	6529	AG	GAC	1	1	1	0	1	0	AAG/+4	1.0000	0;0	17.0	0	90.0	0	60.0	2.000	1.0000	0	0	0	0	0	1	1	AAAATAAATGAAAATTTAAA	ACAGATTTTTTTTTTTAGGA	test:6434-6564	Complex	0	0
 
diff --git a/testdata/integrationtestcases/Simple;mm10.fa;T7-75_S10_chr10_DEL.bam;chr10;11291400-11295300;-f 0.001;.txt b/testdata/integrationtestcases/Simple;mm10.fa;T7-75_S10_chr10_DEL.bam;chr10;11291400-11295300;-f 0.001;.txt
@@ -3944,6 +3944,7 @@ T7-75_S10_chr10_DEL	testbed	10	11295243	11295243	T	T	1203	0	632	571	0	0	T/T	0	2;
 T7-75_S10_chr10_DEL	testbed	10	11295236	11295236	C	T	1165	1	621	543	1	0	C/T	0.0009	2;0	60.0	0	12.0	0	255.0	0	0	0	4	1.000	1	7.0	0	1146	GGTATGGCGATTCAGTACTG	ATTTTGTTCTGTTAATGACT	10:11291401-11295300	SNV	0	0
 T7-75_S10_chr10_DEL	testbed	10	11295237	11295237	A	G	1164	3	620	540	2	1	A/G	0.0026	2;2	27.0	1	17.3	1	255.0	0.500	0.0009	0.0009	3	4.000	1	1.3	1	1133	GTATGGCGATTCAGTACTGC	TTTTGTTCTGTTAATGACTG	10:11291401-11295300	SNV	0	0
 T7-75_S10_chr10_DEL	testbed	10	11295237	11295237	A	T	1164	1	620	540	1	0	A/T	0.0009	2;0	64.0	0	12.0	0	255.0	0	0	0	3	4.000	1	3.0	0	1133	GTATGGCGATTCAGTACTGC	TTTTGTTCTGTTAATGACTG	10:11291401-11295300	SNV	0	0
+T7-75_S10_chr10_DEL	testbed	10	11295238	11295238	T	T	1164	0	623	540	0	0	T/T	0	2;0	39.5	1	39.8	1	254.8	57.150	0.9991	0	0	0	0	0.4	1143	1144	0	0	10:11291401-11295300		0	0
 T7-75_S10_chr10_DEL	testbed	10	11295237	11295238	AT	A	1164	1	623	540	1	0	T/-1	0.0009	2;0	42.0	0	41.0	0	255.0	2.000	0.0009	0	3	4.000	1	0	1	1144	GTATGGCGATTCAGTACTGC	TTTGTTCTGTTAATGACTGG	10:11291401-11295300	Deletion	0	0
 T7-75_S10_chr10_DEL	testbed	10	11295239	11295239	T	T	1165	0	624	541	0	0	T/T	0	2;0	39.6	1	39.8	1	254.8	54.476	1.0000	0	0	0	0	0.4	1144	1144	0	0	10:11291401-11295300		0	0
 T7-75_S10_chr10_DEL	testbed	10	11295232	11295232	A	C	1093	1	614	477	1	0	A/C	0.0009	2;0	21.0	0	12.0	0	255.0	0	0	0	0	1.000	1	2.0	0	1064	TCTAGGTATGGCGATTCAGT	CTGCATTTTGTTCTGTTAAT	10:11291401-11295300	SNV	0	0

Original file line number	Diff line number	Diff line change
`@@ -47,7 +47,7 @@ public void accept(Scope<AlignedVarsData> mapScope) {`
`47`	`47`	`}`
`48`	`48`	`}`
`49`	`49`
`50`		`- if (variantsOnPosition != null && variantsOnPosition.variants.isEmpty()) {`
	`50`	`+ if (variantsOnPosition.variants.isEmpty()) {`
`51`	`51`	`if (!conf.doPileup) {`
`52`	`52`	`continue;`
`53`	`53`	`}`
`@@ -70,6 +70,16 @@ public void accept(Scope<AlignedVarsData> mapScope) {`
`70`	`70`	`continue;`
`71`	`71`	`}`
`72`	`72`	`}`
	`73`	`+ if (vref.startPosition != position && conf.doPileup && vvar.size() == 1) {`
	`74`	`+ Variant refVar = variantsOnPosition.referenceVariant;`
	`75`	`+ if (refVar == null) {`
	`76`	`+ SimpleOutputVariant outputVariant = new SimpleOutputVariant(refVar, mapScope.region, variantsOnPosition.sv, position);`
	`77`	`+ variantPrinter.print(outputVariant);`
	`78`	`+ refVar = new Variant();`
	`79`	`+ }`
	`80`	`+ refVar.vartype = "";`
	`81`	`+ vrefs.add(refVar);`
	`82`	`+ }`
`73`	`83`	`vref.vartype = vref.varType();`
`74`	`84`	`if (!vref.isGoodVar(variantsOnPosition.referenceVariant, vref.vartype, mapScope.splice)) {`
`75`	`85`	`if (!conf.doPileup) {`