Examples with GenotypeBuilder htsjdk.variant.variantcontext.GenotypeBuilder used on opensource projects

Example 41 with GenotypeBuilder

use of htsjdk.variant.variantcontext.GenotypeBuilder in project jvarkit by lindenb.

the class LumpySort method doWork.

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

Example 42 with GenotypeBuilder

use of htsjdk.variant.variantcontext.GenotypeBuilder in project jvarkit by lindenb.

the class FixVcfMissingGenotypes method doVcfToVcf.

@Override
protected int doVcfToVcf(final String inputName, final VcfIterator in, final VariantContextWriter out) {
    final Set<String> bamFiles = IOUtils.unrollFiles(bamList);
    final Map<String, List<SamReader>> sample2bam = new HashMap<>(bamFiles.size());
    try {
        final VCFHeader header = in.getHeader();
        for (final String bamFile : bamFiles) {
            LOG.info("Reading header for " + bamFile);
            final SamReader reader = super.openSamReader(bamFile);
            if (!reader.hasIndex()) {
                LOG.error("No BAM index available for " + bamFile);
                return -1;
            }
            final SAMFileHeader samHeader = reader.getFileHeader();
            for (final SAMReadGroupRecord g : samHeader.getReadGroups()) {
                if (g.getSample() == null)
                    continue;
                final String sample = g.getSample();
                if (StringUtil.isBlank(sample))
                    continue;
                List<SamReader> readers = sample2bam.get(sample);
                if (readers == null) {
                    readers = new ArrayList<>();
                    sample2bam.put(sample, readers);
                }
                readers.add(reader);
            }
        }
        final VCFHeader h2 = new VCFHeader(header);
        if (h2.getFormatHeaderLine(VCFConstants.DEPTH_KEY) == null) {
            h2.addMetaDataLine(VCFStandardHeaderLines.getFormatLine(VCFConstants.DEPTH_KEY));
        }
        if (h2.getFormatHeaderLine(VCFConstants.GENOTYPE_KEY) == null) {
            h2.addMetaDataLine(VCFStandardHeaderLines.getFormatLine(VCFConstants.GENOTYPE_KEY));
        }
        if (h2.getFormatHeaderLine(VCFConstants.GENOTYPE_FILTER_KEY) == null) {
            h2.addMetaDataLine(VCFStandardHeaderLines.getFormatLine(VCFConstants.GENOTYPE_FILTER_KEY));
        }
        h2.addMetaDataLine(new VCFFormatHeaderLine(this.fixedTag, 1, VCFHeaderLineType.Integer, "Genotype was set as homozygous (min depth =" + this.minDepth + ")"));
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(header);
        this.recalculator.setHeader(h2);
        out.writeHeader(h2);
        while (in.hasNext()) {
            final VariantContext ctx = progress.watch(in.next());
            boolean somethingWasChanged = false;
            final List<Genotype> genotypes = new ArrayList<>(ctx.getNSamples());
            for (int i = 0; i < ctx.getNSamples(); ++i) {
                Genotype genotype = ctx.getGenotype(i);
                final String sample = genotype.getSampleName();
                if (!genotype.isCalled() || (!genotype.hasDP() && this.fixDP)) {
                    List<SamReader> samReaders = sample2bam.get(sample);
                    if (samReaders == null)
                        samReaders = Collections.emptyList();
                    final int depth = fetchDP(ctx, sample, samReaders);
                    if (genotype.isCalled() && !genotype.hasDP()) {
                        genotype = new GenotypeBuilder(genotype).DP(depth).make();
                        somethingWasChanged = true;
                    } else // genotype was not called
                    {
                        if (depth >= this.minDepth) {
                            final List<Allele> homozygous = new ArrayList<>(2);
                            homozygous.add(ctx.getReference());
                            homozygous.add(ctx.getReference());
                            final GenotypeBuilder gb = new GenotypeBuilder(genotype);
                            gb.alleles(homozygous);
                            gb.attribute(this.fixedTag, 1);
                            gb.DP(depth);
                            if (!StringUtil.isBlank(this.fixedGenotypesAreFiltered))
                                gb.filter(this.fixedGenotypesAreFiltered);
                            somethingWasChanged = true;
                            genotype = gb.make();
                        } else if (// cannot fix but we can update DP
                        !genotype.hasDP()) {
                            genotype = new GenotypeBuilder(genotype).DP(depth).make();
                            somethingWasChanged = true;
                        }
                    }
                }
                genotypes.add(genotype);
            }
            if (somethingWasChanged) {
                final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
                vcb.genotypes(genotypes);
                out.add(this.recalculator.apply(vcb.make()));
            } else {
                out.add(ctx);
            }
        }
        progress.finish();
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        sample2bam.values().stream().flatMap(L -> L.stream()).forEach(R -> CloserUtil.close(R));
    }
}

Example 43 with GenotypeBuilder

use of htsjdk.variant.variantcontext.GenotypeBuilder in project jvarkit by lindenb.

the class SamScanSplitReads method saveAsVcf.

private void saveAsVcf(Set<String> sampleNames, SAMSequenceDictionary dict) throws IOException {
    final Function<String, Integer> contig2tid = S -> {
        int i = dict.getSequenceIndex(S);
        if (i == -1)
            throw new IllegalArgumentException("cannot find contig " + S + " in dictionary");
        return i;
    };
    final Comparator<Interval> intervalComparator = new Comparator<Interval>() {

        @Override
        public int compare(final Interval r1, final Interval r2) {
            final int cmp = contig2tid.apply(r1.getContig()) - contig2tid.apply(r2.getContig());
            if (cmp != 0)
                return cmp;
            return r1.getStart() - r2.getStart();
        }
    };
    final Allele REF = Allele.create("N", true);
    final Set<VCFHeaderLine> meta = new HashSet<>();
    VCFStandardHeaderLines.addStandardFormatLines(meta, false, VCFConstants.GENOTYPE_KEY, VCFConstants.DEPTH_KEY);
    sampleNames.addAll(this.sample2database.keySet());
    if (sampleNames.isEmpty())
        sampleNames.add(this.defaultSampleName);
    VCFHeader header = new VCFHeader(meta, sampleNames);
    header.setSequenceDictionary(dict);
    VariantContextWriter vcw = super.openVariantContextWriter(outputFile);
    vcw.writeHeader(header);
    final List<Arc> all_arcs = new ArrayList<>();
    for (final String sample : this.sample2database.keySet()) {
        final IntervalTreeMap<Set<Arc>> database = this.sample2database.get(sample);
        for (final Interval interval : database.keySet()) {
            for (final Arc arc : database.get(interval)) {
                all_arcs.add(arc);
            }
        }
    }
    Collections.sort(all_arcs, (A1, A2) -> intervalComparator.compare(A1.intervalFrom, A2.intervalFrom));
    for (final Arc row : all_arcs) {
        final List<Genotype> genotypes = new ArrayList<>();
        final Set<Allele> alleles = new HashSet<>();
        for (final String sample : this.sample2database.keySet()) {
            final IntervalTreeMap<Set<Arc>> database = this.sample2database.get(sample);
            for (final Arc arc : database.get(row.intervalFrom)) {
                if (!arc.equals(row))
                    continue;
                final Allele alt = Allele.create(new StringBuilder().append("<").append(arc.intervalFrom.getContig()).append(":").append(arc.intervalFrom.getStart()).append("-").append(arc.intervalFrom.getEnd()).append("|").append(arc.intervalTo.getContig()).append(":").append(arc.intervalTo.getStart()).append("-").append(arc.intervalTo.getEnd()).append(">").toString(), false);
                alleles.add(alt);
                final Genotype g = new GenotypeBuilder(sample).alleles(Collections.singletonList(alt)).DP(arc.countSupportingReads).make();
                genotypes.add(g);
            }
        }
        alleles.add(REF);
        VariantContextBuilder vcb = new VariantContextBuilder().chr(row.intervalFrom.getContig()).start(row.intervalFrom.getStart()).stop(row.intervalFrom.getStart()).alleles(alleles).genotypes(genotypes);
        vcw.add(vcb.make());
    }
    vcw.close();
}

Example 44 with GenotypeBuilder

use of htsjdk.variant.variantcontext.GenotypeBuilder in project jvarkit by lindenb.

the class VCFReplaceTag method doVcfToVcf.

@Override
protected int doVcfToVcf(final String inputName, final VcfIterator r, final VariantContextWriter w) {
    final VCFHeader header = r.getHeader();
    final HashSet<VCFHeaderLine> copyMeta = new HashSet<>(header.getMetaDataInInputOrder());
    for (final String key : this.transformMap.keySet()) {
        switch(this.replaceTypeNo) {
            case // INFO
            0:
                {
                    final VCFInfoHeaderLine info = header.getInfoHeaderLine(key);
                    if (info != null) {
                        copyMeta.remove(info);
                        copyMeta.add(VCFUtils.renameVCFInfoHeaderLine(info, this.transformMap.get(key)));
                    }
                    break;
                }
            case // FORMAT
            1:
                {
                    final VCFFormatHeaderLine fmt = header.getFormatHeaderLine(key);
                    if (fmt != null) {
                        copyMeta.remove(fmt);
                        copyMeta.add(VCFUtils.renameVCFFormatHeaderLine(fmt, this.transformMap.get(key)));
                    }
                    break;
                }
            case // FILTER
            2:
                {
                    final VCFFilterHeaderLine filter = header.getFilterHeaderLine(key);
                    if (filter != null) {
                        copyMeta.remove(filter);
                        copyMeta.add(VCFUtils.renameVCFFilterHeaderLine(filter, this.transformMap.get(key)));
                    }
                    break;
                }
            default:
                throw new IllegalStateException("" + this.replaceTypeNo);
        }
    }
    final VCFHeader h2 = new VCFHeader(copyMeta, header.getSampleNamesInOrder());
    addMetaData(h2);
    final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(h2);
    w.writeHeader(h2);
    while (r.hasNext()) {
        VariantContext ctx = progress.watch(r.next());
        VariantContextBuilder b = new VariantContextBuilder(ctx);
        switch(this.replaceTypeNo) {
            case // INFO
            0:
                {
                    for (String key : this.transformMap.keySet()) {
                        Object o = ctx.getAttribute(key);
                        if (o != null) {
                            b.rmAttribute(key);
                            b.attribute(this.transformMap.get(key), o);
                        }
                    }
                    break;
                }
            case // FORMAT
            1:
                {
                    List<Genotype> newgenotypes = new ArrayList<>(ctx.getNSamples());
                    for (int i = 0; i < ctx.getNSamples(); ++i) {
                        Genotype g = ctx.getGenotype(i);
                        Map<String, Object> atts = g.getExtendedAttributes();
                        GenotypeBuilder gb = new GenotypeBuilder(g);
                        for (String key : this.transformMap.keySet()) {
                            Object o = atts.get(key);
                            if (o != null) {
                                atts.remove(key);
                                atts.put(this.transformMap.get(key), o);
                            }
                        }
                        gb.attributes(atts);
                        newgenotypes.add(gb.make());
                    }
                    b.genotypes(newgenotypes);
                    break;
                }
            case // FILTER
            2:
                {
                    Set<String> filters = new HashSet<>(ctx.getFilters());
                    for (String key : this.transformMap.keySet()) {
                        if (filters.contains(key)) {
                            filters.remove(key);
                            filters.add(this.transformMap.get(key));
                        }
                    }
                    b.filters(filters);
                    break;
                }
            default:
                throw new IllegalStateException("" + this.replaceTypeNo);
        }
        w.add(b.make());
        if (w.checkError())
            break;
    }
    progress.finish();
    LOG.info("done");
    return 0;
}

Example 45 with GenotypeBuilder

use of htsjdk.variant.variantcontext.GenotypeBuilder 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();
    }
}