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

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

public void addRead (SAMRecord r) {
	// skip secondary of supplemental reads.
	if (r.isSecondaryOrSupplementary()) {
		return;
	}
	
	boolean isDupe = r.getDuplicateReadFlag();
	int mapQuality = r.getMappingQuality();
	boolean unmapped = r.getReadUnmappedFlag();
	if (histogram!=null) {
		histogram.increment(mapQuality);
	}
	
	totalReads++;
	if (!unmapped) {
		mappedReads++;
		if (mapQuality >= this.mapQuality) {
			hqMappedReads++;
			if (!isDupe) {
				hqMappedReadsNoPCRDupes++;					
			}
		}
	}
}
 

BaseDistributionMetricCollection gatherBaseQualities (final File input, final String tag) {
	ProgressLogger pl = new ProgressLogger(this.log);
	SamReader inputSam = SamReaderFactory.makeDefault().open(input);

	BaseDistributionMetricCollection c = new BaseDistributionMetricCollection();
	// MAIN_LOOP:
	for (final SAMRecord samRecord : inputSam) {
		pl.record(samRecord);
		if (samRecord.isSecondaryOrSupplementary()) continue;
		String b = samRecord.getStringAttribute(tag);
		if (b==null) continue;

		byte [] bases = b.getBytes();
		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);

}
 

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

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


}
 

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

SameCoordReadSet getSameCoordReads(PeekIterator<SAMRecord> it, String fileName) throws Exception {
  SameCoordReadSet ret = null;
  try {
    SAMRecord record;
    while (it.hasNext()) {
      record = it.peek();
      if (record.isSecondaryOrSupplementary() ||
          (options.ignoreUnmappedReads && record.getReadUnmappedFlag())) {
        it.next();
        continue;
      }
      if (ret != null) {
        if (record.getAlignmentStart() != ret.coord || !record.getReferenceName().equals(ret.reference)) {
          break;
        }
      } else {
        ret = new SameCoordReadSet();
        ret.map = Maps.newHashMap();
        ret.coord = record.getAlignmentStart();
        ret.reference = record.getReferenceName();
      }
      ret.map.put(record.getReadName(), record);
      it.next();
    }
  } catch (Exception ex) {
    throw new Exception("Error reading from " + fileName + "\n" + ex.getMessage());
  }
  return ret;
}
 

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

public PrimaryAlignmentKey(final SAMRecord rec) {
    if (rec.isSecondaryOrSupplementary()) {
        throw new IllegalArgumentException("PrimaryAligmentKey cannot be a secondary or suplementary alignment");
    }
    this.pairStatus = rec.getReadPairedFlag() ?
            (rec.getSecondOfPairFlag() ? PairStatus.SECOND : PairStatus.FIRST) :
            PairStatus.UNPAIRED;
    this.readName = rec.getReadName();
}
 

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

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
public void acceptRecord(final SAMRecord rec, final ReferenceSequence ref) {
    if (!rec.isSecondaryOrSupplementary()) {
        super.acceptRecord(rec, ref);
    }
}
 

protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    ValidationUtil.assertWritable(OUTPUT, OUTPUT_BY_READGROUP);

    final boolean sanitizing = SANITIZE;
    final SamReader in = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).validationStringency(VALIDATION_STRINGENCY).open(INPUT);
    final SAMFileHeader inHeader = in.getFileHeader();
    ValidationUtil.validateHeaderOverrides(inHeader, SAMPLE_ALIAS, LIBRARY_NAME);

    ////////////////////////////////////////////////////////////////////////////
    // Build the output writer with an appropriate header based on the options
    ////////////////////////////////////////////////////////////////////////////
    final boolean presorted = isPresorted(inHeader, SORT_ORDER, sanitizing);
    if (SAMPLE_ALIAS != null) overwriteSample(inHeader.getReadGroups(), SAMPLE_ALIAS);
    if (LIBRARY_NAME != null) overwriteLibrary(inHeader.getReadGroups(), LIBRARY_NAME);
    final SAMFileHeader singleOutHeader = createOutHeader(inHeader, SORT_ORDER, REMOVE_ALIGNMENT_INFORMATION);
    inHeader.getReadGroups().forEach(readGroup -> singleOutHeader.addReadGroup(readGroup));

    final Map<String, File> outputMap;
    final Map<String, SAMFileHeader> headerMap;
    if (OUTPUT_BY_READGROUP) {
        if (inHeader.getReadGroups().isEmpty()) {
            throw new PicardException(INPUT + " does not contain Read Groups");
        }

        final String defaultExtension;
        if (OUTPUT_BY_READGROUP_FILE_FORMAT == FileType.dynamic) {
            defaultExtension = getDefaultExtension(INPUT.toString());
        } else {
            defaultExtension = "." + OUTPUT_BY_READGROUP_FILE_FORMAT.toString();
        }

        outputMap = createOutputMap(OUTPUT_MAP, OUTPUT, defaultExtension, inHeader.getReadGroups());
        ValidationUtil.assertAllReadGroupsMapped(outputMap, inHeader.getReadGroups());
        headerMap = createHeaderMap(inHeader, SORT_ORDER, REMOVE_ALIGNMENT_INFORMATION);
    } else {
        outputMap = null;
        headerMap = null;
    }

    if (RESTORE_HARDCLIPS && !REMOVE_ALIGNMENT_INFORMATION) {
        throw new PicardException("Cannot revert sam file when RESTORE_HARDCLIPS is true and REMOVE_ALIGNMENT_INFORMATION is false.");
    }

    final SAMFileWriterFactory factory = new SAMFileWriterFactory();
    final RevertSamWriter out = new RevertSamWriter(OUTPUT_BY_READGROUP, headerMap, outputMap, singleOutHeader, OUTPUT, presorted, factory, REFERENCE_SEQUENCE);

    ////////////////////////////////////////////////////////////////////////////
    // Build a sorting collection to use if we are sanitizing
    ////////////////////////////////////////////////////////////////////////////
    final RevertSamSorter sorter;
    if (sanitizing)
        sorter = new RevertSamSorter(OUTPUT_BY_READGROUP, headerMap, singleOutHeader, MAX_RECORDS_IN_RAM);
    else sorter = null;

    final ProgressLogger progress = new ProgressLogger(log, 1000000, "Reverted");
    for (final SAMRecord rec : in) {
        // Weed out non-primary and supplemental read as we don't want duplicates in the reverted file!
        if (rec.isSecondaryOrSupplementary()) continue;

        // log the progress before you revert because otherwise the "last read position" might not be accurate
        progress.record(rec);

        // Actually do the reverting of the remaining records
        revertSamRecord(rec);

        if (sanitizing) sorter.add(rec);
        else out.addAlignment(rec);
    }
    CloserUtil.close(in);

    ////////////////////////////////////////////////////////////////////////////
    // Now if we're sanitizing, clean up the records and write them to the output
    ////////////////////////////////////////////////////////////////////////////
    if (!sanitizing) {
        out.close();
    } else {
        final Map<SAMReadGroupRecord, FastqQualityFormat> readGroupToFormat;
        try {
            readGroupToFormat = createReadGroupFormatMap(inHeader, REFERENCE_SEQUENCE, VALIDATION_STRINGENCY, INPUT, RESTORE_ORIGINAL_QUALITIES);
        } catch (final PicardException e) {
            log.error(e.getMessage());
            return -1;
        }

        final long[] sanitizeResults = sanitize(readGroupToFormat, sorter, out);
        final long discarded = sanitizeResults[0];
        final long total = sanitizeResults[1];
        out.close();

        final double discardRate = discarded / (double) total;
        final NumberFormat fmt = new DecimalFormat("0.000%");
        log.info("Discarded " + discarded + " out of " + total + " (" + fmt.format(discardRate) + ") reads in order to sanitize output.");

        if (discardRate > MAX_DISCARD_FRACTION) {
            throw new PicardException("Discarded " + fmt.format(discardRate) + " which is above MAX_DISCARD_FRACTION of " + fmt.format(MAX_DISCARD_FRACTION));
        }
    }

    return 0;
}
 

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

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

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

public void loadReads(File[] bams) throws IOException{

int count = 0;
int numOp = 0;

for(File bam : bams){
    HLA.log.appendln("Loading reads from:\t" + bam.getName());
    Object2IntOpenHashMap<String> readLoadingSet = new Object2IntOpenHashMap<String>();
    readLoadingSet.defaultReturnValue(0);
    
    final SamReader reader = SamReaderFactory.makeDefault().open(bam);
    
    //Kourami bam checker added
    if(!checkHeader(reader.getFileHeader())){
	HLA.log.appendln("Unexpected BAM :\t"+ bam.getName() 
			+"\nThe input BAM MUST be aligned to the set of IMGT/HLA alleles in " + HLA.MSAFILELOC + "\n" 
			+ "Please use the recommended preprocessing steps explained on the github page:\n"
			+ "https://github.com/Kingsford-Group/kourami");
	System.err.println("Unexpected BAM :\t"+ bam.getName() 
			   +"\nThe input BAM MUST be aligned to the set of IMGT/HLA alleles in " + HLA.MSAFILELOC + "\n" 
			   + "Please use the recommended preprocessing steps explained on the github page:\n"
			   + "https://github.com/Kingsford-Group/kourami");
	HLA.log.outToFile();
	System.exit(1);
    }

    for(final SAMRecord samRecord : reader){
	if(count == 0){
	    HLA.READ_LENGTH = samRecord.getReadLength();
	    HLA.log.appendln("Setting HLA.READ_LEGNTH = " + HLA.READ_LENGTH);
	}
	//added checking to process reads matching to HLA-type sequences
	//discarding decoy hits (DQB2, DQA2)
	boolean qc = false;
	if( (samRecord.getReferenceName().indexOf("*") > -1) 
	    && !samRecord.getReadUnmappedFlag() 
	    && !samRecord.isSecondaryOrSupplementary() 
	    && !this.startWIns(samRecord)){
	    count++;
	    if(samRecord.getReadPairedFlag())
		numOp += processRecord(samRecord, readLoadingSet);
	    else
		numOp += processRecordUnpaired(samRecord);
	}
	
	if(HLA.DEBUG && count%10000 == 0)
	    HLA.log.appendln("Processed 10000 reads...");
    }
    reader.close();
}
HLA.log.appendln("Loaded a total of " + count + " mapped reads.");
HLA.log.appendln("A total of " + numOp + " bases");
   }
 

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

public TagOfTagResults<String,String> getResults (final File inputBAM, final String primaryTag, final String secondaryTag, final boolean filterPCRDuplicates, final Integer readQuality) {
	TagOfTagResults<String,String> result = new TagOfTagResults<>();
	SamReader reader = SamReaderFactory.makeDefault().open(inputBAM);
	CloseableIterator<SAMRecord> iter = CustomBAMIterators.getReadsInTagOrder(reader, primaryTag);
	CloserUtil.close(reader);
	String currentTag="";
	Set<String> otherTagCollection=new HashSet<>();

	ProgressLogger progress = new ProgressLogger(log);

	while (iter.hasNext()) {
		SAMRecord r = iter.next();
		progress.record(r);
		// skip reads that don't pass filters.
		boolean discardResult=(filterPCRDuplicates && r.getDuplicateReadFlag()) || r.getMappingQuality()<readQuality || r.isSecondaryOrSupplementary();
		// short circuit if read is below the quality to be considered.
		if (discardResult) continue;

		Object d = r.getAttribute(secondaryTag);
		// short circuit if there's no tag for this read.
		if (d==null) continue;

		String data=null;

		if (d instanceof String)
			data=r.getStringAttribute(secondaryTag);
		else if (d instanceof Integer)
			data=Integer.toString(r.getIntegerAttribute(secondaryTag));


		String newTag = r.getStringAttribute(primaryTag);
		if (newTag==null)
			newTag="";

		// if you see a new tag.
		if (!currentTag.equals(newTag)) {
			// write out tag results, if any.
			if (!currentTag.equals("") && otherTagCollection.size()>0)
				result.addEntries(currentTag, otherTagCollection);
			currentTag=newTag;
			otherTagCollection.clear();
			if (!discardResult) otherTagCollection=addTagToCollection (data,otherTagCollection);

		} else // gather stats
		if (!discardResult) otherTagCollection=addTagToCollection (data,otherTagCollection);
	}
	if (otherTagCollection.size()>0)
		result.addEntries(currentTag, otherTagCollection);

	return (result);
}