Java Code Examples for htsjdk.samtools.util.SequenceUtil#reverseQualities()

/**
 * Returns input read with alignment-related info cleared
 */
private static GATKRead clearReadAlignment(final GATKRead read, final SAMFileHeader header) {
    final GATKRead newRead = new SAMRecordToGATKReadAdapter(new SAMRecord(header));
    newRead.setName(read.getName());
    newRead.setBases(read.getBases());
    newRead.setBaseQualities(read.getBaseQualities());
    if (read.isReverseStrand()) {
        SequenceUtil.reverseComplement(newRead.getBases());
        SequenceUtil.reverseQualities(newRead.getBaseQualities());
    }
    newRead.setIsUnmapped();
    newRead.setIsPaired(read.isPaired());
    if (read.isPaired()) {
        newRead.setMateIsUnmapped();
        if (read.isFirstOfPair()) {
            newRead.setIsFirstOfPair();
        } else if (read.isSecondOfPair()) {
            newRead.setIsSecondOfPair();
        }
    }
    final String readGroup = read.getReadGroup();
    if (readGroup != null) {
        newRead.setAttribute(SAMTag.RG.name(), readGroup);
    }
    return newRead;
}
 

void reset(EvidenceRecord evidenceRecord) {
	this.evidenceRecord = evidenceRecord;

	negative = "-".equals(evidenceRecord.Strand);

	baseBuf.clear();
	scoreBuf.clear();
	if (evidenceRecord.Strand.equals("+")) {
		baseBuf.put(evidenceRecord.Sequence.getBytes());
		scoreBuf.put(SAMUtils.fastqToPhred(evidenceRecord.Scores));
	} else {
		byte[] bytes = evidenceRecord.Sequence.getBytes();
		SequenceUtil.reverseComplement(bytes);
		baseBuf.put(bytes);
		bytes = SAMUtils.fastqToPhred(evidenceRecord.Scores);
		SequenceUtil.reverseQualities(bytes);
		scoreBuf.put(bytes);
	}
	baseBuf.flip();
	scoreBuf.flip();

	firstHalf.clear();
	secondHalf.clear();
}
 

public FastqRead( final GATKRead read, final boolean includeMappingLocation ) {
    this.header = composeHeaderLine(read, includeMappingLocation);
    this.bases = read.getBases();
    this.quals = read.getBaseQualities();
    if (!read.isUnmapped() && read.isReverseStrand()) {
        SequenceUtil.reverseComplement(this.bases);
        SequenceUtil.reverseQualities(this.quals);
    }
}
 

public void acceptRecord(final SAMRecordAndReference args) {
	mappedRecordCount++;

	final SAMRecord samRecord = args.getSamRecord();
	final ReferenceSequence referenceSequence = args.getReferenceSequence();

	final byte[] readBases = samRecord.getReadBases();
	final byte[] readQualities = samRecord.getBaseQualities();
	final byte[] refBases = referenceSequence.getBases();

	if (samRecord.getReadLength() < minReadLength) {
		smallReadCount++;
		return;
	} else if (SequenceUtil.countMismatches(samRecord, refBases, true) > Math.round(samRecord.getReadLength() * maxMismatchRate)) {
		mismatchCount++;
		return;
	}

	// We only record non-CpG C sites if there was at least one CpG in the read, keep track of
	// the values for this record and then apply to the global totals if valid
	int recordCpgs = 0;

	for (final AlignmentBlock alignmentBlock : samRecord.getAlignmentBlocks()) {
		final int blockLength = alignmentBlock.getLength();
		final int refFragmentStart = alignmentBlock.getReferenceStart() - 1;
		final int readFragmentStart = alignmentBlock.getReadStart() - 1;

		final byte[] refFragment = getFragment(refBases, refFragmentStart, blockLength);
		final byte[] readFragment = getFragment(readBases, readFragmentStart, blockLength);
		final byte[] readQualityFragment = getFragment(readQualities, readFragmentStart, blockLength);

		if (samRecord.getReadNegativeStrandFlag()) {
			// In the case of a negative strand, reverse (and complement for base arrays) the reference,
			// reads & qualities so that it can be treated as a positive strand for the rest of the process
			SequenceUtil.reverseComplement(refFragment);
			SequenceUtil.reverseComplement(readFragment);
			SequenceUtil.reverseQualities(readQualityFragment);
		}

		for (int i=0; i < blockLength-1; i++) {
			final int curRefIndex = getCurRefIndex(refFragmentStart, blockLength, i, samRecord.getReadNegativeStrandFlag());

			// Look at a 2-base window to see if we're on a CpG site, and if so check for conversion
			// (CG -> TG). We do not consider ourselves to be on a CpG site if we're on the last base of a read
			if ((SequenceUtil.basesEqual(refFragment[i], SequenceUtil.C)) &&
					(SequenceUtil.basesEqual(refFragment[i+1], SequenceUtil.G))) {
				// We want to catch the case where there's a CpG in the reference, even if it is not valid
				// to prevent the C showing up as a non-CpG C down below. Otherwise this could have been all
				// in one if statement
				if (isValidCpg(refFragment, readFragment, readQualityFragment, i)) {
					recordCpgs++;
					final CpgLocation curLocation = new CpgLocation(samRecord.getReferenceName(), curRefIndex);
					cpgTotal.increment(curLocation);
					if (SequenceUtil.isBisulfiteConverted(readFragment[i], refFragment[i])) {
						cpgConverted.increment(curLocation);
					}
				}
				i++;
			} else if (isC(refFragment[i], readFragment[i]) && isAboveCytoQcThreshold(readQualities, i) &&
					SequenceUtil.bisulfiteBasesEqual(false, readFragment[i+1], refFragment[i+1])) {
				// C base in the reference that's not associated with a CpG
				nCytoTotal++;
				if (SequenceUtil.isBisulfiteConverted(readFragment[i], refFragment[i])) {
					nCytoConverted++;
				}
			}
		}
	}

	if (recordCpgs == 0) {
		noCpgCount++;
	}
}
 

/**
 * Builds a SAMRecord from unaligned read data and an alignment.
 * qualsArg can be null.
 * readGroup can be null.
 *
 * Also, BWA does this odd thing, supplying AS and XS tags for unaligned reads.
 * To produce SAM records that are byte-identical to the ones created by the bwa mem command line,
 * call this method with alwaysGenerateASandXS=true.
 */
public static SAMRecord applyAlignment( final String readName, final byte[] basesArg, final byte[] qualsArg,
                                       final String readGroup, final BwaMemAlignment alignment,
                                       final List<String> refNames, final SAMFileHeader header,
                                       final boolean softClipAlts, final boolean alwaysGenerateASandXS ) {
    final SAMRecord samRecord = new SAMRecord(header);
    samRecord.setReadName(readName);
    final int samFlag = alignment.getSamFlag();
    samRecord.setFlags(samFlag);
    if ( alignment.getRefId() >= 0 ) samRecord.setReferenceName(refNames.get(alignment.getRefId()));
    else if ( alignment.getMateRefId() >= 0 ) samRecord.setReferenceName(refNames.get(alignment.getMateRefId()));
    if ( alignment.getRefStart() >= 0 ) samRecord.setAlignmentStart(alignment.getRefStart()+1);
    else if ( alignment.getMateRefStart() >= 0 ) samRecord.setAlignmentStart(alignment.getMateRefStart()+1);
    if ( alignment.getMapQual() >= 0 ) samRecord.setMappingQuality(alignment.getMapQual());
    byte[] bases = basesArg;
    byte[] quals = qualsArg == null ? new byte[0] : qualsArg;
    if ( SAMFlag.READ_REVERSE_STRAND.isSet(samFlag) &&
            SAMFlag.SECONDARY_ALIGNMENT.isUnset(samFlag) ) {
        bases = BaseUtils.simpleReverseComplement(bases);
        quals = Arrays.copyOf(quals, quals.length);
        SequenceUtil.reverseQualities(quals);
    }
    if ( alignment.getCigar() != null && !alignment.getCigar().isEmpty() ) {
        final Cigar cigar = TextCigarCodec.decode(alignment.getCigar());
        Cigar tmpCigar = cigar;
        if ( !softClipAlts && SAMFlag.SUPPLEMENTARY_ALIGNMENT.isSet(samFlag) ) {
            if ( tmpCigar.getFirstCigarElement().getOperator() == CigarOperator.S ||
                    tmpCigar.getLastCigarElement().getOperator() == CigarOperator.S ) {
                tmpCigar = new Cigar();
                for ( final CigarElement ele : cigar ) {
                    if ( ele.getOperator() == CigarOperator.S ) {
                        tmpCigar.add(new CigarElement(ele.getLength(), CigarOperator.H));
                    }
                    else {
                        tmpCigar.add(ele);
                    }
                }
            }
            bases = Arrays.copyOfRange(bases, alignment.getSeqStart(), alignment.getSeqEnd());
            if ( quals.length != 0 )
                quals = Arrays.copyOfRange(quals, alignment.getSeqStart(), alignment.getSeqEnd());
        }
        samRecord.setCigar(tmpCigar);
        samRecord.setAttribute("NM", alignment.getNMismatches());
        samRecord.setAttribute("AS", alignment.getAlignerScore());
        samRecord.setAttribute("XS", alignment.getSuboptimalScore());
        samRecord.setAttribute("MD", alignment.getMDTag());
        samRecord.setAttribute("XA", alignment.getXATag());
    }
    else if ( alwaysGenerateASandXS ) {
        samRecord.setAttribute("AS", 0);
        samRecord.setAttribute("XS", 0);
    }
    if ( SAMFlag.READ_PAIRED.isSet(samFlag) ) {
        if ( alignment.getMateRefId() >= 0 ) samRecord.setMateReferenceName(refNames.get(alignment.getMateRefId()));
        else if ( alignment.getRefId() >= 0 ) samRecord.setMateReferenceName(refNames.get(alignment.getRefId()));
        if ( alignment.getMateRefStart() >= 0 ) samRecord.setMateAlignmentStart(alignment.getMateRefStart() + 1);
        else if ( alignment.getRefStart() >= 0 ) samRecord.setMateAlignmentStart(alignment.getRefStart() + 1);
        if ( alignment.getTemplateLen() != 0 ) samRecord.setInferredInsertSize(alignment.getTemplateLen());
    }
    if ( SAMFlag.SECONDARY_ALIGNMENT.isUnset(samFlag) ) {
        samRecord.setReadBases(bases);
        samRecord.setBaseQualities(quals);
    } else {
        samRecord.setReadBases(SAMRecord.NULL_SEQUENCE);
        samRecord.setBaseQualities(SAMRecord.NULL_QUALS);
    }
    //TODO: there ought to be a way to indicate a set of tag names that ought to be copied -- we're just doing RG
    if ( readGroup != null ) samRecord.setAttribute(SAMTag.RG.name(), readGroup);
    return samRecord;
}