/** * Create an artificial read based on the parameters * * @param header the SAM header to associate the read with * @param name the name of the read * @param refIndex the reference index, i.e. what chromosome to associate it with * @param alignmentStart where to start the alignment * @param bases the sequence of the read * @param qual the qualities of the read * @param cigar the cigar string of the read * @return the artificial read */ public static GATKSAMRecord createArtificialRead(SAMFileHeader header, String name, int refIndex, int alignmentStart, byte[] bases, byte[] qual, String cigar) { GATKSAMRecord rec = createArtificialRead(header, name, refIndex, alignmentStart, bases, qual); rec.setCigarString(cigar); return rec; }
private GATKSAMRecord createRead(final boolean useGoodBases) { GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "foo", 0, 10, readBases.getBytes(), useGoodBases ? Arrays.copyOf(goodQuals, goodQuals.length) : Arrays.copyOf(badQuals, badQuals.length)); read.setCigarString("10M"); return read; }
@Test(dataProvider = "ClipUpstreamProvider", enabled = true) public void clipUpstreamTest(final int readStart, final int readLength, final int delLength) { final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "basicRead", 0, readStart, readLength); if ( delLength == 0 ) read.setCigarString(readLength + "M"); else read.setCigarString((readLength / 2) + "M" + delLength + "D" + (readLength / 2) + "M"); final boolean result = PairHMMIndelErrorModel.mustClipUpstream(read, refWindowStart); Assert.assertEquals(result, read.getSoftStart() < refWindowStart && read.getSoftEnd() > refWindowStart); }
@Test(dataProvider = "ClipDownstreamProvider", enabled = true) public void clipDownstreamTest(final int readStart, final int readLength, final int delLength) { final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "basicRead", 0, readStart, readLength); if ( delLength == 0 ) read.setCigarString(readLength + "M"); else read.setCigarString((readLength / 2) + "M" + delLength + "D" + (readLength / 2) + "M"); final boolean result = PairHMMIndelErrorModel.mustClipDownstream(read, refWindowEnd); Assert.assertEquals(result, read.getSoftStart() < refWindowEnd && read.getSoftStart() + readLength > refWindowEnd); }
@Test public void realignAtContigBorderTest() { final int contigEnd = header.getSequence(0).getSequenceLength(); final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "goodRead", 0, contigEnd - 1, 2); read.setCigarString("2M"); Assert.assertEquals(IndelRealigner.realignmentProducesBadAlignment(read, contigEnd), false); read.setCigarString("1M1D1M"); Assert.assertEquals(IndelRealigner.realignmentProducesBadAlignment(read, contigEnd), true); }
@Test public void clipDownstreamAtBorderTest() { final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "basicRead", 0, 5, 10); read.setCigarString("10M"); Assert.assertEquals(PairHMMIndelErrorModel.mustClipDownstream(read, 13), true); Assert.assertEquals(PairHMMIndelErrorModel.mustClipDownstream(read, 14), false); }
private Cigar makeExpectedCigar2(final Cigar originalCigar, final char indelOp, final int indelStart, final int indelSize, final int readLength) { if ( indelStart < 17 || indelStart > (26 - (indelOp == 'D' ? indelSize : 0)) ) return originalCigar; final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "myRead", 0, 1, readLength); if ( indelOp == 'I' && (indelSize == 1 || indelSize == 3) && indelStart % 2 == 1 ) read.setCigarString(buildTestCigarString(indelOp, 0, Math.max(indelStart - indelSize, 16), indelSize, readLength)); else if ( (indelSize == 2 || indelSize == 4) && (indelOp == 'D' || indelStart % 2 == 0) ) read.setCigarString(buildTestCigarString(indelOp, 0, 16, indelSize, readLength)); else return originalCigar; return read.getCigar(); }
private Cigar makeExpectedCigar3(final Cigar originalCigar, final char indelOp, final int indelStart, final int indelSize, final int readLength) { if ( indelStart < 17 || indelStart > (28 - (indelOp == 'D' ? indelSize : 0)) ) return originalCigar; final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "myRead", 0, 1, readLength); if ( indelSize == 3 && (indelOp == 'D' || indelStart % 3 == 1) ) read.setCigarString(buildTestCigarString(indelOp, 0, 16, indelSize, readLength)); else if ( (indelOp == 'I' && indelSize == 4 && indelStart % 3 == 2) || (indelOp == 'I' && indelSize == 2 && indelStart % 3 == 0) || (indelOp == 'I' && indelSize == 1 && indelStart < 28 && indelStart % 3 == 2) ) read.setCigarString(buildTestCigarString(indelOp, 0, Math.max(indelStart - indelSize, 16), indelSize, readLength)); else return originalCigar; return read.getCigar(); }
private Cigar makeExpectedCigar1(final Cigar originalCigar, final char indelOp, final int indelStart, final int indelSize, final int readLength) { if ( indelSize == 0 || indelStart < 17 || indelStart > (26 - (indelOp == 'D' ? indelSize : 0)) ) return originalCigar; final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "myRead", 0, 1, readLength); read.setCigarString(buildTestCigarString(indelOp, 0, 16, indelSize, readLength)); return read.getCigar(); }
public GATKSAMRecord createReadAndAddToBin(String cigarString, int alignmentStart) { final GATKSAMRecord read = ReadUtils.createRandomRead(readLength); read.setCigarString(cigarString); read.setAlignmentStart(alignmentStart); readBin.add(read); return read; } }
@Test (enabled = true) public void testReadWithNsRefAfterDeletion() throws FileNotFoundException { final ReferenceSequenceFile seq = new CachingIndexedFastaSequenceFile(new File(b37KGReference)); final SAMFileHeader header = ArtificialSAMUtils.createArtificialSamHeader(seq.getSequenceDictionary()); final int readLength = 76; final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "myRead", 0, 8975, readLength); read.setReadBases(Utils.dupBytes((byte) 'A', readLength)); read.setBaseQualities(Utils.dupBytes((byte)30, readLength)); read.setCigarString("3M414N1D73M"); final int result = ReadUtils.getReadCoordinateForReferenceCoordinateUpToEndOfRead(read, 9393, ReadUtils.ClippingTail.LEFT_TAIL); Assert.assertEquals(result, 3); }
@Test (enabled = true) public void testReadWithNsRefIndexInDeletion() throws FileNotFoundException { final ReferenceSequenceFile seq = new CachingIndexedFastaSequenceFile(new File(b37KGReference)); final SAMFileHeader header = ArtificialSAMUtils.createArtificialSamHeader(seq.getSequenceDictionary()); final int readLength = 76; final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "myRead", 0, 8975, readLength); read.setReadBases(Utils.dupBytes((byte) 'A', readLength)); read.setBaseQualities(Utils.dupBytes((byte)30, readLength)); read.setCigarString("3M414N1D73M"); final int result = ReadUtils.getReadCoordinateForReferenceCoordinateUpToEndOfRead(read, 9392, ReadUtils.ClippingTail.LEFT_TAIL); Assert.assertEquals(result, 2); }
@Test( ) public void testCreationFromSAMRecordUnmappedButOnGenome() { final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "foo", 0, 1, 5); read.setReadUnmappedFlag(true); read.setCigarString("*"); final GenomeLoc loc = genomeLocParser.createGenomeLoc(read); Assert.assertEquals(loc.getContig(), read.getReferenceName()); Assert.assertEquals(loc.getContigIndex(), (int)read.getReferenceIndex()); Assert.assertEquals(loc.getStart(), read.getAlignmentStart()); Assert.assertEquals(loc.getStop(), read.getAlignmentStart()); }
public GATKSAMRecord makeRead() { GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "read", 0, locus, readLength); read.setReadBases(Utils.dupBytes((byte) 'A', readLength)); final byte[] quals = new byte[readLength]; for ( int i = 0; i < readLength; i++ ) quals[i] = (byte)(i % QualityUtils.MAX_SAM_QUAL_SCORE); read.setBaseQualities(quals); read.setCigarString(cigarString); return read; } }
@Test public void testRBPMappingQuals() { // create a read with high MQ final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "read", 0, 1, 10); read.setReadBases(Utils.dupBytes((byte) 'A', 10)); read.setBaseQualities(Utils.dupBytes((byte) 30, 10)); read.setCigarString("10M"); read.setMappingQuality(200); // set a MQ higher than max signed byte // now create the RBP final List<PileupElement> elts = new LinkedList<>(); elts.add(new PileupElement(read, 0, read.getCigar().getCigarElement(0), 0, 0)); final Map<String, ReadBackedPileupImpl> pileupsBySample = new HashMap<>(); pileupsBySample.put("foo", new ReadBackedPileupImpl(loc, elts)); final ReadBackedPileup pileup = new ReadBackedPileupImpl(loc, pileupsBySample); Assert.assertEquals(pileup.getMappingQuals()[0], 200); } }
@Test(enabled = false && ! DEBUG) public void testWholeIndelReadInIsolation() { final int firstLocus = 44367789; GATKSAMRecord indelOnlyRead = ArtificialSAMUtils.createArtificialRead(header, "indelOnly", 0, firstLocus, 76); indelOnlyRead.setReadBases(Utils.dupBytes((byte)'A',76)); indelOnlyRead.setBaseQualities(Utils.dupBytes((byte) '@', 76)); indelOnlyRead.setCigarString("76I"); List<GATKSAMRecord> reads = Arrays.asList(indelOnlyRead); // create the iterator by state with the fake reads and fake records li = makeLTBS(reads); // Traditionally, reads that end with indels bleed into the pileup at the following locus. Verify that the next pileup contains this read // and considers it to be an indel-containing read. Assert.assertTrue(li.hasNext(),"Should have found a whole-indel read in the normal base pileup without extended events enabled"); AlignmentContext alignmentContext = li.next(); Assert.assertEquals(alignmentContext.getLocation().getStart(), firstLocus, "Base pileup is at incorrect location."); ReadBackedPileup basePileup = alignmentContext.getBasePileup(); Assert.assertEquals(basePileup.getReads().size(),1,"Pileup is of incorrect size"); Assert.assertSame(basePileup.getReads().get(0), indelOnlyRead, "Read in pileup is incorrect"); }
private List<PileupElement> makeReads( final int n, final int mapq, final String op ) { final int readLength = 3; final List<PileupElement> elts = new LinkedList<PileupElement>(); for ( int i = 0; i < n; i++ ) { GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "read", 0, 1, readLength); read.setReadBases(Utils.dupBytes((byte) 'A', readLength)); read.setBaseQualities(Utils.dupBytes((byte) 30, readLength)); read.setCigarString("1M1" + op + "1M"); read.setMappingQuality(mapq); final int baseOffset = op.equals("M") ? 1 : 0; final CigarElement cigarElement = read.getCigar().getCigarElement(1); elts.add(new PileupElement(read, baseOffset, cigarElement, 1, 0)); } return elts; }
@Test public void testUnmappedReadsDoNotFail() { // create an unmapped read final GATKSAMRecord read = new GATKSAMRecord(ArtificialSAMUtils.createArtificialSamHeader()); read.setReadName("foo"); read.setReferenceName("*"); read.setAlignmentStart(100); read.setCigarString("*"); read.setReadUnmappedFlag(true); // try to add it to the manager final OverhangFixingManager manager = new OverhangFixingManager(null, null, null, 100, 1, 30, false); manager.addRead(read); // we just want to make sure that the following call does not fail Assert.assertTrue(true); } }
@Override protected void setUp() { header = ArtificialSAMUtils.createArtificialSamHeader(1, 1, 1000); genomeLocParser = new GenomeLocParser(header.getSequenceDictionary()); for ( int j = 0; j < nReads; j++ ) { GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "read", 0, locus, readLength); read.setReadBases(Utils.dupBytes((byte) 'A', readLength)); final byte[] quals = new byte[readLength]; for ( int i = 0; i < readLength; i++ ) quals[i] = (byte)(i % QualityUtils.MAX_SAM_QUAL_SCORE); read.setBaseQualities(quals); read.setCigarString(cigar); reads.add(read); } }
protected GATKSAMRecord buildSAMRecord(final String readName, final String contig, final int alignmentStart) { GATKSAMRecord record = new GATKSAMRecord(header); record.setReadName(readName); record.setReferenceIndex(dictionary.getSequenceIndex(contig)); record.setAlignmentStart(alignmentStart); record.setCigarString("1M"); record.setReadString("A"); record.setBaseQualityString("A"); record.setReadGroup(readGroup); return record; }