<jupyter_start><jupyter_text>A notebook to compute enrichment analysis using a SPARQL endpoint<jupyter_code>import sys
import rdflib
from IPython.core.display import display, HTML
from SPARQLWrapper import SPARQLWrapper, JSON, XML
import scipy.stats as ss
from decimal import Decimal
import pandas as pd, io
from pandas.io.json import json_normalize
pd.set_option("display.max_colwidth",300)
pd.set_option('colheader_justify', 'left')

def getPrefixDec(prefixes):
    l = ""
    for k,v in prefixes.items():
        l = l + "PREFIX " + k + ": <" + v + ">" + "\r\n"
    return l

def getValuesDec(entities):
    l = ""
    for i in entities:
        if(i[0:4] == "http"):
            l = l + "(<" + i + ">) "
        else:
            l = l + "(" + i + ") "
    return l

def getPopulationCount(endpoint, prefixes, triplepattern):
    prefixDec = getPrefixDec(prefixes)
    sparql = SPARQLWrapper(endpoint)
    query = prefixDec + "SELECT (count(*) AS ?c) {" + triplepattern + " }"
    sparql.setQuery(query)
    sparql.setReturnFormat(JSON)
    results = sparql.query().convert()
    count = int(results["results"]["bindings"][0]["c"]["value"])
    return count

def getProperties(endpoint, prefixes, entities = "", triplepattern = ""):
    prefixDec = getPrefixDec(prefixes)
    eValuesDec = ""
    if len(entities) != 0:
        eValuesDec = " VALUES (?s) {" + getValuesDec(entities) + "}"
        
    sparql = SPARQLWrapper(endpoint)
    query = prefixDec + """ 
    SELECT ?p ?plabel ?c
    {
     {{
        SELECT ?p (count(distinct ?s) AS ?c)
        { """ + eValuesDec  + """
          """ + triplepattern + """
          ?s ?p ?o .
        } GROUP BY ?p
        ORDER BY DESC(?c)
     }}
     OPTIONAL {
        ?p dct:title ?plabel
     }
    }
    """
    #print(query)
    sparql.setQuery(query)
    sparql.setReturnFormat(JSON)
    results = sparql.query().convert()
    return results["results"]["bindings"]

def getFrequencyValuesForSelectedEntitiesAndPredicates(endpoint, prefixes, entities = "", predicates = "", triplepattern = ""):
    prefixDec = getPrefixDec(prefixes)
    eValuesDec = " VALUES (?s) {" + getValuesDec(entities) + "}"
    pValuesDec = " VALUES (?p) {" + getValuesDec(predicates) + "}"
    sparql = SPARQLWrapper(endpoint)
    query = prefixDec + """ 
    SELECT ?p ?o ?olabel ?sc (count(?o) AS ?pc) 
    {
        {{ 
          SELECT ?p ?o (count(?o) AS ?sc) 
          { """ + eValuesDec + pValuesDec + """
            ?s ?p ?o 
          } GROUP BY ?p ?o
        }}
        """ + triplepattern + """
        ?s ?p ?o .
        OPTIONAL {
            ?o dct:title ?olabel
        }
    }
    """
    #print(query)
    
    sparql.setQuery(query)
    sparql.setReturnFormat(JSON)
    results = sparql.query().convert()
    return results["results"]["bindings"]

def makeCURIE(uri,prefixes):
    l = ""
    for prefix,uribase in prefixes.items():
        # now match substr
        size = len(uribase)
        if uribase == uri[0:size]:
            # match
            l = prefix + ":" + uri[size:]
    if l == "":
        l = uri
    return l

def performStatistics(nSamples, nPopulation, fv):
    ret = dict()
    
    meta = dict()
    meta["nSamples"] = nSamples
    meta["nPopulation"] = nPopulation    
    ret["meta"] = meta
    
    results = []
    for i in fv:
        o = dict()
        o["predicate"] = str(i["p"]["value"])
        o["attribute"] = str(i["o"]["value"])
        o['attribute_label'] = str(i["olabel"]["value"])
        o["sample_count"] = int(i["sc"]["value"])
        o["population_count"] = int(i["pc"]["value"])
    
        hpd = ss.hypergeom(nPopulation, o["population_count"], nSamples)
        prob = hpd.pmf(o["sample_count"])
        o["prob"]= prob
        results.append(o)
    ret["results"] = results
    return ret

def printResults(results,prefixes, pfilter = 0.001):
    meta = results["meta"]
    
    print("Sample     size: " + str(meta['nSamples']))
    print("Population size: " + str(meta['nPopulation']))
    for i in results["results"]:
        p = makeCURIE( i['predicate'], prefixes)
        o = makeCURIE( i['attribute'], prefixes)
        ol = i['attribute_label']
        
        if i['prob'] <= pfilter:
            if i['prob'] < 0.0001:
                prob = '{0:.2E}'.format(Decimal(i['prob']))
            else:
                prob = '{0:.5f}'.format(Decimal(i['prob']))
            print(" " + str(i['sample_count']) + " / " + str(i['population_count']) + " p-value: " + str(prob) + " " + str(p) + " " + str(o) + " " + str(ol))

#print(getPrefixDec({ "drugbank":"http://bio2rdf.org/drugbank:","dv":"http://bio2rdf.org/drugbank_vocabulary:"}))
#print(getValuesDec( ["http://bio2rdf.org/drugbank:test", "drugbank:test2"]))
#print(makeCURIE("http://bio2rdf.org/drugbank_vocabulary:category",prefixes))
### drug example
endpoint = "http://bio2rdf.org/sparql"
prefixes = { "dct":"http://purl.org/dc/terms/", "drugbank":"http://bio2rdf.org/drugbank:","dv":"http://bio2rdf.org/drugbank_vocabulary:"}
sample_names  = ["Eletriptan","Zolmitriptan","Dihydroergotamine","Almotriptan","Rizatriptan"]
sample_curies = ["drugbank:DB00216","drugbank:DB00315","drugbank:DB00320","drugbank:DB00918","drugbank:DB00953"]
population_tp = "?s rdf:type dv:Drug ."
attributes    = ["dv:category","dv:group"]

nSamples = len(sample_curies)
nPopulation = getPopulationCount(endpoint, prefixes, population_tp)
print("There are " + str(nSamples) + " samples in a population of " + str(nPopulation))

#fv_test = getProperties(endpoint, prefixes, "", population_tp)
#table = json_normalize(fv_test)
#table
#table[['p.value','plabel.value','c.value']]

fv = getFrequencyValuesForSelectedEntitiesAndPredicates(endpoint, prefixes, sample_curies, attributes, population_tp)
results = performStatistics(nSamples, nPopulation, fv)
printResults(results, prefixes)

endpoint = "http://bio2rdf.org/sparql"
prefixes = { "dct":"http://purl.org/dc/terms/", "sgd":"http://bio2rdf.org/sgd:","sgd_resource":"http://bio2rdf.org/sgd_resource:","sv":"http://bio2rdf.org/sgd_vocabulary:"}
sample_curies = ["sgd_resource:S000004425gp","sgd_resource:S000005376gp","sgd_resource:S000004238gp","sgd_resource:S000003399gp","sgd_resource:S000005853gp"]
population_tp = "?s rdf:type sv:Protein ."
attributes    = ["sv:function"]

nSamples = len(sample_curies)
nPopulation = getPopulationCount(endpoint, prefixes, population_tp)
fv = getFrequencyValuesForSelectedEntitiesAndPredicates(endpoint, prefixes, sample_curies, attributes, population_tp)
results = performStatistics(nSamples, nPopulation, fv)
printResults(results, prefixes)<jupyter_output>Sample     size: 5
Population size: 18303
 5 / 12 p-value: 4.63E-17 sv:function http://bio2rdf.org/go:0005516 calmodulin binding
 1 / 3 p-value: 0.00082 sv:function http://bio2rdf.org/go:0004723 calcium-dependent protein serine/threonine phosphatase activity
 3 / 622 p-value: 0.00036 sv:function http://bio2rdf.org/go:0005524 ATP binding
 2 / 122 p-value: 0.00043 sv:function http://bio2rdf.org/go:0004674 protein serine/threonine kinase activity
 2 / 128 p-value: 0.00048 sv:function http://bio2rdf.org/go:0004672 protein kinase activity
 2 / 141 p-value: 0.00058 sv:function http://bio2rdf.org/go:0016772 transferase activity, transferring phosphorus-containing groups
 1 / 2 p-value: 0.00055 sv:function http://bio2rdf.org/go:0004683 calmodulin-dependent protein kinase activity
 3 / 787 p-value: 0.00073 sv:function http://bio2rdf.org/go:0000166 nucleotide binding
{'fraction_non_alphanumeric': 0.16144158628081456, 'fraction_numerical': 0.033092175777063235, 'mean_word_length': 3.443452380952381, 'pattern_counts': {'":': 12, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '147674', 'license': 'MIT', 'n_tokens_mistral': 2764, 'n_tokens_neox': 2529, 'n_words': 492, 'path': 'notebooks/ldenrich.ipynb', 'repo_name': 'micheldumontier/translator-notebooks', 'size': '7464'}