Search in sources :

Example 6 with AnnPredictionParserFactory

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

the class TestNg01 method testVcfFilterSo.

@Test
public void testVcfFilterSo() throws IOException {
    File output = new File(TEST_RESULTS_DIR, "jeter.filrerso.vcf");
    final AnnPredictionParser parser = new AnnPredictionParserFactory().createDefaultParser();
    final SequenceOntologyTree tree = SequenceOntologyTree.getInstance();
    String acn = "SO:0001583";
    final SequenceOntologyTree.Term term = tree.getTermByAcn(acn);
    final Set<SequenceOntologyTree.Term> terms = term.getAllDescendants();
    Assert.assertNotNull(term);
    Assert.assertTrue(terms.size() > 1);
    Assert.assertTrue(terms.contains(term));
    Assert.assertEquals(0, new VcfFilterSequenceOntology().instanceMain(new String[] { "-o", output.getPath(), "-A", acn, VCF01 }));
    streamVcf(output).forEach(V -> {
        // System.err.println(V.getAttribute("ANN")+" vs "+ terms);
        Assert.assertTrue(parser.getPredictions(V).stream().flatMap(P -> P.getSOTerms().stream()).anyMatch(T -> terms.contains(T)));
    });
    Assert.assertEquals(0, new VcfFilterSequenceOntology().instanceMain(new String[] { "-o", output.getPath(), "-A", acn, "--rmatt", "--invert", VCF01 }));
    streamVcf(output).forEach(V -> {
        Assert.assertFalse(parser.getPredictions(V).stream().flatMap(P -> P.getSOTerms().stream()).anyMatch(T -> terms.contains(T)));
    });
    Assert.assertEquals(0, new VcfFilterSequenceOntology().instanceMain(new String[] { "-o", output.getPath(), "-A", acn, "--rmatt", VCF01 }));
    Assert.assertTrue(streamVcf(output).findAny().isPresent());
    Assert.assertTrue(output.delete());
}
Also used : AnnPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParser) Arrays(java.util.Arrays) VCFFileReader(htsjdk.variant.vcf.VCFFileReader) VCFHeader(htsjdk.variant.vcf.VCFHeader) VcfFilterJdk(com.github.lindenb.jvarkit.tools.vcffilterjs.VcfFilterJdk) Test(org.testng.annotations.Test) VcfBurdenFisherH(com.github.lindenb.jvarkit.tools.burden.VcfBurdenFisherH) Vcf2Xml(com.github.lindenb.jvarkit.tools.vcf2xml.Vcf2Xml) VcfInjectPedigree(com.github.lindenb.jvarkit.tools.burden.VcfInjectPedigree) VcfToTable(com.github.lindenb.jvarkit.tools.misc.VcfToTable) ReferenceGenomeFactory(com.github.lindenb.jvarkit.util.bio.fasta.ReferenceGenomeFactory) VcfMultiToOneAllele(com.github.lindenb.jvarkit.tools.misc.VcfMultiToOneAllele) SAXParser(javax.xml.parsers.SAXParser) VcfNoCallToHomRef(com.github.lindenb.jvarkit.tools.misc.VcfNoCallToHomRef) VcfBurdenFisherV(com.github.lindenb.jvarkit.tools.burden.VcfBurdenFisherV) CloserUtil(htsjdk.samtools.util.CloserUtil) PrintWriter(java.io.PrintWriter) VcfBurdenFilterGenes(com.github.lindenb.jvarkit.tools.burden.VcfBurdenFilterGenes) Set(java.util.Set) PadEmptyFastq(com.github.lindenb.jvarkit.tools.misc.PadEmptyFastq) VcfMultiToOne(com.github.lindenb.jvarkit.tools.onesamplevcf.VcfMultiToOne) VcfOffsetsIndexFactory(com.github.lindenb.jvarkit.tools.vcflist.VcfOffsetsIndexFactory) Stream(java.util.stream.Stream) VcfBurdenFilterExac(com.github.lindenb.jvarkit.tools.burden.VcfBurdenFilterExac) VcfMoveFiltersToInfo(com.github.lindenb.jvarkit.tools.burden.VcfMoveFiltersToInfo) CloseableIterator(htsjdk.samtools.util.CloseableIterator) Bam2Raster(com.github.lindenb.jvarkit.tools.bam2graphics.Bam2Raster) SortVcfOnInfo(com.github.lindenb.jvarkit.tools.sortvcfonref.SortVcfOnInfo) TrapIndexer(com.github.lindenb.jvarkit.tools.trap.TrapIndexer) IterableAdapter(htsjdk.samtools.util.IterableAdapter) VcfTrap(com.github.lindenb.jvarkit.tools.trap.VcfTrap) ArrayList(java.util.ArrayList) FixVcfMissingGenotypes(com.github.lindenb.jvarkit.tools.misc.FixVcfMissingGenotypes) Assert(org.testng.Assert) SequenceOntologyTree(com.github.lindenb.jvarkit.util.so.SequenceOntologyTree) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) VcfRebase(com.github.lindenb.jvarkit.tools.vcfrebase.VcfRebase) StreamSupport(java.util.stream.StreamSupport) Properties(java.util.Properties) Files(java.nio.file.Files) VCFBigWig(com.github.lindenb.jvarkit.tools.vcfbigwig.VCFBigWig) VcfGnomad(com.github.lindenb.jvarkit.tools.gnomad.VcfGnomad) IOException(java.io.IOException) File(java.io.File) DefaultHandler(org.xml.sax.helpers.DefaultHandler) VcfSetSequenceDictionary(com.github.lindenb.jvarkit.tools.misc.VcfSetSequenceDictionary) VcfCreateDictionary(com.github.lindenb.jvarkit.tools.misc.VcfCreateDictionary) NgsFilesSummary(com.github.lindenb.jvarkit.tools.ngsfiles.NgsFilesSummary) VcfBurdenRscriptV(com.github.lindenb.jvarkit.tools.burden.VcfBurdenRscriptV) AnnPredictionParserFactory(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParserFactory) BufferedReader(java.io.BufferedReader) GroupByGene(com.github.lindenb.jvarkit.tools.groupbygene.GroupByGene) VCFFamilies(com.github.lindenb.jvarkit.tools.vcftrios.VCFFamilies) ReferenceGenome(com.github.lindenb.jvarkit.util.bio.fasta.ReferenceGenome) Algorithms(com.github.lindenb.jvarkit.util.Algorithms) Biostar84452(com.github.lindenb.jvarkit.tools.biostar.Biostar84452) VCFTrios(com.github.lindenb.jvarkit.tools.vcftrios.VCFTrios) URL(java.net.URL) LowResBam2Raster(com.github.lindenb.jvarkit.tools.bam2graphics.LowResBam2Raster) VCFBed(com.github.lindenb.jvarkit.tools.vcfbed.VCFBed) Random(java.util.Random) VcfToSql(com.github.lindenb.jvarkit.tools.vcf2sql.VcfToSql) MiniCaller(com.github.lindenb.jvarkit.tools.calling.MiniCaller) GoUtils(com.github.lindenb.jvarkit.tools.misc.GoUtils) FindAVariation(com.github.lindenb.jvarkit.tools.misc.FindAVariation) VcfXmlAmalgamation(com.github.lindenb.jvarkit.tools.vcfamalgation.VcfXmlAmalgamation) ImageIO(javax.imageio.ImageIO) BamToSql(com.github.lindenb.jvarkit.tools.misc.BamToSql) Gff2KnownGene(com.github.lindenb.jvarkit.tools.misc.Gff2KnownGene) VCFFixIndels(com.github.lindenb.jvarkit.tools.vcffixindels.VCFFixIndels) VCFStripAnnotations(com.github.lindenb.jvarkit.tools.vcfstripannot.VCFStripAnnotations) BufferedImage(java.awt.image.BufferedImage) Predicate(java.util.function.Predicate) BeforeClass(org.testng.annotations.BeforeClass) Collectors(java.util.stream.Collectors) IlluminaReadName(com.github.lindenb.jvarkit.tools.misc.IlluminaReadName) List(java.util.List) BackLocate(com.github.lindenb.jvarkit.tools.backlocate.BackLocate) VcfToSvg(com.github.lindenb.jvarkit.tools.misc.VcfToSvg) VariantContext(htsjdk.variant.variantcontext.VariantContext) Pattern(java.util.regex.Pattern) FastaSequenceReader(com.github.lindenb.jvarkit.util.bio.fasta.FastaSequenceReader) VcfList(com.github.lindenb.jvarkit.tools.vcflist.VcfList) SamReaderFactory(htsjdk.samtools.SamReaderFactory) VcfFilterNotInPedigree(com.github.lindenb.jvarkit.tools.burden.VcfFilterNotInPedigree) VcfFilterSequenceOntology(com.github.lindenb.jvarkit.tools.vcffilterso.VcfFilterSequenceOntology) FindAllCoverageAtPosition(com.github.lindenb.jvarkit.tools.misc.FindAllCoverageAtPosition) VcfToRdf(com.github.lindenb.jvarkit.tools.vcf2rdf.VcfToRdf) DataProvider(org.testng.annotations.DataProvider) AnnPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParser) SAXParserFactory(javax.xml.parsers.SAXParserFactory) VcfLoopOverGenes(com.github.lindenb.jvarkit.tools.burden.VcfLoopOverGenes) VcfToHilbert(com.github.lindenb.jvarkit.tools.hilbert.VcfToHilbert) VcfStats(com.github.lindenb.jvarkit.tools.vcfstats.VcfStats) VcfRemoveUnusedAlt(com.github.lindenb.jvarkit.tools.misc.VcfRemoveUnusedAlt) FileInputStream(java.io.FileInputStream) CaseControlCanvas(com.github.lindenb.jvarkit.tools.burden.CaseControlCanvas) VcfCompareCallers(com.github.lindenb.jvarkit.tools.vcfcmp.VcfCompareCallers) ReferenceContig(com.github.lindenb.jvarkit.util.bio.fasta.ReferenceContig) SamReader(htsjdk.samtools.SamReader) FastqShuffle(com.github.lindenb.jvarkit.tools.fastq.FastqShuffle) BamTile(com.github.lindenb.jvarkit.tools.misc.BamTile) VcfBurdenMAF(com.github.lindenb.jvarkit.tools.burden.VcfBurdenMAF) FastaSequence(com.github.lindenb.jvarkit.util.bio.fasta.FastaSequence) ConcatSam(com.github.lindenb.jvarkit.tools.misc.ConcatSam) Bam2Wig(com.github.lindenb.jvarkit.tools.bam2wig.Bam2Wig) FileReader(java.io.FileReader) VcfFilterSequenceOntology(com.github.lindenb.jvarkit.tools.vcffilterso.VcfFilterSequenceOntology) AnnPredictionParserFactory(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParserFactory) SequenceOntologyTree(com.github.lindenb.jvarkit.util.so.SequenceOntologyTree) File(java.io.File) Test(org.testng.annotations.Test)

Example 7 with AnnPredictionParserFactory

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

the class VcfBurdenGoEnrichment method doWork.

@Override
public int doWork(final List<String> args) {
    if (StringUtil.isBlank(this.readingGo.goUri)) {
        LOG.error("Undefined GOs uri.");
        return -1;
    }
    if (this.geneFile == null || !this.geneFile.exists()) {
        LOG.error("Undefined gene file option.");
        return -1;
    }
    try {
        final GoTree gotree = this.readingGo.createParser().setIgnoreDbXRef(true).parse(this.readingGo.goUri);
        List<GoTree.Term> terms = new ArrayList<>(gotree.getTerms());
        final Map<GoTree.Term, Node> term2node = new HashMap<>();
        // build the node TREE
        while (!terms.isEmpty()) {
            int i = 0;
            while (i < terms.size()) {
                final GoTree.Term t = terms.get(i);
                if (!t.hasRelations()) {
                    term2node.put(t, new Node(t));
                    terms.remove(i);
                } else if (t.getRelations().stream().allMatch(L -> term2node.containsKey(L.getTo()))) {
                    final Node n = new Node(t);
                    n.parents.addAll(t.getRelations().stream().map(L -> term2node.get(L.getTo())).collect(Collectors.toSet()));
                    term2node.put(t, n);
                    terms.remove(i);
                } else {
                    i++;
                }
            }
        }
        terms = null;
        final Set<String> unknownAcn = new HashSet<>();
        final Map<String, Set<Node>> gene2node = new HashMap<>();
        final BufferedReader r = IOUtils.openFileForBufferedReading(this.geneFile);
        String line;
        while ((line = r.readLine()) != null) {
            if (line.isEmpty() || line.startsWith("#"))
                continue;
            final int t = line.indexOf('\t');
            if (t == -1) {
                r.close();
                LOG.error("tab missing in " + line + " of " + this.geneFile);
                return -1;
            }
            final String gene = line.substring(0, t).trim();
            if (StringUtil.isBlank(gene)) {
                r.close();
                LOG.error("Emtpy gene in " + line);
                return -1;
            }
            // using getTermByName because found sysnonym in GOA
            final String termAcn = line.substring(t + 1).trim();
            if (unknownAcn.contains(termAcn))
                continue;
            final GoTree.Term term = gotree.getTermByName(termAcn);
            if (term == null && !unknownAcn.contains(termAcn)) {
                unknownAcn.add(termAcn);
                LOG.warning("Don't know this GO term in " + line + " of " + this.geneFile + ". Could be obsolete, synonym, go specific division. Skipping.");
                continue;
            }
            final Node node = term2node.get(term);
            if (node == null) {
                r.close();
                LOG.error("Don't know this node in " + line + " of " + this.geneFile);
                return -1;
            }
            Set<Node> nodes = gene2node.get(gene);
            if (nodes == null) {
                nodes = new HashSet<>();
                gene2node.put(gene, nodes);
            }
            node.numGenes++;
            nodes.add(node);
        }
        ;
        // clean up
        unknownAcn.clear();
        r.close();
        final VcfIterator iter = openVcfIterator(oneFileOrNull(args));
        final VCFHeader header = iter.getHeader();
        final VepPredictionParser vepParser = new VepPredictionParserFactory(header).get();
        final AnnPredictionParser annParser = new AnnPredictionParserFactory(header).get();
        final Set<Pedigree.Person> persons;
        if (this.pedFile != null) {
            final Pedigree pedigree = Pedigree.newParser().parse(this.pedFile);
            persons = new Pedigree.CaseControlExtractor().extract(header, pedigree);
        } else {
            persons = new Pedigree.CaseControlExtractor().extract(header);
        }
        final Set<Pedigree.Person> affected = persons.stream().filter(P -> P.isAffected()).collect(Collectors.toSet());
        final Set<Pedigree.Person> unaffected = persons.stream().filter(P -> P.isUnaffected()).collect(Collectors.toSet());
        if (affected.isEmpty()) {
            LOG.error("No Affected individual");
            return -1;
        }
        if (unaffected.isEmpty()) {
            LOG.error("No unaffected individual");
            return -1;
        }
        final List<String> lookColumns = Arrays.asList("CCDS", "Feature", "ENSP", "Gene", "HGNC", "HGNC_ID", "SYMBOL", "RefSeq");
        final Predicate<Genotype> isWildGenotype = G -> {
            if (G == null)
                return false;
            return G.isHomRef();
        };
        final Predicate<Genotype> isAltGenotype = G -> {
            if (G == null)
                return false;
            return G.isCalled() && !G.isHomRef();
        };
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(header).logger(LOG);
        while (iter.hasNext()) {
            final VariantContext ctx = progress.watch(iter.next());
            if (!this.variantFilter.test(ctx))
                continue;
            final Set<String> genes = new HashSet<>();
            for (final String predStr : ctx.getAttributeAsList(vepParser.getTag()).stream().map(O -> String.class.cast(O)).collect(Collectors.toList())) {
                final VepPredictionParser.VepPrediction pred = vepParser.parseOnePrediction(ctx, predStr);
                for (final String col : lookColumns) {
                    final String token = pred.getByCol(col);
                    if (!StringUtil.isBlank(token)) {
                        genes.add(token);
                    }
                }
            }
            for (final String predStr : ctx.getAttributeAsList(annParser.getTag()).stream().map(O -> String.class.cast(O)).collect(Collectors.toList())) {
                final AnnPredictionParser.AnnPrediction pred = annParser.parseOnePrediction(predStr);
                final String token = pred.getGeneName();
                if (!StringUtil.isBlank(token)) {
                    genes.add(token);
                }
            }
            if (genes.isEmpty())
                continue;
            final Set<Node> nodes = genes.stream().filter(G -> gene2node.containsKey(G)).flatMap(G -> gene2node.get(G).stream()).collect(Collectors.toSet());
            if (nodes.isEmpty())
                continue;
            final long unaffected_alt = unaffected.stream().map(P -> ctx.getGenotype(P.getId())).filter(G -> this.genotypeFilter.test(ctx, G)).filter(isAltGenotype).count();
            final long affected_alt = affected.stream().map(P -> ctx.getGenotype(P.getId())).filter(G -> this.genotypeFilter.test(ctx, G)).filter(isAltGenotype).count();
            /* no informative */
            if (unaffected_alt + affected_alt == 0L) {
                continue;
            }
            final long affected_ref = affected.stream().map(P -> ctx.getGenotype(P.getId())).filter(G -> this.genotypeFilter.test(ctx, G)).filter(isWildGenotype).count();
            final long unaffected_ref = unaffected.stream().map(P -> ctx.getGenotype(P.getId())).filter(G -> this.genotypeFilter.test(ctx, G)).filter(isWildGenotype).count();
            nodes.stream().forEach(N -> N.resetVisitedFlag());
            nodes.stream().forEach(N -> N.visit(unaffected_ref, unaffected_alt, affected_ref, affected_alt));
        }
        iter.close();
        progress.finish();
        LOG.info("Calculating Fisher and dumping.. please wait");
        final PrintWriter pw = super.openFileOrStdoutAsPrintWriter(this.outputFile);
        pw.println("#go_term\tfisher\tname\tgo_term_depth\tcount_genes_in_this_node" + "\tunaffected_ref_gt" + "\tunaffected_alt_gt" + "\taffected_ref_gt" + "\taffected_alt_gt");
        term2node.values().stream().filter(N -> this.show_never_seeen_term || N.sum() > 0L).sorted((n1, n2) -> Double.compare(n1.fisher(), n2.fisher())).forEach(N -> {
            pw.print(N.goTerm.getAcn());
            pw.print('\t');
            pw.print(N.fisher());
            pw.print("\t");
            pw.print(N.goTerm.getName().replaceAll("[ \',\\-]+", "_"));
            pw.print("\t");
            pw.print(N.goTerm.getMinDepth());
            pw.print('\t');
            pw.print(N.numGenes);
            pw.print('\t');
            pw.print(N.unaffected_ref);
            pw.print('\t');
            pw.print(N.unaffected_alt);
            pw.print('\t');
            pw.print(N.affected_ref);
            pw.print('\t');
            pw.print(N.affected_alt);
            pw.println();
        });
        pw.flush();
        pw.close();
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    }
}
Also used : Genotype(htsjdk.variant.variantcontext.Genotype) Arrays(java.util.Arrays) JexlVariantPredicate(com.github.lindenb.jvarkit.util.vcf.JexlVariantPredicate) Program(com.github.lindenb.jvarkit.util.jcommander.Program) Parameter(com.beust.jcommander.Parameter) VCFHeader(htsjdk.variant.vcf.VCFHeader) AnnPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParser) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) HashMap(java.util.HashMap) ParametersDelegate(com.beust.jcommander.ParametersDelegate) ArrayList(java.util.ArrayList) HashSet(java.util.HashSet) BiPredicate(java.util.function.BiPredicate) StringUtil(htsjdk.samtools.util.StringUtil) FisherExactTest(com.github.lindenb.jvarkit.math.stats.FisherExactTest) Pedigree(com.github.lindenb.jvarkit.util.Pedigree) Map(java.util.Map) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) VepPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.VepPredictionParser) VepPredictionParserFactory(com.github.lindenb.jvarkit.util.vcf.predictions.VepPredictionParserFactory) PrintWriter(java.io.PrintWriter) JexlGenotypePredicate(com.github.lindenb.jvarkit.util.vcf.JexlGenotypePredicate) GoTree(com.github.lindenb.jvarkit.util.go.GoTree) Predicate(java.util.function.Predicate) VcfIterator(com.github.lindenb.jvarkit.util.vcf.VcfIterator) Logger(com.github.lindenb.jvarkit.util.log.Logger) Set(java.util.Set) Collectors(java.util.stream.Collectors) File(java.io.File) List(java.util.List) AnnPredictionParserFactory(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParserFactory) VariantContext(htsjdk.variant.variantcontext.VariantContext) BufferedReader(java.io.BufferedReader) AnnPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParser) HashSet(java.util.HashSet) Set(java.util.Set) HashMap(java.util.HashMap) ArrayList(java.util.ArrayList) VariantContext(htsjdk.variant.variantcontext.VariantContext) VcfIterator(com.github.lindenb.jvarkit.util.vcf.VcfIterator) AnnPredictionParserFactory(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParserFactory) VCFHeader(htsjdk.variant.vcf.VCFHeader) HashSet(java.util.HashSet) PrintWriter(java.io.PrintWriter) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) Genotype(htsjdk.variant.variantcontext.Genotype) GoTree(com.github.lindenb.jvarkit.util.go.GoTree) VepPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.VepPredictionParser) Pedigree(com.github.lindenb.jvarkit.util.Pedigree) BufferedReader(java.io.BufferedReader) VepPredictionParserFactory(com.github.lindenb.jvarkit.util.vcf.predictions.VepPredictionParserFactory)

Aggregations

AnnPredictionParserFactory (com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParserFactory)7 VepPredictionParserFactory (com.github.lindenb.jvarkit.util.vcf.predictions.VepPredictionParserFactory)5 VariantContext (htsjdk.variant.variantcontext.VariantContext)5 VCFHeader (htsjdk.variant.vcf.VCFHeader)5 Set (java.util.Set)5 AnnPredictionParser (com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParser)4 File (java.io.File)4 PrintWriter (java.io.PrintWriter)4 ArrayList (java.util.ArrayList)4 HashSet (java.util.HashSet)4 List (java.util.List)4 Collectors (java.util.stream.Collectors)4 Parameter (com.beust.jcommander.Parameter)3 IOUtils (com.github.lindenb.jvarkit.io.IOUtils)3 Launcher (com.github.lindenb.jvarkit.util.jcommander.Launcher)3 Program (com.github.lindenb.jvarkit.util.jcommander.Program)3 Logger (com.github.lindenb.jvarkit.util.log.Logger)3 VepPredictionParser (com.github.lindenb.jvarkit.util.vcf.predictions.VepPredictionParser)3 CloserUtil (htsjdk.samtools.util.CloserUtil)3 IOException (java.io.IOException)3