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

Example 51 with SAMRecordIterator

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

the class SamGrep method doWork.

@Override
public int doWork(List<String> args) {
    readNames.clear();
    if (namefile != null) {
        BufferedReader in = null;
        try {
            in = IOUtils.openFileForBufferedReading(this.namefile);
            String line;
            while ((line = in.readLine()) != null) {
                line = line.trim();
                if (line.isEmpty())
                    continue;
                readNames.put(line, 0);
            }
        } catch (Exception err) {
            LOG.error(err);
            return -1;
        } finally {
            CloserUtil.close(in);
        }
    }
    for (final String line : this.nameStrings) {
        readNames.put(line, 0);
    }
    if (readNames.isEmpty()) {
        LOG.warn("no read found.");
    }
    SAMFileWriter sfw = null;
    SAMFileWriter samStdout = null;
    SamReader sfr = null;
    try {
        sfr = super.openSamReader(oneFileOrNull(args));
        final SAMFileHeader header = sfr.getFileHeader().clone();
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(header);
        if (this.outputFile == null) {
            sfw = this.writingBamArgs.openSAMFileWriter(null, header, true);
        } else {
            samStdout = this.writingBamArgs.openSAMFileWriter(null, header, true);
            sfw = this.writingBamArgs.openSAMFileWriter(outputFile, header, true);
        }
        SAMRecordIterator iter = sfr.iterator();
        while (iter.hasNext()) {
            boolean keep = false;
            final SAMRecord rec = progress.watch(iter.next());
            if (samStdout != null)
                samStdout.addAlignment(rec);
            Integer count = readNames.get(rec.getReadName());
            if (count != null) {
                keep = true;
            }
            if (this.inverse)
                keep = !keep;
            if (keep) {
                sfw.addAlignment(rec);
            }
            if (n_before_remove != -1 && !inverse && keep) {
                count++;
                if (count >= n_before_remove) {
                    readNames.remove(rec.getReadName());
                    if (samStdout == null && readNames.isEmpty())
                        break;
                } else {
                    readNames.put(rec.getReadName(), count);
                }
            }
        }
        progress.finish();
        return RETURN_OK;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(samStdout);
        CloserUtil.close(sfw);
        CloserUtil.close(sfr);
    }
}

Example 52 with SAMRecordIterator

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

the class SamCustomSortJdk method doWork.

@Override
public int doWork(final List<String> args) {
    SAMRecordIterator iter = null;
    SamReader samFileReader = null;
    SAMFileWriter sw = null;
    SortingCollection<SAMRecord> sorter = null;
    CloseableIterator<SAMRecord> iter2 = null;
    try {
        final String code;
        if (this.scriptFile != null) {
            code = IOUtil.slurp(this.scriptFile);
        } else if (!StringUtil.isBlank(this.scriptExpr)) {
            code = this.scriptExpr;
        } else {
            LOG.error("Option -e or -f are required. The content of those empty mut be not empty");
            return -1;
        }
        final Random rand = new Random(System.currentTimeMillis());
        final String javaClassName = SamCustomSortJdk.class.getSimpleName() + "Custom" + Math.abs(rand.nextInt());
        final StringWriter codeWriter = new StringWriter();
        final PrintWriter pw = new PrintWriter(codeWriter);
        pw.println("import java.util.*;");
        pw.println("import java.util.stream.*;");
        pw.println("import java.util.function.*;");
        pw.println("import htsjdk.samtools.*;");
        pw.println("import htsjdk.samtools.util.*;");
        pw.println("import com.github.lindenb.jvarkit.tools.misc.IlluminaReadName;");
        pw.println("import javax.annotation.Generated;");
        pw.println("@Generated(value=\"" + SamCustomSortJdk.class.getSimpleName() + "\",date=\"" + new Iso8601Date(new Date()) + "\")");
        pw.println("public class " + javaClassName + " extends " + AbstractSamComparator.class.getName().replace('$', '.') + " {");
        pw.println("  public " + javaClassName + "(final SAMFileHeader header) {");
        pw.println("  super(header);");
        pw.println("  }");
        if (user_code_is_body) {
            pw.println("   //user's code starts here");
            pw.println(code);
            pw.println("   // user's code ends here");
        } else {
            pw.println("  @Override");
            pw.println("  public int compare(final  SAMRecord R1,final SAMRecord R2) {");
            pw.println("   /** user's code starts here */");
            pw.println(code);
            pw.println("/** user's code ends here */");
            pw.println("   }");
        }
        pw.println("}");
        pw.flush();
        if (!hideGeneratedCode) {
            LOG.debug(" Compiling :\n" + InMemoryCompiler.beautifyCode(codeWriter.toString()));
        }
        if (this.saveCodeInDir != null) {
            PrintWriter cw = null;
            try {
                IOUtil.assertDirectoryIsWritable(this.saveCodeInDir);
                cw = new PrintWriter(new File(this.saveCodeInDir, javaClassName + ".java"));
                cw.write(codeWriter.toString());
                cw.flush();
                cw.close();
                cw = null;
                LOG.info("saved " + javaClassName + ".java in " + this.saveCodeInDir);
            } catch (final Exception err) {
                LOG.error(err);
                return -1;
            } finally {
                CloserUtil.close(cw);
            }
        }
        final InMemoryCompiler inMemoryCompiler = new InMemoryCompiler();
        final Class<?> compiledClass = inMemoryCompiler.compileClass(javaClassName, codeWriter.toString());
        final Constructor<?> ctor = compiledClass.getDeclaredConstructor(SAMFileHeader.class);
        samFileReader = openSamReader(oneFileOrNull(args));
        final SAMFileHeader headerIn = samFileReader.getFileHeader();
        final StableSort customComparator = new StableSort((Comparator<SAMRecord>) ctor.newInstance(headerIn));
        final BAMRecordCodec bamRecordCodec = new BAMRecordCodec(headerIn);
        sorter = SortingCollection.newInstance(SAMRecord.class, bamRecordCodec, customComparator, this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
        sorter.setDestructiveIteration(true);
        SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(headerIn).logger(LOG);
        iter = samFileReader.iterator();
        while (iter.hasNext()) {
            sorter.add(progress.watch(iter.next()));
        }
        samFileReader.close();
        samFileReader = null;
        sorter.doneAdding();
        final SAMFileHeader headerOut = headerIn.clone();
        headerOut.setSortOrder(SAMFileHeader.SortOrder.unsorted);
        headerOut.addComment(getProgramName() + " " + getVersion() + " " + getProgramCommandLine());
        sw = this.writingBamArgs.openSAMFileWriter(this.outputFile, headerOut, false);
        progress = new SAMSequenceDictionaryProgress(headerIn).logger(LOG);
        iter2 = sorter.iterator();
        while (iter2.hasNext()) {
            sw.addAlignment(progress.watch(iter2.next()));
        }
        iter2.close();
        iter2 = null;
        sw.close();
        sw = null;
        progress.finish();
        return RETURN_OK;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        try {
            if (sorter != null)
                sorter.cleanup();
        } catch (Exception e) {
        }
        CloserUtil.close(iter);
        CloserUtil.close(iter2);
        CloserUtil.close(samFileReader);
        CloserUtil.close(sw);
    }
}

Example 53 with SAMRecordIterator

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

the class SamJavascript method doWork.

@Override
public int doWork(final List<String> args) {
    SAMRecordIterator iter = null;
    SamReader samFileReader = null;
    SAMFileWriter sw = null;
    try {
        this.script = super.compileJavascript(this.jsExpression, this.jsFile);
        samFileReader = openSamReader(oneFileOrNull(args));
        final SAMFileHeader header = samFileReader.getFileHeader();
        sw = writingBamArgs.openSAMFileWriter(outputFile, header, true);
        long count = 0L;
        final Bindings bindings = this.script.getEngine().createBindings();
        bindings.put("header", samFileReader.getFileHeader());
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(header).logger(LOG);
        iter = samFileReader.iterator();
        while (iter.hasNext()) {
            final SAMRecord record = iter.next();
            progress.watch(record);
            bindings.put("record", record);
            if (super.evalJavaScriptBoolean(this.script, bindings)) {
                ++count;
                sw.addAlignment(record);
                if (this.LIMIT > 0L && count >= this.LIMIT)
                    break;
            } else {
                failing(record, header);
            }
        }
        sw.close();
        /* create empty if never called */
        openFailing(header);
        return RETURN_OK;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(iter);
        CloserUtil.close(samFileReader);
        CloserUtil.close(sw);
        CloserUtil.close(failingReadsWriter);
    }
}

Example 54 with SAMRecordIterator

use of htsjdk.samtools.SAMRecordIterator in project hmftools by hartwigmedical.

the class MNVValidator method mergeVariants.

@NotNull
public List<VariantContext> mergeVariants(@NotNull final PotentialMNVRegion potentialMnvRegion) {
    if (potentialMnvRegion.potentialMnvs().size() == 0) {
        return potentialMnvRegion.variants();
    } else {
        final SAMRecordIterator samIterator = queryBam(potentialMnvRegion);
        final MNVRegionValidator regionValidator = validateMNVs(samIterator, potentialMnvRegion);
        samIterator.close();
        return outputVariants(regionValidator);
    }
}

Example 55 with SAMRecordIterator

use of htsjdk.samtools.SAMRecordIterator in project hmftools by hartwigmedical.

the class Analysis method determineBreakpointsImprecise.

private BreakpointResult determineBreakpointsImprecise(final HMFVariantContext ctx, final SamReader reader) {
    final Pair<Integer, Integer> ctxOrientation = Pair.of(ctx.OrientationBP1, ctx.OrientationBP2);
    final PairedReads interesting = new PairedReads();
    final PairedReads clipped_proper = new PairedReads();
    final PairedReads secondary_pairs = new PairedReads();
    final List<SAMRecord> currentReads = Lists.newArrayList();
    final SAMRecordIterator iterator = reader.iterator();
    while (iterator.hasNext() || !currentReads.isEmpty()) {
        final SAMRecord record = iterator.hasNext() ? iterator.next() : null;
        if (record != null) {
            if (currentReads.isEmpty() || record.getReadName().equals(currentReads.get(0).getReadName())) {
                currentReads.add(record);
                continue;
            }
        }
        currentReads.sort(comparator);
        final PairedReads pairs = pairs(currentReads);
        currentReads.clear();
        if (record != null) {
            currentReads.add(record);
        }
        for (final Pair<SAMRecord, SAMRecord> pair : pairs) {
            final boolean correctOrientation = orientation(pair).equals(ctxOrientation);
            final boolean correctChromosome = Location.fromSAMRecord(pair.getLeft()).sameChromosomeAs(ctx.MantaBP1) && Location.fromSAMRecord(pair.getRight()).sameChromosomeAs(ctx.MantaBP2);
            final boolean hasExpectedClipping = clippedOnCorrectSide(pair.getLeft(), ctx.OrientationBP1) || clippedOnCorrectSide(pair.getRight(), ctx.OrientationBP2);
            final boolean sameChromosome = pair.getLeft().getReferenceIndex().equals(pair.getRight().getReferenceIndex());
            final boolean potentialSROnly = sameChromosome && Stream.of(ctx.OrientationBP1, ctx.OrientationBP2).anyMatch(orientation -> clippedOnCorrectSide(orientation > 0 ? pair.getRight() : pair.getLeft(), orientation));
            final boolean secondary = stream(pair).anyMatch(SAMRecord::isSecondaryOrSupplementary);
            final boolean proper = stream(pair).anyMatch(SAMRecord::getProperPairFlag) && !secondary;
            LOGGER.trace("determineBreakpoints {} {}->{} {} {}->{} correctOrientation({}) correctChromosome({}) hasExpectedClipping({}) proper({}) potentialSROnly({}) secondary({})", pair.getLeft().getReadName(), pair.getLeft().getAlignmentStart(), pair.getLeft().getMateAlignmentStart(), pair.getRight().getReadName(), pair.getRight().getAlignmentStart(), pair.getRight().getMateAlignmentStart(), correctOrientation, correctChromosome, hasExpectedClipping, proper, potentialSROnly, secondary);
            if (secondary && potentialSROnly) {
                secondary_pairs.add(pair);
            } else if ((!proper || hasExpectedClipping) && correctChromosome && correctOrientation) {
                interesting.add(pair);
            } else if (proper && potentialSROnly) {
                clipped_proper.add(pair);
            }
        }
    }
    iterator.close();
    // load clipping info
    Clipping bp1_clipping = new Clipping();
    Clipping bp2_clipping = new Clipping();
    for (final Pair<SAMRecord, SAMRecord> pair : interesting) {
        if (ctx.OrientationBP1 > 0) {
            bp1_clipping.add(Clipping.getRightClip(pair.getLeft()));
        } else {
            bp1_clipping.add(Clipping.getLeftClip(pair.getLeft()));
        }
        if (ctx.OrientationBP2 > 0) {
            bp2_clipping.add(Clipping.getRightClip(pair.getRight()));
        } else {
            bp2_clipping.add(Clipping.getLeftClip(pair.getRight()));
        }
    }
    for (final Pair<SAMRecord, SAMRecord> pair : clipped_proper) {
        if (stream(pair).allMatch(r -> Location.fromSAMRecord(r).sameChromosomeAs(ctx.MantaBP1))) {
            if (ctx.OrientationBP1 > 0) {
                bp1_clipping.add(Clipping.getRightClip(pair.getRight()));
            } else {
                bp1_clipping.add(Clipping.getLeftClip(pair.getLeft()));
            }
        }
        if (stream(pair).allMatch(r -> Location.fromSAMRecord(r).sameChromosomeAs(ctx.MantaBP2))) {
            if (ctx.OrientationBP2 > 0) {
                bp2_clipping.add(Clipping.getRightClip(pair.getRight()));
            } else {
                bp2_clipping.add(Clipping.getLeftClip(pair.getLeft()));
            }
        }
    }
    for (final Pair<SAMRecord, SAMRecord> pair : secondary_pairs) {
        if (stream(pair).allMatch(r -> Location.fromSAMRecord(r).sameChromosomeAs(ctx.MantaBP1))) {
            if (ctx.OrientationBP1 > 0) {
                bp1_clipping.add(Clipping.getRightClip(pair.getRight()));
            } else {
                bp1_clipping.add(Clipping.getLeftClip(pair.getLeft()));
            }
        }
        if (stream(pair).allMatch(r -> Location.fromSAMRecord(r).sameChromosomeAs(ctx.MantaBP2))) {
            if (ctx.OrientationBP2 > 0) {
                bp2_clipping.add(Clipping.getRightClip(pair.getRight()));
            } else {
                bp2_clipping.add(Clipping.getLeftClip(pair.getLeft()));
            }
        }
    }
    // determinate candidates based on clipping info
    final List<Location> bp1_candidates = bp1_clipping.getSequences().stream().map(c -> c.Alignment).filter(c -> withinRange(c, ctx.MantaBP1, ctx.Uncertainty1)).collect(Collectors.toList());
    if (bp1_candidates.isEmpty()) {
        final Location candidate = interesting.stream().map(Pair::getLeft).map(r -> Location.fromSAMRecord(r, ctx.OrientationBP1 < 0).add(ctx.OrientationBP1 > 0 ? 1 : -1)).filter(l -> withinRange(l, ctx.MantaBP1, ctx.Uncertainty1)).max((a, b) -> ctx.OrientationBP1 > 0 ? a.compareTo(b) : b.compareTo(a)).orElse(null);
        if (candidate != null) {
            bp1_candidates.add(candidate);
        }
    }
    final List<Location> bp2_candidates = bp2_clipping.getSequences().stream().map(c -> c.Alignment).filter(c -> withinRange(c, ctx.MantaBP2, ctx.Uncertainty2)).collect(Collectors.toList());
    if (bp2_candidates.isEmpty()) {
        final Location candidate = interesting.stream().map(Pair::getRight).map(r -> Location.fromSAMRecord(r, ctx.OrientationBP2 < 0).add(ctx.OrientationBP2 > 0 ? 1 : -1)).filter(l -> withinRange(l, ctx.MantaBP2, ctx.Uncertainty2)).max((a, b) -> ctx.OrientationBP2 > 0 ? a.compareTo(b) : b.compareTo(a)).orElse(null);
        if (candidate != null) {
            bp2_candidates.add(candidate);
        }
    }
    // NOTE: we include homology on both sides here and take it out later
    LOGGER.trace("bp1_candidates={} bp2_candidates={}", bp1_candidates, bp2_candidates);
    final Location breakpoint1 = bp1_candidates.isEmpty() ? null : bp1_candidates.get(0).add(-ctx.OrientationBP1);
    final Location breakpoint2 = bp2_candidates.isEmpty() ? null : bp2_candidates.get(0).add(-ctx.OrientationBP2);
    return BreakpointResult.from(Pair.of(breakpoint1, breakpoint2));
}