Examples with GeneType io.repseq.core.GeneType used on opensource projects

Example 16 with GeneType

use of io.repseq.core.GeneType in project mixcr by milaboratory.

the class VDJCAlignmentsFormatter method getTargetAsMultiAlignment.

public static MultiAlignmentHelper getTargetAsMultiAlignment(VDJCObject vdjcObject, int targetId, boolean addHitScore, boolean addReads) {
    if (addReads && !(vdjcObject instanceof VDJCAlignments))
        throw new IllegalArgumentException("Read alignments supported only for VDJCAlignments.");
    NSequenceWithQuality target = vdjcObject.getTarget(targetId);
    NucleotideSequence targetSeq = target.getSequence();
    SequencePartitioning partitioning = vdjcObject.getPartitionedTarget(targetId).getPartitioning();
    List<Alignment<NucleotideSequence>> alignments = new ArrayList<>();
    List<String> alignmentLeftComments = new ArrayList<>();
    List<String> alignmentRightComments = new ArrayList<>();
    for (GeneType gt : GeneType.values()) for (VDJCHit hit : vdjcObject.getHits(gt)) {
        Alignment<NucleotideSequence> alignment = hit.getAlignment(targetId);
        if (alignment == null)
            continue;
        alignment = alignment.invert(targetSeq);
        alignments.add(alignment);
        alignmentLeftComments.add(hit.getGene().getName());
        alignmentRightComments.add(" " + (int) (hit.getAlignment(targetId).getScore()) + (addHitScore ? " (" + (int) (hit.getScore()) + ")" : ""));
    }
    // Adding read information
    if (addReads) {
        VDJCAlignments vdjcAlignments = (VDJCAlignments) vdjcObject;
        SequenceHistory history = vdjcAlignments.getHistory(targetId);
        List<SequenceHistory.RawSequence> reads = history.rawReads();
        for (SequenceHistory.RawSequence read : reads) {
            NucleotideSequence seq = vdjcAlignments.getOriginalSequence(read.index).getSequence();
            int offset = history.offset(read.index);
            Alignment<NucleotideSequence> alignment = Aligner.alignOnlySubstitutions(targetSeq, seq, offset, seq.size(), 0, seq.size(), AffineGapAlignmentScoring.IGBLAST_NUCLEOTIDE_SCORING);
            alignments.add(alignment);
            alignmentLeftComments.add(read.index.toString());
            alignmentRightComments.add("");
        }
    }
    MultiAlignmentHelper helper = MultiAlignmentHelper.build(MultiAlignmentHelper.DEFAULT_SETTINGS, new Range(0, target.size()), targetSeq, alignments.toArray(new Alignment[alignments.size()]));
    if (!alignments.isEmpty())
        drawPoints(helper, partitioning, POINTS_FOR_REARRANGED);
    drawAASequence(helper, partitioning, targetSeq);
    helper.addSubjectQuality("Quality", target.getQuality());
    helper.setSubjectLeftTitle("Target" + targetId);
    helper.setSubjectRightTitle(" Score" + (addHitScore ? " (hit score)" : ""));
    for (int i = 0; i < alignmentLeftComments.size(); i++) {
        helper.setQueryLeftTitle(i, alignmentLeftComments.get(i));
        helper.setQueryRightTitle(i, alignmentRightComments.get(i));
    }
    return helper;
}

Example 17 with GeneType

use of io.repseq.core.GeneType in project mixcr by milaboratory.

the class VDJCAlignmentsWriter method header.

@Override
public void header(VDJCAlignerParameters parameters, List<VDJCGene> genes) {
    if (parameters == null || genes == null)
        throw new IllegalArgumentException();
    if (header)
        throw new IllegalStateException();
    // Writing magic bytes
    assert MAGIC_BYTES.length == MAGIC_LENGTH;
    output.write(MAGIC_BYTES);
    // Writing version information
    output.writeUTF(MiXCRVersionInfo.get().getVersionString(MiXCRVersionInfo.OutputType.ToFile));
    // Writing parameters
    output.writeObject(parameters);
    IOUtil.writeAndRegisterGeneReferences(output, genes, parameters);
    // Registering links to features to align
    for (GeneType gt : GeneType.VDJC_REFERENCE) {
        GeneFeature feature = parameters.getFeatureToAlign(gt);
        output.writeObject(feature);
        if (feature != null)
            output.putKnownReference(feature);
    }
    header = true;
}

Example 18 with GeneType

use of io.repseq.core.GeneType in project mixcr by milaboratory.

the class VDJCAlignmentsReader method init.

void init(Map<GeneFeature, GeneFeature> geneFeatureRefs) {
    if (usedGenes != null)
        return;
    assert MAGIC_BYTES.length == MAGIC_LENGTH;
    byte[] magic = new byte[MAGIC_LENGTH];
    input.readFully(magic);
    String magicString = new String(magic);
    this.magic = magicString;
    SerializersManager serializersManager = input.getSerializersManager();
    switch(magicString) {
        case MAGIC_V9:
            serializersManager.registerCustomSerializer(VDJCAlignments.class, new IO.VDJCAlignmentsSerializer21());
        case MAGIC:
            break;
        default:
            throw new RuntimeException("Unsupported file format; .vdjca file of version " + new String(magic) + " while you are running MiXCR " + MAGIC);
    }
    versionInfo = input.readUTF();
    parameters = input.readObject(VDJCAlignerParameters.class);
    this.usedGenes = IOUtil.readAndRegisterGeneReferences(input, vdjcRegistry, parameters);
    // Registering links to features to align
    for (GeneType gt : GeneType.VDJC_REFERENCE) {
        GeneFeature featureParams = parameters.getFeatureToAlign(gt);
        GeneFeature featureDeserialized = input.readObject(GeneFeature.class);
        if (!Objects.equals(featureDeserialized, featureParams))
            throw new RuntimeException("Wrong format.");
        // Find corresponding reference
        if (geneFeatureRefs != null) {
            featureParams = geneFeatureRefs.get(featureParams);
            if (featureParams == null)
                throw new RuntimeException("Absent record for " + featureDeserialized + " in geneFeatureRefs map.");
        }
        if (featureDeserialized != null)
            input.putKnownReference(featureParams);
    }
}

Example 19 with GeneType

use of io.repseq.core.GeneType 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()]);
    }
}

Example 20 with GeneType

use of io.repseq.core.GeneType in project mixcr by milaboratory.

the class IO method readGF2GTMap.

public static EnumMap<GeneType, GeneFeature> readGF2GTMap(PrimitivI input) {
    int count = input.readInt();
    EnumMap<GeneType, GeneFeature> map = new EnumMap<GeneType, GeneFeature>(GeneType.class);
    for (int i = 0; i < count; i++) {
        GeneType gt = input.readObject(GeneType.class);
        GeneFeature gf = input.readObject(GeneFeature.class);
        map.put(gt, gf);
    }
    return map;
}