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

The following examples show how to use htsjdk.samtools.SAMRecord#getAttribute() . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
Source File: NumberEvents.java    From hmftools with GNU General Public License v3.0 6 votes vote down vote up
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;
}
 
Example 2
Source File: CompareBAMTagValues.java    From Drop-seq with MIT License 6 votes vote down vote up
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;
}
 
Example 3
Source File: BamTagHistogram.java    From Drop-seq with MIT License 6 votes vote down vote up
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);

}
 
Example 4
Source File: IntervalTagComparator.java    From Drop-seq with MIT License 5 votes vote down vote up
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);
    }
 
Example 5
Source File: CollectJumpingLibraryMetrics.java    From picard with MIT License 5 votes vote down vote up
/**
 * 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;
}
 
Example 6
Source File: LibraryIdGenerator.java    From picard with MIT License 5 votes vote down vote up
/**
 * 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;
}
 
Example 7
Source File: ClippingUtilityTest.java    From picard with MIT License 5 votes vote down vote up
@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);
            }
        }
    }
}
 
Example 8
Source File: SortedSAMWriter.java    From abra2 with MIT License 5 votes vote down vote up
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());
				}
			}
		}
	}
}
 
Example 9
Source File: AdapterDescriptor.java    From Drop-seq with MIT License 5 votes vote down vote up
@Override
public String getSequence(SAMRecord rec) {
    final String attribute = (String) rec.getAttribute(binaryTag);
    if (reverseComplement) {
        return SequenceUtil.reverseComplement(attribute);
    } else {
        return attribute;
    }
}
 
Example 10
Source File: FilterBam.java    From Drop-seq with MIT License 5 votes vote down vote up
/**
 * 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);
}
 
Example 11
Source File: TagValueFilteringIterator.java    From Drop-seq with MIT License 5 votes vote down vote up
@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;
 }
 
Example 12
Source File: RequiredTagPredicate.java    From Drop-seq with MIT License 5 votes vote down vote up
@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;
}
 
Example 13
Source File: FilterBamByTag.java    From Drop-seq with MIT License 5 votes vote down vote up
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;
}
 
Example 14
Source File: SAMRecordField.java    From cramtools with Apache License 2.0 5 votes vote down vote up
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());
	}
}
 
Example 15
Source File: AlignmentSummaryMetricsCollector.java    From picard with MIT License 4 votes vote down vote up
private boolean isNoiseRead(final SAMRecord record) {
    final Object noiseAttribute = record.getAttribute(ReservedTagConstants.XN);
    return (noiseAttribute != null && noiseAttribute.equals(1));
}
 
Example 16
Source File: BamToDetailsConversion.java    From chipster with MIT License 4 votes vote down vote up
/**
 * 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();
}
 
Example 17
Source File: AlignmentSummaryMetricsCollector.java    From picard with MIT License 4 votes vote down vote up
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;
                }
            }
        }
    }
}
 
Example 18
Source File: SamToFastq.java    From picard with MIT License 4 votes vote down vote up
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));
}
 
Example 19
Source File: BamTagOfTagCounts.java    From Drop-seq with MIT License 4 votes vote down vote up
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);
}
 
Example 20
Source File: LibraryIdGenerator.java    From gatk with BSD 3-Clause "New" or "Revised" License 2 votes vote down vote up
/**
 * 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);
}