Examples with SamReaderFactory htsjdk.samtools.SamReaderFactory used on opensource projects

Example 16 with SamReaderFactory

use of htsjdk.samtools.SamReaderFactory in project gridss by PapenfussLab.

the class ComputeSamTags method doWork.

@Override
protected int doWork() {
    log.debug("Setting language-neutral locale");
    java.util.Locale.setDefault(Locale.ROOT);
    validateParameters();
    SamReaderFactory readerFactory = SamReaderFactory.make();
    SAMFileWriterFactory writerFactory = new SAMFileWriterFactory();
    try {
        try (SamReader reader = readerFactory.open(INPUT)) {
            SAMFileHeader header = reader.getFileHeader();
            if (!ASSUME_SORTED) {
                if (header.getSortOrder() != SortOrder.queryname) {
                    log.error("INPUT is not sorted by queryname. " + "ComputeSamTags requires that reads with the same name be sorted together. " + "If the input file satisfies this constraint (the output from many aligners do)," + " this check can be disabled with the ASSUME_SORTED option.");
                    return -1;
                }
            }
            try (SAMRecordIterator it = reader.iterator()) {
                File tmpoutput = gridss.Defaults.OUTPUT_TO_TEMP_FILE ? FileSystemContext.getWorkingFileFor(OUTPUT, "gridss.tmp.ComputeSamTags.") : OUTPUT;
                try (SAMFileWriter writer = writerFactory.makeSAMOrBAMWriter(header, true, tmpoutput)) {
                    compute(it, writer, getReference(), TAGS, SOFTEN_HARD_CLIPS, FIX_MATE_INFORMATION, RECALCULATE_SA_SUPPLEMENTARY, INPUT.getName() + "-");
                }
                if (tmpoutput != OUTPUT) {
                    FileHelper.move(tmpoutput, OUTPUT, true);
                }
            }
        }
    } catch (IOException e) {
        log.error(e);
        return -1;
    }
    return 0;
}

Example 17 with SamReaderFactory

use of htsjdk.samtools.SamReaderFactory in project gridss by PapenfussLab.

the class MultipleSamFileCommandLineProgram method ensureDictionariesMatch.

public void ensureDictionariesMatch() throws IOException {
    ReferenceSequenceFile ref = null;
    try {
        ref = ReferenceSequenceFileFactory.getReferenceSequenceFile(REFERENCE_SEQUENCE);
        SAMSequenceDictionary dictionary = ref.getSequenceDictionary();
        final SamReaderFactory samFactory = SamReaderFactory.makeDefault();
        for (File f : INPUT) {
            SamReader reader = null;
            try {
                reader = samFactory.open(f);
                SequenceUtil.assertSequenceDictionariesEqual(reader.getFileHeader().getSequenceDictionary(), dictionary, f, REFERENCE_SEQUENCE);
            } catch (htsjdk.samtools.util.SequenceUtil.SequenceListsDifferException e) {
                log.error("Reference genome used by ", f, " does not match reference genome ", REFERENCE_SEQUENCE, ". ", "The reference supplied must match the reference used for every input.");
                throw e;
            } finally {
                if (reader != null)
                    reader.close();
            }
        }
    } finally {
        if (ref != null)
            ref.close();
    }
}

Example 18 with SamReaderFactory

use of htsjdk.samtools.SamReaderFactory in project gridss by PapenfussLab.

the class BamToBed method extract.

public static void extract(File input, File splitReads, File spanningReads, final double minLengthPortion, final int minMapq, final int minIndelSize, final int maxMappedBases, final int minIndelComponentSize) throws IOException {
    SamReaderFactory factory = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.LENIENT);
    SamReader reader = factory.open(input);
    BufferedWriter srwriter = null;
    if (splitReads != null) {
        srwriter = new BufferedWriter(new FileWriter(splitReads));
    }
    BufferedWriter spwriter = null;
    if (spanningReads != null) {
        spwriter = new BufferedWriter(new FileWriter(spanningReads));
    }
    SAMRecordIterator it = reader.iterator();
    while (it.hasNext()) {
        final SAMRecord r = it.next();
        if (r.isSecondaryOrSupplementary())
            continue;
        if (splitReads != null) {
            for (SimpleBEDFeature f : ChimericAlignment.getChimericAlignments(r).stream().map(ca -> getSplitReadDeletion(r, ca, minLengthPortion, minMapq, minIndelSize)).filter(f -> f != null).collect(Collectors.toList())) {
                writeBed(srwriter, f);
            }
        }
        if (spanningReads != null) {
            for (SimpleBEDFeature f : getSpanningDeletion(r, minMapq, minIndelSize, maxMappedBases, minIndelComponentSize)) {
                writeBed(spwriter, f);
            }
        }
    }
    CloserUtil.close(srwriter);
    CloserUtil.close(spwriter);
}

Example 19 with SamReaderFactory

use of htsjdk.samtools.SamReaderFactory in project jvarkit by lindenb.

the class BamToSVG method doWork.

@Override
public int doWork(List<String> args) {
    /* parse interval */
    if (this.intervalStr == null) {
        LOG.error("bed.interval0.undefined");
        return -1;
    }
    int colon = this.intervalStr.indexOf(':');
    int hyphen = this.intervalStr.indexOf('-', colon + 1);
    if (colon < 1 || hyphen <= colon || hyphen + 1 == intervalStr.length()) {
        LOG.error("Bad interval " + this.intervalStr);
        return -1;
    }
    this.interval = new Interval();
    this.interval.chrom = this.intervalStr.substring(0, colon);
    this.interval.start = Integer.parseInt(this.intervalStr.substring(colon + 1, hyphen)) + 1;
    this.interval.end = Integer.parseInt(this.intervalStr.substring(hyphen + 1));
    this.drawinAreaWidth = Math.max(100, this.drawinAreaWidth);
    SamReader in = null;
    SAMRecordIterator iter = null;
    SamReaderFactory sfrf = SamReaderFactory.makeDefault();
    sfrf.validationStringency(ValidationStringency.SILENT);
    XMLStreamWriter w = null;
    FileOutputStream fout = null;
    try {
        /* get genomic sequence */
        if (this.referenceFile != null) {
            LOG.info("opening " + this.referenceFile);
            this.indexedFastaSequenceFile = new IndexedFastaSequenceFile(this.referenceFile);
            this.genomicSequence = new GenomicSequence(this.indexedFastaSequenceFile, this.interval.chrom);
        }
        for (String vcf : this.vcfFileSet) {
            readVariantFile(vcf);
        }
        /* read SAM data */
        if (args.isEmpty()) {
            LOG.info("Reading from stdin");
            in = sfrf.open(SamInputResource.of(stdin()));
            iter = in.iterator();
            readBamStream(iter);
            iter.close();
            in.close();
        } else {
            for (String arg : args) {
                File filename = new File(arg);
                LOG.info("Reading from " + filename);
                in = sfrf.open(SamInputResource.of(filename));
                if (in.hasIndex()) {
                    iter = in.query(this.interval.getChrom(), this.interval.getStart(), this.interval.getEnd(), false);
                } else {
                    LOG.info("Bam file not indexed !! " + filename);
                    iter = in.iterator();
                }
                readBamStream(iter);
                iter.close();
                in.close();
            }
        }
        this.featureWidth = this.drawinAreaWidth / (double) ((this.interval.end - this.interval.start) + 1);
        this.featureHeight = Math.min(Math.max(5.0, this.featureWidth), 30);
        this.HEIGHT_RULER = (int) (this.niceIntFormat.format(this.interval.end).length() * this.featureHeight + 5);
        LOG.info("Feature height:" + this.featureHeight);
        XMLOutputFactory xof = XMLOutputFactory.newFactory();
        if (this.outputFile == null) {
            w = xof.createXMLStreamWriter(stdout(), "UTF-8");
        } else {
            fout = new FileOutputStream(this.outputFile);
            w = xof.createXMLStreamWriter(fout, "UTF-8");
        }
        w.writeStartDocument("UTF-8", "1.0");
        printDocument(w, intervalStr);
        w.writeEndDocument();
        w.flush();
        w.close();
        return RETURN_OK;
    } catch (Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(iter);
        CloserUtil.close(fout);
        CloserUtil.close(in);
        CloserUtil.close(this.indexedFastaSequenceFile);
        this.indexedFastaSequenceFile = null;
        this.interval = null;
    }
}

Example 20 with SamReaderFactory

use of htsjdk.samtools.SamReaderFactory in project jvarkit by lindenb.

the class Bam2Wig method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.win_shift <= 0) {
        LOG.error("window shift<=0");
        return -1;
    }
    if (this.window_span <= 0) {
        LOG.error("window size<=0");
        return -1;
    }
    final Interval interval;
    PrintWriter pw = null;
    CloseableIterator<SAMRecord> samRecordIterator = null;
    final List<SamReader> samReaders = new ArrayList<>();
    final List<CloseableIterator<SAMRecord>> merginIterators = new ArrayList<>();
    try {
        final SamReaderFactory srf = SamReaderFactory.makeDefault().validationStringency(htsjdk.samtools.ValidationStringency.LENIENT);
        if (args.isEmpty()) {
            if (!StringUtil.isBlank(region_str)) {
                LOG.error("cannot specify region for stdin");
                return -1;
            }
            interval = null;
            samReaders.add(srf.open(SamInputResource.of(stdin())));
            samRecordIterator = samReaders.get(0).iterator();
        } else if (args.size() == 1 && !args.get(0).endsWith(".list")) {
            samReaders.add(srf.open(SamInputResource.of(new File(args.get(0)))));
            if (StringUtil.isBlank(this.region_str)) {
                samRecordIterator = samReaders.get(0).iterator();
                interval = null;
            } else {
                interval = new IntervalParser(samReaders.get(0).getFileHeader().getSequenceDictionary()).setContigNameIsWholeContig(true).parse(region_str);
                if (interval == null) {
                    LOG.error("Cannot parse interval " + this.region_str);
                    return -1;
                }
                LOG.debug("interval " + interval);
                samRecordIterator = samReaders.get(0).query(interval.getContig(), interval.getStart(), interval.getEnd(), false);
            }
        } else {
            final List<File> samFiles;
            if (args.size() == 1 && args.get(0).endsWith(".list")) {
                samFiles = IOUtils.unrollFile(new File(args.get(0)));
            } else {
                samFiles = args.stream().map(S -> new File(S)).collect(Collectors.toList());
            }
            if (samFiles.isEmpty()) {
                LOG.error("No Input SAM file");
                return -1;
            }
            final SAMSequenceDictionary dict0 = SAMSequenceDictionaryExtractor.extractDictionary(samFiles.get(0));
            if (dict0 == null)
                throw new JvarkitException.DictionaryMissing(samFiles.get(0).getPath());
            samFiles.stream().forEach(F -> {
                final SAMSequenceDictionary dicti = SAMSequenceDictionaryExtractor.extractDictionary(F);
                if (dicti == null)
                    throw new JvarkitException.DictionaryMissing(F.getPath());
                if (!SequenceUtil.areSequenceDictionariesEqual(dicti, dict0)) {
                    throw new JvarkitException.DictionariesAreNotTheSame(dict0, dicti);
                }
            });
            for (final File bamFile : samFiles) {
                LOG.info("opening " + bamFile);
                samReaders.add(srf.open(bamFile));
            }
            final SamFileHeaderMerger mergedheader = new SamFileHeaderMerger(SAMFileHeader.SortOrder.coordinate, samReaders.stream().map(SR -> SR.getFileHeader()).collect(Collectors.toList()), false);
            final Map<SamReader, CloseableIterator<SAMRecord>> reader2iter = new HashMap<>();
            if (StringUtil.isBlank(this.region_str)) {
                for (final SamReader sr : samReaders) {
                    reader2iter.put(sr, sr.iterator());
                }
                interval = null;
            } else {
                interval = new IntervalParser(dict0).setContigNameIsWholeContig(true).parse(region_str);
                if (interval == null) {
                    LOG.error("Cannot parse interval " + this.region_str);
                    return -1;
                }
                LOG.info("interval :" + interval);
                for (final SamReader sr : samReaders) {
                    reader2iter.put(sr, sr.query(interval.getContig(), interval.getStart(), interval.getEnd(), false));
                }
            }
            merginIterators.addAll(reader2iter.values());
            samRecordIterator = new MergingSamRecordIterator(mergedheader, reader2iter, true);
        }
        for (final SamReader sr : samReaders) {
            if (sr.getFileHeader().getSortOrder() != SAMFileHeader.SortOrder.coordinate) {
                LOG.error("one of your bam input is not sorted on coordinate");
                return -1;
            }
        }
        pw = openFileOrStdoutAsPrintWriter(this.outputFile);
        run(pw, samRecordIterator, samReaders.get(0).getFileHeader().getSequenceDictionary(), interval);
        samRecordIterator.close();
        samRecordIterator = null;
        CloserUtil.close(samReaders);
        samReaders.clear();
        pw.flush();
        return RETURN_OK;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(merginIterators);
        CloserUtil.close(samRecordIterator);
        CloserUtil.close(samReaders);
        CloserUtil.close(pw);
        pw = null;
    }
}