Search in sources :

Example 21 with CharSplitter

use of com.github.lindenb.jvarkit.lang.CharSplitter in project jvarkit by lindenb.

the class PrettyTable method doWork.

@Override
public int doWork(final List<String> args) {
    final List<List<String>> rows = new ArrayList<>();
    PrintWriter out = null;
    BufferedReader br = null;
    try {
        final CharSplitter delim = CharSplitter.of(this.delimiter);
        int n_columns = -1;
        br = super.openBufferedReader(oneFileOrNull(args));
        String line;
        while ((line = br.readLine()) != null) {
            if (StringUtils.isBlank(line))
                continue;
            final List<String> row = delim.splitAsStringList(line);
            if (rows.isEmpty()) {
                n_columns = row.size();
                LOG.info("n=" + n_columns);
                if (this.noHeader) {
                    final List<String> header = new ArrayList<>(n_columns);
                    while (header.size() < n_columns) {
                        header.add("$" + (1 + header.size()));
                    }
                    rows.add(header);
                }
            } else if (row.size() < n_columns) {
                while (row.size() < n_columns) row.add("");
            } else if (row.size() > n_columns) {
                n_columns = row.size();
                for (int y = 0; y < rows.size(); ++y) {
                    while (rows.get(y).size() < n_columns) rows.get(y).add("");
                }
            } else {
            // ok
            }
            System.err.println(">>>" + row.get(0));
            rows.add(row);
        }
        br.close();
        br = null;
        if (rows.isEmpty()) {
            LOG.error("empty input");
            return -1;
        }
        if (this.transpose) {
            final List<List<String>> rows2 = new ArrayList<>(n_columns);
            for (int x = 0; x < n_columns; ++x) {
                final List<String> row = new ArrayList<>(rows.size());
                for (int y = 0; y < rows.size(); ++y) {
                    row.add(rows.get(y).get(x));
                }
                rows2.add(row);
            }
            rows.clear();
            rows.addAll(rows2);
        }
        final TableFactory tableFactory = new TableFactory();
        final Table table = tableFactory.createTable(rows.get(0));
        for (int y = 1; y < rows.size(); ++y) {
            table.addRow(rows.get(y));
        }
        out = super.openFileOrStdoutAsPrintWriter(this.outputFile);
        if (this.format.equalsIgnoreCase("html")) {
            final HtmlExporter exporter = new HtmlExporter();
            exporter.setRemoveEmptyColumns(!this.keepEmptyComlumns);
            exporter.saveTableTo(table, out);
        } else {
            final TextExporter exporter = new TextExporter();
            exporter.setAllowUnicode(this.withUnicode);
            exporter.setRemoveEmptyColumns(!this.keepEmptyComlumns);
            exporter.saveTableTo(table, out);
        }
        out.flush();
        out.close();
        out = null;
        return 0;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(out);
    }
}
Also used : Table(com.github.lindenb.jvarkit.table.Table) CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter) HtmlExporter(com.github.lindenb.jvarkit.table.HtmlExporter) ArrayList(java.util.ArrayList) TextExporter(com.github.lindenb.jvarkit.table.TextExporter) TableFactory(com.github.lindenb.jvarkit.table.TableFactory) BufferedReader(java.io.BufferedReader) ArrayList(java.util.ArrayList) List(java.util.List) PrintWriter(java.io.PrintWriter)

Example 22 with CharSplitter

use of com.github.lindenb.jvarkit.lang.CharSplitter in project jvarkit by lindenb.

the class IndexCovToSV method doWork.

@Override
public int doWork(final List<String> args) {
    final IndexCovUtils indexCovUtils = new IndexCovUtils(this.treshold);
    final CharSplitter tab = CharSplitter.TAB;
    try {
        final SAMSequenceDictionary dict = SequenceDictionaryUtils.extractRequired(this.refFile);
        final List<String> samples;
        SortingCollection<Cell> sorter1 = SortingCollection.newInstance(Cell.class, new CellCodec(), (A, B) -> A.sortSamplePosition(B), this.sortingDelegate.getMaxRecordsInRam(), this.sortingDelegate.getTmpPaths());
        sorter1.setDestructiveIteration(true);
        try (BufferedReader r = super.openBufferedReader(oneFileOrNull(args))) {
            String line = r.readLine();
            if (line == null) {
                LOG.error("Cannot read first line of input");
                return -1;
            }
            String[] tokens = tab.split(line);
            if (tokens.length < 4 || !tokens[0].equals("#chrom") || !tokens[1].equals("start") || !tokens[2].equals("end")) {
                LOG.error("bad first line " + line);
                return -1;
            }
            samples = Arrays.asList(tokens).subList(3, tokens.length);
            while ((line = r.readLine()) != null) {
                if (StringUtil.isBlank(line))
                    continue;
                tokens = tab.split(line);
                if (tokens.length != 3 + samples.size()) {
                    throw new JvarkitException.TokenErrors(samples.size() + 3, tokens);
                }
                final SAMSequenceRecord ssr = dict.getSequence(tokens[0]);
                if (ssr == null) {
                    LOG.error(JvarkitException.ContigNotFoundInDictionary.getMessage(tokens[0], dict));
                    return -1;
                }
                for (int i = 3; i < tokens.length; i++) {
                    final double score = Double.parseDouble(tokens[i]);
                    if (indexCovUtils.getType(score).isAmbigous())
                        continue;
                    final Cell cell = new Cell();
                    cell.tid = ssr.getSequenceIndex();
                    cell.start = Integer.parseInt(tokens[1]);
                    cell.end = Integer.parseInt(tokens[2]);
                    cell.sample_id = i - 3;
                    cell.score = score;
                    sorter1.add(cell);
                }
            }
            // while readline
            sorter1.doneAdding();
        }
        /* merge adjacent blocks */
        SortingCollection<Cell> sorter2 = SortingCollection.newInstance(Cell.class, new CellCodec(), (A, B) -> A.sortPositionSample(B), this.sortingDelegate.getMaxRecordsInRam(), this.sortingDelegate.getTmpPaths());
        sorter2.setDestructiveIteration(true);
        try (CloseableIterator<Cell> iter1 = sorter1.iterator()) {
            try (PeekableIterator<Cell> peeker1 = new PeekableIterator<>(iter1)) {
                while (peeker1.hasNext()) {
                    final Cell first = peeker1.next();
                    while (peeker1.hasNext()) {
                        final Cell second = peeker1.peek();
                        if (first.sample_id != second.sample_id)
                            break;
                        if (first.tid != second.tid)
                            break;
                        if (first.end + (this.merge_distance + 1) < second.start)
                            break;
                        if (!indexCovUtils.getType(first.score).equals(indexCovUtils.getType(second.score)))
                            break;
                        // extend first with end of second
                        first.end = second.end;
                        // consumme
                        peeker1.next();
                    }
                    sorter2.add(first);
                }
            // while peeker1
            }
        }
        sorter1.cleanup();
        sorter1 = null;
        sorter2.doneAdding();
        final Set<VCFHeaderLine> metaData = new HashSet<>();
        VCFStandardHeaderLines.addStandardFormatLines(metaData, true, VCFConstants.GENOTYPE_KEY, VCFConstants.GENOTYPE_QUALITY_KEY);
        VCFStandardHeaderLines.addStandardInfoLines(metaData, true, VCFConstants.END_KEY, VCFConstants.ALLELE_FREQUENCY_KEY, VCFConstants.ALLELE_COUNT_KEY, VCFConstants.ALLELE_NUMBER_KEY);
        final VCFFilterHeaderLine filterHomDel = new VCFFilterHeaderLine("HOM_DEL", "There is one Homozygous deletion.");
        metaData.add(filterHomDel);
        final VCFFilterHeaderLine filterHomDup = new VCFFilterHeaderLine("HOM_DUP", "There is one Homozygous duplication.");
        metaData.add(filterHomDup);
        /**
         * raw value in indexcov
         */
        final VCFFormatHeaderLine foldHeader = new VCFFormatHeaderLine("F", 1, VCFHeaderLineType.Float, "Relative number of copy: 0.5 deletion 1 normal 2.0 duplication");
        metaData.add(foldHeader);
        final VCFInfoHeaderLine infoSvType = new VCFInfoHeaderLine(VCFConstants.SVTYPE, 1, VCFHeaderLineType.String, "SV type");
        metaData.add(infoSvType);
        final VCFInfoHeaderLine infoSvLen = new VCFInfoHeaderLine("SVLEN", 1, VCFHeaderLineType.Integer, "ABS(SV length)");
        metaData.add(infoSvLen);
        /*
				final VCFInfoHeaderLine infoSvLens = new VCFInfoHeaderLine("SVLENS", 1, VCFHeaderLineType.Integer,"Allele length");
				metaData.add(infoSvLens);
				final VCFInfoHeaderLine infoSvEnds = new VCFInfoHeaderLine("SVENDS", 1, VCFHeaderLineType.Integer,"Allele ends");
				metaData.add(infoSvEnds);
				 */
        final VCFHeader vcfHeader = new VCFHeader(metaData, samples);
        JVarkitVersion.getInstance().addMetaData(this, vcfHeader);
        vcfHeader.setSequenceDictionary(dict);
        try (VariantContextWriter vcw = this.writingDelegate.dictionary(dict).open(outputFile)) {
            vcw.writeHeader(vcfHeader);
            final Allele REF_ALLELE = Allele.create("N", true);
            try (CloseableIterator<Cell> iter2 = sorter2.iterator()) {
                try (EqualRangeIterator<Cell> peeker2 = new EqualRangeIterator<>(iter2, (A, B) -> A.sortChromStart(B))) {
                    while (peeker2.hasNext()) {
                        final List<Cell> array0 = peeker2.next();
                        final Set<Integer> distinct_ends = array0.stream().map(C -> C.end).collect(Collectors.toCollection(TreeSet::new));
                        final List<List<Cell>> array_of_array;
                        if (this.normalize_one_allele && distinct_ends.size() > 1) {
                            array_of_array = new ArrayList<>(distinct_ends.size());
                            for (final int end_variant : distinct_ends) {
                                array_of_array.add(array0.stream().filter(C -> C.end == end_variant).collect(Collectors.toList()));
                            }
                        } else {
                            array_of_array = Collections.singletonList(array0);
                        }
                        for (final List<Cell> array : array_of_array) {
                            int an = 0;
                            if (array.stream().map(C -> indexCovUtils.getType(C.score)).noneMatch(C -> C.isVariant()))
                                continue;
                            final Cell first = array.get(0);
                            final int max_end = array.stream().mapToInt(C -> C.end).max().getAsInt();
                            final List<Genotype> genotypes = new ArrayList<>(array.size());
                            final VariantContextBuilder vcb = new VariantContextBuilder();
                            vcb.chr(dict.getSequence(first.tid).getSequenceName());
                            vcb.start(first.start);
                            vcb.stop(max_end);
                            vcb.attribute(VCFConstants.END_KEY, max_end);
                            final List<AltAllele> altAlleles = new ArrayList<>();
                            boolean has_hom_del = false;
                            boolean has_hom_dup = false;
                            for (int i = 0; i < array.size(); i++) {
                                final Cell cell = array.get(i);
                                final String sampleName = samples.get(cell.sample_id);
                                an += 2;
                                final GenotypeBuilder gb;
                                final IndexCovUtils.SvType type = indexCovUtils.getType(cell.score);
                                final Allele allele;
                                AltAllele altAllele;
                                if (type.isDeletion() || type.isDuplication()) {
                                    altAllele = altAlleles.stream().filter(C -> C.end == cell.end && C.type.equals(type)).findFirst().orElse(null);
                                    if (altAllele == null) {
                                        allele = Allele.create("<" + (type.isDeletion() ? "DEL" : "DUP") + "." + (1 + altAlleles.size()) + ">", false);
                                        altAllele = new AltAllele(type, cell.end, allele);
                                        altAlleles.add(altAllele);
                                    } else {
                                        allele = altAllele.allele;
                                    }
                                } else {
                                    altAllele = null;
                                    allele = null;
                                }
                                switch(type) {
                                    case REF:
                                        gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, REF_ALLELE));
                                        break;
                                    case HET_DEL:
                                        gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, allele));
                                        altAllele.ac++;
                                        break;
                                    case HOM_DEL:
                                        gb = new GenotypeBuilder(sampleName, Arrays.asList(allele, allele));
                                        vcb.filter(filterHomDel.getID());
                                        altAllele.ac += 2;
                                        has_hom_del = true;
                                        break;
                                    case HET_DUP:
                                        gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, allele));
                                        altAllele.ac++;
                                        break;
                                    case HOM_DUP:
                                        gb = new GenotypeBuilder(sampleName, Arrays.asList(allele, allele));
                                        vcb.filter(filterHomDup.getID());
                                        altAllele.ac += 2;
                                        has_hom_dup = true;
                                        break;
                                    default:
                                        gb = new GenotypeBuilder(sampleName, Arrays.asList(Allele.NO_CALL, Allele.NO_CALL));
                                        break;
                                }
                                gb.attribute(foldHeader.getID(), cell.score);
                                if (!type.isAmbigous()) {
                                    gb.GQ(type.getGenotypeQuality(cell.score));
                                }
                                genotypes.add(gb.make());
                            }
                            // end loop over genotypes
                            final List<Allele> alleles = altAlleles.stream().map(A -> A.allele).collect(Collectors.toCollection(ArrayList::new));
                            if (has_hom_del) {
                                vcb.filter(filterHomDel.getID());
                            } else if (has_hom_dup) {
                                vcb.filter(filterHomDup.getID());
                            } else {
                                vcb.passFilters();
                            }
                            alleles.add(0, REF_ALLELE);
                            vcb.alleles(alleles);
                            vcb.genotypes(genotypes);
                            if (altAlleles.stream().allMatch(C -> C.type.isDeletion())) {
                                vcb.attribute(VCFConstants.SVTYPE, "DEL");
                            } else if (altAlleles.stream().allMatch(C -> C.type.isDuplication())) {
                                vcb.attribute(VCFConstants.SVTYPE, "DUP");
                            } else {
                                vcb.attribute(VCFConstants.SVTYPE, "MIXED");
                            }
                            vcb.attribute(infoSvLen.getID(), altAlleles.stream().mapToInt(C -> 1 + C.end - first.start).max().orElse(0));
                            vcb.attribute(VCFConstants.ALLELE_COUNT_KEY, altAlleles.stream().map(A -> A.ac).collect(Collectors.toList()));
                            vcb.attribute(VCFConstants.ALLELE_NUMBER_KEY, an);
                            if (an > 0) {
                                final int final_an = an;
                                vcb.attribute(VCFConstants.ALLELE_FREQUENCY_KEY, altAlleles.stream().map(A -> A.ac / (double) final_an).collect(Collectors.toList()));
                            }
                            vcw.add(vcb.make());
                        }
                    }
                // end loop while iter
                }
            }
            sorter2.cleanup();
            sorter2 = null;
        }
        return 0;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
    }
}
Also used : WritingVariantsDelegate(com.github.lindenb.jvarkit.variant.variantcontext.writer.WritingVariantsDelegate) Genotype(htsjdk.variant.variantcontext.Genotype) VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) DataInputStream(java.io.DataInputStream) CloseableIterator(htsjdk.samtools.util.CloseableIterator) Allele(htsjdk.variant.variantcontext.Allele) Arrays(java.util.Arrays) CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter) Program(com.github.lindenb.jvarkit.util.jcommander.Program) Parameter(com.beust.jcommander.Parameter) VCFStandardHeaderLines(htsjdk.variant.vcf.VCFStandardHeaderLines) VCFHeader(htsjdk.variant.vcf.VCFHeader) TreeSet(java.util.TreeSet) ParametersDelegate(com.beust.jcommander.ParametersDelegate) ArrayList(java.util.ArrayList) HashSet(java.util.HashSet) DistanceParser(com.github.lindenb.jvarkit.util.bio.DistanceParser) NoSplitter(com.github.lindenb.jvarkit.util.jcommander.NoSplitter) DataOutputStream(java.io.DataOutputStream) StringUtil(htsjdk.samtools.util.StringUtil) AbstractDataCodec(com.github.lindenb.jvarkit.util.picard.AbstractDataCodec) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) Path(java.nio.file.Path) PeekableIterator(htsjdk.samtools.util.PeekableIterator) VCFConstants(htsjdk.variant.vcf.VCFConstants) SequenceDictionaryUtils(com.github.lindenb.jvarkit.util.bio.SequenceDictionaryUtils) SortingCollection(htsjdk.samtools.util.SortingCollection) VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine) VCFHeaderLineType(htsjdk.variant.vcf.VCFHeaderLineType) GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder) Logger(com.github.lindenb.jvarkit.util.log.Logger) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) Set(java.util.Set) IOException(java.io.IOException) JVarkitVersion(com.github.lindenb.jvarkit.util.JVarkitVersion) EOFException(java.io.EOFException) Collectors(java.util.stream.Collectors) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) List(java.util.List) FractionConverter(com.github.lindenb.jvarkit.jcommander.converter.FractionConverter) EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine) VCFFormatHeaderLine(htsjdk.variant.vcf.VCFFormatHeaderLine) BufferedReader(java.io.BufferedReader) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) Collections(java.util.Collections) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) ArrayList(java.util.ArrayList) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator) ArrayList(java.util.ArrayList) List(java.util.List) VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VCFHeader(htsjdk.variant.vcf.VCFHeader) HashSet(java.util.HashSet) VCFFormatHeaderLine(htsjdk.variant.vcf.VCFFormatHeaderLine) CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter) Genotype(htsjdk.variant.variantcontext.Genotype) GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder) VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine) Allele(htsjdk.variant.variantcontext.Allele) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) BufferedReader(java.io.BufferedReader) PeekableIterator(htsjdk.samtools.util.PeekableIterator)

Example 23 with CharSplitter

use of com.github.lindenb.jvarkit.lang.CharSplitter in project jvarkit by lindenb.

the class IndexCovToVcf method doWork.

@Override
public int doWork(final List<String> args) {
    final IndexCovUtils indexCovUtils = new IndexCovUtils(this.treshold);
    final CharSplitter tab = CharSplitter.TAB;
    BufferedReader r = null;
    VariantContextWriter vcw = null;
    try {
        final SAMSequenceDictionary dict;
        if (this.refFile == null) {
            dict = null;
        } else {
            dict = SequenceDictionaryUtils.extractRequired(this.refFile);
        }
        r = super.openBufferedReader(oneFileOrNull(args));
        String line = r.readLine();
        if (line == null) {
            LOG.error("Cannot read first line of input");
            return -1;
        }
        String[] tokens = tab.split(line);
        if (tokens.length < 4 || !tokens[0].equals("#chrom") || !tokens[1].equals("start") || !tokens[2].equals("end")) {
            LOG.error("bad first line " + line);
            return -1;
        }
        final Set<VCFHeaderLine> metaData = new HashSet<>();
        VCFStandardHeaderLines.addStandardFormatLines(metaData, true, VCFConstants.GENOTYPE_KEY, VCFConstants.GENOTYPE_QUALITY_KEY);
        VCFStandardHeaderLines.addStandardInfoLines(metaData, true, VCFConstants.END_KEY);
        /**
         * raw value in indexcov
         */
        final VCFFormatHeaderLine foldHeader = new VCFFormatHeaderLine("F", 1, VCFHeaderLineType.Float, "Relative number of copy: 0.5 deletion 1 normal 2.0 duplication");
        metaData.add(foldHeader);
        final VCFFilterHeaderLine filterAllDel = new VCFFilterHeaderLine("ALL_DEL", "number of samples greater than 1 and all are deletions");
        metaData.add(filterAllDel);
        final VCFFilterHeaderLine filterAllDup = new VCFFilterHeaderLine("ALL_DUP", "number of samples  greater than  1 and all are duplication");
        metaData.add(filterAllDup);
        final VCFFilterHeaderLine filterNoSV = new VCFFilterHeaderLine("NO_SV", "There is no DUP or DEL in this variant");
        metaData.add(filterNoSV);
        final VCFFilterHeaderLine filterHomDel = new VCFFilterHeaderLine("HOM_DEL", "There is one Homozygous deletion.");
        metaData.add(filterHomDel);
        final VCFFilterHeaderLine filterHomDup = new VCFFilterHeaderLine("HOM_DUP", "There is one Homozygous duplication.");
        metaData.add(filterHomDup);
        final VCFInfoHeaderLine infoNumDup = new VCFInfoHeaderLine("NDUP", 1, VCFHeaderLineType.Integer, "Number of samples being duplicated");
        metaData.add(infoNumDup);
        final VCFInfoHeaderLine infoNumDel = new VCFInfoHeaderLine("NDEL", 1, VCFHeaderLineType.Integer, "Number of samples being deleted");
        metaData.add(infoNumDel);
        final VCFInfoHeaderLine infoSingleton = new VCFInfoHeaderLine("SINGLETON", 1, VCFHeaderLineType.Flag, "Singleton candidate");
        metaData.add(infoSingleton);
        final VCFInfoHeaderLine infoAllAffected = new VCFInfoHeaderLine("ALL_CASES", 1, VCFHeaderLineType.Flag, "All cases are affected");
        metaData.add(infoAllAffected);
        final List<String> samples = Arrays.asList(tokens).subList(3, tokens.length);
        final Pedigree pedigree;
        final int count_cases_in_pedigree;
        if (this.pedFile == null) {
            pedigree = PedigreeParser.empty();
            count_cases_in_pedigree = 0;
        } else {
            pedigree = new PedigreeParser().parse(this.pedFile);
            final Set<String> set = new HashSet<>(samples);
            count_cases_in_pedigree = (int) pedigree.getAffectedSamples().stream().filter(S -> set.contains(S.getId())).count();
        }
        final VCFHeader vcfHeader = new VCFHeader(metaData, samples);
        JVarkitVersion.getInstance().addMetaData(this, vcfHeader);
        if (dict != null) {
            vcfHeader.setSequenceDictionary(dict);
        }
        vcw = this.writingDelegate.dictionary(dict).open(outputFile);
        vcw.writeHeader(vcfHeader);
        // final List<Allele> NO_CALL_NO_CALL = Arrays.asList(Allele.NO_CALL,Allele.NO_CALL);
        final Allele DUP_ALLELE = Allele.create("<DUP>", false);
        final Allele DEL_ALLELE = Allele.create("<DEL>", false);
        final Allele REF_ALLELE = Allele.create("N", true);
        while ((line = r.readLine()) != null) {
            if (StringUtil.isBlank(line))
                continue;
            tokens = tab.split(line);
            if (tokens.length != 3 + samples.size()) {
                r.close();
                vcw.close();
                throw new JvarkitException.TokenErrors("expected " + (samples.size() + 3) + "columns.", tokens);
            }
            final Set<Allele> alleles = new HashSet<>();
            alleles.add(REF_ALLELE);
            final VariantContextBuilder vcb = new VariantContextBuilder();
            vcb.chr(tokens[0]);
            vcb.start(Integer.parseInt(tokens[1]));
            final int chromEnd = Integer.parseInt(tokens[2]);
            vcb.stop(chromEnd);
            vcb.attribute(VCFConstants.END_KEY, chromEnd);
            if (dict != null) {
                final SAMSequenceRecord ssr = dict.getSequence(tokens[0]);
                if (ssr == null) {
                    LOG.error(JvarkitException.ContigNotFoundInDictionary.getMessage(tokens[0], dict));
                    return -1;
                }
                if (chromEnd > ssr.getSequenceLength()) {
                    LOG.warn("WARNING sequence length in " + line + " is greater than in dictionary ");
                }
            }
            int count_dup = 0;
            int count_del = 0;
            final Map<String, Float> sample2fold = new HashMap<>(samples.size());
            for (int i = 3; i < tokens.length; i++) {
                final String sampleName = samples.get(i - 3);
                final float f = Float.parseFloat(tokens[i]);
                if (f < 0 || Float.isNaN(f) || !Float.isFinite(f)) {
                    LOG.error("Bad fold " + f + " for sample " + sampleName + " in " + line);
                }
                sample2fold.put(sampleName, f);
            }
            final List<Genotype> genotypes = new ArrayList<>(samples.size());
            int count_sv_cases = 0;
            int count_sv_controls = 0;
            int count_ref_cases = 0;
            int count_ref_controls = 0;
            boolean got_sv = false;
            for (final String sampleName : sample2fold.keySet()) {
                final float normDepth = sample2fold.get(sampleName);
                final IndexCovUtils.SvType type = indexCovUtils.getType(normDepth);
                final GenotypeBuilder gb;
                switch(type) {
                    case REF:
                        {
                            gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, REF_ALLELE));
                            break;
                        }
                    case HET_DEL:
                        {
                            gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, DEL_ALLELE));
                            alleles.add(DEL_ALLELE);
                            count_del++;
                            break;
                        }
                    case HOM_DEL:
                        {
                            gb = new GenotypeBuilder(sampleName, Arrays.asList(DEL_ALLELE, DEL_ALLELE));
                            alleles.add(DEL_ALLELE);
                            count_del++;
                            vcb.filter(filterHomDel.getID());
                            break;
                        }
                    case HET_DUP:
                        {
                            gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, DUP_ALLELE));
                            alleles.add(DUP_ALLELE);
                            count_dup++;
                            break;
                        }
                    case HOM_DUP:
                        {
                            gb = new GenotypeBuilder(sampleName, Arrays.asList(DUP_ALLELE, DUP_ALLELE));
                            alleles.add(DUP_ALLELE);
                            vcb.filter(filterHomDup.getID());
                            count_dup++;
                            break;
                        }
                    default:
                        {
                            gb = new GenotypeBuilder(sampleName, Arrays.asList(Allele.NO_CALL, Allele.NO_CALL));
                            break;
                        }
                }
                if (type.isVariant())
                    got_sv = true;
                gb.attribute(foldHeader.getID(), normDepth);
                gb.GQ(type.getGenotypeQuality(normDepth));
                final Sample sn = pedigree.getSampleById(sampleName);
                if (sn != null) {
                    if (type.isVariant()) {
                        if (sn.isAffected()) {
                            count_sv_cases++;
                        } else if (sn.isUnaffected()) {
                            count_sv_controls++;
                        }
                    } else {
                        if (sn.isAffected()) {
                            count_ref_cases++;
                        } else if (sn.isUnaffected()) {
                            count_ref_controls++;
                        }
                    }
                }
                genotypes.add(gb.make());
            }
            vcb.alleles(alleles);
            if (!pedigree.isEmpty() && count_sv_cases == 1 && count_ref_cases > 0 && count_sv_controls == 0 && count_ref_controls > 0) {
                vcb.attribute(infoSingleton.getID(), Boolean.TRUE);
            } else if (!pedigree.isEmpty() && count_sv_cases > 0 && count_sv_cases == count_cases_in_pedigree && count_ref_cases == 0 && count_sv_controls == 0 && count_ref_controls > 0) {
                vcb.attribute(infoAllAffected.getID(), Boolean.TRUE);
            }
            vcb.genotypes(genotypes);
            if (count_dup == samples.size() && samples.size() != 1) {
                vcb.filter(filterAllDup.getID());
            }
            if (count_del == samples.size() && samples.size() != 1) {
                vcb.filter(filterAllDel.getID());
            }
            if (!got_sv) {
                vcb.filter(filterNoSV.getID());
            }
            vcb.attribute(infoNumDel.getID(), count_del);
            vcb.attribute(infoNumDup.getID(), count_dup);
            vcw.add(vcb.make());
        }
        vcw.close();
        vcw = null;
        r.close();
        r = null;
        return 0;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(r);
        CloserUtil.close(vcw);
    }
}
Also used : WritingVariantsDelegate(com.github.lindenb.jvarkit.variant.variantcontext.writer.WritingVariantsDelegate) Genotype(htsjdk.variant.variantcontext.Genotype) VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) Allele(htsjdk.variant.variantcontext.Allele) Arrays(java.util.Arrays) CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter) Program(com.github.lindenb.jvarkit.util.jcommander.Program) Parameter(com.beust.jcommander.Parameter) VCFStandardHeaderLines(htsjdk.variant.vcf.VCFStandardHeaderLines) VCFHeader(htsjdk.variant.vcf.VCFHeader) HashMap(java.util.HashMap) ParametersDelegate(com.beust.jcommander.ParametersDelegate) ArrayList(java.util.ArrayList) HashSet(java.util.HashSet) StringUtil(htsjdk.samtools.util.StringUtil) PedigreeParser(com.github.lindenb.jvarkit.pedigree.PedigreeParser) Map(java.util.Map) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) Path(java.nio.file.Path) CloserUtil(htsjdk.samtools.util.CloserUtil) VCFConstants(htsjdk.variant.vcf.VCFConstants) SequenceDictionaryUtils(com.github.lindenb.jvarkit.util.bio.SequenceDictionaryUtils) VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine) VCFHeaderLineType(htsjdk.variant.vcf.VCFHeaderLineType) Pedigree(com.github.lindenb.jvarkit.pedigree.Pedigree) GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder) Logger(com.github.lindenb.jvarkit.util.log.Logger) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) Set(java.util.Set) JVarkitVersion(com.github.lindenb.jvarkit.util.JVarkitVersion) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) List(java.util.List) FractionConverter(com.github.lindenb.jvarkit.jcommander.converter.FractionConverter) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine) VCFFormatHeaderLine(htsjdk.variant.vcf.VCFFormatHeaderLine) BufferedReader(java.io.BufferedReader) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) Sample(com.github.lindenb.jvarkit.pedigree.Sample) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) HashMap(java.util.HashMap) ArrayList(java.util.ArrayList) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) PedigreeParser(com.github.lindenb.jvarkit.pedigree.PedigreeParser) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine) VCFHeader(htsjdk.variant.vcf.VCFHeader) HashSet(java.util.HashSet) VCFFormatHeaderLine(htsjdk.variant.vcf.VCFFormatHeaderLine) CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter) Sample(com.github.lindenb.jvarkit.pedigree.Sample) Genotype(htsjdk.variant.variantcontext.Genotype) GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder) VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine) Allele(htsjdk.variant.variantcontext.Allele) Pedigree(com.github.lindenb.jvarkit.pedigree.Pedigree) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) BufferedReader(java.io.BufferedReader)

Example 24 with CharSplitter

use of com.github.lindenb.jvarkit.lang.CharSplitter in project jvarkit by lindenb.

the class VcfAncestralAllele method loadManifest.

private void loadManifest() throws IOException {
    BufferedReader r = null;
    try {
        final CharSplitter tab = CharSplitter.TAB;
        final CharSplitter pipe = CharSplitter.PIPE;
        r = IOUtils.openPathForBufferedReading(this.manifestFile);
        String line;
        while ((line = r.readLine()) != null) {
            if (line.isEmpty() || StringUtil.isBlank(line))
                continue;
            final String[] tokens = tab.split(line);
            if (tokens.length < 3)
                throw new JvarkitException.TokenErrors("expected two columns", tokens);
            final String ancestralContig = tokens[1];
            final Path fastaPath = Paths.get(tokens[2]);
            IOUtil.assertFileIsReadable(fastaPath);
            /* no, one fasta can contains more than one sequence
				if(contigName2ancestralFasta.values().contains(fastaPath)) {
					throw new IOException("fasta already defined for  "+fastaPath);
					}
				*/
            final Path faiPath = Paths.get(tokens[2] + ".fai");
            IOUtil.assertFileIsReadable(faiPath);
            for (final String sn : pipe.split(tokens[0])) {
                if (StringUtil.isBlank(line)) {
                    throw new IOException("empty contig name in " + line);
                }
                if (this.contigName2ancestralFasta.containsKey(sn)) {
                    throw new IOException("duplicate contig name in " + line);
                }
                this.contigName2ancestralFasta.put(sn, new Mapping(ancestralContig, fastaPath));
            }
        }
        if (this.contigName2ancestralFasta.isEmpty()) {
            LOG.warn("Manifest is empty " + this.manifestFile);
        }
    } finally {
        CloserUtil.close(r);
    }
}
Also used : JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) Path(java.nio.file.Path) CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter) BufferedReader(java.io.BufferedReader) IOException(java.io.IOException)

Example 25 with CharSplitter

use of com.github.lindenb.jvarkit.lang.CharSplitter in project jvarkit by lindenb.

the class SplitByTile method createKeys.

@Override
protected Set<Integer> createKeys(final SAMRecord rec) {
    final CharSplitter colon = CharSplitter.COLON;
    final String[] tokens = colon.split(rec.getReadName(), 6);
    if (tokens.length < 5) {
        LOG.error("Cannot get the 6th field in " + rec.getReadName());
        return null;
    }
    int tile = -1;
    try {
        tile = Integer.parseInt(tokens[4]);
    } catch (final Throwable e) {
        tile = -1;
    }
    if (tile < 0) {
        LOG.error("Bad tile in read: " + rec.getReadName());
        return null;
    }
    return Collections.singleton(tile);
}
Also used : CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter)

Aggregations

CharSplitter (com.github.lindenb.jvarkit.lang.CharSplitter)28 BufferedReader (java.io.BufferedReader)19 IOException (java.io.IOException)12 List (java.util.List)11 Parameter (com.beust.jcommander.Parameter)10 JvarkitException (com.github.lindenb.jvarkit.lang.JvarkitException)10 Program (com.github.lindenb.jvarkit.util.jcommander.Program)10 Logger (com.github.lindenb.jvarkit.util.log.Logger)10 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)10 Set (java.util.Set)10 Path (java.nio.file.Path)9 ArrayList (java.util.ArrayList)9 HashSet (java.util.HashSet)9 SequenceDictionaryUtils (com.github.lindenb.jvarkit.util.bio.SequenceDictionaryUtils)8 Launcher (com.github.lindenb.jvarkit.util.jcommander.Launcher)7 CloserUtil (htsjdk.samtools.util.CloserUtil)7 VariantContextBuilder (htsjdk.variant.variantcontext.VariantContextBuilder)7 VCFHeader (htsjdk.variant.vcf.VCFHeader)7 VCFInfoHeaderLine (htsjdk.variant.vcf.VCFInfoHeaderLine)7 Map (java.util.Map)7