Examples with VDJCAlignmentsWriterI com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriterI used on opensource projects

Example 1 with VDJCAlignmentsWriterI

use of com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriterI in project mixcr by milaboratory.

the class ActionAlignmentsDiff method go.

@Override
public void go(ActionHelper actionHelper) throws Exception {
    try (VDJCAlignmentsReader reader1 = new VDJCAlignmentsReader(parameters.get1());
        VDJCAlignmentsReader reader2 = new VDJCAlignmentsReader(parameters.get2());
        VDJCAlignmentsWriterI only1 = parameters.onlyFirst == null ? VDJCAlignmentsWriterI.DummyWriter.INSTANCE : new VDJCAlignmentsWriter(parameters.onlyFirst);
        VDJCAlignmentsWriterI only2 = parameters.onlySecond == null ? VDJCAlignmentsWriterI.DummyWriter.INSTANCE : new VDJCAlignmentsWriter(parameters.onlySecond);
        VDJCAlignmentsWriterI diff1 = parameters.diff1 == null ? VDJCAlignmentsWriterI.DummyWriter.INSTANCE : new VDJCAlignmentsWriter(parameters.diff1);
        VDJCAlignmentsWriterI diff2 = parameters.diff1 == null ? VDJCAlignmentsWriterI.DummyWriter.INSTANCE : new VDJCAlignmentsWriter(parameters.diff2);
        PrintStream report = parameters.report().equals(".") ? System.out : new PrintStream(new FileOutputStream(parameters.report()))) {
        if (reader1.getNumberOfReads() > reader2.getNumberOfReads())
            SmartProgressReporter.startProgressReport("Analyzing diff", reader1);
        else
            SmartProgressReporter.startProgressReport("Analyzing diff", reader2);
        long same = 0, onlyIn1 = 0, onlyIn2 = 0, diffFeature = 0, justDiff = 0;
        long[] diffHits = new long[GeneType.NUMBER_OF_TYPES];
        only1.header(reader1.getParameters(), reader1.getUsedGenes());
        diff1.header(reader1.getParameters(), reader1.getUsedGenes());
        only2.header(reader2.getParameters(), reader2.getUsedGenes());
        diff2.header(reader2.getParameters(), reader2.getUsedGenes());
        VDJCAlignmentsDifferenceReader diffReader = new VDJCAlignmentsDifferenceReader(reader1, reader2, parameters.getFeature(), parameters.hitsCompareLevel);
        for (VDJCAlignmentsDifferenceReader.Diff diff : CUtils.it(diffReader)) {
            switch(diff.status) {
                case AlignmentsAreSame:
                    ++same;
                    break;
                case AlignmentPresentOnlyInFirst:
                    ++onlyIn1;
                    only1.write(diff.first);
                    break;
                case AlignmentPresentOnlyInSecond:
                    ++onlyIn2;
                    only2.write(diff.second);
                    break;
                case AlignmentsAreDifferent:
                    ++justDiff;
                    diff1.write(diff.first);
                    diff2.write(diff.second);
                    if (diff.reason.diffGeneFeature)
                        ++diffFeature;
                    for (Map.Entry<GeneType, Boolean> e : diff.reason.diffHits.entrySet()) if (e.getValue())
                        ++diffHits[e.getKey().ordinal()];
            }
        }
        only1.setNumberOfProcessedReads(onlyIn1);
        only2.setNumberOfProcessedReads(onlyIn2);
        diff1.setNumberOfProcessedReads(justDiff);
        diff2.setNumberOfProcessedReads(justDiff);
        report.println("First  file: " + parameters.get1());
        report.println("Second file: " + parameters.get2());
        report.println("Completely same reads: " + same);
        report.println("Aligned reads present only in the FIRST  file: " + onlyIn1 + " (" + Util.PERCENT_FORMAT.format(100. * onlyIn1 / reader1.getNumberOfReads()) + ")%");
        report.println("Aligned reads present only in the SECOND file: " + onlyIn2 + " (" + Util.PERCENT_FORMAT.format(100. * onlyIn2 / reader2.getNumberOfReads()) + ")%");
        report.println("Total number of different reads: " + justDiff);
        report.println("Reads with not same " + parameters.geneFeatureToMatch + ": " + diffFeature);
        for (GeneType geneType : GeneType.VDJC_REFERENCE) report.println("Reads with not same " + geneType.name() + " hits: " + diffHits[geneType.ordinal()]);
    }
}