Examples with Counter com.github.lindenb.jvarkit.util.Counter used on opensource projects

Example 16 with Counter

use of com.github.lindenb.jvarkit.util.Counter in project jvarkit by lindenb.

the class XContaminations method doWork.

@Override
public int doWork(final List<String> args) {
    long last_save_ms = System.currentTimeMillis();
    if (this.output_as_vcf && !this.use_only_sample_name) {
        LOG.error("cannot write vcf if --sample is not set");
        return -1;
    }
    if (args.size() < 2) {
        LOG.error("Illegal Number of args");
        return -1;
    }
    final Set<File> bamFiles = IOUtils.unrollFiles(args.subList(1, args.size())).stream().map(S -> new File(S)).collect(Collectors.toSet());
    if (bamFiles.isEmpty()) {
        LOG.error("Undefined BAM file(s)");
        return -1;
    }
    SAMRecordIterator iter = null;
    VcfIterator in = null;
    Map<String, SamReader> sample2samReader = new HashMap<>();
    VariantContextWriter vcfw = null;
    try {
        final SamReaderFactory srf = super.createSamReaderFactory();
        if (args.get(0).equals("-")) {
            in = super.openVcfIterator(null);
        } else {
            in = super.openVcfIterator(args.get(0));
        }
        VCFHeader vcfHeader = in.getHeader();
        final SAMSequenceDictionary dict1 = vcfHeader.getSequenceDictionary();
        if (dict1 == null) {
            LOG.error(JvarkitException.VcfDictionaryMissing.getMessage(args.get(0)));
            return -1;
        }
        final Set<String> sampleNames = new HashSet<>(vcfHeader.getSampleNamesInOrder());
        if (sampleNames.isEmpty()) {
            LOG.error("VCF contains no sample");
            return -1;
        }
        for (final File bamFile : bamFiles) {
            LOG.info("Opening " + bamFile);
            final SamReader samReader = srf.open(bamFile);
            final SAMFileHeader samHeader = samReader.getFileHeader();
            final SAMSequenceDictionary dict2 = samHeader.getSequenceDictionary();
            if (dict2 == null) {
                samReader.close();
                LOG.error(JvarkitException.BamDictionaryMissing.getMessage(bamFile.getPath()));
                return -1;
            }
            if (!SequenceUtil.areSequenceDictionariesEqual(dict1, dict2)) {
                samReader.close();
                LOG.error(JvarkitException.DictionariesAreNotTheSame.getMessage(dict1, dict2));
                return -1;
            }
            if (!samReader.hasIndex()) {
                samReader.close();
                LOG.error("sam is not indexed : " + bamFile);
                return -1;
            }
            String sampleName = null;
            for (final SAMReadGroupRecord rgr : samHeader.getReadGroups()) {
                final String s = rgr.getSample();
                if (StringUtil.isBlank(s))
                    continue;
                if (sampleName == null) {
                    sampleName = s;
                } else if (!sampleName.equals(s)) {
                    samReader.close();
                    LOG.error("Cannot handle more than one sample/bam  " + bamFile + " " + sampleName);
                    return -1;
                }
            }
            if (sampleName == null) {
                samReader.close();
                LOG.error("No sample in " + bamFile);
                // skip this bam
                continue;
            }
            if (!sampleNames.contains(sampleName)) {
                samReader.close();
                LOG.error("Not in VCF header: sample " + sampleName + " " + bamFile);
                // skip this bam
                continue;
            }
            if (sample2samReader.containsKey(sampleName)) {
                samReader.close();
                LOG.error("Cannot handle more than one bam/sample: " + bamFile + " " + sampleName);
                return -1;
            }
            sample2samReader.put(sampleName, samReader);
        }
        if (sample2samReader.size() < 2) {
            LOG.error("Not engough BAM/samples. Expected at least two valid BAMs");
            return -1;
        }
        sampleNames.retainAll(sample2samReader.keySet());
        /* create a VCF is VCF output asked */
        final List<SamplePair> sampleListForVcf;
        if (this.output_as_vcf) {
            vcfw = super.openVariantContextWriter(outputFile);
            final Set<VCFHeaderLine> metaData = new HashSet<>();
            metaData.add(new VCFFormatHeaderLine("S1S1", 1, VCFHeaderLineType.Integer, "reads sample 1 supporting sample 1"));
            metaData.add(new VCFFormatHeaderLine("S1S2", 1, VCFHeaderLineType.Integer, "reads sample 1 supporting sample 2"));
            metaData.add(new VCFFormatHeaderLine("S1SO", 1, VCFHeaderLineType.Integer, "reads sample 1 supporting others"));
            metaData.add(new VCFFormatHeaderLine("S2S1", 1, VCFHeaderLineType.Integer, "reads sample 2 supporting sample 1"));
            metaData.add(new VCFFormatHeaderLine("S2S2", 1, VCFHeaderLineType.Integer, "reads sample 2 supporting sample 2"));
            metaData.add(new VCFFormatHeaderLine("S2SO", 1, VCFHeaderLineType.Integer, "reads sample 2 supporting others"));
            metaData.add(new VCFFormatHeaderLine("FR", 1, VCFHeaderLineType.Float, "Fraction. '-1' for unavailable."));
            metaData.add(new VCFFormatHeaderLine("S1A", 1, VCFHeaderLineType.Character, "sample 1 allele"));
            metaData.add(new VCFFormatHeaderLine("S2A", 1, VCFHeaderLineType.Character, "sample 2 allele"));
            metaData.add(new VCFFilterHeaderLine("XCONTAMINATION", "Fraction test is > " + fraction_treshold));
            metaData.add(new VCFFilterHeaderLine("BADSAMPLES", "At least one pair of genotype fails the 'LE' test"));
            metaData.add(new VCFInfoHeaderLine("LE", 1, VCFHeaderLineType.Integer, "number of pair of genotypes having (S1S1<=S1S2 or S2S2<=S2S1)."));
            metaData.add(new VCFInfoHeaderLine("BADSAMPLES", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, "Samples founds failing the 'LE' test"));
            sampleListForVcf = new ArrayList<>();
            final List<String> sampleList = new ArrayList<>(sampleNames);
            for (int x = 0; x + 1 < sampleList.size(); ++x) {
                for (int y = x + 1; y < sampleList.size(); ++y) {
                    sampleListForVcf.add(new SamplePair(new SimpleSampleIdenfifier(sampleList.get(x)), new SimpleSampleIdenfifier(sampleList.get(y))));
                }
            }
            final VCFHeader header2 = new VCFHeader(metaData, sampleListForVcf.stream().map(V -> V.getLabel()).sorted().collect(Collectors.toList()));
            header2.setSequenceDictionary(dict1);
            vcfw.writeHeader(header2);
        } else {
            vcfw = null;
            sampleListForVcf = null;
        }
        final Map<SamplePair, SampleAlleles> contaminationTable = new HashMap<>();
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(dict1).logger(LOG);
        while (in.hasNext()) {
            final VariantContext ctx = progress.watch(in.next());
            if (!ctx.isSNP() || ctx.isFiltered() || !ctx.isBiallelic() || ctx.isSymbolic() || !this.variantFilter.test(ctx)) {
                continue;
            }
            int count_homref = 0;
            int count_homvar = 0;
            int count_het = 0;
            final Map<String, Genotype> sample2gt = new HashMap<>();
            for (int gidx = 0; gidx < ctx.getNSamples(); ++gidx) {
                final Genotype G = ctx.getGenotype(gidx);
                if (!G.isCalled())
                    continue;
                if (G.isHet()) {
                    // here because in use_singleton we must be sure that there is only one hom_var
                    count_het++;
                    if (this.use_singleton && count_het > 0)
                        break;
                } else if (G.isHomVar()) {
                    // here because in use_singleton we must be sure that there is only one hom_var
                    count_homvar++;
                    if (this.use_singleton && count_homvar > 1)
                        break;
                }
                if (G.isFiltered())
                    continue;
                if (!sample2samReader.containsKey(G.getSampleName()))
                    continue;
                if (!sampleNames.contains(G.getSampleName()))
                    continue;
                if (!this.genotypeFilter.test(ctx, G))
                    continue;
                sample2gt.put(G.getSampleName(), G);
            }
            if (this.use_singleton && count_het > 0)
                continue;
            if (this.use_singleton && count_homvar > 1)
                continue;
            if (sample2gt.size() < 2)
                continue;
            // reset and recount
            count_homref = 0;
            count_homvar = 0;
            count_het = 0;
            for (final String sampleName : sample2gt.keySet()) {
                final Genotype G = ctx.getGenotype(sampleName);
                switch(G.getType()) {
                    case HOM_REF:
                        count_homref++;
                        break;
                    case HOM_VAR:
                        count_homvar++;
                        break;
                    case HET:
                        count_het++;
                        break;
                    default:
                        break;
                }
            }
            // singleton check
            if (this.use_singleton && (count_het > 0 || count_homvar != 1)) {
                continue;
            }
            // at least one HOM_REF and one HOM_VAR
            if (count_homref == 0)
                continue;
            if (count_homvar == 0)
                continue;
            final Map<SampleIdentifier, Counter<Character>> sample_identifier_2allelesCount = new HashMap<>();
            /* scan Reads for those Genotype/Samples */
            for (final String sampleName : sample2gt.keySet()) {
                if (!sample2samReader.containsKey(sampleName))
                    continue;
                // sample name is not in vcf header
                final SamReader samReader = sample2samReader.get(sampleName);
                if (samReader == null)
                    continue;
                final Genotype genotype = sample2gt.get(sampleName);
                if (genotype == null)
                    continue;
                iter = samReader.query(ctx.getContig(), ctx.getStart(), ctx.getEnd(), false);
                while (iter.hasNext()) {
                    final SAMRecord record = iter.next();
                    if (record.getEnd() < ctx.getStart())
                        continue;
                    if (ctx.getEnd() < record.getStart())
                        continue;
                    if (record.getReadUnmappedFlag())
                        continue;
                    if (this.filter.filterOut(record))
                        continue;
                    final SAMReadGroupRecord srgr = record.getReadGroup();
                    // not current sample
                    if (srgr == null)
                        continue;
                    if (!sampleName.equals(srgr.getSample()))
                        continue;
                    final Cigar cigar = record.getCigar();
                    if (cigar == null || cigar.isEmpty())
                        continue;
                    byte[] readSeq = record.getReadBases();
                    if (readSeq == null || readSeq.length == 0)
                        continue;
                    int readPos = record.getReadPositionAtReferencePosition(ctx.getStart());
                    if (readPos < 1)
                        continue;
                    readPos--;
                    if (readPos >= readSeq.length)
                        continue;
                    final char base = Character.toUpperCase((char) readSeq[readPos]);
                    if (base == 'N')
                        continue;
                    final SampleIdentifier sampleIdentifier;
                    if (this.use_only_sample_name) {
                        sampleIdentifier = new SimpleSampleIdenfifier(sampleName);
                    } else {
                        final ShortReadName readName = ShortReadName.parse(record);
                        if (!readName.isValid()) {
                            LOG.info("No a valid read name " + record.getReadName());
                            continue;
                        }
                        sampleIdentifier = new SequencerFlowCellRunLaneSample(readName, sampleName);
                    }
                    Counter<Character> sampleAlleles = sample_identifier_2allelesCount.get(sampleIdentifier);
                    if (sampleAlleles == null) {
                        sampleAlleles = new Counter<Character>();
                        sample_identifier_2allelesCount.put(sampleIdentifier, sampleAlleles);
                    }
                    sampleAlleles.incr(base);
                }
                iter.close();
                iter = null;
            }
            /* end scan reads for this sample */
            /* sum-up data for this SNP */
            final VariantContextBuilder vcb;
            final List<Genotype> genotypeList;
            if (this.output_as_vcf) {
                vcb = new VariantContextBuilder(args.get(0), ctx.getContig(), ctx.getStart(), ctx.getEnd(), ctx.getAlleles());
                if (ctx.hasID())
                    vcb.id(ctx.getID());
                genotypeList = new ArrayList<>();
            } else {
                vcb = null;
                genotypeList = null;
            }
            for (final String sample1 : sample2gt.keySet()) {
                final Genotype g1 = sample2gt.get(sample1);
                final char a1 = g1.getAllele(0).getBaseString().charAt(0);
                for (final String sample2 : sample2gt.keySet()) {
                    if (sample1.compareTo(sample2) >= 0)
                        continue;
                    final Genotype g2 = sample2gt.get(sample2);
                    if (g2.sameGenotype(g1))
                        continue;
                    final char a2 = g2.getAllele(0).getBaseString().charAt(0);
                    for (final SampleIdentifier sfcr1 : sample_identifier_2allelesCount.keySet()) {
                        if (!sfcr1.getSampleName().equals(sample1))
                            continue;
                        final Counter<Character> counter1 = sample_identifier_2allelesCount.get(sfcr1);
                        if (counter1 == null)
                            continue;
                        for (final SampleIdentifier sfcr2 : sample_identifier_2allelesCount.keySet()) {
                            if (!sfcr2.getSampleName().equals(sample2))
                                continue;
                            final SamplePair samplePair = new SamplePair(sfcr1, sfcr2);
                            final Counter<Character> counter2 = sample_identifier_2allelesCount.get(sfcr2);
                            if (counter2 == null)
                                continue;
                            SampleAlleles sampleAlleles = contaminationTable.get(samplePair);
                            if (sampleAlleles == null) {
                                sampleAlleles = new SampleAlleles();
                                contaminationTable.put(samplePair, sampleAlleles);
                                if (!this.output_as_vcf && contaminationTable.size() % 10000 == 0)
                                    LOG.info("n(pairs)=" + contaminationTable.size());
                            }
                            sampleAlleles.number_of_comparaisons++;
                            for (final Character allele : counter1.keySet()) {
                                final long n = counter1.count(allele);
                                if (allele.equals(a1)) {
                                    sampleAlleles.reads_sample1_supporting_sample1 += n;
                                } else if (allele.equals(a2)) {
                                    sampleAlleles.reads_sample1_supporting_sample2 += n;
                                } else {
                                    sampleAlleles.reads_sample1_supporting_other += n;
                                }
                            }
                            for (final Character allele : counter2.keySet()) {
                                final long n = counter2.count(allele);
                                if (allele.equals(a2)) {
                                    sampleAlleles.reads_sample2_supporting_sample2 += n;
                                } else if (allele.equals(a1)) {
                                    sampleAlleles.reads_sample2_supporting_sample1 += n;
                                } else {
                                    sampleAlleles.reads_sample2_supporting_other += n;
                                }
                            }
                        }
                    }
                }
            }
            if (this.output_as_vcf) {
                final Set<String> bad_samples = new TreeSet<>();
                boolean fraction_flag = false;
                int num_lt = 0;
                for (final SamplePair samplepair : sampleListForVcf) {
                    final GenotypeBuilder gb = new GenotypeBuilder(samplepair.getLabel());
                    final SampleAlleles sampleAlleles = contaminationTable.get(samplepair);
                    if (sampleAlleles != null) {
                        gb.attribute("S1S1", sampleAlleles.reads_sample1_supporting_sample1);
                        gb.attribute("S1S2", sampleAlleles.reads_sample1_supporting_sample2);
                        gb.attribute("S1SO", sampleAlleles.reads_sample1_supporting_other);
                        gb.attribute("S2S1", sampleAlleles.reads_sample2_supporting_sample1);
                        gb.attribute("S2S2", sampleAlleles.reads_sample2_supporting_sample2);
                        gb.attribute("S2SO", sampleAlleles.reads_sample2_supporting_other);
                        gb.attribute("S1A", sample2gt.get(samplepair.sample1.getSampleName()).getAllele(0).getDisplayString().charAt(0));
                        gb.attribute("S2A", sample2gt.get(samplepair.sample2.getSampleName()).getAllele(0).getDisplayString().charAt(0));
                        final double fraction = sampleAlleles.getFraction();
                        gb.attribute("FR", fraction);
                        if (!this.passFractionTreshold.test(fraction)) {
                            fraction_flag = true;
                        }
                        boolean bad_lt_flag = false;
                        if (sampleAlleles.reads_sample1_supporting_sample1 <= this.fail_factor * sampleAlleles.reads_sample1_supporting_sample2) {
                            bad_samples.add(samplepair.sample1.getSampleName());
                            bad_lt_flag = true;
                        }
                        if (sampleAlleles.reads_sample2_supporting_sample2 <= this.fail_factor * sampleAlleles.reads_sample2_supporting_sample1) {
                            bad_samples.add(samplepair.sample2.getSampleName());
                            bad_lt_flag = true;
                        }
                        if (bad_lt_flag) {
                            num_lt++;
                        }
                    } else {
                        gb.attribute("S1S1", -1);
                        gb.attribute("S1S2", -1);
                        gb.attribute("S1SO", -1);
                        gb.attribute("S2S1", -1);
                        gb.attribute("S2S2", -1);
                        gb.attribute("S2SO", -1);
                        gb.attribute("S1A", '.');
                        gb.attribute("S2A", '.');
                        gb.attribute("FR", -1f);
                    }
                    genotypeList.add(gb.make());
                }
                if (!bad_samples.isEmpty()) {
                    vcb.attribute("BADSAMPLES", new ArrayList<>(bad_samples));
                }
                vcb.attribute("LE", num_lt);
                if (fraction_flag || !bad_samples.isEmpty()) {
                    if (fraction_flag)
                        vcb.filter("XCONTAMINATION");
                    if (!bad_samples.isEmpty())
                        vcb.filter("BADSAMPLES");
                } else {
                    vcb.passFilters();
                }
                vcb.genotypes(genotypeList);
                vcfw.add(vcb.make());
                contaminationTable.clear();
            } else {
                final long now = System.currentTimeMillis();
                if (this.outputFile != null && this.save_every_sec > -1L && last_save_ms + (this.save_every_sec * 1000L) > now) {
                    saveToFile(contaminationTable);
                    last_save_ms = now;
                }
            }
        }
        progress.finish();
        if (this.output_as_vcf) {
            vcfw.close();
            vcfw = null;
        } else {
            saveToFile(contaminationTable);
        }
        return 0;
    } catch (final Exception e) {
        LOG.error(e);
        return -1;
    } finally {
        CloserUtil.close(vcfw);
        CloserUtil.close(in);
        CloserUtil.close(iter);
        for (SamReader samReader : sample2samReader.values()) CloserUtil.close(samReader);
        sample2samReader.clear();
    }
}

Example 17 with Counter

use of com.github.lindenb.jvarkit.util.Counter in project jvarkit by lindenb.

the class CompareBams4 method doWork.

@Override
public int doWork(final List<String> args) {
    PrintWriter out = null;
    final SamReader[] samFileReaders = new SamReader[] { null, null };
    final SAMSequenceDictionary[] dicts = new SAMSequenceDictionary[] { null, null };
    @SuppressWarnings("unchecked") final PeekableIterator<SAMRecord>[] iters = new PeekableIterator[] { null, null };
    final List<List<SAMRecord>> recordLists = Arrays.asList(new ArrayList<SAMRecord>(), new ArrayList<SAMRecord>());
    final Counter<Diff> diffs = new Counter<>();
    final short NM_TAG = SAMTagUtil.getSingleton().NM;
    if (args.size() != 2) {
        LOG.error("Expected two and only two bams please, but got " + args.size());
        return -1;
    }
    try {
        final SamReaderFactory srf = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.LENIENT);
        for (int i = 0; i < args.size() && i < samFileReaders.length; ++i) {
            final File samFile = new File(args.get(i));
            LOG.info("opening " + samFile);
            samFileReaders[i] = srf.open(samFile);
            final SAMFileHeader header = samFileReaders[i].getFileHeader();
            if (header.getSortOrder() != SAMFileHeader.SortOrder.queryname) {
                LOG.error("Expected " + samFile + " to be sorted on " + SAMFileHeader.SortOrder.queryname + " but got " + header.getSortOrder());
                return -1;
            }
            dicts[i] = header.getSequenceDictionary();
            if (dicts[i] == null || dicts[i].isEmpty()) {
                LOG.error("In " + samFile + ": No SAMSequenceDictionary in header.");
                return -1;
            }
            iters[i] = new PeekableIterator<>(samFileReaders[i].iterator());
        }
        if (this.chainFile != null) {
            LOG.info("opening chain file " + this.chainFile + ".");
            this.liftOver = new LiftOver(this.chainFile);
            this.liftOver.setLiftOverMinMatch(this.liftOverMismatch <= 0 ? LiftOver.DEFAULT_LIFTOVER_MINMATCH : this.liftOverMismatch);
            if (!this.disableChainValidation) {
                this.liftOver.validateToSequences(dicts[1]);
            }
        }
        final Comparator<SAMRecord> comparator = this.sortMethod.get();
        for (; ; ) {
            for (int i = 0; i < 2; ++i) {
                if (recordLists.get(i).isEmpty()) {
                    while (iters[i].hasNext()) {
                        final SAMRecord rec = iters[i].peek();
                        if (rec.isSecondaryOrSupplementary()) {
                            iters[i].next();
                            continue;
                        }
                        if (!recordLists.get(i).isEmpty() && comparator.compare(recordLists.get(i).get(0), rec) > 0) {
                            throw new SAMException("Something is wrong in sort order of " + args.get(i) + " : got\n\t" + rec + " " + side(rec) + "\nafter\n\t" + recordLists.get(i).get(0) + " " + side(recordLists.get(i).get(0)) + "\nSee also option (samtools querysort)");
                        } else if (recordLists.get(i).isEmpty() || comparator.compare(recordLists.get(i).get(0), rec) == 0) {
                            recordLists.get(i).add(iters[i].next());
                        } else {
                            break;
                        }
                    }
                }
            }
            if (recordLists.get(0).isEmpty() && recordLists.get(1).isEmpty())
                break;
            if (recordLists.get(0).size() > 1) {
                LOG.warn("size>2 for 1/2:" + recordLists.get(0).get(0).getReadName());
                for (final SAMRecord sr : recordLists.get(0)) {
                    LOG.warn(">> " + sr + " flags:" + sr.getFlags() + " pos:" + sr.getReferenceName() + ":" + sr.getStart());
                }
            }
            if (recordLists.get(1).size() > 1) {
                LOG.warn("size>2 for 2/2:" + recordLists.get(1).get(0).getReadName());
                for (final SAMRecord sr : recordLists.get(1)) {
                    LOG.warn(">> " + sr + " flags:" + sr.getFlags() + " pos:" + sr.getReferenceName() + ":" + sr.getStart());
                }
            }
            final Diff diff = new Diff();
            if ((recordLists.get(0).isEmpty() && !recordLists.get(1).isEmpty()) || (!recordLists.get(0).isEmpty() && !recordLists.get(1).isEmpty() && comparator.compare(recordLists.get(0).get(0), recordLists.get(1).get(0)) > 0)) {
                diff.onlyIn = OnlyIn.ONLY_IN_SECOND;
                recordLists.get(1).clear();
            } else if ((!recordLists.get(0).isEmpty() && recordLists.get(1).isEmpty()) || (!recordLists.get(0).isEmpty() && !recordLists.get(1).isEmpty() && comparator.compare(recordLists.get(0).get(0), recordLists.get(1).get(0)) < 0)) {
                diff.onlyIn = OnlyIn.ONLY_IN_FIRST;
                recordLists.get(0).clear();
            } else {
                final SAMRecord rec0 = recordLists.get(0).get(0);
                final SAMRecord rec1 = recordLists.get(1).get(0);
                final Interval rgn0a = interval(rec0);
                final Interval rgn0b;
                diff.onlyIn = OnlyIn.BOTH;
                diff.diffFlags = new Couple<Integer>(rec0.getFlags(), rec1.getFlags());
                diff.diffFlag = (rec0.getFlags() == rec1.getFlags() ? Flag.Same : Flag.Discordant);
                if (this.liftOver == null) {
                    rgn0b = rgn0a;
                    diff.liftover = LiftOverStatus.NoLiftOver;
                } else if (rgn0a == null) {
                    diff.liftover = LiftOverStatus.SourceUnmapped;
                    rgn0b = rgn0a;
                } else {
                    rgn0b = this.liftOver.liftOver(rgn0a);
                    if (rgn0b == null) {
                        diff.liftover = LiftOverStatus.LiftOverFailed;
                    } else if (dicts[1].getSequence(rgn0b.getContig()) == null) {
                        diff.liftover = LiftOverStatus.DestNotInDict;
                    } else if (rgn0a.getContig().equals(rgn0b.getContig())) {
                        diff.liftover = LiftOverStatus.SameChrom;
                    } else {
                        diff.liftover = LiftOverStatus.DiscordantChrom;
                    }
                }
                final Interval rgn1 = interval(rec1);
                if (rgn0b == null && rgn1 == null) {
                    diff.compareContig = CompareContig.BothUnmapped;
                } else if (rgn0b == null && rgn1 != null) {
                    diff.compareContig = CompareContig.GainMapping;
                    diff.diffChrom = new Couple<String>("*", rgn1.getContig());
                } else if (rgn0b != null && rgn1 == null) {
                    diff.compareContig = CompareContig.LostMapping;
                    diff.diffChrom = new Couple<String>(rgn0b.getContig(), "*");
                } else if (rgn0b.getContig().equals(rgn1.getContig())) {
                    diff.compareContig = CompareContig.SameContig;
                    diff.diffChrom = new Couple<String>(rgn0b.getContig(), rgn1.getContig());
                    final int shift = Math.abs(rgn0b.getStart() - rgn1.getStart());
                    if (shift == 0) {
                        diff.shift = Shift.Zero;
                    } else if (shift <= 10) {
                        diff.shift = Shift.Le10;
                    } else if (shift <= 20) {
                        diff.shift = Shift.Le20;
                    } else if (shift <= 100) {
                        diff.shift = Shift.Le100;
                    } else {
                        diff.shift = Shift.Gt100;
                    }
                } else {
                    if (dicts[1].getSequence(rgn0b.getContig()) == null) {
                        diff.diffChrom = new Couple<String>("?", rgn1.getContig());
                    } else {
                        diff.diffChrom = new Couple<String>(rgn0b.getContig(), rgn1.getContig());
                    }
                    diff.compareContig = CompareContig.DiscordantContig;
                }
                if (rgn0b != null && rgn1 != null) {
                    diff.diffCigarOperations = (rec0.getCigar().numCigarElements() - rec1.getCigar().numCigarElements());
                    final Object nm0 = rec0.getAttribute(NM_TAG);
                    final Object nm1 = rec1.getAttribute(NM_TAG);
                    if (nm0 != null && nm1 != null) {
                        diff.diffNM = (Number.class.cast(nm0).intValue() - Number.class.cast(nm1).intValue());
                    }
                }
                recordLists.get(1).clear();
                recordLists.get(0).clear();
            }
            diffs.incr(diff);
        }
        LOG.info("done");
        final StringBuilder sb = new StringBuilder();
        str(sb, "onlyIn");
        str(sb, "liftover");
        str(sb, "compareContig");
        str(sb, "shift");
        str(sb, "diffCigarOperations");
        str(sb, "diffNM");
        str(sb, "diffFlags");
        str(sb, "diffChroms");
        str(sb, "diffFlag");
        sb.append("Count");
        out = super.openFileOrStdoutAsPrintWriter(outputFile);
        out.println(sb);
        for (final Diff key : diffs.keySet()) {
            out.print(key.toString());
            out.println(diffs.count(key));
        }
        out.flush();
        out.close();
        return RETURN_OK;
    } catch (Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        this.liftOver = null;
        for (int i = 0; i < samFileReaders.length; ++i) {
            CloserUtil.close(iters[i]);
            CloserUtil.close(samFileReaders[i]);
        }
        CloserUtil.close(out);
    }
}