Search in sources :

Example 11 with BedLineCodec

use of com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec in project jvarkit by lindenb.

the class BedLiftOver method scan.

private void scan(BufferedReader r, PrintWriter out, PrintWriter failed) throws IOException {
    String line;
    final BedLineCodec bedCodec = new BedLineCodec();
    while ((line = r.readLine()) != null) {
        if (line.startsWith("#") || line.trim().isEmpty())
            continue;
        final BedLine bedLine = bedCodec.decode(line);
        if (bedLine == null)
            continue;
        final Interval srcInterval = bedLine.toInterval();
        Interval dest = this.liftOver.liftOver(srcInterval);
        if (dest != null) {
            out.print(dest.getContig());
            out.print('\t');
            out.print(dest.getStart() - 1);
            out.print('\t');
            out.print(dest.getEnd());
            for (int i = 3; i < bedLine.getColumnCount(); ++i) {
                out.print('\t');
                out.print(bedLine.get(i));
            }
            out.println();
        } else if (failed != null) {
            failed.println(line);
        }
    }
}
Also used : BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) BedLine(com.github.lindenb.jvarkit.util.bio.bed.BedLine) Interval(htsjdk.samtools.util.Interval)

Example 12 with BedLineCodec

use of com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec in project jvarkit by lindenb.

the class LumpySort method doWork.

@Override
public int doWork(final List<String> args) {
    VariantContextWriter vcw = null;
    LineIterator vcfIn = null;
    Environment environment = null;
    Database variantsDb1 = null;
    final List<File> inputs = IOUtil.unrollFiles(args.stream().map(S -> new File(S)).collect(Collectors.toList()), ".vcf", ".vcf.gz");
    if (inputs.isEmpty()) {
        LOG.error("empty vcf list");
        return -1;
    }
    try {
        IOUtil.assertDirectoryIsWritable(this.bdbHomeDir);
        final Set<VCFHeaderLine> metaData = new HashSet<>();
        final Set<String> sampleNames = new TreeSet<>();
        final IntervalTreeMap<Boolean> intervalTreeMapBed;
        if (this.bedFile != null) {
            intervalTreeMapBed = new IntervalTreeMap<>();
            final BedLineCodec bedLineCodec = new BedLineCodec();
            final BufferedReader br = IOUtils.openFileForBufferedReading(this.bedFile);
            br.lines().map(L -> bedLineCodec.decode(L)).filter(L -> L != null).forEach(B -> intervalTreeMapBed.put(B.toInterval(), true));
            br.close();
        } else {
            intervalTreeMapBed = null;
        }
        for (int idx = 0; idx < inputs.size(); ++idx) {
            final File vcfFile = inputs.get(idx);
            LOG.info("Read header " + (idx + 1) + "/" + inputs.size());
            final VCFFileReader r = new VCFFileReader(vcfFile, false);
            final VCFHeader header = r.getFileHeader();
            if (!LumpyConstants.isLumpyHeader(header)) {
                LOG.error("doesn't look like a Lumpy-SV vcf header " + vcfFile);
                r.close();
                return -1;
            }
            if (!header.hasGenotypingData()) {
                LOG.error("No sample in " + vcfFile);
                r.close();
                return -1;
            }
            for (final String sampleName : header.getSampleNamesInOrder()) {
                if (sampleNames.contains(sampleName)) {
                    LOG.error("Sample found twice " + sampleName + " in " + vcfFile);
                    r.close();
                    return -1;
                }
                sampleNames.add(sampleName);
            }
            metaData.addAll(header.getMetaDataInInputOrder().stream().filter(H -> !H.getKey().equals("fileDate")).collect(Collectors.toSet()));
            r.close();
        }
        final VCFInfoHeaderLine nSampleInfoHeaderLine = new VCFInfoHeaderLine("NSAMPLES", 1, VCFHeaderLineType.Integer, "Number of affected samples.");
        metaData.add(nSampleInfoHeaderLine);
        final VCFFormatHeaderLine chromStartFormatHeaderLine = new VCFFormatHeaderLine("CB", 1, VCFHeaderLineType.Integer, "Original Variant POS");
        metaData.add(chromStartFormatHeaderLine);
        final VCFFormatHeaderLine chromEndFormatHeaderLine = new VCFFormatHeaderLine("CE", 1, VCFHeaderLineType.Integer, "Original Variant END");
        metaData.add(chromEndFormatHeaderLine);
        final VCFHeader outHeader = new VCFHeader(metaData, sampleNames);
        final VCFHeaderVersion[] versions = VCFHeaderVersion.values();
        this.vcfEncoder = new VCFEncoder(outHeader, false, true);
        this.vcfCodec.setVCFHeader(outHeader, versions[versions.length - 1]);
        /* open BDB env */
        final Transaction txn = null;
        environment = new Environment(this.bdbHomeDir, new EnvironmentConfig().setAllowCreate(true).setReadOnly(false));
        variantsDb1 = environment.openDatabase(txn, "variants1", new DatabaseConfig().setBtreeComparator(KeySorterComparator.class).setAllowCreate(true).setReadOnly(false).setTemporary(true));
        long total_variants = 0L;
        final LumpyVarBinding lumpVarBinding = new LumpyVarBinding();
        final KeySorterBinding keySorterBinding = new KeySorterBinding();
        for (int idx = 0; idx < inputs.size(); ++idx) {
            final long millisecstart = System.currentTimeMillis();
            final File vcfFile = inputs.get(idx);
            int nVariant = 0;
            final VCFFileReader r = new VCFFileReader(vcfFile, false);
            final List<Genotype> missing = new ArrayList<>(sampleNames.size());
            for (final String sn : sampleNames) {
                if (r.getFileHeader().getSampleNamesInOrder().contains(sn))
                    continue;
                missing.add(GenotypeBuilder.createMissing(sn, 2));
            }
            final CloseableIterator<VariantContext> iter = r.iterator();
            while (iter.hasNext()) {
                VariantContext ctx = iter.next();
                if (!this.keep_secondary) {
                    if (ctx.hasAttribute("SECONDARY"))
                        continue;
                }
                if (!this.variantFilter.test(ctx))
                    continue;
                if (intervalTreeMapBed != null && !intervalTreeMapBed.containsOverlapping(ctx))
                    continue;
                final List<Genotype> gtList = new ArrayList<>(ctx.getGenotypes());
                for (int gi = 0; gi < gtList.size(); gi++) {
                    Genotype g = gtList.get(gi);
                    final GenotypeBuilder gb;
                    if (this.do_genotype && isAvailableGenotype(g)) {
                        gb = new GenotypeBuilder(g.getSampleName(), ctx.getAlternateAlleles());
                        gb.attributes(g.getExtendedAttributes());
                    } else {
                        gb = new GenotypeBuilder(g);
                    }
                    gb.attribute(chromStartFormatHeaderLine.getID(), ctx.getStart());
                    gb.attribute(chromEndFormatHeaderLine.getID(), ctx.getEnd());
                    gtList.set(gi, gb.make());
                }
                gtList.addAll(missing);
                ctx = new VariantContextBuilder(ctx).genotypes(gtList).rmAttribute("PRPOS").make();
                final LumpyVar lvar = new LumpyVar(ctx, total_variants);
                final DatabaseEntry key = new DatabaseEntry();
                final DatabaseEntry data = new DatabaseEntry();
                lumpVarBinding.objectToEntry(lvar, data);
                keySorterBinding.objectToEntry(lvar.getSortKey(), key);
                if (variantsDb1.put(txn, key, data) != OperationStatus.SUCCESS) {
                    r.close();
                    LOG.error("insertion failed");
                    return -1;
                }
                nVariant++;
                total_variants++;
            }
            iter.close();
            r.close();
            LOG.info("Read  " + (idx + 1) + "/" + inputs.size() + " variants of " + vcfFile + " N=" + nVariant + " Total:" + total_variants + " That took: " + Duration.ofMillis(System.currentTimeMillis() - millisecstart));
            System.gc();
        }
        if (intervalTreeMapBed != null)
            intervalTreeMapBed.clear();
        System.gc();
        LOG.info("Writing output");
        final List<Allele> ALLELES_NO_CALLS = this.do_genotype ? Collections.singletonList(Allele.NO_CALL) : Arrays.asList(Allele.NO_CALL, Allele.NO_CALL);
        final Cursor cursor = variantsDb1.openCursor(txn, null);
        vcw = super.openVariantContextWriter(this.outputFile);
        vcw.writeHeader(outHeader);
        for (; ; ) {
            final DatabaseEntry key = new DatabaseEntry();
            final DatabaseEntry data = new DatabaseEntry();
            OperationStatus status = cursor.getNext(key, data, LockMode.DEFAULT);
            if (!status.equals(OperationStatus.SUCCESS))
                break;
            final LumpyVar first = lumpVarBinding.entryToObject(data);
            if (this.do_not_merge_ctx) {
                vcw.add(first.ctx);
                continue;
            }
            final KeySorter keySorter1 = keySorterBinding.entryToObject(key);
            final List<LumpyVar> buffer = new ArrayList<>();
            buffer.add(first);
            final DatabaseEntry key2 = new DatabaseEntry();
            final DatabaseEntry data2 = new DatabaseEntry();
            final Cursor cursor2 = cursor.dup(true);
            for (; ; ) {
                status = cursor2.getNext(key2, data2, LockMode.DEFAULT);
                if (!status.equals(OperationStatus.SUCCESS))
                    break;
                final KeySorter keySorter2 = keySorterBinding.entryToObject(key2);
                if (keySorter1.compare1(keySorter2) != 0) {
                    break;
                }
                final LumpyVar lv = lumpVarBinding.entryToObject(data2);
                if (lv.ctx.getStart() > first.ctx.getEnd()) {
                    break;
                }
                if (first.canMerge(lv)) {
                    buffer.add(lv);
                    cursor2.delete();
                }
            }
            cursor2.close();
            // delete 'first'
            cursor.delete();
            final int variantStartA = buffer.stream().mapToInt(V -> V.ctx.getStart()).min().getAsInt();
            final int variantStartB = (int) buffer.stream().mapToInt(V -> V.ctx.getStart()).average().getAsDouble();
            final int variantStartC = buffer.stream().mapToInt(V -> V.ctx.getStart()).max().getAsInt();
            final int variantEndA = buffer.stream().mapToInt(V -> V.ctx.getEnd()).min().getAsInt();
            final int variantEndB = (int) buffer.stream().mapToInt(V -> V.ctx.getEnd()).average().getAsDouble();
            final int variantEndC = buffer.stream().mapToInt(V -> V.ctx.getEnd()).max().getAsInt();
            final VariantContextBuilder vcb = new VariantContextBuilder("lumpymerge", first.ctx.getContig(), variantStartB, variantEndB, first.ctx.getAlleles());
            vcb.attribute("END", variantEndB);
            vcb.attribute("SVTYPE", first.ctx.getAttribute("SVTYPE"));
            vcb.attribute("SVLEN", (int) Percentile.median().evaluate(buffer.stream().mapToInt(V -> V.ctx.getEnd() - V.ctx.getStart())));
            vcb.attribute("CIPOS", Arrays.asList(variantStartB - variantStartA, variantStartC - variantStartB));
            vcb.attribute("CIEND", Arrays.asList(variantEndB - variantEndA, variantEndC - variantEndB));
            vcb.attribute("SU", buffer.stream().flatMap(V -> V.ctx.getGenotypes().stream()).mapToInt(G -> G.getAttributeAsInt("SU", 0)).sum());
            vcb.attribute("SR", buffer.stream().flatMap(V -> V.ctx.getGenotypes().stream()).mapToInt(G -> G.getAttributeAsInt("SR", 0)).sum());
            vcb.attribute("PE", buffer.stream().flatMap(V -> V.ctx.getGenotypes().stream()).mapToInt(G -> G.getAttributeAsInt("PE", 0)).sum());
            final Map<String, Genotype> sample2genotype = new HashMap<>(sampleNames.size());
            buffer.stream().flatMap(V -> V.ctx.getGenotypes().stream()).filter(G -> isAvailableGenotype(G)).forEach(G -> {
                sample2genotype.put(G.getSampleName(), G);
            });
            vcb.attribute(nSampleInfoHeaderLine.getID(), sample2genotype.size());
            for (final String sn : sampleNames) {
                if (!sample2genotype.containsKey(sn)) {
                    sample2genotype.put(sn, new GenotypeBuilder(sn, ALLELES_NO_CALLS).attribute("SU", 0).attribute("SR", 0).attribute("PE", 0).make());
                }
            }
            vcb.genotypes(sample2genotype.values());
            vcw.add(vcb.make());
        }
        cursor.close();
        vcw.close();
        vcw = null;
        variantsDb1.close();
        variantsDb1 = null;
        environment.close();
        environment = null;
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(vcfIn);
        CloserUtil.close(vcw);
        CloserUtil.close(variantsDb1);
        CloserUtil.close(environment);
    }
}
Also used : Allele(htsjdk.variant.variantcontext.Allele) Arrays(java.util.Arrays) JexlVariantPredicate(com.github.lindenb.jvarkit.util.vcf.JexlVariantPredicate) Program(com.github.lindenb.jvarkit.util.jcommander.Program) LineIterator(htsjdk.tribble.readers.LineIterator) IOUtil(htsjdk.samtools.util.IOUtil) VCFFileReader(htsjdk.variant.vcf.VCFFileReader) VCFHeader(htsjdk.variant.vcf.VCFHeader) LockMode(com.sleepycat.je.LockMode) VCFEncoder(htsjdk.variant.vcf.VCFEncoder) VCFHeaderVersion(htsjdk.variant.vcf.VCFHeaderVersion) Duration(java.time.Duration) Map(java.util.Map) CloserUtil(htsjdk.samtools.util.CloserUtil) TupleInput(com.sleepycat.bind.tuple.TupleInput) IntervalTreeMap(htsjdk.samtools.util.IntervalTreeMap) GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder) Predicate(java.util.function.Predicate) Logger(com.github.lindenb.jvarkit.util.log.Logger) DatabaseEntry(com.sleepycat.je.DatabaseEntry) Set(java.util.Set) Collectors(java.util.stream.Collectors) Percentile(com.github.lindenb.jvarkit.math.stats.Percentile) List(java.util.List) DatabaseConfig(com.sleepycat.je.DatabaseConfig) StructuralVariantType(htsjdk.variant.variantcontext.StructuralVariantType) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine) VariantContext(htsjdk.variant.variantcontext.VariantContext) EnvironmentConfig(com.sleepycat.je.EnvironmentConfig) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) Genotype(htsjdk.variant.variantcontext.Genotype) VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) OperationStatus(com.sleepycat.je.OperationStatus) CloseableIterator(htsjdk.samtools.util.CloseableIterator) Parameter(com.beust.jcommander.Parameter) TupleOutput(com.sleepycat.bind.tuple.TupleOutput) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) HashMap(java.util.HashMap) Function(java.util.function.Function) TreeSet(java.util.TreeSet) TupleBinding(com.sleepycat.bind.tuple.TupleBinding) ArrayList(java.util.ArrayList) HashSet(java.util.HashSet) Interval(htsjdk.samtools.util.Interval) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) Database(com.sleepycat.je.Database) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) VCFCodec(htsjdk.variant.vcf.VCFCodec) VCFHeaderLineType(htsjdk.variant.vcf.VCFHeaderLineType) Environment(com.sleepycat.je.Environment) File(java.io.File) Cursor(com.sleepycat.je.Cursor) VCFFormatHeaderLine(htsjdk.variant.vcf.VCFFormatHeaderLine) BufferedReader(java.io.BufferedReader) Comparator(java.util.Comparator) Transaction(com.sleepycat.je.Transaction) Collections(java.util.Collections) VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) HashMap(java.util.HashMap) ArrayList(java.util.ArrayList) Cursor(com.sleepycat.je.Cursor) TreeSet(java.util.TreeSet) Database(com.sleepycat.je.Database) HashSet(java.util.HashSet) VCFHeaderVersion(htsjdk.variant.vcf.VCFHeaderVersion) EnvironmentConfig(com.sleepycat.je.EnvironmentConfig) Genotype(htsjdk.variant.variantcontext.Genotype) GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) VCFEncoder(htsjdk.variant.vcf.VCFEncoder) Transaction(com.sleepycat.je.Transaction) Environment(com.sleepycat.je.Environment) File(java.io.File) VCFFileReader(htsjdk.variant.vcf.VCFFileReader) VariantContext(htsjdk.variant.variantcontext.VariantContext) DatabaseEntry(com.sleepycat.je.DatabaseEntry) LineIterator(htsjdk.tribble.readers.LineIterator) OperationStatus(com.sleepycat.je.OperationStatus) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VCFHeader(htsjdk.variant.vcf.VCFHeader) VCFFormatHeaderLine(htsjdk.variant.vcf.VCFFormatHeaderLine) DatabaseConfig(com.sleepycat.je.DatabaseConfig) VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine) Allele(htsjdk.variant.variantcontext.Allele) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) BufferedReader(java.io.BufferedReader)

Example 13 with BedLineCodec

use of com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec in project jvarkit by lindenb.

the class GcPercentAndDepth method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.windowSize <= 0) {
        LOG.error("Bad window size.");
        return -1;
    }
    if (this.windowStep <= 0) {
        LOG.error("Bad window step.");
        return -1;
    }
    if (this.refFile == null) {
        LOG.error("Undefined REF File");
        return -1;
    }
    if (args.isEmpty()) {
        LOG.error("Illegal Number of arguments.");
        return -1;
    }
    ReferenceGenome indexedFastaSequenceFile = null;
    List<SamReader> readers = new ArrayList<SamReader>();
    PrintWriter out = null;
    try {
        LOG.info("Loading " + this.refFile);
        indexedFastaSequenceFile = new ReferenceGenomeFactory().openFastaFile(this.refFile);
        this.samSequenceDictionary = indexedFastaSequenceFile.getDictionary();
        if (this.samSequenceDictionary == null) {
            LOG.error("Cannot get sequence dictionary for " + this.refFile);
            return -1;
        }
        out = super.openFileOrStdoutAsPrintWriter(outPutFile);
        Set<String> all_samples = new TreeSet<String>();
        /* create input, collect sample names */
        for (int optind = 0; optind < args.size(); ++optind) {
            LOG.info("Opening " + args.get(optind));
            final SamReader samFileReaderScan = super.openSamReader(args.get(optind));
            readers.add(samFileReaderScan);
            final SAMFileHeader header = samFileReaderScan.getFileHeader();
            if (!SequenceUtil.areSequenceDictionariesEqual(this.samSequenceDictionary, header.getSequenceDictionary())) {
                LOG.error(JvarkitException.DictionariesAreNotTheSame.getMessage(this.samSequenceDictionary, header.getSequenceDictionary()));
                return -1;
            }
            for (final SAMReadGroupRecord g : header.getReadGroups()) {
                final String sample = this.partition.apply(g);
                if (StringUtil.isBlank(sample)) {
                    LOG.warning("Read group " + g.getId() + " has no sample in merged dictionary");
                    continue;
                }
                all_samples.add(sample);
            }
        }
        LOG.info("N " + this.partition.name() + "=" + all_samples.size());
        /* print header */
        out.print("#");
        if (!this.hide_genomic_index) {
            out.print("id");
            out.print("\t");
        }
        out.print("chrom");
        out.print("\t");
        out.print("start");
        out.print("\t");
        out.print("end");
        out.print("\t");
        out.print("GCPercent");
        for (final String sample : all_samples) {
            out.print("\t");
            out.print(sample);
        }
        out.println();
        final List<RegionCaptured> regionsCaptured = new ArrayList<RegionCaptured>();
        if (bedFile != null) {
            LOG.info("Reading BED:" + bedFile);
            final BedLineCodec bedLineCodec = new BedLineCodec();
            BufferedReader r = IOUtils.openFileForBufferedReading(bedFile);
            r.lines().filter(L -> !L.startsWith("#")).filter(L -> !StringUtil.isBlank(L)).map(L -> bedLineCodec.decode(L)).filter(B -> B != null).forEach(B -> {
                final SAMSequenceRecord ssr = this.samSequenceDictionary.getSequence(B.getContig());
                if (ssr == null) {
                    LOG.warning("Cannot resolve " + B.getContig());
                    return;
                }
                final RegionCaptured roi = new RegionCaptured(ssr, B.getStart() - 1, B.getEnd());
                regionsCaptured.add(roi);
            });
            CloserUtil.close(r);
            LOG.info("end Reading BED:" + bedFile);
            Collections.sort(regionsCaptured);
        } else {
            LOG.info("No capture, peeking everything");
            for (final SAMSequenceRecord ssr : this.samSequenceDictionary.getSequences()) {
                final RegionCaptured roi = new RegionCaptured(ssr, 0, ssr.getSequenceLength());
                regionsCaptured.add(roi);
            }
        }
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(this.samSequenceDictionary).logger(LOG);
        ReferenceContig genomicSequence = null;
        for (final RegionCaptured roi : regionsCaptured) {
            if (genomicSequence == null || !genomicSequence.hasName(roi.getContig())) {
                genomicSequence = indexedFastaSequenceFile.getContig(roi.getContig());
                if (genomicSequence == null) {
                    LOG.error(JvarkitException.ContigNotFoundInDictionary.getMessage(roi.getContig(), this.samSequenceDictionary));
                    return -1;
                }
            }
            Map<String, int[]> sample2depth = new HashMap<String, int[]>();
            Map<String, Double> sample2meanDepth = new HashMap<String, Double>();
            for (final String sample : all_samples) {
                int[] depth = new int[roi.length()];
                Arrays.fill(depth, 0);
                sample2depth.put(sample, depth);
            }
            List<CloseableIterator<SAMRecord>> iterators = new ArrayList<CloseableIterator<SAMRecord>>();
            for (final SamReader r : readers) {
                iterators.add(r.query(roi.getContig(), roi.getStart(), roi.getEnd(), false));
            }
            final MergingIterator<SAMRecord> merginIter = new MergingIterator<>(new SAMRecordCoordinateComparator(), iterators);
            while (merginIter.hasNext()) {
                final SAMRecord rec = merginIter.next();
                if (rec.getReadUnmappedFlag())
                    continue;
                if (this.filter.filterOut(rec))
                    continue;
                final String sample = this.partition.getPartion(rec, null);
                if (sample == null)
                    continue;
                final int[] depth = sample2depth.get(sample);
                if (depth == null)
                    continue;
                final Cigar cigar = rec.getCigar();
                if (cigar == null)
                    continue;
                int refpos1 = rec.getAlignmentStart();
                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) {
                            if (refpos1 + i < roi.getStart())
                                continue;
                            if (refpos1 + i > roi.getEnd())
                                break;
                            depth[refpos1 + i - roi.getStart()]++;
                        }
                    }
                    refpos1 += ce.getLength();
                }
            }
            merginIter.close();
            for (final RegionCaptured.SlidingWindow win : roi) {
                double total = 0f;
                int countN = 0;
                for (int pos1 = win.getStart(); pos1 <= win.getEnd(); ++pos1) {
                    switch(genomicSequence.charAt(pos1 - 1)) {
                        case 'c':
                        case 'C':
                        case 'g':
                        case 'G':
                        case 's':
                        case 'S':
                            {
                                total++;
                                break;
                            }
                        case 'n':
                        case 'N':
                            countN++;
                            break;
                        default:
                            break;
                    }
                }
                if (skip_if_contains_N && countN > 0)
                    continue;
                double GCPercent = total / (double) win.length();
                int max_depth_for_win = 0;
                sample2meanDepth.clear();
                for (final String sample : all_samples) {
                    int[] depth = sample2depth.get(sample);
                    double sum = 0;
                    for (int pos = win.getStart(); pos < win.getEnd() && (pos - roi.getStart()) < depth.length; ++pos) {
                        sum += depth[pos - roi.getStart()];
                    }
                    double mean = (sum / (double) depth.length);
                    max_depth_for_win = Math.max(max_depth_for_win, (int) mean);
                    sample2meanDepth.put(sample, mean);
                }
                if (max_depth_for_win < this.min_depth)
                    continue;
                if (!this.hide_genomic_index) {
                    out.print(win.getGenomicIndex());
                    out.print("\t");
                }
                out.print(win.getContig());
                out.print("\t");
                out.print(win.getStart() - 1);
                out.print("\t");
                out.print(win.getEnd());
                out.print("\t");
                out.printf("%.2f", GCPercent);
                for (String sample : all_samples) {
                    out.print("\t");
                    out.printf("%.2f", (double) sample2meanDepth.get(sample));
                }
                out.println();
            }
        }
        progress.finish();
        out.flush();
        return 0;
    } catch (Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        for (SamReader r : readers) CloserUtil.close(r);
        CloserUtil.close(indexedFastaSequenceFile);
        CloserUtil.close(out);
    }
}
Also used : Cigar(htsjdk.samtools.Cigar) CloseableIterator(htsjdk.samtools.util.CloseableIterator) Arrays(java.util.Arrays) SequenceUtil(htsjdk.samtools.util.SequenceUtil) MergingIterator(htsjdk.samtools.util.MergingIterator) Program(com.github.lindenb.jvarkit.util.jcommander.Program) Parameter(com.beust.jcommander.Parameter) CigarElement(htsjdk.samtools.CigarElement) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) CigarOperator(htsjdk.samtools.CigarOperator) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) HashMap(java.util.HashMap) SAMRecordPartition(com.github.lindenb.jvarkit.util.samtools.SAMRecordPartition) SAMFileHeader(htsjdk.samtools.SAMFileHeader) TreeSet(java.util.TreeSet) ArrayList(java.util.ArrayList) StringUtil(htsjdk.samtools.util.StringUtil) ReferenceGenomeFactory(com.github.lindenb.jvarkit.util.bio.fasta.ReferenceGenomeFactory) Map(java.util.Map) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) CloserUtil(htsjdk.samtools.util.CloserUtil) PrintWriter(java.io.PrintWriter) AbstractIterator(htsjdk.samtools.util.AbstractIterator) Locatable(htsjdk.samtools.util.Locatable) Iterator(java.util.Iterator) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) Logger(com.github.lindenb.jvarkit.util.log.Logger) Set(java.util.Set) SamReader(htsjdk.samtools.SamReader) ReferenceContig(com.github.lindenb.jvarkit.util.bio.fasta.ReferenceContig) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) File(java.io.File) SAMRecord(htsjdk.samtools.SAMRecord) SamRecordFilter(htsjdk.samtools.filter.SamRecordFilter) List(java.util.List) SamRecordJEXLFilter(com.github.lindenb.jvarkit.util.samtools.SamRecordJEXLFilter) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) SAMRecordCoordinateComparator(htsjdk.samtools.SAMRecordCoordinateComparator) BufferedReader(java.io.BufferedReader) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) ReferenceGenome(com.github.lindenb.jvarkit.util.bio.fasta.ReferenceGenome) Collections(java.util.Collections) ReferenceContig(com.github.lindenb.jvarkit.util.bio.fasta.ReferenceContig) HashMap(java.util.HashMap) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) ArrayList(java.util.ArrayList) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) SamReader(htsjdk.samtools.SamReader) SAMRecordCoordinateComparator(htsjdk.samtools.SAMRecordCoordinateComparator) TreeSet(java.util.TreeSet) PrintWriter(java.io.PrintWriter) CloseableIterator(htsjdk.samtools.util.CloseableIterator) ReferenceGenome(com.github.lindenb.jvarkit.util.bio.fasta.ReferenceGenome) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) ReferenceGenomeFactory(com.github.lindenb.jvarkit.util.bio.fasta.ReferenceGenomeFactory) CigarOperator(htsjdk.samtools.CigarOperator) CigarElement(htsjdk.samtools.CigarElement) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) MergingIterator(htsjdk.samtools.util.MergingIterator) Cigar(htsjdk.samtools.Cigar) SAMRecord(htsjdk.samtools.SAMRecord) BufferedReader(java.io.BufferedReader) SAMFileHeader(htsjdk.samtools.SAMFileHeader)

Example 14 with BedLineCodec

use of com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec in project jvarkit by lindenb.

the class VCFAnnoBam method doVcfToVcf.

@Override
protected int doVcfToVcf(final String inputName, final VcfIterator r, final VariantContextWriter w) {
    BufferedReader bedIn = null;
    List<SamReader> samReaders = new ArrayList<SamReader>();
    IntervalTreeMap<Rgn> capture = new IntervalTreeMap<Rgn>();
    try {
        SAMFileHeader firstHeader = null;
        for (final File samFile : new HashSet<File>(BAMFILE)) {
            LOG.info("open bam " + samFile);
            final SamReader samReader = super.openSamReader(samFile.getPath());
            final SAMFileHeader samHeader = samReader.getFileHeader();
            samReaders.add(samReader);
            if (firstHeader == null) {
                firstHeader = samHeader;
            } else if (!SequenceUtil.areSequenceDictionariesEqual(firstHeader.getSequenceDictionary(), samHeader.getSequenceDictionary())) {
                throw new JvarkitException.DictionariesAreNotTheSame(firstHeader.getSequenceDictionary(), samHeader.getSequenceDictionary());
            }
        }
        IntervalList intervalList = new IntervalList(firstHeader);
        LOG.info("read bed " + BEDILE);
        bedIn = IOUtils.openFileForBufferedReading(BEDILE);
        String line;
        final BedLineCodec bedCodec = new BedLineCodec();
        while ((line = bedIn.readLine()) != null) {
            if (line.isEmpty() || line.startsWith("#"))
                continue;
            final BedLine bed = bedCodec.decode(line);
            if (bed == null)
                continue;
            if (firstHeader.getSequenceDictionary().getSequence(bed.getContig()) == null) {
                LOG.error("error in BED +" + BEDILE + " : " + line + " chromosome is not in sequence dict of " + BAMFILE);
                continue;
            }
            intervalList.add(bed.toInterval());
        }
        bedIn.close();
        bedIn = null;
        intervalList = intervalList.sorted();
        for (final Interval interval : intervalList.uniqued()) {
            final Rgn rgn = new Rgn();
            rgn.interval = interval;
            capture.put(rgn.interval, rgn);
        }
        intervalList = null;
        VCFHeader header = r.getHeader();
        VCFHeader h2 = new VCFHeader(header.getMetaDataInInputOrder(), header.getSampleNamesInOrder());
        h2.addMetaDataLine(new VCFInfoHeaderLine(this.capture_tag, 1, VCFHeaderLineType.String, "Capture stats: Format is (start|end|mean|min|max|length|not_covered|percent_covered) BAM files: " + BAMFILE + " CAPTURE:" + BEDILE));
        w.writeHeader(h2);
        while (r.hasNext()) {
            final VariantContext ctx = r.next();
            Interval interval = new Interval(ctx.getContig(), ctx.getStart(), ctx.getEnd());
            Collection<Rgn> rgns = capture.getOverlapping(interval);
            Iterator<Rgn> it = rgns.iterator();
            if (!it.hasNext()) {
                w.add(ctx);
                continue;
            }
            final Rgn rgn = it.next();
            if (!rgn.processed) {
                // LOG.info("processing "+rgn.interval);
                process(rgn, samReaders);
            }
            final VariantContextBuilder b = new VariantContextBuilder(ctx);
            b.attribute(this.capture_tag, rgn.toString());
            w.add(b.make());
        }
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        for (final SamReader samReader : samReaders) CloserUtil.close(samReader);
    }
}
Also used : ArrayList(java.util.ArrayList) VariantContext(htsjdk.variant.variantcontext.VariantContext) VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) BedLine(com.github.lindenb.jvarkit.util.bio.bed.BedLine) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) SamReader(htsjdk.samtools.SamReader) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) IntervalList(htsjdk.samtools.util.IntervalList) BufferedReader(java.io.BufferedReader) SAMFileHeader(htsjdk.samtools.SAMFileHeader) VCFHeader(htsjdk.variant.vcf.VCFHeader) IntervalTreeMap(htsjdk.samtools.util.IntervalTreeMap) File(java.io.File) HashSet(java.util.HashSet) Interval(htsjdk.samtools.util.Interval)

Example 15 with BedLineCodec

use of com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec in project jvarkit by lindenb.

the class IndexCovJfx method filterBed.

private void filterBed(boolean filterOut) {
    final FileChooser fc = new FileChooser();
    fc.setSelectedExtensionFilter(new FileChooser.ExtensionFilter("Bed file", "bed", "bed.gz"));
    final File bed = fc.showOpenDialog(null);
    if (bed == null)
        return;
    BufferedReader r = null;
    final BedLineCodec bedLineCodec = new BedLineCodec();
    try {
        r = IOUtils.openFileForBufferedReading(bed);
        final List<Interval> intervals = r.lines().filter(L -> !StringUtil.isBlank(L)).map(L -> bedLineCodec.decode(L)).filter(B -> B != null).map(B -> B.toInterval()).collect(Collectors.toList());
        r.close();
        r = null;
        this.visibleIndexCovRows.removeIf(R -> intervals.stream().anyMatch(B -> R.overlaps(B.getContig(), B.getStart(), B.getEnd()) == filterOut));
        repaintCanvas();
    } catch (final Exception err) {
        LOG.error(err);
    } finally {
        CloserUtil.close(r);
    }
}
Also used : Button(javafx.scene.control.Button) Arrays(java.util.Arrays) Program(com.github.lindenb.jvarkit.util.jcommander.Program) AbstractList(java.util.AbstractList) VBox(javafx.scene.layout.VBox) MultipleSelectionModel(javafx.scene.control.MultipleSelectionModel) KeyCombination(javafx.scene.input.KeyCombination) Application(javafx.application.Application) ScrollPane(javafx.scene.control.ScrollPane) StringUtil(htsjdk.samtools.util.StringUtil) ListChangeListener(javafx.collections.ListChangeListener) AlertType(javafx.scene.control.Alert.AlertType) Locale(java.util.Locale) Map(java.util.Map) Hershey(com.github.lindenb.jvarkit.util.Hershey) CloserUtil(htsjdk.samtools.util.CloserUtil) Alert(javafx.scene.control.Alert) HBox(javafx.scene.layout.HBox) Rectangle2D(javafx.geometry.Rectangle2D) TextField(javafx.scene.control.TextField) MenuItem(javafx.scene.control.MenuItem) GraphicsContext(javafx.scene.canvas.GraphicsContext) JfxLauncher(com.github.lindenb.jvarkit.util.jcommander.JfxLauncher) Logger(com.github.lindenb.jvarkit.util.log.Logger) Set(java.util.Set) Canvas(javafx.scene.canvas.Canvas) Spinner(javafx.scene.control.Spinner) Screen(javafx.stage.Screen) Collectors(java.util.stream.Collectors) Platform(javafx.application.Platform) SeparatorMenuItem(javafx.scene.control.SeparatorMenuItem) List(java.util.List) SmartComparator(com.github.lindenb.jvarkit.lang.SmartComparator) ScrollBarPolicy(javafx.scene.control.ScrollPane.ScrollBarPolicy) TextInputDialog(javafx.scene.control.TextInputDialog) Optional(java.util.Optional) Pattern(java.util.regex.Pattern) ObservableList(javafx.collections.ObservableList) IntStream(java.util.stream.IntStream) Scene(javafx.scene.Scene) ListView(javafx.scene.control.ListView) ButtonType(javafx.scene.control.ButtonType) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) FXCollections(javafx.collections.FXCollections) HashMap(java.util.HashMap) NumberFormat(java.text.NumberFormat) ArrayList(java.util.ArrayList) Interval(htsjdk.samtools.util.Interval) Insets(javafx.geometry.Insets) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) GridPane(javafx.scene.layout.GridPane) Locatable(htsjdk.samtools.util.Locatable) KeyCode(javafx.scene.input.KeyCode) Color(javafx.scene.paint.Color) Label(javafx.scene.control.Label) MenuBar(javafx.scene.control.MenuBar) File(java.io.File) Menu(javafx.scene.control.Menu) KeyCodeCombination(javafx.scene.input.KeyCodeCombination) FileChooser(javafx.stage.FileChooser) SelectionMode(javafx.scene.control.SelectionMode) Stage(javafx.stage.Stage) BufferedReader(java.io.BufferedReader) FileChooser(javafx.stage.FileChooser) BufferedReader(java.io.BufferedReader) File(java.io.File) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) Interval(htsjdk.samtools.util.Interval)

Aggregations

BedLineCodec (com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec)22 BufferedReader (java.io.BufferedReader)17 BedLine (com.github.lindenb.jvarkit.util.bio.bed.BedLine)16 Interval (htsjdk.samtools.util.Interval)12 File (java.io.File)12 ArrayList (java.util.ArrayList)11 List (java.util.List)10 Logger (com.github.lindenb.jvarkit.util.log.Logger)9 CloserUtil (htsjdk.samtools.util.CloserUtil)9 Parameter (com.beust.jcommander.Parameter)8 IOUtils (com.github.lindenb.jvarkit.io.IOUtils)8 Program (com.github.lindenb.jvarkit.util.jcommander.Program)8 SAMFileHeader (htsjdk.samtools.SAMFileHeader)8 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)8 SamReader (htsjdk.samtools.SamReader)8 HashMap (java.util.HashMap)8 Map (java.util.Map)8 Launcher (com.github.lindenb.jvarkit.util.jcommander.Launcher)7 StringUtil (htsjdk.samtools.util.StringUtil)7 VCFHeader (htsjdk.variant.vcf.VCFHeader)7