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

public Map<String, ReadQualityMetrics> gatherMetrics(final File inputSamOrBamFile) {
	ProgressLogger p = new ProgressLogger(this.log);
       Map<String, ReadQualityMetrics> result = new TreeMap<>();

	ReadQualityMetrics globalMetrics = new ReadQualityMetrics(this.MINIMUM_MAPPING_QUALITY, this.GLOBAL, true);

	SamReader in = SamReaderFactory.makeDefault().open(INPUT);

	for (final SAMRecord r : in)
		if (!r.getReadFailsVendorQualityCheckFlag() || INCLUDE_NON_PF_READS) {
               p.record(r);

               globalMetrics.addRead(r);
               // gather per tag metrics if required.
               result = addMetricsPerTag(r, result);
           }

	CloserUtil.close(in);
	result.put(this.GLOBAL, globalMetrics);
	return (result);
}
 

@Override
protected void acceptRead(final SAMRecord rec, final ReferenceSequence ref) {
    // Skip unwanted records
    if (PF_READS_ONLY && rec.getReadFailsVendorQualityCheckFlag()) return;
    if (ALIGNED_READS_ONLY && rec.getReadUnmappedFlag()) return;
    if (rec.isSecondaryOrSupplementary()) return;

    final byte[] bases = rec.getReadBases();
    final byte[] quals = rec.getBaseQualities();
    final byte[] oq    = rec.getOriginalBaseQualities();

    final int length = quals.length;

    for (int i=0; i<length; ++i) {
        if (INCLUDE_NO_CALLS || !SequenceUtil.isNoCall(bases[i])) {
            qCounts[quals[i]]++;
            if (oq != null) oqCounts[oq[i]]++;
        }
    }
}
 

/**
 * Return the next usable read in the iterator
 *
 * @param iterator the iterator to pull from
 * @param justPeek if true, just peek the next usable read rather than pulling it (note: it may remove unusable reads from the iterator)
 * @return the next read or null if none are left
 */
private static SAMRecord getNextUsableRead(final PeekableIterator<SAMRecord> iterator, final boolean justPeek) {

    while (iterator.hasNext()) {
        // trash the next read if it fails PF, is secondary, or is supplementary
        final SAMRecord nextRead = iterator.peek();
        if (nextRead.getReadFailsVendorQualityCheckFlag() || nextRead.isSecondaryOrSupplementary()) {
            iterator.next();
        }
        // otherwise, return it
        else {
            return justPeek ? nextRead : iterator.next();
        }
    }

    // no good reads left
    return null;
}
 

@Override
       public void acceptRecord(final SAMRecord rec) {
           if (MT_SEQUENCE!=null &&  MT_SEQUENCE.contains(rec.getReferenceName()) && !rec.getReadFailsVendorQualityCheckFlag() &&
                   !rec.isSecondaryOrSupplementary() && !rec.getReadUnmappedFlag()) {

               metrics.PF_BASES += rec.getReadLength();
               final int numAlignedBases = getNumAlignedBases(rec);
               castMetrics().MT_BASES += numAlignedBases;
               metrics.PF_ALIGNED_BASES += numAlignedBases;
           } else
super.acceptRecord(rec);
       }
 

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
}
 

private void handleRecord(final SAMRecord currentRecord, final Map<SAMReadGroupRecord, FastqWriters> writers,
                          final Map<SAMReadGroupRecord, List<FastqWriter>> additionalWriters,
                          final Map<String, SAMRecord> firstSeenMates) {
    if (currentRecord.isSecondaryOrSupplementary() && !INCLUDE_NON_PRIMARY_ALIGNMENTS) {
        return;
    }

    // Skip non-PF reads as necessary
    if (currentRecord.getReadFailsVendorQualityCheckFlag() && !INCLUDE_NON_PF_READS) {
        return;
    }

    final FastqWriters fq = writers.get(currentRecord.getReadGroup());
    SAMRecord read1 = null;
    SAMRecord read2 = null;
    if (currentRecord.getReadPairedFlag()) {
        final String currentReadName = currentRecord.getReadName();
        final SAMRecord firstRecord = firstSeenMates.remove(currentReadName);
        if (firstRecord == null) {
            firstSeenMates.put(currentReadName, currentRecord);
        } else {
            assertPairedMates(firstRecord, currentRecord);

            read1 = currentRecord.getFirstOfPairFlag() ? currentRecord : firstRecord;
            read2 = currentRecord.getFirstOfPairFlag() ? firstRecord : currentRecord;
            writeRecord(read1, 1, fq.getFirstOfPair(), READ1_TRIM, READ1_MAX_BASES_TO_WRITE);
            final FastqWriter secondOfPairWriter = fq.getSecondOfPair();
            if (secondOfPairWriter == null) {
                throw new PicardException("Input contains paired reads but no SECOND_END_FASTQ specified.");
            }
            writeRecord(read2, 2, secondOfPairWriter, READ2_TRIM, READ2_MAX_BASES_TO_WRITE);
        }
    } else {
        writeRecord(currentRecord, null, fq.getUnpaired(), READ1_TRIM, READ1_MAX_BASES_TO_WRITE);
    }

    handleAdditionalRecords(currentRecord, additionalWriters, read1, read2);
}
 

@Override
protected void acceptRead(final SAMRecord rec, final ReferenceSequence ref) {
    if ((PF_READS_ONLY) && (rec.getReadFailsVendorQualityCheckFlag())) {
        return;
    }
    if ((ALIGNED_READS_ONLY) && (rec.getReadUnmappedFlag())) {
        return;
    }
    if (rec.isSecondaryOrSupplementary()) {
        return;
    }
    hist.addRecord(rec);
}
 

@Override
protected void acceptRead(final SAMRecord rec, final ReferenceSequence ref) {
    // Skip unwanted records
    if (PF_READS_ONLY && rec.getReadFailsVendorQualityCheckFlag()) return;
    if (ALIGNED_READS_ONLY && rec.getReadUnmappedFlag()) return;
    if (rec.isSecondaryOrSupplementary()) return;

    q.addRecord(rec);
    oq.addRecord(rec);
}
 

/**
 * Count the number of records in the bamFile that could potentially be written
 *
 * @return  the number of records in the Bam file
 */
private int countWritableRecords() {
    int count = 0;

    final SamReader reader = SamReaderFactory.makeDefault().open(this.bamFile);
    if(!reader.getFileHeader().getSortOrder().equals(SAMFileHeader.SortOrder.queryname)) {
    	//this is a fix for issue PIC-274: It looks like BamToBfqWriter requires that the input BAM is queryname sorted, 
    	//but it doesn't check this early, nor produce an understandable error message."
    	throw new PicardException("Input file (" + this.bamFile.getAbsolutePath() +") needs to be sorted by queryname.");
    }
    final PeekableIterator<SAMRecord> it = new PeekableIterator<SAMRecord>(reader.iterator());
    if (!this.pairedReads) {
        // Filter out noise reads and reads that fail the quality filter
        final List<SamRecordFilter> filters = new ArrayList<SamRecordFilter>();
        filters.add(new TagFilter(ReservedTagConstants.XN, 1));
        if (!this.includeNonPfReads) {
            filters.add(new FailsVendorReadQualityFilter());
        }
        final FilteringSamIterator itr = new FilteringSamIterator(it, new AggregateFilter(filters));
        while (itr.hasNext()) {
            itr.next();
            count++;
        }
    }
    else {
        while (it.hasNext()) {
            final SAMRecord first = it.next();
            final SAMRecord second = it.next();
            // If both are noise reads, filter them out
            if (first.getAttribute(ReservedTagConstants.XN) != null &&
                second.getAttribute(ReservedTagConstants.XN) != null)  {
                // skip it
            }
            // If either fails to pass filter, then exclude them as well
            else if (!this.includeNonPfReads && (first.getReadFailsVendorQualityCheckFlag() || second.getReadFailsVendorQualityCheckFlag()) ) {
                // skip it
            }
            // Otherwise, write them out
            else {
                count++;
            }
        }
    }
    it.close();
    CloserUtil.close(reader);
    return count;
}
 

/**
 * This is factored out of doWork only for unit testing.
 */
int writeSamText(PrintStream printStream) {
    try {
        final CloseableIterator<SAMRecord> samRecordsIterator;
        final SamReader samReader = SamReaderFactory.makeDefault()
                .referenceSequence(REFERENCE_SEQUENCE)
                .open(SamInputResource.of(INPUT));

        // if we are only using the header or we aren't using intervals, then use the reader as the iterator.
        // otherwise use the SamRecordIntervalIteratorFactory to make an interval-ed iterator
        if (HEADER_ONLY || INTERVAL_LIST == null) {
            samRecordsIterator = samReader.iterator();
        } else {
            IOUtil.assertFileIsReadable(INTERVAL_LIST);

            final List<Interval> intervals = IntervalList.fromFile(INTERVAL_LIST).uniqued().getIntervals();
            samRecordsIterator = new SamRecordIntervalIteratorFactory().makeSamRecordIntervalIterator(samReader, intervals, samReader.hasIndex());
        }
        final AsciiWriter writer = new AsciiWriter(printStream);
        final SAMFileHeader header = samReader.getFileHeader();
        if (!RECORDS_ONLY) {
            if (header.getTextHeader() != null) {
                writer.write(header.getTextHeader());
            } else {
                // Headers that are too large are not retained as text, so need to regenerate text
                new SAMTextHeaderCodec().encode(writer, header, true);
            }
        }
        if (!HEADER_ONLY) {
            while (samRecordsIterator.hasNext()) {
                final SAMRecord rec = samRecordsIterator.next();

                if (printStream.checkError()) {
                    return 1;
                }

                if (this.ALIGNMENT_STATUS == AlignmentStatus.Aligned && rec.getReadUnmappedFlag()) continue;
                if (this.ALIGNMENT_STATUS == AlignmentStatus.Unaligned && !rec.getReadUnmappedFlag()) continue;

                if (this.PF_STATUS == PfStatus.PF && rec.getReadFailsVendorQualityCheckFlag()) continue;
                if (this.PF_STATUS == PfStatus.NonPF && !rec.getReadFailsVendorQualityCheckFlag()) continue;
                writer.write(rec.getSAMString());
            }
        }
        writer.flush();
        if (printStream.checkError()) {
            return 1;
        }
        CloserUtil.close(writer);
        CloserUtil.close(samRecordsIterator);
        return 0;
    } catch (IOException e) {
        throw new PicardException("Exception writing SAM text", e);
    }
}
 

public void acceptRecord(final SAMRecord rec, final ReferenceSequence ref) {
    if (!this.includeSecondaryAlignments    && rec.getNotPrimaryAlignmentFlag()) return;
    if (!this.includeSupplementalAlignments && rec.getSupplementaryAlignmentFlag()) return;

    final int length = rec.getReadLength();
    metrics.TOTAL_READS++;
    metrics.TOTAL_BASES += length;

    final boolean isPfRead = !rec.getReadFailsVendorQualityCheckFlag();
    if (isPfRead) {
        metrics.PF_READS++;
        metrics.PF_BASES += length;
    }

    final byte[] quals;
    if (this.useOriginalQualities) {
        byte[] tmp = rec.getOriginalBaseQualities();
        if (tmp == null) tmp = rec.getBaseQualities();
        quals = tmp;
    } else {
        quals = rec.getBaseQualities();
    }

    // add up quals, and quals >= 20
    for (final int qual : quals) {
        metrics.Q20_EQUIVALENT_YIELD += qual;

        if (qual >= 30) {
            metrics.Q20_BASES++;
            metrics.Q30_BASES++;
        } else if (qual >= 20) {
            metrics.Q20_BASES++;
        }

        if (isPfRead) {
            metrics.PF_Q20_EQUIVALENT_YIELD += qual;
            if (qual >= 30) {
                metrics.PF_Q20_BASES++;
                metrics.PF_Q30_BASES++;
            } else if (qual >= 20) {
                metrics.PF_Q20_BASES++;
            }
        }
    }
}
 

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 boolean isHighQualityMapping(final SAMRecord record) {
    return !record.getReadFailsVendorQualityCheckFlag() &&
            record.getMappingQuality() >= MAPPING_QUALITY_THRESHOLD;
}