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

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

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 static void addDuplicateReadToMetrics(final SAMRecord rec, final DuplicationMetrics metrics) {
    // only update duplicate counts for "decider" reads, not tag-a-long reads
    if (!rec.isSecondaryOrSupplementary() && !rec.getReadUnmappedFlag()) {
        // Update the duplication metrics
        if (!rec.getReadPairedFlag() || rec.getMateUnmappedFlag()) {
            ++metrics.UNPAIRED_READ_DUPLICATES;
        } else {
            ++metrics.READ_PAIR_DUPLICATES;// will need to be divided by 2 at the end
        }
    }
}
 

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

@Override
public void acceptRecord(final SAMRecord record, final ReferenceSequence refSeq) {
    if (!record.getReadPairedFlag() ||
            record.getReadUnmappedFlag() ||
            record.getMateUnmappedFlag() ||
            record.getFirstOfPairFlag() ||
            record.isSecondaryOrSupplementary() ||
            (record.getDuplicateReadFlag() && !this.includeDuplicates) ||
            record.getInferredInsertSize() == 0) {
        return;
    }

    super.acceptRecord(record, refSeq);
}
 

/**
 * Read the data and create a list of lengths
 */
private void read(){
	int counter = 0;
	SAMFileReader reader = new SAMFileReader(bam,index);
	ArrayList<Double> temp = new ArrayList<Double>();
	for (int i = 0; i < genome.size();i++){
		String chr = genome.get(i).getChrom();
		int start = genome.get(i).getStart();
		int stop = genome.get(i).getStop();
		CloseableIterator<SAMRecord> iter = reader.query(chr,start,stop,false);
		while (iter.hasNext()){
			SAMRecord record = null;
			try{
				record = iter.next();
			}
			catch(SAMFormatException ex){
				System.out.println("SAM Record is problematic. Has mapQ != 0 for unmapped read. Will continue anyway");
			}
			if(record != null){
			if(!record.getReadUnmappedFlag() && !record.getMateUnmappedFlag() && record.getFirstOfPairFlag()) {
				if (record.getMappingQuality() >= minQ){
					
					if (Math.abs(record.getInferredInsertSize()) > 100 && Math.abs(record.getInferredInsertSize())
							< 1000){
						counter+=1;
						temp.add((double)Math.abs(record.getInferredInsertSize()));
					}
				}
			}
			}
		}
		iter.close();
	}
	reader.close();
	lengths = new double[counter];
	for (int i = 0;i < temp.size();i++){
		if (temp.get(i) > 100){
			lengths[i] = temp.get(i);
		}
	}
	
}
 

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

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

/**
 * Checks if we can skip the SAM/BAM file when reverting origin base qualities and adding mate cigars.
 *
 * @param inputFile                   the SAM/BAM input file
 * @param maxRecordsToExamine         the maximum number of records to examine before quitting
 * @param revertOriginalBaseQualities true if we are to revert original base qualities, false otherwise
 * @return whether we can skip or not, and the explanation why.
 */
public static CanSkipSamFile canSkipSAMFile(final File inputFile, final int maxRecordsToExamine, boolean revertOriginalBaseQualities,
                                            final File referenceFasta) {
    final SamReader in = SamReaderFactory.makeDefault().referenceSequence(referenceFasta).enable(SamReaderFactory.Option.EAGERLY_DECODE).open(inputFile);
    final Iterator<SAMRecord> iterator = in.iterator();
    int numRecordsExamined = 0;
    CanSkipSamFile returnType = CanSkipSamFile.FOUND_NO_EVIDENCE;

    while (iterator.hasNext() && numRecordsExamined < maxRecordsToExamine) {
        final SAMRecord record = iterator.next();

        if (revertOriginalBaseQualities && null != record.getOriginalBaseQualities()) {
            // has OQ, break and return case #2
            returnType = CanSkipSamFile.CANNOT_SKIP_FOUND_OQ;
            break;
        }

        // check if mate pair and its mate is mapped
        if (record.getReadPairedFlag() && !record.getMateUnmappedFlag()) {
            if (null == SAMUtils.getMateCigar(record)) {
                // has no MC, break and return case #2
                returnType = CanSkipSamFile.CANNOT_SKIP_FOUND_NO_MC;
                break;
            } else {
                // has MC, previously checked that it does not have OQ, break and return case #1
                returnType = CanSkipSamFile.CAN_SKIP;
                break;
            }
        }

        numRecordsExamined++;
    }

    // no more records anyhow, so we can skip
    if (!iterator.hasNext() && CanSkipSamFile.FOUND_NO_EVIDENCE == returnType) {
        returnType = CanSkipSamFile.CAN_SKIP;
    }

    CloserUtil.close(in);

    return returnType;
}
 

private static boolean isPairedAndBothMapped(final SAMRecord record) {
    return record.getReadPairedFlag() &&
            !record.getReadUnmappedFlag() &&
            !record.getMateUnmappedFlag();
}
 

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 static boolean isMappedPair(final SAMRecord rec) {
    return rec.getReadPairedFlag() && !rec.getReadUnmappedFlag() && !rec.getMateUnmappedFlag();
}
 

@Override
public boolean reallyFilterOut(final SAMRecord record) { return !record.getReadPairedFlag() || record.getMateUnmappedFlag(); }
 

@ProcessElement
public void processElement(DoFn<Read, String>.ProcessContext c, BoundedWindow window)
    throws Exception {

  this.window = window;

  if (headerInfo == null) {
    headerInfo = c.sideInput(headerView);
  }
  final Read read = c.element();

  if (readCount == 0) {

    shardContig = KeyReadsFn.shardKeyForRead(read, 1);
    sequenceIndex = headerInfo.header.getSequenceIndex(shardContig.referenceName);
    final boolean isFirstShard = headerInfo.shardHasFirstRead(shardContig);
    final String outputFileName = options.getOutput();
    shardName = outputFileName + "-" + String.format("%012d", sequenceIndex) + "-"
        + shardContig.referenceName
        + ":" + String.format("%012d", shardContig.start);
    LOG.info("Writing shard file " + shardName);
    final OutputStream outputStream =
        Channels.newOutputStream(
            new GcsUtil.GcsUtilFactory().create(options)
              .create(GcsPath.fromUri(shardName),
                  BAMIO.BAM_INDEX_FILE_MIME_TYPE));
    ts = new TruncatedOutputStream(
        outputStream, BlockCompressedStreamConstants.EMPTY_GZIP_BLOCK.length);
    bw = new BAMBlockWriter(ts, null /*file*/);
    bw.setSortOrder(headerInfo.header.getSortOrder(), true);
    bw.setHeader(headerInfo.header);
    if (isFirstShard) {
      LOG.info("First shard - writing header to " + shardName);
      bw.writeHeader(headerInfo.header);
    }
  }
  SAMRecord samRecord = ReadUtils.makeSAMRecord(read, headerInfo.header);
  if (prevRead != null && prevRead.getAlignmentStart() > samRecord.getAlignmentStart()) {
    LOG.info("Out of order read " + prevRead.getAlignmentStart() + " " +
        samRecord.getAlignmentStart() + " during writing of shard " + shardName +
        " after processing " + readCount + " reads, min seen alignment is " +
        minAlignment + " and max is " + maxAlignment + ", this read is " +
        (samRecord.getReadUnmappedFlag() ? "unmapped" : "mapped") + " and its mate is " +
        (samRecord.getMateUnmappedFlag() ? "unmapped" : "mapped"));
    Metrics.counter(WriteBAMFn.class, "Out of order reads").inc();
    readCount++;
    hadOutOfOrder = true;
    return;
  }
  minAlignment = Math.min(minAlignment, samRecord.getAlignmentStart());
  maxAlignment = Math.max(maxAlignment, samRecord.getAlignmentStart());
  prevRead = samRecord;
  if (samRecord.getReadUnmappedFlag()) {
    if (!samRecord.getMateUnmappedFlag()) {
      samRecord.setReferenceName(samRecord.getMateReferenceName());
      samRecord.setAlignmentStart(samRecord.getMateAlignmentStart());
    }
    unmappedReadCount++;
  }
  bw.addAlignment(samRecord);
  readCount++;
}