annotacell.py

import pandas as pd
import numpy as np
import os, requests, json, csv
import psycopg2
import os
from openai import OpenAI
from tavily import TavilyClient
from dotenv import load_dotenv

#Process the Data
def fetch_expression_data(cluster, ensembl_ids):
    payload = {
        "filter": {
            "dataset_ids": [],
            "development_stage_ontology_term_ids": [],
            "disease_ontology_term_ids": [],
            "gene_ontology_term_ids": ensembl_ids,
            "organism_ontology_term_id": "NCBITaxon:9606",
            "self_reported_ethnicity_ontology_term_ids": [],
            "sex_ontology_term_ids": [],
            "publication_citations": [],
        },
        "is_rollup": True
    }

    # URL for the POST request
    API_URL = "https://api.cellxgene.cziscience.com/wmg/v2/query"

    try:
        response = requests.post(API_URL, json=payload)
        response.raise_for_status()  # Raise an exception for HTTP errors
        return response.json()['expression_summary']
    except requests.RequestException as e:
        print(f"Error fetching data for Cluster {cluster}: {e}")
        return None
    
def expression_data_to_df(data):
    flattened_data = []
    for gene_id, anatomical_structures in data.items():
        for anatomical_id, cell_types in anatomical_structures.items():
            for cell_type_id, aggregated_data in cell_types.items():
                metrics = aggregated_data['aggregated']
                flattened_data.append({
                    'gene': gene_id,
                    'tissue': anatomical_id,
                    'cell': cell_type_id,
                    'expression': metrics.get('me', None),
                    'cell count': metrics.get('n', None),
                    'cell percentage': metrics.get('pc', None),
                    'tissue composition': metrics.get('tpc', None)
                })
    return pd.DataFrame(flattened_data)

def filter_expression_data(response_df, TARGET_UBERON_ID):
    tissue_df = response_df[response_df['tissue'] == TARGET_UBERON_ID].drop(columns=['tissue'])
    return tissue_df[~tissue_df['cell'].isin(['tissue_stats', 'CL:0000000'])]

def translate_ontology(df):
    df['cell'] = df['cell'].apply(lambda x: ontl.cell(x))
    df['gene'] = df['gene'].apply(lambda x: ontl.gene(x)).fillna(df['gene'])

    return df

def calculate_cell_score(
    results,
    w_e=1.0, w_p=1.5, w_ct=2.5,             # base score weights
    w_m=7.0, w_f=3.0,                       # mean vs flatness weights
    alpha=1.0, beta=1.0, gamma=1.0,         # flatness breakdown
    w_c=1.0                                 # count weight
):

    cell_df = results.copy()

    # Calculate total cell count per cell
    sum_cell_count = cell_df.groupby('cell')['cell count'].sum().reset_index()
    sum_cell_count.columns = ['cell', 'total cell count']
    cell_df = cell_df.merge(sum_cell_count, on='cell')

    # Adaptive cutoff: max(μ - 0.5σ, min(total count))
    mu = cell_df['total cell count'].mean()
    sigma = cell_df['total cell count'].std()
    minimum = cell_df['total cell count'].min()
    min_count = max(mu - 0.5 * sigma, minimum)
    cell_df = cell_df[cell_df['total cell count'] >= min_count]

    # Normalize expression and cell count
    norm_expr = (cell_df['expression'] - cell_df['expression'].min()) / (cell_df['expression'].max() - cell_df['expression'].min())
    norm_count = (cell_df['cell count'] - cell_df['cell count'].min()) / (cell_df['cell count'].max() - cell_df['cell count'].min())

    # Compute base score
    cell_df['score'] = (
        w_e * norm_expr +
        w_p * cell_df['cell percentage'] +
        w_ct * norm_count
    )

    # Compute flatness components
    std_scores = cell_df.groupby('cell')['score'].std().reset_index().rename(columns={'score': 'std'})
    mean_scores = cell_df.groupby('cell')['score'].mean().reset_index().rename(columns={'score': 'mean'})
    kurtosis_scores = cell_df.groupby('cell')['score'].apply(pd.Series.kurt).reset_index(name='kurtosis')

    cell_df = cell_df.merge(std_scores, on='cell')
    cell_df = cell_df.merge(mean_scores, on='cell')
    cell_df = cell_df.merge(kurtosis_scores, on='cell')

    cell_df['CV'] = cell_df['std'] / (cell_df['mean'].replace(0, np.nan))

    # Normalize flatness metrics
    for col in ['std', 'CV', 'kurtosis', 'mean']:
        col_min, col_max = cell_df[col].min(), cell_df[col].max()
        cell_df[f'{col} norm'] = (cell_df[col] - col_min) / (col_max - col_min + 1e-9)

    # Stability (flatness) score
    stability = (
        alpha * (1 - cell_df['std norm'].fillna(0)) +
        beta * (1 - cell_df['CV norm'].fillna(0)) +
        gamma * (1 - cell_df['kurtosis norm'].fillna(0))
    )
    flatness_weight = stability / (alpha + beta + gamma)

    # Count-based log-scaled soft weight
    cell_counts = cell_df['cell'].value_counts().reset_index()
    cell_counts.columns = ['cell', 'count']
    cell_df = cell_df.merge(cell_counts, on='cell')

    log_scaled_count = np.log1p(cell_df['count']) / np.log1p(cell_df['count'].max())

    # Final cell score
    cell_df['cell score'] = (
        (w_m * cell_df['mean norm'] + w_f * flatness_weight) *
        (1 + w_c * log_scaled_count)
    )

    return cell_df.sort_values(by=['count', 'cell score', 'cell', 'gene'], ascending=[False, False, False, True]).reset_index(drop=True)

def run_data_processing(cluster, ensembl_ids, target_uberon_id):
    try:
        data = fetch_expression_data(cluster, ensembl_ids)

        response_df = expression_data_to_df(data)

        filtered_df = filter_expression_data(response_df, target_uberon_id)

        translated_df = translate_ontology(filtered_df)

        ranked_df = calculate_cell_score(translated_df)

        return ranked_df

    except Exception as e:
        print(f"An error occurred during the pipeline execution: {e}")
        return pd.DataFrame()
    
def main_data_analysis(cluster_dict, target_uberon_id):
    results = {}

    # Iterate over each cluster (column) in the DataFrame
    for cluster, name_list in cluster_dict.items():
        # Extract Ensembl IDs for the current cluster, dropping any NaN values
        ensembl_ids = name_list.dropna()
        
        if not ensembl_ids:
            continue
        
        # Run the expression pipeline for the current cluster
        result_df = run_data_processing(cluster, ensembl_ids, target_uberon_id)
        
        # Store the result in the dictionary
        results[cluster] = result_df.reset_index(drop=True)

    return results

def export_analysis_dataframe(results):
    ex_df = results.copy()
    for cluster in ex_df:
        unique_cells = ex_df[cluster]['cell'].unique()[:10]
        ex_df[cluster] = ex_df[cluster][['gene', 'cell', 'cell score']][ex_df[cluster]['cell'].isin(unique_cells)]
    return ex_df


#Validation
def initialize_tavily():
    # Initialize Tavily and Gemini clients
    return TavilyClient(api_key=os.getenv("TAVILY_API_KEY"))

def initialize_openai():
    return OpenAI(api_key=os.getenv("OPENAI_API_KEY"))

# Connect to PostgreSQL
def get_db_connection():
    return psycopg2.connect(os.getenv("DATABASE_URL"))

def check_database(tissue_name, cell_name, gene_name):
    query = """
        SELECT relationcertainty FROM gene_expression
        WHERE uberon = (
            SELECT uberon 
            FROM tissue 
            WHERE tissuename = %s
        )
        AND cl = (
            SELECT cl 
            FROM cell 
            WHERE cellname = %s
        )
        AND ensembl = (
            SELECT ensembl 
            FROM gene 
            WHERE genename = %s
        )
    """
    try:
        with get_db_connection() as conn:
            with conn.cursor() as cur:
                cur.execute(query, (tissue_name, cell_name, gene_name))
                result = cur.fetchone()
                return result[0] if result else None
    except Exception as e:
        print(f"Database error: {e}")
        return None
    
# Function to query Tavily using SDK
def query_tavily(tissue_name, cell_name, gene_name, tavily_client = initialize_tavily()):
    query = f"Does {cell_name} express {gene_name} in {tissue_name}?"
    try:
        response = tavily_client.search(query= query, search_depth= "basic")
        sources = [result["url"] for result in response.get("results", [])]
        return sources if sources else None
    except Exception as e:
        print(f"Tavily Error: {e}")
        return None
    

# Function to summarize using OpenAI
def summarize_with_openai(tissue_name, cell_name, gene_name, sources, client = initialize_openai()):
    prompt = f"""
    Given the following sources, determine the confidence level that {cell_name} expresses {gene_name} in {tissue_name}.
    Return one of these values in valid JSON format:
    {{"status": 1.0}} → Yes
    {{"status": 0.67}} → Mostly Yes
    {{"status": 0.33}} → Mostly No
    {{"status": 0.0}} → No
    never return with markdown formatting.
    Sources: {sources}
    """

    try:
        response = client.chat.completions.create(
            model="gpt-4o-mini",
            messages=[{"role": "system", "content": "You are an expert in cellular biology."},
                      {"role": "user", "content": prompt}]
        )
        output = response.choices[0].message.content
        return json.loads(output)["status"]
    except Exception as e:
        print(f"OpenAI Error: {e}")
        return None
    
def insert_into_database(tissue_name, cell_name, gene_name, source_urls, certainty):
    conn = psycopg2.connect(os.getenv("DATABASE_URL"))
    cur = conn.cursor()

    try:
        # Lookup IDs
        cur.execute("SELECT uberon FROM tissue WHERE tissuename = %s", (tissue_name,))
        uberon = cur.fetchone()
        if not uberon:
            raise ValueError(f"Tissue '{tissue_name}' not found.")
        uberon = uberon[0]

        cur.execute("SELECT cl FROM cell WHERE cellname = %s", (cell_name,))
        cl = cur.fetchone()
        if not cl:
            raise ValueError(f"Cell '{cell_name}' not found.")
        cl = cl[0]

        cur.execute("SELECT ensembl FROM gene WHERE genename = %s", (gene_name,))
        ensembl = cur.fetchone()
        if not ensembl:
            raise ValueError(f"Gene '{gene_name}' not found.")
        ensembl = ensembl[0]

        # Insert or update gene_expression
        cur.execute("""
            INSERT INTO gene_expression (uberon, cl, ensembl, relationcertainty)
            VALUES (%s, %s, %s, %s)
            ON CONFLICT (uberon, cl, ensembl)
            DO UPDATE SET relationcertainty = EXCLUDED.relationcertainty
            RETURNING expressionid
        """, (uberon, cl, ensembl, certainty))
        expressionid = cur.fetchone()[0]

        # Insert sources and validate links
        for url in source_urls:
            # Check or insert source
            cur.execute("SELECT sourceid FROM source WHERE url = %s", (url,))
            row = cur.fetchone()
            if row:
                sourceid = row[0]
            else:
                cur.execute("INSERT INTO source (url) VALUES (%s) RETURNING sourceid", (url,))
                sourceid = cur.fetchone()[0]

            # Insert validate link if not exists
            cur.execute("""
                INSERT INTO validate (expressionid, sourceid)
                VALUES (%s, %s)
                ON CONFLICT DO NOTHING
            """, (expressionid, sourceid))

        conn.commit()
        print(f"✅ Inserted/Updated expressionid {expressionid} with sources.")

    except Exception as e:
        conn.rollback()
        print(f"❌ Error: {e}")

    finally:
        cur.close()
        conn.close()


def validate_cell_gene_relation(tissue_name, cell_name, gene_name):
    print(f"\n[Validation Start] Processing: {tissue_name} - {cell_name} - {gene_name}")

    # Step 2: Check database first
    stored_status = check_database(tissue_name, cell_name, gene_name)
    if stored_status is not None:
        print(f"[Validation] Relation already exists in DB with status {stored_status}. Skipping process.")
        return float(stored_status)  # Ensure status is returned as a float

    # Step 3: Query Tavily
    try:
        sources = query_tavily(tissue_name, cell_name, gene_name)
    except Exception as e:
        print(f"Tavily Error: {e}")
        return None

    # Step 3a: If no sources found, set status to 0.0 and insert into database
    if not sources:
        status = 0.0
        print(f"[Validation] No references found in Tavily. Setting status to {status}.")
        insert_into_database(tissue_name, cell_name, gene_name, [], status)
        return float(status)  # Ensure status is a float

    # Step 4: Summarize with OpenAI
    status = summarize_with_openai(tissue_name, cell_name, gene_name, sources)
    status = float(status)  # Ensure status is a float
    print(f"[Validation] OpenAI summary status: {status}")

    # Step 5: Insert into Database
    insert_into_database(tissue_name, cell_name, gene_name, sources, status)
    print(f"[Validation] Inserted/Updated in DB with status: {status}")

    return status  # Return the validation score as a float

def main_data_validation(results, target_uberon_id):
    ex_df = results.copy()
    for cluster in ex_df:
        unique_cells = ex_df[cluster]['cell'].unique()[:10]
        ex_df[cluster] = ex_df[cluster][['gene', 'cell', 'cell score']][ex_df[cluster]['cell'].isin(unique_cells)]
    tissue_name = ontl.tissue(target_uberon_id)

    for cluster in ex_df:
        ex_df[cluster]["certainty score"] = ex_df[cluster].apply(lambda row: validate_cell_gene_relation(tissue_name, row["cell"], row["gene"]), axis=1)

    return ex_df

def export_validation_dataframe(results):
    ex_df = results.copy()
    for cluster in ex_df:
        unique_cells = ex_df[cluster]['cell'].unique()[:10]
        ex_df[cluster] = ex_df[cluster][['gene', 'cell', 'cell score', 'certainty score']][ex_df[cluster]['cell'].isin(unique_cells)]
    return ex_df