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

Example 6 with CloneSet

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

the class FullSeqAssemblerTest method testLargeCloneNoMismatches.

@Test
public void testLargeCloneNoMismatches() throws Exception {
    MasterSequence master = FullSeqAssemblerTest.masterSeq1WT;
    NSequenceWithQuality seq = new NSequenceWithQuality(master.getRange(-master.vPart + 10, 80), SequenceQuality.GOOD_QUALITY_VALUE);
    RunMiXCR.RunMiXCRAnalysis params0 = new RunMiXCR.RunMiXCRAnalysis(new SingleReadImpl(0, seq, ""));
    params0.cloneAssemblerParameters.setAssemblingFeatures(new GeneFeature[] { GeneFeature.VDJRegion });
    Clone largeClone = RunMiXCR.assemble(RunMiXCR.align(params0)).cloneSet.get(0);
    // ActionExportClonesPretty.outputCompact(System.out, largeClone);
    // System.exit(0);
    Well44497b rnd = new Well44497b(1234567889L);
    int nReads = 100_000;
    int readLength = 75, readGap = 150;
    // slice seq randomly
    PairedRead[] slicedReads = new PairedRead[nReads];
    for (int i = 0; i < nReads; ++i) {
        int r1from = rnd.nextInt(seq.size() - readLength - 1), r1to = r1from + readLength, r2from = r1from + 1 + rnd.nextInt(seq.size() - r1from - readLength - 1), r2to = r2from + readLength;
        assert r2from > r1from;
        slicedReads[i] = new PairedRead(new SingleReadImpl(i, seq.getRange(r1from, r1to), "" + i), new SingleReadImpl(i, seq.getRange(r2from, r2to).getReverseComplement(), "" + i));
    }
    RunMiXCR.RunMiXCRAnalysis params = new RunMiXCR.RunMiXCRAnalysis(slicedReads);
    // params.alignerParameters = VDJCParametersPresets.getByName("rna-seq");
    params.alignerParameters.setSaveOriginalReads(true);
    RunMiXCR.AlignResult align = RunMiXCR.align(params);
    RunMiXCR.AssembleResult assemble = RunMiXCR.assemble(align);
    // for (VDJCAlignments al : align.alignments) {
    // for (int i = 0; i < al.numberOfTargets(); i++) {
    // System.out.println(VDJCAlignmentsFormatter.getTargetAsMultiAlignment(al, i));
    // System.out.println();
    // }
    // System.out.println();
    // System.out.println(" ================================================ ");
    // System.out.println();
    // }
    Assert.assertEquals(1, assemble.cloneSet.size());
    Clone initialClone = assemble.cloneSet.get(0);
    NSequenceWithQuality cdr3 = initialClone.getFeature(GeneFeature.CDR3);
    List<VDJCAlignments> alignments = align.alignments.stream().filter(al -> cdr3.equals(al.getFeature(GeneFeature.CDR3))).collect(Collectors.toList());
    alignments.stream().filter(al -> Arrays.stream(al.getBestHit(GeneType.Variable).getAlignments()).filter(Objects::nonNull).anyMatch(a -> !a.getAbsoluteMutations().isEmpty())).filter(al -> al.getBestHit(GeneType.Variable).getGene().getName().contains("3-74")).forEach(al -> {
        for (int i = 0; i < al.numberOfTargets(); i++) {
            System.out.println(VDJCAlignmentsFormatter.getTargetAsMultiAlignment(al, i));
            System.out.println();
        }
        System.out.println();
        System.out.println(" ================================================ ");
        System.out.println();
    });
    // System.exit(0);
    System.out.println("=> Agg");
    CloneFactory cloneFactory = new CloneFactory(align.parameters.cloneAssemblerParameters.getCloneFactoryParameters(), align.parameters.cloneAssemblerParameters.getAssemblingFeatures(), align.usedGenes, align.parameters.alignerParameters.getFeaturesToAlignMap());
    FullSeqAssembler agg = new FullSeqAssembler(cloneFactory, DEFAULT_PARAMETERS, initialClone, align.parameters.alignerParameters);
    FullSeqAssembler.RawVariantsData prep = agg.calculateRawData(() -> CUtils.asOutputPort(alignments));
    List<Clone> clones = new ArrayList<>(new CloneSet(Arrays.asList(agg.callVariants(prep))).getClones());
    clones.sort(Comparator.comparingDouble(Clone::getCount).reversed());
    for (Clone clone : clones) {
        ActionExportClonesPretty.outputCompact(System.out, clone);
        System.out.println();
        System.out.println(" ================================================ ");
        System.out.println();
    }
}

Example 7 with CloneSet

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

the class RunMiXCR method assemble.

public static AssembleResult assemble(final AlignResult align, boolean close) {
    CloneAssembler assembler = null;
    try {
        RunMiXCRAnalysis parameters = align.parameters;
        assembler = new CloneAssembler(parameters.cloneAssemblerParameters, false, align.usedGenes, align.parameters.alignerParameters);
        CloneAssemblerReport report = new CloneAssemblerReport();
        assembler.setListener(report);
        CloneAssemblerRunner assemblerRunner = new CloneAssemblerRunner(new AlignmentsProvider() {

            @Override
            public OutputPortCloseable<VDJCAlignments> create() {
                return opCloseable(CUtils.asOutputPort(align.alignments));
            }

            @Override
            public long getTotalNumberOfReads() {
                return align.alignments.size();
            }
        }, assembler, parameters.threads);
        // start progress reporting
        SmartProgressReporter.startProgressReport(assemblerRunner);
        assemblerRunner.run();
        CloneSet cloneSet = assemblerRunner.getCloneSet(align.parameters.alignerParameters);
        return new AssembleResult(cloneSet, report, assembler);
    } finally {
        if (close)
            assembler.close();
    }
}

Example 8 with CloneSet

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

the class VersionInfoAction method go.

@Override
public void go(ActionHelper helper) throws Exception {
    String inputFile = parameters.getInputFile();
    String i = inputFile.toLowerCase();
    if (i.endsWith(".vdjca.gz") || i.endsWith(".vdjca")) {
        try (VDJCAlignmentsReader reader = new VDJCAlignmentsReader(inputFile)) {
            reader.init();
            System.out.println("MagicBytes = " + reader.getMagic());
            System.out.println(reader.getVersionInfo());
        }
    } else if (i.endsWith(".clns.gz") || i.endsWith(".clns")) {
        CloneSet cs = CloneSetIO.read(inputFile);
        System.out.println(cs.getVersionInfo());
    } else if (i.endsWith(".clna")) {
        try (ClnAReader reader = new ClnAReader(inputFile, VDJCLibraryRegistry.createDefaultRegistry())) {
            System.out.println(reader.getVersionInfo());
        }
    } else
        throw new ParameterException("Wrong file type.");
}

Example 9 with CloneSet

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

the class ActionAssemble method go.

@Override
public void go(ActionHelper helper) throws Exception {
    // Saving initial timestamp
    long beginTimestamp = System.currentTimeMillis();
    // Checking consistency between actionParameters.doWriteClnA() value and file extension
    if ((actionParameters.getOutputFileName().toLowerCase().endsWith(".clna") && !actionParameters.doWriteClnA()) || (actionParameters.getOutputFileName().toLowerCase().endsWith(".clns") && actionParameters.doWriteClnA()))
        System.out.println("WARNING: Unexpected file extension, use .clns extension for clones-only (normal) output and\n" + ".clna if -a / --write-alignments options specified.");
    // Extracting V/D/J/C gene list from input vdjca file
    final List<VDJCGene> genes;
    final VDJCAlignerParameters alignerParameters;
    try (VDJCAlignmentsReader reader = new VDJCAlignmentsReader(actionParameters.getInputFileName(), VDJCLibraryRegistry.getDefault())) {
        genes = reader.getUsedGenes();
        // Saving aligner parameters to correct assembler parameters
        alignerParameters = reader.getParameters();
    }
    AlignmentsProvider alignmentsProvider = AlignmentsProvider.Util.createProvider(actionParameters.getInputFileName(), VDJCLibraryRegistry.getDefault());
    CloneAssemblerParameters assemblerParameters = actionParameters.getCloneAssemblerParameters();
    // set aligner parameters
    assemblerParameters.updateFrom(alignerParameters);
    // Overriding JSON parameters
    if (!actionParameters.overrides.isEmpty()) {
        assemblerParameters = JsonOverrider.override(assemblerParameters, CloneAssemblerParameters.class, actionParameters.overrides);
        if (assemblerParameters == null) {
            System.err.println("Failed to override some parameter.");
            System.exit(1);
        }
    }
    // Performing assembly
    try (CloneAssembler assembler = new CloneAssembler(assemblerParameters, actionParameters.doWriteClnA() || actionParameters.events != null, genes, alignerParameters.getFeaturesToAlignMap())) {
        // Creating event listener to collect run statistics
        CloneAssemblerReport report = new CloneAssemblerReport();
        report.setStartMillis(beginTimestamp);
        report.setInputFiles(new String[] { actionParameters.getInputFileName() });
        report.setOutputFiles(new String[] { actionParameters.getOutputFileName() });
        report.setCommandLine(helper.getCommandLineArguments());
        assembler.setListener(report);
        // Running assembler
        CloneAssemblerRunner assemblerRunner = new CloneAssemblerRunner(alignmentsProvider, assembler, actionParameters.threads);
        SmartProgressReporter.startProgressReport(assemblerRunner);
        assemblerRunner.run();
        // Getting results
        final CloneSet cloneSet = assemblerRunner.getCloneSet(alignerParameters);
        // Passing final cloneset to assemble last pieces of statistics for report
        report.onClonesetFinished(cloneSet);
        // Writing results
        if (actionParameters.doWriteClnA())
            try (ClnAWriter writer = new ClnAWriter(actionParameters.getOutputFileName())) {
                // writer will supply current stage and completion percent to the progress reporter
                SmartProgressReporter.startProgressReport(writer);
                // Writing clone block
                writer.writeClones(cloneSet);
                // Pre-soring alignments
                try (AlignmentsMappingMerger merged = new AlignmentsMappingMerger(alignmentsProvider.create(), assembler.getAssembledReadsPort())) {
                    writer.sortAlignments(merged, assembler.getAlignmentsCount());
                }
                writer.writeAlignmentsAndIndex();
            }
        else
            try (CloneSetIO.CloneSetWriter writer = new CloneSetIO.CloneSetWriter(cloneSet, actionParameters.getOutputFileName())) {
                SmartProgressReporter.startProgressReport(writer);
                writer.write();
            }
        // Writing report
        report.setFinishMillis(System.currentTimeMillis());
        assert cloneSet.getClones().size() == report.getCloneCount();
        report.setTotalReads(alignmentsProvider.getTotalNumberOfReads());
        // Writing report to stout
        System.out.println("============= Report ==============");
        Util.writeReportToStdout(report);
        if (actionParameters.report != null)
            Util.writeReport(actionParameters.report, report);
        if (actionParameters.jsonReport != null)
            Util.writeJsonReport(actionParameters.jsonReport, report);
        // Writing raw events (not documented feature)
        if (actionParameters.events != null)
            try (PipeWriter<ReadToCloneMapping> writer = new PipeWriter<>(actionParameters.events)) {
                CUtils.drain(assembler.getAssembledReadsPort(), writer);
            }
    }
}

Example 10 with CloneSet

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

the class ActionClonesDiff method go.

@Override
public void go(ActionHelper helper) throws Exception {
    try (InputStream is1 = IOUtil.createIS(params.get1());
        InputStream is2 = IOUtil.createIS(params.get2());
        PrintStream report = params.report().equals(".") ? System.out : new PrintStream(new FileOutputStream(params.report()))) {
        CloneSet cs1 = CloneSetIO.readClns(is1);
        CloneSet cs2 = CloneSetIO.readClns(is2);
        HashMap<CKey, CRec> recs = new HashMap<>();
        populate(recs, cs1, 0);
        populate(recs, cs2, 1);
        int newClones1 = 0, newClones2 = 0;
        long newClones1Reads = 0, newClones2Reads = 0;
        for (CRec cRec : recs.values()) {
            if (cRec.clones[0] == null) {
                newClones2++;
                newClones2Reads += cRec.clones[1].getCount();
            }
            if (cRec.clones[1] == null) {
                newClones1++;
                newClones1Reads += cRec.clones[0].getCount();
            }
        }
        report.println("Unique clones in cloneset 1: " + newClones1 + " (" + Util.PERCENT_FORMAT.format(100.0 * newClones1 / cs1.size()) + "%)");
        report.println("Reads in unique clones in cloneset 1: " + newClones1Reads + " (" + Util.PERCENT_FORMAT.format(100.0 * newClones1Reads / cs1.getTotalCount()) + "%)");
        report.println("Unique clones in cloneset 2: " + newClones2 + " (" + Util.PERCENT_FORMAT.format(100.0 * newClones2 / cs2.size()) + "%)");
        report.println("Reads in unique clones in cloneset 2: " + newClones2Reads + " (" + Util.PERCENT_FORMAT.format(100.0 * newClones2Reads / cs2.getTotalCount()) + "%)");
    }
}