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

public ObjectCounter<String> getBamTagCounts (final Iterator<SAMRecord> iterator, final String tag, final int readQuality, final boolean filterPCRDuplicates) {
    ProgressLogger pl = new ProgressLogger(log, 10000000);

    ObjectCounter<String> counter = new ObjectCounter<>();

    for (final SAMRecord r : new IterableAdapter<>(iterator)) {
        pl.record(r);
        if (filterPCRDuplicates && r.getDuplicateReadFlag()) continue;
        if (r.getMappingQuality()<readQuality) continue;
        if (r.isSecondaryOrSupplementary()) continue;
        //String s1 = r.getStringAttribute(tag);
        String s1 = getAnyTagAsString(r, tag);
        if (s1!=null && s1!="") counter.increment(s1); // if the tag doesn't have a value, don't increment it.

    }
    return (counter);
}
 

/**
 * Checks the first ADAPTER_MATCH_LENGTH bases of the read against known adapter sequences and returns
 * true if the read matches an adapter sequence with MAX_ADAPTER_ERRORS mismsatches or fewer.
 * <p>
 * Only unmapped reads and reads with MQ=0 are considers eligible for being adapter
 */

 public boolean isAdapter(final SAMRecord record) {

    // If the read mapped with mapping quality > 0 we do not consider it an adapter read.
    if (!record.getReadUnmappedFlag() && record.getMappingQuality() != 0) {
        return false;
    }

    // if the read is mapped and is aligned to the reverse strand, it needs to be
    // rev-comp'ed before calling isAdapterSequence.
    final boolean revComp = !record.getReadUnmappedFlag() &&
            record.getReadNegativeStrandFlag();

    final byte[] readBases = record.getReadBases();
    return isAdapterSequence(readBases, revComp);
}
 

public void considerBest(final SAMRecord rec) {
    if (bestMapq == -1) {
        bestMapq = rec.getMappingQuality();
        bestAlignments.add(rec);
    } else {
        final int cmp = SAMUtils.compareMapqs(bestMapq, rec.getMappingQuality());
        if (cmp < 0) {
            bestMapq = rec.getMappingQuality();
            bestAlignments.clear();
            bestAlignments.add(rec);
        } else if (cmp == 0) {
            bestAlignments.add(rec);
        }
    }
}
 

private static boolean testReadCategory(ReadCategory c, SAMRecord record) {
	switch (c.type) {
	case ALL:
		return true;
	case HIGHER_MAPPING_SCORE:
		return record.getMappingQuality() > c.param;
	case LOWER_MAPPING_SCORE:
		return record.getMappingQuality() < c.param;
	case UNPLACED:
		return record.getReadUnmappedFlag();

	default:
		throw new RuntimeException("Unknown read category: " + c.type.name());
	}
}
 

/**
 * Do the work after command line has been parsed.
 * RuntimeException may be thrown by this method, and are reported appropriately.
 *
 * @return program exit status.
 */
@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsWritable(OUTPUT);
    final SamReaderFactory factory = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE);
    if (VALIDATION_STRINGENCY == ValidationStringency.STRICT) {
        factory.validationStringency(ValidationStringency.LENIENT);
    }
    final SamReader reader = factory.open(INPUT);
    final SAMFileWriter writer = new SAMFileWriterFactory().makeSAMOrBAMWriter(reader.getFileHeader(), true, OUTPUT);
    final CloseableIterator<SAMRecord> it = reader.iterator();
    final ProgressLogger progress = new ProgressLogger(Log.getInstance(CleanSam.class));

    // If the read (or its mate) maps off the end of the alignment, clip it
    while (it.hasNext()) {
        final SAMRecord rec = it.next();

        // If the read (or its mate) maps off the end of the alignment, clip it
        AbstractAlignmentMerger.createNewCigarsIfMapsOffEndOfReference(rec);

        // check the read's mapping quality
        if (rec.getReadUnmappedFlag() && 0 != rec.getMappingQuality()) {
            rec.setMappingQuality(0);
        }

        writer.addAlignment(rec);
        progress.record(rec);
    }

    writer.close();
    it.close();
    CloserUtil.close(reader);
    return 0;
}
 

public void pickPrimaryAlignment(final HitsForInsert hitsForInsert) {

        if (hitsForInsert.numHits() == 0) throw new IllegalArgumentException("No alignments to pick from");

        // Gather the earliest alignment(s) with best MAPQ
        final List<Integer> earliestAlignments = new ArrayList<Integer>();
        int earliestMappedBase = Integer.MAX_VALUE;
        int bestMapQ = -1;
        for (int i = 0; i < hitsForInsert.numHits(); ++i) {
            final SAMRecord rec = hitsForInsert.getFragment(i);
            if (rec.getReadUnmappedFlag()) continue;
            final int thisFirstMappedBase = getIndexOfFirstAlignedBase(rec);
            final int thisMapQ = rec.getMappingQuality();
            if (thisFirstMappedBase < earliestMappedBase ||
                    (thisFirstMappedBase == earliestMappedBase && thisMapQ > bestMapQ)) {
                earliestAlignments.clear();
                earliestAlignments.add(i);
                earliestMappedBase = thisFirstMappedBase;
                bestMapQ = thisMapQ;
            } else if (thisFirstMappedBase == earliestMappedBase && thisMapQ == bestMapQ) {
                earliestAlignments.add(i);
            } // else it is not the earliest or the best so skip it.
        }


        if (earliestAlignments.size() == 1) {
            // If only one best, pick it.
            hitsForInsert.setPrimaryAlignment(earliestAlignments.get(0));
        } else {
            // Arbitrarily select one of the best
            hitsForInsert.setPrimaryAlignment(earliestAlignments.get(random.nextInt(earliestAlignments.size())));
        }
    }
 

private double calculateQualityScore(int readBaseIndex, final SAMRecord record, final QualityConfig qualityConfig, int numberOfEvents) {
    final double baseQuality = baseQuality(readBaseIndex, record);
    final int distanceFromReadEdge = readDistanceFromEdge(readBaseIndex, record);

    final int mapQuality = record.getMappingQuality();
    int modifiedMapQuality = qualityConfig.modifiedMapQuality(variant, mapQuality, numberOfEvents, record.getProperPairFlag());
    double modifiedBaseQuality = qualityConfig.modifiedBaseQuality(baseQuality, distanceFromReadEdge);

    return Math.max(0, Math.min(modifiedMapQuality, modifiedBaseQuality));
}
 

private boolean passesFilters(@NotNull final SAMRecord record)
{
    if(record.getMappingQuality() < mMinMappingQuality || record.getReadUnmappedFlag())
        return false;

    if(mDropSupplementaries && (record.getSupplementaryAlignmentFlag() || record.isSecondaryAlignment()))
        return false;

    if(mDropDuplicates && record.getDuplicateReadFlag())
        return false;

    return true;
}
 

@Override
public void pickPrimaryAlignment(final HitsForInsert hitsForInsert) {

    if (hitsForInsert.numHits() == 0) throw new IllegalArgumentException("No alignments to pick from");

    // Gather the earliest alignment(s) with best MAPQ
    final List<Integer> earliestAlignments = new ArrayList<>();
    int earliestMappedBase = Integer.MAX_VALUE;
    int bestMapQ = -1;
    for (int i = 0; i < hitsForInsert.numHits(); ++i) {
        final SAMRecord rec = hitsForInsert.getFragment(i);
        if (rec.getReadUnmappedFlag()) continue;
        final int thisFirstMappedBase = getIndexOfFirstAlignedBase(rec);
        final int thisMapQ = rec.getMappingQuality();
        if (thisFirstMappedBase < earliestMappedBase ||
                (thisFirstMappedBase == earliestMappedBase && thisMapQ > bestMapQ)) {
            earliestAlignments.clear();
            earliestAlignments.add(i);
            earliestMappedBase = thisFirstMappedBase;
            bestMapQ = thisMapQ;
        } else if (thisFirstMappedBase == earliestMappedBase && thisMapQ == bestMapQ) {
            earliestAlignments.add(i);
        } // else it is not the earliest or the best so skip it.
    }


    if (earliestAlignments.size() == 1) {
        // If only one best, pick it.
        hitsForInsert.setPrimaryAlignment(earliestAlignments.get(0));
    } else {
        // Arbitrarily select one of the best
        hitsForInsert.setPrimaryAlignment(earliestAlignments.get(random.nextInt(earliestAlignments.size())));
    }
}
 

private boolean samRecordMeetsQualityRequirements(@NotNull final SAMRecord record) {
    return record.getMappingQuality() >= minMappingQuality && !record.getReadUnmappedFlag()
            && (!record.getDuplicateReadFlag() || !dropDuplicates) && !record.isSecondaryOrSupplementary();
}
 

private boolean samRecordMeetsQualityRequirements(@NotNull final SAMRecord record) {
    return record.getMappingQuality() >= minMappingQuality && !record.getReadUnmappedFlag() && !record.getDuplicateReadFlag() && !record
            .isSecondaryOrSupplementary();
}
 

private boolean rejectOnMapQuality (final SAMRecord r) {
	if (this.MINIMUM_MAPPING_QUALITY==null) return (false);
	if (r.getMappingQuality() < this.MINIMUM_MAPPING_QUALITY) return (true);
	return (false);
}
 

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 isEligible(@NotNull final SAMRecord record) {
    return record.getMappingQuality() >= MIN_MAPPING_QUALITY && !(record.getReadUnmappedFlag() || record.getDuplicateReadFlag()
            || record.isSecondaryOrSupplementary());
}
 

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

@Override
public boolean test(SAMRecord r) {
    boolean flag=(rejectNonPrimaryReads && r.isSecondaryOrSupplementary()) ||
            (this.mapQuality!=null && r.getMappingQuality() < this.mapQuality);
    return !flag;
}
 

private void testBestFragmentMapqStrategy(final String testName, final int[] firstMapQs, final int[] secondMapQs,
                                          final boolean includeSecondary, final int expectedFirstMapq,
                                          final int expectedSecondMapq) throws Exception {
    final File unmappedSam = File.createTempFile("unmapped.", ".sam");
    unmappedSam.deleteOnExit();
    final SAMFileWriterFactory factory = new SAMFileWriterFactory();
    final SAMFileHeader header = new SAMFileHeader();
    header.setSortOrder(SAMFileHeader.SortOrder.queryname);

    final String readName = "theRead";
    final SAMRecord firstUnmappedRead = new SAMRecord(header);
    firstUnmappedRead.setReadName(readName);
    firstUnmappedRead.setReadString("ACGTACGTACGTACGT");
    firstUnmappedRead.setBaseQualityString("5555555555555555");
    firstUnmappedRead.setReadUnmappedFlag(true);
    firstUnmappedRead.setMateUnmappedFlag(true);
    firstUnmappedRead.setReadPairedFlag(true);
    firstUnmappedRead.setFirstOfPairFlag(true);

    final SAMRecord secondUnmappedRead = new SAMRecord(header);
    secondUnmappedRead.setReadName(readName);
    secondUnmappedRead.setReadString("TCGAACGTTCGAACTG");
    secondUnmappedRead.setBaseQualityString("6666666666666666");
    secondUnmappedRead.setReadUnmappedFlag(true);
    secondUnmappedRead.setMateUnmappedFlag(true);
    secondUnmappedRead.setReadPairedFlag(true);
    secondUnmappedRead.setSecondOfPairFlag(true);

    final SAMFileWriter unmappedWriter = factory.makeSAMWriter(header, false, unmappedSam);
    unmappedWriter.addAlignment(firstUnmappedRead);
    unmappedWriter.addAlignment(secondUnmappedRead);
    unmappedWriter.close();

    final File alignedSam = File.createTempFile("aligned.", ".sam");
    alignedSam.deleteOnExit();

    final String sequence = "chr1";
    // Populate the header with SAMSequenceRecords
    header.setSequenceDictionary(SAMSequenceDictionaryExtractor.extractDictionary(sequenceDict2.toPath()));

    final SAMFileWriter alignedWriter = factory.makeSAMWriter(header, false, alignedSam);

    addAlignmentsForBestFragmentMapqStrategy(alignedWriter, firstUnmappedRead, sequence, firstMapQs);
    addAlignmentsForBestFragmentMapqStrategy(alignedWriter, secondUnmappedRead, sequence, secondMapQs);
    alignedWriter.close();

    final File output = File.createTempFile("testBestFragmentMapqStrategy." + testName, ".sam");
    output.deleteOnExit();
    doMergeAlignment(unmappedSam, Collections.singletonList(alignedSam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true,
            fasta, output,
            SamPairUtil.PairOrientation.FR,
            MergeBamAlignment.PrimaryAlignmentStrategy.BestEndMapq,
            null, includeSecondary, null, null);
    final SamReader reader = SamReaderFactory.makeDefault().open(output);

    int numFirstRecords = 0;
    int numSecondRecords = 0;
    int firstPrimaryMapq = -1;
    int secondPrimaryMapq = -1;
    for (final SAMRecord rec: reader) {
        Assert.assertTrue(rec.getReadPairedFlag());
        if (rec.getFirstOfPairFlag()) ++numFirstRecords;
        else if (rec.getSecondOfPairFlag()) ++ numSecondRecords;
        else Assert.fail("unpossible!");
        if (!rec.getReadUnmappedFlag() && !rec.getNotPrimaryAlignmentFlag()) {
            if (rec.getFirstOfPairFlag()) {
                Assert.assertEquals(firstPrimaryMapq, -1);
                firstPrimaryMapq = rec.getMappingQuality();
            } else {
                Assert.assertEquals(secondPrimaryMapq, -1);
                secondPrimaryMapq = rec.getMappingQuality();
            }
        } else if (rec.getNotPrimaryAlignmentFlag()) {
            Assert.assertTrue(rec.getMateUnmappedFlag());
        }
    }
    reader.close();
    Assert.assertEquals(firstPrimaryMapq, expectedFirstMapq);
    Assert.assertEquals(secondPrimaryMapq, expectedSecondMapq);
    if (!includeSecondary) {
        Assert.assertEquals(numFirstRecords, 1);
        Assert.assertEquals(numSecondRecords, 1);
    } else {
        // If no alignments for an end, there will be a single unmapped record
        Assert.assertEquals(numFirstRecords, Math.max(1, firstMapQs.length));
        Assert.assertEquals(numSecondRecords, Math.max(1, secondMapQs.length));
    }
}
 

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

/**
 * Read a single fragment read from a file, and create one or more aligned records for the read pair based on
 * the lists, merge with the original read, and assert expected results.
 * @param description
 * @param hitSpecs List that describes the aligned SAMRecords to create.
 * @param expectedPrimaryHitIndex Expected value for the HI tag in the primary alignment in the merged output.
 * @param expectedNumReads Expected number of SAMRecords in the merged output.
 * @param expectedPrimaryMapq Mapq of both ends of primary alignment in the merged output.
 * @throws Exception
 */
@Test(dataProvider = "testFragmentMultiHitWithFilteringTestCases")
public void testFragmentMultiHitWithFiltering(final String description, final List<HitSpec> hitSpecs,
                                            final Integer expectedPrimaryHitIndex, final int expectedNumReads,
                                            final int expectedPrimaryMapq) throws Exception {

    // Create the aligned file by copying bases, quals, readname from the unmapped read, and conforming to each HitSpec.
    final File unmappedSam = new File(TEST_DATA_DIR, "multihit.filter.fragment.unmapped.sam");
    final SAMRecordIterator unmappedSamFileIterator = SamReaderFactory.makeDefault().open(unmappedSam).iterator();
    final SAMRecord unmappedRec = unmappedSamFileIterator.next();
    unmappedSamFileIterator.close();
    final File alignedSam = File.createTempFile("aligned.", ".sam");
    alignedSam.deleteOnExit();
    final SAMFileHeader alignedHeader = new SAMFileHeader();
    alignedHeader.setSortOrder(SAMFileHeader.SortOrder.queryname);
    alignedHeader.setSequenceDictionary(SamReaderFactory.makeDefault().getFileHeader(sequenceDict).getSequenceDictionary());
    final SAMFileWriter alignedWriter = new SAMFileWriterFactory().makeSAMWriter(alignedHeader, true, alignedSam);
    for (int i = 0; i < hitSpecs.size(); ++i) {
        final HitSpec hitSpec = hitSpecs.get(i);
        final SAMRecord mappedRec = makeRead(alignedHeader, unmappedRec, hitSpec, i);
        if (mappedRec != null) {
            alignedWriter.addAlignment(mappedRec);
        }
    }
    alignedWriter.close();

    // Merge aligned file with original unmapped file.
    final File mergedSam = File.createTempFile("merged.", ".sam");
    mergedSam.deleteOnExit();

    doMergeAlignment(unmappedSam, Collections.singletonList(alignedSam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            false, fasta, mergedSam,
            null, null, null, null, null, null);

    assertSamValid(mergedSam);

    // Tally metrics and check for agreement with expected.
    final SamReader mergedReader = SamReaderFactory.makeDefault().open(mergedSam);
    int numReads = 0;
    Integer primaryHitIndex = null;
    int primaryMapq = 0;
    for (final SAMRecord rec : mergedReader) {
        ++numReads;
        if (!rec.getNotPrimaryAlignmentFlag()  && !rec.getReadUnmappedFlag()) {
            final Integer hitIndex = rec.getIntegerAttribute(SAMTag.HI.name());
            final int newHitIndex = (hitIndex == null? -1: hitIndex);
            Assert.assertNull(primaryHitIndex);
            primaryHitIndex = newHitIndex;
            primaryMapq = rec.getMappingQuality();
        }
    }
    Assert.assertEquals(numReads, expectedNumReads);
    Assert.assertEquals(primaryHitIndex, expectedPrimaryHitIndex);
    Assert.assertEquals(primaryMapq, expectedPrimaryMapq);
}
 

private boolean isHighQualityMapping(final SAMRecord record) {
    return !record.getReadFailsVendorQualityCheckFlag() &&
            record.getMappingQuality() >= MAPPING_QUALITY_THRESHOLD;
}