@Test(enabled = !DEBUG, dataProvider = "NumHardClipped") public void testNumHardClipped(final Cigar cigar, final int expected) { final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "myRead", 0, 1, cigar == null ? 10 : cigar.getReadLength()); read.setCigar(cigar); Assert.assertEquals(AlignmentUtils.getNumHardClippedBases(read), expected, "Cigar " + cigar + " failed num hard clips"); }
@Test(enabled = !DEBUG, dataProvider = "NumAlignedBasesCountingSoftClips") public void testNumAlignedBasesCountingSoftClips(final Cigar cigar, final int expected) { final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "myRead", 0, 1, cigar == null ? 10 : cigar.getReadLength()); read.setCigar(cigar); Assert.assertEquals(AlignmentUtils.getNumAlignedBasesCountingSoftClips(read), expected, "Cigar " + cigar + " failed NumAlignedBasesCountingSoftClips"); }
@Test(enabled = !DEBUG, dataProvider = "NumAlignedBlocks") public void testNumAlignedBlocks(final Cigar cigar, final int expected) { final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "myRead", 0, 1, cigar == null ? 10 : cigar.getReadLength()); read.setCigar(cigar); Assert.assertEquals(AlignmentUtils.getNumAlignmentBlocks(read), expected, "Cigar " + cigar + " failed NumAlignedBlocks"); }
public boolean finalizeUpdate() { // if we haven't made any changes, don't do anything if ( newCigar == null ) return false; if ( newStart == -1 ) newStart = read.getAlignmentStart(); else if ( Math.abs(newStart - read.getAlignmentStart()) > MAX_POS_MOVE_ALLOWED ) { logger.debug(String.format("Attempting to realign read %s at %d more than %d bases to %d.", read.getReadName(), read.getAlignmentStart(), MAX_POS_MOVE_ALLOWED, newStart)); return false; } // store the old CIGAR and start in case we need to back out final Cigar oldCigar = read.getCigar(); final int oldStart = read.getAlignmentStart(); // try updating the read with the new CIGAR and start read.setCigar(newCigar); read.setAlignmentStart(newStart); // back out if necessary if ( realignmentProducesBadAlignment(read) ) { read.setCigar(oldCigar); read.setAlignmentStart(oldStart); return false; } // annotate the record with the original cigar and start (if it changed) if ( !NO_ORIGINAL_ALIGNMENT_TAGS ) { read.setAttribute(ORIGINAL_CIGAR_TAG, oldCigar.toString()); if ( newStart != oldStart ) read.setAttribute(ORIGINAL_POSITION_TAG, oldStart); } return true; }
@Override public GATKSAMRecord apply(final GATKSAMRecord read) { if(read == null) throw new UserException.BadInput("try to transform a null GATKSAMRecord"); final Cigar originalCigar = read.getCigar(); if (originalCigar.isValid(read.getReadName(),-1) != null) throw new UserException.BadInput("try to transform a read with non-valid cigar string: readName: "+read.getReadName()+" Cigar String: "+originalCigar); read.setCigar(refactorNDNtoN(originalCigar)); return read; }
List<CigarElement> elements = new ArrayList<CigarElement>(); elements.add(new CigarElement(length, CigarOperator.characterToEnum('M'))); record.setCigar(new Cigar(elements)); record.setProperPairFlag(false);
hardClippedRead.setCigar(cigarShift.cigar); if (start == 0) hardClippedRead.setAlignmentStart(read.getAlignmentStart() + calculateAlignmentStartShift(read.getCigar(), cigarShift.cigar));
public Integer map(ReferenceContext ref, GATKSAMRecord read, RefMetaDataTracker metaDataTracker) { // we can not deal with screwy records if ( read.getReadUnmappedFlag() || read.getCigar().numCigarElements() == 0 ) { emit(read); return 0; } // move existing indels (for 1 indel reads only) to leftmost position within identical sequence int numBlocks = AlignmentUtils.getNumAlignmentBlocks(read); if ( numBlocks == 2 ) { Cigar newCigar = AlignmentUtils.leftAlignIndel(IndelRealigner.unclipCigar(read.getCigar()), ref.getBases(), read.getReadBases(), 0, 0, true); newCigar = IndelRealigner.reclipCigar(newCigar, read); read.setCigar(newCigar); } emit(read); return 1; }
originalRead.getReadBases(), swPairwiseAlignment.getAlignmentStart2wrt1(), 0, true); read.setCigar(readToRefCigar);
read.setCigar(newCigar);
final ArrayList<CigarElement> cList = new ArrayList<CigarElement>(); cList.add(c); returnRead.setCigar( new Cigar( cList )); returnRead.setMappingQuality( firstRead.getMappingQuality() );
protected GATKSAMRecord buildSAMRecord(String readName, String contig, int alignmentStart, int alignmentEnd) { SAMFileHeader header = ArtificialSAMUtils.createDefaultReadGroup(new SAMFileHeader(), "test", "test"); header.setSequenceDictionary(dictionary); header.setSortOrder(SAMFileHeader.SortOrder.coordinate); GATKSAMRecord record = new GATKSAMRecord(header); record.setReadName(readName); record.setReferenceIndex(dictionary.getSequenceIndex(contig)); record.setAlignmentStart(alignmentStart); Cigar cigar = new Cigar(); int len = alignmentEnd - alignmentStart + 1; cigar.add(new CigarElement(len, CigarOperator.M)); record.setCigar(cigar); record.setReadString(new String(new char[len]).replace("\0", "A")); record.setBaseQualities(new byte[len]); record.setReadGroup(new GATKSAMReadGroupRecord(header.getReadGroup("test"))); return record; }
private GATKSAMRecord revertSoftClippedBases(GATKSAMRecord read) { GATKSAMRecord unclipped = (GATKSAMRecord) read.clone(); Cigar unclippedCigar = new Cigar(); int matchesCount = 0; for (CigarElement element : read.getCigar().getCigarElements()) { if (element.getOperator() == CigarOperator.SOFT_CLIP || element.getOperator() == CigarOperator.MATCH_OR_MISMATCH) matchesCount += element.getLength(); else if (matchesCount > 0) { unclippedCigar.add(new CigarElement(matchesCount, CigarOperator.MATCH_OR_MISMATCH)); matchesCount = 0; unclippedCigar.add(element); } else unclippedCigar.add(element); } if (matchesCount > 0) unclippedCigar.add(new CigarElement(matchesCount, CigarOperator.MATCH_OR_MISMATCH)); unclipped.setCigar(unclippedCigar); final int newStart = read.getAlignmentStart() + calculateAlignmentStartShift(read.getCigar(), unclippedCigar); unclipped.setAlignmentStart(newStart); if ( newStart <= 0 ) { // if the start of the unclipped read occurs before the contig, // we must hard clip away the bases since we cannot represent reads with // negative or 0 alignment start values in the SAMRecord (e.g., 0 means unaligned) return hardClip(unclipped, 0, - newStart); } else { return unclipped; } }
/** * Write out a representation of this haplotype as a read * * @param haplotype a haplotype to write out. Cannot be null * @param paddedRefLoc the reference location. Cannot be null * @param isAmongBestHaplotypes true if among the best haplotypes, false if it was just one possible but not so good */ private void writeHaplotype(final Haplotype haplotype, final GenomeLoc paddedRefLoc, final boolean isAmongBestHaplotypes) { final GATKSAMRecord record = new GATKSAMRecord(output.getHeader()); record.setReadBases(haplotype.getBases()); record.setAlignmentStart(paddedRefLoc.getStart() + haplotype.getAlignmentStartHapwrtRef()); record.setBaseQualities(Utils.dupBytes((byte) '!', haplotype.getBases().length)); record.setCigar(AlignmentUtils.consolidateCigar(haplotype.getCigar())); record.setMappingQuality(isAmongBestHaplotypes ? 60 : 0); record.setReadName("HC" + uniqueNameCounter++); record.setAttribute(AlignmentUtils.HAPLOTYPE_TAG,haplotype.hashCode()); record.setReadUnmappedFlag(false); record.setReferenceIndex(paddedRefLoc.getContigIndex()); record.setAttribute(SAMTag.RG.toString(), READ_GROUP_ID); record.setFlags(16); output.add(record); }
/** * Build a SAM record featuring the absolute minimum required dataset. * * @param contig Contig to populate. * @param alignmentStart start of alignment * @param alignmentEnd end of alignment * * @return New SAM Record */ protected GATKSAMRecord buildSAMRecord(String readName, String contig, int alignmentStart, int alignmentEnd) { SAMFileHeader header = new SAMFileHeader(); header.setSequenceDictionary(sequenceSourceFile.getSequenceDictionary()); GATKSAMRecord record = new GATKSAMRecord(header); record.setReadName(readName); record.setReferenceIndex(sequenceSourceFile.getSequenceDictionary().getSequenceIndex(contig)); record.setAlignmentStart(alignmentStart); Cigar cigar = new Cigar(); int len = alignmentEnd - alignmentStart + 1; cigar.add(new CigarElement(len, CigarOperator.M)); record.setCigar(cigar); record.setReadBases(new byte[len]); record.setBaseQualities(new byte[len]); return record; }
private GATKSAMRecord createReadOffContig(final SAMFileHeader header, final boolean negStrand, final int pre, final int post) { final int contigLen = header.getSequence(0).getSequenceLength(); final int readLen = pre + contigLen + post; final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "read1", 0, 1, readLen); read.setAlignmentStart(1); read.setCigar(TextCigarCodec.decode(pre + "S" + contigLen + "M" + post + "S")); read.setBaseQualities(Utils.dupBytes((byte) 30, readLen)); read.setReadBases(Utils.dupBytes((byte)'A', readLen)); read.setMappingQuality(60); read.setMateAlignmentStart(1); read.setProperPairFlag(true); read.setReadPairedFlag(true); read.setInferredInsertSize(30); read.setReadNegativeStrandFlag(negStrand); read.setMateNegativeStrandFlag(! negStrand); read.setReadGroup(new GATKSAMReadGroupRecord("foo")); return read; }