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

Example 76 with SAMFileWriter

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

the class BamQueryReadNames method doWork.

@Override
public int doWork(final List<String> args) {
    PrintWriter notFoundStream = new PrintWriter(new NullOuputStream());
    SamReader sfr = null;
    SAMFileWriter bamw = null;
    try {
        if (!(2 == args.size() || 1 == args.size())) {
            LOG.error("illegal.number.of.arguments");
            return -1;
        }
        if (this.notFoundFile == null) {
            notFoundStream.close();
            notFoundStream = openFileOrStdoutAsPrintWriter(notFoundFile);
        }
        File bamFile = new File(args.get(0));
        sfr = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.SILENT).open(bamFile);
        File nameIdxFile = new File(bamFile.getParentFile(), bamFile.getName() + NAME_IDX_EXTENSION);
        this.indexDef = new NameIndexDef();
        this.raf = new RandomAccessFile(nameIdxFile, "r");
        indexDef.countReads = raf.readLong();
        indexDef.maxNameLengt = raf.readInt();
        LineIterator r = null;
        if (args.size() == 2) {
            r = IOUtils.openURIForLineIterator(args.get(1));
        } else {
            r = IOUtils.openStdinForLineIterator();
        }
        SAMFileHeader header = sfr.getFileHeader().clone();
        bamw = writingBamArgs.openSAMFileWriter(this.outputFile, header, true);
        long iter_start = 0L;
        while (r.hasNext()) {
            String line = r.next();
            String searchRead = null;
            int side = -1;
            if (line.isEmpty() || line.startsWith("#"))
                continue;
            /* forward or reverse is specified ? */
            if (line.endsWith("/1")) {
                side = 1;
                searchRead = line.substring(0, line.length() - 2);
            } else if (line.endsWith("/2")) {
                side = 2;
                searchRead = line.substring(0, line.length() - 2);
            } else {
                side = -1;
                searchRead = line;
            }
            long index = lower_bound(iter_start, this.indexDef.countReads, searchRead);
            if (index >= this.indexDef.countReads) {
                notFoundStream.println(line);
                continue;
            }
            if (query_reads_is_sorted) {
                iter_start = index;
            }
            Set<SAMRecord> found = new LinkedHashSet<SAMRecord>();
            while (index < this.indexDef.countReads) {
                NameAndPos nap = getNameAndPosAt(index);
                if (nap.name.compareTo(searchRead) < 0) {
                    ++index;
                    continue;
                } else if (nap.name.compareTo(searchRead) > 0) {
                    break;
                }
                SAMRecordIterator iter;
                if (nap.tid < 0) {
                    iter = sfr.queryUnmapped();
                } else {
                    iter = sfr.query(header.getSequence(nap.tid).getSequenceName(), nap.pos, 0, true);
                }
                while (iter.hasNext()) {
                    SAMRecord rec = iter.next();
                    if (nap.tid >= 0) {
                        if (nap.tid != rec.getReferenceIndex())
                            throw new IllegalStateException();
                        if (rec.getAlignmentStart() < nap.pos) {
                            continue;
                        }
                        if (rec.getAlignmentStart() > nap.pos) {
                            break;
                        }
                    }
                    if (rec.getReadName().equals(searchRead)) {
                        if (side == 1 && !(rec.getReadPairedFlag() && rec.getFirstOfPairFlag())) {
                            continue;
                        } else if (side == 2 && !(rec.getReadPairedFlag() && rec.getSecondOfPairFlag())) {
                            continue;
                        }
                        found.add(rec);
                    }
                }
                iter.close();
                ++index;
            }
            if (found.isEmpty()) {
                notFoundStream.println(line);
            } else {
                for (SAMRecord rec : found) {
                    bamw.addAlignment(rec);
                }
            }
        }
        CloserUtil.close(r);
        notFoundStream.flush();
        notFoundStream.close();
        notFoundStream = null;
        return 0;
    } catch (Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(notFoundStream);
        CloserUtil.close(raf);
        CloserUtil.close(sfr);
        CloserUtil.close(bamw);
    }
}