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

public void addQualityScores(final SAMRecord s, final CramCompressionRecord r, final ReferenceTracks t) {
	if (s.getBaseQualities() == SAMRecord.NULL_QUALS) {
		r.qualityScores = SAMRecord.NULL_QUALS;
		r.setForcePreserveQualityScores(false);
		return;
	}

	final byte[] scores = new byte[s.getReadLength()];
	Arrays.fill(scores, (byte) -1);
	for (final PreservationPolicy p : policyList)
		addQS(s, r, scores, t, p);

	if (!r.isForcePreserveQualityScores()) {
		for (int i = 0; i < scores.length; i++) {
			if (scores[i] > -1) {
				if (r.readFeatures == null || r.readFeatures == Collections.EMPTY_LIST)
					r.readFeatures = new LinkedList<ReadFeature>();

				r.readFeatures.add(new BaseQualityScore(i + 1, scores[i]));
			}
		}
		if (r.readFeatures != null)
			Collections.sort(r.readFeatures, readFeaturePositionComparator);
	}
	r.qualityScores = scores;
}
 

/**
 * When an adapter is matched in only one end of a pair, we check it again with
 * stricter thresholds.  If it still matches, then we trim both ends of the read
 * at the same location.
  */
private static boolean attemptOneSidedMatch(final SAMRecord read1,
                                         final SAMRecord read2,
                                         final int index1,
                                         final int index2,
                                         final int stricterMinMatchBases) {

    // Save all the data about the read where we found the adapter match
    final int matchedIndex = index1 == NO_MATCH ? index2 : index1;
    final SAMRecord matchedRead = index1 == NO_MATCH ? read2 : read1;


    // If it still matches with a stricter minimum matched bases, then
    // clip both reads
    if (matchedRead.getReadLength() - matchedIndex >= stricterMinMatchBases) {
        if (read1.getReadBases().length > matchedIndex) {
            read1.setAttribute(ReservedTagConstants.XT, matchedIndex + 1);
        }
        if (read2.getReadBases().length > matchedIndex) {
            read2.setAttribute(ReservedTagConstants.XT, matchedIndex + 1);
        }
        return true;
    }
    return false;
}
 

/**
 *  Returns original edit distance as set in YX tag.
 */
public static int getOrigEditDistance(SAMRecord read) {
	
	Integer distance = null;
	
	if (read.getReadUnmappedFlag()) {
		distance = read.getReadLength();
	} else {
		distance = read.getIntegerAttribute("YX");
		
		if (distance == null) {
			distance = read.getReadLength();
		}
	}
			
	return distance;
}
 

private void addRead(final GcObject gcObj, final SAMRecord rec, final String group, final byte[] gc, final byte[] refBases) {
    if (!rec.getReadPairedFlag() || rec.getFirstOfPairFlag()) ++gcObj.totalClusters;
    final int pos = rec.getReadNegativeStrandFlag() ? rec.getAlignmentEnd() - scanWindowSize : rec.getAlignmentStart();
    ++gcObj.totalAlignedReads;
    if (pos > 0) {
        final int windowGc = gc[pos];
        if (windowGc >= 0) {
            ++gcObj.readsByGc[windowGc];
            gcObj.basesByGc[windowGc] += rec.getReadLength();
            gcObj.errorsByGc[windowGc] +=
                    SequenceUtil.countMismatches(rec, refBases, bisulfite) +
                            SequenceUtil.countInsertedBases(rec) + SequenceUtil.countDeletedBases(rec);
        }
    }
    if (gcObj.group == null) {
        gcObj.group = group;
    }
}
 

private boolean shouldAssemble(SAMRecord read) {
	return (
		(!read.getDuplicateReadFlag()) && 
		(!read.getReadFailsVendorQualityCheckFlag()) &&
		read.getReadLength() > 0 &&
		(read.getMappingQuality() >= this.minMapqForAssembly || read.getReadUnmappedFlag()) &&
		SAMRecordUtils.isPrimary(read));  // Was previously an id check, so supplemental / secondary alignments could be included
}
 

public static int getUnclippedLength(SAMRecord read) {
	int softClipLen = 0;
	for (CigarElement elem : read.getCigar().getCigarElements()) {
		if (elem.getOperator() == CigarOperator.S) {
			softClipLen += elem.getLength();
		}
	}
	
	return read.getReadLength() - softClipLen;
}
 

public static int getMappedLength(SAMRecord read) {
	int length = read.getReadLength();
	for (CigarElement elem : read.getCigar().getCigarElements()) {
		if (elem.getOperator() == CigarOperator.S) {
			length -= elem.getLength();
		}
	}
	
	return length;
}
 

/**
 * Method returns reverse complemented sequence for the part of the record. Can work with 3' and 5' ends
 * (if start index < 0, then it will found the index in the end of sequence by adding the length of record).
 * @param record read from SAM file to process
 * @param startIndex index where start the sequence
 * @param length length of pert of sequence
 * @return reverse complemented part of record
 */
public static String getReverseComplementedSequence(SAMRecord record, int startIndex, int length) {
    if (startIndex < 0) {
        startIndex = record.getReadLength() + startIndex;
    }
    byte[] rangeBytes = Arrays.copyOfRange(record.getReadBases(), startIndex, startIndex + length);
    SequenceUtil.reverseComplement(rangeBytes);
    return new String(rangeBytes);
}
 

/**
 * Returns 1-based index of first base in read that corresponds to M in CIGAR string.
 * Note that first is relative to 5' end, so that for reverse-strand alignment, the index of
 * the last base aligned is computed relative to the end of the read.
 */
int getIndexOfFirstAlignedBase(final SAMRecord rec) {
    final List<AlignmentBlock> alignmentBlocks = rec.getAlignmentBlocks();
    if (rec.getReadNegativeStrandFlag()) {
        final AlignmentBlock alignmentBlock = alignmentBlocks.get(alignmentBlocks.size() - 1);
        return rec.getReadLength() - CoordMath.getEnd(alignmentBlock.getReadStart(), alignmentBlock.getLength()) + 1;
    } else {
        return alignmentBlocks.get(0).getReadStart();
    }
}
 

/**
 * Returns 1-based index of first base in read that corresponds to M in CIGAR string.
 * Note that first is relative to 5' end, so that for reverse-strand alignment, the index of
 * the last base aligned is computed relative to the end of the read.
 */
int getIndexOfFirstAlignedBase(final SAMRecord rec) {
    final List<AlignmentBlock> alignmentBlocks = rec.getAlignmentBlocks();
    if (rec.getReadNegativeStrandFlag()) {
        final AlignmentBlock alignmentBlock = alignmentBlocks.get(alignmentBlocks.size() - 1);
        return rec.getReadLength() - CoordMath.getEnd(alignmentBlock.getReadStart(), alignmentBlock.getLength()) + 1;
    } else {
        return alignmentBlocks.get(0).getReadStart();
    }
}
 

private Character getBaseAtPosition(SAMRecord read, int refPos) {
	int readPos = 0;
	int refPosInRead = read.getAlignmentStart();
	int cigarElementIdx = 0;
	
	while (refPosInRead <= refPos && cigarElementIdx < read.getCigar().numCigarElements() && readPos < read.getReadLength()) {
		CigarElement elem = read.getCigar().getCigarElement(cigarElementIdx++);
		
		switch(elem.getOperator()) {
			case H: //NOOP
				break;
			case S:
			case I:
				readPos += elem.getLength();
				break;
			case D:
			case N:
				refPosInRead += elem.getLength();
				break;
			case M:
				if (refPos < (refPosInRead + elem.getLength())) {
					readPos += refPos - refPosInRead;
					if (readPos < read.getReadLength()) {
						// Found the base.  Return it
						return read.getReadString().charAt(readPos);
					}
				} else {
					readPos += elem.getLength();
					refPosInRead += elem.getLength();
				}
				break;
			default:
				throw new IllegalArgumentException("Invalid Cigar Operator: " + elem.getOperator() + " for read: " + read.getSAMString());					
		}
	}
	
	return null;
}
 

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;
}
 

public static int basesDeletedAfterPosition(@NotNull final SAMRecord record, int position) {
    int startReadPosition = record.getReadPositionAtReferencePosition(position);
    assert (startReadPosition != 0);

    int nextReferencePosition = record.getReferencePositionAtReadPosition(startReadPosition + 1);
    return nextReferencePosition == 0 && startReadPosition == record.getReadLength()
            ? record.getAlignmentEnd() - position
            : Math.max(0, nextReferencePosition - position - 1);
}
 

ModifiableVariantHotspotEvidence findEvidenceOfDelete(@NotNull final ModifiableVariantHotspotEvidence builder,
        @NotNull final VariantHotspot hotspot, @NotNull final SAMRecord record) {
    assert (hotspot.isSimpleDelete());

    int hotspotStartPosition = (int) hotspot.position();
    int recordStartPosition = record.getReadPositionAtReferencePosition(hotspotStartPosition);
    if (recordStartPosition == 0) {
        return builder;
    }

    int recordStartQuality = getBaseQuality(record, recordStartPosition);

    if (containsDelete(record, hotspotStartPosition, hotspot.ref())) {
        int quality =
                record.getReadLength() > recordStartPosition ? getAvgBaseQuality(record, recordStartPosition, 2) : recordStartQuality;
        if (quality < minBaseQuality) {
            return builder;
        }
        return builder.setReadDepth(builder.readDepth() + 1)
                .setIndelSupport(builder.indelSupport() + 1)
                .setAltQuality(builder.altQuality() + quality)
                .setAltSupport(builder.altSupport() + 1);
    }

    if (recordStartQuality < minBaseQuality) {
        return builder;
    }

    int insertedBases = basesInsertedAfterPosition(record, hotspotStartPosition);
    int deletedBases = basesDeletedAfterPosition(record, hotspotStartPosition);
    if (insertedBases == 0 && deletedBases == 0 && record.getReadString().charAt(recordStartPosition - 1) == hotspot.ref().charAt(0)) {
        builder.setRefSupport(builder.refSupport() + 1);
    } else if (insertedBases > 0 || deletedBases > 0) {
        builder.setIndelSupport(builder.indelSupport() + 1);
    }

    return builder.setReadDepth(builder.readDepth() + 1);
}
 

private void populateNext(boolean force) throws IOException {
  final SAMRecord previousRecord = mNextRecord;
  mNextRecord = null;
  if (force) {
    advanceSubIterator();
  }
  while (mCurrentIt != null) {
    if (!mCurrentIt.hasNext()) {  // Only happens when stream is exhausted, so effectively just closes things out.
      advanceSubIterator();
    } else {
      final SAMRecord rec = nextOrBuffered();
      final int refId = rec.getReferenceIndex();

      if (refId > mCurrentTemplate) { // current offset has exceeded region and block overlapped next template
        recordPreviousAndAdvance(previousRecord, rec);
      } else {

        if (refId < mCurrentTemplate) { // Current block may occasionally return records from the previous template if the block overlaps
          continue;
        }

        final int alignmentStart = rec.getAlignmentStart() - 1; // to 0-based
        int alignmentEnd = rec.getAlignmentEnd(); // to 0-based exclusive = noop
        if (alignmentEnd == 0) {
          // Use the read length to get an end point for unmapped reads
          alignmentEnd = rec.getAlignmentStart() + rec.getReadLength();
        }

        if (alignmentEnd <= mCurrentRegion.getStart()) {  // before region
          continue;
        }

        if (alignmentStart <= mPreviousAlignmentStart) { // this record would have been already returned by an earlier region
          //Diagnostic.developerLog("Ignoring record from earlier block at " + rec.getReferenceName() + ":" + rec.getAlignmentStart());
          ++mDoubleFetched;
          if (mDoubleFetched % 100000 == 0) {
            Diagnostic.developerLog("Many double-fetched records for source " + mLabel + " noticed at " + rec.getReferenceName() + ":" + rec.getAlignmentStart() + " in region " + mCurrentRegion + " (skipping through to " + mPreviousAlignmentStart + ")");
          }
          continue;
        }

        if (alignmentStart >= mCurrentRegion.getEnd()) {  // past current region, advance the iterator and record the furtherest we got
          recordPreviousAndAdvance(previousRecord, rec);
          continue;
        }

        mNextRecord = rec;
        break;
      }
    }
  }
}
 

private IndelInfo checkForIndelAtLocus(SAMRecord read, int refPos) {
		IndelInfo elem = null;
		
//		if (refPos == 105243047 && read.getReadName().equals("D7T4KXP1:400:C5F94ACXX:5:2302:20513:30410")) {
//			System.out.println("bar");
//		}
		
		String contigInfo = read.getStringAttribute("YA");
		if (contigInfo != null) {
			// Get assembled contig info.
			String[] fields = contigInfo.split(":");
			int contigPos = Integer.parseInt(fields[1]);
			
			Cigar contigCigar = TextCigarCodec.decode(fields[2]);
			
			// Check to see if contig contains indel at current locus
			elem = checkForIndelAtLocus(contigPos, contigCigar, refPos);
			
			if (elem != null) {
				// Now check to see if this read supports the indel
				IndelInfo readElem = checkForIndelAtLocus(read.getAlignmentStart(),
						read.getCigar(), refPos);
				
				// Allow partially overlapping indels to support contig
				// (Should only matter for inserts)
				if (readElem == null || readElem.getCigarElement().getOperator() != elem.getCigarElement().getOperator()) {
					// Read element doesn't match contig indel
					elem = null;
				} else {
					elem.setReadIndex(readElem.getReadIndex());
					
					// If this read overlaps the entire insert, capture the bases.
					if (elem.getCigarElement().getOperator() == CigarOperator.I) {

						if (elem.getCigarElement().getLength() == readElem.getCigarElement().getLength()) {
					
							String insertBases = read.getReadString().substring(readElem.getReadIndex(), readElem.getReadIndex()+readElem.getCigarElement().getLength());
							elem.setInsertBases(insertBases);
						} else if (readElem.getCigarElement().getLength() < elem.getCigarElement().getLength()) {
							
							int lengthDiff = elem.getCigarElement().getLength() - readElem.getCigarElement().getLength();
							
							if (readElem.getReadIndex() == 0) {
								elem.setReadIndex(readElem.getReadIndex() - lengthDiff);
							} else if (readElem.getReadIndex() == read.getReadLength()-1) {
								elem.setReadIndex(readElem.getReadIndex() + lengthDiff);
							}
						}
					}
				}
			}
		}
		
		return elem;
	}
 

private boolean isAssemblyTriggerCandidate(SAMRecord read, CompareToReference2 c2r, Feature region) {
		
		// High quality unmapped read anchored by mate
		if (!isSkipUnmappedTrigger && read.getReadUnmappedFlag() && !read.getMateUnmappedFlag() &&
			read.getReadLength() >= readLength * .9 && 
			SAMRecordUtils.getNumHighQualBases(read, MIN_CANDIDATE_BASE_QUALITY) >= readLength * .9) {
			return true;
		}
		
		int numGaps = countGaps(read.getCigar());
		
		int insertBases = getInsertBases(read);
		
		// More than one indel and/or splice in read
//		if (numGaps > 1) {
//			return true;
//		}
		
//		if (numGaps > 0) {
//			int qualAdjustedEditDist = c2r.numHighQualityMismatches(read, MIN_CANDIDATE_BASE_QUALITY) + SAMRecordUtils.getNumIndelBases(read);
//			if (qualAdjustedEditDist > (readLength * .25)) {
//				return true;
//			}			
//		}
		
		if (insertBases > readLength * .15) {
			return true;
		}
		
		if (maxSoftClipLength(read, region) > readLength * .25) {
			if (SAMRecordUtils.getNumHighQualBases(read, MIN_CANDIDATE_BASE_QUALITY) >= readLength * .9) {
				return true;
			}
		}
		
		// Increment candidate count for substantial high quality soft clipping
		// TODO: Higher base quality threshold?
//		if (read.getCigarString().contains("S") && (c2r.numHighQualityMismatches(read, MIN_CANDIDATE_BASE_QUALITY) > (readLength/10))) {
//			return true;
//		}


		
		// if indel is at least 10% of read length
		// TODO: Longer threshold (especially for deletions)?
//		if (numGaps == 1 && firstIndelLength(read.getCigar()) >= (readLength/10)) {
//			return true;
//		}
		
		// Read contains indel / intron and SNV (inclusive of soft clip mismatch)
		// TODO: Higher base qual threshold?
		// TODO: Minimum edit dist (1 base indel + 1 mismatch insufficient)
//		if (numGaps > 0 && c2r.numHighQualityMismatches(read, MIN_CANDIDATE_BASE_QUALITY) > 0) {
//			return true;
//		}
		
		// Increment candidate count for indels
//		if (read.getCigarString().contains("I") || read.getCigarString().contains("D") || read.getCigarString().contains("N")) {
//			isCandidate = true;
//		}
		
		// Increment candidate count if read contains at least 3 high quality mismatches
//		if (SAMRecordUtils.getIntAttribute(read, "NM") >= 3) {
//			if (c2r.numHighQualityMismatches(read, MIN_CANDIDATE_BASE_QUALITY) > 3) {
//				isCandidate = true;
//			}
//		}

		return false;
	}
 

private boolean findReadContext(int readIndex, @NotNull final SAMRecord record) {
    return readIndex >= config.readContextFlankSize() && readIndex < record.getReadLength() - config.readContextFlankSize();
}
 

/**
 * Filter a SAM record based on specified filtering parameters.
 * Does not consider position based filters, or inversion status.
 *
 * @param params parameters
 * @param rec record
 * @return true if record should be accepted for processing
 */
private static boolean checkRecordProperties(final SamFilterParams params, final SAMRecord rec) {
  final int flags = rec.getFlags();
  if ((flags & params.requireUnsetFlags()) != 0) {
    return false;
  }
  if ((flags & params.requireSetFlags()) != params.requireSetFlags()) {
    return false;
  }
  if (rec.getAlignmentStart() == 0 && params.excludeUnplaced()) {
    return false;
  }

  final boolean mated = (flags & SamBamConstants.SAM_READ_IS_MAPPED_IN_PROPER_PAIR) != 0;
  final boolean unmapped = rec.getReadUnmappedFlag();
  if (!mated && !unmapped && params.excludeUnmated()) {
    return false;
  }

  final int minMapQ = params.minMapQ();
  if (minMapQ >= 0 && rec.getMappingQuality() < minMapQ) {
    return false;
  }

  final Integer nh = SamUtils.getNHOrIH(rec);
  final int maxNH = params.maxAlignmentCount();
  if (maxNH >= 0 && nh != null && nh > maxNH) {
      return false;
  }
  if (params.excludeVariantInvalid()) {
    if (nh != null && nh == 0) {
      return false;
    }
    if (!rec.getReadUnmappedFlag() && rec.getAlignmentStart() <= 0) {
      return false;
    }
  }

  final IntegerOrPercentage maxAS = mated ? params.maxMatedAlignmentScore() : params.maxUnmatedAlignmentScore();
  if (maxAS != null) {
    final Integer as = rec.getIntegerAttribute(SamUtils.ATTRIBUTE_ALIGNMENT_SCORE);
    if (as != null && as > maxAS.getValue(rec.getReadLength())) {
      return false;
    }
  }
  final int minReadLength = params.minReadLength();
  if (minReadLength >= 0 && rec.getReadLength() < minReadLength) {
    return false;
  }

  if (params.subsampleFraction() != null) {
    final double sFrac;
    if (params.subsampleRampFraction() != null) { // Use subsample ramping
      if (rec.getAlignmentStart() == 0) {
        sFrac = (params.subsampleFraction() + params.subsampleRampFraction()) / 2;
      } else {
        final int pos = rec.getAlignmentStart();
        final int refLength = rec.getHeader().getSequence(rec.getReferenceIndex()).getSequenceLength();
        final double lengthFrac = Math.max(0, Math.min(1.0, (double) pos / refLength));
        sFrac = params.subsampleFraction() + lengthFrac * (params.subsampleRampFraction() - params.subsampleFraction());
      }
    } else {
      sFrac = params.subsampleFraction();
    }
    // Subsample using hash of read name, ensures pairs are kept together
    return (internalHash(rec.getReadName(), params.subsampleSeed()) & SUBSAMPLE_MASK) < sFrac * SUBSAMPLE_MAX;
  }

  if (params.selectReadGroups() != null) {
    final SAMReadGroupRecord srg = rec.getReadGroup();
    if (srg == null || !params.selectReadGroups().contains(srg.getReadGroupId())) {
      return false;
    }
  }

  return true;
}
 

public boolean qcCheck(SAMRecord sr){
Cigar cigar = sr.getCigar();
int rLen = sr.getReadLength();
int effectiveLen = 0;
if(cigar==null) 
    return false;
else{
    for(final CigarElement ce : cigar.getCigarElements()){
	CigarOperator op = ce.getOperator();
	int cigarLen = ce.getLength();
	switch(op)
	    {
	    case M:
		{
		    effectiveLen += cigarLen;
		    break;
		}
	    case I:
		{
		    effectiveLen += cigarLen;
		    break;
		}
	    default: 
		break;
	    }
    }
}
boolean readdebug = false;
if(readdebug){
    HLA.log.appendln(sr.getSAMString());
    HLA.log.appendln("EffectiveLen:\t" + effectiveLen);
    HLA.log.appendln("ReadLen:\t" + rLen);
}
Integer i = sr.getIntegerAttribute("NM");
int nm = 0;
if(i!=null)
    nm = i.intValue();
if(readdebug)
    HLA.log.appendln("NM=\t" + nm);
if(nm < 16){
    if(readdebug)
	HLA.log.appendln("PASSWED QC");
    return true;
}
if(readdebug){
    HLA.log.appendln("FAILED QC");
    HLA.log.appendln(sr.getSAMString());
}
return false;
   }