import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
@Test @Transactional(propagation = Propagation.REQUIRES_NEW, rollbackFor = Exception.class)
public void nibFileAndFastaFileIdenticalTest() throws ReferenceReadingException, IOException {
    ReferenceSequenceFile referenceSequenceFile =
            ReferenceSequenceFileFactory.getReferenceSequenceFile(resource.getFile());
    //Go through all chromosomes if FASTA-file
    Track<Sequence> track = new Track<>();
    track.setType(TrackType.REF);
    track.setScaleFactor(SCALE_FACTOR_4_BASE);
    track.setId(reference.getId());

    byte[] arrayOfNucleicAcids;

    //Nucleicide frome FASTA-file
    arrayOfNucleicAcids = referenceSequenceFile.nextSequence().getBases();
    Chromosome chromosome = reference.getChromosomes().get(0);
    track.setChromosome(chromosome);
    track.setStartIndex(START_INDEX);
    track.setEndIndex(END_INDEX);
    //Nucleicide frome nib-file
    track = referenceManager.getNucleotidesResultFromNib(track);
    for (Sequence sequence : track.getBlocks()) {
        assertTrue(String.valueOf((char) (arrayOfNucleicAcids[sequence.getStartIndex() - 1]))
                .equals(sequence.getText()));
    }
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
public static ReferenceSequenceFile findReferenceFile(String fileName)
		throws FileNotFoundException {
	ReferenceSequenceFile rsf;

	File refFile = new File(fileName);

	try {
		rsf = ReferenceSequenceFileFactory
				.getReferenceSequenceFile(refFile);
	} catch (SAMException ex) {
		System.err.println("Could not load reference sequence file \""
				+ refFile + "\".");
		throw new FileNotFoundException(refFile.toString());
	}
	return rsf;
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的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.ReferenceSequenceFile; //导入依赖的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.ReferenceSequenceFile; //导入依赖的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.ReferenceSequenceFile; //导入依赖的package包/类
public HsMetricCollector(final Set<MetricAccumulationLevel> accumulationLevels,
                         final List<SAMReadGroupRecord> samRgRecords,
                         final ReferenceSequenceFile refFile,
                         final File perTargetCoverage,
                         final File perBaseCoverage,
                         final IntervalList targetIntervals,
                         final IntervalList probeIntervals,
                         final String probeSetName,
                         final int nearProbeDistance,
                         final int minimumMappingQuality,
                         final int minimumBaseQuality,
                         final boolean clipOverlappingReads,
                         final int coverageCap,
                         final int sampleSize) {
    super(accumulationLevels, samRgRecords, refFile, perTargetCoverage, perBaseCoverage, targetIntervals, probeIntervals, probeSetName, nearProbeDistance, minimumMappingQuality, minimumBaseQuality, clipOverlappingReads, coverageCap, sampleSize);
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
public TargetMetricsCollector(final Set<MetricAccumulationLevel> accumulationLevels,
                              final List<SAMReadGroupRecord> samRgRecords,
                              final ReferenceSequenceFile refFile,
                              final File perTargetCoverage,
                              final File perBaseCoverage,
                              final IntervalList targetIntervals,
                              final IntervalList probeIntervals,
                              final String probeSetName,
                              final int nearProbeDistance,
                              final int minimumMappingQuality,
                              final int minimumBaseQuality,
                              final boolean clipOverlappingReads,
                              final int coverageCap,
                              final int sampleSize) {
    this(accumulationLevels, samRgRecords, refFile, perTargetCoverage, perBaseCoverage, targetIntervals, probeIntervals, probeSetName, nearProbeDistance, minimumMappingQuality, minimumBaseQuality, clipOverlappingReads, false, coverageCap, sampleSize);
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
public TargetedPcrMetricsCollector(final Set<MetricAccumulationLevel> accumulationLevels,
                                   final List<SAMReadGroupRecord> samRgRecords,
                                   final ReferenceSequenceFile refFile,
                                   final File perTargetCoverage,
                                   final File perBaseCoverage,
                                   final IntervalList targetIntervals,
                                   final IntervalList probeIntervals,
                                   final String probeSetName,
                                   final int nearProbeDistance,
                                   final int minimumMappingQuality,
                                   final int minimumBaseQuality,
                                   final boolean clipOverlappingReads,
                                   final int coverageCap,
                                   final int sampleSize) {
    super(accumulationLevels, samRgRecords, refFile, perTargetCoverage, perBaseCoverage, targetIntervals, probeIntervals, probeSetName, nearProbeDistance, minimumMappingQuality, minimumBaseQuality, clipOverlappingReads, coverageCap, sampleSize);
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
@Override
public List<Entry> call() throws Exception {
	List<Entry> list = new ArrayList<RefRepo.Entry>();

	ReferenceSequenceFile rsFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(file);

	ReferenceSequence sequence = null;
	while ((sequence = rsFile.nextSequence()) != null) {
		sequence.getBases();

		Entry e = new Entry();
		e.md5 = Utils.calculateMD5String(sequence.getBases());
		e.file = "file://" + file.getAbsolutePath();
		e.name = sequence.getName();
		e.length = sequence.length();
		log.info(String.format("New entry: %s", e.toString()));
		list.add(e);
	}
	return list;
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
/**
 * @param header SECRAM file header
 * @param inputStream The input stream where we read SECRAM records.
 * @param referenceFile Reference sequence file.
 * @param filter Security filter for decryption.
 */
public SECRAMIterator(SecramHeader header, InputStream inputStream,
		ReferenceSequenceFile referenceFile, SECRAMSecurityFilter filter) {
	this.secramHeader = header;
	this.mRsf = referenceFile;
	this.filter = filter;
	this.containerIterator = new SECRAMContainerIterator(inputStream,
			filter);
	this.parser = new SecramContainerParser();
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
@Override
public void onTraversalStart() {
    final ReferenceSequenceFile referenceReader = getReferenceReader(referenceArguments);
    hcEngine = new HaplotypeCallerEngine(hcArgs, getHeaderForReads(), referenceReader);

    // The HC engine will make the right kind (VCF or GVCF) of writer for us
    final SAMSequenceDictionary sequenceDictionary = getHeaderForReads().getSequenceDictionary();
    vcfWriter = hcEngine.makeVCFWriter(outputVCF, sequenceDictionary);
    hcEngine.writeHeader(vcfWriter, sequenceDictionary, getDefaultToolVCFHeaderLines());
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
public void writeAsVcf(final File output, final File refFile) throws FileNotFoundException {
    ReferenceSequenceFile ref = new IndexedFastaSequenceFile(refFile);
    try (VariantContextWriter writer = new VariantContextWriterBuilder()
            .setOutputFile(output)
            .setReferenceDictionary(ref.getSequenceDictionary())
            .build()) {

        final VCFHeader vcfHeader = new VCFHeader(
                VCFUtils.withUpdatedContigsAsLines(Collections.emptySet(), refFile, header.getSequenceDictionary(), false),
                Collections.singleton(HET_GENOTYPE_FOR_PHASING));

        VCFUtils.withUpdatedContigsAsLines(Collections.emptySet(), refFile, header.getSequenceDictionary(), false);

        vcfHeader.addMetaDataLine(new VCFHeaderLine(VCFHeaderVersion.VCF4_2.getFormatString(), VCFHeaderVersion.VCF4_2.getVersionString()));
        vcfHeader.addMetaDataLine(new VCFInfoHeaderLine(VCFConstants.ALLELE_FREQUENCY_KEY, VCFHeaderLineCount.A, VCFHeaderLineType.Float, "Allele Frequency, for each ALT allele, in the same order as listed"));
        vcfHeader.addMetaDataLine(new VCFFormatHeaderLine(VCFConstants.GENOTYPE_KEY, 1, VCFHeaderLineType.String, "Genotype"));
        vcfHeader.addMetaDataLine(new VCFFormatHeaderLine(VCFConstants.PHASE_SET_KEY, 1, VCFHeaderLineType.String, "Phase-set identifier for phased genotypes."));
        vcfHeader.addMetaDataLine(new VCFHeaderLine(VCFHeader.SOURCE_KEY,"HaplotypeMap::writeAsVcf"));
        vcfHeader.addMetaDataLine(new VCFHeaderLine("reference","HaplotypeMap::writeAsVcf"));


      //  vcfHeader.addMetaDataLine(new VCFHeaderLine());
        writer.writeHeader(vcfHeader);
        final LinkedList<VariantContext> variants = new LinkedList<>(this.asVcf(ref));
        variants.sort(vcfHeader.getVCFRecordComparator());
        variants.forEach(writer::add);
    }
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
/**
 * A factory method for the TargetMetricsCollector to use this time.  Examples of TargetMetricsCollector:
 * (TargetedPcrMetricsCollector, HsMetricsCalculator)
 *
 * @return A TargetMetricsCollector to which we will pass SAMRecords
 */
protected abstract COLLECTOR makeCollector(final Set<MetricAccumulationLevel> accumulationLevels,
                                           final List<SAMReadGroupRecord> samRgRecords,
                                           final ReferenceSequenceFile refFile,
                                           final File perTargetCoverage,
                                           final File perBaseCoverage,
                                           final IntervalList targetIntervals,
                                           final IntervalList probeIntervals,
                                           final String probeSetName,
                                           final int nearProbeDistance);
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
@Override
protected TargetedPcrMetricsCollector makeCollector(final Set<MetricAccumulationLevel> accumulationLevels,
                                                    final List<SAMReadGroupRecord> samRgRecords,
                                                    final ReferenceSequenceFile refFile,
                                                    final File perTargetCoverage,
                                                    final File perBaseCoverage,
                                                    final IntervalList targetIntervals,
                                                    final IntervalList probeIntervals,
                                                    final String probeSetName,
                                                    final int nearProbeDistance) {
    return new TargetedPcrMetricsCollector(accumulationLevels, samRgRecords, refFile, perTargetCoverage, perBaseCoverage, targetIntervals, probeIntervals, probeSetName, nearProbeDistance,
            MINIMUM_MAPPING_QUALITY, MINIMUM_BASE_QUALITY, CLIP_OVERLAPPING_READS, true, COVERAGE_CAP, SAMPLE_SIZE);
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
@Override
protected HsMetricCollector makeCollector(final Set<MetricAccumulationLevel> accumulationLevels,
                                          final List<SAMReadGroupRecord> samRgRecords,
                                          final ReferenceSequenceFile refFile,
                                          final File perTargetCoverage,
                                          final File perBaseCoverage,
                                          final IntervalList targetIntervals,
                                          final IntervalList probeIntervals,
                                          final String probeSetName,
                                          final int nearProbeDistance) {
    return new HsMetricCollector(accumulationLevels, samRgRecords, refFile, perTargetCoverage, perBaseCoverage, targetIntervals, probeIntervals, probeSetName, nearProbeDistance,
            MINIMUM_MAPPING_QUALITY, MINIMUM_BASE_QUALITY, CLIP_OVERLAPPING_READS, true, COVERAGE_CAP, SAMPLE_SIZE);
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(REFERENCE_SEQUENCE);
    IOUtil.assertFileIsWritable(OUTPUT);

    final ReferenceSequenceFile refFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(REFERENCE_SEQUENCE, true);
    if (!refFile.isIndexed()) {
        throw new IllegalStateException("Reference file must be indexed, but no index file was found");
    }
    if (refFile.getSequenceDictionary() == null) {
        throw new IllegalStateException("Reference file must include a dictionary, but no dictionary file was found");
    }

    // get the intervals
    final IntervalList intervals = segregateReference(refFile, MAX_TO_MERGE);

    log.info(String.format("Found %d intervals in %d loci during %s seconds", intervalProgress.getCount(), locusProgress.getCount(), locusProgress.getElapsedSeconds()));

    /**********************************
     * Now output regions for calling *
     **********************************/

    final IntervalList outputIntervals = new IntervalList(intervals.getHeader().clone());
    log.info(String.format("Collecting requested type of intervals (%s)", OUTPUT_TYPE));

    intervals.getIntervals().stream().filter(i -> OUTPUT_TYPE.accepts(i.getName())).forEach(outputIntervals::add);

    log.info("Writing Intervals.");
    outputIntervals.write(OUTPUT);

    log.info(String.format("Execution ending. Total time %d seconds", locusProgress.getElapsedSeconds()));

    return 0;
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INTERVAL_LIST);
    IOUtil.assertFileIsReadable(REFERENCE_SEQUENCE);
    IOUtil.assertFileIsWritable(OUTPUT);

    final IntervalList intervals = IntervalList.fromFile(INTERVAL_LIST);
    final ReferenceSequenceFile ref = ReferenceSequenceFileFactory.getReferenceSequenceFile(REFERENCE_SEQUENCE);
    SequenceUtil.assertSequenceDictionariesEqual(intervals.getHeader().getSequenceDictionary(), ref.getSequenceDictionary());

    final BufferedWriter out = IOUtil.openFileForBufferedWriting(OUTPUT);

    for (final Interval interval : intervals) {
        final ReferenceSequence seq = ref.getSubsequenceAt(interval.getContig(), interval.getStart(), interval.getEnd());
        final byte[] bases = seq.getBases();
        if (interval.isNegativeStrand()) SequenceUtil.reverseComplement(bases);

        try {
            out.write(">");
            out.write(interval.getName());
            out.write("\n");

            for (int i=0; i<bases.length; ++i) {
                if (i > 0 && i % LINE_LENGTH == 0) out.write("\n");
                out.write(bases[i]);
            }

            out.write("\n");
        }
        catch (IOException ioe) {
            throw new PicardException("Error writing to file " + OUTPUT.getAbsolutePath(), ioe);

        }
    }

    CloserUtil.close(out);

    return 0;
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
@Test
public void testForFilteredBases(){
    AbstractLocusIterator iterator = createReadEndsIterator(exampleSamOneRead);
    CollectWgsMetrics collectWgsMetrics = new CollectWgsMetrics();
    FastWgsMetricsCollector collector =  new FastWgsMetricsCollector(collectWgsMetrics, 100, createIntervalList());
    String secondReferenceString = ">ref\nNNNNNNNNNNAATATTCTTC";
    ReferenceSequenceFile referenceSequenceFile = createReferenceSequenceFile(secondReferenceString);
    ReferenceSequenceFileWalker refWalker =  new ReferenceSequenceFileWalker(referenceSequenceFile);
    WgsMetricsProcessorImpl wgsMetricsProcessor =
            new WgsMetricsProcessorImpl(iterator, refWalker, collector, progress);
    wgsMetricsProcessor.processFile();
    assertEquals(10, collector.counter);
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
@Override
public ReferenceBases getReferenceBases(final SimpleInterval interval) throws IOException {
    try ( ReferenceSequenceFile referenceSequenceFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(getReferencePath()) ) {
        ReferenceSequence sequence = referenceSequenceFile.getSubsequenceAt(interval.getContig(), interval.getStart(), interval.getEnd());
        return new ReferenceBases(sequence.getBases(), interval);
    }
}
 

import htsjdk.samtools.reference.ReferenceSequenceFile; //导入依赖的package包/类
public Map<String, ReferenceBases> getAllReferenceBases() throws IOException {
    try ( ReferenceSequenceFile referenceSequenceFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(getReferencePath()) ) {
        Map<String, ReferenceBases> bases = new LinkedHashMap<>();
        ReferenceSequence seq;
        while ( (seq = referenceSequenceFile.nextSequence()) != null ) {
            String name = seq.getName();
            bases.put(name, new ReferenceBases(seq.getBases(), new SimpleInterval(name, 1, seq.length())));
        }
        return bases;
    }
}