Search in sources :

Example 1 with VDJCAlignerParameters

use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters in project mixcr by milaboratory.

the class ActionAlign method go.

@Override
@SuppressWarnings("unchecked")
public void go(ActionHelper helper) throws Exception {
    // FIXME remove in 2.2
    if (actionParameters.printNonFunctionalWarnings())
        System.out.println("WARNING: -wf / --non-functional-warnings option is deprecated, will be removed in 2.2 " + "release. Use --verbose instead.");
    // Saving initial timestamp
    long beginTimestamp = System.currentTimeMillis();
    // Getting aligner parameters
    VDJCAlignerParameters alignerParameters = actionParameters.getAlignerParameters();
    // FIXME remove in 2.3
    if (actionParameters.getSaveOriginalReads()) {
        System.out.println("WARNING: -g / --save-reads option is deprecated, will be removed in 2.3 " + "release. Use -OsaveOriginalReads=true.");
        alignerParameters.setSaveOriginalReads(true);
    }
    // FIXME remove in 2.3
    if (actionParameters.getSaveReadDescription()) {
        System.out.println("WARNING: -a / --save-description option is deprecated, will be removed in 2.3 " + "release. Use -OsaveOriginalReads=true.");
        alignerParameters.setSaveOriginalReads(true);
    }
    if (!actionParameters.overrides.isEmpty()) {
        // Perform parameters overriding
        alignerParameters = JsonOverrider.override(alignerParameters, VDJCAlignerParameters.class, actionParameters.overrides);
        if (alignerParameters == null)
            throw new ProcessException("Failed to override some parameter.");
    }
    // FIXME remove in 2.2
    if (actionParameters.allowDifferentVJLoci != null && actionParameters.allowDifferentVJLoci) {
        System.out.println("Warning: usage of --diff-loci is deprecated. Use -OallowChimeras=true instead.");
        alignerParameters.setAllowChimeras(true);
    }
    // Creating aligner
    VDJCAligner aligner = VDJCAligner.createAligner(alignerParameters, actionParameters.isInputPaired(), !actionParameters.getNoMerge());
    // Detect if automatic featureToAlign correction is required
    int totalV = 0, totalVErrors = 0, hasVRegion = 0;
    GeneFeature correctingFeature = alignerParameters.getVAlignerParameters().getGeneFeatureToAlign().hasReversedRegions() ? GeneFeature.VRegionWithP : GeneFeature.VRegion;
    VDJCLibrary library = VDJCLibraryRegistry.getDefault().getLibrary(actionParameters.library, actionParameters.species);
    System.out.println("Reference library: " + library.getLibraryId());
    // Printing library level warnings, if specified for the library
    if (!library.getWarnings().isEmpty()) {
        System.out.println("Library warnings:");
        for (String l : library.getWarnings()) System.out.println(l);
    }
    // Printing citation notice, if specified for the library
    if (!library.getCitations().isEmpty()) {
        System.out.println("Please cite:");
        for (String l : library.getCitations()) System.out.println(l);
    }
    for (VDJCGene gene : library.getGenes(actionParameters.getChains())) {
        if (gene.getGeneType() == GeneType.Variable) {
            totalV++;
            if (!alignerParameters.containsRequiredFeature(gene)) {
                totalVErrors++;
                if (gene.getPartitioning().isAvailable(correctingFeature))
                    hasVRegion++;
            }
        }
    }
    // Performing V featureToAlign correction if needed
    if (totalVErrors > totalV * 0.9 && hasVRegion > totalVErrors * 0.8) {
        System.out.println("WARNING: forcing -OvParameters.geneFeatureToAlign=" + GeneFeature.encode(correctingFeature) + " since current gene feature (" + GeneFeature.encode(alignerParameters.getVAlignerParameters().getGeneFeatureToAlign()) + ") is absent in " + Util.PERCENT_FORMAT.format(100.0 * totalVErrors / totalV) + "% of V genes.");
        alignerParameters.getVAlignerParameters().setGeneFeatureToAlign(correctingFeature);
    }
    int numberOfExcludedNFGenes = 0;
    int numberOfExcludedFGenes = 0;
    for (VDJCGene gene : library.getGenes(actionParameters.getChains())) {
        NucleotideSequence featureSequence = alignerParameters.extractFeatureToAlign(gene);
        // exclusionReason is null ==> gene is not excluded
        String exclusionReason = null;
        if (featureSequence == null)
            exclusionReason = "absent " + GeneFeature.encode(alignerParameters.getFeatureToAlign(gene.getGeneType()));
        else if (featureSequence.containsWildcards())
            exclusionReason = "wildcard symbols in " + GeneFeature.encode(alignerParameters.getFeatureToAlign(gene.getGeneType()));
        if (exclusionReason == null)
            // If there are no reasons to exclude the gene, adding it to aligner
            aligner.addGene(gene);
        else {
            if (gene.isFunctional()) {
                ++numberOfExcludedFGenes;
                if (actionParameters.verbose())
                    System.out.println("WARNING: Functional gene " + gene.getName() + " excluded due to " + exclusionReason);
            } else
                ++numberOfExcludedNFGenes;
        }
    }
    if (actionParameters.printWarnings() && numberOfExcludedFGenes > 0)
        System.out.println("WARNING: " + numberOfExcludedFGenes + " functional genes were excluded, re-run " + "with --verbose option to see the list of excluded genes and exclusion reason.");
    if (actionParameters.verbose() && numberOfExcludedNFGenes > 0)
        System.out.println("WARNING: " + numberOfExcludedNFGenes + " non-functional genes excluded.");
    if (aligner.getVGenesToAlign().isEmpty())
        throw new ProcessException("No V genes to align. Aborting execution. See warnings for more info " + "(turn on verbose warnings by adding --verbose option).");
    if (aligner.getJGenesToAlign().isEmpty())
        throw new ProcessException("No J genes to align. Aborting execution. See warnings for more info " + "(turn on verbose warnings by adding --verbose option).");
    AlignerReport report = new AlignerReport();
    report.setStartMillis(beginTimestamp);
    report.setInputFiles(actionParameters.getInputsForReport());
    report.setOutputFiles(actionParameters.getOutputsForReport());
    report.setCommandLine(helper.getCommandLineArguments());
    // Attaching report to aligner
    aligner.setEventsListener(report);
    try (SequenceReaderCloseable<? extends SequenceRead> reader = actionParameters.createReader();
        VDJCAlignmentsWriter writer = actionParameters.getOutputName().equals(".") ? null : new VDJCAlignmentsWriter(actionParameters.getOutputName());
        SequenceWriter notAlignedWriter = actionParameters.failedReadsR1 == null ? null : (actionParameters.isInputPaired() ? new PairedFastqWriter(actionParameters.failedReadsR1, actionParameters.failedReadsR2) : new SingleFastqWriter(actionParameters.failedReadsR1))) {
        if (writer != null)
            writer.header(aligner);
        OutputPort<? extends SequenceRead> sReads = reader;
        CanReportProgress progress = (CanReportProgress) reader;
        if (actionParameters.limit != 0) {
            sReads = new CountLimitingOutputPort<>(sReads, actionParameters.limit);
            progress = SmartProgressReporter.extractProgress((CountLimitingOutputPort<?>) sReads);
        }
        final boolean writeAllResults = actionParameters.getWriteAllResults();
        EnumMap<GeneType, VDJCHit[]> emptyHits = new EnumMap<>(GeneType.class);
        for (GeneType gt : GeneType.values()) if (alignerParameters.getGeneAlignerParameters(gt) != null)
            emptyHits.put(gt, new VDJCHit[0]);
        final PairedEndReadsLayout readsLayout = alignerParameters.getReadsLayout();
        SmartProgressReporter.startProgressReport("Alignment", progress);
        OutputPort<Chunk<? extends SequenceRead>> mainInputReads = CUtils.buffered((OutputPort) chunked(sReads, 64), 16);
        OutputPort<VDJCAlignmentResult> alignments = unchunked(new ParallelProcessor(mainInputReads, chunked(aligner), actionParameters.threads));
        for (VDJCAlignmentResult result : CUtils.it(new OrderedOutputPort<>(alignments, new Indexer<VDJCAlignmentResult>() {

            @Override
            public long getIndex(VDJCAlignmentResult o) {
                return o.read.getId();
            }
        }))) {
            VDJCAlignments alignment = result.alignment;
            SequenceRead read = result.read;
            if (alignment == null) {
                if (writeAllResults) // Creating empty alignment object if alignment for current read failed
                {
                    Target target = readsLayout.createTargets(read)[0];
                    alignment = new VDJCAlignments(emptyHits, target.targets, SequenceHistory.RawSequence.of(read.getId(), target), alignerParameters.isSaveOriginalReads() ? new SequenceRead[] { read } : null);
                } else {
                    if (notAlignedWriter != null)
                        notAlignedWriter.write(result.read);
                    continue;
                }
            }
            if (alignment.isChimera())
                report.onChimera();
            if (writer != null)
                writer.write(alignment);
        }
        if (writer != null)
            writer.setNumberOfProcessedReads(reader.getNumberOfReads());
    }
    report.setFinishMillis(System.currentTimeMillis());
    // 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);
}
Also used : VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) VDJCAlignmentsWriter(com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriter) VDJCAligner(com.milaboratory.mixcr.vdjaligners.VDJCAligner) ParallelProcessor(cc.redberry.pipe.blocks.ParallelProcessor) Target(com.milaboratory.core.Target) Indexer(cc.redberry.pipe.util.Indexer) CanReportProgress(com.milaboratory.util.CanReportProgress) SingleFastqWriter(com.milaboratory.core.io.sequence.fastq.SingleFastqWriter) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments) PairedEndReadsLayout(com.milaboratory.core.PairedEndReadsLayout) VDJCAlignmentResult(com.milaboratory.mixcr.vdjaligners.VDJCAlignmentResult) PairedFastqWriter(com.milaboratory.core.io.sequence.fastq.PairedFastqWriter) Chunk(cc.redberry.pipe.util.Chunk) ProcessException(com.milaboratory.cli.ProcessException) CountLimitingOutputPort(cc.redberry.pipe.util.CountLimitingOutputPort) SequenceWriter(com.milaboratory.core.io.sequence.SequenceWriter) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) SequenceRead(com.milaboratory.core.io.sequence.SequenceRead)

Example 2 with VDJCAlignerParameters

use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters in project mixcr by milaboratory.

the class ActionAssembleContigs method go.

@Override
public void go(ActionHelper helper) throws Exception {
    // TODO FIX!!!!!!!!!!!!!
    if (parameters.threads > 1)
        throw new ParameterException("Multithreaded processing is not supported yet.");
    long beginTimestamp = System.currentTimeMillis();
    FullSeqAssemblerParameters p = FullSeqAssemblerParameters.getByName("default");
    if (!parameters.overrides.isEmpty()) {
        // Perform parameters overriding
        p = JsonOverrider.override(p, FullSeqAssemblerParameters.class, parameters.overrides);
        if (p == null)
            throw new ProcessException("Failed to override some parameter.");
    }
    final FullSeqAssemblerReport report = new FullSeqAssemblerReport();
    FullSeqAssemblerParameters assemblerParameters = p;
    int totalClonesCount = 0;
    List<VDJCGene> genes;
    VDJCAlignerParameters alignerParameters;
    CloneAssemblerParameters cloneAssemblerParameters;
    try (ClnAReader reader = new ClnAReader(parameters.getInputFileName(), VDJCLibraryRegistry.getDefault());
        PrimitivO tmpOut = new PrimitivO(new BufferedOutputStream(new FileOutputStream(parameters.getOutputFileName())))) {
        final CloneFactory cloneFactory = new CloneFactory(reader.getAssemblerParameters().getCloneFactoryParameters(), reader.getAssemblingFeatures(), reader.getGenes(), reader.getAlignerParameters().getFeaturesToAlignMap());
        alignerParameters = reader.getAlignerParameters();
        cloneAssemblerParameters = reader.getAssemblerParameters();
        genes = reader.getGenes();
        IOUtil.registerGeneReferences(tmpOut, genes, alignerParameters);
        ClnAReader.CloneAlignmentsPort cloneAlignmentsPort = reader.clonesAndAlignments();
        SmartProgressReporter.startProgressReport("Assembling", cloneAlignmentsPort);
        OutputPort<Clone[]> parallelProcessor = new ParallelProcessor<>(cloneAlignmentsPort, cloneAlignments -> {
            FullSeqAssembler fullSeqAssembler = new FullSeqAssembler(cloneFactory, assemblerParameters, cloneAlignments.clone, alignerParameters);
            fullSeqAssembler.setReport(report);
            FullSeqAssembler.RawVariantsData rawVariantsData = fullSeqAssembler.calculateRawData(() -> {
                try {
                    return cloneAlignments.alignments();
                } catch (IOException e) {
                    throw new RuntimeException(e);
                }
            });
            return fullSeqAssembler.callVariants(rawVariantsData);
        }, parameters.threads);
        for (Clone[] clones : CUtils.it(parallelProcessor)) {
            totalClonesCount += clones.length;
            for (Clone cl : clones) tmpOut.writeObject(cl);
        }
        assert report.getInitialCloneCount() == reader.numberOfClones();
    }
    assert report.getFinalCloneCount() == totalClonesCount;
    assert report.getFinalCloneCount() >= report.getInitialCloneCount();
    Clone[] clones = new Clone[totalClonesCount];
    try (PrimitivI tmpIn = new PrimitivI(new BufferedInputStream(new FileInputStream(parameters.getOutputFileName())))) {
        IOUtil.registerGeneReferences(tmpIn, genes, alignerParameters);
        int i = 0;
        for (Clone clone : CUtils.it(new PipeDataInputReader<>(Clone.class, tmpIn, totalClonesCount))) clones[i++] = clone;
    }
    Arrays.sort(clones, Comparator.comparingDouble(c -> -c.getCount()));
    for (int i = 0; i < clones.length; i++) clones[i] = clones[i].setId(i);
    CloneSet cloneSet = new CloneSet(Arrays.asList(clones), genes, alignerParameters.getFeaturesToAlignMap(), alignerParameters, cloneAssemblerParameters);
    try (CloneSetIO.CloneSetWriter writer = new CloneSetIO.CloneSetWriter(cloneSet, parameters.getOutputFileName())) {
        SmartProgressReporter.startProgressReport(writer);
        writer.write();
    }
    ReportWrapper reportWrapper = new ReportWrapper(command(), report);
    reportWrapper.setStartMillis(beginTimestamp);
    reportWrapper.setInputFiles(parameters.getInputFileName());
    reportWrapper.setOutputFiles(parameters.getOutputFileName());
    reportWrapper.setCommandLine(helper.getCommandLineArguments());
    reportWrapper.setFinishMillis(System.currentTimeMillis());
    // Writing report to stout
    System.out.println("============= Report ==============");
    Util.writeReportToStdout(report);
    if (parameters.report != null)
        Util.writeReport(parameters.report, reportWrapper);
    if (parameters.jsonReport != null)
        Util.writeJsonReport(parameters.jsonReport, reportWrapper);
}
Also used : Parameters(com.beust.jcommander.Parameters) java.util(java.util) ParameterException(com.beust.jcommander.ParameterException) Action(com.milaboratory.cli.Action) Parameter(com.beust.jcommander.Parameter) CloneAssemblerParameters(com.milaboratory.mixcr.assembler.CloneAssemblerParameters) com.milaboratory.mixcr.basictypes(com.milaboratory.mixcr.basictypes) FullSeqAssemblerParameters(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssemblerParameters) ParallelProcessor(cc.redberry.pipe.blocks.ParallelProcessor) VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) PositiveInteger(com.beust.jcommander.validators.PositiveInteger) CloneFactory(com.milaboratory.mixcr.assembler.CloneFactory) FullSeqAssembler(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssembler) ActionHelper(com.milaboratory.cli.ActionHelper) PrimitivO(com.milaboratory.primitivio.PrimitivO) CUtils(cc.redberry.pipe.CUtils) FullSeqAssemblerReport(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssemblerReport) PrimitivI(com.milaboratory.primitivio.PrimitivI) OutputPort(cc.redberry.pipe.OutputPort) PipeDataInputReader(com.milaboratory.primitivio.PipeDataInputReader) java.io(java.io) VDJCGene(io.repseq.core.VDJCGene) SmartProgressReporter(com.milaboratory.util.SmartProgressReporter) DynamicParameter(com.beust.jcommander.DynamicParameter) VDJCLibraryRegistry(io.repseq.core.VDJCLibraryRegistry) ProcessException(com.milaboratory.cli.ProcessException) ActionParametersWithOutput(com.milaboratory.cli.ActionParametersWithOutput) VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) FullSeqAssemblerParameters(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssemblerParameters) CloneAssemblerParameters(com.milaboratory.mixcr.assembler.CloneAssemblerParameters) ParallelProcessor(cc.redberry.pipe.blocks.ParallelProcessor) ParameterException(com.beust.jcommander.ParameterException) PrimitivO(com.milaboratory.primitivio.PrimitivO) FullSeqAssembler(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssembler) FullSeqAssemblerReport(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssemblerReport) ProcessException(com.milaboratory.cli.ProcessException) VDJCGene(io.repseq.core.VDJCGene) CloneFactory(com.milaboratory.mixcr.assembler.CloneFactory) PrimitivI(com.milaboratory.primitivio.PrimitivI)

Example 3 with VDJCAlignerParameters

use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters in project mixcr by milaboratory.

the class CloneSetIO method readClns.

public static CloneSet readClns(InputStream inputStream, VDJCLibraryRegistry libraryRegistry) {
    PrimitivI input = new PrimitivI(inputStream);
    // Registering custom serializer
    input.getSerializersManager().registerCustomSerializer(GeneFeature.class, new GeneFeatureSerializer(true));
    byte[] magicBytes = new byte[MAGIC_LENGTH];
    input.readFully(magicBytes);
    String magicString = new String(magicBytes);
    SerializersManager serializersManager = input.getSerializersManager();
    switch(magicString) {
        case MAGIC_V5:
            serializersManager.registerCustomSerializer(Clone.class, new CompatibilityIO.CloneSerializerV5());
            break;
        case MAGIC:
            break;
        default:
            throw new RuntimeException("Unsupported file format; .clns file of version " + magicString + " while you are running MiXCR " + MAGIC);
    }
    String versionInfo = input.readUTF();
    if (!magicString.equals(MAGIC))
        // Dropping this field for v5 files
        input.readObject(GeneFeature[].class);
    VDJCAlignerParameters alignerParameters;
    CloneAssemblerParameters assemblerParameters;
    if (magicString.equals(MAGIC)) {
        alignerParameters = input.readObject(VDJCAlignerParameters.class);
        assemblerParameters = input.readObject(CloneAssemblerParameters.class);
    } else {
        alignerParameters = null;
        assemblerParameters = null;
    }
    EnumMap<GeneType, GeneFeature> alignedFeatures = IO.readGF2GTMap(input);
    List<VDJCGene> genes = IOUtil.readAndRegisterGeneReferences(input, libraryRegistry, new GT2GFAdapter(alignedFeatures));
    int count = input.readInt();
    List<Clone> clones = new ArrayList<>(count);
    for (int i = 0; i < count; i++) clones.add(input.readObject(Clone.class));
    CloneSet cloneSet = new CloneSet(clones, genes, alignedFeatures, alignerParameters, assemblerParameters);
    cloneSet.versionInfo = versionInfo;
    return cloneSet;
}
Also used : VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) SerializersManager(com.milaboratory.primitivio.SerializersManager) ArrayList(java.util.ArrayList) CloneAssemblerParameters(com.milaboratory.mixcr.assembler.CloneAssemblerParameters) PrimitivI(com.milaboratory.primitivio.PrimitivI)

Example 4 with VDJCAlignerParameters

use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters in project mixcr by milaboratory.

the class VDJCObjectExtenderTest method testTripleRead.

@Test
public void testTripleRead() throws Exception {
    final boolean print = true;
    Well44497b rg = new Well44497b(12312);
    final VDJCAlignerParameters rnaSeqParams = VDJCParametersPresets.getByName("rna-seq");
    final PartialAlignmentsAssemblerAligner aligner = new PartialAlignmentsAssemblerAligner(rnaSeqParams);
    final VDJCLibrary lib = VDJCLibraryRegistry.getDefault().getLibrary("default", "hs");
    for (VDJCGene gene : VDJCLibraryRegistry.getDefault().getLibrary("default", "hs").getGenes()) if (gene.isFunctional())
        aligner.addGene(gene);
    TargetBuilder.VDJCGenes genes = new TargetBuilder.VDJCGenes(lib, "TRBV12-3*00", "TRBD1*00", "TRBJ1-3*00", "TRBC2*00");
    // | 310  | 338   | 438
    // 250V + 60CDR3 (20V 7N 10D 3N 20J) + 28J + 100C + 100N
    final NucleotideSequence baseSeq = TargetBuilder.generateSequence(genes, "{CDR3Begin(-250)}V*270 NNNNNNN {DBegin(0)}D*10 NNN {CDR3End(-20):FR4End} {CBegin}C*100 N*100", rg);
    F4 goAssert = new F4() {

        @Override
        public VDJCAlignments go(int len, int offset1, int offset2, int offset3) {
            NucleotideSequence seq1 = baseSeq.getRange(offset1, offset1 + len);
            NucleotideSequence seq2 = baseSeq.getRange(offset2, offset2 + len);
            NucleotideSequence seq3 = offset3 == -1 ? null : baseSeq.getRange(offset3, offset3 + len);
            VDJCAlignmentResult<VDJCMultiRead> alignment = offset3 == -1 ? aligner.process(MiXCRTestUtils.createMultiRead(seq1, seq2)) : aligner.process(MiXCRTestUtils.createMultiRead(seq1, seq2, seq3));
            VDJCAlignments al = alignment.alignment;
            Assert.assertNotNull(al);
            VDJCObjectExtender<VDJCAlignments> extender = new VDJCObjectExtender<>(Chains.TCR, (byte) 35, rnaSeqParams.getVAlignerParameters().getParameters().getScoring(), rnaSeqParams.getJAlignerParameters().getParameters().getScoring(), 100, 70, ReferencePoint.CDR3Begin, ReferencePoint.CDR3End);
            MiXCRTestUtils.assertAlignments(al);
            if (print) {
                MiXCRTestUtils.printAlignment(al);
                System.out.println();
                System.out.println("-------------------------------------------");
                System.out.println();
            }
            VDJCAlignments processed = extender.process(al);
            if (print) {
                MiXCRTestUtils.printAlignment(processed);
                System.out.println();
                System.out.println("===========================================");
                System.out.println();
                System.out.println();
            }
            MiXCRTestUtils.assertAlignments(processed);
            Assert.assertEquals(al.getFeature(GeneFeature.VDJunction), processed.getFeature(GeneFeature.VDJunction));
            Assert.assertEquals(al.getFeature(GeneFeature.DJJunction), processed.getFeature(GeneFeature.DJJunction));
            Assert.assertEquals(al.getFeature(GeneFeature.VJJunction), processed.getFeature(GeneFeature.VJJunction));
            return processed;
        }
    };
    VDJCAlignments a1 = goAssert.go(60, 245, 307, 450);
    Assert.assertEquals(2, a1.numberOfTargets());
    VDJCAlignments a2 = goAssert.go(60, 245, 315, 450);
    Assert.assertEquals(3, a2.numberOfTargets());
    VDJCAlignments a3 = goAssert.go(60, 245, 315, -1);
    Assert.assertEquals(2, a3.numberOfTargets());
    VDJCAlignments a4 = goAssert.go(60, 245, 307, -1);
    Assert.assertEquals(1, a4.numberOfTargets());
    VDJCAlignments a5 = goAssert.go(53, 252, 307, -1);
    Assert.assertEquals(1, a5.numberOfTargets());
    VDJCAlignments a6 = goAssert.go(53, 252, 307, 450);
    Assert.assertEquals(2, a6.numberOfTargets());
}
Also used : VDJCMultiRead(com.milaboratory.mixcr.partialassembler.VDJCMultiRead) VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) PartialAlignmentsAssemblerAligner(com.milaboratory.mixcr.partialassembler.PartialAlignmentsAssemblerAligner) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) Well44497b(org.apache.commons.math3.random.Well44497b) TargetBuilder(com.milaboratory.mixcr.tests.TargetBuilder) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments) Test(org.junit.Test)

Example 5 with VDJCAlignerParameters

use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters in project mixcr by milaboratory.

the class PartialAlignmentsAssembler method searchOverlaps.

public void searchOverlaps(VDJCAlignmentsReader reader) {
    final VDJCAlignerParameters alignerParameters = reader.getParameters();
    PartialAlignmentsAssemblerAligner aligner = new PartialAlignmentsAssemblerAligner(alignerParameters);
    targetMerger.setAlignerParameters(alignerParameters);
    for (VDJCGene gene : reader.getUsedGenes()) aligner.addGene(gene);
    for (final VDJCAlignments alignment : CUtils.it(reader)) {
        total.incrementAndGet();
        if (alignment.getFeature(GeneFeature.CDR3) != null) {
            containsCDR3.incrementAndGet();
            if (!overlappedOnly) {
                totalWritten.incrementAndGet();
                writer.write(alignment);
            }
            continue;
        }
        if (alreadyMergedIds.contains(alignment.getAlignmentsIndex()))
            continue;
        final OverlapSearchResult searchResult = searchOverlaps(alignment, alignerParameters.isAllowChimeras());
        // Common procedure to cancel processing of current input alignment if it fails to pass some good
        // overlap filtering criterion
        final Runnable cancelCurrentResult = new Runnable() {

            @Override
            public void run() {
                if (searchResult != null)
                    searchResult.cancel();
                if (writePartial && !overlappedOnly && notInLeftIndexIds.contains(alignment.getAlignmentsIndex())) {
                    totalWritten.incrementAndGet();
                    partialAsIs.incrementAndGet();
                    writer.write(alignment);
                }
            }
        };
        if (searchResult == null) {
            cancelCurrentResult.run();
            continue;
        }
        List<AlignedTarget> mergedTargets = searchResult.result;
        VDJCMultiRead mRead = new VDJCMultiRead(mergedTargets);
        final VDJCAlignments mAlignment = aligner.process(mRead).alignment;
        // Checking number of overlapped non-template (NRegion) letters
        int overlapTargetId = -1;
        Range overlapRange = null;
        for (int i = 0; i < mergedTargets.size(); i++) {
            overlapRange = AlignedTarget.getOverlapRange(mergedTargets.get(i));
            if (overlapRange != null) {
                overlapTargetId = i;
                break;
            }
        }
        if (overlapTargetId == -1) {
            // No alignments for Best V Hit and Best J Hit in central (overlapped) target
            cancelCurrentResult.run();
            continue;
        }
        int targetLength = mergedTargets.get(overlapTargetId).getTarget().size();
        VDJCHit bestVHit = mAlignment.getBestHit(GeneType.Variable), bestJHit = mAlignment.getBestHit(GeneType.Joining);
        if (bestVHit == null || bestJHit == null || bestVHit.getAlignment(overlapTargetId) == null || bestJHit.getAlignment(overlapTargetId) == null) {
            cancelCurrentResult.run();
            continue;
        }
        int ndnRegionBegin = 0;
        int ndnRegionEnd = targetLength;
        Alignment<NucleotideSequence> vAlignment = bestVHit.getAlignment(overlapTargetId);
        if (vAlignment != null)
            ndnRegionBegin = vAlignment.getSequence2Range().getTo();
        Alignment<NucleotideSequence> jAlignment = bestJHit.getAlignment(overlapTargetId);
        if (jAlignment != null)
            ndnRegionEnd = jAlignment.getSequence2Range().getFrom();
        RangeSet nRegion = ndnRegionBegin >= ndnRegionEnd ? RangeSet.EMPTY : RangeSet.create(ndnRegionBegin, ndnRegionEnd);
        Range dRange = mAlignment.getPartitionedTarget(overlapTargetId).getPartitioning().getRange(GeneFeature.DRegionTrimmed);
        if (dRange != null)
            nRegion = nRegion.subtract(dRange);
        RangeSet nRegionInOverlap = nRegion.intersection(overlapRange);
        int actualNRegionLength = nRegion.totalLength();
        int minimalN = Math.min(minimalNOverlap, actualNRegionLength);
        if (nRegionInOverlap.totalLength() < minimalN) {
            droppedSmallOverlapNRegion.incrementAndGet();
            cancelCurrentResult.run();
            continue;
        }
        // DDDDDDDDDJJJJJJJJJJJJJJJ
        if (minimalN == 0 && (!overlapRange.contains(ndnRegionBegin - 1) || !overlapRange.contains(ndnRegionEnd))) {
            droppedNoNRegion.incrementAndGet();
            cancelCurrentResult.run();
            continue;
        }
        overlapped.incrementAndGet();
        totalWritten.incrementAndGet();
        writer.write(mAlignment);
        // Saving alignment that where merge to prevent it's use as left part
        alreadyMergedIds.add(alignment.getAlignmentsIndex());
        alreadyMergedIds.add(searchResult.KMerInfo.alignments.getAlignmentsIndex());
    }
    if (writePartial && !overlappedOnly)
        for (List<KMerInfo> kMerInfos : kToIndexLeft.valueCollection()) for (KMerInfo kMerInfo : kMerInfos) if (alreadyMergedIds.add(kMerInfo.alignments.getAlignmentsIndex())) {
            totalWritten.incrementAndGet();
            partialAsIs.incrementAndGet();
            writer.write(kMerInfo.getAlignments());
        }
    writer.setNumberOfProcessedReads(reader.getNumberOfReads() - overlapped.get());
}
Also used : VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) Range(com.milaboratory.core.Range) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) ArrayList(java.util.ArrayList) List(java.util.List)

Aggregations

VDJCAlignerParameters (com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters)12 NucleotideSequence (com.milaboratory.core.sequence.NucleotideSequence)6 VDJCGene (io.repseq.core.VDJCGene)6 VDJCAlignments (com.milaboratory.mixcr.basictypes.VDJCAlignments)5 Test (org.junit.Test)5 ArrayList (java.util.ArrayList)4 ParallelProcessor (cc.redberry.pipe.blocks.ParallelProcessor)3 TargetBuilder (com.milaboratory.mixcr.tests.TargetBuilder)3 Well44497b (org.apache.commons.math3.random.Well44497b)3 Chunk (cc.redberry.pipe.util.Chunk)2 Indexer (cc.redberry.pipe.util.Indexer)2 ProcessException (com.milaboratory.cli.ProcessException)2 CloneAssemblerParameters (com.milaboratory.mixcr.assembler.CloneAssemblerParameters)2 PartialAlignmentsAssemblerAligner (com.milaboratory.mixcr.partialassembler.PartialAlignmentsAssemblerAligner)2 VDJCMultiRead (com.milaboratory.mixcr.partialassembler.VDJCMultiRead)2 PrimitivI (com.milaboratory.primitivio.PrimitivI)2 CUtils (cc.redberry.pipe.CUtils)1 OutputPort (cc.redberry.pipe.OutputPort)1 CountLimitingOutputPort (cc.redberry.pipe.util.CountLimitingOutputPort)1 DynamicParameter (com.beust.jcommander.DynamicParameter)1