Java ReferenceSequenceFileFactory.getReferenceSequenceFile方法代码示例

import htsjdk.samtools.reference.ReferenceSequenceFileFactory; //导入方法依赖的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.ReferenceSequenceFileFactory; //导入方法依赖的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.ReferenceSequenceFileFactory; //导入方法依赖的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.ReferenceSequenceFileFactory; //导入方法依赖的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.ReferenceSequenceFileFactory; //导入方法依赖的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.ReferenceSequenceFileFactory; //导入方法依赖的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.ReferenceSequenceFileFactory; //导入方法依赖的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.ReferenceSequenceFileFactory; //导入方法依赖的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.ReferenceSequenceFileFactory; //导入方法依赖的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.ReferenceSequenceFileFactory; //导入方法依赖的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;
    }
}
 

import htsjdk.samtools.reference.ReferenceSequenceFileFactory; //导入方法依赖的package包/类
public ReferenceSource(File file) {
	if (file != null) {
		rsFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(file);

		File indexFile = new File(file.getAbsoluteFile() + ".fai");
		if (indexFile.exists())
			fastaSequenceIndex = new FastaSequenceIndex(indexFile);
	}
}
 

import htsjdk.samtools.reference.ReferenceSequenceFileFactory; //导入方法依赖的package包/类
public static byte[] getBasesFromReferenceFile(String referenceFilePath, String seqName, int from, int length) {
	ReferenceSequenceFile referenceSequenceFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(new File(
			referenceFilePath));
	ReferenceSequence sequence = referenceSequenceFile.getSequence(seqName);
	byte[] bases = referenceSequenceFile.getSubsequenceAt(sequence.getName(), from, from + length).getBases();
	return bases;
}
 

import htsjdk.samtools.reference.ReferenceSequenceFileFactory; //导入方法依赖的package包/类
public static void main(String[] args) {
	ReferenceSequenceFile rsf = ReferenceSequenceFileFactory
			.getReferenceSequenceFile(new File("data/hs37d5.fa"));
	SAMSequenceDictionary dict = rsf.getSequenceDictionary();
	ReferenceSequence rs1 = rsf.nextSequence(), rs2 = rsf.nextSequence();
	ReferenceSequence rs = rsf.getSequence("1"); // get reference sequence
													// by sequence name
	System.out.println(rs.getName());
	System.out.println(rs.getContigIndex());
	System.out.println(rs.length());
	System.out.println(rs.getBases().length);
	System.out.println("===============");
	byte[] bases = rs.getBases();
	int count[] = new int[256];
	for (int i = 0; i < bases.length; i++) {
		count[bases[i]]++;
	}
	for (int i = 0; i < count.length; i++) {
		System.out.println((char) i + " " + count[i]);
	}

	rs = rsf.getSequence("22");
	System.out.println(rs.getName());
	System.out.println(rs.getContigIndex());
	System.out.println(rs.length());
	System.out.println(rs.getBases().length);
	System.out.println("===============");

	rs = rsf.getSequence("X");
	System.out.println(rs.getName());
	System.out.println(rs.getContigIndex());
	System.out.println(rs.length());
	System.out.println(rs.getBases().length);
	System.out.println(rs.getBases()[rs.length() - 1]);

	rs = rsf.getSubsequenceAt("X", 0, 2);
	System.out.println(rs.getName());
	System.out.println(rs.getContigIndex());
	System.out.println(rs.length());
	System.out.println(rs.getBases().length);
	System.out.println(String.format("%02x", rs.getBases()[0]));

	// for (int i=0;i<bases.length;++i) {
	// int unsigned = bases[i]&0xff;
	//
	// System.out.println("\""+(char)(unsigned)+"\" ("+unsigned+")");
	// }
}
 

import htsjdk.samtools.reference.ReferenceSequenceFileFactory; //导入方法依赖的package包/类
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsReadable(REFERENCE_SEQUENCE);
    IOUtil.assertFileIsWritable(OUTPUT);

    final SamReader in = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).open(INPUT);

    final ReferenceSequenceFile reference = ReferenceSequenceFileFactory.getReferenceSequenceFile(REFERENCE_SEQUENCE);
    final SAMSequenceDictionary refDict = reference.getSequenceDictionary();

    if (refDict == null) {
        log.error("No reference sequence dictionary found. Aborting.  You can create a sequence dictionary for the reference fasta using CreateSequenceDictionary.jar.");
        CloserUtil.close(in);
        return 1;
    }

    printDictionary("SAM/BAM file", in.getFileHeader().getSequenceDictionary());
    printDictionary("Reference", refDict);
    final Map<Integer, Integer> newOrder = buildSequenceDictionaryMap(refDict, in.getFileHeader().getSequenceDictionary());

    // has to be after we create the newOrder
    final SAMFileHeader outHeader = in.getFileHeader().clone();
    outHeader.setSequenceDictionary(refDict);

    log.info("Writing reads...");
    if (in.hasIndex()) {
        try( final SAMFileWriter out = new SAMFileWriterFactory().makeSAMOrBAMWriter(outHeader, true, OUTPUT)) {

            // write the reads in contig order
            for (final SAMSequenceRecord contig : refDict.getSequences()) {
                final SAMRecordIterator it = in.query(contig.getSequenceName(), 0, 0, false);
                writeReads(out, it, newOrder, contig.getSequenceName());
            }
            // don't forget the unmapped reads
            writeReads(out, in.queryUnmapped(), newOrder, "unmapped");
        }
    } else {
        try (final SAMFileWriter out = new SAMFileWriterFactory().makeSAMOrBAMWriter(outHeader, false, OUTPUT)) {
            writeReads(out, in.iterator(), newOrder, "All reads");
        }
    }

    // cleanup
    CloserUtil.close(in);
    return 0;
}
 

import htsjdk.samtools.reference.ReferenceSequenceFileFactory; //导入方法依赖的package包/类
/**
 * Asserts that files are readable and writable and then fires off an
 * HsMetricsCalculator instance to do the real work.
 */
protected int doWork() {
    for (final File targetInterval : TARGET_INTERVALS) IOUtil.assertFileIsReadable(targetInterval);
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsWritable(OUTPUT);
    if (PER_TARGET_COVERAGE != null) IOUtil.assertFileIsWritable(PER_TARGET_COVERAGE);

    final SamReader reader = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).open(INPUT);
    final IntervalList targetIntervals = IntervalList.fromFiles(TARGET_INTERVALS);

    // Validate that the targets and baits have the same references as the reads file
    SequenceUtil.assertSequenceDictionariesEqual(
            reader.getFileHeader().getSequenceDictionary(),
            targetIntervals.getHeader().getSequenceDictionary());
    SequenceUtil.assertSequenceDictionariesEqual(
            reader.getFileHeader().getSequenceDictionary(),
            getProbeIntervals().getHeader().getSequenceDictionary()
    );

    ReferenceSequenceFile ref = null;
    if (REFERENCE_SEQUENCE != null) {
        IOUtil.assertFileIsReadable(REFERENCE_SEQUENCE);
        ref = ReferenceSequenceFileFactory.getReferenceSequenceFile(REFERENCE_SEQUENCE);
        SequenceUtil.assertSequenceDictionariesEqual(
                reader.getFileHeader().getSequenceDictionary(), ref.getSequenceDictionary(),
                INPUT, REFERENCE_SEQUENCE
        );
    }

    final COLLECTOR collector = makeCollector(
            METRIC_ACCUMULATION_LEVEL,
            reader.getFileHeader().getReadGroups(),
            ref,
            PER_TARGET_COVERAGE,
            PER_BASE_COVERAGE,
            targetIntervals,
            getProbeIntervals(),
            getProbeSetName(),
            NEAR_DISTANCE
    );

    final ProgressLogger progress = new ProgressLogger(log);
    for (final SAMRecord record : reader) {
        collector.acceptRecord(record, null);
        progress.record(record);
    }

    // Write the output file
    final MetricsFile<METRIC, Integer> metrics = getMetricsFile();
    collector.finish();

    collector.addAllLevelsToFile(metrics);

    metrics.write(OUTPUT);

    CloserUtil.close(reader);
    return 0;
}