Java ReferenceSequence类代码示例

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
private char getReferenceBase(long pos)
		throws ArrayIndexOutOfBoundsException, IOException {
	int refID = (int) (pos >> 32);
	if (refID == cachedRefID) {
		return (char) cachedRefSequence[(int) pos];
	}
	SAMSequenceRecord seq = secramHeader.getSamFileHeader().getSequence(
			refID);
	ReferenceSequence rs = mRsf.getSequence(seq.getSequenceName());

	if (rs == null || rs.length() != seq.getSequenceLength()) {
		System.err.println("Could not find the reference sequence "
				+ seq.getSequenceName() + " in the file");
		throw new IOException("No such sequence in file");
	}

	cachedRefID = refID;
	cachedRefSequence = rs.getBases();

	return (char) cachedRefSequence[(int) pos];
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
/**
 * Get the reference base of the position.
 * @param pos
 *            The absolute position. Chromosome ID is in the higher 32 bits,
 *            and chromosome position is in the lower 32 bits.
 * @throws ArrayIndexOutOfBoundsException
 * @throws IOException
 */
private char getReferenceBase(long pos)
		throws ArrayIndexOutOfBoundsException, IOException {
	int refID = (int) (pos >> 32);
	if (refID == cachedRefID) {
		return (char) cachedRefSequence[(int) pos];
	}
	SAMSequenceRecord seq = mSAMFileHeader.getSequence(refID);
	ReferenceSequence rs = mRsf.getSequence(seq.getSequenceName());

	if (rs == null || rs.length() != seq.getSequenceLength()) {
		System.err.println("Could not find the reference sequence "
				+ seq.getSequenceName() + " in the file");
		throw new IOException("No such sequence in file");
	}

	cachedRefID = refID;
	cachedRefSequence = rs.getBases();

	return (char) cachedRefSequence[(int) pos];
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
public static void main(String[] args) {
	String bam = "./data/HG00115.chrom11.ILLUMINA.bwa.GBR.exome.20130415.bam";
	SamReader reader = SamReaderFactory.makeDefault()
			.validationStringency(ValidationStringency.SILENT)
			.open(new File(bam));
	int length = 0;
	for (final SAMRecord record : reader) {
		if (record.getReadUnmappedFlag())
			continue;
		length += record.getReadLength();
	}

	ReferenceSequenceFile rsf = ReferenceSequenceFileFactory
			.getReferenceSequenceFile(new File("data/hs37d5.fa"));
	ReferenceSequence rs = rsf.getSequence("11");

	System.out.println(length * 1.0 / rs.length());
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
private void checkOutputGCContent(final ReferenceFileSource reference, final SimpleInterval interval, final double observed) {
    final ReferenceSequence sequence = reference.queryAndPrefetch(interval);
    int total = 0;
    int gc = 0;
    for (final byte base : sequence.getBases()) {
        switch (Character.toLowerCase(base)) {
            case 'g':
            case 'c':
                gc++; total++; break;
            case 'a':
            case 't':
            case 'u':
                total++; break;
            default:
        }
    }
    final double expectedGCContent = total == 0 ? Double.NaN : ((double) gc / (double) total);
    if (Double.isNaN(expectedGCContent)) {
        Assert.assertTrue(Double.isNaN(observed));
    } else {
        Assert.assertEquals(observed, expectedGCContent, 0.0001);
    }
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
/**
 * Read all the sequences from the given reference file, and convert into SAMSequenceRecords
 * @param referenceFile fasta or fasta.gz
 * @return SAMSequenceRecords containing info from the fasta, plus from cmd-line arguments.
 */
@Deprecated
public SAMSequenceDictionary makeSequenceDictionary(final File referenceFile) {
    final ReferenceSequenceFile refSeqFile =
            ReferenceSequenceFileFactory.getReferenceSequenceFile(referenceFile, TRUNCATE_NAMES_AT_WHITESPACE);
    ReferenceSequence refSeq;
    final List<SAMSequenceRecord> ret = new ArrayList<>();
    final Set<String> sequenceNames = new HashSet<>();
    for (int numSequences = 0; numSequences < NUM_SEQUENCES && (refSeq = refSeqFile.nextSequence()) != null; ++numSequences) {
        if (sequenceNames.contains(refSeq.getName())) {
            throw new PicardException("Sequence name appears more than once in reference: " + refSeq.getName());
        }
        sequenceNames.add(refSeq.getName());
        ret.add(makeSequenceRecord(refSeq));
    }
    return new SAMSequenceDictionary(ret);
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
/**
 * Create one SAMSequenceRecord from a single fasta sequence
 */
private SAMSequenceRecord makeSequenceRecord(final ReferenceSequence refSeq) {
    final SAMSequenceRecord ret = new SAMSequenceRecord(refSeq.getName(), refSeq.length());

    // Compute MD5 of upcased bases
    final byte[] bases = refSeq.getBases();
    for (int i = 0; i < bases.length; ++i) {
            bases[i] = StringUtil.toUpperCase(bases[i]);
        }

    ret.setAttribute(SAMSequenceRecord.MD5_TAG, md5Hash(bases));
    if (GENOME_ASSEMBLY != null) {
        ret.setAttribute(SAMSequenceRecord.ASSEMBLY_TAG, GENOME_ASSEMBLY);
    }
    ret.setAttribute(SAMSequenceRecord.URI_TAG, URI);
    if (SPECIES != null) {
            ret.setAttribute(SAMSequenceRecord.SPECIES_TAG, SPECIES);
        }
    return ret;
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
public static int[] calculateRefWindowsByGc(final int windows, final File referenceSequence, final int windowSize) {
    final ReferenceSequenceFile refFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(referenceSequence);
    ReferenceSequence ref;

    final int [] windowsByGc = new int [windows];

    while ((ref = refFile.nextSequence()) != null) {
        final byte[] refBases = ref.getBases();
        StringUtil.toUpperCase(refBases);
        final int refLength = refBases.length;
        final int lastWindowStart = refLength - windowSize;

        final CalculateGcState state = new GcBiasUtils().new CalculateGcState();

        for (int i = 1; i < lastWindowStart; ++i) {
            final int windowEnd = i + windowSize;
            final int gcBin = calculateGc(refBases, i, windowEnd, state);
            if (gcBin != -1) windowsByGc[gcBin]++;
        }
    }

    return windowsByGc;
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
/**
 * This loads a FASTA file into the memory, can load multi-FASTA files, can call multiple times
 *
 * @param filename FASTA file with fai index
 */
public void addReference(String filename) {
    File f = new File(filename);
    FastaSequenceFile fa = new FastaSequenceFile(f, true);
    ReferenceSequence seq = fa.nextSequence();
    while (seq != null) {
        ChrString name = new ChrString(seq.getName());
        byte[] seq_bytes = seq.getBases();

        if (seq_bytes == null) {
            log.error("Contig error: " + seq.getName());
        } else {
            log.info("Read ref: " + seq.getName() + ":" + name);
            if (!data.containsKey(name)) {
                Sequence contig = new Sequence(seq.getName(),
                        seq_bytes, seq_bytes.length);
                data.put(name, contig);
            } else {
                log.warn("Duplicate Key!");
            }
        }
        seq = fa.nextSequence();
    }
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
@Override
public void addInfo(final AbstractLocusInfo<EdgingRecordAndOffset> info, final ReferenceSequence ref, boolean referenceBaseN) {
    prepareCollector(info);
    for (final EdgingRecordAndOffset record : info.getRecordAndOffsets()) {
        final String readName = record.getReadName();
        Optional<Set<EdgingRecordAndOffset>> recordsAndOffsetsForName = Optional.ofNullable((Set<EdgingRecordAndOffset>)readsNames.get(readName));
        if (record.getType() == EdgingRecordAndOffset.Type.BEGIN) {
            processRecord(info.getPosition(), ref, record, recordsAndOffsetsForName.orElse(new HashSet<>()));
        } else {
            recordsAndOffsetsForName.ifPresent(
                    edgingRecordAndOffsets -> removeRecordFromMap(record,
                            edgingRecordAndOffsets));
        }
    }
    if (!referenceBaseN) {
        final int readNamesSize = pileupSize.get(info.getPosition());
        final int highQualityDepth = Math.min(readNamesSize, coverageCap);
        if (highQualityDepth < readNamesSize) {
            basesExcludedByCapping += readNamesSize - coverageCap;
        }
        highQualityDepthHistogramArray[highQualityDepth]++;
        unfilteredDepthHistogramArray[unfilteredDepthSize.get(info.getPosition())]++;
    }
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
@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]]++;
        }
    }
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
@Override
List<Bait> design(final BaitDesigner designer, final Interval target, final ReferenceSequence reference) {
    final List<Bait> baits = new LinkedList<Bait>();
    final int baitSize = designer.BAIT_SIZE;
    final int baitOffset = designer.BAIT_OFFSET;
    final int lastPossibleBaitStart = Math.min(target.getEnd(), reference.length() - baitSize);
    final int baitCount = 1 + (int) Math.floor((lastPossibleBaitStart - target.getStart()) / (double) baitOffset);

    int i = 0;
    for (int start = target.getStart(); start < lastPossibleBaitStart; start += baitOffset) {
        final Bait bait = new Bait(target.getContig(),
                start,
                CoordMath.getEnd(start, baitSize),
                target.isNegativeStrand(),
                designer.makeBaitName(target.getName(), ++i, baitCount));
        bait.addBases(reference, designer.DESIGN_ON_TARGET_STRAND);
        baits.add(bait);
    }
    return baits;
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
private static Allele reverseComplement(final Allele oldAllele, final Interval target, final ReferenceSequence referenceSequence, final boolean isBiAllelicIndel, final boolean addToStart){

        if (oldAllele.isSymbolic() || oldAllele.isNoCall()) {
            return oldAllele;
        }
        else if (isBiAllelicIndel) {
            // target.getStart is 1-based, reference bases are 0-based
            final StringBuilder alleleBuilder = new StringBuilder(target.getEnd() - target.getStart() + 1);

            if (addToStart) {
                alleleBuilder.append((char) referenceSequence.getBases()[target.getStart() - 2]);
            }
            alleleBuilder.append(SequenceUtil.reverseComplement(oldAllele.getBaseString().substring(1, oldAllele.length())));
            if (!addToStart) {
                alleleBuilder.append((char) referenceSequence.getBases()[target.getEnd() - 1]);
            }

            return Allele.create(alleleBuilder.toString(), oldAllele.isReference());
        } else {
            return Allele.create(SequenceUtil.reverseComplement(oldAllele.getBaseString()), oldAllele.isReference());
        }
    }
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
@Test
public void testForCollectorWithoutData(){
    long[] templateQualHistogram =  new long[127];
    long[] templateHistogramArray = new long[11];
    CollectWgsMetrics collectWgsMetrics = new CollectWgsMetrics();
    AbstractWgsMetricsCollector collector = new AbstractWgsMetricsCollector(collectWgsMetrics,
            10, createIntervalList()) {
        @Override
        public void addInfo(AbstractLocusInfo info, ReferenceSequence ref, boolean referenceBaseN) {
        }
    };
    assertEquals(templateHistogramArray, collector.highQualityDepthHistogramArray);
    assertEquals(templateQualHistogram, collector.unfilteredBaseQHistogramArray);
    assertEquals(0, collector.basesExcludedByCapping);
    assertEquals(0, collector.basesExcludedByOverlap);
    assertEquals(0, collector.basesExcludedByBaseq);
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
/**
 * Query a specific interval on this reference, and get back all bases spanning that interval at once.
 * Call getBases() on the returned ReferenceSequence to get the actual reference bases. See the BaseUtils
 * class for guidance on how to work with bases in this format.
 *
 * @param contig query interval contig
 * @param start query interval start
 * @param stop query interval stop (included)
 * @return a ReferenceSequence containing all bases spanning the query interval, prefetched
 */
@Override
public ReferenceSequence queryAndPrefetch( final String contig, final long start , final long stop) {
    final int contigIndex = sequenceDictionary.getSequenceIndex(contig);
    int startIndex = (int)(start - bases.getInterval().getStart());
    int length = (int)(stop - start + 1);
    byte[] basesBytes = bases.getBases();
    if (startIndex==0 && length==basesBytes.length) {
        // special case: no need to make a copy
        return new ReferenceSequence(contig, contigIndex, basesBytes);
    }
    Utils.validIndex(startIndex, basesBytes.length);
    Utils.validateArg(startIndex+length <= basesBytes.length, () -> String.format("Asking for stop %d on contig %s but the ReferenceData only has data until %d.", stop, contig, bases.getInterval().getEnd()));
    Utils.validateArg(length >= 0, () -> String.format("Asking for stop<start (%d < %d)", stop, start));
    return new ReferenceSequence(contig, contigIndex, Arrays.copyOfRange(basesBytes, startIndex, startIndex+length));
}
 

import htsjdk.samtools.reference.ReferenceSequence; //导入依赖的package包/类
/**
 * Gets the next complete in-frame codon from the given {@link ReferenceSequence} according to the current codon position and strand.
 * @param referenceSequence The {@link ReferenceSequence} for the current codon.  Must not be {@code null}.
 * @param currentAlignedCodingSequenceAlleleStart The starting position (1-based, inclusive) of the current codon.  Must be > 0.
 * @param currentAlignedCodingSequenceAlleleStop The ending position (1-based, inclusive) of the current codon.  Must be > 0.
 * @param strand The {@link Strand} on which the current codon resides.  Must not be {@code null}.  Must not be {@link Strand#NONE}.
 * @return The next codon in frame with the current codon as specified by the given current codon positions.
 */
private static String getNextReferenceCodon(final ReferenceSequence referenceSequence,
                                            final int currentAlignedCodingSequenceAlleleStart,
                                            final int currentAlignedCodingSequenceAlleleStop,
                                            final Strand strand) {

    Utils.nonNull( referenceSequence );
    ParamUtils.isPositive(currentAlignedCodingSequenceAlleleStart, "Genomic positions must be > 0.");
    ParamUtils.isPositive(currentAlignedCodingSequenceAlleleStop, "Genomic positions must be > 0.");
    assertValidStrand(strand);

    final String nextRefCodon;
    if ( strand == Strand.POSITIVE ) {
        nextRefCodon = referenceSequence.getBaseString().substring(currentAlignedCodingSequenceAlleleStop, currentAlignedCodingSequenceAlleleStop + 3 );
    }
    else {
        nextRefCodon = ReadUtils.getBasesReverseComplement(
                referenceSequence.getBaseString().substring(currentAlignedCodingSequenceAlleleStart - 3, currentAlignedCodingSequenceAlleleStart ).getBytes()
        );
    }
    return nextRefCodon;
}