Examples with GranthamScore com.github.lindenb.jvarkit.util.bio.GranthamScore used on opensource projects

Example 1 with GranthamScore

use of com.github.lindenb.jvarkit.util.bio.GranthamScore in project jvarkit by lindenb.

the class VCFCombineTwoSnvs method doVcfToVcf.

@Override
protected int doVcfToVcf(final String inputName, File saveAs) {
    BufferedReader bufferedReader = null;
    htsjdk.variant.variantcontext.writer.VariantContextWriter w = null;
    SortingCollection<CombinedMutation> mutations = null;
    CloseableIterator<Variant> varIter = null;
    CloseableIterator<CombinedMutation> mutIter = null;
    final Map<String, SamReader> sample2samReader = new HashMap<>();
    PrintWriter bedPeReport = null;
    try {
        bufferedReader = inputName == null ? IOUtils.openStreamForBufferedReader(stdin()) : IOUtils.openURIForBufferedReading(inputName);
        final VCFUtils.CodecAndHeader cah = VCFUtils.parseHeader(bufferedReader);
        /* get VCF header */
        final VCFHeader header = cah.header;
        final List<String> sampleList = header.getSampleNamesInOrder();
        this.indexedFastaSequenceFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(this.referencePath);
        final SAMSequenceDictionary dict = SequenceDictionaryUtils.extractRequired(this.indexedFastaSequenceFile);
        this.rnaSequenceFactory.setContigToGenomicSequence(C -> getGenomicSequenceForContig(C));
        if (this.bamIn != null) {
            final Set<String> sampleSet = new HashSet<>(sampleList);
            /**
             * unroll and open bam file
             */
            for (final Path bamFile : IOUtils.unrollPaths(Collections.singletonList(this.bamIn.toString()))) {
                LOG.info("opening BAM :" + this.bamIn);
                final SamReader samReader = SamReaderFactory.makeDefault().referenceSequence(this.referencePath).validationStringency(ValidationStringency.LENIENT).open(this.bamIn);
                if (!samReader.hasIndex()) {
                    samReader.close();
                    throw new IOException("Sam file is NOT indexed: " + bamFile);
                }
                final SAMFileHeader samHeader = samReader.getFileHeader();
                if (samHeader.getSequenceDictionary() == null || !SequenceUtil.areSequenceDictionariesEqual(dict, samReader.getFileHeader().getSequenceDictionary())) {
                    samReader.close();
                    throw new JvarkitException.DictionariesAreNotTheSame(dict, samReader.getFileHeader().getSequenceDictionary());
                }
                /* get sample name */
                String sampleName = null;
                for (final SAMReadGroupRecord rg : samHeader.getReadGroups()) {
                    if (rg.getSample() == null)
                        continue;
                    if (sampleName != null && !sampleName.equals(rg.getSample())) {
                        samReader.close();
                        throw new IOException(bamFile + " Contains two samples " + sampleName + " " + rg.getSample());
                    }
                    sampleName = rg.getSample();
                }
                if (sampleName == null) {
                    samReader.close();
                    LOG.warn("no sample in " + bamFile);
                    continue;
                }
                if (!sampleSet.contains(sampleName)) {
                    samReader.close();
                    LOG.warn("no sample " + sampleName + " in vcf. Ignoring " + bamFile);
                    continue;
                }
                sample2samReader.put(sampleName, samReader);
            }
        }
        loadTranscripts();
        this.variants = SortingCollection.newInstance(Variant.class, new VariantCodec(), new VariantComparatorTwo(dict), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
        this.variants.setDestructiveIteration(true);
        ProgressFactory.Watcher<VariantContext> progress1 = ProgressFactory.newInstance().dictionary(header).logger(LOG).build();
        String vcfLine = null;
        while ((vcfLine = bufferedReader.readLine()) != null) {
            final VariantContext ctx = progress1.apply(cah.codec.decode(vcfLine));
            /* discard non SNV variant */
            if (!ctx.isVariant() || ctx.isIndel()) {
                continue;
            }
            /* find the overlapping genes : extend the interval of the variant to include the stop codon */
            final Collection<Transcript> genes = this.knownGenes.getOverlapping(new Interval(ctx.getContig(), Math.max(1, ctx.getStart() - 3), ctx.getEnd() + 3)).stream().flatMap(L -> L.stream()).collect(Collectors.toList());
            final List<Allele> alternateAlleles = ctx.getAlternateAlleles();
            /* loop over overlapping genes */
            for (final Transcript kg : genes) {
                /* loop over available alleles */
                for (int allele_idx = 0; allele_idx < alternateAlleles.size(); ++allele_idx) {
                    final Allele alt = alternateAlleles.get(allele_idx);
                    challenge(ctx, alt, kg, vcfLine);
                }
            }
        }
        progress1.close();
        this.variants.doneAdding();
        bedPeReport = this.bedPePath == null ? new PrintWriter(new NullOuputStream()) : IOUtils.openPathForPrintWriter(this.bedPePath);
        mutations = SortingCollection.newInstance(CombinedMutation.class, new MutationCodec(), new MutationComparatorTwo(dict), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
        mutations.setDestructiveIteration(true);
        final VCFFilterHeaderLine vcfFilterHeaderLine = new VCFFilterHeaderLine("TwoHaplotypes", "(number of reads carrying both mutation) < (reads carrying variant 1 + reads carrying variant 2) ");
        varIter = this.variants.iterator();
        @SuppressWarnings("resource") EqualRangeIterator<Variant> eqVarIter = new EqualRangeIterator<>(varIter, new VariantComparatorOne(dict));
        ProgressFactory.Watcher<Variant> progress2 = ProgressFactory.newInstance().dictionary(header).logger(LOG).build();
        while (eqVarIter.hasNext()) {
            final List<Variant> buffer = eqVarIter.next();
            if (buffer.size() < 2)
                continue;
            for (int i = 0; i + 1 < buffer.size(); ++i) {
                final Variant v1 = buffer.get(i);
                for (int j = i + 1; j < buffer.size(); ++j) {
                    final Variant v2 = buffer.get(j);
                    if (v1.codonStart() != v2.codonStart())
                        continue;
                    if (v1.positionInCodon() == v2.positionInCodon())
                        continue;
                    if (!v1.wildCodon.equals(v2.wildCodon)) {
                        throw new IllegalStateException();
                    }
                    // no sample share the two variants
                    final Set<Integer> sharedSamplesIdx = v1.getSharedSampleIndexes(v2);
                    if (sharedSamplesIdx.isEmpty() && !sampleList.isEmpty())
                        continue;
                    final StringBuilder combinedCodon = new StringBuilder(v1.wildCodon);
                    combinedCodon.setCharAt(v1.positionInCodon(), v1.mutCodon.charAt(v1.positionInCodon()));
                    combinedCodon.setCharAt(v2.positionInCodon(), v2.mutCodon.charAt(v2.positionInCodon()));
                    final String pwild = PeptideSequence.of(v1.wildCodon).toString();
                    final String p1 = PeptideSequence.of(v1.mutCodon).toString();
                    final String p2 = PeptideSequence.of(v2.mutCodon).toString();
                    final String pCombined = PeptideSequence.of(combinedCodon).toString();
                    final String combinedSO;
                    final String combinedType;
                    /* both AA are synonymous, while combined is not */
                    if (!pCombined.equals(pwild) && p1.equals(pwild) && p2.equals(pwild)) {
                        combinedType = "combined_is_nonsynonymous";
                        if (pCombined.equals("*")) {
                            /* http://www.sequenceontology.org/browser/current_svn/term/SO:0001587 */
                            combinedSO = "stop_gained";
                        } else if (pwild.equals("*")) {
                            /* http://www.sequenceontology.org/browser/current_svn/term/SO:0002012 */
                            combinedSO = "stop_lost";
                        } else {
                            /* http://www.sequenceontology.org/miso/current_svn/term/SO:0001992 */
                            combinedSO = "nonsynonymous_variant";
                        }
                    } else if (!pCombined.equals(p1) && !pCombined.equals(p2) && !pCombined.equals(pwild)) {
                        combinedType = "combined_is_new";
                        if (pCombined.equals("*")) {
                            /* http://www.sequenceontology.org/browser/current_svn/term/SO:0001587 */
                            combinedSO = "stop_gained";
                        } else {
                            /* http://www.sequenceontology.org/miso/current_svn/term/SO:0001992 */
                            combinedSO = "nonsynonymous_variant";
                        }
                    } else {
                        combinedType = null;
                        combinedSO = null;
                    }
                    /**
                     * ok, there is something interesting here ,
                     * create two new Mutations carrying the
                     * two variants
                     */
                    if (combinedSO != null) {
                        /**
                         * grantham score is max found combined vs (p1/p2/wild)
                         */
                        int grantham_score = GranthamScore.score(pCombined.charAt(0), pwild.charAt(0));
                        grantham_score = Math.max(grantham_score, GranthamScore.score(pCombined.charAt(0), p1.charAt(0)));
                        grantham_score = Math.max(grantham_score, GranthamScore.score(pCombined.charAt(0), p2.charAt(0)));
                        /**
                         * info that will be displayed in the vcf
                         */
                        final Map<String, Object> info1 = v1.getInfo(v2);
                        final Map<String, Object> info2 = v2.getInfo(v1);
                        // filter for this combined: default it fails the filter
                        String filter = vcfFilterHeaderLine.getID();
                        final Map<String, Object> combinedMap = new LinkedHashMap<>();
                        combinedMap.put("CombinedCodon", combinedCodon);
                        combinedMap.put("CombinedAA", pCombined);
                        combinedMap.put("CombinedSO", combinedSO);
                        combinedMap.put("CombinedType", combinedType);
                        combinedMap.put("GranthamScore", grantham_score);
                        info1.putAll(combinedMap);
                        info2.putAll(combinedMap);
                        final Map<String, CoverageInfo> sample2coverageInfo = new HashMap<>(sample2samReader.size());
                        final int chromStart = Math.min(v1.genomicPosition1, v2.genomicPosition1);
                        final int chromEnd = Math.max(v1.genomicPosition1, v2.genomicPosition1);
                        /* get phasing info for each sample*/
                        for (final String sampleName : sample2samReader.keySet()) {
                            final SamReader samReader = sample2samReader.get(sampleName);
                            final CoverageInfo covInfo = new CoverageInfo();
                            sample2coverageInfo.put(sampleName, covInfo);
                            SAMRecordIterator iter = null;
                            try {
                                iter = samReader.query(v1.contig, chromStart, chromEnd, false);
                                while (iter.hasNext()) {
                                    final SAMRecord rec = iter.next();
                                    if (rec.getReadUnmappedFlag())
                                        continue;
                                    if (rec.isSecondaryOrSupplementary())
                                        continue;
                                    if (rec.getDuplicateReadFlag())
                                        continue;
                                    if (rec.getReadFailsVendorQualityCheckFlag())
                                        continue;
                                    // get DEPTh for variant 1
                                    if (rec.getAlignmentStart() <= v1.genomicPosition1 && v1.genomicPosition1 <= rec.getAlignmentEnd()) {
                                        covInfo.depth1++;
                                    }
                                    // get DEPTh for variant 2
                                    if (rec.getAlignmentStart() <= v2.genomicPosition1 && v2.genomicPosition1 <= rec.getAlignmentEnd()) {
                                        covInfo.depth2++;
                                    }
                                    if (rec.getAlignmentEnd() < chromEnd)
                                        continue;
                                    if (rec.getAlignmentStart() > chromStart)
                                        continue;
                                    final Cigar cigar = rec.getCigar();
                                    if (cigar == null)
                                        continue;
                                    final byte[] bases = rec.getReadBases();
                                    if (bases == null)
                                        continue;
                                    int refpos1 = rec.getAlignmentStart();
                                    int readpos = 0;
                                    boolean found_variant1_on_this_read = false;
                                    boolean found_variant2_on_this_read = false;
                                    /**
                                     * loop over cigar
                                     */
                                    for (final CigarElement ce : cigar.getCigarElements()) {
                                        final CigarOperator op = ce.getOperator();
                                        switch(op) {
                                            case P:
                                                continue;
                                            case S:
                                            case I:
                                                readpos += ce.getLength();
                                                break;
                                            case D:
                                            case N:
                                                refpos1 += ce.getLength();
                                                break;
                                            case H:
                                                continue;
                                            case EQ:
                                            case M:
                                            case X:
                                                for (int x = 0; x < ce.getLength(); ++x) {
                                                    if (refpos1 == v1.genomicPosition1 && same(bases[readpos], v1.altAllele)) {
                                                        found_variant1_on_this_read = true;
                                                    } else if (refpos1 == v2.genomicPosition1 && same(bases[readpos], v2.altAllele)) {
                                                        found_variant2_on_this_read = true;
                                                    }
                                                    refpos1++;
                                                    readpos++;
                                                }
                                                break;
                                            default:
                                                throw new IllegalStateException(op.name());
                                        }
                                        /* skip remaining bases after last variant */
                                        if (refpos1 > chromEnd)
                                            break;
                                    }
                                    /* sum-up what we found */
                                    if (found_variant1_on_this_read && found_variant2_on_this_read) {
                                        covInfo.count_reads_having_both_variants++;
                                    } else if (!found_variant1_on_this_read && !found_variant2_on_this_read) {
                                        covInfo.count_reads_having_no_variants++;
                                    } else if (found_variant1_on_this_read) {
                                        covInfo.count_reads_having_variant1++;
                                    } else if (found_variant2_on_this_read) {
                                        covInfo.count_reads_having_variant2++;
                                    }
                                }
                            /* end of loop over reads */
                            } finally {
                                iter.close();
                                iter = null;
                            }
                            info1.put("N_READS_BOTH_VARIANTS_" + sampleName, covInfo.count_reads_having_both_variants);
                            info2.put("N_READS_BOTH_VARIANTS_" + sampleName, covInfo.count_reads_having_both_variants);
                            info1.put("N_READS_NO_VARIANTS_" + sampleName, covInfo.count_reads_having_no_variants);
                            info2.put("N_READS_NO_VARIANTS_" + sampleName, covInfo.count_reads_having_no_variants);
                            info1.put("N_READS_TOTAL_" + sampleName, covInfo.count_reads_having_both_variants + covInfo.count_reads_having_no_variants + covInfo.count_reads_having_variant1 + covInfo.count_reads_having_variant2);
                            info2.put("N_READS_TOTAL_" + sampleName, covInfo.count_reads_having_both_variants + covInfo.count_reads_having_no_variants + covInfo.count_reads_having_variant1 + covInfo.count_reads_having_variant2);
                            // count for variant 1
                            info1.put("N_READS_ONLY_1_" + sampleName, covInfo.count_reads_having_variant1);
                            info1.put("N_READS_ONLY_2_" + sampleName, covInfo.count_reads_having_variant2);
                            info1.put("DEPTH_1_" + sampleName, covInfo.depth1);
                            // inverse previous count
                            info2.put("N_READS_ONLY_1_" + sampleName, covInfo.count_reads_having_variant2);
                            info2.put("N_READS_ONLY_2_" + sampleName, covInfo.count_reads_having_variant1);
                            info2.put("DEPTH_2_" + sampleName, covInfo.depth2);
                            /* number of reads with both variant is greater than
								 * reads carrying only one variant: reset the filter 
								 */
                            if (2 * covInfo.count_reads_having_both_variants > (covInfo.count_reads_having_variant1 + covInfo.count_reads_having_variant2)) {
                                /* reset filter */
                                filter = VCFConstants.UNFILTERED;
                                info1.put("FILTER_1_" + sampleName, ".");
                                info2.put("FILTER_2_" + sampleName, ".");
                            } else {
                                info1.put("FILTER_1_" + sampleName, vcfFilterHeaderLine.getID());
                                info2.put("FILTER_2_" + sampleName, vcfFilterHeaderLine.getID());
                            }
                        }
                        /* end of loop over bams */
                        final CombinedMutation m1 = new CombinedMutation();
                        m1.contig = v1.contig;
                        m1.genomicPosition1 = v1.genomicPosition1;
                        m1.id = v1.id;
                        m1.refAllele = v1.refAllele;
                        m1.altAllele = v1.altAllele;
                        m1.vcfLine = v1.vcfLine;
                        m1.info = mapToString(info1);
                        m1.filter = filter;
                        m1.grantham_score = grantham_score;
                        m1.sampleIndexes.addAll(sharedSamplesIdx);
                        m1.sorting_id = ID_GENERATOR++;
                        mutations.add(m1);
                        final CombinedMutation m2 = new CombinedMutation();
                        m2.contig = v2.contig;
                        m2.genomicPosition1 = v2.genomicPosition1;
                        m2.id = v2.id;
                        m2.refAllele = v2.refAllele;
                        m2.altAllele = v2.altAllele;
                        m2.vcfLine = v2.vcfLine;
                        m2.info = mapToString(info2);
                        m2.filter = filter;
                        m2.grantham_score = grantham_score;
                        m2.sampleIndexes.addAll(sharedSamplesIdx);
                        m2.sorting_id = ID_GENERATOR++;
                        mutations.add(m2);
                        bedPeReport.print(m1.contig);
                        bedPeReport.print('\t');
                        bedPeReport.print(m1.genomicPosition1 - 1);
                        bedPeReport.print('\t');
                        bedPeReport.print(m1.genomicPosition1);
                        bedPeReport.print('\t');
                        bedPeReport.print(m2.contig);
                        bedPeReport.print('\t');
                        bedPeReport.print(m2.genomicPosition1 - 1);
                        bedPeReport.print('\t');
                        bedPeReport.print(m2.genomicPosition1);
                        bedPeReport.print('\t');
                        // name
                        bedPeReport.print(v1.transcriptId);
                        bedPeReport.print('\t');
                        // score
                        bedPeReport.print(grantham_score == GranthamScore.getDefaultScore() ? 0 : (int) ((grantham_score / 255.0) * 1000.0));
                        bedPeReport.print('\t');
                        final Transcript kg = this.knownGenes.getOverlapping(new Interval(v1.getContig(), v1.genomicPosition1 - 1, v1.genomicPosition1 + 1)).stream().flatMap(L -> L.stream()).filter(P -> P.getContig().equals(v1.contig) && P.getId().equals(v1.transcriptId)).findFirst().orElseThrow(IllegalStateException::new);
                        // strand1
                        bedPeReport.print(kg.isNegativeStrand() ? "-" : "+");
                        bedPeReport.print('\t');
                        // strand2
                        bedPeReport.print(kg.isNegativeStrand() ? "-" : "+");
                        bedPeReport.print('\t');
                        if (sharedSamplesIdx.isEmpty()) {
                            bedPeReport.print('.');
                        } else {
                            bedPeReport.print(sharedSamplesIdx.stream().map(I -> sampleList.get(I.intValue())).collect(Collectors.joining(";")));
                        }
                        bedPeReport.print('\t');
                        bedPeReport.print(combinedSO);
                        bedPeReport.print('\t');
                        bedPeReport.print(String.join(":", pwild, p1, p2, pCombined));
                        bedPeReport.println();
                    }
                }
            }
        }
        progress2.close();
        mutations.doneAdding();
        eqVarIter.close();
        eqVarIter = null;
        varIter.close();
        varIter = null;
        variants.cleanup();
        variants = null;
        bedPeReport.flush();
        bedPeReport.close();
        bedPeReport = null;
        final VCFHeader header2 = new VCFHeader(header);
        header2.addMetaDataLine(new VCFHeaderLine(getProgramName() + "AboutQUAL", "QUAL is filled with Grantham Score  http://www.ncbi.nlm.nih.gov/pubmed/4843792"));
        final StringBuilder infoDesc = new StringBuilder("Variant affected by two distinct mutation. Format is defined in the INFO column. ");
        final VCFInfoHeaderLine CodonVariantHeader = new VCFInfoHeaderLine("CodonVariant", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, infoDesc.toString());
        header2.addMetaDataLine(CodonVariantHeader);
        final VCFInfoHeaderLine CodonSampleHeader = new VCFInfoHeaderLine("Samples", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, "Samples that could be affected");
        header2.addMetaDataLine(CodonSampleHeader);
        JVarkitVersion.getInstance().addMetaData(this, header2);
        if (!sample2samReader.isEmpty()) {
            header2.addMetaDataLine(vcfFilterHeaderLine);
        }
        w = this.writingVariantsDelegate.dictionary(dict).open(IOUtil.toPath(saveAs));
        w.writeHeader(header2);
        ProgressFactory.Watcher<CombinedMutation> progress3 = ProgressFactory.newInstance().dictionary(header).logger(LOG).build();
        mutIter = mutations.iterator();
        EqualRangeIterator<CombinedMutation> eqRangeMutIter = new EqualRangeIterator<>(mutIter, new MutationComparatorOne(dict));
        while (eqRangeMutIter.hasNext()) {
            final List<CombinedMutation> mBuffer = eqRangeMutIter.next();
            if (mBuffer.isEmpty())
                break;
            progress3.apply(mBuffer.get(0));
            // default grantham score used in QUAL
            int grantham_score = -1;
            // default filter fails
            String filter = vcfFilterHeaderLine.getID();
            final CombinedMutation first = mBuffer.get(0);
            final Set<String> info = new HashSet<>();
            final VariantContext ctx = cah.codec.decode(first.vcfLine);
            final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
            vcb.chr(first.contig);
            vcb.start(first.genomicPosition1);
            vcb.stop(first.genomicPosition1 + first.refAllele.length() - 1);
            if (!first.id.equals(VCFConstants.EMPTY_ID_FIELD))
                vcb.id(first.id);
            for (final CombinedMutation m : mBuffer) {
                info.add(m.info);
                grantham_score = Math.max(grantham_score, m.grantham_score);
                if (VCFConstants.UNFILTERED.equals(m.filter)) {
                    // at least one SNP is ok one this line
                    filter = null;
                }
            }
            if (!sampleList.isEmpty()) {
                vcb.attribute(CodonSampleHeader.getID(), new ArrayList<>(mBuffer.stream().flatMap(S -> S.sampleIndexes.stream()).map(IDX -> sampleList.get(IDX)).collect(Collectors.toSet())));
            }
            vcb.unfiltered();
            if (filter != null && !sample2samReader.isEmpty()) {
                vcb.filter(filter);
            } else {
                vcb.passFilters();
            }
            vcb.attribute(CodonVariantHeader.getID(), new ArrayList<String>(info));
            if (grantham_score > 0) {
                vcb.log10PError(grantham_score / -10.0);
            } else {
                vcb.log10PError(VariantContext.NO_LOG10_PERROR);
            }
            w.add(vcb.make());
        }
        progress3.close();
        eqRangeMutIter.close();
        mutIter.close();
        mutations.cleanup();
        mutations = null;
        return RETURN_OK;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(this.indexedFastaSequenceFile);
        CloserUtil.close(mutIter);
        CloserUtil.close(varIter);
        CloserUtil.close(bedPeReport);
        if (this.variants != null)
            this.variants.cleanup();
        if (mutations != null)
            mutations.cleanup();
        this.variants = null;
        for (SamReader r : sample2samReader.values()) CloserUtil.close(r);
        CloserUtil.close(w);
        CloserUtil.close(bufferedReader);
    }
}