Search in sources :

Example 6 with IntervalParser

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

the class VCFMerge method workUsingSortingCollection.

private int workUsingSortingCollection() {
    VariantContextWriter w = null;
    SortingCollection<VariantOfFile> array = null;
    InputStream in = null;
    CloseableIterator<VariantOfFile> iter = null;
    try {
        final List<String> IN = new ArrayList<String>(this.userVcfFiles);
        final Set<String> genotypeSampleNames = new TreeSet<String>();
        final Set<VCFHeaderLine> metaData = new HashSet<VCFHeaderLine>();
        array = SortingCollection.newInstance(VariantOfFile.class, new VariantCodec(), new VariantComparator(), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
        array.setDestructiveIteration(true);
        for (int fileIndex = 0; fileIndex < IN.size(); ++fileIndex) {
            final String vcfFile = IN.get(fileIndex);
            LOG.info("reading from " + vcfFile + " " + (fileIndex + 1) + "/" + IN.size());
            final VCFHandler handler = new VCFHandler(vcfFile);
            vcfHandlers.add(handler);
            in = IOUtils.openURIForReading(vcfFile);
            final LineReader lr = new SynchronousLineReader(in);
            final LineIterator lit = new LineIteratorImpl(lr);
            handler.header = (VCFHeader) handler.vcfCodec.readActualHeader(lit);
            final SAMSequenceDictionary dict1 = handler.header.getSequenceDictionary();
            if (dict1 == null)
                throw new RuntimeException("dictionary missing in " + vcfFile);
            if (dict1.isEmpty())
                throw new RuntimeException("dictionary is Empty in " + vcfFile);
            genotypeSampleNames.addAll(handler.header.getSampleNamesInOrder());
            metaData.addAll(handler.header.getMetaDataInInputOrder());
            if (fileIndex == 0) {
                this.global_dictionary = dict1;
            } else if (!SequenceUtil.areSequenceDictionariesEqual(global_dictionary, dict1)) {
                throw new JvarkitException.DictionariesAreNotTheSame(global_dictionary, dict1);
            }
            final Predicate<VariantOfFile> accept;
            if (!StringUtil.isBlank(VCFMerge.this.regionStr)) {
                final IntervalParser intervalParser = new IntervalParser(dict1);
                intervalParser.setContigNameIsWholeContig(true);
                final Interval rgn = intervalParser.parse(VCFMerge.this.regionStr);
                accept = (VOL) -> {
                    final VariantContext ctx = VOL.parse();
                    return rgn.intersects(new Interval(ctx.getContig(), ctx.getStart(), ctx.getEnd()));
                };
            } else {
                accept = (VOL) -> true;
            }
            while (lit.hasNext()) {
                final VariantOfFile vof = new VariantOfFile();
                vof.fileIndex = fileIndex;
                vof.line = lit.next();
                if (!accept.test(vof))
                    continue;
                array.add(vof);
            }
            in.close();
            in = null;
        }
        array.doneAdding();
        LOG.info("merging..." + vcfHandlers.size() + " vcfs");
        /* CREATE THE NEW VCH Header */
        VCFHeader mergeHeader = null;
        if (this.doNotMergeRowLines) {
            metaData.add(NO_MERGE_INFO_HEADER);
        }
        mergeHeader = new VCFHeader(metaData, genotypeSampleNames);
        // create the context writer
        w = super.openVariantContextWriter(outputFile);
        w.writeHeader(mergeHeader);
        iter = array.iterator();
        final List<VariantOfFile> row = new ArrayList<VariantOfFile>();
        for (; ; ) {
            VariantOfFile var = null;
            if (iter.hasNext()) {
                var = iter.next();
            } else {
                LOG.info("end of iteration");
                if (!row.isEmpty()) {
                    for (final VariantContext merged : buildContextFromVariantOfFiles(mergeHeader, row)) {
                        w.add(merged);
                    }
                }
                break;
            }
            if (!row.isEmpty() && !row.get(0).same(var)) {
                for (final VariantContext merged : buildContextFromVariantOfFiles(mergeHeader, row)) {
                    w.add(merged);
                }
                row.clear();
            }
            row.add(var);
        }
        CloserUtil.close(w);
        w = null;
        array.cleanup();
        array = null;
        CloserUtil.close(iter);
        iter = null;
        LOG.info("done");
        return RETURN_OK;
    } catch (Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        this.userVcfFiles.clear();
        CloserUtil.close(w);
        CloserUtil.close(in);
        CloserUtil.close(iter);
        if (array != null)
            array.cleanup();
    }
}
Also used : VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) ArrayList(java.util.ArrayList) VariantContext(htsjdk.variant.variantcontext.VariantContext) LineIterator(htsjdk.tribble.readers.LineIterator) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) TreeSet(java.util.TreeSet) SynchronousLineReader(htsjdk.tribble.readers.SynchronousLineReader) LineReader(htsjdk.tribble.readers.LineReader) SynchronousLineReader(htsjdk.tribble.readers.SynchronousLineReader) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) LineIteratorImpl(htsjdk.tribble.readers.LineIteratorImpl) VCFHeader(htsjdk.variant.vcf.VCFHeader) HashSet(java.util.HashSet) IntervalParser(com.github.lindenb.jvarkit.util.bio.IntervalParser) DataInputStream(java.io.DataInputStream) InputStream(java.io.InputStream) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) IOException(java.io.IOException) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) Interval(htsjdk.samtools.util.Interval)

Example 7 with IntervalParser

use of com.github.lindenb.jvarkit.util.bio.IntervalParser 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;
    }
}
Also used : SamReader(htsjdk.samtools.SamReader) IntervalParser(com.github.lindenb.jvarkit.util.bio.IntervalParser) SamReaderFactory(htsjdk.samtools.SamReaderFactory) SAMFileHeader(htsjdk.samtools.SAMFileHeader) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile)

Example 8 with IntervalParser

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

the class KnimeVariantHelper method parseVariantIntervalFilters.

public Predicate<VariantContext> parseVariantIntervalFilters(final String... array) {
    final IntervalParser parser = new IntervalParser();
    parser.setRaiseExceptionOnError(true);
    Predicate<VariantContext> filter = V -> false;
    for (final String str : array) {
        final Interval interval = parser.parse(str);
        filter = filter.or(V -> (V.getContig().equals(interval.getContig()) && !(interval.getEnd() < V.getStart() || V.getEnd() < interval.getStart())));
    }
    return filter;
}
Also used : Allele(htsjdk.variant.variantcontext.Allele) Arrays(java.util.Arrays) Program(com.github.lindenb.jvarkit.util.jcommander.Program) IOUtil(htsjdk.samtools.util.IOUtil) VCFFileReader(htsjdk.variant.vcf.VCFFileReader) VCFHeader(htsjdk.variant.vcf.VCFHeader) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) ParsingUtils(htsjdk.tribble.util.ParsingUtils) VariantContextWriterBuilder(htsjdk.variant.variantcontext.writer.VariantContextWriterBuilder) VCFEncoder(htsjdk.variant.vcf.VCFEncoder) IntervalParser(com.github.lindenb.jvarkit.util.bio.IntervalParser) TabixUtils(htsjdk.tribble.util.TabixUtils) VcfToTable(com.github.lindenb.jvarkit.tools.misc.VcfToTable) IndexFactory(htsjdk.tribble.index.IndexFactory) Map(java.util.Map) CloserUtil(htsjdk.samtools.util.CloserUtil) IntervalTreeMap(htsjdk.samtools.util.IntervalTreeMap) GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder) Predicate(java.util.function.Predicate) Logger(com.github.lindenb.jvarkit.util.log.Logger) Set(java.util.Set) Collectors(java.util.stream.Collectors) List(java.util.List) Stream(java.util.stream.Stream) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VcfTools(com.github.lindenb.jvarkit.util.vcf.VcfTools) VariantContext(htsjdk.variant.variantcontext.VariantContext) Pattern(java.util.regex.Pattern) VcfPeekVcf(com.github.lindenb.jvarkit.tools.vcfvcf.VcfPeekVcf) BedLine(com.github.lindenb.jvarkit.util.bio.bed.BedLine) TabixIndex(htsjdk.tribble.index.tabix.TabixIndex) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) Genotype(htsjdk.variant.variantcontext.Genotype) VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) VCFUtils(com.github.lindenb.jvarkit.util.vcf.VCFUtils) CloseableIterator(htsjdk.samtools.util.CloseableIterator) GenotypeLikelihoods(htsjdk.variant.variantcontext.GenotypeLikelihoods) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) HashMap(java.util.HashMap) Function(java.util.function.Function) ArrayList(java.util.ArrayList) HashSet(java.util.HashSet) Interval(htsjdk.samtools.util.Interval) StringTokenizer(java.util.StringTokenizer) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) StreamSupport(java.util.stream.StreamSupport) VCFCodec(htsjdk.variant.vcf.VCFCodec) Index(htsjdk.tribble.index.Index) VCFConstants(htsjdk.variant.vcf.VCFConstants) IndexedVcfFileReader(com.github.lindenb.jvarkit.util.vcf.IndexedVcfFileReader) VcfIterator(com.github.lindenb.jvarkit.util.vcf.VcfIterator) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) IOException(java.io.IOException) File(java.io.File) IndexedBedReader(com.github.lindenb.jvarkit.util.bio.bed.IndexedBedReader) IndexType(htsjdk.tribble.index.IndexFactory.IndexType) BufferedReader(java.io.BufferedReader) Collections(java.util.Collections) Tribble(htsjdk.tribble.Tribble) IntervalParser(com.github.lindenb.jvarkit.util.bio.IntervalParser) VariantContext(htsjdk.variant.variantcontext.VariantContext) Interval(htsjdk.samtools.util.Interval)

Example 9 with IntervalParser

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

the class VcfLoopOverGenes method doWork.

@SuppressWarnings("resource")
@Override
public int doWork(final List<String> args) {
    PrintWriter pw = null;
    VCFFileReader vcfFileReader = null;
    CloseableIterator<VariantContext> iter = null;
    CloseableIterator<GeneLoc> iter2 = null;
    BufferedReader br = null;
    ArchiveFactory archive = null;
    try {
        final File vcf = new File(oneAndOnlyOneFile(args));
        vcfFileReader = new VCFFileReader(vcf, (this.geneFile != null || !StringUtil.isBlank(this.regionStr)));
        this.dictionary = vcfFileReader.getFileHeader().getSequenceDictionary();
        if (this.dictionary == null) {
            throw new JvarkitException.VcfDictionaryMissing(vcf);
        }
        final VcfTools tools = new VcfTools(vcfFileReader.getFileHeader());
        if (!this.prefix.isEmpty() && !this.prefix.endsWith(".")) {
            this.prefix += ".";
        }
        if (this.geneFile == null) {
            final SortingCollection<GeneLoc> sortingCollection = SortingCollection.newInstance(GeneLoc.class, new GeneLocCodec(), (A, B) -> A.compareTo(B), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
            sortingCollection.setDestructiveIteration(true);
            if (StringUtil.isBlank(this.regionStr)) {
                iter = vcfFileReader.iterator();
            } else {
                final IntervalParser parser = new IntervalParser(this.dictionary);
                parser.setContigNameIsWholeContig(true);
                final Interval interval = parser.parse(this.regionStr);
                if (interval == null) {
                    LOG.error("Cannot parse interval " + this.regionStr);
                    return -1;
                }
                iter = vcfFileReader.query(interval.getContig(), interval.getStart(), interval.getEnd());
            }
            final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(vcfFileReader.getFileHeader()).logger(LOG);
            if (this.splitMethod.equals(SplitMethod.Annotations)) {
                while (iter.hasNext()) {
                    final VariantContext ctx = progress.watch(iter.next());
                    for (final AnnPredictionParser.AnnPrediction pred : tools.getAnnPredictionParser().getPredictions(ctx)) {
                        if (this.snpEffNoIntergenic && pred.isIntergenicRegion()) {
                            continue;
                        }
                        if (!StringUtil.isBlank(pred.getGeneName())) {
                            sortingCollection.add(create(ctx, pred.getGeneName(), SourceType.ANN_GeneName));
                        }
                        if (!StringUtil.isBlank(pred.getGeneId())) {
                            sortingCollection.add(create(ctx, pred.getGeneId(), SourceType.ANN_GeneID));
                        }
                        if (!StringUtil.isBlank(pred.getFeatureId())) {
                            sortingCollection.add(create(ctx, pred.getFeatureId(), SourceType.ANN_FeatureID));
                        }
                    }
                    for (final VepPredictionParser.VepPrediction pred : tools.getVepPredictionParser().getPredictions(ctx)) {
                        if (!StringUtil.isBlank(pred.getGene())) {
                            sortingCollection.add(create(ctx, pred.getGene(), SourceType.VEP_Gene));
                        }
                        if (!StringUtil.isBlank(pred.getFeature())) {
                            sortingCollection.add(create(ctx, pred.getFeature(), SourceType.VEP_Feature));
                        }
                        if (!StringUtil.isBlank(pred.getSymbol())) {
                            sortingCollection.add(create(ctx, pred.getSymbol(), SourceType.VEP_Symbol));
                        }
                        if (!StringUtil.isBlank(pred.getHgncId())) {
                            sortingCollection.add(create(ctx, pred.getHgncId(), SourceType.VEP_HgncId));
                        }
                    }
                }
            } else /**
             * split VCF per sliding window of variants
             */
            if (this.splitMethod.equals(SplitMethod.VariantSlidingWindow)) {
                if (this.variantsWinCount < 1) {
                    LOG.error("Bad value for variantsWinCount");
                    return -1;
                }
                if (this.variantsWinShift < 1 || this.variantsWinShift > this.variantsWinCount) {
                    LOG.error("Bad value for variantsWinShift");
                    return -1;
                }
                final List<VariantContext> buffer = new ArrayList<>(this.variantsWinCount);
                /**
                 * routine to dump buffer into sorting collection
                 */
                final Runnable dumpBuffer = () -> {
                    if (buffer.isEmpty())
                        return;
                    final String contig = buffer.get(0).getContig();
                    final int chromStart = buffer.stream().mapToInt(CTX -> CTX.getStart()).min().getAsInt();
                    // use last of start too
                    final int chromEnd0 = buffer.stream().mapToInt(CTX -> CTX.getStart()).max().getAsInt();
                    // final int chromEnd1 = buffer.stream().mapToInt(CTX->CTX.getEnd()).max().getAsInt();
                    final String identifier = contig + "_" + String.format(NUM_FORMAT, chromStart) + "_" + String.format(NUM_FORMAT, chromEnd0);
                    for (final VariantContext ctx : buffer) {
                        sortingCollection.add(create(ctx, identifier, SourceType.SlidingVariants));
                    }
                };
                while (iter.hasNext()) {
                    VariantContext ctx = progress.watch(iter.next());
                    /* reduce the memory footprint for this context */
                    ctx = new VariantContextBuilder(ctx).genotypes(Collections.emptyList()).unfiltered().rmAttributes(new ArrayList<>(ctx.getAttributes().keySet())).make();
                    if (!buffer.isEmpty() && !buffer.get(0).getContig().equals(ctx.getContig())) {
                        dumpBuffer.run();
                        buffer.clear();
                    }
                    buffer.add(ctx);
                    if (buffer.size() >= this.variantsWinCount) {
                        dumpBuffer.run();
                        final int fromIndex = Math.min(this.variantsWinShift, buffer.size());
                        buffer.subList(0, fromIndex).clear();
                    }
                }
                dumpBuffer.run();
                buffer.clear();
            } else if (this.splitMethod.equals(SplitMethod.ContigSlidingWindow)) {
                if (this.contigWinLength < 1) {
                    LOG.error("Bad value for contigWinCount");
                    return -1;
                }
                if (this.contigWinShift < 1 || this.contigWinShift > this.contigWinLength) {
                    LOG.error("Bad value for contigWinShift");
                    return -1;
                }
                while (iter.hasNext()) {
                    VariantContext ctx = progress.watch(iter.next());
                    /* reduce the memory footprint for this context */
                    ctx = new VariantContextBuilder(ctx).genotypes(Collections.emptyList()).unfiltered().rmAttributes(new ArrayList<>(ctx.getAttributes().keySet())).make();
                    int start = 0;
                    while (start <= ctx.getStart()) {
                        if (start + this.contigWinLength >= ctx.getStart()) {
                            final int chromStart = start;
                            final int chromEnd0 = start + this.contigWinLength;
                            final String identifier = ctx.getContig() + "_" + String.format(NUM_FORMAT, chromStart) + "_" + String.format(NUM_FORMAT, chromEnd0);
                            sortingCollection.add(create(ctx, identifier, SourceType.SlidingContig));
                        }
                        start += this.contigWinShift;
                    }
                }
            } else {
                throw new IllegalStateException("No such method: " + this.splitMethod);
            }
            sortingCollection.doneAdding();
            progress.finish();
            iter.close();
            iter = null;
            pw = super.openFileOrStdoutAsPrintWriter(this.outputFile);
            iter2 = sortingCollection.iterator();
            final EqualRangeIterator<GeneLoc> eqiter = new EqualRangeIterator<>(iter2, this.compareGeneName);
            int geneIdentifierId = 0;
            while (eqiter.hasNext()) {
                final List<GeneLoc> gene = eqiter.next();
                pw.print(gene.get(0).contig);
                pw.print('\t');
                // -1 for BED
                pw.print(gene.stream().mapToInt(G -> G.start).min().getAsInt() - 1);
                pw.print('\t');
                pw.print(gene.stream().mapToInt(G -> G.end).max().getAsInt());
                pw.print('\t');
                pw.print(this.prefix + String.format("%09d", ++geneIdentifierId));
                pw.print('\t');
                pw.print(gene.get(0).geneName);
                pw.print('\t');
                pw.print(gene.get(0).sourceType);
                pw.print('\t');
                pw.print(gene.size());
                pw.println();
            }
            pw.flush();
            pw.close();
            pw = null;
            eqiter.close();
            iter2.close();
            iter2 = null;
            sortingCollection.cleanup();
        } else {
            if (this.nJobs < 1) {
                this.nJobs = Math.max(1, Runtime.getRuntime().availableProcessors());
                LOG.info("setting njobs to " + this.nJobs);
            }
            final ExecutorService executorService;
            final List<Future<Integer>> futureResults;
            if (this.nJobs > 1) {
                executorService = new ThreadPoolExecutor(this.nJobs, this.nJobs, 0L, TimeUnit.MILLISECONDS, new LinkedBlockingQueue<Runnable>());
                futureResults = new ArrayList<>();
            } else {
                executorService = null;
                futureResults = Collections.emptyList();
            }
            if (this.outputFile == null) {
                LOG.error("When scanning a VCF with " + this.geneFile + ". Output file must be defined");
            }
            if (!this.exec.isEmpty()) {
                if (this.outputFile.getName().endsWith(".zip")) {
                    LOG.error("Cannot execute " + this.exec + " when saving to a zip.");
                    return -1;
                }
            }
            archive = ArchiveFactory.open(this.outputFile);
            PrintWriter manifest = this.deleteAfterCommand && !this.exec.isEmpty() ? // all files will be deleted, no manifest needed
            new PrintWriter(new NullOuputStream()) : archive.openWriter(this.prefix + "manifest.txt");
            br = IOUtils.openFileForBufferedReading(this.geneFile);
            final BedLineCodec bedCodec = new BedLineCodec();
            for (; ; ) {
                if (!futureResults.isEmpty()) {
                    int i = 0;
                    while (i < futureResults.size()) {
                        final Future<Integer> r = futureResults.get(i);
                        if (r.isCancelled()) {
                            LOG.error("Task was canceled. Break.");
                            return -1;
                        } else if (r.isDone()) {
                            futureResults.remove(i);
                            int rez = r.get();
                            if (rez != 0) {
                                LOG.error("Task Failed (" + rez + "). Break");
                            }
                        } else {
                            i++;
                        }
                    }
                }
                final String line = br.readLine();
                if (line == null)
                    break;
                if (line.startsWith("#") || line.isEmpty())
                    continue;
                final BedLine bedLine = bedCodec.decode(line);
                if (bedLine == null)
                    continue;
                // ID
                final String geneIdentifier = bedLine.get(3);
                // name
                final String geneName = bedLine.get(4);
                final SourceType sourceType = SourceType.valueOf(bedLine.get(5));
                final String filename = geneIdentifier;
                final String outputVcfName = (filename.startsWith(this.prefix) ? "" : this.prefix) + filename + ".vcf" + (this.compress ? ".gz" : "");
                LOG.info(bedLine.getContig() + ":" + bedLine.getStart() + "-" + bedLine.getEnd() + " length :" + (bedLine.getEnd() - bedLine.getStart()));
                if (bedLine.getEnd() - bedLine.getStart() > 1E6) {
                    LOG.warn("That's a large region ! " + bedLine);
                }
                OutputStream vcfOutputStream = null;
                VariantContextWriter vw = null;
                int countVariants = 0;
                final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(vcfFileReader.getFileHeader()).logger(LOG).prefix(geneName + " " + bedLine.getContig() + ":" + bedLine.getStart() + "-" + bedLine.getEnd());
                iter = vcfFileReader.query(bedLine.getContig(), bedLine.getStart(), bedLine.getEnd());
                while (iter.hasNext()) {
                    VariantContext ctx = progress.watch(iter.next());
                    switch(sourceType) {
                        case SlidingVariants:
                            {
                                // nothing
                                break;
                            }
                        case SlidingContig:
                            {
                                // nothing
                                break;
                            }
                        case ANN_GeneName:
                        case ANN_FeatureID:
                        case ANN_GeneID:
                            {
                                final List<String> preds = new ArrayList<>();
                                for (final AnnPredictionParser.AnnPrediction pred : tools.getAnnPredictionParser().getPredictions(ctx)) {
                                    final String predictionIdentifier;
                                    switch(sourceType) {
                                        case ANN_GeneName:
                                            predictionIdentifier = pred.getGeneName();
                                            break;
                                        case ANN_FeatureID:
                                            predictionIdentifier = pred.getFeatureId();
                                            break;
                                        case ANN_GeneID:
                                            predictionIdentifier = pred.getGeneId();
                                            break;
                                        default:
                                            throw new IllegalStateException(bedLine.toString());
                                    }
                                    if (StringUtil.isBlank(predictionIdentifier))
                                        continue;
                                    if (!geneName.equals(predictionIdentifier))
                                        continue;
                                    preds.add(pred.getOriginalAttributeAsString());
                                }
                                if (preds.isEmpty()) {
                                    ctx = null;
                                } else {
                                    ctx = new VariantContextBuilder(ctx).rmAttribute(tools.getAnnPredictionParser().getTag()).attribute(tools.getAnnPredictionParser().getTag(), preds).make();
                                }
                                break;
                            }
                        case VEP_Gene:
                        case VEP_Feature:
                        case VEP_Symbol:
                        case VEP_HgncId:
                            {
                                final List<String> preds = new ArrayList<>();
                                for (final VepPredictionParser.VepPrediction pred : tools.getVepPredictions(ctx)) {
                                    final String predictionIdentifier;
                                    switch(sourceType) {
                                        case VEP_Gene:
                                            predictionIdentifier = pred.getGene();
                                            break;
                                        case VEP_Feature:
                                            predictionIdentifier = pred.getFeature();
                                            break;
                                        case VEP_Symbol:
                                            predictionIdentifier = pred.getSymbol();
                                            break;
                                        case VEP_HgncId:
                                            predictionIdentifier = pred.getHgncId();
                                            break;
                                        default:
                                            throw new IllegalStateException(bedLine.toString());
                                    }
                                    if (StringUtil.isBlank(predictionIdentifier))
                                        continue;
                                    if (!geneName.equals(predictionIdentifier))
                                        continue;
                                    preds.add(pred.getOriginalAttributeAsString());
                                }
                                if (preds.isEmpty()) {
                                    ctx = null;
                                } else {
                                    ctx = new VariantContextBuilder(ctx).rmAttribute(tools.getVepPredictionParser().getTag()).attribute(tools.getVepPredictionParser().getTag(), preds).make();
                                }
                                break;
                            }
                        default:
                            throw new IllegalStateException(bedLine.toString());
                    }
                    if (ctx == null)
                        continue;
                    if (vcfOutputStream == null) {
                        LOG.info(filename);
                        manifest.println(outputVcfName);
                        final VCFHeader header = new VCFHeader(vcfFileReader.getFileHeader());
                        header.addMetaDataLine(new VCFHeaderLine(VCF_HEADER_SPLITKEY, filename));
                        vcfOutputStream = archive.openOuputStream(outputVcfName);
                        vw = VCFUtils.createVariantContextWriterToOutputStream(vcfOutputStream);
                        vw.writeHeader(header);
                    }
                    countVariants++;
                    vw.add(ctx);
                    if (countVariants % 1000 == 0) {
                        LOG.info("Loading : " + geneIdentifier + " N=" + countVariants);
                    }
                }
                progress.finish();
                LOG.info(geneIdentifier + " N=" + countVariants);
                if (vcfOutputStream != null) {
                    vw.close();
                    vcfOutputStream.flush();
                    vcfOutputStream.close();
                    vw = null;
                    if (!this.exec.isEmpty()) {
                        final Callable<Integer> callable = () -> {
                            final File vcfOutFile = new File(this.outputFile, outputVcfName);
                            IOUtil.assertFileIsReadable(vcfOutFile);
                            final String vcfPath = vcfOutFile.getPath();
                            final StringTokenizer st = new StringTokenizer(this.exec);
                            final List<String> command = new ArrayList<>(1 + st.countTokens());
                            while (st.hasMoreTokens()) {
                                String token = st.nextToken().replaceAll("__PREFIX__", this.prefix).replaceAll("__CONTIG__", bedLine.getContig()).replaceAll("__CHROM__", bedLine.getContig()).replaceAll("__ID__", geneIdentifier).replaceAll("__NAME__", geneName).replaceAll("__START__", String.valueOf(bedLine.getStart())).replaceAll("__END__", String.valueOf(bedLine.getEnd())).replaceAll("__SOURCE__", sourceType.name()).replaceAll("__VCF__", vcfPath);
                                command.add(token);
                            }
                            LOG.info(command.stream().map(S -> "'" + S + "'").collect(Collectors.joining(" ")));
                            final ProcessBuilder pb = new ProcessBuilder(command);
                            pb.redirectErrorStream(true);
                            final Process p = pb.start();
                            final Thread stdoutThread = new Thread(() -> {
                                try {
                                    InputStream in = p.getInputStream();
                                    IOUtils.copyTo(in, stdout());
                                } catch (Exception err) {
                                    LOG.error(err);
                                }
                            });
                            stdoutThread.start();
                            int exitValue = p.waitFor();
                            if (exitValue != 0) {
                                LOG.error("Command failed (" + exitValue + "):" + String.join(" ", command));
                                return -1;
                            } else {
                                if (deleteAfterCommand) {
                                    if (!vcfOutFile.delete()) {
                                        LOG.warn("Cannot delete " + vcfOutFile);
                                    }
                                }
                                return 0;
                            }
                        };
                        if (executorService != null) {
                            final Future<Integer> rez = executorService.submit(callable);
                            futureResults.add(rez);
                        } else {
                            final int ret = callable.call();
                            if (ret != 0) {
                                LOG.error("Error with process (" + ret + ")");
                                return ret;
                            }
                        }
                    }
                } else {
                    manifest.println("#" + filename);
                    LOG.warn("No Variant Found for " + line);
                }
                iter.close();
            }
            ;
            if (executorService != null) {
                LOG.info("shutdown");
                executorService.shutdown();
                executorService.awaitTermination(365, TimeUnit.DAYS);
            }
            br.close();
            br = null;
            manifest.close();
            archive.close();
            archive = null;
            LOG.info("Done");
        }
        vcfFileReader.close();
        vcfFileReader = null;
        return 0;
    } catch (Exception e) {
        LOG.error(e);
        return -1;
    } finally {
        {
            CloserUtil.close(iter2);
            CloserUtil.close(iter);
            CloserUtil.close(pw);
            CloserUtil.close(vcfFileReader);
            CloserUtil.close(br);
            CloserUtil.close(archive);
        }
    }
}
Also used : Program(com.github.lindenb.jvarkit.util.jcommander.Program) IOUtil(htsjdk.samtools.util.IOUtil) VCFFileReader(htsjdk.variant.vcf.VCFFileReader) VCFHeader(htsjdk.variant.vcf.VCFHeader) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) IntervalParser(com.github.lindenb.jvarkit.util.bio.IntervalParser) Future(java.util.concurrent.Future) DataOutputStream(java.io.DataOutputStream) StringUtil(htsjdk.samtools.util.StringUtil) AbstractDataCodec(com.github.lindenb.jvarkit.util.picard.AbstractDataCodec) CloserUtil(htsjdk.samtools.util.CloserUtil) PrintWriter(java.io.PrintWriter) Logger(com.github.lindenb.jvarkit.util.log.Logger) LinkedBlockingQueue(java.util.concurrent.LinkedBlockingQueue) Collectors(java.util.stream.Collectors) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) List(java.util.List) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VcfTools(com.github.lindenb.jvarkit.util.vcf.VcfTools) VariantContext(htsjdk.variant.variantcontext.VariantContext) BedLine(com.github.lindenb.jvarkit.util.bio.bed.BedLine) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) DataInputStream(java.io.DataInputStream) VCFUtils(com.github.lindenb.jvarkit.util.vcf.VCFUtils) CloseableIterator(htsjdk.samtools.util.CloseableIterator) ThreadPoolExecutor(java.util.concurrent.ThreadPoolExecutor) Parameter(com.beust.jcommander.Parameter) NullOuputStream(com.github.lindenb.jvarkit.io.NullOuputStream) AnnPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParser) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) Callable(java.util.concurrent.Callable) Function(java.util.function.Function) ParametersDelegate(com.beust.jcommander.ParametersDelegate) ArrayList(java.util.ArrayList) Interval(htsjdk.samtools.util.Interval) StringTokenizer(java.util.StringTokenizer) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) VepPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.VepPredictionParser) ExecutorService(java.util.concurrent.ExecutorService) OutputStream(java.io.OutputStream) SortingCollection(htsjdk.samtools.util.SortingCollection) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) IOException(java.io.IOException) File(java.io.File) TimeUnit(java.util.concurrent.TimeUnit) EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator) BufferedReader(java.io.BufferedReader) Comparator(java.util.Comparator) Collections(java.util.Collections) ArchiveFactory(com.github.lindenb.jvarkit.io.ArchiveFactory) InputStream(java.io.InputStream) VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) DataOutputStream(java.io.DataOutputStream) OutputStream(java.io.OutputStream) ArrayList(java.util.ArrayList) LinkedBlockingQueue(java.util.concurrent.LinkedBlockingQueue) EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator) List(java.util.List) ArrayList(java.util.ArrayList) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) BedLine(com.github.lindenb.jvarkit.util.bio.bed.BedLine) VcfTools(com.github.lindenb.jvarkit.util.vcf.VcfTools) File(java.io.File) Interval(htsjdk.samtools.util.Interval) AnnPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParser) VCFFileReader(htsjdk.variant.vcf.VCFFileReader) VariantContext(htsjdk.variant.variantcontext.VariantContext) NullOuputStream(com.github.lindenb.jvarkit.io.NullOuputStream) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VCFHeader(htsjdk.variant.vcf.VCFHeader) PrintWriter(java.io.PrintWriter) ArchiveFactory(com.github.lindenb.jvarkit.io.ArchiveFactory) IntervalParser(com.github.lindenb.jvarkit.util.bio.IntervalParser) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) DataInputStream(java.io.DataInputStream) InputStream(java.io.InputStream) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) IOException(java.io.IOException) StringTokenizer(java.util.StringTokenizer) VepPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.VepPredictionParser) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) BufferedReader(java.io.BufferedReader) ExecutorService(java.util.concurrent.ExecutorService) Future(java.util.concurrent.Future) ThreadPoolExecutor(java.util.concurrent.ThreadPoolExecutor)

Example 10 with IntervalParser

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

the class LowResBam2Raster 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;
    }
    if (this.gcWinSize <= 0) {
        LOG.info("adjusting GC win size to 5");
        this.gcWinSize = 5;
    }
    SamReader samFileReader = null;
    try {
        if (this.referenceFile != null) {
            LOG.info("loading reference");
            this.indexedFastaSequenceFile = new IndexedFastaSequenceFile(this.referenceFile);
        }
        final SamReaderFactory srf = super.createSamReaderFactory();
        this.interval = new IntervalParser(this.indexedFastaSequenceFile == null ? null : this.indexedFastaSequenceFile.getSequenceDictionary()).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();
        if (this.knownGeneUrl != null) {
            IntervalTreeMap<List<KnownGene>> map = KnownGene.loadUriAsIntervalTreeMap(this.knownGeneUrl, (KG) -> (KG.getContig().equals(this.interval.getContig()) && !(KG.getEnd() < this.interval.getStart() || KG.getStart() + 1 > this.interval.getEnd())));
            this.knownGenes.addAll(map.values().stream().flatMap(L -> L.stream()).collect(Collectors.toList()));
        }
        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. no Read-Group specified ? no data in that region ?");
            return -1;
        }
        this.key2partition.values().stream().forEach(P -> P.make());
        saveImages(this.key2partition.values().stream().map(P -> P.image).collect(Collectors.toList()));
        return RETURN_OK;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(samFileReader);
    }
}
Also used : SamReader(htsjdk.samtools.SamReader) IntervalParser(com.github.lindenb.jvarkit.util.bio.IntervalParser) SamReaderFactory(htsjdk.samtools.SamReaderFactory) ArrayList(java.util.ArrayList) List(java.util.List) SAMFileHeader(htsjdk.samtools.SAMFileHeader) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile)

Aggregations

IntervalParser (com.github.lindenb.jvarkit.util.bio.IntervalParser)10 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)9 Interval (htsjdk.samtools.util.Interval)8 SamReader (htsjdk.samtools.SamReader)7 ArrayList (java.util.ArrayList)7 SAMSequenceDictionaryProgress (com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress)6 IOException (java.io.IOException)6 File (java.io.File)5 PrintWriter (java.io.PrintWriter)5 JvarkitException (com.github.lindenb.jvarkit.lang.JvarkitException)4 SAMFileHeader (htsjdk.samtools.SAMFileHeader)4 SAMRecord (htsjdk.samtools.SAMRecord)4 IOUtils (com.github.lindenb.jvarkit.io.IOUtils)3 Program (com.github.lindenb.jvarkit.util.jcommander.Program)3 Logger (com.github.lindenb.jvarkit.util.log.Logger)3 SAMRecordIterator (htsjdk.samtools.SAMRecordIterator)3 SamReaderFactory (htsjdk.samtools.SamReaderFactory)3 IndexedFastaSequenceFile (htsjdk.samtools.reference.IndexedFastaSequenceFile)3 CloseableIterator (htsjdk.samtools.util.CloseableIterator)3 VariantContext (htsjdk.variant.variantcontext.VariantContext)3