Java SAMUtils类代码示例

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
/**
 * Setters and Accessors for base insertion and base deletion quality scores
 */
public void setBaseQualities( final byte[] quals, final EventType errorModel ) {
    switch( errorModel ) {
        case BASE_SUBSTITUTION:
            setBaseQualities(quals);
            break;
        case BASE_INSERTION:
            setAttribute( GATKSAMRecord.BQSR_BASE_INSERTION_QUALITIES, quals == null ? null : SAMUtils.phredToFastq(quals) );
            break;
        case BASE_DELETION:
            setAttribute( GATKSAMRecord.BQSR_BASE_DELETION_QUALITIES, quals == null ? null : SAMUtils.phredToFastq(quals) );
            break;
        default:
            throw new ReviewedGATKException("Unrecognized Base Recalibration type: " + errorModel );
    }
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
private void initializeCachedData() {
    for ( int i = 0; i < 256; i++ )
        for ( int j = 0; j < 256; j++ )
            for ( int q = 0; q <= SAMUtils.MAX_PHRED_SCORE; q++ ) {
                EPSILONS[i][j][q] = 1.0;
            }

    for ( char b1 : "ACGTacgt".toCharArray() ) {
        for ( char b2 : "ACGTacgt".toCharArray() ) {
            for ( int q = 0; q <= SAMUtils.MAX_PHRED_SCORE; q++ ) {
                double qual = qual2prob[q < minBaseQual ? minBaseQual : q];
                double e = Character.toLowerCase(b1) == Character.toLowerCase(b2) ? 1 - qual : qual * EM;
                EPSILONS[(byte)b1][(byte)b2][q] = e;
            }
        }
    }
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
/**
 * Helper function that returns the phred-scaled base quality score we should use for calculating
 * likelihoods for a pileup element.  May return 0 to indicate that the observation is bad, and may
 * cap the quality score by the mapping quality of the read itself.
 *
 * @param p                           Pileup element
 * @param ignoreBadBases              Should we ignore bad bases?
 * @param capBaseQualsAtMappingQual   Should we cap the base qualities at the mapping quality of the read?
 * @param minBaseQual                 Minimum allowed base quality
 * @return the actual base quality to use
 */
private static byte qualToUse(PileupElement p, boolean ignoreBadBases, boolean capBaseQualsAtMappingQual, int minBaseQual) {
    if ( ignoreBadBases && !BaseUtils.isRegularBase( p.getBase() ) )
        return 0;

    byte qual = p.getQual();

    if ( qual > SAMUtils.MAX_PHRED_SCORE )
        throw new UserException.MisencodedBAM(p.getRead(), "we encountered an extremely high quality score (" + (int)qual + ")");
    if ( capBaseQualsAtMappingQual )
        qual = (byte) Math.min( 0xff & qual, p.getMappingQual());
    if ( (int)qual < minBaseQual )
        qual = (byte)0;

    return qual;
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
private void initializeCachedData() {
	for (int i = 0; i < 256; i++)
		for (int j = 0; j < 256; j++)
			for (int q = 0; q <= SAMUtils.MAX_PHRED_SCORE; q++) {
				EPSILONS[i][j][q] = 1.0;
			}

	for (char b1 : "ACGTacgt".toCharArray()) {
		for (char b2 : "ACGTacgt".toCharArray()) {
			for (int q = 0; q <= SAMUtils.MAX_PHRED_SCORE; q++) {
				double qual = qual2prob[q < minBaseQual ? minBaseQual : q];
				double e = Character.toLowerCase(b1) == Character
						.toLowerCase(b2) ? 1 - qual : qual * EM;
				EPSILONS[(byte) b1][(byte) b2][q] = e;
			}
		}
	}
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
@Override
public GATKRead apply(final GATKRead read) {
    // get the qualities
    final byte[] quals = read.getBaseQualities();
    try {
        // transform them in place
        for (int i = 0; i < quals.length; ++i) {
            if (quals[i] < MIN_SOLEXA_BEFORE_CONVERSION) {
                // throw an exception if they are missencoded
                throwException();
            }
            // convert to Solexa
            quals[i] = SOLEXA_QUALITY_CONVERTER
                    .solexaCharToPhredBinary((byte) SAMUtils.phredToFastq(quals[i]));
        }
    } catch (IndexOutOfBoundsException e) {
        // if there is an index exception, that means that the qualities are not correctly encoded
        throwException();
    }
    // set the base qualities
    read.setBaseQualities(quals);
    return read;
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
@Override
protected int doWork() {
    final File output =
            OUTPUT == null
                    ? new File(INPUT.getParentFile(), getOutputFileName(INPUT))
                    : OUTPUT;

    IOUtil.assertFileIsWritable(output);
    final String hashText = SAMUtils.calculateReadGroupRecordChecksum(INPUT, REFERENCE_SEQUENCE);

    try {
        final FileWriter outputWriter = new FileWriter(output);
        outputWriter.write(hashText);
        outputWriter.close();
    } catch (final IOException ioe) {
        throw new PicardException(
                "Could not write the computed hash (" + hashText + ") to the output file: " + ioe.getMessage(), ioe);
    }
    return 0;
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
/**
 *
 * @param querynameOrderIterator
 * @param primaryAlignmentSelectionStrategy Algorithm for selecting primary alignment when it is not clear from
 *                                          the input what should be primary.
 */
MultiHitAlignedReadIterator(final CloseableIterator<SAMRecord> querynameOrderIterator,
                            final PrimaryAlignmentSelectionStrategy primaryAlignmentSelectionStrategy) {
    this.primaryAlignmentSelectionStrategy = primaryAlignmentSelectionStrategy;
    peekIterator = new PeekableIterator<SAMRecord>(new FilteringSamIterator(querynameOrderIterator,
            new SamRecordFilter() {
                // Filter unmapped reads.
                public boolean filterOut(final SAMRecord record) {
                    return record.getReadUnmappedFlag() || SAMUtils.cigarMapsNoBasesToRef(record.getCigar());
                }
                public boolean filterOut(final SAMRecord first, final SAMRecord second) {
                    return ((first.getReadUnmappedFlag() || SAMUtils.cigarMapsNoBasesToRef(first.getCigar()))
                            && (second.getReadUnmappedFlag() || SAMUtils.cigarMapsNoBasesToRef(second.getCigar())));
                }
            }));


    advance();
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
@Test
public void testUsingOriginalQualities() throws Exception {
    final MeanQualityByCycleSpark.HistogramGenerator hg = new MeanQualityByCycleSpark.HistogramGenerator(true);
    Assert.assertEquals(hg.useOriginalQualities, true);

    GATKRead read1 = ArtificialReadUtils.createArtificialRead("aa".getBytes(), new byte[]{50, 50}, "2M");
    hg.addRead(read1);
    assertEqualsLongArray(hg.firstReadCountsByCycle, new long[0]);
    assertEqualsDoubleArray(hg.firstReadTotalsByCycle, new double[0], 1e-05);
    assertEqualsLongArray(hg.secondReadCountsByCycle, new long[0]);
    assertEqualsDoubleArray(hg.secondReadTotalsByCycle, new double[0], 1e-05);

    GATKRead read2 = ArtificialReadUtils.createArtificialRead("aa".getBytes(), new byte[]{50, 50}, "2M");
    read2.setAttribute(SAMTag.OQ.name(), SAMUtils.phredToFastq(new byte[]{30, 40}));
    hg.addRead(read2);

    assertEqualsLongArray(hg.firstReadCountsByCycle, new long[]{0, 1, 1});
    assertEqualsDoubleArray(hg.firstReadTotalsByCycle, new double[]{0, 30, 40}, 1e-05);
    assertEqualsLongArray(hg.secondReadCountsByCycle, new long[]{0, 0, 0});
    assertEqualsDoubleArray(hg.secondReadTotalsByCycle, new double[]{0, 0, 0}, 1e-05);
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
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();
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
SAMRecord firstSAMRecord(SAMFileHeader header) {
	SAMRecord r = new SAMRecord(header);
	r.setReadName(evidenceRecord.getReadName());
	r.setReferenceName(evidenceRecord.Chromosome);
	r.setAlignmentStart(Integer.valueOf(evidenceRecord.OffsetInReference) + 1);
	r.setMappingQuality(Integer.valueOf(evidenceRecord.ScoreAllele0));
	r.setReadPairedFlag(true);
	r.setReadUnmappedFlag(false);
	r.setReadNegativeStrandFlag(negative);
	r.setFirstOfPairFlag(evidenceRecord.side == 0);
	r.setSecondOfPairFlag(!r.getFirstOfPairFlag());

	r.setCigar(new Cigar(Utils.toCigarOperatorList(firstHalf.samCigarElements)));

	r.setReadBases(Utils.toByteArray(firstHalf.readBasesBuf));
	r.setBaseQualities(Utils.toByteArray(firstHalf.readScoresBuf));

	r.setAttribute("GC", Utils.toString(firstHalf.gcList));
	r.setAttribute("GS", Utils.toString(firstHalf.gsBuf));
	r.setAttribute("GQ", SAMUtils.phredToFastq(Utils.toByteArray(firstHalf.gqBuf)));

	return r;
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
SAMRecord secondSAMRecord(SAMFileHeader header) {
	SAMRecord r = new SAMRecord(header);
	r.setReadName(evidenceRecord.getReadName());
	r.setReferenceName(evidenceRecord.Chromosome);
	r.setAlignmentStart(Integer.valueOf(evidenceRecord.MateOffsetInReference) + 1);
	r.setMappingQuality(Integer.valueOf(evidenceRecord.ScoreAllele0));
	r.setReadPairedFlag(true);
	r.setReadUnmappedFlag(false);
	r.setReadNegativeStrandFlag(negative);
	r.setFirstOfPairFlag(evidenceRecord.side == 1);
	r.setSecondOfPairFlag(!r.getFirstOfPairFlag());

	r.setCigar(new Cigar(Utils.toCigarOperatorList(secondHalf.samCigarElements)));

	r.setReadBases(Utils.toByteArray(secondHalf.readBasesBuf));
	r.setBaseQualities(Utils.toByteArray(secondHalf.readScoresBuf));

	r.setAttribute("GC", Utils.toString(secondHalf.gcList));
	r.setAttribute("GS", Utils.toString(secondHalf.gsBuf));
	r.setAttribute("GQ", SAMUtils.phredToFastq(Utils.toByteArray(secondHalf.gqBuf)));

	return r;
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
private void analyseSamRecord(final SAMRecord rec) {
if(rec.getReadUnmappedFlag()) return;
if(rec.getReadFailsVendorQualityCheckFlag()) return;
if(rec.isSecondaryOrSupplementary()) return;
if(rec.getDuplicateReadFlag()) return;

final List<SAMRecord> others= SAMUtils.getOtherCanonicalAlignments(rec);
if(others.isEmpty()) return;
String sample=this.defaultSampleName;
final SAMReadGroupRecord g=rec.getReadGroup();
if(g!=null) {
	final String sa = g.getSample();
	if(sa!=null) sample=sa;
	}

IntervalTreeMap<Set<Arc>> database = this.sample2database.get(sample);
if(database==null) {
	database=new IntervalTreeMap<Set<Arc>>();
	this.sample2database.put(sample, database);
	}

for(final SAMRecord other:others)
	{
	analyseSamPair(database,rec,other);
	}
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
private void fillCache() {
    // cache Q(j,k) = log10(j/2N*(1-ek) + (2N-j)/2N*ek) for j = 0:2N

    logMismatchProbabilityArray = new double[1+numChromosomes][1+SAMUtils.MAX_PHRED_SCORE];
    for (int i=0; i <= numChromosomes; i++) {
        for (int j=0; j <= SAMUtils.MAX_PHRED_SCORE; j++) {
            double phi = (double)i/numChromosomes;
            logMismatchProbabilityArray[i][j] = Math.log10(phi * (1.0-qualVec[j]) + qualVec[j]/3.0 * (1.0-phi));
        }
    }
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
private void writeDebugLikelihoods(final GATKSAMRecord processedRead, final Haplotype haplotype, final double log10l) {
    likelihoodsStream.printf("%s %s %s %s %s %s %f%n",
            haplotype.getBaseString(),
            new String(processedRead.getReadBases()),
            SAMUtils.phredToFastq(processedRead.getBaseQualities()),
            SAMUtils.phredToFastq(processedRead.getBaseInsertionQualities()),
            SAMUtils.phredToFastq(processedRead.getBaseDeletionQualities()),
            SAMUtils.phredToFastq(constantGCP),
            log10l);
}
 

import htsjdk.samtools.SAMUtils; //导入依赖的package包/类
private void writeDebugLikelihoods(final GATKRead processedRead, final Haplotype haplotype, final double log10l){
    // Note: the precision of log10l in the debug output is only ~6 digits (ie., not all digits are necessarily printed)
    likelihoodsStream.printf("%s %s %s %s %s %s %f%n",
            haplotype.getBaseString(),
            new String(processedRead.getBases()),
            SAMUtils.phredToFastq(processedRead.getBaseQualities()),
            SAMUtils.phredToFastq(ReadUtils.getBaseInsertionQualities(processedRead)),
            SAMUtils.phredToFastq(ReadUtils.getBaseDeletionQualities(processedRead)),
            SAMUtils.phredToFastq(constantGCP),
            log10l);
}