Java Code Examples for htsjdk.samtools.reference.ReferenceSequenceFileFactory#getReferenceSequenceFile()

@Override
protected int doWork() {
	IOUtil.assertFileIsReadable(this.REFERENCE_SEQUENCE);
	IOUtil.assertFileIsWritable(this.OUTPUT);
	// validate that an index is present for the reference sequence, since it's required.
	final ReferenceSequenceFile ref = ReferenceSequenceFileFactory.getReferenceSequenceFile(REFERENCE_SEQUENCE, true, true);
	if (!ref.isIndexed())
		throw new IllegalStateException ("Input fasta must be indexed.  You can do this by using samtools faidx to create an index");

	FastaSequenceFileWriter writer = new FastaSequenceFileWriter(OUTPUT, OUTPUT_LINE_LENGTH);
	if (this.CONTIG_PATTERN_TO_IGNORE!=null && !this.CONTIG_PATTERN_TO_IGNORE.isEmpty()) processByWholeContig(ref, writer, this.CONTIG_PATTERN_TO_IGNORE);
	if (this.INTERVALS!=null) processByPartialContig(ref, writer, this.INTERVALS);

	CloserUtil.close(ref);
	CloserUtil.close(writer);
	return 0;

}
 

@Override
protected String[] customCommandLineValidation() {

    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsReadable(EXTENDED_ILLUMINA_MANIFEST);
    IOUtil.assertFileIsReadable(ILLUMINA_BEAD_POOL_MANIFEST_FILE);
    IOUtil.assertFileIsWritable(OUTPUT);
    refSeq = ReferenceSequenceFileFactory.getReferenceSequenceFile(REFERENCE_SEQUENCE);
    final SAMSequenceDictionary sequenceDictionary = refSeq.getSequenceDictionary();
    final String assembly = sequenceDictionary.getSequence(0).getAssembly();
    if (!assembly.equals("GRCh37")) {
        return new String[]{"The selected reference sequence ('" + assembly + "') is not supported.  This tool is currently only implemented to support NCBI Build 37 / HG19 Reference Sequence."};
    }

    if (FINGERPRINT_GENOTYPES_VCF_FILE != null) {
        IOUtil.assertFileIsReadable(FINGERPRINT_GENOTYPES_VCF_FILE);
    }
    if (GENDER_GTC != null) {
        IOUtil.assertFileIsReadable(GENDER_GTC);
    }

    return super.customCommandLineValidation();
}
 

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;
}
 

@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;
}
 

private void validateReferenceBases(File referenceFile) {
    final ReferenceSequenceFile refSeqFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(referenceFile, true);
    ReferenceSequence sequence;
    while ((sequence = refSeqFile.nextSequence()) != null) {
        for (final byte base: sequence.getBases()) {
            if (!IUPAC_TABLE[base]) {
                messages.baseErrors = String.format("WARNING: AT least one invalid base '%c' (decimal %d) in reference sequence named %s",
                        StringUtil.byteToChar(base), base, sequence.getName());
                break;
            }
        }
    }
}
 

/**
 * A function that takes a Fingerprint and writes it as a VCF to a file
 *
 * @param fingerprint               the fingerprint to write
 * @param outputFile                the file to write to
 * @param referenceSequenceFileName the reference sequence (file)
 * @param sample                    the sample name to use in the vcf
 * @param source                    a "source" comment to use in the VCF
 * @throws IOException
 */
public static void writeFingerPrint(final Fingerprint fingerprint,
                                    final File outputFile,
                                    final File referenceSequenceFileName,
                                    final String sample,
                                    final String source) throws IOException {

    try (final ReferenceSequenceFile ref = ReferenceSequenceFileFactory.getReferenceSequenceFile(referenceSequenceFileName);
         final VariantContextWriter variantContextWriter = getVariantContextWriter(outputFile, referenceSequenceFileName, sample, source, ref)) {

        createVCSetFromFingerprint(fingerprint, ref, sample).forEach(variantContextWriter::add);
    }
}
 

private SAMSequenceDictionary makeSequenceDictionary(final File referenceFile) {
    final ReferenceSequenceFile refSeqFile =
            ReferenceSequenceFileFactory.getReferenceSequenceFile(referenceFile, true);
    ReferenceSequence refSeq;
    final List<SAMSequenceRecord> ret = new ArrayList<>();
    final Set<String> sequenceNames = new HashSet<>();
    while ((refSeq = refSeqFile.nextSequence()) != null) {
        if (sequenceNames.contains(refSeq.getName())) {
            throw new PicardException("Sequence name appears more than once in reference: " + refSeq.getName());
        }
        sequenceNames.add(refSeq.getName());
        ret.add(new SAMSequenceRecord(refSeq.getName(), refSeq.length()));
    }
    return new SAMSequenceDictionary(ret);
}
 

public SamSequenceRecordsIterator(File referenceSequence, boolean truncateNamesAtWhitespace) {
    this.truncateNamesAtWhitespace = truncateNamesAtWhitespace;
    this.refSeqFile = ReferenceSequenceFileFactory.
            getReferenceSequenceFile(referenceSequence, truncateNamesAtWhitespace);

    this.nextRefSeq = refSeqFile.nextSequence();
    try {
        md5 = MessageDigest.getInstance("MD5");
    } catch (NoSuchAlgorithmException e) {
        throw new PicardException("MD5 algorithm not found", e);
    }
}
 

@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;
}
 

@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;
}
 

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);
	}
}
 

@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);
    }
}
 

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;
    }
}
 

public static void main(String[] args) throws IOException, IllegalArgumentException, IllegalAccessException {
	Params params = new Params();
	JCommander jc = new JCommander(params);
	try {
		jc.parse(args);
	} catch (Exception e) {
		System.out.println("Failed to parse parameteres, detailed message below: ");
		System.out.println(e.getMessage());
		System.out.println();
		System.out.println("See usage: -h");
		System.exit(1);
	}

	if (args.length == 0 || params.help) {
		printUsage(jc);
		System.exit(1);
	}

	if (params.reference == null) {
		System.out.println("A reference fasta file is required.");
		System.exit(1);
	}

	if (params.cramFile == null) {
		System.out.println("A CRAM input file is required. ");
		System.exit(1);
	}

	Log.setGlobalLogLevel(Log.LogLevel.INFO);

	ReferenceSequenceFile referenceSequenceFile = ReferenceSequenceFileFactory
			.getReferenceSequenceFile(params.reference);

	FileInputStream fis = new FileInputStream(params.cramFile);
	BufferedInputStream bis = new BufferedInputStream(fis);

	CRAMIterator iterator = new CRAMIterator(bis, new ReferenceSource(params.reference),
			ValidationStringency.STRICT);
	CramHeader cramHeader = iterator.getCramHeader();

	iterator.close();

	ProgressLogger progress = new ProgressLogger(log, 100000, "Validated Read");
	SamFileValidator v = new SamFileValidator(new PrintWriter(System.out), 1);
	final SamReader reader = SamReaderFactory.make().referenceSequence(params.reference).open(params.cramFile);
	List<SAMValidationError.Type> errors = new ArrayList<SAMValidationError.Type>();
	errors.add(SAMValidationError.Type.MATE_NOT_FOUND);
	// errors.add(Type.MISSING_TAG_NM);
	v.setErrorsToIgnore(errors);
	v.validateSamFileSummary(reader, ReferenceSequenceFileFactory.getReferenceSequenceFile(params.reference));
	log.info("Elapsed seconds: " + progress.getElapsedSeconds());
}
 

@Override
public SAMSequenceDictionary getReferenceSequenceDictionary(final SAMSequenceDictionary optReadSequenceDictionaryToMatch) {
    ReferenceSequenceFile referenceSequenceFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(new GATKPath(referenceURI.toString()).toPath());
    return referenceSequenceFile.getSequenceDictionary();
}
 

@Override
public ReferenceBases getReferenceBases(final SimpleInterval interval) {
    ReferenceSequenceFile referenceSequenceFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(new GATKPath(referenceURI.toString()).toPath());
    ReferenceSequence sequence = referenceSequenceFile.getSubsequenceAt(interval.getContig(), interval.getStart(), interval.getEnd());
    return new ReferenceBases(sequence.getBases(), interval);
}
 

@Override
public SAMSequenceDictionary getReferenceSequenceDictionary(final SAMSequenceDictionary optReadSequenceDictionaryToMatch) throws IOException {
    try ( ReferenceSequenceFile referenceSequenceFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(getReferencePath()) ) {
        return referenceSequenceFile.getSequenceDictionary();
    }
}
 

@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsWritable(OUTPUT);

    // set up the reference and a mask so that we only count the positions requested by the user
    final ReferenceSequenceFile ref = ReferenceSequenceFileFactory.getReferenceSequenceFile(INPUT);
    final ReferenceSequenceMask referenceSequenceMask;
    if (INTERVALS != null) {
        IOUtil.assertFileIsReadable(INTERVALS);
        final IntervalList intervalList = IntervalList.fromFile(INTERVALS);
        referenceSequenceMask = new IntervalListReferenceSequenceMask(intervalList);
    } else {
        final SAMFileHeader header = new SAMFileHeader();
        header.setSequenceDictionary(ref.getSequenceDictionary());
        referenceSequenceMask = new WholeGenomeReferenceSequenceMask(header);
    }

    long nonNbases = 0L;

    for (final SAMSequenceRecord rec : ref.getSequenceDictionary().getSequences()) {
        // pull out the contig and set up the bases
        final ReferenceSequence sequence = ref.getSequence(rec.getSequenceName());
        final byte[] bases = sequence.getBases();
        StringUtil.toUpperCase(bases);

        for (int i = 0; i < bases.length; i++) {
            // only investigate this position if it's within our mask
            if (referenceSequenceMask.get(sequence.getContigIndex(), i+1)) {
                nonNbases += bases[i] == SequenceUtil.N ? 0 : 1;
            }
        }
    }

    try {
        final BufferedWriter out = IOUtil.openFileForBufferedWriting(OUTPUT);
        out.write(nonNbases + "\n");
        out.close();
    }
    catch (IOException ioe) {
        throw new PicardException("Error writing to file " + OUTPUT.getAbsolutePath(), ioe);
    }

    return 0;
}
 

@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsWritable(OUTPUT);

    if (INPUT.getAbsoluteFile().equals(OUTPUT.getAbsoluteFile())) {
        throw new IllegalArgumentException("Input and output cannot be the same file.");
    }

    final ReferenceSequenceFile ref = ReferenceSequenceFileFactory.getReferenceSequenceFile(INPUT, TRUNCATE_SEQUENCE_NAMES_AT_WHITESPACE);
    final BufferedWriter out = IOUtil.openFileForBufferedWriting(OUTPUT);

    ReferenceSequence seq = null;
    while ((seq = ref.nextSequence()) != null) {
        final String name  = seq.getName();
        final byte[] bases = seq.getBases();

        try {
            out.write(">");
            out.write(name);
            out.newLine();

            if (bases.length == 0) {
                log.warn("Sequence " + name + " contains 0 bases.");
            }
            else {
                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);

        }
    }
    try {
        out.close();
    } catch (IOException e) {
        throw new RuntimeIOException(e);
    }
    return 0;
}
 

/**
 * 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);

    if (THEORETICAL_SENSITIVITY_OUTPUT != null) {
        // Write out theoretical sensitivity results.
        final MetricsFile<TheoreticalSensitivityMetrics, ?> theoreticalSensitivityMetrics = getMetricsFile();
        log.info("Calculating theoretical sentitivity at " + ALLELE_FRACTION.size() + " allele fractions.");
        List<TheoreticalSensitivityMetrics> tsm = TheoreticalSensitivity.calculateSensitivities(SAMPLE_SIZE, collector.getDepthHistogram(), collector.getBaseQualityHistogram(), ALLELE_FRACTION);
        theoreticalSensitivityMetrics.addAllMetrics(tsm);
        theoreticalSensitivityMetrics.write(THEORETICAL_SENSITIVITY_OUTPUT);
    }

    CloserUtil.close(reader);
    return 0;
}