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

/**
 * Determines if the read represented by a SAM record belongs to the top or bottom strand
 * or if it cannot determine strand position due to one of the reads being unmapped.
 * Top strand is defined as having a read 1 unclipped 5' coordinate
 * less than the read 2 unclipped 5' coordinate.  If a read is unmapped
 * we do not attempt to determine the strand to which the read or its mate belongs.
 * If the mate belongs to a different contig from the read, then the reference
 * index for the read and its mate is used in leu of the unclipped 5' coordinate.
 * @param rec Record to determine top or bottom strand
 * @return Top or bottom strand, unknown if it cannot be determined.
 */
static ReadStrand getStrand(final SAMRecord rec) {
    if (rec.getReadUnmappedFlag() || rec.getMateUnmappedFlag()) {
        return ReadStrand.UNKNOWN;
    }

    // If the read pair are aligned to different contigs we use
    // the reference index to determine relative 5' coordinate ordering.
    // Both the read and its mate should not have their unmapped flag set to true.
    if (!rec.getReferenceIndex().equals(rec.getMateReferenceIndex())) {
        if (rec.getFirstOfPairFlag() == (rec.getReferenceIndex() < rec.getMateReferenceIndex())) {
            return ReadStrand.TOP;
        } else {
            return ReadStrand.BOTTOM;
        }
    }

    final int read5PrimeStart = (rec.getReadNegativeStrandFlag()) ? rec.getUnclippedEnd() : rec.getUnclippedStart();
    final int mate5PrimeStart = (rec.getMateNegativeStrandFlag()) ? SAMUtils.getMateUnclippedEnd(rec) : SAMUtils.getMateUnclippedStart(rec);

    if (rec.getFirstOfPairFlag() == (read5PrimeStart <= mate5PrimeStart)) {
        return ReadStrand.TOP;
    } else {
        return ReadStrand.BOTTOM;
    }
}
 

/**
 * This is supposed to check if the mates have valid pairing flags.
 * 
 * @param r1
 * @param r2
 * @return
 */
private boolean checkMateFlags(SAMRecord r1, SAMRecord r2) {
	if (!r1.getReadPairedFlag() || !r2.getReadPairedFlag())
		return false;

	if (r1.getReadUnmappedFlag() != r2.getMateUnmappedFlag())
		return false;
	if (r1.getReadNegativeStrandFlag() != r2.getMateNegativeStrandFlag())
		return false;
	if (r1.getProperPairFlag() != r2.getProperPairFlag())
		return false;
	if (r1.getFirstOfPairFlag() && r2.getFirstOfPairFlag())
		return false;
	if (r1.getSecondOfPairFlag() && r2.getSecondOfPairFlag())
		return false;

	if (r2.getReadUnmappedFlag() != r1.getMateUnmappedFlag())
		return false;
	if (r2.getReadNegativeStrandFlag() != r1.getMateNegativeStrandFlag())
		return false;

	return true;
}
 

private boolean isCandidateRecord(final SAMRecord record)
{
    int fragmentLength = abs(record.getInferredInsertSize());
    if(fragmentLength > MAX_FRAGMENT_LENGTH)
        return false;

    // ignore translocations and inversions
    if(!record.getMateReferenceName().equals(record.getReferenceName()) || record.getMateNegativeStrandFlag() == record.getReadNegativeStrandFlag())
        return false;

    // ignore split and soft-clipped reads above the read length
    if(record.getCigar().containsOperator(CigarOperator.N) || !record.getCigar().containsOperator(CigarOperator.M))
        return false;

    int readLength = max(mMaxReadLength, mConfig.ReadLength);

    if(readLength > 0 && fragmentLength > readLength && record.getCigar().containsOperator(CigarOperator.S))
    {
        // permit small soft-clips up to a point and then none for longer fragments
        if(record.getCigar().getCigarElements().stream().anyMatch(x -> x.getOperator() == CigarOperator.S && x.getLength() > 2))
            return false;
    }

    // both reads must fall in the current gene
    int otherStartPos = record.getMateAlignmentStart();
    if(!positionWithin(otherStartPos, mCurrentGenesRange[SE_START], mCurrentGenesRange[SE_END]))
        return false;

    // reads cannot cover any part of an exon
    int posStart = record.getStart();
    int posEnd = record.getEnd();

    for(final TranscriptData transData : mCurrentTransDataList)
    {
        if(transData.exons().stream().anyMatch(x -> positionsOverlap(posStart, posEnd, x.ExonStart, x.ExonEnd)))
            return false;
    }

    return true;
}
 

public MateKey getMateKey(SAMRecord read) {
	
	// If mate is mapped, use read flag.
	// If mate is not mapped, use opposite of this read's RC flag
	boolean isMateRevOrientation = read.getMateUnmappedFlag() ? !read.getReadNegativeStrandFlag() : read.getMateNegativeStrandFlag();
	int matePos = read.getMateUnmappedFlag() ? -1 : read.getMateAlignmentStart();
	
	int mateNum = read.getFirstOfPairFlag() ? 2 : 1;
	
	return new MateKey(read.getReadName(), matePos,
			read.getMateUnmappedFlag(), isMateRevOrientation, mateNum, read.getAlignmentStart());
}
 

public MateKey getOriginalReadInfo(SAMRecord read) {
		int pos = read.getAlignmentStart();
		boolean isUnmapped = read.getReadUnmappedFlag();
		boolean isRc = read.getReadNegativeStrandFlag();
		
		String yo = read.getStringAttribute("YO");
		
		if (yo != null) {
			if (yo.startsWith("N/A")) {
				// Original alignment was unmapped
				isUnmapped = true;
//				isRc = false;
				// Orientation is forced to be opposite of mate during realignment
				// regardless of the original alignment.
				isRc = !read.getMateNegativeStrandFlag();
			} else {
				String[] fields = yo.split(":");
				pos = Integer.parseInt(fields[1]);
				isUnmapped = false;
				isRc = fields[2].equals("-") ? true : false;
			}
		}
		
		int readNum = read.getFirstOfPairFlag() ? 1 : 2;
		
		return new MateKey(read.getReadName(), pos, isUnmapped, isRc, readNum, read.getAlignmentStart());
	}
 

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

public boolean checkDuplicates(final SAMRecord record)
{
    if(record.getDuplicateReadFlag())
        return true;

    if(!mMarkDuplicates)
        return false;

    if(mDuplicateReadIds.contains(record.getReadName()))
    {
        mDuplicateReadIds.remove(record.getReadName());
        return true;
    }

    if(!record.getReferenceName().equals(record.getMateReferenceName()) || record.getReadNegativeStrandFlag() == record.getMateNegativeStrandFlag())
        return false;

    int firstStartPos = record.getFirstOfPairFlag() ? record.getStart() : record.getMateAlignmentStart();
    int secondStartPos = record.getFirstOfPairFlag() ? record.getMateAlignmentStart() : record.getStart();
    int readLength = record.getReadLength();
    int insertSize = record.getInferredInsertSize();

    List<int[]> dupDataList = mDuplicateCache.get(firstStartPos);

    if(dupDataList == null)
    {
        dupDataList = Lists.newArrayList();
        mDuplicateCache.put(firstStartPos, dupDataList);
    }
    else
    {
        // search for a match
        if(dupDataList.stream().anyMatch(x -> x[DUP_DATA_SECOND_START] == secondStartPos
                && x[DUP_DATA_READ_LEN] == readLength && insertSize == x[DUP_DATA_INSERT_SIZE]))
        {
            ISF_LOGGER.trace("duplicate fragment: id({}) chr({}) pos({}->{}) otherReadStart({}) insertSize({})",
                    record.getReadName(), record.getReferenceName(), firstStartPos, record.getEnd(), secondStartPos, insertSize);

            // cache so the second read can be identified immediately
            mDuplicateReadIds.add(record.getReadName());
            return true;
        }
    }

    int[] dupData = {secondStartPos, readLength, insertSize};
    dupDataList.add(dupData);

    return false;
}
 

public static int mergeReadPair(SAMRecordWrapper readWrapper, Map<String, SAMRecordWrapper> firstReads, Map<String, SAMRecordWrapper> secondReads) {
	
	int alignmentStart = -1;
	SAMRecord read = readWrapper.getSamRecord();
	
	if (read.getReadPairedFlag() && !read.getReadUnmappedFlag() && !read.getMateUnmappedFlag() &&
			read.getReadNegativeStrandFlag() != read.getMateNegativeStrandFlag()) {
		
		SAMRecordWrapper pair = null;
		
		if (read.getFirstOfPairFlag()) {
			pair = secondReads.get(read.getReadName() + "_" + read.getMateAlignmentStart());
		} else {
			pair = firstReads.get(read.getReadName() + "_" + read.getMateAlignmentStart());
		}
		
		if (pair != null) {
			SAMRecordWrapper first = null;
			SAMRecordWrapper second = null;
			if (read.getReadNegativeStrandFlag()) {
				first = pair;
				second = readWrapper;
			} else {
				first = readWrapper;
				second = pair;
			}

			
			if (first.getAdjustedAlignmentStart() < second.getAdjustedAlignmentStart() &&
					first.getAdjustedAlignmentEnd() > second.getAdjustedAlignmentStart() &&
					first.getSamRecord().getReadLength() > MIN_OVERLAP && 
					second.getSamRecord().getReadLength() > MIN_OVERLAP) {

				Pair<String, String> merged = mergeSequences(first.getSamRecord().getReadString(), second.getSamRecord().getReadString(),
						first.getSamRecord().getBaseQualityString(), second.getSamRecord().getBaseQualityString());
				
				if (merged != null) {
					readWrapper.setMerged(merged.getFirst(), merged.getSecond(), first.getAdjustedAlignmentStart(), second.getAdjustedAlignmentEnd());
					pair.setMerged(merged.getFirst(), merged.getSecond(), first.getAdjustedAlignmentStart(), second.getAdjustedAlignmentEnd());
					
					alignmentStart = first.getAdjustedAlignmentStart();
				}
			}
		}
	}
	
	return alignmentStart;
}
 

public static boolean hasPossibleAdapterReadThrough(SAMRecord read, Map<String, SAMRecordWrapper> firstReads, Map<String, SAMRecordWrapper> secondReads) {
	
	boolean hasReadThrough = false;
	
	// Check for fragment read through in paired end
	if (read.getReadPairedFlag() && !read.getReadUnmappedFlag() && !read.getMateUnmappedFlag() &&
			read.getAlignmentStart() == read.getMateAlignmentStart() &&
			read.getReadNegativeStrandFlag() != read.getMateNegativeStrandFlag()) {
		
		SAMRecordWrapper pair = null;
		
		if (read.getFirstOfPairFlag()) {
			pair = secondReads.get(read.getReadName() + "_" + read.getMateAlignmentStart());
		} else {
			pair = firstReads.get(read.getReadName() + "_" + read.getMateAlignmentStart());
		}
		
		if (pair != null && read.getCigar().getCigarElements().size() > 0 && pair.getSamRecord().getCigar().getCigarElements().size() > 0) {
			
			// Looking for something like:
			//     -------->
			//  <--------
			SAMRecord first = null;
			SAMRecord second = null;
			if (read.getReadNegativeStrandFlag()) {
				first = read;
				second = pair.getSamRecord();
			} else {
				first = pair.getSamRecord();
				second = read;
			}
			
			CigarElement firstElement = first.getCigar().getFirstCigarElement();
			CigarElement lastElement = second.getCigar().getLastCigarElement();
			
			if (firstElement.getOperator() == CigarOperator.S && lastElement.getOperator() == CigarOperator.S) {
				// We likely have read through into adapter here.
				hasReadThrough = true;
			}
		}
	}

	return hasReadThrough;
}