Examples with ArchiveFactory com.github.lindenb.jvarkit.io.ArchiveFactory used on opensource projects

Example 1 with ArchiveFactory

use of com.github.lindenb.jvarkit.io.ArchiveFactory in project jvarkit by lindenb.

the class CaseControlPlot method doWork.

@Override
public int doWork(final List<String> args) {
    ArchiveFactory archiveFactory = null;
    VcfIterator in = null;
    VariantContextWriter teeVariantWriter = null;
    final List<CaseControlExtractor> excractors = new ArrayList<>();
    try {
        in = super.openVcfIterator(oneFileOrNull(args));
        final VCFHeader header = in.getHeader();
        excractors.addAll(parseConfigFile(header));
        final Pedigree pedigree;
        if (this.pedigreeFile != null) {
            pedigree = Pedigree.newParser().parse(this.pedigreeFile);
        } else {
            pedigree = Pedigree.newParser().parse(header);
        }
        if (pedigree == null || pedigree.isEmpty()) {
            LOG.error("No pedigree defined , or it is empty");
            return -1;
        }
        final Set<Pedigree.Person> casepersons = pedigree.getPersons().stream().filter(F -> F.isAffected() && header.getSampleNamesInOrder().indexOf(F.getId()) != -1).collect(Collectors.toSet());
        if (casepersons.isEmpty()) {
            LOG.error("No Affected individuals in pedigree/header");
            return -1;
        }
        final Set<Pedigree.Person> controlpersons = pedigree.getPersons().stream().filter(F -> F.isUnaffected() && header.getSampleNamesInOrder().indexOf(F.getId()) != -1).collect(Collectors.toSet());
        if (controlpersons.isEmpty()) {
            LOG.error("No Unaffected individuals in pedigree/header");
            return -1;
        }
        if (this.teeToStdout) {
            teeVariantWriter = super.openVariantContextWriter(null);
            teeVariantWriter.writeHeader(header);
        }
        archiveFactory = ArchiveFactory.open(this.outputDirOrZip);
        for (final CaseControlExtractor extractor : excractors) {
            extractor.pw = archiveFactory.openWriter(extractor.name);
        }
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(header);
        while (in.hasNext()) {
            final VariantContext ctx = progress.watch(in.next());
            if (teeVariantWriter != null)
                teeVariantWriter.add(ctx);
            for (final CaseControlExtractor handler : excractors) {
                handler.visit(ctx, casepersons, controlpersons);
            }
        }
        for (final CaseControlExtractor handler : excractors) {
            handler.close();
        }
        progress.finish();
        if (teeVariantWriter != null) {
            teeVariantWriter.close();
            teeVariantWriter = null;
        }
        in.close();
        in = null;
        archiveFactory.close();
        archiveFactory = null;
        return RETURN_OK;
    } catch (final Exception e) {
        LOG.error(e);
        return -1;
    } finally {
        CloserUtil.close(in);
        CloserUtil.close(archiveFactory);
        CloserUtil.close(teeVariantWriter);
        for (final CaseControlExtractor handler : excractors) {
            handler.close();
        }
    }
}

Example 2 with ArchiveFactory

use of com.github.lindenb.jvarkit.io.ArchiveFactory in project jvarkit by lindenb.

the class CopyNumber01 method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.refFile == null) {
        LOG.error("Undefined REF file");
        return -1;
    }
    if (this.archiveFile == null) {
        LOG.error("Undefined output file.");
        return -1;
    }
    if (!StringUtil.isBlank(prefix)) {
        if (!(prefix.endsWith(".") || prefix.endsWith("_"))) {
            prefix = prefix + ".";
        }
    }
    SamReader samReader = null;
    ArchiveFactory archive = null;
    try {
        final String input = oneAndOnlyOneFile(args);
        samReader = super.openSamReader(input);
        if (!samReader.hasIndex()) {
            LOG.error("file is not indexed " + input);
            return -1;
        }
        final String sampleName = samReader.getFileHeader().getReadGroups().stream().map(RG -> RG.getSample()).filter(S -> !StringUtil.isBlank(S)).findFirst().orElse("SAMPLE");
        final SAMFileHeader header = samReader.getFileHeader();
        this.samDictionary = header.getSequenceDictionary();
        if (this.samDictionary == null || this.samDictionary.isEmpty()) {
            throw new JvarkitException.DictionaryMissing(input);
        }
        /* loading REF Reference */
        this.indexedFastaSequenceFile = new IndexedFastaSequenceFile(refFile);
        final SAMSequenceDictionary dict = this.indexedFastaSequenceFile.getSequenceDictionary();
        if (dict == null) {
            throw new JvarkitException.DictionaryMissing(refFile.getPath());
        }
        this.sam2faiContigNameConverter = ContigNameConverter.fromDictionaries(this.samDictionary, dict);
        this.sam2faiContigNameConverter.setOnNotFound(OnNotFound.SKIP);
        archive = ArchiveFactory.open(this.archiveFile);
        PrintWriter mkw = archive.openWriter(this.prefix + "cnv01.mk");
        mkw.println("## " + this.getProgramCommandLine());
        mkw.println(".PHONY=all all2");
        mkw.println("TARGETS=");
        mkw.println("SCREEN_WIDTH?=2600");
        mkw.println("SCREEN_HEIGHT?=1000");
        mkw.println("GNUPLOT?=gnuplot");
        mkw.println("all : all2");
        final List<String> all_chrom_files = new ArrayList<>();
        double maxDepth = 0.0;
        for (final SAMSequenceRecord ssr : this.samDictionary.getSequences()) {
            if (!StringUtil.isBlank(this.limitToChrom) && !this.limitToChrom.equals(ssr.getSequenceName())) {
                LOG.info("Skipping " + ssr.getSequenceName() + "....");
                continue;
            }
            if (this.sam2faiContigNameConverter.apply(ssr.getSequenceName()) == null) {
                LOG.info("Ignoring " + ssr.getSequenceName() + " because it's not in REF");
                continue;
            }
            if (ignoreChromosomeName(ssr.getSequenceName())) {
                LOG.info("Ignoring " + ssr.getSequenceName());
                continue;
            }
            LOG.info("processing chromosome " + ssr.getSequenceName());
            ContigProcessor proc = new ContigProcessor(samReader, ssr, sampleName);
            proc.run();
            final String tsvFilename = this.prefix + proc.getContig() + "." + proc.sampleName + ".bed";
            final String pngFilename = "$(addsuffix .png,$(basename " + tsvFilename + "))";
            final double depth = proc.items.stream().mapToDouble(I -> I.norm_depth).max().orElse(4);
            maxDepth = Math.max(maxDepth, depth);
            all_chrom_files.add(tsvFilename);
            mkw.println("TARGETS+=" + pngFilename);
            mkw.println(pngFilename + ":" + tsvFilename);
            mkw.println("\techo 'set terminal png truecolor size ${SCREEN_WIDTH},${SCREEN_HEIGHT};" + "set key autotitle columnhead;" + "set datafile separator \"\t\";" + "set title \"" + proc.sampleName + " " + ssr.getSequenceName() + "\";" + "set ylabel \"Normalized Depth - Number of copies\";" + "set xlabel \"Position on " + ssr.getSequenceName() + "\";" + // + "set style circle radius 0.02;"
            "set nokey;" + "set yrange [0:" + Math.min(4, Math.ceil(depth)) + "];" + "set xrange [1:" + ssr.getSequenceLength() + "];" + "set xtic rotate by 90 right;" + "set output \"$@\";" + // + "plot \"$<\" using 1:2:1 w points linecolor variable "
            "plot \"$<\" using 2:7 w lines " + "' | ${GNUPLOT}");
            PrintWriter pw = archive.openWriter(tsvFilename);
            proc.saveCoverage(pw);
            pw.flush();
            pw.close();
            proc = null;
        }
        samReader.close();
        final String pngFilename = this.prefix + "wholeGenome.png";
        mkw.println("TARGETS+=" + pngFilename);
        mkw.println(pngFilename + ":" + String.join(" ", all_chrom_files));
        mkw.println("\tgrep -v '^#' $^ | cut -f 1,4,7 | sed " + samDictionary.getSequences().stream().map(SSR -> " -e 's/^" + SSR.getSequenceName() + "\t/" + SSR.getSequenceIndex() + "\t/' ").collect(Collectors.joining()) + "| sort -t '\t' -k2,2n > $(addsuffix .tmp.tsv,$@)");
        mkw.println("\techo 'set terminal png truecolor size ${SCREEN_WIDTH},${SCREEN_HEIGHT};" + "set datafile separator \"\t\";" + "set title \"" + sampleName + " (Whole Genome)\";" + "set ylabel \"Normalized Depth - Number of copies\";" + "set xlabel \"Genomic Position\";" + "set yrange [0:" + Math.min(4, Math.ceil(maxDepth)) + "];" + "set xrange [1:" + this.samDictionary.getReferenceLength() + ".0];" + "set xtic rotate by 90 right;" + "set nokey;" + "set output \"$@\";" + "plot \"$(addsuffix .tmp.tsv,$@)\" using 2:3:1 w points linecolor variable " + "' | ${GNUPLOT}");
        mkw.println("\trm -f $(addsuffix .tmp.tsv,$@)");
        mkw.println("all2: ${TARGETS}");
        mkw.flush();
        mkw.close();
        mkw = null;
        archive.close();
        archive = null;
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(samReader);
        CloserUtil.close(archive);
        CloserUtil.close(this.indexedFastaSequenceFile);
    }
}

Example 3 with ArchiveFactory

use of com.github.lindenb.jvarkit.io.ArchiveFactory in project jvarkit by lindenb.

the class VcfCompareCallers method doWork.

@SuppressWarnings("resource")
@Override
public int doWork(final List<String> args) {
    if (this.archiveFile == null) {
        LOG.error("undefined output");
        return -1;
    }
    if (!this.archivePrefix.isEmpty() && !archivePrefix.endsWith(".")) {
        this.archivePrefix = this.archivePrefix + ".";
    }
    if (!this.archiveFile.getName().endsWith(".zip")) {
        IOUtil.assertDirectoryIsWritable(this.archiveFile);
    }
    PeekVCF vcfIterator1 = null;
    PeekVCF vcfIterator2 = null;
    ArchiveFactory archiveFactory = null;
    IntervalTreeMap<Boolean> capture = null;
    PrintWriter makefileWriter = null;
    try {
        if (args.size() == 1) {
            LOG.info("Reading from VCF1=stdin and VCF2=" + args.get(0));
            vcfIterator1 = new PeekVCF(VCFUtils.createVcfIteratorFromInputStream(stdin()), "<STDIN>");
            vcfIterator2 = new PeekVCF(VCFUtils.createVcfIterator(args.get(0)), args.get(0));
        } else if (args.size() == 2) {
            LOG.info("Reading from VCF1=" + args.get(0) + " and VCF2=" + args.get(1));
            vcfIterator1 = new PeekVCF(VCFUtils.createVcfIterator(args.get(0)), args.get(0));
            vcfIterator2 = new PeekVCF(VCFUtils.createVcfIterator(args.get(1)), args.get(1));
        } else {
            LOG.error("illegal number of arguments");
            return -1;
        }
        if (this.captureFile != null) {
            LOG.info("Reading " + this.captureFile);
            capture = super.readBedFileAsBooleanIntervalTreeMap(this.captureFile);
        }
        this.global_dictionary = vcfIterator1.dict;
        if (!SequenceUtil.areSequenceDictionariesEqual(vcfIterator1.dict, vcfIterator2.dict)) {
            throw new JvarkitException.DictionariesAreNotTheSame(vcfIterator1.dict, vcfIterator2.dict);
        }
        for (int side = 0; side < 2; ++side) {
            final List<String> jexlExprStrings = (side == 0 ? this.jexlExprStrings1 : this.jexlExprStrings2);
            final PeekVCF peek = (side == 0 ? vcfIterator1 : vcfIterator2);
            // initialize JEXL map
            if (!jexlExprStrings.isEmpty())
                continue;
            final ArrayList<String> dummyNames = new ArrayList<String>(jexlExprStrings.size());
            for (int expCount = 1; expCount < jexlExprStrings.size(); ++expCount) {
                dummyNames.add(String.format("vce%d", expCount));
            }
            peek.jexlVCMatchExps.addAll(VariantContextUtils.initializeMatchExps(dummyNames, jexlExprStrings));
        }
        /* samples */
        final Set<String> samples0 = new HashSet<>(vcfIterator1.header.getSampleNamesInOrder());
        final Set<String> samples1 = new HashSet<>(vcfIterator2.header.getSampleNamesInOrder());
        final Set<String> commonSamples = new TreeSet<>(samples0);
        commonSamples.retainAll(samples1);
        if (commonSamples.size() != samples0.size() || commonSamples.size() != samples1.size()) {
            LOG.warn("Warning: Not the same samples set. Using intersection of both lists.");
        }
        if (commonSamples.isEmpty()) {
            LOG.error("No common samples");
            return -1;
        }
        final Map<String, SampleInfo> sample2info = new HashMap<>(commonSamples.size());
        for (final String sampleName : commonSamples) {
            sample2info.put(sampleName, new SampleInfo(sampleName));
        }
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(this.global_dictionary);
        for (; ; ) {
            VariantContext ctx0 = vcfIterator1.peek();
            VariantContext ctx1 = vcfIterator2.peek();
            final VariantContext smallest;
            if (ctx0 == null && ctx1 == null) {
                smallest = null;
                break;
            } else if (ctx0 == null && ctx1 != null) {
                smallest = ctx1;
            } else if (ctx0 != null && ctx1 == null) {
                smallest = ctx0;
            } else {
                final int diff = this.compareChromPosRef.compare(ctx0, ctx1);
                if (diff < 0) {
                    smallest = ctx0;
                    ctx1 = null;
                } else if (diff > 0) {
                    smallest = ctx1;
                    ctx0 = null;
                } else {
                    smallest = ctx0;
                }
            }
            progress.watch(smallest);
            vcfIterator1.reset(smallest);
            vcfIterator2.reset(smallest);
            if (capture != null) {
                final Interval interval = new Interval(smallest.getContig(), smallest.getStart(), smallest.getEnd());
                if (!capture.containsOverlapping(interval))
                    continue;
            }
            for (final String sampleName : sample2info.keySet()) {
                sample2info.get(sampleName).visit(ctx0, ctx1);
            }
        }
        progress.finish();
        vcfIterator1.close();
        vcfIterator2.close();
        archiveFactory = ArchiveFactory.open(this.archiveFile);
        makefileWriter = archiveFactory.openWriter(this.archivePrefix + "Makefile");
        makefileWriter.println(".PHONY:all all2");
        makefileWriter.println("ALL_TARGETS=");
        makefileWriter.println("all:all2");
        for (final SampleInfo sampleInfo : sample2info.values()) {
            for (final SampleCategory sampleCat : sampleInfo.sampleCat.values()) {
                final String basename = this.archivePrefix + sampleInfo.sampleName + "." + sampleCat.variantCatName.replace(' ', '_');
                final String tsv = basename + ".tsv";
                PrintWriter dataW = archiveFactory.openWriter(tsv);
                for (final String sk : new TreeSet<String>(sampleCat.counter.keySet())) {
                    dataW.print(escapeUnderscore(sk));
                    dataW.print("\t");
                    dataW.print(sampleCat.counter.count(sk));
                    dataW.println();
                }
                dataW.flush();
                dataW.close();
                final String png = "$(addsuffix .png," + basename + ")";
                makefileWriter.println("ALL_TARGETS+=" + png);
                makefileWriter.println(png + ":" + tsv);
                makefileWriter.println("\techo '" + "set ylabel \"Number of Genotypes  " + escapeUnderscore(sampleInfo.sampleName) + "\";" + "set yrange [0:];" + "set xlabel \"Category " + escapeUnderscore(vcf1Name) + ": " + vcfIterator1.count + ", " + escapeUnderscore(vcf2Name) + ": " + vcfIterator2.count + " variants  \";" + "set xtic rotate by 90 right;" + "set size  ratio 0.618;" + "set title \"" + escapeUnderscore(vcf1Name) + " vs " + escapeUnderscore(vcf2Name) + " : Genotypes " + escapeUnderscore(sampleInfo.sampleName) + " / Variants: " + escapeUnderscore(sampleCat.variantCatName) + " \";" + "set key  off;" + "set datafile separator \"\t\";" + "set auto x;" + "set style histogram;" + "set style data histogram;" + "set style fill solid border -1;" + "set terminal png truecolor  size " + (500 + sampleCat.counter.getCountCategories() * 50) + ", 1500;" + "set output \"$@\";" + "plot \"$<\" using 2:xtic(1) ti \"\";' | " + "gnuplot");
            }
        }
        makefileWriter.println("all2:${ALL_TARGETS}");
        makefileWriter.close();
        makefileWriter = null;
        archiveFactory.close();
        archiveFactory = null;
        return RETURN_OK;
    } catch (Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(makefileWriter);
        CloserUtil.close(archiveFactory);
        CloserUtil.close(vcfIterator1);
        CloserUtil.close(vcfIterator1);
    }
}

Example 4 with ArchiveFactory

use of com.github.lindenb.jvarkit.io.ArchiveFactory in project jvarkit by lindenb.

the class LumpyVcfToCircos method doWork.

@Override
public int doWork(List<String> args) {
    VcfIterator r = null;
    ArchiveFactory archiveFactory = null;
    PrintWriter pw = null;
    PrintWriter conf = null;
    try {
        r = super.openVcfIterator(oneFileOrNull(args));
        final VCFHeader header = r.getHeader();
        final SAMSequenceDictionary dict = header.getSequenceDictionary();
        if (dict == null) {
            LOG.error("No dictionary in vcf");
            return -1;
        }
        archiveFactory = ArchiveFactory.open(outputFile);
        conf = archiveFactory.openWriter(createFilename("ideogram.conf"));
        conf.println("<ideogram>");
        conf.println("show = yes");
        conf.println("<spacing>");
        conf.println("default = 1u");
        conf.println("<pairwise \"/hs/ /hs/\">");
        conf.println("spacing = 0.5r");
        conf.println("</pairwise>");
        conf.println("");
        conf.println("</spacing>");
        conf.println("");
        conf.println("thickness         = 25p");
        conf.println("");
        conf.println("#stroke_thickness = 1");
        conf.println("#stroke_color     = black");
        conf.println("");
        conf.println("fill           = yes");
        conf.println("fill_color     = black");
        conf.println("");
        conf.println("radius         = 1r");
        conf.println("show_label     = yes");
        conf.println("label_font     = default");
        conf.println("label_radius   = dims(ideogram,radius_outer) + 10p");
        conf.println("label_size     = 12p");
        conf.println("label_parallel = yes");
        conf.println("");
        conf.println("show_bands            = yes");
        conf.println("fill_bands            = yes");
        conf.println("band_stroke_thickness = 0");
        conf.println("band_stroke_color     = black");
        conf.println("band_transparency     = 4");
        conf.println("");
        conf.println("</ideogram>");
        conf.close();
        conf = archiveFactory.openWriter(createFilename("ticks.conf"));
        conf.println("");
        conf.println("show_ticks        = no");
        conf.println("show_tick_labels  = no");
        conf.println("show_grid         = no");
        conf.println("");
        conf.println("<ticks>");
        conf.println("tick_label_font  = light");
        conf.println("radius           = dims(ideogram,radius_outer) + 180p");
        conf.println("label_offset     = 5p");
        conf.println("label_size       = 16p");
        conf.println("multiplier       = 1e-6");
        conf.println("color            = black");
        conf.println("thickness        = 1p");
        conf.println("");
        conf.println("# 25 Mb ticks, all chromosomes, with labels");
        conf.println("");
        conf.println("<tick>");
        conf.println("spacing        = 25u");
        conf.println("size           = 12p");
        conf.println("show_label     = yes");
        conf.println("format         = %d");
        conf.println("</tick>");
        conf.println("");
        conf.println("# 5 Mb ticks, all chromosomes, with labels");
        conf.println("# labels must be separated by at least 1px, which");
        conf.println("#   avoids overlap on human chrs");
        conf.println("<tick>");
        conf.println("label_separation = 1p");
        conf.println("spacing          = 5u");
        conf.println("size             = 7p");
        conf.println("show_label       = yes");
        conf.println("format           = %d");
        conf.println("</tick>");
        conf.println("");
        conf.println("# 5 Mb ticks with grid on rn1 and mm1, drawn only");
        conf.println("# for their grid (note size=0p)");
        conf.println("<tick>");
        conf.println("chromosomes_display_default = no");
        conf.println("chromosomes    = rn1;mm1");
        conf.println("spacing        = 5u");
        conf.println("size           = 0p");
        conf.println("force_display  = yes");
        conf.println("grid_start     = 0.45r");
        conf.println("grid_end       = dims(ideogram,radius_outer) + 180p");
        conf.println("grid_color     = grey");
        conf.println("grid_thickness = 1p");
        conf.println("grid           = yes");
        conf.println("</tick>");
        conf.println("");
        conf.println("# 20% relative ticks on human chromosomes");
        conf.println("");
        conf.println("<tick>");
        conf.println("chromosomes    = -rn1;-mm1");
        conf.println("radius         = 0.95r");
        conf.println("spacing_type   = relative");
        conf.println("rspacing       = 0.20");
        conf.println("size           = 6p");
        conf.println("show_label     = yes");
        conf.println("label_relative = yes");
        conf.println("rmultiplier    = 100");
        conf.println("format         = %d");
        conf.println("suffix         = %");
        conf.println("");
        conf.println("skip_last_label= yes");
        conf.println("");
        conf.println("grid_start     = 0.885r");
        conf.println("grid_end       = 0.95r");
        conf.println("grid_color     = grey");
        conf.println("grid_thickness = 1p");
        conf.println("grid           = yes");
        conf.println("</tick>");
        conf.println("");
        conf.println("# relative ticks at start and end of human chromosomes");
        conf.println("# with grid");
        conf.println("");
        conf.println("<tick>");
        conf.println("spacing_type   = relative");
        conf.println("rspacing       = 1");
        conf.println("size           = 6p");
        conf.println("grid_start     = 0.45r");
        conf.println("grid_end       = dims(ideogram,radius_outer) + 180p");
        conf.println("grid_color     = grey");
        conf.println("grid_thickness = 1p");
        conf.println("grid           = yes");
        conf.println("</tick>");
        conf.println("");
        conf.println("<tick>");
        conf.println("chromosomes    = -rn1;-mm1");
        conf.println("radius         = 0.95r");
        conf.println("spacing_type   = relative");
        conf.println("rspacing       = 0.10");
        conf.println("size           = 3p");
        conf.println("show_label     = no");
        conf.println("");
        conf.println("grid_start     = 0.885r");
        conf.println("grid_end       = 0.95r");
        conf.println("grid_color     = lgrey");
        conf.println("grid_thickness = 1p");
        conf.println("grid           = yes");
        conf.println("</tick>");
        conf.println("");
        conf.println("# 25% relative ticks on human chromosomes");
        conf.println("");
        conf.println("<tick>");
        conf.println("chromosomes    = -rn1;-mm1");
        conf.println("radius         = 0.82r");
        conf.println("spacing_type   = relative");
        conf.println("rspacing       = 0.25");
        conf.println("size           = 6p");
        conf.println("show_label     = yes");
        conf.println("label_relative = yes");
        conf.println("rmultiplier    = 100");
        conf.println("format         = %d");
        conf.println("");
        conf.println("skip_last_label= yes");
        conf.println("");
        conf.println("grid_start     = 0.755r");
        conf.println("grid_end       = 0.82r");
        conf.println("grid_color     = grey");
        conf.println("grid_thickness = 1p");
        conf.println("grid           = yes");
        conf.println("</tick>");
        conf.println("");
        conf.println("</ticks>");
        conf.close();
        conf = archiveFactory.openWriter(createFilename("circos.conf"));
        conf.println("karyotype=" + createFilename("karyotype.txt"));
        conf.println("chromosomes_units=1000000");
        conf.println("<<include " + createFilename("ideogram.conf") + ">>");
        conf.println("<<include " + createFilename("ticks.conf") + ">>");
        final Map<String, SampleWriters> sample2writers = new HashMap<>(header.getNGenotypeSamples());
        for (final String sampleName : header.getSampleNamesInOrder()) {
            sample2writers.put(sampleName, new SampleWriters(sampleName));
        }
        while (r.hasNext()) {
            final VariantContext vc = r.next();
            for (final SampleWriters sw : sample2writers.values()) {
                sw.visit(vc);
            }
        }
        r.close();
        r = null;
        LOG.info("writing conf");
        final Set<String> seen_chromosomes = new HashSet<>(dict.size());
        for (final SampleWriters sw : sample2writers.values()) {
            seen_chromosomes.addAll(sw.bdns.stream().map(SV -> SV.contig).collect(Collectors.toSet()));
            seen_chromosomes.addAll(sw.bdns.stream().map(SV -> SV.bnb_contig).collect(Collectors.toSet()));
            seen_chromosomes.addAll(sw.duplications.stream().map(SV -> SV.contig).collect(Collectors.toSet()));
            seen_chromosomes.addAll(sw.invertions.stream().map(SV -> SV.contig).collect(Collectors.toSet()));
            seen_chromosomes.addAll(sw.deletions.stream().map(SV -> SV.contig).collect(Collectors.toSet()));
        }
        pw = archiveFactory.openWriter(createFilename("karyotype.txt"));
        for (int i = 0; i < dict.size(); ++i) {
            final SAMSequenceRecord ssr = dict.getSequence(i);
            if (!seen_chromosomes.contains(ssr.getSequenceName()))
                continue;
            pw.print("chr - ");
            pw.print(ssr.getSequenceName());
            pw.print(" ");
            pw.print(ssr.getSequenceName());
            pw.print(" 0 ");
            pw.print(ssr.getSequenceLength());
            pw.print(" ");
            pw.println(i % 2 == 0 ? "blue" : "grey");
        }
        pw.flush();
        pw.close();
        Function<String, String> sample2color = S -> {
            int i = header.getSampleNamesInOrder().indexOf(S);
            if (i != -1) {
                if (i % 2 == 0)
                    i = ((header.getSampleNamesInOrder().size() - 1) - i);
                int g = 30 + (int) ((i / (double) header.getSampleNamesInOrder().size()) * 50.0);
                return "gray" + (int) g;
            }
            return (i % 2 == 0 ? "green" : "red");
        };
        final double minR = 0.1;
        final double maxR = 0.95;
        final double numberOfTracks = sample2writers.values().stream().mapToInt(X -> X.getNumberOfTracks()).sum();
        // final double maxSu = sample2writers.values().stream().mapToInt(X->X.getMaxSupportingEvidences()).max().orElse(1);
        double radius = minR;
        final double dr = (maxR - minR) / numberOfTracks;
        LOG.info("N tracks = " + numberOfTracks + " dr=" + dr);
        conf.println("<plots>");
        int maxSu = sample2writers.values().stream().flatMap(S -> S.deletions.stream()).mapToInt(SV -> SV.su).max().orElse(0);
        for (final SampleWriters sw : sample2writers.values()) {
            if (maxSu == 0 || sw.deletions.isEmpty())
                continue;
            String file = sw.getDeletionFilename();
            pw = archiveFactory.openWriter(file);
            for (final SVDel sv : sw.deletions) sv.print(pw);
            pw.flush();
            pw.close();
            conf.println("<plot>");
            conf.println("type = histogram");
            conf.println("file = " + file);
            conf.println("min = 0");
            conf.println("max = " + maxSu);
            conf.println("orientation = out");
            conf.println("layers = 1");
            conf.println("margin = 0.02u");
            conf.println("color = black_a4");
            conf.println("r0 = " + radius + "r");
            conf.println("r1 = " + (radius + dr) + "r");
            conf.println("fill_color = " + sample2color.apply(sw.sampleName));
            // +sample2color.apply(sw.sampleName));
            conf.println("stroke_color = grey");
            conf.println("thickness = 1");
            conf.println("extend_bin  = no");
            // conf.println("<backgrounds>\n<background>\ncolor ="+(zz%2==0?"pink":"yellow")+"\n</background>\n</backgrounds>");
            // conf.println("<axes>\n<axis>\nspacing   = "+(dr/10.0)+"\ncolor=lgrey\nthickness=2\n</axis>\n</axes>");
            conf.println("</plot>");
            radius += dr;
        }
        maxSu = sample2writers.values().stream().flatMap(S -> S.invertions.stream()).mapToInt(SV -> SV.su).max().orElse(0);
        for (final SampleWriters sw : sample2writers.values()) {
            if (maxSu == 0 || sw.invertions.isEmpty())
                continue;
            String file = sw.getInvertionFilename();
            pw = archiveFactory.openWriter(file);
            for (final SVInv sv : sw.invertions) sv.print(pw);
            pw.flush();
            pw.close();
            conf.println("<plot>");
            conf.println("type = histogram");
            conf.println("file = " + file);
            conf.println("min = 0");
            conf.println("max = " + maxSu);
            conf.println("color = black_a4");
            conf.println("r0 = " + radius + "r");
            conf.println("r1 = " + (radius + dr) + "r");
            conf.println("fill_color = " + sample2color.apply(sw.sampleName));
            conf.println("stroke_color = " + sample2color.apply(sw.sampleName));
            conf.println("thickness = 1");
            conf.println("extend_bin  = no");
            conf.println("</plot>");
            radius += dr;
        }
        maxSu = sample2writers.values().stream().flatMap(S -> S.duplications.stream()).mapToInt(SV -> SV.su).max().orElse(0);
        for (final SampleWriters sw : sample2writers.values()) {
            if (maxSu == 0 || sw.duplications.isEmpty())
                continue;
            String file = sw.getDuplicationFilename();
            pw = archiveFactory.openWriter(file);
            for (final SVDup sv : sw.duplications) sv.print(pw);
            pw.flush();
            pw.close();
            conf.println("<plot>");
            conf.println("type = histogram");
            conf.println("file = " + file);
            conf.println("min = 0");
            conf.println("max = " + maxSu);
            conf.println("color = black_a4");
            conf.println("r0 = " + radius + "r");
            conf.println("r1 = " + (radius + dr) + "r");
            conf.println("fill_color = " + sample2color.apply(sw.sampleName));
            conf.println("stroke_color = " + sample2color.apply(sw.sampleName));
            conf.println("thickness = 1");
            conf.println("extend_bin  = no");
            conf.println("</plot>");
            radius += dr;
        }
        conf.println("</plots>");
        conf.println("<links>");
        for (final SampleWriters sw : sample2writers.values()) {
            if (sw.bdns.isEmpty())
                continue;
            pw = archiveFactory.openWriter(sw.getBnbFilename());
            for (final SVBnB sv : sw.bdns) sv.print(pw);
            pw.flush();
            pw.close();
            conf.println("<link>");
            conf.println("file=" + sw.getBnbFilename());
            conf.println("radius=" + radius + "r");
            conf.println("bezier_radius = " + (radius + dr) + "r");
            conf.println("fill_color = " + sample2color.apply(sw.sampleName));
            conf.println("stroke_color = " + sample2color.apply(sw.sampleName));
            conf.println("</link>");
            radius += dr;
        }
        conf.println("</links>");
        sample2writers.clear();
        conf.println("");
        conf.println("<image>");
        conf.println("# Included from Circos distribution.");
        conf.println("<<include etc/image.conf>>");
        conf.println("</image>");
        conf.println("");
        conf.println("<<include etc/colors_fonts_patterns.conf>>");
        conf.println("");
        conf.println("# Debugging, I/O an dother system parameters");
        conf.println("# Included from Circos distribution.");
        conf.println("<<include etc/housekeeping.conf>>");
        conf.println("");
        conf.flush();
        conf.close();
        conf = null;
        archiveFactory.close();
        archiveFactory = null;
        LOG.info("done");
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(conf);
        CloserUtil.close(archiveFactory);
        CloserUtil.close(r);
    }
}

Example 5 with ArchiveFactory

use of com.github.lindenb.jvarkit.io.ArchiveFactory in project jvarkit by lindenb.

the class VcfStats method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.binSize <= 0) {
        LOG.error("binSize < 0");
        return -1;
    }
    VariantContextWriter teeOut = null;
    VcfIterator iter = null;
    final Map<String, VariantStats> category2stats = new HashMap<>();
    PrintWriter makefileWriter = null;
    try {
        this.archiveFactory = ArchiveFactory.open(this.outputFile);
        if (this.tee)
            teeOut = super.openVariantContextWriter(null);
        iter = super.openVcfIterator(oneFileOrNull(args));
        final VCFHeader header = iter.getHeader();
        this.sampleNamesInOrder = Collections.unmodifiableList(header.getSampleNamesInOrder());
        final SAMSequenceDictionary dict = header.getSequenceDictionary();
        if (dict != null && !dict.isEmpty()) {
            this.the_dictionary = dict;
        }
        if (this.kgFile != null) {
            LOG.info("load " + kgFile);
            this.knownGeneTreeMap = KnownGene.loadUriAsIntervalTreeMap(this.kgFile, KG -> (dict == null || dict.getSequence(KG.getContig()) != null));
        } else {
            this.knownGeneTreeMap = null;
        }
        if (this.pedigreeFile != null) {
            this.pedigree = Pedigree.newParser().parse(this.pedigreeFile);
        } else {
            Pedigree tmpPed = null;
            try {
                tmpPed = Pedigree.newParser().parse(header);
            } catch (Exception err) {
                tmpPed = Pedigree.createEmptyPedigree();
            }
            this.pedigree = tmpPed;
        }
        makefileWriter = this.archiveFactory.openWriter(this.prefix + "Makefile");
        makefileWriter.println(".PHONY: all all_targets ");
        makefileWriter.println("SCREEN_WIDTH?=2600");
        makefileWriter.println("SCREEN_HEIGHT?=1000");
        makefileWriter.println("ALL_TARGETS=");
        makefileWriter.println("all: all_targets");
        if (teeOut != null)
            teeOut.writeHeader(header);
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(header).logger(LOG);
        while (iter.hasNext()) {
            final VariantContext ctx = progress.watch(iter.next());
            if (teeOut != null)
                teeOut.add(ctx);
            for (final String category : this.variantToCategoryKeys.apply(ctx)) {
                VariantStats vcstat = category2stats.get(category);
                if (vcstat == null) {
                    vcstat = new VariantStats(category, header);
                    category2stats.put(category, vcstat);
                }
                vcstat.visit(ctx);
            }
        }
        for (final String category : category2stats.keySet()) {
            final VariantStats vcstats = category2stats.get(category);
            vcstats.finish(makefileWriter);
        }
        progress.finish();
        makefileWriter.println("all_targets : ${ALL_TARGETS}");
        makefileWriter.flush();
        makefileWriter.close();
        makefileWriter = null;
        iter.close();
        iter = null;
        this.archiveFactory.close();
        archiveFactory = null;
        if (teeOut != null)
            teeOut.close();
        teeOut = null;
        return 0;
    } catch (Exception e) {
        LOG.error(e);
        return -1;
    } finally {
        knownGeneTreeMap = null;
        CloserUtil.close(archiveFactory);
        CloserUtil.close(teeOut);
        CloserUtil.close(iter);
        CloserUtil.close(makefileWriter);
    }
}