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

Example 51 with SamReaderFactory

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

the class CompareBams4 method doWork.

@Override
public int doWork(final List<String> args) {
    PrintWriter out = null;
    final SamReader[] samFileReaders = new SamReader[] { null, null };
    final SAMSequenceDictionary[] dicts = new SAMSequenceDictionary[] { null, null };
    @SuppressWarnings("unchecked") final PeekableIterator<SAMRecord>[] iters = new PeekableIterator[] { null, null };
    final List<List<SAMRecord>> recordLists = Arrays.asList(new ArrayList<SAMRecord>(), new ArrayList<SAMRecord>());
    final Counter<Diff> diffs = new Counter<>();
    final short NM_TAG = SAMTagUtil.getSingleton().NM;
    if (args.size() != 2) {
        LOG.error("Expected two and only two bams please, but got " + args.size());
        return -1;
    }
    try {
        final SamReaderFactory srf = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.LENIENT);
        for (int i = 0; i < args.size() && i < samFileReaders.length; ++i) {
            final File samFile = new File(args.get(i));
            LOG.info("opening " + samFile);
            samFileReaders[i] = srf.open(samFile);
            final SAMFileHeader header = samFileReaders[i].getFileHeader();
            if (header.getSortOrder() != SAMFileHeader.SortOrder.queryname) {
                LOG.error("Expected " + samFile + " to be sorted on " + SAMFileHeader.SortOrder.queryname + " but got " + header.getSortOrder());
                return -1;
            }
            dicts[i] = header.getSequenceDictionary();
            if (dicts[i] == null || dicts[i].isEmpty()) {
                LOG.error("In " + samFile + ": No SAMSequenceDictionary in header.");
                return -1;
            }
            iters[i] = new PeekableIterator<>(samFileReaders[i].iterator());
        }
        if (this.chainFile != null) {
            LOG.info("opening chain file " + this.chainFile + ".");
            this.liftOver = new LiftOver(this.chainFile);
            this.liftOver.setLiftOverMinMatch(this.liftOverMismatch <= 0 ? LiftOver.DEFAULT_LIFTOVER_MINMATCH : this.liftOverMismatch);
            if (!this.disableChainValidation) {
                this.liftOver.validateToSequences(dicts[1]);
            }
        }
        final Comparator<SAMRecord> comparator = this.sortMethod.get();
        for (; ; ) {
            for (int i = 0; i < 2; ++i) {
                if (recordLists.get(i).isEmpty()) {
                    while (iters[i].hasNext()) {
                        final SAMRecord rec = iters[i].peek();
                        if (rec.isSecondaryOrSupplementary()) {
                            iters[i].next();
                            continue;
                        }
                        if (!recordLists.get(i).isEmpty() && comparator.compare(recordLists.get(i).get(0), rec) > 0) {
                            throw new SAMException("Something is wrong in sort order of " + args.get(i) + " : got\n\t" + rec + " " + side(rec) + "\nafter\n\t" + recordLists.get(i).get(0) + " " + side(recordLists.get(i).get(0)) + "\nSee also option (samtools querysort)");
                        } else if (recordLists.get(i).isEmpty() || comparator.compare(recordLists.get(i).get(0), rec) == 0) {
                            recordLists.get(i).add(iters[i].next());
                        } else {
                            break;
                        }
                    }
                }
            }
            if (recordLists.get(0).isEmpty() && recordLists.get(1).isEmpty())
                break;
            if (recordLists.get(0).size() > 1) {
                LOG.warn("size>2 for 1/2:" + recordLists.get(0).get(0).getReadName());
                for (final SAMRecord sr : recordLists.get(0)) {
                    LOG.warn(">> " + sr + " flags:" + sr.getFlags() + " pos:" + sr.getReferenceName() + ":" + sr.getStart());
                }
            }
            if (recordLists.get(1).size() > 1) {
                LOG.warn("size>2 for 2/2:" + recordLists.get(1).get(0).getReadName());
                for (final SAMRecord sr : recordLists.get(1)) {
                    LOG.warn(">> " + sr + " flags:" + sr.getFlags() + " pos:" + sr.getReferenceName() + ":" + sr.getStart());
                }
            }
            final Diff diff = new Diff();
            if ((recordLists.get(0).isEmpty() && !recordLists.get(1).isEmpty()) || (!recordLists.get(0).isEmpty() && !recordLists.get(1).isEmpty() && comparator.compare(recordLists.get(0).get(0), recordLists.get(1).get(0)) > 0)) {
                diff.onlyIn = OnlyIn.ONLY_IN_SECOND;
                recordLists.get(1).clear();
            } else if ((!recordLists.get(0).isEmpty() && recordLists.get(1).isEmpty()) || (!recordLists.get(0).isEmpty() && !recordLists.get(1).isEmpty() && comparator.compare(recordLists.get(0).get(0), recordLists.get(1).get(0)) < 0)) {
                diff.onlyIn = OnlyIn.ONLY_IN_FIRST;
                recordLists.get(0).clear();
            } else {
                final SAMRecord rec0 = recordLists.get(0).get(0);
                final SAMRecord rec1 = recordLists.get(1).get(0);
                final Interval rgn0a = interval(rec0);
                final Interval rgn0b;
                diff.onlyIn = OnlyIn.BOTH;
                diff.diffFlags = new Couple<Integer>(rec0.getFlags(), rec1.getFlags());
                diff.diffFlag = (rec0.getFlags() == rec1.getFlags() ? Flag.Same : Flag.Discordant);
                if (this.liftOver == null) {
                    rgn0b = rgn0a;
                    diff.liftover = LiftOverStatus.NoLiftOver;
                } else if (rgn0a == null) {
                    diff.liftover = LiftOverStatus.SourceUnmapped;
                    rgn0b = rgn0a;
                } else {
                    rgn0b = this.liftOver.liftOver(rgn0a);
                    if (rgn0b == null) {
                        diff.liftover = LiftOverStatus.LiftOverFailed;
                    } else if (dicts[1].getSequence(rgn0b.getContig()) == null) {
                        diff.liftover = LiftOverStatus.DestNotInDict;
                    } else if (rgn0a.getContig().equals(rgn0b.getContig())) {
                        diff.liftover = LiftOverStatus.SameChrom;
                    } else {
                        diff.liftover = LiftOverStatus.DiscordantChrom;
                    }
                }
                final Interval rgn1 = interval(rec1);
                if (rgn0b == null && rgn1 == null) {
                    diff.compareContig = CompareContig.BothUnmapped;
                } else if (rgn0b == null && rgn1 != null) {
                    diff.compareContig = CompareContig.GainMapping;
                    diff.diffChrom = new Couple<String>("*", rgn1.getContig());
                } else if (rgn0b != null && rgn1 == null) {
                    diff.compareContig = CompareContig.LostMapping;
                    diff.diffChrom = new Couple<String>(rgn0b.getContig(), "*");
                } else if (rgn0b.getContig().equals(rgn1.getContig())) {
                    diff.compareContig = CompareContig.SameContig;
                    diff.diffChrom = new Couple<String>(rgn0b.getContig(), rgn1.getContig());
                    final int shift = Math.abs(rgn0b.getStart() - rgn1.getStart());
                    if (shift == 0) {
                        diff.shift = Shift.Zero;
                    } else if (shift <= 10) {
                        diff.shift = Shift.Le10;
                    } else if (shift <= 20) {
                        diff.shift = Shift.Le20;
                    } else if (shift <= 100) {
                        diff.shift = Shift.Le100;
                    } else {
                        diff.shift = Shift.Gt100;
                    }
                } else {
                    if (dicts[1].getSequence(rgn0b.getContig()) == null) {
                        diff.diffChrom = new Couple<String>("?", rgn1.getContig());
                    } else {
                        diff.diffChrom = new Couple<String>(rgn0b.getContig(), rgn1.getContig());
                    }
                    diff.compareContig = CompareContig.DiscordantContig;
                }
                if (rgn0b != null && rgn1 != null) {
                    diff.diffCigarOperations = (rec0.getCigar().numCigarElements() - rec1.getCigar().numCigarElements());
                    final Object nm0 = rec0.getAttribute(NM_TAG);
                    final Object nm1 = rec1.getAttribute(NM_TAG);
                    if (nm0 != null && nm1 != null) {
                        diff.diffNM = (Number.class.cast(nm0).intValue() - Number.class.cast(nm1).intValue());
                    }
                }
                recordLists.get(1).clear();
                recordLists.get(0).clear();
            }
            diffs.incr(diff);
        }
        LOG.info("done");
        final StringBuilder sb = new StringBuilder();
        str(sb, "onlyIn");
        str(sb, "liftover");
        str(sb, "compareContig");
        str(sb, "shift");
        str(sb, "diffCigarOperations");
        str(sb, "diffNM");
        str(sb, "diffFlags");
        str(sb, "diffChroms");
        str(sb, "diffFlag");
        sb.append("Count");
        out = super.openFileOrStdoutAsPrintWriter(outputFile);
        out.println(sb);
        for (final Diff key : diffs.keySet()) {
            out.print(key.toString());
            out.println(diffs.count(key));
        }
        out.flush();
        out.close();
        return RETURN_OK;
    } catch (Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        this.liftOver = null;
        for (int i = 0; i < samFileReaders.length; ++i) {
            CloserUtil.close(iters[i]);
            CloserUtil.close(samFileReaders[i]);
        }
        CloserUtil.close(out);
    }
}

Example 52 with SamReaderFactory

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

the class Biostar78285 method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.gc_percent_window < 1) {
        LOG.error("Bad GC% window size:" + this.gc_percent_window);
        return -1;
    }
    final List<File> bamFiles = IOUtil.unrollFiles(args.stream().map(F -> new File(F)).collect(Collectors.toCollection(HashSet::new)), ".bam");
    SAMSequenceDictionary dict = null;
    final List<SamReader> samReaders = new ArrayList<>();
    final List<CloseableIterator<SAMRecord>> samIterators = new ArrayList<>();
    final TreeSet<String> samples = new TreeSet<>();
    final String DEFAULT_PARTITION = "UNDEFINED_PARTITION";
    IndexedFastaSequenceFile indexedFastaSequenceFile = null;
    VariantContextWriter out = null;
    try {
        final SamReaderFactory samReaderFactory = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.LENIENT);
        for (final File bamFile : bamFiles) {
            LOG.info("Opening " + bamFile);
            final SamReader samReader = samReaderFactory.open(bamFile);
            samReaders.add(samReader);
            final SAMFileHeader header = samReader.getFileHeader();
            if (header == null) {
                LOG.error("No header in " + bamFile);
                return -1;
            }
            if (header.getSortOrder() != SortOrder.coordinate) {
                LOG.error("Sam file " + bamFile + " is not sorted on coordinate :" + header.getSortOrder());
                return -1;
            }
            samples.addAll(header.getReadGroups().stream().map(RG -> this.partition.apply(RG, DEFAULT_PARTITION)).collect(Collectors.toSet()));
            final SAMSequenceDictionary currDict = header.getSequenceDictionary();
            if (currDict == null) {
                LOG.error("SamFile doesn't contain a SAMSequenceDictionary : " + bamFile);
                return -1;
            }
            if (dict == null) {
                dict = currDict;
            } else if (!SequenceUtil.areSequenceDictionariesEqual(dict, currDict)) {
                LOG.error(JvarkitException.DictionariesAreNotTheSame.getMessage(dict, currDict));
                return -1;
            }
        }
        if (samReaders.isEmpty()) {
            LOG.error("no bam");
            return -1;
        }
        if (dict == null) {
            LOG.error("no dictionary");
            return -1;
        }
        final QueryInterval[] intervals;
        if (this.captureBed != null) {
            LOG.info("Opening " + this.captureBed);
            ContigNameConverter.setDefaultAliases(dict);
            final List<QueryInterval> L = new ArrayList<>();
            final BedLineCodec codec = new BedLineCodec();
            final LineIterator li = IOUtils.openFileForLineIterator(this.captureBed);
            while (li.hasNext()) {
                final BedLine bed = codec.decode(li.next());
                if (bed == null)
                    continue;
                final QueryInterval q = bed.toQueryInterval(dict);
                L.add(q);
            }
            CloserUtil.close(li);
            intervals = QueryInterval.optimizeIntervals(L.toArray(new QueryInterval[L.size()]));
        } else {
            intervals = null;
        }
        for (final SamReader samReader : samReaders) {
            LOG.info("querying " + samReader.getResourceDescription());
            final CloseableIterator<SAMRecord> iter;
            if (intervals == null) {
                iter = samReader.iterator();
            } else {
                iter = samReader.queryOverlapping(intervals);
            }
            samIterators.add(new FilterIterator<SAMRecord>(iter, R -> !R.getReadUnmappedFlag() && !filter.filterOut(R)));
        }
        if (this.refFile != null) {
            LOG.info("opening " + refFile);
            indexedFastaSequenceFile = new IndexedFastaSequenceFile(this.refFile);
            final SAMSequenceDictionary refdict = indexedFastaSequenceFile.getSequenceDictionary();
            ContigNameConverter.setDefaultAliases(refdict);
            if (refdict == null) {
                throw new JvarkitException.FastaDictionaryMissing(this.refFile);
            }
            if (!SequenceUtil.areSequenceDictionariesEqual(dict, refdict)) {
                LOG.error(JvarkitException.DictionariesAreNotTheSame.getMessage(dict, refdict));
                return -1;
            }
        }
        out = openVariantContextWriter(this.outputFile);
        final Set<VCFHeaderLine> metaData = new HashSet<>();
        VCFStandardHeaderLines.addStandardFormatLines(metaData, true, VCFConstants.DEPTH_KEY, VCFConstants.GENOTYPE_KEY);
        VCFStandardHeaderLines.addStandardInfoLines(metaData, true, VCFConstants.DEPTH_KEY);
        metaData.add(new VCFFormatHeaderLine("DF", 1, VCFHeaderLineType.Integer, "Number of Reads on plus strand"));
        metaData.add(new VCFFormatHeaderLine("DR", 1, VCFHeaderLineType.Integer, "Number of Reads on minus strand"));
        metaData.add(new VCFInfoHeaderLine("AVG_DP", 1, VCFHeaderLineType.Float, "Mean depth"));
        metaData.add(new VCFInfoHeaderLine("MEDIAN_DP", 1, VCFHeaderLineType.Float, "Median depth"));
        metaData.add(new VCFInfoHeaderLine("MIN_DP", 1, VCFHeaderLineType.Integer, "Min depth"));
        metaData.add(new VCFInfoHeaderLine("MAX_DP", 1, VCFHeaderLineType.Integer, "Max depth"));
        metaData.add(new VCFHeaderLine(Biostar78285.class.getSimpleName() + ".SamFilter", this.filter.toString()));
        for (final Integer treshold : this.minDepthTresholds) {
            metaData.add(new VCFFilterHeaderLine("DP_LT_" + treshold, "All  genotypes have DP< " + treshold));
            metaData.add(new VCFInfoHeaderLine("NUM_DP_LT_" + treshold, 1, VCFHeaderLineType.Integer, "Number of genotypes having DP< " + treshold));
            metaData.add(new VCFInfoHeaderLine("FRACT_DP_LT_" + treshold, 1, VCFHeaderLineType.Float, "Fraction of genotypes having DP< " + treshold));
        }
        if (indexedFastaSequenceFile != null) {
            metaData.add(new VCFInfoHeaderLine("GC_PERCENT", 1, VCFHeaderLineType.Integer, "GC% window_size:" + this.gc_percent_window));
        }
        final List<Allele> refAlleles = Collections.singletonList(Allele.create("N", true));
        final List<Allele> NO_CALLS = Arrays.asList(Allele.NO_CALL, Allele.NO_CALL);
        final VCFHeader vcfHeader = new VCFHeader(metaData, samples);
        vcfHeader.setSequenceDictionary(dict);
        out.writeHeader(vcfHeader);
        final SAMRecordCoordinateComparator samRecordCoordinateComparator = new SAMRecordCoordinateComparator();
        final PeekableIterator<SAMRecord> peekIter = new PeekableIterator<>(new MergingIterator<>((R1, R2) -> samRecordCoordinateComparator.fileOrderCompare(R1, R2), samIterators));
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(dict);
        for (final SAMSequenceRecord ssr : dict.getSequences()) {
            final IntervalTree<Boolean> capturePos;
            if (intervals != null) {
                if (!Arrays.stream(intervals).anyMatch(I -> I.referenceIndex == ssr.getSequenceIndex())) {
                    continue;
                }
                capturePos = new IntervalTree<>();
                Arrays.stream(intervals).filter(I -> I.referenceIndex == ssr.getSequenceIndex()).forEach(I -> capturePos.put(I.start, I.end, true));
                ;
            } else {
                capturePos = null;
            }
            final GenomicSequence genomicSequence;
            if (indexedFastaSequenceFile != null && indexedFastaSequenceFile.getSequenceDictionary().getSequence(ssr.getSequenceName()) != null) {
                genomicSequence = new GenomicSequence(indexedFastaSequenceFile, ssr.getSequenceName());
            } else {
                genomicSequence = null;
            }
            final List<SAMRecord> buffer = new ArrayList<>();
            for (int ssr_pos = 1; ssr_pos <= ssr.getSequenceLength(); ++ssr_pos) {
                if (capturePos != null && !capturePos.overlappers(ssr_pos, ssr_pos).hasNext())
                    continue;
                progress.watch(ssr.getSequenceName(), ssr_pos);
                while (peekIter.hasNext()) {
                    final SAMRecord rec = peekIter.peek();
                    if (rec.getReadUnmappedFlag()) {
                        // consumme
                        peekIter.next();
                        continue;
                    }
                    if (this.filter.filterOut(rec)) {
                        // consumme
                        peekIter.next();
                        continue;
                    }
                    if (rec.getReferenceIndex() < ssr.getSequenceIndex()) {
                        throw new IllegalStateException("should not happen");
                    }
                    if (rec.getReferenceIndex() > ssr.getSequenceIndex()) {
                        break;
                    }
                    if (rec.getAlignmentEnd() < ssr_pos) {
                        throw new IllegalStateException("should not happen");
                    }
                    if (rec.getAlignmentStart() > ssr_pos) {
                        break;
                    }
                    buffer.add(peekIter.next());
                }
                int x = 0;
                while (x < buffer.size()) {
                    final SAMRecord R = buffer.get(x);
                    if (R.getReferenceIndex() != ssr.getSequenceIndex() || R.getAlignmentEnd() < ssr_pos) {
                        buffer.remove(x);
                    } else {
                        x++;
                    }
                }
                final Map<String, PosInfo> count = samples.stream().map(S -> new PosInfo(S)).collect(Collectors.toMap(P -> P.sample, Function.identity()));
                for (final SAMRecord rec : buffer) {
                    if (rec.getReferenceIndex() != ssr.getSequenceIndex())
                        throw new IllegalStateException("should not happen");
                    if (rec.getAlignmentEnd() < ssr_pos)
                        continue;
                    if (rec.getAlignmentStart() > ssr_pos)
                        continue;
                    final Cigar cigar = rec.getCigar();
                    if (cigar == null)
                        continue;
                    int refpos = rec.getAlignmentStart();
                    final String sample = this.partition.getPartion(rec, DEFAULT_PARTITION);
                    for (final CigarElement ce : cigar.getCigarElements()) {
                        if (refpos > ssr_pos)
                            break;
                        final CigarOperator op = ce.getOperator();
                        if (op.consumesReferenceBases()) {
                            if (op.consumesReadBases()) {
                                if (refpos <= ssr_pos && ssr_pos <= refpos + ce.getLength()) {
                                    final PosInfo posInfo = count.get(sample);
                                    if (posInfo != null) {
                                        posInfo.dp++;
                                        if (rec.getReadNegativeStrandFlag()) {
                                            posInfo.negative_strand++;
                                        }
                                    }
                                    break;
                                }
                            }
                            refpos += ce.getLength();
                        }
                    }
                }
                final VariantContextBuilder vcb = new VariantContextBuilder();
                final Set<String> filters = new HashSet<>();
                vcb.chr(ssr.getSequenceName());
                vcb.start(ssr_pos);
                vcb.stop(ssr_pos);
                if (genomicSequence == null) {
                    vcb.alleles(refAlleles);
                } else {
                    vcb.alleles(Collections.singletonList(Allele.create((byte) genomicSequence.charAt(ssr_pos - 1), true)));
                    final GenomicSequence.GCPercent gcp = genomicSequence.getGCPercent(Math.max((ssr_pos - 1) - this.gc_percent_window, 0), Math.min(ssr_pos + this.gc_percent_window, ssr.getSequenceLength()));
                    if (!gcp.isEmpty()) {
                        vcb.attribute("GC_PERCENT", gcp.getGCPercentAsInteger());
                    }
                }
                vcb.attribute(VCFConstants.DEPTH_KEY, (int) count.values().stream().mapToInt(S -> S.dp).sum());
                vcb.genotypes(count.values().stream().map(C -> new GenotypeBuilder(C.sample, NO_CALLS).DP((int) C.dp).attribute("DR", C.negative_strand).attribute("DF", C.dp - C.negative_strand).make()).collect(Collectors.toList()));
                for (final Integer treshold : this.minDepthTresholds) {
                    final int count_lt = (int) count.values().stream().filter(S -> S.dp < treshold).count();
                    if (count_lt == samples.size()) {
                        filters.add("DP_LT_" + treshold);
                    }
                    vcb.attribute("NUM_DP_LT_" + treshold, count_lt);
                    if (!samples.isEmpty()) {
                        vcb.attribute("FRACT_DP_LT_" + treshold, count_lt / (float) samples.size());
                    }
                }
                if (!samples.isEmpty()) {
                    final int[] array = count.values().stream().mapToInt(S -> S.dp).toArray();
                    vcb.attribute("AVG_DP", Percentile.average().evaluate(array));
                    vcb.attribute("MEDIAN_DP", Percentile.median().evaluate(array));
                    vcb.attribute("MIN_DP", (int) Percentile.min().evaluate(array));
                    vcb.attribute("MAX_DP", (int) Percentile.max().evaluate(array));
                }
                if (filters.isEmpty()) {
                    vcb.passFilters();
                } else {
                    vcb.filters(filters);
                }
                out.add(vcb.make());
            }
        }
        progress.finish();
        peekIter.close();
        out.close();
        out = null;
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(out);
        CloserUtil.close(samIterators);
        CloserUtil.close(samReaders);
        CloserUtil.close(indexedFastaSequenceFile);
    }
}

Example 53 with SamReaderFactory

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

the class BamCmpCoverage method doWork.

@Override
public int doWork(final List<String> args) {
    if (outputFile == null) {
        LOG.error("output image file not defined");
        return -1;
    }
    if (this.imgageSize < 1) {
        LOG.error("Bad image size:" + this.imgageSize);
        return -1;
    }
    if (this.minDepth < 0) {
        LOG.error("Bad min depth : " + this.minDepth);
        return -1;
    }
    if (this.minDepth >= this.maxDepth) {
        LOG.error("Bad min<max depth : " + this.minDepth + "<" + this.maxDepth);
        return 1;
    }
    if (this.bedFile != null) {
        readBedFile(this.bedFile);
    }
    if (regionStr != null && this.intervals != null) {
        LOG.error("bed and interval both defined.");
        return -1;
    }
    try {
        final ConcatSam.Factory concatSamFactory = new ConcatSam.Factory();
        final SamReaderFactory srf = concatSamFactory.getSamReaderFactory();
        srf.disable(SamReaderFactory.Option.EAGERLY_DECODE);
        srf.disable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS);
        srf.disable(SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS);
        if (this.regionStr != null) {
            concatSamFactory.addInterval(this.regionStr);
        }
        ConcatSam.ConcatSamIterator concatIter = concatSamFactory.open(args);
        final SAMSequenceDictionary dict = concatIter.getFileHeader().getSequenceDictionary();
        final Set<String> samples = concatIter.getFileHeader().getReadGroups().stream().map(RG -> this.samRecordPartition.apply(RG, "N/A")).collect(Collectors.toSet());
        LOG.info("Samples:" + samples.size());
        for (String sample : samples) {
            this.sample2column.put(sample, this.sample2column.size());
        }
        // create image
        LOG.info("Creating image " + this.imgageSize + "x" + this.imgageSize);
        this.image = new BufferedImage(this.imgageSize, this.imgageSize, BufferedImage.TYPE_INT_RGB);
        Graphics2D g = this.image.createGraphics();
        this.marginWidth = this.imgageSize * 0.05;
        double drawingWidth = (this.imgageSize - 1) - marginWidth;
        this.sampleWidth = drawingWidth / samples.size();
        // g.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, this.imgageSize, this.imgageSize);
        g.setColor(Color.BLACK);
        Hershey hershey = new Hershey();
        for (final String sample_x : samples) {
            double labelHeight = marginWidth;
            if (labelHeight > 50)
                labelHeight = 50;
            g.setColor(Color.BLACK);
            hershey.paint(g, sample_x, marginWidth + sample2column.get(sample_x) * sampleWidth, marginWidth - labelHeight, sampleWidth * 0.9, labelHeight * 0.9);
            AffineTransform old = g.getTransform();
            AffineTransform tr = AffineTransform.getTranslateInstance(marginWidth, marginWidth + sample2column.get(sample_x) * sampleWidth);
            tr.rotate(Math.PI / 2);
            g.setTransform(tr);
            hershey.paint(g, sample_x, 0.0, 0.0, sampleWidth * 0.9, labelHeight * 0.9);
            // g.drawString(this.tabixFile.getFile().getName(),0,0);
            g.setTransform(old);
            for (String sample_y : samples) {
                Rectangle2D rect = new Rectangle2D.Double(marginWidth + sample2column.get(sample_x) * sampleWidth, marginWidth + sample2column.get(sample_y) * sampleWidth, sampleWidth, sampleWidth);
                g.setColor(Color.BLUE);
                g.draw(new Line2D.Double(rect.getMinX(), rect.getMinY(), rect.getMaxX(), rect.getMaxY()));
                g.setColor(Color.BLACK);
                g.draw(rect);
            }
        }
        // ceate bit-array
        BitSampleMatrix bitMatrix = new BitSampleMatrix(samples.size());
        // preivous chrom
        // int prev_tid=-1;
        BufferedList<Depth> depthList = new BufferedList<Depth>();
        g.setColor(Color.BLACK);
        SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(dict).logger(LOG);
        LOG.info("Scanning bams...");
        while (concatIter.hasNext()) {
            final SAMRecord rec = progress.watch(concatIter.next());
            if (this.samRecordFilter.filterOut(rec))
                continue;
            final String sample = this.samRecordPartition.getPartion(rec, "N/A");
            final int sample_id = this.sample2column.get(sample);
            final Cigar cigar = rec.getCigar();
            if (cigar == null || cigar.isEmpty())
                continue;
            int refPos = rec.getAlignmentStart();
            /* cleanup front pos */
            while (!depthList.isEmpty()) {
                final Depth front = depthList.getFirst();
                if (front.tid != rec.getReferenceIndex().intValue() || front.pos < refPos) {
                    paint(bitMatrix, front);
                    depthList.removeFirst();
                    continue;
                } else {
                    break;
                }
            }
            for (final CigarElement ce : cigar.getCigarElements()) {
                final CigarOperator op = ce.getOperator();
                if (!op.consumesReferenceBases())
                    continue;
                if (op.consumesReadBases()) {
                    for (int i = 0; i < ce.getLength(); ++i) {
                        Depth depth = null;
                        int pos = refPos + i;
                        // ignore non-overlapping BED
                        if (this.intervals != null && !this.intervals.containsOverlapping(new Interval(rec.getReferenceName(), pos, pos))) {
                            continue;
                        } else if (depthList.isEmpty()) {
                            depth = new Depth();
                            depth.pos = pos;
                            depth.tid = rec.getReferenceIndex();
                            depthList.add(depth);
                        } else if (depthList.getLast().pos < pos) {
                            Depth prev = depthList.getLast();
                            while (prev.pos < pos) {
                                depth = new Depth();
                                depth.pos = prev.pos + 1;
                                depth.tid = rec.getReferenceIndex();
                                depthList.add(depth);
                                prev = depth;
                            }
                            depth = prev;
                        } else {
                            int lastPos = depthList.get(depthList.size() - 1).pos;
                            int distance = lastPos - pos;
                            int indexInList = (depthList.size() - 1) - (distance);
                            if (indexInList < 0) {
                                // can appen when BED declared and partially overlap the read
                                continue;
                            }
                            depth = depthList.get((depthList.size() - 1) - (distance));
                            if (depth.pos != pos) {
                                LOG.error(" " + pos + " vs " + depth.pos + " " + lastPos);
                                return -1;
                            }
                        }
                        depth.depths[sample_id]++;
                    }
                }
                refPos += ce.getLength();
            }
        }
        while (!depthList.isEmpty()) {
            // paint(g,depthList.remove(0));
            paint(bitMatrix, depthList.remove(0));
        }
        progress.finish();
        concatIter.close();
        concatIter = null;
        for (int x = 0; x < bitMatrix.n_samples; ++x) {
            for (int y = 0; y < bitMatrix.n_samples; ++y) {
                LOG.info("Painting...(" + x + "/" + y + ")");
                paint(g, bitMatrix.get(x, y));
            }
        }
        g.dispose();
        // save file
        LOG.info("saving " + this.outputFile);
        if (this.outputFile.getName().toLowerCase().endsWith(".png")) {
            ImageIO.write(this.image, "PNG", this.outputFile);
        } else {
            ImageIO.write(this.image, "JPG", this.outputFile);
        }
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
    }
}

Example 54 with SamReaderFactory

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

the class Bam2Raster method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.regionStr == null) {
        LOG.error("Region was not defined.");
        return -1;
    }
    if (this.WIDTH < 100) {
        LOG.info("adjusting WIDTH to 100");
        this.WIDTH = 100;
    }
    SamReader samFileReader = null;
    try {
        final SamReaderFactory srf = super.createSamReaderFactory();
        if (this.referenceFile != null) {
            LOG.info("loading reference");
            this.indexedFastaSequenceFile = new IndexedFastaSequenceFile(this.referenceFile);
            srf.referenceSequence(this.referenceFile);
        }
        final IntervalParser intervalParser = new IntervalParser(this.indexedFastaSequenceFile == null ? null : this.indexedFastaSequenceFile.getSequenceDictionary()).setFixContigName(true);
        this.interval = intervalParser.parse(this.regionStr);
        if (this.interval == null) {
            LOG.error("Cannot parse interval " + regionStr + " or chrom doesn't exists in sam dictionary.");
            return -1;
        }
        LOG.info("Interval is " + this.interval);
        loadVCFs();
        for (final String bamFile : IOUtils.unrollFiles(args)) {
            samFileReader = srf.open(SamInputResource.of(bamFile));
            final SAMFileHeader header = samFileReader.getFileHeader();
            final SAMSequenceDictionary dict = header.getSequenceDictionary();
            if (dict == null) {
                LOG.error("no dict in " + bamFile);
                return -1;
            }
            if (dict.getSequence(this.interval.getContig()) == null) {
                LOG.error("no such chromosome in " + bamFile + " " + this.interval);
                return -1;
            }
            scan(samFileReader);
            samFileReader.close();
            samFileReader = null;
        }
        if (this.key2partition.isEmpty()) {
            LOG.error("No data was found.(not Read-Group specified ?");
            return -1;
        }
        this.key2partition.values().stream().forEach(P -> P.build());
        saveImages(this.key2partition.values().stream().map(P -> P.image).collect(Collectors.toList()));
        return RETURN_OK;
    } catch (Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(indexedFastaSequenceFile);
        CloserUtil.close(samFileReader);
        indexedFastaSequenceFile = null;
    }
}

Example 55 with SamReaderFactory

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

the class BamStats05 method doWork.

@Override
public int doWork(List<String> args) {
    if (BEDILE == null) {
        LOG.error("missing bed file");
        return -1;
    }
    SamReader in = null;
    BufferedReader r = null;
    PrintWriter pw = null;
    try {
        Map<String, List<Interval>> gene2interval = readBedFile(BEDILE);
        pw = super.openFileOrStdoutAsPrintWriter(this.outputFile);
        pw.println("#chrom\tstart\tend\tgene\tsample\tlength\tmincov\tmaxcov\tavg\tnocoverage.bp\tpercentcovered");
        final SamReaderFactory srf = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.SILENT);
        final Set<String> files = new HashSet<>();
        if (args.isEmpty()) {
            LOG.info("reading BAM path from stdin");
            r = new BufferedReader(new InputStreamReader(stdin()));
            String line;
            while ((line = r.readLine()) != null) {
                if (line.startsWith("#") || line.trim().isEmpty())
                    continue;
                if (!line.endsWith(".bam")) {
                    LOG.error("line should end with .bam :" + line);
                    return -1;
                }
            }
            CloserUtil.close(r);
        } else {
            for (final String fname : args) {
                files.addAll(IOUtils.unrollFiles(Collections.singletonList(fname)));
            }
        }
        for (final String f : files) {
            in = srf.open(new File(f));
            int tl = doWork(pw, gene2interval, f, in);
            CloserUtil.close(in);
            in = null;
            if (tl != 0)
                return tl;
        }
        pw.flush();
        pw.close();
        pw = null;
        return RETURN_OK;
    } catch (Exception e) {
        LOG.error(e);
        return -1;
    } finally {
        CloserUtil.close(in);
        CloserUtil.close(r);
        CloserUtil.close(pw);
    }
}