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

Example 51 with Genotype

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

the class IgvReview method saveVariants.

private File saveVariants(Window owner, final File f) {
    if (f == null)
        return null;
    VariantContextWriterBuilder vcb = new VariantContextWriterBuilder();
    VariantContextWriter w = null;
    try {
        SAMSequenceDictionary dict = this.vcfHeader.getSequenceDictionary();
        if (dict != null) {
            vcb.setReferenceDictionary(dict);
        }
        vcb.setOutputFile(f);
        final VCFHeader header2 = new VCFHeader(this.vcfHeader);
        if (header2.getFormatHeaderLine(this.reviewFormat.getID()) == null) {
            header2.addMetaDataLine(this.reviewFormat);
        }
        w = vcb.build();
        w.writeHeader(header2);
        int x = 0;
        while (x < this.genotypeTable.getItems().size()) {
            VariantContext v1 = this.genotypeTable.getItems().get(x).ctx;
            List<Genotype> genotypes = new ArrayList<>();
            genotypes.add(this.genotypeTable.getItems().get(x).makeGenotype());
            int y = x + 1;
            while (y < this.genotypeTable.getItems().size()) {
                VariantContext v2 = this.genotypeTable.getItems().get(y).ctx;
                // yes '!=' is enough
                if (v2 != v1)
                    break;
                genotypes.add(this.genotypeTable.getItems().get(y).makeGenotype());
                y++;
            }
            VariantContextBuilder vb = new VariantContextBuilder(v1);
            vb.genotypes(genotypes);
            w.add(vb.make());
            x = y;
        }
        w.close();
        return f;
    } catch (final Exception err) {
        JfxUtils.dialog().cause(err).show(owner);
        return null;
    } finally {
        CloserUtil.close(w);
    }
}

Example 52 with Genotype

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

the class VcfDoest method run.

private void run(final LineIterator lr, final PrintWriter pw) throws IOException {
    SortingCollection<TranscriptInfo> sorting = null;
    CloseableIterator<TranscriptInfo> iter2 = null;
    try {
        while (lr.hasNext()) {
            VcfIterator in = VCFUtils.createVcfIteratorFromLineIterator(lr, true);
            final VCFHeader header = in.getHeader();
            final Pedigree pedigree = Pedigree.newParser().parse(header);
            if (pedigree.isEmpty()) {
                throw new IOException("No pedigree found in header VCF header. use VcfInjectPedigree to add it");
            }
            final SortedSet<Pedigree.Person> individuals = new TreeSet<>();
            for (final Pedigree.Person individual : pedigree.getPersons()) {
                if (individual.isAffected() || individual.isUnaffected()) {
                    individuals.add(individual);
                }
            }
            boolean first = true;
            pw.println("# samples ( 0: unaffected 1:affected)");
            pw.print("population <- data.frame(family=c(");
            first = true;
            for (final Pedigree.Person person : individuals) {
                if (!first)
                    pw.print(",");
                pw.print("\"" + person.getFamily().getId() + "\"");
                first = false;
            }
            pw.print("),name=c(");
            first = true;
            for (final Pedigree.Person person : individuals) {
                if (!first)
                    pw.print(",");
                pw.print("\"" + person.getId() + "\"");
                first = false;
            }
            pw.print("),status=c(");
            first = true;
            for (final Pedigree.Person person : individuals) {
                if (!first)
                    pw.print(",");
                pw.print(person.isUnaffected() ? 0 : 1);
                first = false;
            }
            pw.println("))");
            sorting = SortingCollection.newInstance(TranscriptInfo.class, new TranscriptInfoCodec(), new TranscriptInfoCmp(), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
            sorting.setDestructiveIteration(true);
            final SAMSequenceDictionaryProgress progess = new SAMSequenceDictionaryProgress(header.getSequenceDictionary());
            /**
             * loop over variants
             */
            while (in.hasNext() && !pw.checkError()) {
                final VariantContext ctx = progess.watch(in.next());
                if (ctx.isFiltered())
                    continue;
                if (ctx.getAlternateAlleles().isEmpty())
                    continue;
                final Allele altAllele = ctx.getAltAlleleWithHighestAlleleCount();
                final MafCalculator mafCalculator = new MafCalculator(altAllele, ctx.getContig());
                boolean genotyped = false;
                for (final Pedigree.Person p : pedigree.getPersons()) {
                    if (!(p.isAffected() || p.isUnaffected()))
                        continue;
                    final Genotype g = ctx.getGenotype(p.getId());
                    if (g == null)
                        throw new IOException("Strange I cannot find individual " + p + " in the pedigree. Aborting.");
                    if (g.isCalled()) {
                        mafCalculator.add(g, p.isMale());
                    }
                    if (g.isHet() || g.isHomVar()) {
                        if (!g.getAlleles().contains(altAllele))
                            continue;
                        genotyped = true;
                        break;
                    }
                }
                if (!genotyped)
                    continue;
                final Interval interval = new Interval(ctx.getContig(), ctx.getStart(), ctx.getEnd());
                final List<KnownGene> genes = this.overlap(interval);
                if (genes.isEmpty())
                    continue;
                for (final KnownGene kg : genes) {
                    final TranscriptInfo trInfo = new TranscriptInfo();
                    trInfo.contig = kg.getContig();
                    trInfo.txStart = kg.getTxStart();
                    trInfo.txEnd = kg.getTxEnd();
                    trInfo.transcriptName = kg.getName();
                    trInfo.strand = (byte) (kg.isPositiveStrand() ? '+' : '-');
                    trInfo.exonCount = kg.getExonCount();
                    trInfo.transcriptLength = kg.getTranscriptLength();
                    trInfo.ctxStart = ctx.getStart();
                    trInfo.ref = ctx.getReference();
                    trInfo.alt = altAllele;
                    trInfo.maf = mafCalculator.getMaf();
                    trInfo.genotypes = new byte[individuals.size()];
                    int idx = 0;
                    for (final Pedigree.Person individual : individuals) {
                        final Genotype genotype = ctx.getGenotype(individual.getId());
                        final byte b;
                        if (genotype.isHomRef()) {
                            b = 0;
                        } else if (genotype.isHomVar() && genotype.getAlleles().contains(altAllele)) {
                            b = 2;
                        } else if (genotype.isHet() && genotype.getAlleles().contains(altAllele) && genotype.getAlleles().contains(ctx.getReference())) {
                            b = 1;
                        } else /* we treat 0/2 has hom-ref */
                        if (genotype.isHet() && !genotype.getAlleles().contains(altAllele) && genotype.getAlleles().contains(ctx.getReference())) {
                            LOG.warn("Treating " + genotype + " as hom-ref (0) alt=" + altAllele);
                            b = 0;
                        } else /* we treat 2/2 has hom-ref */
                        if (genotype.isHomVar() && !genotype.getAlleles().contains(altAllele)) {
                            LOG.warn("Treating " + genotype + " as hom-ref (0) alt=" + altAllele);
                            b = 0;
                        } else {
                            b = -9;
                        }
                        trInfo.genotypes[idx] = b;
                        ++idx;
                    }
                    KnownGene archetype = kg;
                    /* find gene archetype = longest overlapping */
                    for (final KnownGene kg2 : genes) {
                        if (kg2 == kg)
                            continue;
                        if (archetype.getStrand().equals(kg2.getStrand()) && archetype.getTranscriptLength() < kg2.getTranscriptLength()) {
                            archetype = kg2;
                        }
                    }
                    trInfo.archetypeName = archetype.getName();
                    trInfo.archetypeLength = archetype.getTranscriptLength();
                    boolean ctxWasFoundInExon = false;
                    final int ctxPos0 = ctx.getStart() - 1;
                    int indexInTranscript0 = 0;
                    for (final KnownGene.Exon exon : kg.getExons()) {
                        // variant in exon ?
                        if (!(exon.getStart() > (ctx.getEnd() - 1) || (ctx.getStart() - 1) >= exon.getEnd())) {
                            ctxWasFoundInExon = true;
                            indexInTranscript0 += (ctxPos0 - exon.getStart());
                            if (kg.isNegativeStrand()) {
                                indexInTranscript0 = (kg.getTranscriptLength() - 1) - indexInTranscript0;
                            }
                            trInfo.indexInTranscript0 = indexInTranscript0;
                            trInfo.overlapName = exon.getName();
                            sorting.add(trInfo);
                            break;
                        } else {
                            indexInTranscript0 += (exon.getEnd() - exon.getStart());
                        }
                    }
                    if (ctxWasFoundInExon) {
                        continue;
                    }
                    indexInTranscript0 = 0;
                    // search closest intron/exon junction
                    for (int ex = 0; ex + 1 < kg.getExonCount(); ++ex) {
                        final KnownGene.Exon exon1 = kg.getExon(ex);
                        indexInTranscript0 += (exon1.getEnd() - exon1.getStart());
                        final KnownGene.Exon exon2 = kg.getExon(ex + 1);
                        if (exon1.getEnd() <= ctxPos0 && ctxPos0 < exon2.getStart()) {
                            final int dist_to_exon1 = ctxPos0 - exon1.getEnd();
                            final int dist_to_exon2 = exon2.getStart() - ctxPos0;
                            if (dist_to_exon2 < dist_to_exon1) {
                                indexInTranscript0++;
                            }
                            if (kg.isNegativeStrand()) {
                                indexInTranscript0 = (kg.getTranscriptLength() - 1) - indexInTranscript0;
                            }
                            trInfo.indexInTranscript0 = indexInTranscript0;
                            trInfo.overlapName = exon1.getNextIntron().getName();
                            sorting.add(trInfo);
                            break;
                        }
                    }
                }
            // end loop over genes
            }
            // end while loop over variants
            progess.finish();
            sorting.doneAdding();
            LOG.info("done adding");
            iter2 = sorting.iterator();
            final EqualRangeIterator<TranscriptInfo> eqiter = new EqualRangeIterator<TranscriptInfo>(iter2, new Comparator<TranscriptInfo>() {

                @Override
                public int compare(final TranscriptInfo o1, final TranscriptInfo o2) {
                    int i = o1.contig.compareTo(o2.contig);
                    if (i != 0)
                        return i;
                    i = o1.transcriptName.compareTo(o2.transcriptName);
                    return i;
                }
            });
            while (eqiter.hasNext()) {
                final List<TranscriptInfo> list = eqiter.next();
                final TranscriptInfo front = list.get(0);
                pw.println("# BEGIN TRANSCRIPT " + front.transcriptName + " ##########################################");
                pw.println("transcript.chrom <- \"" + front.contig + "\"");
                pw.println("transcript.txStart0 <- " + front.txStart + "");
                pw.println("transcript.txEnd0 <- " + front.txEnd + "");
                pw.println("transcript.name <- \"" + front.transcriptName + "\"");
                pw.println("transcript.strand <- \"" + ((char) front.strand) + "\"");
                pw.println("transcript.length <- " + front.transcriptLength + "");
                pw.println("transcript.exonCount <- " + front.exonCount + "");
                pw.println("archetype.name <- \"" + front.archetypeName + "\"");
                pw.println("archetype.length <- " + front.archetypeLength + "");
                pw.print("variants <- data.frame(chrom=c(");
                first = true;
                for (final TranscriptInfo v : list) {
                    if (!first)
                        pw.print(",");
                    pw.print("\"" + v.contig + "\"");
                    first = false;
                }
                pw.print("),chromStart=c(");
                first = true;
                for (final TranscriptInfo v : list) {
                    if (!first)
                        pw.print(",");
                    pw.print(v.ctxStart);
                    first = false;
                }
                pw.print("),chromEnd=c(");
                first = true;
                for (final TranscriptInfo v : list) {
                    if (!first)
                        pw.print(",");
                    pw.print(v.ctxStart + v.ref.length() - 1);
                    first = false;
                }
                pw.print("),refAllele=c(");
                first = true;
                for (final TranscriptInfo v : list) {
                    if (!first)
                        pw.print(",");
                    pw.print("\"" + v.ref.getDisplayString() + "\"");
                    first = false;
                }
                pw.print("),altAllele=c(");
                first = true;
                for (final TranscriptInfo v : list) {
                    if (!first)
                        pw.print(",");
                    pw.print("\"" + v.alt.getDisplayString() + "\"");
                    first = false;
                }
                pw.print("),positionInTranscript1=c(");
                first = true;
                for (final TranscriptInfo v : list) {
                    if (!first)
                        pw.print(",");
                    pw.print(v.indexInTranscript0 + 1);
                    first = false;
                }
                pw.print("),maf=c(");
                first = true;
                for (final TranscriptInfo v : list) {
                    if (!first)
                        pw.print(",");
                    pw.print(v.maf);
                    first = false;
                }
                pw.print("),overlapName=c(");
                first = true;
                for (final TranscriptInfo v : list) {
                    if (!first)
                        pw.print(",");
                    pw.print("\"" + v.overlapName + "\"");
                    first = false;
                }
                pw.println("))");
                pw.println("# genotypes as a list. Should be a multiple of length(samples).");
                pw.println("# 0 is homref (0/0), 1 is het (0/1), 2 is homvar (1/1)");
                pw.println("# if the variant contains another ALT allele: (0/2) and (2/2) are considered 0 (homref)");
                pw.print("genotypes <- c(");
                first = true;
                for (final TranscriptInfo tr : list) {
                    for (byte g : tr.genotypes) {
                        if (!first)
                            pw.print(",");
                        first = false;
                        pw.print((int) g);
                    }
                }
                pw.println(")");
                pw.println("stopifnot(NROW(variants) * NROW(population) == length(genotypes) )");
                if (this.userDefinedFunName == null || this.userDefinedFunName.trim().isEmpty()) {
                    pw.println("## WARNING not user-defined R function was defined");
                } else {
                    pw.println("# consumme data with user-defined R function ");
                    pw.println(this.userDefinedFunName + "()");
                }
                pw.println("# END TRANSCRIPT " + front.transcriptName + " ##########################################");
            }
            // end while eqiter
            eqiter.close();
            iter2.close();
            iter2 = null;
            sorting.cleanup();
            sorting = null;
        }
    } finally {
        CloserUtil.close(iter2);
        if (sorting != null)
            sorting.cleanup();
    }
}

Example 53 with Genotype

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

the class GroupByGenotypes method map.

@Override
public Map<Category, Long> map(final RefMetaDataTracker tracker, final ReferenceContext refctx, final AlignmentContext context) {
    if (tracker == null)
        return Collections.emptyMap();
    final Map<Category, Long> counts = new HashMap<>();
    for (final VariantContext ctx : tracker.getValues(this.variants, context.getLocation())) {
        int index_singleton = -1;
        if (onlysingletons) {
            for (int i = 0; i < ctx.getNSamples(); ++i) {
                final Genotype g = ctx.getGenotype(i);
                if (g == null || !g.isCalled() || g.isNoCall() || g.isHomRef())
                    continue;
                if (index_singleton != -1) {
                    // not anymore a singleton
                    index_singleton = -1;
                    break;
                }
                index_singleton = i;
            }
        }
        for (int i = 0; i < ctx.getNSamples(); ++i) {
            if (onlysingletons && index_singleton != i) {
                continue;
            }
            final Genotype genotype = ctx.getGenotype(i);
            final List<Object> labels = new ArrayList<>();
            labels.add(genotype.getSampleName());
            if (bychrom)
                labels.add(ctx.getContig());
            if (byID)
                labels.add(ctx.hasID());
            if (byType)
                labels.add(ctx.getType().name());
            if (byGenotypeType)
                labels.add(genotype.getType());
            if (byFilter)
                labels.add(ctx.isFiltered());
            if (byGFilter)
                labels.add(genotype.isFiltered());
            if (minGenotypeQuality >= 0) {
                labels.add(genotype.hasGQ() && genotype.getGQ() >= this.minGenotypeQuality ? "." : "LOWQUAL");
            }
            if (byImpact) {
                AnnPredictionParser.Impact impact = null;
                for (final AnnPredictionParser.AnnPrediction pred : super.annParser.getPredictions(ctx)) {
                    // see http://stackoverflow.com/questions/41678374/
                    final Predicate<Allele> afilter = new Predicate<Allele>() {

                        @Override
                        public boolean test(final Allele A) {
                            return A.getDisplayString().equals(pred.getAllele());
                        }
                    };
                    if (genotype.getAlleles().stream().filter(afilter).findAny().isPresent() == false)
                        continue;
                    final AnnPredictionParser.Impact currImpact = pred.getPutativeImpact();
                    if (impact != null && currImpact.compareTo(impact) < 0)
                        continue;
                    impact = currImpact;
                }
                if (byImpact)
                    labels.add(impact == null ? "." : impact.name());
            }
            final Category cat = new Category(labels);
            Long n = counts.get(cat);
            counts.put(cat, n == null ? 1L : n + 1);
        }
    }
    return counts;
}

Example 54 with Genotype

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

the class MiniCaller method doWork.

@Override
public int doWork(final List<String> args) {
    ConcatSam.ConcatSamIterator iter = null;
    try {
        if (this.fastaFile == null) {
            LOG.error("no REF");
            return -1;
        }
        /* load faid */
        final ReferenceGenomeFactory referenceGenomeFactory = new ReferenceGenomeFactory();
        this.referenceGenome = referenceGenomeFactory.openFastaFile(this.fastaFile);
        this.dictionary = this.referenceGenome.getDictionary();
        if (this.dictionary == null) {
            LOG.error(JvarkitException.FastaDictionaryMissing.getMessage(this.fastaFile.getPath()));
        }
        /* create sam record iterator */
        iter = new ConcatSam.Factory().addInterval(this.rgnStr).setEnableUnrollList(true).open(args);
        final SAMFileHeader samFileheader = iter.getFileHeader();
        final SAMSequenceDictionary dict = samFileheader.getSequenceDictionary();
        if (dict == null) {
            LOG.error(JvarkitException.BamDictionaryMissing.getMessage(String.join(", ", args)));
            return -1;
        }
        if (!SequenceUtil.areSequenceDictionariesEqual(dict, this.dictionary)) {
            LOG.error(JvarkitException.DictionariesAreNotTheSame.getMessage(dict, this.dictionary));
            return -1;
        }
        final List<SAMReadGroupRecord> groups = samFileheader.getReadGroups();
        if (groups == null || groups.isEmpty()) {
            LOG.error("No group defined in input");
            return -1;
        }
        final Set<String> sampleSet = groups.stream().map(srgr -> this.samRecordPartition.apply(srgr, samRecordPartition.name())).collect(Collectors.toSet());
        /* create VCF metadata */
        final Set<VCFHeaderLine> metaData = new HashSet<VCFHeaderLine>();
        metaData.add(VCFStandardHeaderLines.getFormatLine(VCFConstants.GENOTYPE_KEY));
        metaData.add(VCFStandardHeaderLines.getFormatLine(VCFConstants.DEPTH_KEY));
        metaData.add(VCFStandardHeaderLines.getInfoLine(VCFConstants.DEPTH_KEY));
        metaData.add(VCFStandardHeaderLines.getInfoLine(VCFConstants.ALLELE_COUNT_KEY));
        metaData.add(VCFStandardHeaderLines.getInfoLine(VCFConstants.ALLELE_NUMBER_KEY));
        metaData.add(VCFStandardHeaderLines.getInfoLine(VCFConstants.ALLELE_FREQUENCY_KEY));
        metaData.add(new VCFFormatHeaderLine("DPG", // one value of each genotype
        VCFHeaderLineCount.G, VCFHeaderLineType.Integer, "Depth for each allele"));
        metaData.add(new VCFFormatHeaderLine("DP4", 4, VCFHeaderLineType.Integer, "Depth ReforAlt|Strand : RF,RR,AF,AR"));
        metaData.add(new VCFInfoHeaderLine("INDEL", 1, VCFHeaderLineType.Flag, "Variant is indel"));
        // addMetaData(metaData);
        final VCFHeader vcfHeader = new VCFHeader(metaData, sampleSet);
        vcfHeader.setSequenceDictionary(this.dictionary);
        /* create variant context */
        this.variantContextWriter = super.openVariantContextWriter(outputFile);
        this.variantContextWriter.writeHeader(vcfHeader);
        ReferenceContig genomicSeq = null;
        SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(this.dictionary);
        for (; ; ) {
            SAMRecord rec = null;
            if (iter.hasNext()) {
                rec = progress.watch(iter.next());
                if (rec.getReadUnmappedFlag())
                    continue;
                if (this.readFilter.filterOut(rec))
                    continue;
                /* flush buffer if needed */
                while (!this.buffer.isEmpty() && (this.buffer.get(0).tid < rec.getReferenceIndex() || (this.buffer.get(0).tid == rec.getReferenceIndex() && (this.buffer.get(0).getEnd()) < rec.getAlignmentStart()))) {
                    this.buffer.remove(0).print();
                }
                /* get genomic sequence at this position */
                if (genomicSeq == null || !genomicSeq.getContig().equals(rec.getContig())) {
                    genomicSeq = this.referenceGenome.getContig(rec.getContig());
                }
                final Cigar cigar = rec.getCigar();
                if (cigar == null)
                    continue;
                int readPos = 0;
                // 0 based-reference
                int refPos0 = rec.getAlignmentStart() - 1;
                final byte[] bases = rec.getReadBases();
                final byte[] quals = rec.getBaseQualities();
                final String sampleName = this.samRecordPartition.getPartion(rec, samRecordPartition.name());
                for (final CigarElement ce : cigar.getCigarElements()) {
                    final CigarOperator op = ce.getOperator();
                    switch(op) {
                        case P:
                            break;
                        case H:
                            break;
                        case S:
                            readPos += ce.getLength();
                            break;
                        // go
                        case N:
                        case D:
                            {
                                if (// we need base before deletion
                                refPos0 > 0) {
                                    char refBase = genomicSeq.charAt(refPos0 - 1);
                                    /* we use base before deletion */
                                    final StringBuilder sb = new StringBuilder(ce.getLength());
                                    sb.append(refBase);
                                    for (int i = 0; i < ce.getLength(); ++i) {
                                        sb.append(genomicSeq.charAt(refPos0 + i));
                                    }
                                    findContext(rec.getReferenceIndex(), // we use base *before deletion */
                                    refPos0 - 1, Allele.create(sb.toString(), true)).getSample(sampleName).getAllele(Allele.create(String.valueOf(refBase), false)).incr(rec.getReadNegativeStrandFlag());
                                }
                                refPos0 += ce.getLength();
                                break;
                            }
                        case I:
                            {
                                if (refPos0 > 0) {
                                    // float qual=0;
                                    char refBase = Character.toUpperCase(genomicSeq.charAt(refPos0 - 1));
                                    final StringBuilder sb = new StringBuilder(1 + ce.getLength());
                                    sb.append(refBase);
                                    for (int i = 0; i < ce.getLength(); ++i) {
                                        sb.append((char) bases[readPos + i]);
                                    // qual+=(readPos + i < quals.length?quals[ readPos + i]:0);
                                    }
                                    findContext(rec.getReferenceIndex(), // we use base *before deletion */
                                    refPos0 - 1, Allele.create(String.valueOf(refBase), true)).getSample(sampleName).getAllele(Allele.create(sb.toString().toUpperCase(), false)).incr(rec.getReadNegativeStrandFlag());
                                }
                                readPos += ce.getLength();
                                break;
                            }
                        case EQ:
                        case M:
                        case X:
                            {
                                for (int i = 0; i < ce.getLength(); ++i) {
                                    findContext(rec.getReferenceIndex(), refPos0 + i, Allele.create(String.valueOf(genomicSeq.charAt(refPos0 + i)), true)).getSample(sampleName).getAllele(Allele.create(String.valueOf((char) bases[readPos + i]), false)).incr(rec.getReadNegativeStrandFlag());
                                }
                                readPos += ce.getLength();
                                refPos0 += ce.getLength();
                                break;
                            }
                        default:
                            throw new IllegalStateException("Case statement didn't deal with cigar op: " + op);
                    }
                }
            } else {
                break;
            }
        }
        while (!buffer.isEmpty()) buffer.remove(0).print();
        progress.finish();
        iter.close();
        iter = null;
        this.variantContextWriter.close();
        this.variantContextWriter = null;
        return RETURN_OK;
    } catch (Exception e) {
        LOG.error(e);
        return -1;
    } finally {
        CloserUtil.close(iter);
        CloserUtil.close(this.referenceGenome);
        CloserUtil.close(this.variantContextWriter);
    }
}

Example 55 with Genotype

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

the class VcfRegistryCGI method doWork.

private void doWork(XMLStreamWriter w, final GroupFile gf) throws XMLStreamException {
    Position pos = parsePosition();
    if (pos == null)
        return;
    w.writeStartElement("div");
    w.writeStartElement("h2");
    w.writeCharacters(pos.chrom + ":" + pos.pos);
    w.writeEndElement();
    w.writeStartElement("table");
    w.writeStartElement("thead");
    w.writeStartElement("tr");
    for (String header : new String[] { "CHROM", "POS", "ID", "REF", "QUAL", "Sample", "Alleles", "DP", "GQ", "File" }) {
        w.writeStartElement("th");
        w.writeCharacters(header);
        // td
        w.writeEndElement();
    }
    // tr
    w.writeEndElement();
    // thead
    w.writeEndElement();
    w.writeStartElement("tbody");
    Set<String> samplesWithGenotypes = new HashSet<String>();
    Set<String> allSamples = new HashSet<String>();
    for (VcfFile f : getVcfFiles(gf)) {
        TabixReader tabixReader = null;
        TabixReader.Iterator iter = null;
        BlockCompressedInputStream bgzin = null;
        VCFHeader header = null;
        AbstractVCFCodec vcfCodec = VCFUtils.createDefaultVCFCodec();
        LineIterator lineIterator = null;
        for (int i = 0; i < 2; i++) {
            try {
                if (i == 0) {
                    bgzin = new BlockCompressedInputStream(f.file);
                    lineIterator = new LineIteratorImpl(new SynchronousLineReader(bgzin));
                    header = (VCFHeader) vcfCodec.readActualHeader(lineIterator);
                    allSamples.addAll(header.getGenotypeSamples());
                } else {
                    tabixReader = new TabixReader(f.file.getPath());
                    String line;
                    int[] x = tabixReader.parseReg(pos.chrom + ":" + pos.pos + "-" + (pos.pos));
                    if (x != null && x.length > 2 && x[0] != -1) {
                        iter = tabixReader.query(x[0], x[1], x[2]);
                    } else {
                    }
                    while (iter != null && (line = iter.next()) != null) {
                        VariantContext var = vcfCodec.decode(line);
                        for (String sample : header.getSampleNamesInOrder()) {
                            final Genotype genotype = var.getGenotype(sample);
                            if (genotype == null || !genotype.isCalled())
                                continue;
                            if (!genotype.isAvailable())
                                continue;
                            samplesWithGenotypes.add(sample);
                            w.writeStartElement("tr");
                            w.writeStartElement("td");
                            w.writeCharacters(var.getContig());
                            w.writeEndElement();
                            w.writeStartElement("td");
                            w.writeCharacters(String.valueOf(var.getStart()));
                            w.writeEndElement();
                            if (var.hasID()) {
                                w.writeStartElement("td");
                                if (var.getID().matches("rs[0-9]+")) {
                                    w.writeStartElement("a");
                                    w.writeAttribute("href", "http://www.ncbi.nlm.nih.gov/snp/" + var.getID().substring(2));
                                    w.writeCharacters(var.getID());
                                    // a
                                    w.writeEndElement();
                                } else {
                                    w.writeCharacters(var.getID());
                                }
                                // td
                                w.writeEndElement();
                            } else {
                                w.writeEmptyElement("td");
                            }
                            if (var.getReference() != null) {
                                w.writeStartElement("td");
                                w.writeCharacters(var.getReference().getBaseString());
                                w.writeEndElement();
                            } else {
                                w.writeEmptyElement("td");
                            }
                            if (var.hasLog10PError()) {
                                w.writeStartElement("td");
                                w.writeCharacters(String.valueOf((int) var.getPhredScaledQual()));
                                w.writeEndElement();
                            } else {
                                w.writeEmptyElement("td");
                            }
                            w.writeStartElement("td");
                            w.writeCharacters(sample);
                            w.writeEndElement();
                            List<Allele> alleles = genotype.getAlleles();
                            w.writeStartElement("td");
                            w.writeStartElement("span");
                            if (genotype.isHomRef()) {
                                w.writeAttribute("style", "color:green;");
                            } else if (genotype.isHomVar()) {
                                w.writeAttribute("style", "color:red;");
                            } else if (genotype.isHet()) {
                                w.writeAttribute("style", "color:blue;");
                            }
                            for (int j = 0; j < alleles.size(); ++j) {
                                if (j > 0)
                                    w.writeCharacters(genotype.isPhased() ? "|" : "/");
                                w.writeCharacters(alleles.get(j).getBaseString());
                            }
                            // span
                            w.writeEndElement();
                            w.writeEndElement();
                            if (genotype.hasDP()) {
                                w.writeStartElement("td");
                                w.writeCharacters(String.valueOf(genotype.getDP()));
                                w.writeEndElement();
                            } else {
                                w.writeEmptyElement("td");
                            }
                            if (genotype.hasGQ()) {
                                w.writeStartElement("td");
                                w.writeCharacters(String.valueOf(genotype.getGQ()));
                                w.writeEndElement();
                            } else {
                                w.writeEmptyElement("td");
                            }
                            w.writeStartElement("td");
                            w.writeCharacters(f.file.getName());
                            w.writeEndElement();
                            // tr
                            w.writeEndElement();
                            w.flush();
                        }
                    }
                }
            } catch (Exception err) {
                w.writeComment("BOUM " + err);
                header = null;
                lastException = err;
            } finally {
                CloserUtil.close(lineIterator);
                CloserUtil.close(bgzin);
                CloserUtil.close(tabixReader);
                CloserUtil.close(iter);
            }
            if (i == 0 && header == null)
                break;
        }
        w.flush();
    }
    // tbody
    w.writeEndElement();
    // table
    w.writeEndElement();
    allSamples.removeAll(samplesWithGenotypes);
    if (!allSamples.isEmpty()) {
        w.writeStartElement("h3");
        w.writeCharacters("Samples not found");
        w.writeEndElement();
        w.writeStartElement("ol");
        for (String sample : new TreeSet<String>(allSamples)) {
            w.writeStartElement("li");
            w.writeCharacters(sample);
            w.writeEndElement();
        }
        w.writeEndElement();
    }
    writeHTMLException(w, this.lastException);
    // div
    w.writeEndElement();
}