FilterConsensus implements Callable<Integer> { @Override public Integer call() throws Exception { PrintWriter writer = null; try { writer = writer(outputFile); Map<Integer, Allele> exons = readBedFile(bedFile); ListMultimap<Integer, Allele> regions = ArrayListMultimap.create(); for (SAMRecord record : new SAMFileReader(bamFile)) { List<Edit> edits = cigarToEditList(record); String referenceName = record.getReferenceName(); int start = record.getAlignmentStart(); int end = record.getAlignmentEnd(); Locus alignment = new Locus(referenceName, start, end); SymbolList sequence = DNATools.createDNA(record.getReadString()); for (Edit edit : edits) { sequence.edit(edit); } String name = record.getReadName(); Allele contig = Allele.builder() .withContig(referenceName) .withStart(start) .withEnd(end) .withSequence(sequence) .build(); for (Map.Entry<Integer, Allele> entry : exons.entrySet()) { int index = entry.getKey(); Allele exon = entry.getValue(); String chr = exon.getContig(); int min = exon.getMin(); int max = exon.getMax(); String range = chr + ":" + min + "-" + max; if (alignment.overlaps(exon)) { Allele xover = exon.doubleCrossover(contig); Allele clipped = xover.leftHardClip("-").rightHardClip("-"); clipped.setName(">" + name + "|gene=" + gene + "|exon=" + index + "|location=" + range + "|" + (max - min)); regions.put(index, clipped); int offset = 0; for (Edit edit : edits) { if (edit.replacement.equals(SymbolList.EMPTY_LIST) && exon.contains(edit.pos)) { exon.sequence.edit(new Edit(edit.pos + offset, 0, edit.replacement)); offset += edit.length; } } } } } Map<String, ListMultimap<Integer, Allele>> contigs = new HashMap<String, ListMultimap<Integer, Allele>>(); for (int index : regions.keySet()) { for (Allele allele : regions.get(index)) { int sequenceLength = allele.sequence.seqString().length(); List<String> fields = Splitter.on("|").splitToList(allele.getName()); Double locusLength = Double.parseDouble(fields.get(fields.size() - 1)); if (sequenceLength / locusLength >= minimumBreadth) { if (!contigs.containsKey(fields.get(0))) { contigs.put(fields.get(0), ArrayListMultimap.<Integer, Allele>create()); } contigs.get(fields.get(0)).put(index, allele); if (!cdna) { writer.println(allele.getName() + "|" + sequenceLength); writer.println(allele.sequence.seqString().toUpperCase()); } } } } if (cdna) { Map<String, String> cdnas = new HashMap<String, String>(); for (String contig : contigs.keySet()) { StringBuilder sb = new StringBuilder(); for (int index : contigs.get(contig).keySet()) { List<Allele> list = contigs.get(contig).get(index); Collections.sort(list, new Comparator<Allele>() { @Override public int compare(final Allele first, final Allele second) { return first.sequence.seqString().length() - second.sequence.seqString().length(); } }); Allele best = list.get(0); sb.append(best.sequence.seqString()); } String cdnaSequence = sb.toString(); if (removeGaps) { cdnaSequence = cdnaSequence.replaceAll("-", ""); } if (!(gene.equals("HLA-A") || gene.equals("HLA-DPB1"))) { cdnaSequence = DNATools.reverseComplement(DNATools.createDNA(cdnaSequence)).seqString(); } cdnas.put(contig, cdnaSequence.toUpperCase()); } List<Map.Entry<String, String>> list = new ArrayList<Map.Entry<String, String>>(cdnas.entrySet()); Collections.sort(list, new Comparator<Map.Entry<String, String>>() { @Override public int compare(final Map.Entry<String, String> first, final Map.Entry<String, String> second) { return second.getValue().length() - first.getValue().length(); } }); for (Map.Entry<String, String> entry : list.subList(0, expectedPloidy > list.size() ? list.size() : expectedPloidy)) { writer.println(entry.getKey() + "\n" + entry.getValue()); } } return 0; } finally { try { writer.close(); } catch (Exception e) { } } } FilterConsensus(final File bamFile,
                           final File bedFile,
                           final File outputFile,
                           final String gene,
                           final boolean cdna,
                           final boolean removeGaps,
                           final double minimumBreadth,
                           final int expectedPloidy); @Override Integer call(); static void main(final String[] args);  }
@Test public void testCall() throws Exception { copyResource("hla-a.bam", bamFile); copyResource("hla-a.bed", bedFile); new FilterConsensus(bamFile, bedFile, outputFile, gene, cdna, removeGaps, minimumBreadth, expectedPloidy).call(); }
@Test public void testCallPloidyDefault() throws Exception { Files.write(Resources.toByteArray(getClass().getResource("testfile2.bam")), bamFile); Files.write(Resources.toByteArray(getClass().getResource("hla-a.bed")), bedFile); new FilterConsensus(bamFile, bedFile, outputFile, gene, cdna, removeGaps, minimumBreadth, expectedPloidy).call(); assertEquals(4, countLines(outputFile)); }
@Test public void testCallPloidyThree() throws Exception { Files.write(Resources.toByteArray(getClass().getResource("testfile2.bam")), bamFile); Files.write(Resources.toByteArray(getClass().getResource("hla-a.bed")), bedFile); new FilterConsensus(bamFile, bedFile, outputFile, gene, cdna, removeGaps, minimumBreadth, 3).call(); assertEquals(6, countLines(outputFile)); }
@Test public void testCallKirExons() throws Exception { Files.write(Resources.toByteArray(getClass().getResource("2DL1_0020101.bwa.sorted.bam")), bamFile); Files.write(Resources.toByteArray(getClass().getResource("kir-2dl1.exons.bed")), bedFile); cdna = false; removeGaps = false; new FilterConsensus(bamFile, bedFile, outputFile, gene, cdna, removeGaps, minimumBreadth, expectedPloidy).call(); assertEquals(16, countLines(outputFile)); }