Java Code Examples for htsjdk.samtools.SAMRecord#setCigar()

SAMRecord firstSAMRecord(SAMFileHeader header) {
	SAMRecord r = new SAMRecord(header);
	r.setReadName(evidenceRecord.getReadName());
	r.setReferenceName(evidenceRecord.Chromosome);
	r.setAlignmentStart(Integer.valueOf(evidenceRecord.OffsetInReference) + 1);
	r.setMappingQuality(Integer.valueOf(evidenceRecord.ScoreAllele0));
	r.setReadPairedFlag(true);
	r.setReadUnmappedFlag(false);
	r.setReadNegativeStrandFlag(negative);
	r.setFirstOfPairFlag(evidenceRecord.side == 0);
	r.setSecondOfPairFlag(!r.getFirstOfPairFlag());

	r.setCigar(new Cigar(Utils.toCigarOperatorList(firstHalf.samCigarElements)));

	r.setReadBases(Utils.toByteArray(firstHalf.readBasesBuf));
	r.setBaseQualities(Utils.toByteArray(firstHalf.readScoresBuf));

	r.setAttribute("GC", Utils.toString(firstHalf.gcList));
	r.setAttribute("GS", Utils.toString(firstHalf.gsBuf));
	r.setAttribute("GQ", SAMUtils.phredToFastq(Utils.toByteArray(firstHalf.gqBuf)));

	return r;
}
 

SAMRecord secondSAMRecord(SAMFileHeader header) {
	SAMRecord r = new SAMRecord(header);
	r.setReadName(evidenceRecord.getReadName());
	r.setReferenceName(evidenceRecord.Chromosome);
	r.setAlignmentStart(Integer.valueOf(evidenceRecord.MateOffsetInReference) + 1);
	r.setMappingQuality(Integer.valueOf(evidenceRecord.ScoreAllele0));
	r.setReadPairedFlag(true);
	r.setReadUnmappedFlag(false);
	r.setReadNegativeStrandFlag(negative);
	r.setFirstOfPairFlag(evidenceRecord.side == 1);
	r.setSecondOfPairFlag(!r.getFirstOfPairFlag());

	r.setCigar(new Cigar(Utils.toCigarOperatorList(secondHalf.samCigarElements)));

	r.setReadBases(Utils.toByteArray(secondHalf.readBasesBuf));
	r.setBaseQualities(Utils.toByteArray(secondHalf.readScoresBuf));

	r.setAttribute("GC", Utils.toString(secondHalf.gcList));
	r.setAttribute("GS", Utils.toString(secondHalf.gsBuf));
	r.setAttribute("GQ", SAMUtils.phredToFastq(Utils.toByteArray(secondHalf.gqBuf)));

	return r;
}
 

@Test
public void testGetConsistentExons () {
	// need to build a new read that falls into 2 genes to check if it's inconsistent.
	AnnotationUtils u = AnnotationUtils.getInstance();
	OverlapDetector<Gene> geneOverlapDetector = getGeneOverlapDetector(BAM, GTF);
	Map<String, Gene> map = getGeneMap(geneOverlapDetector);
	// expected UTR and coding.
	SAMRecord rec = getReadByName("000000000-AMY9M:1:2104:20987:12124", BAM);
	// change the cigar so the read maps to two genes exons.
	// rec.setCigar(cigar);
	rec.setAlignmentStart(6269430);
	Cigar c = new Cigar ();
	c.add(new CigarElement(20, CigarOperator.M));
	c.add(new CigarElement(9741377, CigarOperator.N));
	c.add(new CigarElement(30, CigarOperator.M));
	rec.setCigar(c);
	Set<Gene> genes = new HashSet<>();
	genes.add(map.get("RPL22"));
	genes.add(map.get("PLEKHM2"));

	// This read spans 2 genes, so it's not consistent under the strict test.
	Set<Gene> actual = u.getConsistentExons(rec, genes, false);
	Set<Gene> expected = new HashSet<>();
	Assert.assertEquals(actual, expected);

	// under the less strict option, we allow the read to span 2 different genes.
	actual = u.getConsistentExons(rec, genes, true);
	expected = genes;
	Assert.assertEquals(actual, expected);

}
 

/**
 * Replace hard clips with soft clips.
 */
public static void replaceHardClips(SAMRecord read) {
	Cigar cigar = read.getCigar();
	
	if (cigar.getCigarElements().size() > 0) {
		CigarElement firstElement = cigar.getCigarElement(0);
		CigarElement lastElement  = cigar.getCigarElement(cigar.numCigarElements()-1);
		
		if ((firstElement.getOperator() == CigarOperator.H) ||
			(lastElement.getOperator() == CigarOperator.H)) {
			
			Cigar newCigar = new Cigar();
			
			boolean isFirst = true;
			
			for (CigarElement element : cigar.getCigarElements()) {
				if (element.getOperator() != CigarOperator.H) {
					newCigar.add(element);
				} else {
					CigarElement softClip = new CigarElement(element.getLength(), CigarOperator.S);
					newCigar.add(softClip);
					
					if (isFirst) {
						read.setReadString(padBases(element.getLength()) + read.getReadString());
					} else {
						read.setReadString(read.getReadString() + padBases(element.getLength()));							
					}
				}
				
				isFirst = false;
			}
			
			read.setCigar(newCigar);
		}
	}
}
 

/**
 * Remove leading or trailing soft clips from the input read.
 * Does not modify a read entirely comprised of soft clips.
 */
public static void removeSoftClips(SAMRecord read) {
	
	Cigar cigar = read.getCigar();
	
	CigarElement firstElement = cigar.getCigarElement(0);
	CigarElement lastElement  = cigar.getCigarElement(cigar.numCigarElements()-1);
	
	if ((firstElement.getOperator() == CigarOperator.S) ||
		(lastElement.getOperator() == CigarOperator.S)) {
	
		Cigar newCigar = new Cigar();
		
		String bases = read.getReadString();
		//String qualities = read.getBaseQualityString();
				
		if (firstElement.getOperator() == CigarOperator.S) {
			bases = bases.substring(firstElement.getLength(), bases.length());
			//qualities = qualities.substring(firstElement.getLength(), qualities.length()-1);
		} else {
			newCigar.add(firstElement);
		}
		
		for (int i=1; i<cigar.numCigarElements()-1; i++) {
			newCigar.add(cigar.getCigarElement(i));
		}
		
		if (lastElement.getOperator() == CigarOperator.S) {
			bases = bases.substring(0, bases.length() - lastElement.getLength());
			//qualities = qualities.substring(0, qualities.length() - lastElement.getLength() - 1);
		} else {
			newCigar.add(lastElement);
		}
		
		read.setCigar(newCigar);
		read.setReadString(bases);
		//read.setBaseQualityString(qualities);
	}
}
 

@Test (groups = "unit")
public void testGetLeadingClips() {
	SAMRecord read = new SAMRecord(null);
	Cigar cigar = getCigar("10S15M8S");
	read.setCigar(cigar);
	String leading = SAMRecordUtils.getLeadingClips(read);
	Assert.assertEquals(leading, "10S");
}
 

@Test (groups = "unit")
public void testGetLeadingClips_empty() {
	SAMRecord read = new SAMRecord(null);
	Cigar cigar = getCigar("15M8S");
	read.setCigar(cigar);
	String leading = SAMRecordUtils.getLeadingClips(read);
	Assert.assertEquals(leading, "");
}
 

@Test (groups = "unit")
public void testGetTrailingClips() {
	SAMRecord read = new SAMRecord(null);
	Cigar cigar = getCigar("10S15M8S");
	read.setCigar(cigar);
	String leading = SAMRecordUtils.getTrailingClips(read);
	Assert.assertEquals(leading, "8S");
}
 

@Test (groups = "unit")
public void testGetTrailingClips_empty() {
	SAMRecord read = new SAMRecord(null);
	Cigar cigar = getCigar("10S15M");
	read.setCigar(cigar);
	String leading = SAMRecordUtils.getTrailingClips(read);
	Assert.assertEquals(leading, "");
}
 

@Test (groups = "unit")
public void testGetMappedReadPortion() {
	SAMRecord read = new SAMRecord(null);
	Cigar cigar = getCigar("4S5M7S");
	String seq = "ATCGCGATACCCCCCC";
	read.setCigar(cigar);
	read.setReadString(seq);
	
	String unclipped = SAMRecordUtils.getMappedReadPortion(read);
	Assert.assertEquals(unclipped, "CGATA");
}
 

/** Replaces hard clips with soft clips and fills in bases and qualities with dummy values as needed. */
 private void replaceHardWithSoftClips(final SAMRecord rec) {
     if (rec.getReadUnmappedFlag()) return;
     if (rec.getCigar().isEmpty()) return;

     List<CigarElement> elements = rec.getCigar().getCigarElements();
     final CigarElement first = elements.get(0);
     final CigarElement last  = elements.size() == 1 ? null : elements.get(elements.size()-1);
     final int startHardClip = first.getOperator() == CigarOperator.H ? first.getLength() : 0;
     final int endHardClip   = (last != null && last.getOperator() == CigarOperator.H) ? last.getLength() : 0;

     if (startHardClip + endHardClip > 0) {
         final int len = rec.getReadBases().length + startHardClip + endHardClip;

         // Fix the basecalls
         final byte[] bases = new byte[len];
         Arrays.fill(bases, (byte) 'N');
         System.arraycopy(rec.getReadBases(), 0, bases, startHardClip, rec.getReadBases().length);

         // Fix the quality scores
         final byte[] quals = new byte[len];
         Arrays.fill(quals, (byte) 2  );
         System.arraycopy(rec.getBaseQualities(), 0, quals, startHardClip, rec.getBaseQualities().length);

         // Fix the cigar!
         elements = new ArrayList<CigarElement>(elements); // make it modifiable
         if (startHardClip > 0) elements.set(0, new CigarElement(first.getLength(), CigarOperator.S));
         if (endHardClip   > 0) elements.set(elements.size()-1, new CigarElement(last.getLength(), CigarOperator.S));

         // Set the update structures on the new record
         rec.setReadBases(bases);
         rec.setBaseQualities(quals);
         rec.setCigar(new Cigar(elements));
     }
}
 

private SAMRecord makeRead(final SAMFileHeader alignedHeader, final SAMRecord unmappedRec, final HitSpec hitSpec, final int hitIndex) {
    final SAMRecord rec = new SAMRecord(alignedHeader);
    rec.setReadName(unmappedRec.getReadName());
    rec.setReadBases(unmappedRec.getReadBases());
    rec.setBaseQualities(unmappedRec.getBaseQualities());
    rec.setMappingQuality(hitSpec.mapq);
    if (!hitSpec.primary) rec.setNotPrimaryAlignmentFlag(true);
    final Cigar cigar = new Cigar();
    final int readLength = rec.getReadLength();
    if (hitSpec.filtered) {
        // Add two insertions so alignment is filtered.
        cigar.add(new CigarElement(readLength-4, CigarOperator.M));
        cigar.add(new CigarElement(1, CigarOperator.I));
        cigar.add(new CigarElement(1, CigarOperator.M));
        cigar.add(new CigarElement(1, CigarOperator.I));
        cigar.add(new CigarElement(1, CigarOperator.M));
    } else {
        cigar.add(new CigarElement(readLength, CigarOperator.M));
    }
    rec.setCigar(cigar);

    rec.setReferenceName(bigSequenceName);
    rec.setAttribute(SAMTag.HI.name(), hitIndex);
    rec.setAlignmentStart(hitIndex + 1);

    return rec;
}
 

/**
 * Write out a representation of this haplotype as a read
 *
 * @param haplotype a haplotype to write out, must not be null
 * @param paddedRefLoc the reference location, must not be null
 * @param isAmongBestHaplotypes true if among the best haplotypes, false if it was just one possible haplotype
 * @param callableRegion the region over which variants are being called
 */
private void writeHaplotype(final Haplotype haplotype,
                            final Locatable paddedRefLoc,
                            final boolean isAmongBestHaplotypes,
                            final Locatable callableRegion) {
    Utils.nonNull(haplotype, "haplotype cannot be null");
    Utils.nonNull(paddedRefLoc, "paddedRefLoc cannot be null");

    final SAMRecord record = new SAMRecord(output.getBAMOutputHeader());
    record.setReadBases(haplotype.getBases());
    record.setAlignmentStart(paddedRefLoc.getStart() + haplotype.getAlignmentStartHapwrtRef());
    // Use a base quality value "!" for it's display value (quality value is not meaningful)
    record.setBaseQualities(Utils.dupBytes((byte) '!', haplotype.getBases().length));
    record.setCigar(haplotype.getCigar());
    record.setMappingQuality(isAmongBestHaplotypes ? bestHaplotypeMQ : otherMQ);
    record.setReadName(output.getHaplotypeSampleTag() + uniqueNameCounter++);
    record.setAttribute(output.getHaplotypeSampleTag(), haplotype.hashCode());
    record.setReadUnmappedFlag(false);
    record.setReferenceIndex(output.getBAMOutputHeader().getSequenceIndex(paddedRefLoc.getContig()));
    record.setAttribute(SAMTag.RG.toString(), output.getHaplotypeReadGroupID());
    record.setFlags(SAMFlag.READ_REVERSE_STRAND.intValue());
    if (callableRegion != null) {
        record.setAttribute(AssemblyBasedCallerUtils.CALLABLE_REGION_TAG, callableRegion.toString());
    }

    output.add(new SAMRecordToGATKReadAdapter(record));
}
 

@Test
public void test() {
	SAMFileHeader header = new SAMFileHeader();
	SAMRecord r1 = new SAMRecord(header);
	r1.setReadName("readName");
	r1.setFlags(4);
	r1.setReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
	r1.setAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
	r1.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);
	r1.setCigar(new Cigar());
	r1.setMateReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
	r1.setMateAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
	r1.setReadBases("A".getBytes());
	r1.setBaseQualityString("!");

	BAMRecordView view = new BAMRecordView(new byte[1024]);
	translate(r1, view);
	r1.setReadName("2");
	translate(r1, view);

	List<SAMRecord> list = toSAMRecord(view, header);
	assertEquals(2, list.size());

	Iterator<SAMRecord> iterator = list.iterator();
	SAMRecord r2 = iterator.next();
	r1.setReadName("readName");
	compare(r1, r2);

	r1.setReadName("2");
	r2 = iterator.next();
	compare(r1, r2);
}
 

public static void changeReadLength(SAMRecord record, int newLength) {
	if (newLength == record.getReadLength())
		return;
	if (newLength < 1 || newLength >= record.getReadLength())
		throw new IllegalArgumentException("Cannot change read length to " + newLength);

	List<CigarElement> newCigarElements = new ArrayList<CigarElement>();
	int len = 0;
	for (CigarElement ce : record.getCigar().getCigarElements()) {
		switch (ce.getOperator()) {
		case D:
			break;
		case S:
			// dump = true;
			// len -= ce.getLength();
			// break;
		case M:
		case I:
		case X:
			len += ce.getLength();
			break;

		default:
			throw new IllegalArgumentException("Unexpected cigar operator: " + ce.getOperator() + " in cigar "
					+ record.getCigarString());
		}

		if (len <= newLength) {
			newCigarElements.add(ce);
			continue;
		}
		CigarElement newCe = new CigarElement(ce.getLength() - (record.getReadLength() - newLength),
				ce.getOperator());
		if (newCe.getLength() > 0)
			newCigarElements.add(newCe);
		break;
	}

	byte[] newBases = new byte[newLength];
	System.arraycopy(record.getReadBases(), 0, newBases, 0, newLength);
	record.setReadBases(newBases);

	byte[] newScores = new byte[newLength];
	System.arraycopy(record.getBaseQualities(), 0, newScores, 0, newLength);

	record.setCigar(new Cigar(newCigarElements));
}
 

/**
 * Updates the CIGAR in the supplied record to be legacy-compatible.
 * @param record the <code>SAMRecord</code> to update.
 */
public static void convertToLegacyCigar(SAMRecord record) {
  record.setCigar(convertToLegacyCigar(record.getCigar()));
}