Skip to content

Commit 101dbcc

Browse files
bbimberbbimber
committed
Debug merge lofreq VCFs
1 parent c2147e2 commit 101dbcc

File tree

1 file changed

+24
-25
lines changed

1 file changed

+24
-25
lines changed

‎SequenceAnalysis/src/org/labkey/sequenceanalysis/run/analysis/MergeLoFreqVcfHandler.java‎

Lines changed: 24 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -110,8 +110,8 @@ private class SiteAndAlleles
110110
private final int _start;
111111
private final Allele _ref;
112112

113-
Map<String, Map<Allele, Allele>> _encounteredAlleles = new HashMap<>();
114-
List<Allele> _alternates = new ArrayList<>();
113+
Map<String, Map<Allele, String>> _encounteredAlleles = new HashMap<>();
114+
List<String> _alternates = new ArrayList<>();
115115

116116
public SiteAndAlleles(String contig, int start, Allele ref)
117117
{
@@ -120,17 +120,17 @@ public SiteAndAlleles(String contig, int start, Allele ref)
120120
_ref = ref;
121121
}
122122

123-
public Allele getRenamedAllele(VariantContext vc, Allele toCheck)
123+
public String getRenamedAllele(VariantContext vc, Allele toCheck)
124124
{
125-
Map<Allele, Allele> map = _encounteredAlleles.get(getEncounteredKey(vc.getStart(), vc.getReference()));
125+
Map<Allele, String> map = _encounteredAlleles.get(getEncounteredKey(vc.getStart(), vc.getReference()));
126126
if (map == null)
127127
{
128128
throw new IllegalArgumentException("Ref not found at pos " + _start + ": " + vc.toStringWithoutGenotypes() + ", existing: " + StringUtils.join(_encounteredAlleles.keySet(), ";"));
129129
}
130130

131131
if (!map.containsKey(toCheck))
132132
{
133-
throw new IllegalArgumentException("Allele not found at pos " + _start + ": " + vc.toStringWithoutGenotypes() + " for allele: " + toCheck.getBaseString() + ", existing: " + map.keySet().stream().map(x -> x.getBaseString() + "-" + map.get(x).getBaseString()).collect(Collectors.joining(";")));
133+
throw new IllegalArgumentException("Allele not found at pos " + _start + ": " + vc.toStringWithoutGenotypes() + " for allele: " + toCheck.getBaseString() + ", existing: " + map.keySet().stream().map(x -> x.getBaseString() + "-" + map.get(x)).collect(Collectors.joining(";")));
134134
}
135135

136136
return map.get(toCheck);
@@ -143,12 +143,12 @@ private String getEncounteredKey(int start, Allele ref)
143143

144144
public void addSite(VariantContext vc, Logger log)
145145
{
146-
Map<Allele, Allele> alleles = _encounteredAlleles.getOrDefault(getEncounteredKey(vc.getStart(), vc.getReference()), new HashMap<>());
146+
Map<Allele, String> alleles = _encounteredAlleles.getOrDefault(getEncounteredKey(vc.getStart(), vc.getReference()), new HashMap<>());
147147
vc.getAlternateAlleles().forEach(a -> {
148-
Allele translated = extractAlleleForPosition(vc, a, log);
148+
String translated = extractAlleleForPosition(vc, a, log);
149149
if (alleles.containsKey(a) && alleles.get(a) != null && !alleles.get(a).equals(translated))
150150
{
151-
throw new IllegalArgumentException("Translated allele does not match at pos " + _start + ": " + vc.toStringWithoutGenotypes() + " for allele: " + a.getBaseString() + ", existing: " + alleles.keySet().stream().map(x -> x.getBaseString() + "-" + alleles.get(x).getBaseString()).collect(Collectors.joining(";")));
151+
throw new IllegalArgumentException("Translated allele does not match at pos " + _start + ": " + vc.toStringWithoutGenotypes() + " for allele: " + a.getBaseString() + ", existing: " + alleles.keySet().stream().map(x -> x.getBaseString() + "-" + alleles.get(x)).collect(Collectors.joining(";")));
152152
}
153153

154154
alleles.put(a, translated);
@@ -159,7 +159,7 @@ public void addSite(VariantContext vc, Logger log)
159159
_encounteredAlleles.put(getEncounteredKey(vc.getStart(), vc.getReference()), alleles);
160160
}
161161

162-
protected Allele extractAlleleForPosition(VariantContext vc, Allele a, Logger log)
162+
protected String extractAlleleForPosition(VariantContext vc, Allele a, Logger log)
163163
{
164164
int offset = _start - vc.getStart();
165165
if (offset < 0)
@@ -172,17 +172,17 @@ protected Allele extractAlleleForPosition(VariantContext vc, Allele a, Logger lo
172172
//deletion
173173
if (a.length() <= offset)
174174
{
175-
return Allele.SPAN_DEL;
175+
return Allele.SPAN_DEL.getBaseString();
176176
}
177177
else if (offset == 0 && a.length() == 1)
178178
{
179-
return a;
179+
return a.getBaseString();
180180
}
181181
else
182182
{
183183
Allele ret = Allele.create(a.getBaseString().substring(offset, offset + 1));
184184

185-
return _ref.equals(ret, true) ? null: ret;
185+
return _ref.equals(ret, true) ? null: ret.getBaseString();
186186
}
187187
}
188188
}
@@ -294,7 +294,7 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
294294
int idx = 0;
295295
try (CSVWriter writer = new CSVWriter(IOUtil.openFileForBufferedUtf8Writing(output), '\t', CSVWriter.NO_QUOTE_CHARACTER))
296296
{
297-
writer.writeNext(new String[]{"OutputFileId", "ReadsetId", "Contig", "Start", "End", "Ref", "AltAlleles", "Depth", "RefAF", "AltAFs", "NonRefCount", "AltCounts"});
297+
writer.writeNext(new String[]{"ReadsetName", "OutputFileId", "ReadsetId", "Contig", "Start", "End", "Ref", "AltAlleles", "Depth", "RefAF", "AltAFs", "NonRefCount", "AltCounts"});
298298

299299
for (Pair<String, Integer> site : whitelistSites)
300300
{
@@ -306,6 +306,7 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
306306
try (CloseableIterator<VariantContext> it = reader.query(site.getLeft(), site.getRight(), site.getRight()))
307307
{
308308
List<String> line = new ArrayList<>();
309+
line.add(ctx.getSequenceSupport().getCachedReadset(so.getReadset()).getName());
309310
line.add(String.valueOf(so.getRowid()));
310311
line.add(String.valueOf(so.getReadset()));
311312

@@ -321,7 +322,7 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
321322
line.add(String.valueOf(site.getRight() + siteDef._ref.length() - 1));
322323

323324
line.add(siteDef._ref.getBaseString());
324-
line.add(siteDef._alternates.stream().map(Allele::getBaseString).collect(Collectors.joining(";")));
325+
line.add(StringUtils.join(siteDef._alternates, ";"));
325326

326327
if (!it.hasNext())
327328
{
@@ -353,8 +354,8 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
353354
Integer siteDepth = null;
354355
Double totalAltAf = 0.0;
355356
int totalAltDepth = 0;
356-
Map<Allele, Double> alleleToAf = new HashMap<>();
357-
Map<Allele, Integer> alleleToDp = new HashMap<>();
357+
Map<String, Double> alleleToAf = new HashMap<>();
358+
Map<String, Integer> alleleToDp = new HashMap<>();
358359

359360
while (it.hasNext())
360361
{
@@ -390,23 +391,21 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
390391
if (siteDepth < minDepth)
391392
{
392393
vc.getAlternateAlleles().forEach(a -> {
393-
Allele translatedAllele = siteDef.getRenamedAllele(vc, a);
394+
String translatedAllele = siteDef.getRenamedAllele(vc, a);
394395
if (translatedAllele != null)
395396
{
396-
double val = alleleToAf.getOrDefault(a, NO_DATA_VAL);
397-
alleleToAf.put(a, val);
397+
double val = alleleToAf.getOrDefault(translatedAllele, NO_DATA_VAL);
398+
alleleToAf.put(translatedAllele, val);
398399

399-
int val1 = alleleToDp.getOrDefault(a, (int)NO_DATA_VAL);
400-
alleleToDp.put(a, val1);
400+
int val1 = alleleToDp.getOrDefault(translatedAllele, (int)NO_DATA_VAL);
401+
alleleToDp.put(translatedAllele, val1);
401402
}
402403
});
403404
}
404405
else
405406
{
406407
Double af = vc.getAttributeAsDouble("AF", 0.0);
407-
408-
Allele a = vc.getAlternateAlleles().get(0);
409-
a = siteDef.getRenamedAllele(vc, a);
408+
String a = siteDef.getRenamedAllele(vc, vc.getAlternateAlleles().get(0));
410409
if (a != null)
411410
{
412411
totalAltAf += af;
@@ -440,7 +439,7 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
440439
//Add AFs in order:
441440
List<Object> toAdd = new ArrayList<>();
442441
List<Object> toAddDp = new ArrayList<>();
443-
for (Allele a : siteDef._alternates)
442+
for (String a : siteDef._alternates)
444443
{
445444
double af = alleleToAf.getOrDefault(a, 0.0);
446445
toAdd.add(af == NO_DATA_VAL ? "ND" : af);

0 commit comments

Comments
 (0)