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

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

/**
 * 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 areReadsMates(final SAMRecord read1, final SAMRecord read2) {
    // must have same name
    return (read1.getReadName().equals(read2.getReadName()) &&
            // must be paired
            read1.getReadPairedFlag() &&
            // one must be first while the other is not
            read1.getFirstOfPairFlag() != read2.getFirstOfPairFlag() &&
            // one must be second while the other is not
            read1.getSecondOfPairFlag() != read2.getSecondOfPairFlag() &&
            // read1 must be mapped
            !read1.getReadUnmappedFlag() &&
            // read2 must be mapped
            !read2.getReadUnmappedFlag() &&
            // read1 must be non-secondary
            !read1.isSecondaryAlignment() &&
            // read2 must be non-secondary
            !read2.isSecondaryAlignment() &&
            // read1 mate position must agree with read2's position
            read1.getMateAlignmentStart() == read2.getAlignmentStart() &&
            // read1 mate reference must agree with read2's reference
            read1.getMateReferenceIndex().equals(read2.getReferenceIndex())
    );
}
 

static byte getFlags(SAMRecord record) {
  byte flags = 0;
  if (record.getReadPairedFlag()) {
    flags ^= READ_PAIRED_FLAG;
    if (record.getFirstOfPairFlag()) {
      flags ^= FIRST_OF_PAIR_FLAG;
    }
  }
  if (record.getReadNegativeStrandFlag()) {
    flags ^= READ_STRAND_FLAG;
  }
  if (record.isSecondaryAlignment()) {
    flags ^= NOT_PRIMARY_ALIGNMENT_FLAG;
  }
  return flags;
}
 

boolean retainByReadNumber (final SAMRecord r, final int desiredReadNumber) {
	boolean readPaired = r.getReadPairedFlag();
	// if the read isn't paired, then there's no read number, so keep it.
	if (!readPaired)
		return true;
	// read is paired, is it the first of pair?
	boolean firstRead = r.getFirstOfPairFlag();
	// if you're looking for the first read and this isn't, or looking for the 2nd read and this isn't, then go to the next read.
	if ((firstRead && desiredReadNumber!=1) || (!firstRead && desiredReadNumber==1)) return (false);
	return true;
}
 

private void countRecord(final SAMRecord rec) {
    final boolean isNegativeStrand = rec.getReadNegativeStrandFlag();
    final boolean isReadTwo = rec.getReadPairedFlag() && rec.getSecondOfPairFlag();
    if (isReadTwo) {
        if (isNegativeStrand) this.R2_NEG++;
        else this.R2_POS++;
    } else {
        if (isNegativeStrand) this.R1_NEG++;
        else this.R1_POS++;
    }
}
 

public int getStreamIndex(SAMRecord record) {
	if (!record.getReadPairedFlag())
		return 0;
	if (record.getFirstOfPairFlag())
		return 1;
	if (record.getSecondOfPairFlag())
		return 2;

	return 0;
}
 

@Override
public void addAlignment(SAMRecord alignment) {
	PrintStream ps = s1;
	if (s2 != null && alignment.getReadPairedFlag() && alignment.getFirstOfPairFlag())
		ps = s2;

	printFastq(ps, alignment);
}
 

BaseDistributionMetricCollection gatherBaseQualities (final File input, final int readNumber) {
	ProgressLogger p = new ProgressLogger(this.log);

	SamReader inputSam = SamReaderFactory.makeDefault().open(input);
	BaseDistributionMetricCollection c = new BaseDistributionMetricCollection();

	MAIN_LOOP:
	for (final SAMRecord samRecord : inputSam) {

		p.record(samRecord);
		if (samRecord.isSecondaryOrSupplementary()) continue;
		boolean readPaired = samRecord.getReadPairedFlag();

		boolean firstRead=false;
		if (!readPaired & readNumber==2)
			continue;
		else if (!readPaired & readNumber==1)
			firstRead=true;
		else
			firstRead = samRecord.getFirstOfPairFlag();

		// if you're looking for the first read and this isn't, or looking for the 2nd read and this isn't, then go to the next read.
		if ((firstRead && readNumber!=1) || (!firstRead && readNumber==1)) continue MAIN_LOOP;

		byte [] bases = samRecord.getReadBases();

		for (int i=0; i<bases.length; i++) {
			char base = (char) (bases[i]);
			int idx=i+1;
			c.addBase(base, idx);
		}
	}

	CloserUtil.close(inputSam);
	return (c);


}
 

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

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

@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsWritable(METRICS);

    final SamReader in = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).open(INPUT);
    final SAMFileHeader.SortOrder order = in.getFileHeader().getSortOrder();
    SAMFileWriter out = null;
    if (OUTPUT != null) {
        IOUtil.assertFileIsWritable(OUTPUT);
        out = new SAMFileWriterFactory().makeSAMOrBAMWriter(in.getFileHeader(), true, OUTPUT);
    }

    final Histogram<Integer> histo = new Histogram<Integer>("clipped_bases", "read_count");

    // Combine any adapters and custom adapter pairs from the command line into an array for use in clipping
    final AdapterPair[] adapters;
    {
        final List<AdapterPair> tmp = new ArrayList<AdapterPair>();
        tmp.addAll(ADAPTERS);
        if (FIVE_PRIME_ADAPTER != null && THREE_PRIME_ADAPTER != null) {
            tmp.add(new CustomAdapterPair(FIVE_PRIME_ADAPTER, THREE_PRIME_ADAPTER));
        }
        adapters = tmp.toArray(new AdapterPair[tmp.size()]);
    }

    ////////////////////////////////////////////////////////////////////////
    // Main loop that consumes reads, clips them and writes them to the output
    ////////////////////////////////////////////////////////////////////////
    final ProgressLogger progress = new ProgressLogger(log, 1000000, "Read");
    final SAMRecordIterator iterator = in.iterator();

    final AdapterMarker adapterMarker = new AdapterMarker(ADAPTER_TRUNCATION_LENGTH, adapters).
            setMaxPairErrorRate(MAX_ERROR_RATE_PE).setMinPairMatchBases(MIN_MATCH_BASES_PE).
            setMaxSingleEndErrorRate(MAX_ERROR_RATE_SE).setMinSingleEndMatchBases(MIN_MATCH_BASES_SE).
            setNumAdaptersToKeep(NUM_ADAPTERS_TO_KEEP).
            setThresholdForSelectingAdaptersToKeep(PRUNE_ADAPTER_LIST_AFTER_THIS_MANY_ADAPTERS_SEEN);

    while (iterator.hasNext()) {
        final SAMRecord rec = iterator.next();
        final SAMRecord rec2 = rec.getReadPairedFlag() && iterator.hasNext() ? iterator.next() : null;
        rec.setAttribute(ReservedTagConstants.XT, null);

        // Do the clipping one way for PE and another for SE reads
        if (rec.getReadPairedFlag()) {
            // Assert that the input file is in query name order only if we see some PE reads
            if (order != SAMFileHeader.SortOrder.queryname) {
                throw new PicardException("Input BAM file must be sorted by queryname");
            }

            if (rec2 == null) throw new PicardException("Missing mate pair for paired read: " + rec.getReadName());
            rec2.setAttribute(ReservedTagConstants.XT, null);

            // Assert that we did in fact just get two mate pairs
            if (!rec.getReadName().equals(rec2.getReadName())) {
                throw new PicardException("Adjacent reads expected to be mate-pairs have different names: " +
                        rec.getReadName() + ", " + rec2.getReadName());
            }

            // establish which of pair is first and which second
            final SAMRecord first, second;

            if (rec.getFirstOfPairFlag() && rec2.getSecondOfPairFlag()) {
                first = rec;
                second = rec2;
            } else if (rec.getSecondOfPairFlag() && rec2.getFirstOfPairFlag()) {
                first = rec2;
                second = rec;
            } else {
                throw new PicardException("Two reads with same name but not correctly marked as 1st/2nd of pair: " + rec.getReadName());
            }

            adapterMarker.adapterTrimIlluminaPairedReads(first, second);
        } else {
            adapterMarker.adapterTrimIlluminaSingleRead(rec);
        }

        // Then output the records, update progress and metrics
        for (final SAMRecord r : new SAMRecord[]{rec, rec2}) {
            if (r != null) {
                progress.record(r);
                if (out != null) out.addAlignment(r);

                final Integer clip = r.getIntegerAttribute(ReservedTagConstants.XT);
                if (clip != null) histo.increment(r.getReadLength() - clip + 1);
            }
        }
    }

    if (out != null) out.close();

    // Lastly output the metrics to file
    final MetricsFile<?, Integer> metricsFile = getMetricsFile();
    metricsFile.setHistogram(histo);
    metricsFile.write(METRICS);

    CloserUtil.close(in);
    return 0;
}
 

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

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

private void updateTotalClusters(final SAMRecord rec, final Map<String, GcObject> gcData) {
    for (final Map.Entry<String, GcObject> entry : gcData.entrySet()) {
        final GcObject gcCur = entry.getValue();
        if (!rec.getReadPairedFlag() || rec.getFirstOfPairFlag()) ++gcCur.totalClusters;
    }
}
 

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

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

/**
 * @return The ith hit for a un-paired read.  Never returns null.
 * Do not call if paired read.
 */
public SAMRecord getFragment(final int i) {
    final SAMRecord samRecord = firstOfPairOrFragment.get(i);
    if (samRecord.getReadPairedFlag()) throw new UnsupportedOperationException("getFragment called for paired read");
    return samRecord;
}
 

private HitsForInsert nextMaybeEmpty() {
    if (!peekIterator.hasNext()) throw new IllegalStateException();
    final String readName = peekIterator.peek().getReadName();
    final HitsForInsert hits = new HitsForInsert();

    Boolean isPaired = null;

    // Accumulate the alignments matching readName.
    do {
        final SAMRecord rec = peekIterator.next();
        replaceHardWithSoftClips(rec);
        // It is critical to do this here, because SamAlignmentMerger uses this exception to determine
        // if the aligned input needs to be sorted.
        if (peekIterator.hasNext() && queryNameComparator.fileOrderCompare(rec, peekIterator.peek()) > 0) {
            throw new IllegalStateException("Underlying iterator is not queryname sorted: " +
            rec + " > " + peekIterator.peek());
        }

        if (isPaired == null) {
            isPaired = rec.getReadPairedFlag();
        } else if (isPaired != rec.getReadPairedFlag()) {
            throw new PicardException("Got a mix of paired and unpaired alignments for read " + readName);
        }

        // Records w/ a supplemental flag are stashed to the side until the primary alignment has
        // been determined, and then re-added into the process later
        if (!rec.getReadPairedFlag() || rec.getFirstOfPairFlag()) {
            if (rec.getSupplementaryAlignmentFlag()) {
                hits.addSupplementalFirstOfPairOrFragment(rec);
            } else {
                hits.addFirstOfPairOrFragment(rec);
            }
        } else if (rec.getSecondOfPairFlag()) {
            if (rec.getSupplementaryAlignmentFlag()) {
                hits.addSupplementalSecondOfPair(rec);
            } else {
                hits.addSecondOfPair(rec);
            }
        } else throw new PicardException("Read is marked as pair but neither first or second: " + readName);
    } while (peekIterator.hasNext() && peekIterator.peek().getReadName().equals(readName));

    // If there is no more than one alignment for each end, no need to do any coordination.
    if (hits.numHits() <= 1) {
        // No HI tags needed if only a single hit
        if (hits.getFirstOfPair(0) != null) {
            hits.getFirstOfPair(0).setAttribute(SAMTag.HI.name(), null);
            hits.getFirstOfPair(0).setNotPrimaryAlignmentFlag(false);
        }
        if (hits.getSecondOfPair(0) != null) {
            hits.getSecondOfPair(0).setAttribute(SAMTag.HI.name(), null);
            hits.getSecondOfPair(0).setNotPrimaryAlignmentFlag(false);
        }
    } else {
        primaryAlignmentSelectionStrategy.pickPrimaryAlignment(hits);
    }

    // Used to check that alignments for first and second were correlated, but this is no longer required.
    return hits;
}
 

public PrimaryAlignmentKey(final SAMRecord rec) {
    this.pairStatus = rec.getReadPairedFlag() ?
            (rec.getSecondOfPairFlag() ? PairStatus.SECOND : PairStatus.FIRST) :
            PairStatus.UNPAIRED;
    this.readName = rec.getReadName();
}