Examples with VcfIterator com.github.lindenb.jvarkit.util.vcf.VcfIterator used on opensource projects

Example 46 with VcfIterator

use of com.github.lindenb.jvarkit.util.vcf.VcfIterator in project jvarkit by lindenb.

the class VCFMerge method copyTo.

private void copyTo(final InputStream in) throws IOException {
    final VcfIterator iter = VCFUtils.createVcfIteratorFromInputStream(in);
    final VariantContextWriter out = this.openVariantContextWriter(outputFile);
    VCFUtils.copyHeaderAndVariantsTo(iter, out);
    CloserUtil.close(out);
    CloserUtil.close(iter);
}

Example 47 with VcfIterator

use of com.github.lindenb.jvarkit.util.vcf.VcfIterator in project jvarkit by lindenb.

the class VcfPeekVcf method doVcfToVcf.

/**
 * public for knime
 */
@Override
public int doVcfToVcf(final String inputName, final VcfIterator vcfIn, final VariantContextWriter out) {
    try {
        final Set<String> unmatchedcontigs = new HashSet<>();
        final VCFHeader h = vcfIn.getHeader();
        final VCFHeader h2 = new VCFHeader(h);
        super.addMetaData(h2);
        final Map<String, VCFInfoHeaderLine> databaseTags = new HashMap<String, VCFInfoHeaderLine>();
        final VCFHeader databaseHeader = this.indexedVcfFileReader.getFileHeader();
        final ContigNameConverter nameConverter = (h.getSequenceDictionary() != null && !h.getSequenceDictionary().isEmpty() && databaseHeader.getSequenceDictionary() != null && !databaseHeader.getSequenceDictionary().isEmpty() ? ContigNameConverter.fromDictionaries(h.getSequenceDictionary(), databaseHeader.getSequenceDictionary()) : ContigNameConverter.getIdentity()).setOnNotFound(this.onContigNotFound);
        ;
        for (final String key : this.peek_info_tags) {
            VCFInfoHeaderLine hinfo = databaseHeader.getInfoHeaderLine(key);
            if (hinfo == null) {
                final String msg = "INFO name=" + key + " missing in " + this.resourceVcfFile;
                if (this.missingIdIsError) {
                    LOG.warn(msg);
                    continue;
                } else {
                    LOG.error(msg);
                    return -1;
                }
            }
            switch(hinfo.getCountType()) {
                case G:
                    throw new JvarkitException.UserError("Cannot handle VCFHeaderLineCount.G for " + hinfo.getID());
                default:
                    databaseTags.put(hinfo.getID(), hinfo);
                    break;
            }
            hinfo = VCFUtils.renameVCFInfoHeaderLine(hinfo, this.peekTagPrefix + key);
            if (h2.getInfoHeaderLine(hinfo.getID()) != null) {
                throw new JvarkitException.UserError("key " + this.peekTagPrefix + key + " already defined in VCF header");
            }
            h2.addMetaDataLine(hinfo);
            ;
        }
        out.writeHeader(h2);
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(h).logger(LOG);
        while (vcfIn.hasNext()) {
            final VariantContext ctx = progress.watch(vcfIn.next());
            final String outContig = nameConverter.apply(ctx.getContig());
            if (outContig == null) {
                unmatchedcontigs.add(ctx.getContig());
                continue;
            }
            final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
            CloseableIterator<VariantContext> iter = this.indexedVcfFileReader.query(outContig, Math.max(0, ctx.getStart() - 1), (ctx.getEnd() + 1));
            while (iter.hasNext()) {
                final VariantContext ctx2 = iter.next();
                if (!outContig.equals(ctx2.getContig()))
                    continue;
                if (ctx.getStart() != ctx2.getStart())
                    continue;
                if (!ctx.getReference().equals(ctx2.getReference()))
                    continue;
                final boolean okAllele;
                switch(altAlleleMatcher) {
                    case all:
                        {
                            okAllele = ctx.getAlternateAlleles().stream().filter(A -> ctx2.hasAlternateAllele(A)).count() == ctx.getAlternateAlleles().size();
                            break;
                        }
                    case at_least_one:
                        {
                            okAllele = ctx.getAlternateAlleles().stream().filter(A -> ctx2.hasAlternateAllele(A)).findAny().isPresent();
                            break;
                        }
                    case none:
                        okAllele = true;
                        break;
                    default:
                        throw new IllegalStateException(altAlleleMatcher.name());
                }
                if (!okAllele)
                    continue;
                if (this.peekId && ctx2.hasID()) {
                    vcb.id(ctx2.getID());
                }
                boolean somethingWasChanged = false;
                for (final String key : databaseTags.keySet()) {
                    if (!ctx2.hasAttribute(key))
                        continue;
                    final VCFInfoHeaderLine dbHeader = databaseTags.get(key);
                    switch(dbHeader.getCountType()) {
                        case A:
                            {
                                final List<Object> newatt = new ArrayList<>();
                                final List<Object> ctx2att = ctx2.getAttributeAsList(key);
                                for (int i = 0; i < ctx.getAlternateAlleles().size(); ++i) {
                                    final Allele ctxalt = ctx.getAlternateAllele(i);
                                    int index2 = ctx2.getAlternateAlleles().indexOf(ctxalt);
                                    if (index2 == -1 || index2 >= ctx2att.size()) {
                                        newatt.add(null);
                                    } else {
                                        newatt.add(ctx2att.get(index2));
                                    }
                                }
                                if (newatt.stream().filter(Obj -> !(Obj == null || VCFConstants.EMPTY_INFO_FIELD.equals(Obj))).count() > 0) {
                                    vcb.attribute(this.peekTagPrefix + key, newatt);
                                    somethingWasChanged = true;
                                }
                                break;
                            }
                        case R:
                            {
                                final List<Object> newatt = new ArrayList<>();
                                final List<Object> ctx2att = ctx2.getAttributeAsList(key);
                                for (int i = 0; i < ctx.getAlleles().size(); ++i) {
                                    final Allele ctxalt = ctx.getAlleles().get(i);
                                    int index2 = ctx2.getAlleleIndex(ctxalt);
                                    if (index2 == -1 || index2 >= ctx2att.size()) {
                                        newatt.add(null);
                                    } else {
                                        newatt.add(ctx2att.get(index2));
                                    }
                                }
                                if (newatt.stream().filter(Obj -> !(Obj == null || VCFConstants.EMPTY_INFO_FIELD.equals(Obj))).count() > 0) {
                                    vcb.attribute(this.peekTagPrefix + key, newatt);
                                    somethingWasChanged = true;
                                }
                                break;
                            }
                        default:
                            {
                                final Object o = ctx2.getAttribute(key);
                                vcb.attribute(this.peekTagPrefix + key, o);
                                somethingWasChanged = true;
                                break;
                            }
                    }
                }
                if (somethingWasChanged)
                    break;
            }
            iter.close();
            iter = null;
            out.add(vcb.make());
            if (out.checkError())
                break;
        }
        progress.finish();
        if (!unmatchedcontigs.isEmpty()) {
            LOG.debug("Unmatched contigs: " + unmatchedcontigs.stream().collect(Collectors.joining("; ")));
        }
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    }
}

Example 48 with VcfIterator

use of com.github.lindenb.jvarkit.util.vcf.VcfIterator 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);
    }
}

Example 49 with VcfIterator

use of com.github.lindenb.jvarkit.util.vcf.VcfIterator in project jvarkit by lindenb.

the class XContaminations method doWork.

@Override
public int doWork(final List<String> args) {
    long last_save_ms = System.currentTimeMillis();
    if (this.output_as_vcf && !this.use_only_sample_name) {
        LOG.error("cannot write vcf if --sample is not set");
        return -1;
    }
    if (args.size() < 2) {
        LOG.error("Illegal Number of args");
        return -1;
    }
    final Set<File> bamFiles = IOUtils.unrollFiles(args.subList(1, args.size())).stream().map(S -> new File(S)).collect(Collectors.toSet());
    if (bamFiles.isEmpty()) {
        LOG.error("Undefined BAM file(s)");
        return -1;
    }
    SAMRecordIterator iter = null;
    VcfIterator in = null;
    Map<String, SamReader> sample2samReader = new HashMap<>();
    VariantContextWriter vcfw = null;
    try {
        final SamReaderFactory srf = super.createSamReaderFactory();
        if (args.get(0).equals("-")) {
            in = super.openVcfIterator(null);
        } else {
            in = super.openVcfIterator(args.get(0));
        }
        VCFHeader vcfHeader = in.getHeader();
        final SAMSequenceDictionary dict1 = vcfHeader.getSequenceDictionary();
        if (dict1 == null) {
            LOG.error(JvarkitException.VcfDictionaryMissing.getMessage(args.get(0)));
            return -1;
        }
        final Set<String> sampleNames = new HashSet<>(vcfHeader.getSampleNamesInOrder());
        if (sampleNames.isEmpty()) {
            LOG.error("VCF contains no sample");
            return -1;
        }
        for (final File bamFile : bamFiles) {
            LOG.info("Opening " + bamFile);
            final SamReader samReader = srf.open(bamFile);
            final SAMFileHeader samHeader = samReader.getFileHeader();
            final SAMSequenceDictionary dict2 = samHeader.getSequenceDictionary();
            if (dict2 == null) {
                samReader.close();
                LOG.error(JvarkitException.BamDictionaryMissing.getMessage(bamFile.getPath()));
                return -1;
            }
            if (!SequenceUtil.areSequenceDictionariesEqual(dict1, dict2)) {
                samReader.close();
                LOG.error(JvarkitException.DictionariesAreNotTheSame.getMessage(dict1, dict2));
                return -1;
            }
            if (!samReader.hasIndex()) {
                samReader.close();
                LOG.error("sam is not indexed : " + bamFile);
                return -1;
            }
            String sampleName = null;
            for (final SAMReadGroupRecord rgr : samHeader.getReadGroups()) {
                final String s = rgr.getSample();
                if (StringUtil.isBlank(s))
                    continue;
                if (sampleName == null) {
                    sampleName = s;
                } else if (!sampleName.equals(s)) {
                    samReader.close();
                    LOG.error("Cannot handle more than one sample/bam  " + bamFile + " " + sampleName);
                    return -1;
                }
            }
            if (sampleName == null) {
                samReader.close();
                LOG.error("No sample in " + bamFile);
                // skip this bam
                continue;
            }
            if (!sampleNames.contains(sampleName)) {
                samReader.close();
                LOG.error("Not in VCF header: sample " + sampleName + " " + bamFile);
                // skip this bam
                continue;
            }
            if (sample2samReader.containsKey(sampleName)) {
                samReader.close();
                LOG.error("Cannot handle more than one bam/sample: " + bamFile + " " + sampleName);
                return -1;
            }
            sample2samReader.put(sampleName, samReader);
        }
        if (sample2samReader.size() < 2) {
            LOG.error("Not engough BAM/samples. Expected at least two valid BAMs");
            return -1;
        }
        sampleNames.retainAll(sample2samReader.keySet());
        /* create a VCF is VCF output asked */
        final List<SamplePair> sampleListForVcf;
        if (this.output_as_vcf) {
            vcfw = super.openVariantContextWriter(outputFile);
            final Set<VCFHeaderLine> metaData = new HashSet<>();
            metaData.add(new VCFFormatHeaderLine("S1S1", 1, VCFHeaderLineType.Integer, "reads sample 1 supporting sample 1"));
            metaData.add(new VCFFormatHeaderLine("S1S2", 1, VCFHeaderLineType.Integer, "reads sample 1 supporting sample 2"));
            metaData.add(new VCFFormatHeaderLine("S1SO", 1, VCFHeaderLineType.Integer, "reads sample 1 supporting others"));
            metaData.add(new VCFFormatHeaderLine("S2S1", 1, VCFHeaderLineType.Integer, "reads sample 2 supporting sample 1"));
            metaData.add(new VCFFormatHeaderLine("S2S2", 1, VCFHeaderLineType.Integer, "reads sample 2 supporting sample 2"));
            metaData.add(new VCFFormatHeaderLine("S2SO", 1, VCFHeaderLineType.Integer, "reads sample 2 supporting others"));
            metaData.add(new VCFFormatHeaderLine("FR", 1, VCFHeaderLineType.Float, "Fraction. '-1' for unavailable."));
            metaData.add(new VCFFormatHeaderLine("S1A", 1, VCFHeaderLineType.Character, "sample 1 allele"));
            metaData.add(new VCFFormatHeaderLine("S2A", 1, VCFHeaderLineType.Character, "sample 2 allele"));
            metaData.add(new VCFFilterHeaderLine("XCONTAMINATION", "Fraction test is > " + fraction_treshold));
            metaData.add(new VCFFilterHeaderLine("BADSAMPLES", "At least one pair of genotype fails the 'LE' test"));
            metaData.add(new VCFInfoHeaderLine("LE", 1, VCFHeaderLineType.Integer, "number of pair of genotypes having (S1S1<=S1S2 or S2S2<=S2S1)."));
            metaData.add(new VCFInfoHeaderLine("BADSAMPLES", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, "Samples founds failing the 'LE' test"));
            sampleListForVcf = new ArrayList<>();
            final List<String> sampleList = new ArrayList<>(sampleNames);
            for (int x = 0; x + 1 < sampleList.size(); ++x) {
                for (int y = x + 1; y < sampleList.size(); ++y) {
                    sampleListForVcf.add(new SamplePair(new SimpleSampleIdenfifier(sampleList.get(x)), new SimpleSampleIdenfifier(sampleList.get(y))));
                }
            }
            final VCFHeader header2 = new VCFHeader(metaData, sampleListForVcf.stream().map(V -> V.getLabel()).sorted().collect(Collectors.toList()));
            header2.setSequenceDictionary(dict1);
            vcfw.writeHeader(header2);
        } else {
            vcfw = null;
            sampleListForVcf = null;
        }
        final Map<SamplePair, SampleAlleles> contaminationTable = new HashMap<>();
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(dict1).logger(LOG);
        while (in.hasNext()) {
            final VariantContext ctx = progress.watch(in.next());
            if (!ctx.isSNP() || ctx.isFiltered() || !ctx.isBiallelic() || ctx.isSymbolic() || !this.variantFilter.test(ctx)) {
                continue;
            }
            int count_homref = 0;
            int count_homvar = 0;
            int count_het = 0;
            final Map<String, Genotype> sample2gt = new HashMap<>();
            for (int gidx = 0; gidx < ctx.getNSamples(); ++gidx) {
                final Genotype G = ctx.getGenotype(gidx);
                if (!G.isCalled())
                    continue;
                if (G.isHet()) {
                    // here because in use_singleton we must be sure that there is only one hom_var
                    count_het++;
                    if (this.use_singleton && count_het > 0)
                        break;
                } else if (G.isHomVar()) {
                    // here because in use_singleton we must be sure that there is only one hom_var
                    count_homvar++;
                    if (this.use_singleton && count_homvar > 1)
                        break;
                }
                if (G.isFiltered())
                    continue;
                if (!sample2samReader.containsKey(G.getSampleName()))
                    continue;
                if (!sampleNames.contains(G.getSampleName()))
                    continue;
                if (!this.genotypeFilter.test(ctx, G))
                    continue;
                sample2gt.put(G.getSampleName(), G);
            }
            if (this.use_singleton && count_het > 0)
                continue;
            if (this.use_singleton && count_homvar > 1)
                continue;
            if (sample2gt.size() < 2)
                continue;
            // reset and recount
            count_homref = 0;
            count_homvar = 0;
            count_het = 0;
            for (final String sampleName : sample2gt.keySet()) {
                final Genotype G = ctx.getGenotype(sampleName);
                switch(G.getType()) {
                    case HOM_REF:
                        count_homref++;
                        break;
                    case HOM_VAR:
                        count_homvar++;
                        break;
                    case HET:
                        count_het++;
                        break;
                    default:
                        break;
                }
            }
            // singleton check
            if (this.use_singleton && (count_het > 0 || count_homvar != 1)) {
                continue;
            }
            // at least one HOM_REF and one HOM_VAR
            if (count_homref == 0)
                continue;
            if (count_homvar == 0)
                continue;
            final Map<SampleIdentifier, Counter<Character>> sample_identifier_2allelesCount = new HashMap<>();
            /* scan Reads for those Genotype/Samples */
            for (final String sampleName : sample2gt.keySet()) {
                if (!sample2samReader.containsKey(sampleName))
                    continue;
                // sample name is not in vcf header
                final SamReader samReader = sample2samReader.get(sampleName);
                if (samReader == null)
                    continue;
                final Genotype genotype = sample2gt.get(sampleName);
                if (genotype == null)
                    continue;
                iter = samReader.query(ctx.getContig(), ctx.getStart(), ctx.getEnd(), false);
                while (iter.hasNext()) {
                    final SAMRecord record = iter.next();
                    if (record.getEnd() < ctx.getStart())
                        continue;
                    if (ctx.getEnd() < record.getStart())
                        continue;
                    if (record.getReadUnmappedFlag())
                        continue;
                    if (this.filter.filterOut(record))
                        continue;
                    final SAMReadGroupRecord srgr = record.getReadGroup();
                    // not current sample
                    if (srgr == null)
                        continue;
                    if (!sampleName.equals(srgr.getSample()))
                        continue;
                    final Cigar cigar = record.getCigar();
                    if (cigar == null || cigar.isEmpty())
                        continue;
                    byte[] readSeq = record.getReadBases();
                    if (readSeq == null || readSeq.length == 0)
                        continue;
                    int readPos = record.getReadPositionAtReferencePosition(ctx.getStart());
                    if (readPos < 1)
                        continue;
                    readPos--;
                    if (readPos >= readSeq.length)
                        continue;
                    final char base = Character.toUpperCase((char) readSeq[readPos]);
                    if (base == 'N')
                        continue;
                    final SampleIdentifier sampleIdentifier;
                    if (this.use_only_sample_name) {
                        sampleIdentifier = new SimpleSampleIdenfifier(sampleName);
                    } else {
                        final ShortReadName readName = ShortReadName.parse(record);
                        if (!readName.isValid()) {
                            LOG.info("No a valid read name " + record.getReadName());
                            continue;
                        }
                        sampleIdentifier = new SequencerFlowCellRunLaneSample(readName, sampleName);
                    }
                    Counter<Character> sampleAlleles = sample_identifier_2allelesCount.get(sampleIdentifier);
                    if (sampleAlleles == null) {
                        sampleAlleles = new Counter<Character>();
                        sample_identifier_2allelesCount.put(sampleIdentifier, sampleAlleles);
                    }
                    sampleAlleles.incr(base);
                }
                iter.close();
                iter = null;
            }
            /* end scan reads for this sample */
            /* sum-up data for this SNP */
            final VariantContextBuilder vcb;
            final List<Genotype> genotypeList;
            if (this.output_as_vcf) {
                vcb = new VariantContextBuilder(args.get(0), ctx.getContig(), ctx.getStart(), ctx.getEnd(), ctx.getAlleles());
                if (ctx.hasID())
                    vcb.id(ctx.getID());
                genotypeList = new ArrayList<>();
            } else {
                vcb = null;
                genotypeList = null;
            }
            for (final String sample1 : sample2gt.keySet()) {
                final Genotype g1 = sample2gt.get(sample1);
                final char a1 = g1.getAllele(0).getBaseString().charAt(0);
                for (final String sample2 : sample2gt.keySet()) {
                    if (sample1.compareTo(sample2) >= 0)
                        continue;
                    final Genotype g2 = sample2gt.get(sample2);
                    if (g2.sameGenotype(g1))
                        continue;
                    final char a2 = g2.getAllele(0).getBaseString().charAt(0);
                    for (final SampleIdentifier sfcr1 : sample_identifier_2allelesCount.keySet()) {
                        if (!sfcr1.getSampleName().equals(sample1))
                            continue;
                        final Counter<Character> counter1 = sample_identifier_2allelesCount.get(sfcr1);
                        if (counter1 == null)
                            continue;
                        for (final SampleIdentifier sfcr2 : sample_identifier_2allelesCount.keySet()) {
                            if (!sfcr2.getSampleName().equals(sample2))
                                continue;
                            final SamplePair samplePair = new SamplePair(sfcr1, sfcr2);
                            final Counter<Character> counter2 = sample_identifier_2allelesCount.get(sfcr2);
                            if (counter2 == null)
                                continue;
                            SampleAlleles sampleAlleles = contaminationTable.get(samplePair);
                            if (sampleAlleles == null) {
                                sampleAlleles = new SampleAlleles();
                                contaminationTable.put(samplePair, sampleAlleles);
                                if (!this.output_as_vcf && contaminationTable.size() % 10000 == 0)
                                    LOG.info("n(pairs)=" + contaminationTable.size());
                            }
                            sampleAlleles.number_of_comparaisons++;
                            for (final Character allele : counter1.keySet()) {
                                final long n = counter1.count(allele);
                                if (allele.equals(a1)) {
                                    sampleAlleles.reads_sample1_supporting_sample1 += n;
                                } else if (allele.equals(a2)) {
                                    sampleAlleles.reads_sample1_supporting_sample2 += n;
                                } else {
                                    sampleAlleles.reads_sample1_supporting_other += n;
                                }
                            }
                            for (final Character allele : counter2.keySet()) {
                                final long n = counter2.count(allele);
                                if (allele.equals(a2)) {
                                    sampleAlleles.reads_sample2_supporting_sample2 += n;
                                } else if (allele.equals(a1)) {
                                    sampleAlleles.reads_sample2_supporting_sample1 += n;
                                } else {
                                    sampleAlleles.reads_sample2_supporting_other += n;
                                }
                            }
                        }
                    }
                }
            }
            if (this.output_as_vcf) {
                final Set<String> bad_samples = new TreeSet<>();
                boolean fraction_flag = false;
                int num_lt = 0;
                for (final SamplePair samplepair : sampleListForVcf) {
                    final GenotypeBuilder gb = new GenotypeBuilder(samplepair.getLabel());
                    final SampleAlleles sampleAlleles = contaminationTable.get(samplepair);
                    if (sampleAlleles != null) {
                        gb.attribute("S1S1", sampleAlleles.reads_sample1_supporting_sample1);
                        gb.attribute("S1S2", sampleAlleles.reads_sample1_supporting_sample2);
                        gb.attribute("S1SO", sampleAlleles.reads_sample1_supporting_other);
                        gb.attribute("S2S1", sampleAlleles.reads_sample2_supporting_sample1);
                        gb.attribute("S2S2", sampleAlleles.reads_sample2_supporting_sample2);
                        gb.attribute("S2SO", sampleAlleles.reads_sample2_supporting_other);
                        gb.attribute("S1A", sample2gt.get(samplepair.sample1.getSampleName()).getAllele(0).getDisplayString().charAt(0));
                        gb.attribute("S2A", sample2gt.get(samplepair.sample2.getSampleName()).getAllele(0).getDisplayString().charAt(0));
                        final double fraction = sampleAlleles.getFraction();
                        gb.attribute("FR", fraction);
                        if (!this.passFractionTreshold.test(fraction)) {
                            fraction_flag = true;
                        }
                        boolean bad_lt_flag = false;
                        if (sampleAlleles.reads_sample1_supporting_sample1 <= this.fail_factor * sampleAlleles.reads_sample1_supporting_sample2) {
                            bad_samples.add(samplepair.sample1.getSampleName());
                            bad_lt_flag = true;
                        }
                        if (sampleAlleles.reads_sample2_supporting_sample2 <= this.fail_factor * sampleAlleles.reads_sample2_supporting_sample1) {
                            bad_samples.add(samplepair.sample2.getSampleName());
                            bad_lt_flag = true;
                        }
                        if (bad_lt_flag) {
                            num_lt++;
                        }
                    } else {
                        gb.attribute("S1S1", -1);
                        gb.attribute("S1S2", -1);
                        gb.attribute("S1SO", -1);
                        gb.attribute("S2S1", -1);
                        gb.attribute("S2S2", -1);
                        gb.attribute("S2SO", -1);
                        gb.attribute("S1A", '.');
                        gb.attribute("S2A", '.');
                        gb.attribute("FR", -1f);
                    }
                    genotypeList.add(gb.make());
                }
                if (!bad_samples.isEmpty()) {
                    vcb.attribute("BADSAMPLES", new ArrayList<>(bad_samples));
                }
                vcb.attribute("LE", num_lt);
                if (fraction_flag || !bad_samples.isEmpty()) {
                    if (fraction_flag)
                        vcb.filter("XCONTAMINATION");
                    if (!bad_samples.isEmpty())
                        vcb.filter("BADSAMPLES");
                } else {
                    vcb.passFilters();
                }
                vcb.genotypes(genotypeList);
                vcfw.add(vcb.make());
                contaminationTable.clear();
            } else {
                final long now = System.currentTimeMillis();
                if (this.outputFile != null && this.save_every_sec > -1L && last_save_ms + (this.save_every_sec * 1000L) > now) {
                    saveToFile(contaminationTable);
                    last_save_ms = now;
                }
            }
        }
        progress.finish();
        if (this.output_as_vcf) {
            vcfw.close();
            vcfw = null;
        } else {
            saveToFile(contaminationTable);
        }
        return 0;
    } catch (final Exception e) {
        LOG.error(e);
        return -1;
    } finally {
        CloserUtil.close(vcfw);
        CloserUtil.close(in);
        CloserUtil.close(iter);
        for (SamReader samReader : sample2samReader.values()) CloserUtil.close(samReader);
        sample2samReader.clear();
    }
}

Example 50 with VcfIterator

use of com.github.lindenb.jvarkit.util.vcf.VcfIterator in project jvarkit by lindenb.

the class BioAlcidae method executeAsVcf.

/**
 * moved to public for knime
 */
public int executeAsVcf(final String source) throws IOException {
    LOG.info("source: " + source);
    VcfIterator in = null;
    try {
        in = VCFUtils.createVcfIterator(source);
        return executeAsVcf(in);
    } catch (Exception e) {
        LOG.error(e);
        return -1;
    } finally {
        CloserUtil.close(in);
    }
}