Search in sources :

Example 76 with NucleotideSequence

use of com.milaboratory.core.sequence.NucleotideSequence in project mixcr by milaboratory.

the class PartialAlignmentsAssemblerTest method test2a.

@Test
public void test2a() throws Exception {
    RandomUtil.reseedThreadLocal(47);
    final InputTestData input = createTestData(47);
    final NucleotideSequence reference = input.reference;
    final EnumMap<GeneType, int[]> refPositions = input.refPositions;
    NucleotideSequence left1 = reference.getRange(refPositions.get(Diversity)[0] - 85, refPositions.get(Diversity)[0] + 10);
    final char[] chars = left1.toString().toCharArray();
    chars[3] = 'a';
    left1 = new NucleotideSequence(new String(chars));
    PairedRead[] data = { createPair(0, left1, reference.getRange(refPositions.get(Diversity)[1], refPositions.get(Diversity)[1] + 85).getReverseComplement()), createPair(1, reference.getRange(refPositions.get(Diversity)[0] - 135, refPositions.get(Diversity)[0] - 30), reference.getRange(refPositions.get(Diversity)[0] - 8, refPositions.get(Diversity)[0] + 85).getReverseComplement()) };
    final TestResult testResult = processData(data, input);
    for (VDJCAlignments al : testResult.assembled) {
        MiXCRTestUtils.printAlignment(al);
    // System.out.println(input.VJJunction);
    // System.out.println(al.getFeature(GeneFeature.VJJunction).getSequence());
    // Assert.assertTrue(input.VJJunction.toString().contains(al.getFeature(GeneFeature.VJJunction).getSequence().toString()));
    }
}
Also used : NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) GeneType(io.repseq.core.GeneType) PairedRead(com.milaboratory.core.io.sequence.PairedRead) Test(org.junit.Test)

Example 77 with NucleotideSequence

use of com.milaboratory.core.sequence.NucleotideSequence in project mixcr by milaboratory.

the class PartialAlignmentsAssemblerTest method createTestData.

public static InputTestData createTestData(long seed) throws Exception {
    EnumMap<GeneType, String> geneNames = new EnumMap<GeneType, String>(GeneType.class) {

        {
            put(Variable, "TRBV20-1*00");
            put(Diversity, "TRBD2*00");
            put(Joining, "TRBJ2-6*00");
            put(Constant, "TRBC2*00");
        }
    };
    // config
    RandomGenerator rnd = RandomUtil.getThreadLocalRandom();
    rnd.setSeed(seed);
    final VDJCAlignerParameters defaultFeatures = VDJCParametersPresets.getByName("default");
    defaultFeatures.getVAlignerParameters().setGeneFeatureToAlign(VRegion);
    defaultFeatures.getDAlignerParameters().setGeneFeatureToAlign(DRegion);
    defaultFeatures.getJAlignerParameters().setGeneFeatureToAlign(JRegion);
    // used alleles
    EnumMap<GeneType, VDJCGene> genes = new EnumMap<>(GeneType.class);
    // germline parts of sequences
    EnumMap<GeneType, NucleotideSequence> germlineRegions = gtMap();
    // left, right cut of germline
    EnumMap<GeneType, int[]> germlineCuts = gtMap();
    // begin, end positions in assembled sequence
    EnumMap<GeneType, int[]> refPositions = gtMap();
    // single assembled sequence
    SequenceBuilder<NucleotideSequence> referenceBuilder = NucleotideSequence.ALPHABET.createBuilder();
    NucleotideSequence VDJunction = TestUtil.randomSequence(NucleotideSequence.ALPHABET, 3, 10);
    NucleotideSequence DJJunction = TestUtil.randomSequence(NucleotideSequence.ALPHABET, 3, 10);
    for (GeneType gt : GeneType.VDJC_REFERENCE) {
        VDJCGene gene = VDJCLibraryRegistry.getDefault().getLibrary("default", "hs").get(geneNames.get(gt));
        NucleotideSequence seq = gene.getFeature(defaultFeatures.getFeatureToAlign(gt));
        int[] cuts = null;
        switch(gt) {
            case Variable:
                cuts = new int[] { 0, rnd.nextInt(gene.getFeature(GermlineVCDR3Part).size() - 5) };
                break;
            case Diversity:
                cuts = new int[] { rnd.nextInt(seq.size() / 3), rnd.nextInt(seq.size() / 3) };
                break;
            case Joining:
                cuts = new int[] { rnd.nextInt(gene.getFeature(GermlineJCDR3Part).size() - 5), 0 };
                break;
            case Constant:
                cuts = new int[] { 0, rnd.nextInt(seq.size() / 2) };
                break;
        }
        NucleotideSequence gSeq = seq.getRange(cuts[0], seq.size() - cuts[1]);
        int[] positions = new int[2];
        positions[0] = referenceBuilder.size();
        referenceBuilder.append(gSeq);
        positions[1] = referenceBuilder.size();
        if (gt == Variable)
            referenceBuilder.append(VDJunction);
        if (gt == Diversity)
            referenceBuilder.append(DJJunction);
        genes.put(gt, gene);
        germlineCuts.put(gt, cuts);
        germlineRegions.put(gt, gSeq);
        refPositions.put(gt, positions);
    }
    NucleotideSequence VJJunction = NucleotideSequence.ALPHABET.createBuilder().append(VDJunction).append(germlineRegions.get(Diversity)).append(DJJunction).createAndDestroy();
    NucleotideSequence reference = referenceBuilder.createAndDestroy();
    return new InputTestData(genes, germlineRegions, germlineCuts, refPositions, VDJunction, DJJunction, reference, VJJunction);
}
Also used : VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) RandomGenerator(org.apache.commons.math3.random.RandomGenerator) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) VDJCGene(io.repseq.core.VDJCGene) GeneType(io.repseq.core.GeneType) EnumMap(java.util.EnumMap)

Example 78 with NucleotideSequence

use of com.milaboratory.core.sequence.NucleotideSequence in project mixcr by milaboratory.

the class PartialAlignmentsAssemblerTest method test1.

// @Test
// public void testMaxAllele() throws Exception {
// final LociLibrary ll = LociLibraryManager.getDefault().getLibrary("mi");
// final Locus locus = Locus.TRB;
// 
// for (GeneFeature feature : new GeneFeature[]{VRegionWithP, DRegion, JRegionWithP, CRegion}) {
// Allele maxAllele = null;
// for (Allele allele : ll.getAllAlleles()) {
// if (allele.getLocusContainer().getSpeciesAndLocus().taxonId != Species.HomoSapiens)
// continue;
// if (!allele.isFunctional())
// continue;
// if (allele.getLocus() != locus)
// continue;
// if (maxAllele == null && allele.getFeature(feature) != null)
// maxAllele = allele;
// if (allele.getFeature(feature) != null && allele.getFeature(feature).size() > maxAllele.getFeature(feature).size())
// maxAllele = allele;
// }
// 
// System.out.println(maxAllele.getName() + "    Size: " + maxAllele.getFeature(feature).size());
// }
// }
@Test
public void test1() throws Exception {
    final InputTestData input = createTestData();
    final NucleotideSequence reference = input.reference;
    final EnumMap<GeneType, int[]> refPositions = input.refPositions;
    PairedRead[] data = { createPair(0, reference.getRange(refPositions.get(Variable)[1] - 85, refPositions.get(Variable)[1] + 15), reference.getRange(refPositions.get(Joining)[1] - 20, refPositions.get(Joining)[1] + 80).getReverseComplement()), createPair(1, reference.getRange(refPositions.get(Variable)[1] - 186, refPositions.get(Variable)[1] - 86), reference.getRange(refPositions.get(Variable)[1] - 10, refPositions.get(Variable)[1] + 102).getReverseComplement()) };
    final TestResult testResult = processData(data, input);
    for (VDJCAlignments al : testResult.assembled) {
        MiXCRTestUtils.printAlignment(al);
    }
}
Also used : NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) GeneType(io.repseq.core.GeneType) PairedRead(com.milaboratory.core.io.sequence.PairedRead) Test(org.junit.Test)

Example 79 with NucleotideSequence

use of com.milaboratory.core.sequence.NucleotideSequence in project mixcr by milaboratory.

the class TargetMergerTest method test1.

@Test
public void test1() throws Exception {
    // acgatgggcgcaa     atatagggagctccgatcgacatcgagTTTTTtatcgccctggtacgatcccggtgacaaagcgttc   ggacctgtctggacgctagaacgcag
    // acgatgggcgcaa                          atcgggtatcgccctggtacgatAccggtga aaagcgttc   ggacctgtctggacgctagaacgcag
    // acgatgggcgcaa     atatagg agctAcgatcgacatcgggtatcgccc                              ggacctgtctggacgctagaacgcag
    final NucleotideSequence seq1 = new NucleotideSequence("atatagg agctAcgatcgacatcgAgtatcgccc                           ".replace(" ", ""));
    final NucleotideSequence seq2 = new NucleotideSequence("                     atcgggtatcgccctggtacgatAccggtga aaagcgttc".replace(" ", ""));
    final NucleotideSequence merg = new NucleotideSequence("atatagg agctAcgatcgacatcgggtatcgccctggtacgatAccggtga aaagcgttc".replace(" ", ""));
    final NucleotideSequence reference = new NucleotideSequence("acgatgggcgcaa     atatagggagctccgatcgacatcgagTTTTTtatcgccctggtacgatcccggtgacaaagcgttc   ggacctgtctggacgctagaacgcag".replace(" ", ""));
    AffineGapAlignmentScoring<NucleotideSequence> scoring = new AffineGapAlignmentScoring<>(NucleotideSequence.ALPHABET, 8, -5, -8, -1);
    final Alignment<NucleotideSequence> al1 = Aligner.alignLocal(scoring, reference, seq1);
    final Alignment<NucleotideSequence> al2 = Aligner.alignLocal(scoring, reference, seq2);
    final Alignment<NucleotideSequence> merge = TargetMerger.merge(scoring, 10, merg, 20, al1, al2);
    System.out.println(al1 + "\n\n" + al2 + "\n\n" + merge);
    Assert.assertEquals(merg.getRange(merge.getSequence2Range()), AlignmentUtils.getAlignedSequence2Part(merge));
}
Also used : NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) AffineGapAlignmentScoring(com.milaboratory.core.alignment.AffineGapAlignmentScoring) Test(org.junit.Test)

Example 80 with NucleotideSequence

use of com.milaboratory.core.sequence.NucleotideSequence in project mixcr by milaboratory.

the class MiXCRTestUtils method assertAlignments.

public static void assertAlignments(VDJCAlignments alignments) {
    for (GeneType gt : GeneType.VDJC_REFERENCE) {
        for (VDJCHit hit : alignments.getHits(gt)) {
            for (int targetIndex = 0; targetIndex < alignments.numberOfTargets(); targetIndex++) {
                Alignment<NucleotideSequence> al = hit.getAlignment(targetIndex);
                if (al == null)
                    continue;
                NucleotideSequence sequence = alignments.getTarget(targetIndex).getSequence();
                AlignerTest.assertAlignment(al, sequence);
            }
        }
    }
}
Also used : NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) GeneType(io.repseq.core.GeneType) VDJCHit(com.milaboratory.mixcr.basictypes.VDJCHit)

Aggregations

NucleotideSequence (com.milaboratory.core.sequence.NucleotideSequence)81 Test (org.junit.Test)32 Range (com.milaboratory.core.Range)19 VDJCGene (io.repseq.core.VDJCGene)15 GeneType (io.repseq.core.GeneType)14 NSequenceWithQuality (com.milaboratory.core.sequence.NSequenceWithQuality)13 GeneFeature (io.repseq.core.GeneFeature)9 VDJCAlignments (com.milaboratory.mixcr.basictypes.VDJCAlignments)8 VDJCHit (com.milaboratory.mixcr.basictypes.VDJCHit)8 ReferencePoint (io.repseq.core.ReferencePoint)7 VDJCLibrary (io.repseq.core.VDJCLibrary)7 ArrayList (java.util.ArrayList)7 VDJCAlignerParameters (com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters)6 AlignmentHit (com.milaboratory.core.alignment.batch.AlignmentHit)5 PairedRead (com.milaboratory.core.io.sequence.PairedRead)5 AminoAcidSequence (com.milaboratory.core.sequence.AminoAcidSequence)5 Path (java.nio.file.Path)5 Alignment (com.milaboratory.core.alignment.Alignment)4 Pattern (java.util.regex.Pattern)4 AlignmentHelper (com.milaboratory.core.alignment.AlignmentHelper)3