Examples with ReadClipper com.github.lindenb.jvarkit.tools.pcr.ReadClipper used on opensource projects

Example 1 with ReadClipper

use of com.github.lindenb.jvarkit.tools.pcr.ReadClipper in project jvarkit by lindenb.

the class Biostar480685 method doWork.

@Override
public int doWork(final List<String> args) {
    SamReader in = null;
    SAMFileWriter out = null;
    try {
        final SamReaderFactory srf = super.createSamReaderFactory();
        if (this.faidx != null) {
            srf.referenceSequence(this.faidx);
            writingBamArgs.setReferencePath(this.faidx);
        }
        final String input = oneFileOrNull(args);
        if (input == null) {
            in = srf.open(SamInputResource.of(stdin()));
        } else {
            in = srf.open(SamInputResource.of(input));
        }
        final SAMFileHeader header = in.getFileHeader();
        if (!(header.getSortOrder().equals(SAMFileHeader.SortOrder.unsorted) || header.getSortOrder().equals(SAMFileHeader.SortOrder.queryname))) {
            LOG.error("input should be sorted with 'samtools sort -n' or 'samtools collate' but got " + header.getSortOrder());
            return -1;
        }
        final ReadClipper clipper = new ReadClipper();
        header.setSortOrder(SAMFileHeader.SortOrder.unsorted);
        final SAMProgramRecord prg = header.createProgramRecord();
        prg.setCommandLine(this.getProgramCommandLine());
        prg.setProgramName(this.getProgramName());
        prg.setProgramVersion(this.getGitHash());
        JVarkitVersion.getInstance().addMetaData(this, header);
        out = this.writingBamArgs.openSamWriter(this.outputFile, header, true);
        try (CloseableIterator<List<SAMRecord>> iter = new EqualIterator<>(in.iterator(), (A, B) -> A.getReadName().compareTo(B.getReadName()))) {
            while (iter.hasNext()) {
                final List<SAMRecord> buffer = iter.next();
                int read1_idx = -1;
                int read2_idx = -1;
                for (int i = 0; i < buffer.size(); i++) {
                    final SAMRecord rec = buffer.get(i);
                    if (!rec.getReadPairedFlag())
                        continue;
                    if (rec.getReadUnmappedFlag())
                        continue;
                    if (rec.getMateUnmappedFlag())
                        continue;
                    if (rec.isSecondaryOrSupplementary())
                        continue;
                    if (rec.getFirstOfPairFlag()) {
                        read1_idx = i;
                    } else if (rec.getSecondOfPairFlag()) {
                        read2_idx = i;
                    }
                }
                if (read1_idx == -1 || read2_idx == -1 || read1_idx == read2_idx)
                    continue;
                final SAMRecord rec1a = buffer.get(read1_idx);
                final SAMRecord rec2a = buffer.get(read2_idx);
                if (!rec1a.overlaps(rec2a))
                    continue;
                final int chromStart = Math.max(rec1a.getStart(), rec2a.getStart());
                final int chromEnd = Math.min(rec1a.getEnd(), rec2a.getEnd());
                if (chromStart > chromEnd)
                    continue;
                final SimpleInterval rgn = new SimpleInterval(rec1a.getContig(), chromStart, chromEnd);
                final SAMRecord rec1b = clipper.clip(rec1a, rgn);
                if (rec1b == null || rec1b.getReadUnmappedFlag())
                    continue;
                final SAMRecord rec2b = clipper.clip(rec2a, rgn);
                if (rec2b == null || rec2b.getReadUnmappedFlag())
                    continue;
                rec1b.setAttribute("PG", prg.getId());
                rec2b.setAttribute("PG", prg.getId());
                rec1b.setAlignmentStart(chromStart);
                rec1b.setMateAlignmentStart(rec2b.getAlignmentStart());
                rec2b.setAlignmentStart(chromStart);
                rec2b.setMateAlignmentStart(rec1b.getAlignmentStart());
                rec1b.setAttribute("MC", rec2b.getCigarString());
                rec2b.setAttribute("MC", rec1b.getCigarString());
                rec1b.setAttribute("NM", null);
                rec2b.setAttribute("NM", null);
                buffer.set(read1_idx, rec1b);
                buffer.set(read2_idx, rec2b);
                for (SAMRecord rec : buffer) {
                    out.addAlignment(rec);
                }
            }
        }
        in.close();
        in = null;
        out.close();
        out = null;
        return 0;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(in);
        CloserUtil.close(out);
    }
}