Search in sources :

Example 51 with SAMReadGroupRecord

use of htsjdk.samtools.SAMReadGroupRecord in project gatk by broadinstitute.

the class ReadPileupUnitTest method testSimplePileupWithIndelsOffset.

@Test
public void testSimplePileupWithIndelsOffset() {
    final int readlength = 10;
    final byte[] bases1 = Utils.repeatChars('A', readlength);
    final byte[] quals1 = Utils.repeatBytes((byte) 10, readlength);
    final String cigar1 = "10M";
    final byte[] bases2 = Utils.repeatChars('C', readlength);
    final byte[] quals2 = Utils.repeatBytes((byte) 20, readlength);
    final String cigar2 = "5M3I2M";
    final SAMFileHeader header = ArtificialReadUtils.createArtificialSamHeader();
    final GATKRead read1 = ArtificialReadUtils.createArtificialRead(header, bases1, quals1, cigar1);
    read1.setName("read1");
    final GATKRead read2 = ArtificialReadUtils.createArtificialRead(header, bases2, quals2, cigar2);
    read1.setName("read2");
    final List<GATKRead> reads = Arrays.asList(read1, read2);
    final int offset = 4;
    final ReadPileup pu = new ReadPileup(loc, reads, offset);
    final PileupElement firstElem = pu.getElementForRead(read1);
    Assert.assertNull(firstElem.getAdjacentOperator(PileupElement.Direction.NEXT));
    Assert.assertNull(firstElem.getAdjacentOperator(PileupElement.Direction.PREV));
    final PileupElement secondElem = pu.getElementForRead(read2);
    Assert.assertEquals(secondElem.getAdjacentOperator(PileupElement.Direction.NEXT), CigarOperator.I);
    Assert.assertNull(secondElem.getAdjacentOperator(PileupElement.Direction.PREV));
    //checking non-blowup. We're not making any claims about the format of toString
    Assert.assertNotNull(pu.toString());
    Assert.assertEquals(pu.getPileupString('N'), String.format("%s %s N %s%s %s%s", // the position
    loc.getContig(), // the position
    loc.getStart(), // the bases
    (char) read1.getBase(0), // the bases
    (char) read2.getBase(0), // base quality in FASTQ format
    SAMUtils.phredToFastq(read1.getBaseQuality(0)), // base quality in FASTQ format
    SAMUtils.phredToFastq(read2.getBaseQuality(0))));
    Assert.assertEquals(pu.getBases(), new byte[] { (byte) 'A', (byte) 'C' });
    Assert.assertFalse(pu.isEmpty());
    Assert.assertEquals(pu.size(), 2, "size");
    Assert.assertEquals(pu.getBaseCounts(), new int[] { 1, 1, 0, 0 });
    Assert.assertEquals(pu.getBaseQuals(), new byte[] { 10, 20 });
    Assert.assertEquals(pu.getLocation(), loc);
    Assert.assertEquals(pu.getMappingQuals(), new int[] { NO_MAPPING_QUALITY, NO_MAPPING_QUALITY });
    Assert.assertEquals(pu.makeFilteredPileup(p -> p.getQual() >= 12).makeFilteredPileup(p -> p.getMappingQual() >= 0).size(), 1, "getBaseAndMappingFilteredPileup");
    Assert.assertEquals(pu.makeFilteredPileup(r -> r.getQual() >= 12).size(), 1, "getBaseFilteredPileup");
    Assert.assertEquals(pu.getNumberOfElements(p -> true), 2);
    Assert.assertEquals(pu.getNumberOfElements(p -> false), 0);
    Assert.assertEquals(pu.getNumberOfElements(p -> p.isDeletion()), 0);
    Assert.assertEquals(pu.getNumberOfElements(p -> p.isBeforeDeletionStart()), 0);
    Assert.assertEquals(pu.getNumberOfElements(p -> p.isBeforeInsertion()), 1);
    Assert.assertEquals(pu.getNumberOfElements(p -> p.getRead().getMappingQuality() == 0), 2);
    Assert.assertEquals(pu.getOffsets(), Arrays.asList(offset, offset), "getOffsets");
    Assert.assertEquals(pu.getReadGroupIDs(), Arrays.asList(new SAMReadGroupRecord[] { null }), "getReadGroups");
    Assert.assertEquals(pu.getReads(), Arrays.asList(read1, read2), "getReads");
    Assert.assertEquals(pu.getSamples(header), Arrays.asList(new String[] { null }), "getSamples");
    Assert.assertTrue(pu.getPileupForLane("fred").isEmpty());
    Assert.assertTrue(pu.makeFilteredPileup(r -> r.getMappingQual() >= 10).isEmpty());
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) Locatable(htsjdk.samtools.util.Locatable) java.util(java.util) SAMUtils(htsjdk.samtools.SAMUtils) DataProvider(org.testng.annotations.DataProvider) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) CigarElement(htsjdk.samtools.CigarElement) BeforeClass(org.testng.annotations.BeforeClass) CigarOperator(htsjdk.samtools.CigarOperator) QualityUtils(org.broadinstitute.hellbender.utils.QualityUtils) Test(org.testng.annotations.Test) GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) SAMFileHeader(htsjdk.samtools.SAMFileHeader) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) UserException(org.broadinstitute.hellbender.exceptions.UserException) ArtificialReadUtils(org.broadinstitute.hellbender.utils.read.ArtificialReadUtils) Assert(org.testng.Assert) NO_MAPPING_QUALITY(org.broadinstitute.hellbender.utils.read.ReadConstants.NO_MAPPING_QUALITY) Utils(org.broadinstitute.hellbender.utils.Utils) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) SAMFileHeader(htsjdk.samtools.SAMFileHeader) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 52 with SAMReadGroupRecord

use of htsjdk.samtools.SAMReadGroupRecord in project gatk by broadinstitute.

the class ReadPileupUnitTest method testOverlappingFragmentFilter.

@Test
public void testOverlappingFragmentFilter() throws Exception {
    final String rgID = "MY.ID";
    final SAMReadGroupRecord rg = new SAMReadGroupRecord(rgID);
    rg.setSample("MySample");
    final SAMFileHeader header = ArtificialReadUtils.createArtificialSamHeaderWithReadGroup(rg);
    final int readlength = 10;
    final String cigar = "10M";
    final byte[] lowQuals = Utils.repeatBytes((byte) 10, readlength);
    final byte[] highQuals = Utils.repeatBytes((byte) 20, readlength);
    final byte[] basesAllA = Utils.repeatChars('A', readlength);
    final byte[] basesAllC = Utils.repeatChars('C', readlength);
    final byte[] basesAllT = Utils.repeatChars('T', readlength);
    GATKRead read1BasesDisagree = ArtificialReadUtils.createArtificialRead(header, basesAllA, lowQuals, cigar);
    GATKRead read2BasesDisagree = ArtificialReadUtils.createArtificialRead(header, basesAllC, highQuals, cigar);
    GATKRead read3BasesDisagree = ArtificialReadUtils.createArtificialRead(header, basesAllT, lowQuals, cigar);
    read1BasesDisagree.setName("BasesDisagree");
    read2BasesDisagree.setName("BasesDisagree");
    read3BasesDisagree.setName("BasesDisagree");
    read1BasesDisagree.setReadGroup(rgID);
    read2BasesDisagree.setReadGroup(rgID);
    read3BasesDisagree.setReadGroup(rgID);
    GATKRead read1BasesAgree = ArtificialReadUtils.createArtificialRead(header, basesAllA, lowQuals, cigar);
    GATKRead read2BasesAgree = ArtificialReadUtils.createArtificialRead(header, basesAllA, highQuals, cigar);
    read1BasesAgree.setName("BasesAgree");
    read2BasesAgree.setName("BasesAgree");
    read1BasesAgree.setReadGroup(rgID);
    read2BasesAgree.setReadGroup(rgID);
    final List<GATKRead> reads = Arrays.asList(read1BasesDisagree, read2BasesDisagree, read1BasesAgree, read2BasesAgree);
    final int off = 0;
    final ReadPileup pu = new ReadPileup(loc, reads, off);
    final ReadPileup filteredPileupDiscardDiscordant = pu.getOverlappingFragmentFilteredPileup(header);
    final List<GATKRead> filteredReadsDiscardDiscordant = filteredPileupDiscardDiscordant.getReads();
    Assert.assertFalse(filteredReadsDiscardDiscordant.contains(read1BasesDisagree), "Reads with disagreeing bases were kept.");
    Assert.assertFalse(filteredReadsDiscardDiscordant.contains(read2BasesDisagree), "Reads with disagreeing bases were kept.");
    Assert.assertFalse(filteredReadsDiscardDiscordant.contains(read3BasesDisagree), "Reads with disagreeing bases were kept.");
    Assert.assertFalse(filteredReadsDiscardDiscordant.contains(read1BasesAgree), "The lower quality base was kept.");
    Assert.assertTrue(filteredReadsDiscardDiscordant.contains(read2BasesAgree), "The higher quality base is missing.");
    final ReadPileup filteredPileupKeepDiscordant = pu.getOverlappingFragmentFilteredPileup(false, ReadPileup.baseQualTieBreaker, header);
    final List<GATKRead> filteredReadsKeepDiscordant = filteredPileupKeepDiscordant.getReads();
    Assert.assertFalse(filteredReadsKeepDiscordant.contains(read1BasesDisagree), "Low quality base was kept.");
    Assert.assertTrue(filteredReadsKeepDiscordant.contains(read2BasesDisagree), "High quality base is missing.");
    Assert.assertFalse(filteredReadsKeepDiscordant.contains(read3BasesDisagree), "Low quality base was kept.");
    Assert.assertFalse(filteredReadsKeepDiscordant.contains(read1BasesAgree), "The lower quality base was kept.");
    Assert.assertTrue(filteredReadsKeepDiscordant.contains(read2BasesAgree), "The higher quality base is missing.");
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) SAMFileHeader(htsjdk.samtools.SAMFileHeader) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 53 with SAMReadGroupRecord

use of htsjdk.samtools.SAMReadGroupRecord in project gatk by broadinstitute.

the class ReadPileupUnitTest method testSplitByReadGroup.

/**
     * Ensure that basic read group splitting works.
     */
@Test
public void testSplitByReadGroup() {
    final SAMReadGroupRecord readGroupOne = new SAMReadGroupRecord("rg1");
    final SAMReadGroupRecord readGroupTwo = new SAMReadGroupRecord("rg2");
    final SAMFileHeader header = ArtificialReadUtils.createArtificialSamHeader(1, 1, 1000);
    header.addReadGroup(readGroupOne);
    header.addReadGroup(readGroupTwo);
    final GATKRead read1 = ArtificialReadUtils.createArtificialRead(header, "read1", 0, 1, 10);
    read1.setReadGroup(readGroupOne.getId());
    final GATKRead read2 = ArtificialReadUtils.createArtificialRead(header, "read2", 0, 1, 10);
    read2.setReadGroup(readGroupTwo.getId());
    final GATKRead read3 = ArtificialReadUtils.createArtificialRead(header, "read3", 0, 1, 10);
    read3.setReadGroup(readGroupOne.getId());
    final GATKRead read4 = ArtificialReadUtils.createArtificialRead(header, "read4", 0, 1, 10);
    read4.setReadGroup(readGroupTwo.getId());
    final GATKRead read5 = ArtificialReadUtils.createArtificialRead(header, "read5", 0, 1, 10);
    read5.setReadGroup(readGroupTwo.getId());
    final GATKRead read6 = ArtificialReadUtils.createArtificialRead(header, "read6", 0, 1, 10);
    read6.setReadGroup(readGroupOne.getId());
    final GATKRead read7 = ArtificialReadUtils.createArtificialRead(header, "read7", 0, 1, 10);
    read7.setReadGroup(readGroupOne.getId());
    final ReadPileup pileup = new ReadPileup(loc, Arrays.asList(read1, read2, read3, read4, read5, read6, read7), 1);
    final ReadPileup rg1Pileup = pileup.makeFilteredPileup(pe -> "rg1".equals(pe.getRead().getReadGroup()));
    final List<GATKRead> rg1Reads = rg1Pileup.getReads();
    Assert.assertEquals(rg1Reads.size(), 4, "Wrong number of reads in read group rg1");
    Assert.assertEquals(rg1Reads.get(0), read1, "Read " + read1.getName() + " should be in rg1 but isn't");
    Assert.assertEquals(rg1Reads.get(1), read3, "Read " + read3.getName() + " should be in rg1 but isn't");
    Assert.assertEquals(rg1Reads.get(2), read6, "Read " + read6.getName() + " should be in rg1 but isn't");
    Assert.assertEquals(rg1Reads.get(3), read7, "Read " + read7.getName() + " should be in rg1 but isn't");
    final ReadPileup rg2Pileup = pileup.makeFilteredPileup(pe -> "rg2".equals(pe.getRead().getReadGroup()));
    final List<GATKRead> rg2Reads = rg2Pileup.getReads();
    Assert.assertEquals(rg2Reads.size(), 3, "Wrong number of reads in read group rg2");
    Assert.assertEquals(rg2Reads.get(0), read2, "Read " + read2.getName() + " should be in rg2 but isn't");
    Assert.assertEquals(rg2Reads.get(1), read4, "Read " + read4.getName() + " should be in rg2 but isn't");
    Assert.assertEquals(rg2Reads.get(2), read5, "Read " + read5.getName() + " should be in rg2 but isn't");
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) SAMFileHeader(htsjdk.samtools.SAMFileHeader) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 54 with SAMReadGroupRecord

use of htsjdk.samtools.SAMReadGroupRecord in project gatk by broadinstitute.

the class NGSPlatformUnitTest method testPLFromReadWithRG.

/**
     * A unit test that creates an artificial read for testing some code that uses reads
     */
@Test(dataProvider = "TestMappings")
public void testPLFromReadWithRG(final String plField, final NGSPlatform expected) {
    final SAMFileHeader header = ArtificialReadUtils.createArtificialSamHeader(seq.getSequenceDictionary());
    final String rgID = "ID";
    final SAMReadGroupRecord rg = new SAMReadGroupRecord(rgID);
    if (plField != null)
        rg.setPlatform(plField);
    header.addReadGroup(rg);
    final GATKRead read = ArtificialReadUtils.createArtificialRead(header, "myRead", 0, 1, 10);
    read.setAttribute("RG", rgID);
    Assert.assertEquals(NGSPlatform.fromRead(read, header), expected);
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) SAMFileHeader(htsjdk.samtools.SAMFileHeader) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 55 with SAMReadGroupRecord

use of htsjdk.samtools.SAMReadGroupRecord in project gatk by broadinstitute.

the class ReadGroupBlackListReadFilterUnitTest method testFilterOutByListFile.

@Test
public void testFilterOutByListFile() throws IOException {
    int recordsPerGroup = 3;
    List<GATKRead> records = new ArrayList<>();
    int alignmentStart = 0;
    for (int x = 0; x < GROUP_COUNT; x++) {
        SAMReadGroupRecord groupRecord = header.getReadGroup(getReadGroupId(x));
        for (int y = 1; y <= recordsPerGroup; y++) {
            GATKRead record = ArtificialReadUtils.createArtificialRead(header, "readUno", 0, ++alignmentStart, 20);
            record.setReadGroup(groupRecord.getReadGroupId());
            records.add(record);
        }
    }
    List<String> filterList = new ArrayList<>();
    filterList.add(publicTestDir + "readgroupblacklisttestlist.txt");
    ReadGroupBlackListReadFilter filter = new ReadGroupBlackListReadFilter(filterList, header);
    int filtered = 0;
    int unfiltered = 0;
    for (GATKRead record : records) {
        String readGroup = record.getReadGroup();
        if (!filter.test(record)) {
            if (!("ReadGroup3".equals(readGroup) || "ReadGroup4".equals(readGroup)))
                Assert.fail("Read group " + readGroup + " was filtered");
            filtered++;
        } else {
            if ("ReadGroup3".equals(readGroup) || "ReadGroup4".equals(readGroup))
                Assert.fail("Read group " + readGroup + " was not filtered");
            unfiltered++;
        }
    }
    int filteredExpected = recordsPerGroup * 2;
    int unfilteredExpected = recordsPerGroup * (GROUP_COUNT - 2);
    Assert.assertEquals(filtered, filteredExpected, "Filtered");
    Assert.assertEquals(unfiltered, unfilteredExpected, "Unfiltered");
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) ArrayList(java.util.ArrayList) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Aggregations

SAMReadGroupRecord (htsjdk.samtools.SAMReadGroupRecord)81 SAMFileHeader (htsjdk.samtools.SAMFileHeader)48 SAMRecord (htsjdk.samtools.SAMRecord)33 Test (org.testng.annotations.Test)31 SamReader (htsjdk.samtools.SamReader)29 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)26 File (java.io.File)23 ArrayList (java.util.ArrayList)22 SAMRecordIterator (htsjdk.samtools.SAMRecordIterator)20 GATKRead (org.broadinstitute.hellbender.utils.read.GATKRead)20 HashMap (java.util.HashMap)18 CigarElement (htsjdk.samtools.CigarElement)17 Cigar (htsjdk.samtools.Cigar)16 HashSet (java.util.HashSet)16 SAMFileWriter (htsjdk.samtools.SAMFileWriter)15 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)15 CigarOperator (htsjdk.samtools.CigarOperator)14 IOException (java.io.IOException)14 SAMSequenceDictionaryProgress (com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress)13 List (java.util.List)12