Merge pull request #85 from daisybio/dev

expression clustering, model filtering, new stringSearchTranscript

Author: strasserle

app/controllers/dataset.py

@@ -3,6 +3,7 @@
 import app.models as models
 from app.config import LATEST, db
 from typing import List
+import re
 
 
 def _dataset_query(query = None, sponge_db_version = LATEST, **kwargs):
@@ -174,3 +175,60 @@ def read_spongeRunInformation(dataset_ID: int = None, disease_name: str = None,
             "type": "about:blank",
             "data": []
         }), 200
+    
+def _extract_tss_code(sample_id):
+    """
+    Extracts the Tissue Source Site (TSS) code from a TCGA sample ID.
+
+    Args:
+        sample_id (str): The TCGA sample ID.
+
+    Returns:
+        str: The TSS code (2 characters) or None if not found.
+    """
+    match = re.match(r"TCGA-([A-Z0-9]{2})-[A-Z0-9]{4}-[A-Z0-9]{2}.*", sample_id)
+    return match.group(1) if match else None
+
+
+def get_disease_from_sample(sample_ID: str = None):
+    """
+    Get the disease name from a sample ID. This handles the api route /sponge/get_disease_from_sample?sample_ID={sample_ID}&sponge_db_version={sponge_db_version} 
+    Args:
+        sample_ID (str): The TCGA sample ID.
+        sponge_db_version (int): The version of the database.
+    Returns:
+        str: The disease name or None if not found.
+    """
+    # Extract the TSS code from the sample ID
+    
+    query = db.select(models.TissueSourceSite)
+    if sample_ID: 
+        tss_code = _extract_tss_code(sample_ID)
+
+        if tss_code is None:
+            return jsonify({
+                "detail": 'No valid sample ID.',
+                "status": 200,
+                "title": "No Content",
+                "type": "about:blank",
+                "data": []
+            }), 200
+        
+        query = query.where(models.TissueSourceSite.tissue_source_site_code == tss_code)
+
+    data = db.session.execute(query).scalars().all()
+
+    if len(data) == 0:
+        return jsonify({
+            "detail": 'Issue with sample ID: {sample_ID}'.format(sample_ID=sample_ID),
+            "status": 200,
+            "title": "No Content",
+            "type": "about:blank",
+            "data": []
+        }), 200
+
+    if sample_ID: 
+        return {sample_ID: data[0].disease_name}
+    else:
+        return {d.tissue_source_site_code: d.disease_name for d in data}   
+

app/controllers/diffExpression.py

@@ -4,17 +4,15 @@
 from app.controllers.dataset import _dataset_query
 from app.controllers.comparison import _comparison_query
 
-def get_diff_expr(dataset_ID_1: str = None, dataset_ID_2: int = None, disease_name_1=None, disease_name_2=None, disease_subtype_1=None, disease_subtype_2=None, condition_1=None, 
-                  condition_2=None, ensg_number=None, gene_symbol=None, sponge_db_version: int = LATEST):
+def get_diff_expr(dataset_ID_1: str = None, dataset_ID_2: int = None, 
+                  condition_1=None, condition_2=None, 
+                  ensg_number=None, gene_symbol=None, sponge_db_version: int = LATEST, 
+                  limit: int = None, offset: int = None):
     """
     API call /differentialExpression,
     get differential expression results between genes
     :param dataset_ID_1: dataset_ID of the first dataset of interest
     :param dataset_ID_2: dataset_ID of the second dataset of interest
-    :param disease_name_1: disease name of the first part of comparison (e.g. Sarcoma)
-    :param disease_name_2: disease name of the second part of comparison (e.g. Sarcoma)
-    :param disease_subtype_1: subtype of first part of comparison, overtype if none is provided (e.g. LMS) 
-    :param disease_subtype_2: subtype of second part of comparison, overtype if none is provided (e.g. LMS) 
     :param condition_1: condition of first part of comparison (e.g. disease, normal)
     :param condition_2: condition of second part of comparison (e.g. disease, normal)
     :param ensg_number: esng number of the gene(s) of interest
@@ -45,14 +43,17 @@ def get_diff_expr(dataset_ID_1: str = None, dataset_ID_2: int = None, disease_na
     if len(gene) > 0:
         gene_IDs = [i.gene_ID for i in gene]
 
-    dataset_1 = _dataset_query(sponge_db_version=sponge_db_version, disease_name=disease_name_1, dataset_ID=dataset_ID_1, disease_subtype=disease_subtype_1)
+    dataset_1 = _dataset_query(sponge_db_version=sponge_db_version, dataset_ID=dataset_ID_1)
     dataset_1 = [x.dataset_ID for x in dataset_1]
 
-    dataset_2 = _dataset_query(sponge_db_version=sponge_db_version, disease_name=disease_name_2, dataset_ID=dataset_ID_2, disease_subtype=disease_subtype_2)
+    dataset_2 = _dataset_query(sponge_db_version=sponge_db_version, dataset_ID=dataset_ID_2)
     dataset_2 = [x.dataset_ID for x in dataset_2]
 
     comparisons, reverse = _comparison_query(dataset_1, dataset_2, condition_1, condition_2, "gene")
 
+    if type(comparisons) is not list: 
+        return comparisons
+
     comparison_ID = comparisons[0].comparison_ID
 
     result = models.DifferentialExpression.query \
@@ -61,6 +62,10 @@ def get_diff_expr(dataset_ID_1: str = None, dataset_ID_2: int = None, disease_na
     if len(gene) > 0:
         result = result.filter(models.DifferentialExpression.gene_ID.in_(gene_IDs))
 
+    if limit is not None:
+        result = result.limit(limit)
+    if offset is not None:
+        result = result.offset(offset)
     result = result.all()
 
     if len(result) > 0:
@@ -81,17 +86,15 @@ def get_diff_expr(dataset_ID_1: str = None, dataset_ID_2: int = None, disease_na
         }), 200
 
 
-def get_diff_expr_transcript(dataset_ID_1: int = None, dataset_ID_2: int = None, disease_name_1=None, disease_name_2=None, disease_subtype_1=None, disease_subtype_2=None, 
-                             condition_1=None, condition_2=None, enst_number=None, sponge_db_version: int = LATEST):
+def get_diff_expr_transcript(dataset_ID_1: int = None, dataset_ID_2: int = None, 
+                             condition_1=None, condition_2=None, 
+                             enst_number=None, sponge_db_version: int = LATEST,
+                             limit: int = None, offset: int = None):
     """
     API call /differentialExpressionTranscript,
     get differential expression results between transcripts.
     :param dataset_ID_1: dataset_ID of the first dataset of interest
     :param dataset_ID_2: dataset_ID of the second dataset of interest
-    :param disease_name_1: disease name of the first part of comparison (e.g. Sarcoma)
-    :param disease_name_2: disease name of the second part of comparison (e.g. Sarcoma)
-    :param disease_subtype_1: subtype of first part of comparison, overtype if none is provided (e.g. LMS) 
-    :param disease_subtype_2: subtype of second part of comparison, overtype if none is provided (e.g. LMS) 
     :param condition_1: condition of first part of comparison (e.g. disease, normal)
     :param condition_2: condition of second part of comparison (e.g. disease, normal)
     :param enst_number: esng number of the transcript(s) of interest
@@ -109,10 +112,10 @@ def get_diff_expr_transcript(dataset_ID_1: int = None, dataset_ID_2: int = None,
     if len(transcript) > 0:
         transcript_IDs = [i.transcript_ID for i in transcript]
 
-    dataset_1 = _dataset_query(sponge_db_version=sponge_db_version, disease_name=disease_name_1, dataset_ID=dataset_ID_1, disease_subtype=disease_subtype_1)
+    dataset_1 = _dataset_query(sponge_db_version=sponge_db_version, dataset_ID=dataset_ID_1)
     dataset_1 = [x.dataset_ID for x in dataset_1]
 
-    dataset_2 = _dataset_query(sponge_db_version=sponge_db_version, disease_name=disease_name_2, dataset_ID=dataset_ID_2, disease_subtype=disease_subtype_1)
+    dataset_2 = _dataset_query(sponge_db_version=sponge_db_version, dataset_ID=dataset_ID_2)
     dataset_2 = [x.dataset_ID for x in dataset_2]
 
     comparisons, reverse = _comparison_query(dataset_1, dataset_2, condition_1, condition_2, "transcript")
@@ -124,6 +127,11 @@ def get_diff_expr_transcript(dataset_ID_1: int = None, dataset_ID_2: int = None,
     if len(transcript) > 0:
         result = result.filter(models.DifferentialExpression.transcript_ID.in_(transcript_IDs))
 
+    if limit is not None:
+        result = result.limit(limit)
+    if offset is not None:  
+        result = result.offset(offset)
+
     result = result.all()
 
     if len(result) > 0:

app/controllers/expressionValues.py

@@ -1,18 +1,21 @@
-from flask import jsonify
+from flask import Response, jsonify, stream_with_context
+from scipy.cluster.hierarchy import linkage, dendrogram
+import pandas as pd
 import app.models as models
 from app.config import LATEST, db
 
-def get_gene_expr(dataset_ID: int = None, disease_name=None, ensg_number=None, gene_symbol=None, sponge_db_version: int = LATEST):
+def get_gene_expr(dataset_ID: int = None, disease_name=None, ensg_number=None, gene_symbol=None, cluster: bool = False, limit: int = None, offset: int = None, sponge_db_version: int = LATEST):
     """
     Handles API call /exprValue/getceRNA to get gene expression values
     :param dataset_ID: dataset_ID of interest
     :param disease_name: disease_name of interest
     :param ensg_number: esng number of the gene of interest
     :param gene_symbol: gene symbol of the gene of interest
     :param sponge_db_version: version of the database
+    :param cluster: whether to cluster the gene expression (rows and columns)
+    :param limit: limit the number of results
     :return: all expression values for the genes of interest
     """
-
     # test if any of the two identification possibilities is given
     if ensg_number is None and gene_symbol is None:
         return jsonify({
@@ -82,9 +85,70 @@ def get_gene_expr(dataset_ID: int = None, disease_name=None, ensg_number=None, g
     result = models.GeneExpressionValues.query \
         .filter(*queries) \
         .all()
-
+    
     if len(result) > 0:
-        return models.geneExpressionSchema(many=True).dump(result)
+        # perform hierarchical clustering on rows and columns
+        if cluster:
+            # Convert result to a DataFrame for clustering
+            data = pd.DataFrame([{
+                "gene_ID": r.gene.gene_symbol if r.gene.gene_symbol else r.gene.ensg_number,
+                "sample_ID": r.sample_ID + "___" + (
+                    str(r.dataset.disease_name) if disease_name == "pancancer" else
+                    str(r.dataset.disease_subtype)
+                    ),
+                "expression_value": r.expr_value,
+            } for r in result])
+
+            # Pivot the data to create a matrix for clustering
+            expression_matrix = data.pivot(index="gene_ID", columns="sample_ID", values="expression_value").fillna(0)
+
+            # Perform hierarchical clustering on rows (genes) and columns (datasets)
+            try:
+                row_linkage = linkage(expression_matrix, method='ward', optimal_ordering=False)
+                col_linkage = linkage(expression_matrix.T, method='ward', optimal_ordering=False)
+            except ValueError as e:
+                # Handle the case where the data is not suitable for clustering
+                return jsonify({
+                    "detail": str(e),
+                    "status": 400,
+                    "title": "Bad Request",
+                    "type": "about:blank"
+                }), 400
+
+            # Add clustering results to the response
+            row_order = dendrogram(row_linkage, labels=expression_matrix.index, no_plot=True).get('leaves')
+            col_order = dendrogram(col_linkage, labels=expression_matrix.columns, no_plot=True).get('leaves')
+            expression_matrix = expression_matrix.iloc[row_order, col_order]
+
+            result = expression_matrix.reset_index().melt(id_vars='gene_ID', var_name='sample_ID', value_name='expression_value')
+            result = [models.GeneExpressionValues(gene={"gene_symbol": row['gene_ID'], "ensg_number": None}, 
+                                                    expr_value=row['expression_value'],
+                                                    sample_ID=row['sample_ID'], #.split('___')[0],
+                                                    # note that this is 'pancancer' if the disease is 'pancancer'
+                                                    dataset={"disease_subtype": row['sample_ID'].split('___')[1]},
+                                                    )
+                for _, row in result.iterrows()]
+            
+
+        # Limit the number of results if specified
+        if offset is not None:
+            result = result[offset:]
+        if limit is not None:
+            result = result[:limit]
+
+        def _generate():
+            yield "["
+            first = True
+            for r in result:
+                if not first:
+                    yield ","
+                yield models.geneExpressionSchema().dumps(r)
+                first = False
+            yield "]"
+            
+        return Response(stream_with_context(_generate()), content_type='application/json')
+        
+        # return models.geneExpressionSchema(many=True).dump(result)
     else:
         return jsonify({
             "detail": "No results.",
@@ -93,16 +157,21 @@ def get_gene_expr(dataset_ID: int = None, disease_name=None, ensg_number=None, g
             "type": "about:blank",
             "data": []
         }), 200
+        
+    
+
 
 
-def get_transcript_expression(dataset_ID: int = None, disease_name: str = None, enst_number: str = None, ensg_number: str = None, gene_symbol: str = None, sponge_db_version: int = LATEST):
+def get_transcript_expression(dataset_ID: int = None, disease_name: str = None, enst_number: str = None, ensg_number: str = None, gene_symbol: str = None, cluster: bool = False, limit: int = None, offset: int = None, sponge_db_version: int = LATEST):
     """
     Handles API call /exprValue/getTranscriptExpr to return transcript expressions
     :param dataset_ID: dataset_ID of interest
     :param disease_name: Name of the disease
     :param enst_number: Ensembl transcript ID
     :param ensg_number: Ensembl gene ID
     :param gene_symbol: gene symbol
+    :param limit: limit the number of results
+    :param cluster: whether to cluster the gene expression (rows and columns)
     :param sponge_db_version: version of the database
     :return: expression values for given search parameters
     """
@@ -187,7 +256,59 @@ def get_transcript_expression(dataset_ID: int = None, disease_name: str = None,
     result = db.session.execute(query).scalars().all()
 
     if len(result) > 0:
-        return models.ExpressionDataTranscriptSchema(many=True).dump(result)
+                # perform hierarchical clustering on rows and columns
+        if cluster:
+            # Convert result to a DataFrame for clustering
+            data = pd.DataFrame([{
+                "ensembl_ID": r.transcript.enst_number + "___" + (r.transcript.gene.gene_symbol if r.transcript.gene.gene_symbol else r.transcript.gene.ensg_number),
+                "sample_ID": r.sample_ID + "___" + (
+                    str(r.dataset.disease_name) if disease_name == "pancancer" else
+                    str(r.dataset.disease_subtype)
+                    ),
+                "expression_value": r.expr_value,
+                "expression_value": r.expr_value,
+            } for r in result])
+
+            # Pivot the data to create a matrix for clustering
+            expression_matrix = data.pivot(index="ensembl_ID", columns="sample_ID", values="expression_value").fillna(0)
+
+            # Perform hierarchical clustering on rows (genes) and columns (datasets)
+            row_linkage = linkage(expression_matrix, method='ward')
+            col_linkage = linkage(expression_matrix.T, method='ward')
+
+            # Add clustering results to the response
+            row_order = dendrogram(row_linkage, labels=expression_matrix.index, no_plot=True).get('leaves')
+            col_order = dendrogram(col_linkage, labels=expression_matrix.columns, no_plot=True).get('leaves')
+            expression_matrix = expression_matrix.iloc[row_order, col_order]
+
+            result = expression_matrix.reset_index().melt(id_vars='ensembl_ID', var_name='sample_ID', value_name='expression_value')
+            result = [models.ExpressionDataTranscript(
+                transcript={"enst_number": row['ensembl_ID'].split("___")[0], "gene": {"gene_symbol": row['ensembl_ID'].split('___')[1]}}, 
+                expr_value=row['expression_value'],
+                sample_ID=row['sample_ID'], #.split('___')[0],
+                # note that this is 'pancancer' if the disease is 'pancancer'
+                dataset={"disease_subtype": row['sample_ID'].split('___')[1]},
+            )for _, row in result.iterrows()]
+
+        # Limit the number of results if specified
+        if offset is not None:
+            result = result[offset:]
+        if limit is not None:
+            result = result[:limit]
+
+        def _generate():
+            yield "["
+            first = True
+            for r in result:
+                if not first:
+                    yield ","
+                yield models.ExpressionDataTranscriptSchema().dumps(r)
+                first = False
+            yield "]"
+            
+        return Response(stream_with_context(_generate()), content_type='application/json')
+    
+        # return models.ExpressionDataTranscriptSchema(many=True).dump(result)
     else:
         return jsonify({
             "detail": "No transcript expression data found for the given filters.",