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

public static int numberOfEvents(@NotNull final SAMRecord record) {
    Object nm = record.getAttribute("NM");
    if (!(nm instanceof Integer)) {
        return 0;
    }

    int additionalIndels = 0;
    int softClips = 0;
    for (CigarElement cigarElement : record.getCigar()) {
        switch (cigarElement.getOperator()) {
            case D:
            case I:
                additionalIndels += (cigarElement.getLength() - 1);
                break;
            case S:
                softClips++;
                break;
        }
    }

    return (Integer) nm - additionalIndels + softClips;
}
 

boolean compareTags (final SAMRecord r1, final SAMRecord r2, final List<String> tags) {
	boolean sameRead = r1.getReadName().equals(r2.getReadName());
	if (!sameRead)
		log.error("Read names differ at iteration.  R1: "+ r1.toString(), "R2: ", r2.toString());

	for (String tag: tags) {
		Object o1 = r1.getAttribute(tag);
		Object o2 = r2.getAttribute(tag);
		if (o1==null && o2==null) return true;
		if ((o1==null && o2!=null ) || (o1!=null && o2==null)) {
			log.error("Read tag values differ for tag: ["+ tag.toString()+ "] R1 is null [" + String.valueOf(o1==null) + "] R2 is null [" + String.valueOf(o2==null)+"]");
			return false;
		}


		if (!o1.equals(o2)) {
			log.error("Read tag values differ for tag: ["+ tag.toString()+ "] "  + r1.toString()+ "[" +o1.toString()+ "] "+ r2.toString() +" [" + o2.toString()+"]");
			return false;
		}

	}
	return true;
}
 

public String getAnyTagAsString (final SAMRecord r, final String tag) {
	String s = null;
	Object o = r.getAttribute(tag);
	if (o==null) return (s);
	if (o instanceof String) {
		s = (String) o;
		return (s);
	}
	if (o instanceof Integer) {
		Integer i = (Integer) o;
		s= i.toString();
		return (s);
	}
	return (s);

}
 

private String [] getFirstSplit(final SAMRecord rec) {

    	Object strIntervalRec = rec.getAttribute(tag);
    	if (!(strIntervalRec instanceof String))
			throw new IllegalArgumentException(SAMTagUtil.getSingleton().makeStringTag(this.tag) + " does not have a String value");
    	String intervalString = (String) strIntervalRec;
    	String [] result = new String [4];
    	StringUtil.splitConcatenateExcessTokens(intervalString, result, ENCODE_DELIMITER);
    	return (result);
    }
 

/**
 * Calculates the mode for outward-facing pairs, using the first SAMPLE_FOR_MODE
 * outward-facing pairs found in INPUT
 */
private double getOutieMode() {

    int samplePerFile = SAMPLE_FOR_MODE / INPUT.size();

    Histogram<Integer> histo = new Histogram<Integer>();

    for (File f : INPUT) {
        SamReader reader = SamReaderFactory.makeDefault().open(f);
        int sampled = 0;
        for (Iterator<SAMRecord> it = reader.iterator(); it.hasNext() && sampled < samplePerFile; ) {
            SAMRecord sam = it.next();
            if (!sam.getFirstOfPairFlag()) {
                continue;
            }
            // If we get here we've hit the end of the aligned reads
            if (sam.getReadUnmappedFlag() && sam.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
                break;
            } else if (sam.getReadUnmappedFlag() || sam.getMateUnmappedFlag()) {
                continue;
            } else if ((sam.getAttribute(SAMTag.MQ.name()) == null ||
                    sam.getIntegerAttribute(SAMTag.MQ.name()) >= MINIMUM_MAPPING_QUALITY) &&
                    sam.getMappingQuality() >= MINIMUM_MAPPING_QUALITY &&
                    sam.getMateNegativeStrandFlag() != sam.getReadNegativeStrandFlag() &&
                    sam.getMateReferenceIndex().equals(sam.getReferenceIndex()) &&
                    SamPairUtil.getPairOrientation(sam) == PairOrientation.RF) {
                histo.increment(Math.abs(sam.getInferredInsertSize()));
                sampled++;
            }
        }
        CloserUtil.close(reader);
    }

    return histo.size() > 0 ? histo.getMode() : 0;
}
 

/**
 * Gets the library name from the header for the record. If the RG tag is not present on
 * the record, or the library isn't denoted on the read group, a constant string is
 * returned.
 */
public static String getLibraryName(final SAMFileHeader header, final SAMRecord rec) {
    final String readGroupId = (String) rec.getAttribute(ReservedTagConstants.READ_GROUP_ID);

    if (readGroupId != null) {
        final SAMReadGroupRecord rg = header.getReadGroup(readGroupId);
        if (rg != null) {
            final String libraryName = rg.getLibrary();
            if (null != libraryName) return libraryName;
        }
    }

    return UNKNOWN_LIBRARY;
}
 

@Test(dataProvider = "testAdapterInAllReadPositionsDataProvider")
public void testAdapterInAllReadPositions(final int readLength) {
    final int minAdapterLength = 6;
    for (final IlluminaAdapterPair adapterPair : IlluminaAdapterPair.values()) {
        final AdapterMarker marker = new AdapterMarker(adapterPair);
        for (int adapterPosition = 0; adapterPosition < readLength; ++adapterPosition) {
            final SAMRecord rec = createSamRecordWithAdapterSequence(readLength, adapterPair, adapterPosition);
            final AdapterPair matchedAdapter = ClippingUtility.adapterTrimIlluminaSingleRead(rec, minAdapterLength, 0.1, adapterPair);
            final Object xt = rec.getAttribute(ReservedTagConstants.XT);

            rec.setAttribute(ReservedTagConstants.XT, null);
            final AdapterPair truncatedAdapter = marker.adapterTrimIlluminaSingleRead(rec, minAdapterLength, 0.1);
            final Object xtFromMarker = rec.getAttribute(ReservedTagConstants.XT);

            if (adapterPosition <= readLength - minAdapterLength) {
                Assert.assertEquals(matchedAdapter, adapterPair, "Failed at position " + adapterPosition);
                Assert.assertEquals((Integer)xt, Integer.valueOf(adapterPosition + 1));
                Assert.assertNotNull(truncatedAdapter);
                Assert.assertEquals((Integer)xtFromMarker, Integer.valueOf(adapterPosition + 1));
            } else {
                Assert.assertNull(matchedAdapter);
                Assert.assertNull(xt);
                Assert.assertNull(truncatedAdapter);
                Assert.assertNull(xtFromMarker);
            }
        }
    }
}
 

private void setMateInfo(SAMRecord read, Map<MateKey, SAMRecord> mates) {
	
	if (read.getReadPairedFlag()) {
		SAMRecord mate = mates.get(getMateKey(read));
		if (mate != null) {
			// Only update mate info if a read has been modified
			if (read.getAttribute("YO") != null || mate.getAttribute("YO") != null) {
				read.setMateAlignmentStart(mate.getAlignmentStart());
				read.setMateUnmappedFlag(mate.getReadUnmappedFlag());
				read.setMateNegativeStrandFlag(mate.getReadNegativeStrandFlag());
				 
				int start = read.getAlignmentStart() < mate.getAlignmentStart() ? read.getAlignmentStart() : mate.getAlignmentStart();
				int stop  = read.getAlignmentEnd() > mate.getAlignmentEnd() ? read.getAlignmentEnd() : mate.getAlignmentEnd();
				
				int insert = stop-start+1;
				
				if (read.getAlignmentStart() > mate.getAlignmentStart()) {
					insert *= -1;
				} else if (read.getAlignmentStart() == mate.getAlignmentStart() && mate.getFirstOfPairFlag()) {
					insert *= -1;
				}
				
				read.setInferredInsertSize(insert);
				
				if (read.getStringAttribute("MC") != null) {
					read.setAttribute("MC", mate.getCigarString());
				}
				
				if (!mate.getReadUnmappedFlag()) {
					read.setMateReferenceName(mate.getReferenceName());
				}
			}
		}
	}
}
 

@Override
public String getSequence(SAMRecord rec) {
    final String attribute = (String) rec.getAttribute(binaryTag);
    if (reverseComplement) {
        return SequenceUtil.reverseComplement(attribute);
    } else {
        return attribute;
    }
}
 

/**
 * Reject on any of the tags if they are set.
 */
boolean rejectOnTags (final List<String> tags, final SAMRecord r) {
	if (tags==null || tags.isEmpty()) return(false);

	for (String tag: tags) {
		Object v = r.getAttribute(tag);
		// if not null and union, return true if any are true.
		if (v!=null && this.TAG_REJECT_COMBINE_FLAG!=null && this.TAG_REJECT_COMBINE_FLAG.equals(UNION)) return true;
		// if null and intersect, return false
		if (v==null && this.TAG_REJECT_COMBINE_FLAG!=null && this.TAG_REJECT_COMBINE_FLAG.equals(INTERSECT)) return false;
		// if not null and single iteration, return true
		if (v!=null && this.TAG_REJECT_COMBINE_FLAG==null) return true;
	}
	return (false);
}
 

@Override
 public boolean filterOut(final SAMRecord rec) {
 	Object value = rec.getAttribute(requiredTag);
 	if (value == null)
return true;

 	if (this.expectedValues.contains(value)) return false;
     return true;
 }
 

@Override
public boolean test(SAMRecord rec) {
    for (final short tag : requiredTags)
        // String strTag=SAMTagUtil.getSingleton().makeStringTag(tag);
        // List<SAMTagAndValue> vals = rec.getAttributes();
        if (rec.getAttribute(tag) == null)
            return false;
    return true;
}
 

boolean filterRead(final SAMRecord r, final String tag, final Set<String> values,
		final boolean acceptFlag, final Integer mapQuality, boolean allowPartial) {

	// quickly filter on map quality if provided.
	if (mapQuality!=null && r.getMappingQuality() < mapQuality) return true;
	
	Object v = r.getAttribute(tag);
	// if the tag is not set, and you need a TAG set, then filter the read.
	if (v == null && acceptFlag)
		return true;
	// if the tag is not set, and you don't need a TAG set, then keep the
	// read.
	if (v == null && !acceptFlag)
		return false;

	final String vv =  getTagValue(v);

	// if there are no values to scan, it's a match. Start with that.
	boolean hasElement = true;
	
	// if there are values, check to see if this tag matches one.
	if (values != null && values.size() > 0)
		if (!allowPartial)
			hasElement = values.contains(vv);
		else {
			hasElement = values.stream().anyMatch((x -> vv.contains(x)));
		}

	if ((hasElement & acceptFlag)
			| (hasElement == false & acceptFlag == false))
		return false;
	if ((hasElement == false & acceptFlag)
			| (hasElement & acceptFlag == false))
		return true;

	return false;
}
 

public Object getValue(SAMRecord record) {
	switch (type) {
	case QNAME:
		return record.getReadName();
	case FLAG:
		return Integer.toString(record.getFlags());
	case RNAME:
		return record.getReferenceName();
	case POS:
		return Integer.toString(record.getAlignmentStart());
	case MAPQ:
		return Integer.toString(record.getMappingQuality());
	case CIGAR:
		return record.getCigarString();
	case RNEXT:
		return record.getMateReferenceName();
	case PNEXT:
		return Integer.toString(record.getMateAlignmentStart());
	case TLEN:
		return Integer.toString(record.getInferredInsertSize());
	case SEQ:
		return record.getReadString();
	case QUAL:
		return record.getBaseQualityString();

	case TAG:
		if (tagId == null)
			throw new IllegalArgumentException("Tag mismatch reqiues tag id. ");
		return record.getAttribute(tagId);

	default:
		throw new IllegalArgumentException("Unknown record field: " + type.name());
	}
}
 

private boolean isNoiseRead(final SAMRecord record) {
    final Object noiseAttribute = record.getAttribute(ReservedTagConstants.XN);
    return (noiseAttribute != null && noiseAttribute.equals(1));
}
 

/**
 * Find reads in a given range.
 * 
 * <p>
 * TODO add cigar to the list of returned values
 * <p>
 * TODO add pair information to the list of returned values
 * 
 * @param request
 * @return
 * @throws InterruptedException 
 */
@Override
protected void processDataRequest(DataRequest request) throws InterruptedException {
	
	// Read the given region
	CloseableIterator<SAMRecord> iterator = dataSource.query(request.start.chr, request.start.bp.intValue(), request.end.bp.intValue());
	
	// Produce results
	while (iterator.hasNext()) {

		List<Feature> responseList = new LinkedList<Feature>();

		// Split results into chunks
		for (int c = 0; c < RESULT_CHUNK_SIZE && iterator.hasNext(); c++) {
			SAMRecord record = iterator.next();

			// Region for this read
			Region recordRegion = new Region((long) record.getAlignmentStart(), (long) record.getAlignmentEnd(), request.start.chr);

			// Values for this read
			LinkedHashMap<DataType, Object> values = new LinkedHashMap<DataType, Object>();

			Feature read = new Feature(recordRegion, values);

			if (request.getRequestedContents().contains(DataType.ID)) {
				values.put(DataType.ID, record.getReadName());
			}

			if (request.getRequestedContents().contains(DataType.STRAND)) {
				values.put(DataType.STRAND, BamUtils.getStrand(record, coverageType));					
			}

			if (request.getRequestedContents().contains(DataType.QUALITY)) {
				values.put(DataType.QUALITY, record.getBaseQualityString());
			}

			if (request.getRequestedContents().contains(DataType.CIGAR)) {
				Cigar cigar = new Cigar(read, record.getCigar());
				values.put(DataType.CIGAR, cigar);
			}

			// TODO Deal with "=" and "N" in read string
			if (request.getRequestedContents().contains(DataType.SEQUENCE)) {
				String seq = record.getReadString();
				values.put(DataType.SEQUENCE, seq);
			}

			if (request.getRequestedContents().contains(DataType.MATE_POSITION)) {
				
				BpCoord mate = new BpCoord((Long)(long)record.getMateAlignmentStart(),
						new Chromosome(record.getMateReferenceName()));
				
				values.put(DataType.MATE_POSITION, mate);
			}
			
			if (request.getRequestedContents().contains(DataType.BAM_TAG_NH)) {
				Object ng = record.getAttribute("NH");
				if (ng != null) {
					values.put(DataType.BAM_TAG_NH, (Integer)record.getAttribute("NH"));
				}
			}
			
			/*
			 * NOTE! RegionContents created from the same read area has to be equal in methods equals, hash and compareTo. Primary types
			 * should be ok, but objects (including tables) has to be handled in those methods separately. Otherwise tracks keep adding
			 * the same reads to their read sets again and again.
			 */
			responseList.add(read);
		}

		// Send result			
		super.createDataResult(new DataResult(request.getStatus(), responseList));			
	}

	// We are done
	iterator.close();
}
 

private void collectReadData(final SAMRecord record) {
    // NB: for read count metrics, do not include supplementary records, but for base count metrics, do include supplementary records.
    if (record.getSupplementaryAlignmentFlag()) return;

    metrics.TOTAL_READS++;
    readLengthHistogram.increment(record.getReadBases().length);

    if (!record.getReadFailsVendorQualityCheckFlag()) {
        metrics.PF_READS++;
        if (isNoiseRead(record)) metrics.PF_NOISE_READS++;

        // See if the read is an adapter sequence
        if (adapterUtility.isAdapter(record)) {
            this.adapterReads++;
        }

        if (!record.getReadUnmappedFlag() && doRefMetrics) {
            metrics.PF_READS_ALIGNED++;
            if (record.getReadPairedFlag() && !record.getProperPairFlag()) metrics.PF_READS_IMPROPER_PAIRS++;
            if (!record.getReadNegativeStrandFlag()) numPositiveStrand++;
            if (record.getReadPairedFlag() && !record.getMateUnmappedFlag()) {
                metrics.READS_ALIGNED_IN_PAIRS++;

                // Check that both ends have mapq > minimum
                final Integer mateMq = record.getIntegerAttribute(SAMTag.MQ.toString());
                if (mateMq == null || mateMq >= MAPPING_QUALITY_THRESHOLD && record.getMappingQuality() >= MAPPING_QUALITY_THRESHOLD) {
                    ++this.chimerasDenominator;

                    // With both reads mapped we can see if this pair is chimeric
                    if (ChimeraUtil.isChimeric(record, maxInsertSize, expectedOrientations)) {
                        ++this.chimeras;
                    }
                }
            } else { // fragment reads or read pairs with one end that maps
                // Consider chimeras that occur *within* the read using the SA tag
                if (record.getMappingQuality() >= MAPPING_QUALITY_THRESHOLD) {
                    ++this.chimerasDenominator;
                    if (record.getAttribute(SAMTag.SA.toString()) != null) ++this.chimeras;
                }
            }
        }
    }
}
 

private void writeRecord(final SAMRecord read, final Integer mateNumber, final FastqWriter writer,
                         final int basesToTrim, final Integer maxBasesToWrite) {
    final String seqHeader = mateNumber == null ? read.getReadName() : read.getReadName() + "/" + mateNumber;
    String readString = read.getReadString();
    String baseQualities = read.getBaseQualityString();

    // If we're clipping, do the right thing to the bases or qualities
    if (CLIPPING_ATTRIBUTE != null) {
        Integer clipPoint = (Integer) read.getAttribute(CLIPPING_ATTRIBUTE);
        if (clipPoint != null && clipPoint < CLIPPING_MIN_LENGTH) {
            clipPoint = Math.min(readString.length(), CLIPPING_MIN_LENGTH);
        }

        if (clipPoint != null) {
            if (CLIPPING_ACTION.equalsIgnoreCase(CLIP_TRIM)) {
                readString = clip(readString, clipPoint, null, !read.getReadNegativeStrandFlag());
                baseQualities = clip(baseQualities, clipPoint, null, !read.getReadNegativeStrandFlag());
            } else if (CLIPPING_ACTION.equalsIgnoreCase(CLIP_TO_N)) {
                readString = clip(readString, clipPoint, CLIP_TO_N.charAt(0), !read.getReadNegativeStrandFlag());
            } else {
                final char newQual = SAMUtils.phredToFastq(new byte[]{(byte) Integer.parseInt(CLIPPING_ACTION)}).charAt(0);
                baseQualities = clip(baseQualities, clipPoint, newQual, !read.getReadNegativeStrandFlag());
            }
        }
    }

    if (RE_REVERSE && read.getReadNegativeStrandFlag()) {
        readString = SequenceUtil.reverseComplement(readString);
        baseQualities = StringUtil.reverseString(baseQualities);
    }

    if (basesToTrim > 0) {
        readString = readString.substring(basesToTrim);
        baseQualities = baseQualities.substring(basesToTrim);
    }

    // Perform quality trimming if desired, making sure to leave at least one base!
    if (QUALITY != null) {
        final byte[] quals = SAMUtils.fastqToPhred(baseQualities);
        final int qualityTrimIndex = Math.max(1, TrimmingUtil.findQualityTrimPoint(quals, QUALITY));
        if (qualityTrimIndex < quals.length) {
            readString = readString.substring(0, qualityTrimIndex);
            baseQualities = baseQualities.substring(0, qualityTrimIndex);
        }
    }

    if (maxBasesToWrite != null && maxBasesToWrite < readString.length()) {
        readString = readString.substring(0, maxBasesToWrite);
        baseQualities = baseQualities.substring(0, maxBasesToWrite);
    }

    writer.write(new FastqRecord(seqHeader, readString, "", baseQualities));
}
 

public TagOfTagResults<String,String> getResults (final File inputBAM, final String primaryTag, final String secondaryTag, final boolean filterPCRDuplicates, final Integer readQuality) {
	TagOfTagResults<String,String> result = new TagOfTagResults<>();
	SamReader reader = SamReaderFactory.makeDefault().open(inputBAM);
	CloseableIterator<SAMRecord> iter = CustomBAMIterators.getReadsInTagOrder(reader, primaryTag);
	CloserUtil.close(reader);
	String currentTag="";
	Set<String> otherTagCollection=new HashSet<>();

	ProgressLogger progress = new ProgressLogger(log);

	while (iter.hasNext()) {
		SAMRecord r = iter.next();
		progress.record(r);
		// skip reads that don't pass filters.
		boolean discardResult=(filterPCRDuplicates && r.getDuplicateReadFlag()) || r.getMappingQuality()<readQuality || r.isSecondaryOrSupplementary();
		// short circuit if read is below the quality to be considered.
		if (discardResult) continue;

		Object d = r.getAttribute(secondaryTag);
		// short circuit if there's no tag for this read.
		if (d==null) continue;

		String data=null;

		if (d instanceof String)
			data=r.getStringAttribute(secondaryTag);
		else if (d instanceof Integer)
			data=Integer.toString(r.getIntegerAttribute(secondaryTag));


		String newTag = r.getStringAttribute(primaryTag);
		if (newTag==null)
			newTag="";

		// if you see a new tag.
		if (!currentTag.equals(newTag)) {
			// write out tag results, if any.
			if (!currentTag.equals("") && otherTagCollection.size()>0)
				result.addEntries(currentTag, otherTagCollection);
			currentTag=newTag;
			otherTagCollection.clear();
			if (!discardResult) otherTagCollection=addTagToCollection (data,otherTagCollection);

		} else // gather stats
		if (!discardResult) otherTagCollection=addTagToCollection (data,otherTagCollection);
	}
	if (otherTagCollection.size()>0)
		result.addEntries(currentTag, otherTagCollection);

	return (result);
}
 

/**
 * Gets the library name from the header for the record. If the RG tag is not present on
 * the record, or the library isn't denoted on the read group, a constant string is
 * returned.
 */
public static String getLibraryName(final SAMFileHeader header, final SAMRecord rec) {
    final String readGroupId = (String) rec.getAttribute(ReservedTagConstants.READ_GROUP_ID);
    return getLibraryName(header, readGroupId);
}