Search in sources :

Example 11 with SAMRecordIterator

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

the class SequentialCoverageAnnotator method createLookup.

private List<ReferenceCoverageLookup> createLookup(ProcessingContext context, List<SAMEvidenceSource> sources, int windowSize) {
    List<ReferenceCoverageLookup> result = new ArrayList<>();
    for (SAMEvidenceSource ses : sources) {
        assert (ses.getSourceCategory() >= 0);
        assert (ses.getSourceCategory() < context.getCategoryCount());
        // one read-ahead thread per input file
        SamReader reader = SamReaderFactory.makeDefault().open(ses.getFile());
        SAMRecordIterator rawIterator = reader.iterator();
        rawIterator.assertSorted(SortOrder.coordinate);
        CloseableIterator<SAMRecord> sit = new AsyncBufferedIterator<SAMRecord>(rawIterator, ses.getFile().getName() + "-Coverage");
        // close the async iterator first to prevent aysnc reading from a closed stream
        toclose.add(sit);
        toclose.add(rawIterator);
        toclose.add(reader);
        sit = new ProgressLoggingSAMRecordIterator(sit, new ProgressLogger(log, 10000000));
        SequentialReferenceCoverageLookup sourceLookup = new SequentialReferenceCoverageLookup(sit, ses.getMetrics().getIdsvMetrics(), ses.getReadPairConcordanceCalculator(), windowSize, ses.getSourceCategory(), context.isFilterDuplicates());
        context.registerBuffer(ses.getFile().getName(), sourceLookup);
        result.add(sourceLookup);
    }
    return result;
}
Also used : SamReader(htsjdk.samtools.SamReader) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) SAMRecord(htsjdk.samtools.SAMRecord) ArrayList(java.util.ArrayList) AsyncBufferedIterator(au.edu.wehi.idsv.util.AsyncBufferedIterator) ProgressLogger(htsjdk.samtools.util.ProgressLogger)

Example 12 with SAMRecordIterator

use of htsjdk.samtools.SAMRecordIterator 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;
}
Also used : SamReader(htsjdk.samtools.SamReader) SamReaderFactory(htsjdk.samtools.SamReaderFactory) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) SAMFileWriter(htsjdk.samtools.SAMFileWriter) SAMFileWriterFactory(htsjdk.samtools.SAMFileWriterFactory) IOException(java.io.IOException) SAMFileHeader(htsjdk.samtools.SAMFileHeader) File(java.io.File)

Example 13 with SAMRecordIterator

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

the class ByReadNameSinglePassSamProgram method makeItSo.

public static void makeItSo(final File input, final File referenceSequence, final boolean assumeSorted, final long stopAfter, final Collection<ByReadNameSinglePassSamProgram> programs) throws FileNotFoundException {
    // Setup the standard inputs
    IOUtil.assertFileIsReadable(input);
    SamReader in = SamReaderFactory.makeDefault().referenceSequence(referenceSequence).open(input);
    // Optionally load up the reference sequence and double check sequence dictionaries
    final ReferenceLookup lookup;
    if (referenceSequence == null) {
        lookup = null;
    } else {
        IOUtil.assertFileIsReadable(referenceSequence);
        lookup = new TwoBitBufferedReferenceSequenceFile(new IndexedFastaSequenceFile(referenceSequence));
        if (!in.getFileHeader().getSequenceDictionary().isEmpty()) {
            SequenceUtil.assertSequenceDictionariesEqual(in.getFileHeader().getSequenceDictionary(), lookup.getSequenceDictionary());
        }
    }
    // Check on the sort order of the BAM file
    final SortOrder sort = in.getFileHeader().getSortOrder();
    if (sort != SortOrder.queryname) {
        if (assumeSorted) {
            log.warn("File reports sort order '" + sort + "', assuming it's queryname sorted anyway.");
        } else {
            throw new PicardException("File " + input.getAbsolutePath() + " should be queryname sorted but " + "the header says the sort order is " + sort + ". If you believe the file " + "to be queryname sorted you may pass ASSUME_SORTED=true");
        }
    }
    for (final ByReadNameSinglePassSamProgram program : programs) {
        program.setReference(lookup);
        program.setup(in.getFileHeader(), input);
    }
    final ProgressLogger progress = new ProgressLogger(log);
    final SAMRecordIterator rawit = in.iterator();
    final CloseableIterator<SAMRecord> it = new AsyncBufferedIterator<SAMRecord>(rawit, "ByReadNameSinglePassSamProgram " + input.getName());
    try {
        List<SAMRecord> currentRecords = new ArrayList<>();
        String currentReadName = null;
        while (it.hasNext()) {
            SAMRecord r = it.next();
            String readname = r.getReadName();
            // if read name we have to just treat it as a single read
            if (readname == null || !readname.equals(currentReadName)) {
                if (currentRecords.size() > 0) {
                    for (final ByReadNameSinglePassSamProgram program : programs) {
                        program.acceptFragment(currentRecords, lookup);
                    }
                }
                currentRecords.clear();
                currentReadName = readname;
                if (stopAfter > 0 && progress.getCount() >= stopAfter) {
                    break;
                }
            }
            currentRecords.add(r);
            progress.record(r);
        }
        if (currentRecords.size() > 0) {
            for (final ByReadNameSinglePassSamProgram program : programs) {
                program.acceptFragment(currentRecords, lookup);
            }
        }
    } finally {
        CloserUtil.close(it);
        CloserUtil.close(rawit);
        CloserUtil.close(in);
    }
    for (final ByReadNameSinglePassSamProgram program : programs) {
        program.finish();
    }
}
Also used : TwoBitBufferedReferenceSequenceFile(au.edu.wehi.idsv.picard.TwoBitBufferedReferenceSequenceFile) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) ArrayList(java.util.ArrayList) SortOrder(htsjdk.samtools.SAMFileHeader.SortOrder) AsyncBufferedIterator(au.edu.wehi.idsv.util.AsyncBufferedIterator) ProgressLogger(htsjdk.samtools.util.ProgressLogger) PicardException(picard.PicardException) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile) SamReader(htsjdk.samtools.SamReader) ReferenceLookup(au.edu.wehi.idsv.picard.ReferenceLookup) SAMRecord(htsjdk.samtools.SAMRecord)

Example 14 with SAMRecordIterator

use of htsjdk.samtools.SAMRecordIterator 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);
}
Also used : CommandLineProgramProperties(org.broadinstitute.barclay.argparser.CommandLineProgramProperties) BufferedWriter(java.io.BufferedWriter) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) Argument(org.broadinstitute.barclay.argparser.Argument) CigarElement(htsjdk.samtools.CigarElement) FileWriter(java.io.FileWriter) CigarOperator(htsjdk.samtools.CigarOperator) IOException(java.io.IOException) SamReader(htsjdk.samtools.SamReader) Collectors(java.util.stream.Collectors) ValidationStringency(htsjdk.samtools.ValidationStringency) File(java.io.File) SAMRecord(htsjdk.samtools.SAMRecord) StandardOptionDefinitions(picard.cmdline.StandardOptionDefinitions) ArrayList(java.util.ArrayList) SimpleBEDFeature(htsjdk.tribble.bed.SimpleBEDFeature) CommandLineProgram(picard.cmdline.CommandLineProgram) List(java.util.List) Writer(java.io.Writer) CigarUtil(au.edu.wehi.idsv.sam.CigarUtil) SamReaderFactory(htsjdk.samtools.SamReaderFactory) ChimericAlignment(au.edu.wehi.idsv.sam.ChimericAlignment) CloserUtil(htsjdk.samtools.util.CloserUtil) SamReader(htsjdk.samtools.SamReader) SamReaderFactory(htsjdk.samtools.SamReaderFactory) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) SimpleBEDFeature(htsjdk.tribble.bed.SimpleBEDFeature) FileWriter(java.io.FileWriter) SAMRecord(htsjdk.samtools.SAMRecord) BufferedWriter(java.io.BufferedWriter)

Example 15 with SAMRecordIterator

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

the class LongReadSupportFinder method evaluateDeletion.

/**
 * evaluate the long read support for the given deletion event
 * @param chr
 * @param start1
 * @param end1
 * @param start2
 * @param end2
 * @param length length of deletion
 * @param softClipMargin error margin allowed for soft clip location
 * @param spanningWindowMargin
 * @param minDeletionSize
 * @return
 */
public LongReadSupportLevel evaluateDeletion(String chr, int start1, int end1, int start2, int end2, int softClipMargin, int minSoftClipLength, int spanningWindowMargin, int minDeletionSize) {
    chr = translateReference(chr);
    if (chr == null)
        return null;
    SAMRecordIterator it = null;
    try {
        it = reader.query(chr, start1, end1, false);
        LongReadSupportLevel support = new LongReadSupportLevel();
        while (it.hasNext()) {
            SAMRecord r = it.next();
            if (SAMRecordUtil.getEndSoftClipLength(r) >= minSoftClipLength && IntervalUtil.overlapsClosed(r.getAlignmentEnd(), r.getAlignmentEnd(), start1 - softClipMargin, end1 + softClipMargin)) {
                support.startClipLocations.add(r.getAlignmentEnd());
            }
            int deletions = countDeletions(r, minDeletionSize, start1, end2);
            if (deletions > 0) {
                support.spanningAlignments.add(deletions);
            }
        }
        it.close();
        it = null;
        it = reader.query(chr, start2, end2, false);
        while (it.hasNext()) {
            SAMRecord r = it.next();
            if (SAMRecordUtil.getStartSoftClipLength(r) >= minSoftClipLength && IntervalUtil.overlapsClosed(r.getAlignmentStart(), r.getAlignmentStart(), start2 - softClipMargin, end2 + softClipMargin)) {
                support.endClipLocations.add(r.getAlignmentStart());
            }
        }
        it.close();
        it = null;
        return support;
    } catch (Exception e) {
        e.printStackTrace();
        return null;
    } finally {
        if (it != null) {
            it.close();
        }
    }
}
Also used : SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) SAMRecord(htsjdk.samtools.SAMRecord)

Aggregations

SAMRecordIterator (htsjdk.samtools.SAMRecordIterator)107 SAMRecord (htsjdk.samtools.SAMRecord)92 SamReader (htsjdk.samtools.SamReader)83 SAMFileHeader (htsjdk.samtools.SAMFileHeader)49 SAMSequenceDictionaryProgress (com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress)47 File (java.io.File)47 SAMFileWriter (htsjdk.samtools.SAMFileWriter)45 IOException (java.io.IOException)41 ArrayList (java.util.ArrayList)34 CigarElement (htsjdk.samtools.CigarElement)30 Cigar (htsjdk.samtools.Cigar)26 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)24 SamReaderFactory (htsjdk.samtools.SamReaderFactory)21 SAMReadGroupRecord (htsjdk.samtools.SAMReadGroupRecord)18 CigarOperator (htsjdk.samtools.CigarOperator)16 Interval (htsjdk.samtools.util.Interval)16 PrintWriter (java.io.PrintWriter)15 HashMap (java.util.HashMap)15 SAMFileWriterFactory (htsjdk.samtools.SAMFileWriterFactory)14 SAMSequenceRecord (htsjdk.samtools.SAMSequenceRecord)14