alethiotx.artemis.clinical

397def compute(mesh_headings: list, chembl_version: str = '36', trials_only_last_n_years: int = None, filter_families: bool = True) -> dict:
398   """
399   Compute clinical validation scores for drug targets across multiple disease indications.
400   
401   This is the core function of the Artemis clinical scoring module. It processes ChEMBL drug-target-disease
402   relationships combined with clinical trial data to calculate comprehensive validation scores for each
403   target gene. These scores quantify the strength of clinical evidence supporting each target, enabling
404   data-driven target prioritization for drug discovery.
405   
406   **What This Function Does:**
407   
408   - Loads drug molecule, target, and clinical trial data from ChEMBL
409   - Expands diseases using MeSH hierarchies (e.g., "Lung Neoplasms" → all lung cancer subtypes)
410   - Counts clinical trials at each phase (0-4) for each drug-target pair
411   - Optionally filters out over-represented gene families to reduce bias
412   - Aggregates evidence across all drugs targeting the same gene
413   - Returns comprehensive target-level scores for each disease
414   
415   **Processing Workflow:**
416   
417   1. **Data Loading** (once for all diseases):
418      
419      - Downloads pre-computed ChEMBL molecules.csv from S3
420      - Contains drug-target-indication-trial relationships
421      - Optional temporal filtering to recent trials
422   
423   2. **Gene Family Analysis** (optional, once for all diseases):
424      
425      - Downloads HGNC gene family data
426      - Creates lookup table of drug-family representation percentages
427      - Identifies multi-target drugs (e.g., multi-kinase inhibitors)
428   
429   3. **Per-Disease Processing** (for each disease in mesh_headings):
430      
431   - **MeSH Expansion**: Retrieves all descendant disease terms
432      (e.g., "Lung Neoplasms" includes "Adenocarcinoma of Lung", "Small Cell Lung Carcinoma", etc.)
433   
434   - **Trial De-duplication**: Removes duplicate trials from:
435
436      - Same trial appearing in multiple mesh descendants
437      - Same trial with different molecular forms (salts) of the same parent drug
438   
439   - **Family Filtering** (optional): For each drug-disease combination:
440
441      - Identifies over-represented gene families (>40% of family targeted)
442      - Keeps one representative target per family (first alphabetically)
443      - Eliminates bias from multi-target drugs
444   
445   - **Trial Counting**: Counts trials per phase for each target:
446
447      - Phase 0-3: Experimental/developmental trials
448      - Phase 4: Approved drugs (post-marketing surveillance)
449   
450   - **Score Aggregation**: For each target gene:
451
452      - Sum all trial phases (0-3) across all drugs → phase_scores
453      - Check if any drug is approved (phase 4) → approved (boolean)
454      - Calculate final score: clinical_scores = phase_scores + (20 if approved else 0)
455   
456   - **Output**: DataFrame with one row per target, sorted by clinical_scores (descending)
457   
458   **Scoring System Details:**
459   
460   Each target receives a **clinical_scores** value based on:
461   
462   1. **Phase Scores** (0-3 points per trial):
463      
464   - Phase 0: 0 points (exploratory studies)
465   - Phase 1: 1 point (safety/dosing in small groups)
466   - Phase 2: 2 points (efficacy in larger groups)
467   - Phase 3: 3 points (large-scale confirmation)
468   - Phase 4: 0 points (counted separately as "approved")
469   
470   2. **Approval Bonus**: +20 points if target has any approved drug
471   
472   3. **Multi-Drug Aggregation**: If multiple drugs target the same gene, all their phase scores sum
473   
474   4. **Multi-Trial Aggregation**: Each distinct clinical trial contributes its phase score
475   
476   **Example Scoring:**
477   
478   - Target A: 2 Phase 2 trials + 1 Phase 3 trial = 2+2+3 = 7 points
479   - Target B: 1 Phase 3 trial + 1 approved drug = 3+20 = 23 points
480   - Target C: 3 drugs with Phase 2 trials + approved = (2+2+2)+20 = 26 points
481   
482   **Why Family Filtering Matters:**
483   
484   Multi-target drugs (e.g., sunitinib targeting 8 RTKs) can artificially inflate scores for
485   specific gene families. Without filtering, all 8 RTKs receive high scores from the same drug,
486   biasing results toward that family. With filtering (default), only 1 RTK representative is kept,
487   producing unbiased scores that better reflect biological diversity.
488   
489   :param mesh_headings: List of MeSH disease headings to analyze. Standard MeSH terms like
490   - 'Breast Neoplasms' (breast cancer)
491   - 'Lung Neoplasms' (lung cancer)
492   - 'Diabetes Mellitus, Type 2' (type 2 diabetes)
493   - 'Cardiovascular Diseases' (heart disease)
494   
495   Each heading is automatically expanded to include all MeSH descendant terms.
496   Use exact MeSH terminology for best results.
497   :type mesh_headings: list[str]
498
499   :param chembl_version: ChEMBL database version to use. Versions correspond to data release dates
500   - '36' (default): Data through 2024
501   - '35': Data through 2023
502   
503   See https://www.ebi.ac.uk/chembl/ for version details.
504   
505   :type chembl_version: str, optional
506
507   :param trials_only_last_n_years: Temporal filter for recent trials only. If provided
508   - Filters to trials registered in last N years
509   - Year inferred from ClinicalTrials.gov NCT IDs
510   - Useful for identifying emerging targets
511   
512   Examples:
513   - ``6``: Last 6 years (captures recent development)
514   - ``10``: Last decade
515   - ``None`` (default): All historical trials
516   :type trials_only_last_n_years: int or None, optional
517
518   :param filter_families: Whether to apply gene family filtering to reduce bias
519                           
520   - ``True`` (default): Filter over-represented families
521      (recommended for unbiased target prioritization)
522   - ``False``: Include all targets without filtering
523      (use when family representation is biologically relevant)
524   :type filter_families: bool, optional
525   :return: Dictionary mapping each MeSH heading to its results DataFrame.
526            
527   **Dictionary Structure:**
528   
529   - Keys: MeSH disease headings (str) - same as input ``mesh_headings``
530   - Values: DataFrames with target-level clinical scores
531   
532   **DataFrame Columns:**
533   
534   - ``target_gene_name`` (str): HGNC gene symbol (e.g., 'EGFR', 'TP53')
535   - ``phase_0`` (int): Count of phase 0 trials for this target
536   - ``phase_1`` (int): Count of phase 1 trials for this target
537   - ``phase_2`` (int): Count of phase 2 trials for this target
538   - ``phase_3`` (int): Count of phase 3 trials for this target
539   - ``phase_4`` (int): Count of phase 4 (approved) drugs for this target
540   - ``phase_scores`` (int): Sum of phase_1 + 2×phase_2 + 3×phase_3
541   - ``approved`` (bool): True if phase_4 > 0 (any approved drug exists)
542   - ``clinical_scores`` (int): phase_scores + (20 if approved else 0)
543   
544   DataFrames are sorted by ``clinical_scores`` descending (highest validation first).
545      
546   :rtype: dict[str, DataFrame]
547   
548   .. note::
549      **Efficiency Optimization**: ChEMBL data and gene family lookup tables are loaded once and
550      shared across all diseases. Processing 7 diseases takes ~5-10 minutes total, not per disease.
551   
552   .. note::
553      **MeSH Hierarchy Expansion**: Diseases are automatically expanded to include all subtypes.
554      For example, "Breast Neoplasms" includes "Ductal Carcinoma", "Lobular Carcinoma", etc.
555      This ensures comprehensive coverage of clinical evidence.
556   
557   .. note::
558      **Trial De-duplication**: The function carefully de-duplicates clinical trials to avoid
559      double-counting. Trials are unique by (drug, NCT ID, target, disease) to ensure accurate counts.
560   
561   .. warning::
562      **AWS S3 Access Required**: This function reads ChEMBL and MeSH data from the
563      ``alethiotx-artemis`` S3 bucket. Configure AWS credentials via environment variables
564      or ``aws configure`` before calling.
565   
566   .. warning::
567      **Memory Usage**: ChEMBL v36 molecules.csv is ~500MB-1GB in memory. Processing multiple
568      diseases simultaneously requires sufficient RAM. For memory-constrained environments,
569      process diseases in smaller batches.
570   
571   .. warning::
572      **Processing Time**: Full computation takes 5-10 minutes for 7 diseases with family filtering.
573      Set ``filter_families=False`` to reduce runtime by ~50% (but scores may be biased).
574   
575   .. warning::
576      **Family Filtering Off**: When ``filter_families=False``, results may show artificial
577      enrichment for gene families with multi-target drugs. Use this option only when
578      family-level representation is biologically meaningful for your analysis.
579   
580   **Examples**
581   
582   Basic usage - single disease::
583   
584   >>> from alethiotx.artemis.clinical import compute
585   >>> results = compute(['Lung Neoplasms'])
586   >>> lung_targets = results['Lung Neoplasms']
587   >>> print(f"Found {len(lung_targets)} lung cancer targets")
588   >>> print(lung_targets.head())
589      target_gene_name  phase_0  phase_1  phase_2  phase_3  phase_4  phase_scores  approved  clinical_scores
590   0  EGFR             0        5        15       8        3         45           True      65
591   1  ALK              0        2        8        6        2         32           True      52
592      
593   Multiple diseases simultaneously::
594   
595   >>> diseases = ['Breast Neoplasms', 'Lung Neoplasms', 'Prostatic Neoplasms']
596   >>> results = compute(diseases, chembl_version='36', filter_families=True)
597   >>> for disease, df in results.items():
598   ...     approved_count = df['approved'].sum()
599   ...     print(f"{disease}: {len(df)} total targets, {approved_count} approved")
600   
601   Recent trials only (emerging targets)::
602   
603   >>> # Focus on targets with recent clinical activity (last 6 years)
604   >>> results = compute(
605   ...     ['Lung Neoplasms'],
606   ...     trials_only_last_n_years=6,
607   ...     filter_families=True
608   ... )
609   >>> recent_lung = results['Lung Neoplasms']
610   >>> print("Emerging lung cancer targets:", recent_lung['target_gene_name'].tolist()[:10])
611   
612   Disable family filtering for comprehensive view::
613   
614   >>> # Include all targets without filtering (may show family bias)
615   >>> results = compute(['Lung Neoplasms'], filter_families=False)
616   >>> lung_all = results['Lung Neoplasms']
617   >>> print(f"All targets (no filtering): {len(lung_all)}")
618   
619   Filter to approved targets only::
620   
621   >>> results = compute(['Diabetes Mellitus, Type 2'])
622   >>> diabetes = results['Diabetes Mellitus, Type 2']
623   >>> approved_diabetes = diabetes[diabetes['approved'] == True]
624   >>> print(f"Approved diabetes targets: {len(approved_diabetes)}")
625   >>> print(approved_diabetes[['target_gene_name', 'clinical_scores']].head())
626   
627   Compare targets across diseases::
628   
629   >>> diseases = ['Breast Neoplasms', 'Lung Neoplasms', 'Prostatic Neoplasms']
630   >>> results = compute(diseases)
631   >>> target = 'EGFR'
632   >>> for disease in diseases:
633   ...     df = results[disease]
634   ...     if target in df['target_gene_name'].values:
635   ...         score = df[df['target_gene_name'] == target]['clinical_scores'].iloc[0]
636   ...         print(f"{disease}: {target} = {score} points")
637   
638   Save results to CSV files::
639   
640   >>> results = compute(['Breast Neoplasms', 'Lung Neoplasms'])
641   >>> for disease, df in results.items():
642   ...     filename = disease.replace(' ', '_').replace(',', '') + '.csv'
643   ...     df.to_csv(filename, index=False)
644   ...     print(f"Saved {filename}")
645   
646   Batch processing for all major diseases::
647   
648   >>> all_diseases = [
649   ...     'Breast Neoplasms',
650   ...     'Lung Neoplasms', 
651   ...     'Prostatic Neoplasms',
652   ...     'Skin Neoplasms',
653   ...     'Intestinal Neoplasms',
654   ...     'Diabetes Mellitus, Type 2',
655   ...     'Cardiovascular Diseases'
656   ... ]
657   >>> results = compute(all_diseases, trials_only_last_n_years=6)
658   >>> print(f"Processed {len(results)} diseases in one batch")
659   
660   .. seealso::
661      - :func:`lookup_drug_family_representation`: Creates family lookup table used when ``filter_families=True``
662      - :func:`filter_overrepresented_families`: Filters individual drug-disease combinations
663      - :func:`load`: Load pre-computed scores from S3 (faster than computing from scratch)
664      - :func:`~alethiotx.artemis.mesh.descendants`: Expands MeSH terms to include all subtypes
665      - :func:`~alethiotx.artemis.chembl.molecules`: Retrieves raw ChEMBL molecule data
666   """
667   print("="*80)
668   print("Loading ChEMBL data...")
669   print("="*80)
670   
671   # Load ChEMBL clinical trials data once for all diseases (efficiency optimization)
672   chembl = read_csv(f's3://alethiotx-artemis/data/chembl/{chembl_version}/molecules.csv')
673   
674   # Apply temporal filter if requested (e.g., only trials from last 6 years)
675   if trials_only_last_n_years is not None:
676      current_year = int(date.today().strftime("%Y"))
677      cutoff_year = current_year - trials_only_last_n_years
678      chembl = chembl[chembl['trial_year'] >= cutoff_year]
679      print(f"Filtered to trials from last {trials_only_last_n_years} years (>= {cutoff_year})")
680   
681   # Create gene family representation lookup table once (shared across all diseases)
682   # This identifies which drugs over-represent specific gene families
683   lookup_table = None
684   if filter_families:
685      print("\nCreating drug-family representation lookup table...")
686      lookup_table = lookup_drug_family_representation(chembl)
687      print(f"Lookup table created with {len(lookup_table)} drug-mesh-family combinations")
688   else:
689      print("\nFamily filtering disabled - all targets will be included")
690   
691   # Process each disease indication independently
692   results = {}
693   
694   for idx, mesh_heading in enumerate(mesh_headings, 1):
695      print("\n" + "="*80)
696      print(f"Processing disease {idx}/{len(mesh_headings)}: {mesh_heading}")
697      print("="*80)
698      
699      # Step 1: Get all MeSH descendants first
700      # Example: "Lung Neoplasms" includes "Adenocarcinoma of Lung", "Small Cell Lung Carcinoma", etc.
701      mesh_descendants_list = mesh_descendants(
702         mesh_heading,
703         s3_base='s3://alethiotx-artemis/data/mesh/',
704         file_base='d2025'
705      )
706      print(f"Including {len(mesh_descendants_list)} MeSH headings (including descendants)")
707      
708      # Filter ChEMBL data to this disease and its descendants
709      disease_data = chembl[chembl['mesh_heading'].isin(mesh_descendants_list)].copy()
710      
711      # Step 2: De-duplicate clinical trials
712      # Each row in ChEMBL represents one clinical trial. We need to avoid counting:
713      # - Same trial with different child molecules (salt forms) of the same parent
714      # - Same trial appearing in multiple mesh descendants
715      # De-duplicate by (parent_molregno, clinical_trial_id, mesh_heading) first
716      if 'clinical_trial_id' in disease_data.columns:
717         disease_data = disease_data.drop_duplicates(
718            subset=['parent_molregno', 'clinical_trial_id', 'mesh_heading', 'target_gene_name']
719         )
720      
721      # Step 3: Group by mesh_heading + parent_molregno + target to aggregate trials per mesh
722      # This intermediate step is needed for family filtering (which operates per mesh heading)
723      parent_per_mesh = disease_data.groupby(
724         ['mesh_heading', 'parent_molregno', 'target_gene_name'], 
725         as_index=False
726      ).agg(
727         max_phase=('phase', 'max'),  # Keep max phase for this parent in this mesh
728         chembl_id=('chembl_id', 'first')
729      )
730      
731      # Step 4: Identify approved (per mesh)
732      parent_per_mesh['approved'] = (parent_per_mesh['max_phase'] == 4)
733      
734      # Step 5: Calculate phase scores (per mesh)
735      parent_per_mesh['phase_scores'] = parent_per_mesh['max_phase'].apply(
736         lambda p: p if p != 4 else 0
737      )
738      
739      res = parent_per_mesh
740      
741      # Handle empty result case (no trials for this disease)
742      if len(res) == 0:
743         print(f"Warning: No data found for {mesh_heading}")
744         results[mesh_heading] = DataFrame(
745            columns=['target_gene_name', 'phase_scores', 'approved', 'clinical_scores']
746         )
747         continue
748
749      # Step 5: Filter over-represented gene families per drug across ALL mesh descendants (optional)
750      # This prevents multi-target drugs (e.g., multi-kinase inhibitors) from biasing results
751      # IMPORTANT: Filter at the parent mesh_heading level (not per descendant) to avoid
752      # keeping multiple representatives from different sub-indications
753      if filter_families:
754         unique_drugs = res['chembl_id'].drop_duplicates()
755         print(f"Filtering over-represented families for {len(unique_drugs)} drugs across all descendants...")
756         
757         filtered_results = []
758         for drug_id in unique_drugs:
759            # Extract all targets for this drug across ALL mesh descendants
760            drug_data = res[res['chembl_id'] == drug_id].copy()
761            
762            # Get unique mesh headings for this drug
763            unique_mesh = drug_data['mesh_heading'].unique()
764            
765            # Apply family filtering FOR EACH mesh heading separately
766            # This is critical because we need to filter each mesh heading's records independently
767            drug_filtered = []
768            for mesh in unique_mesh:
769               mesh_data = drug_data[drug_data['mesh_heading'] == mesh].copy()
770               
771               # Try filtering with the specific mesh first, fall back to parent if no data
772               filtered_mesh = filter_overrepresented_families(
773                  mesh_data, 
774                  drug_chembl_id=drug_id,
775                  mesh_heading=mesh,  # Use specific mesh heading
776                  lookup_table=lookup_table
777               )
778               
779               # If no filtering happened (no family data for specific mesh), try parent heading
780               if len(filtered_mesh) == len(mesh_data) and mesh != mesh_heading:
781                  filtered_mesh = filter_overrepresented_families(
782                     mesh_data, 
783                     drug_chembl_id=drug_id,
784                     mesh_heading=mesh_heading,  # Fall back to parent heading
785                     lookup_table=lookup_table
786                  )
787               
788               drug_filtered.append(filtered_mesh)
789            
790            # Combine all mesh headings for this drug
791            if len(drug_filtered) > 0:
792               filtered_results.append(pd_concat(drug_filtered, ignore_index=True))
793         
794         # Combine all filtered drug data
795         res = pd_concat(filtered_results, ignore_index=True)
796         print(f"After filtering: {len(res)} records, {res['target_gene_name'].nunique()} unique targets")
797      else:
798         print("Skipping family filtering (filter_families=False)")
799      
800      # Step 6: Use FILTERED data and go back to trial level for counting
801      # Important: Use 'res' (after family filtering) to preserve filtering results
802      # Get the filtered parent-target combinations
803      filtered_combinations = res[['parent_molregno', 'target_gene_name', 'chembl_id']].drop_duplicates()
804      
805      # Now go back to original trial data but only keep filtered combinations
806      # This preserves family filtering while counting ALL trials
807      if 'clinical_trial_id' in disease_data.columns:
808         # Merge to keep only filtered parent-target combinations
809         trials_filtered = disease_data.merge(
810            filtered_combinations,
811            on=['parent_molregno', 'target_gene_name'],
812            how='inner'
813         )
814         
815         # De-duplicate by clinical trial ID globally across mesh descendants
816         trials_global = trials_filtered.drop_duplicates(
817            subset=['parent_molregno', 'clinical_trial_id', 'target_gene_name']
818         ).copy()  # Create explicit copy to avoid SettingWithCopyWarning
819      else:
820         # Fallback: if no trial IDs, use the aggregated filtered data
821         trials_global = res.drop_duplicates(
822            subset=['parent_molregno', 'target_gene_name', 'max_phase']
823         ).copy()
824      
825      # Now aggregate all trials per target
826      trials_global['approved'] = (trials_global['phase'] == 4)
827      trials_global['phase_scores'] = trials_global['phase'].apply(
828         lambda p: p if p != 4 else 0
829      )
830      
831      print(f"After de-duplicating trials globally: {len(trials_global)} trial records, {trials_global['target_gene_name'].nunique()} unique targets")
832      
833      # Step 7: Count TRIALS per phase for each target
834      # Each clinical trial in a given phase counts towards that phase's total
835      phase_counts = trials_global.groupby('target_gene_name', as_index=False).agg(
836         phase_0=('phase', lambda x: (x == 0).sum()),
837         phase_1=('phase', lambda x: (x == 1).sum()),
838         phase_2=('phase', lambda x: (x == 2).sum()),
839         phase_3=('phase', lambda x: (x == 3).sum()),
840         phase_4=('phase', lambda x: (x == 4).sum())
841      )
842      
843      # Step 8: Aggregate ALL trial scores per target
844      # Sum all clinical trial phases (0-3 for non-approved, 0 for approved)
845      # Each trial contributes its phase score to the cumulative total
846      phase_scores_final = trials_global.groupby('target_gene_name', as_index=False).agg(
847         phase_scores=('phase_scores', 'sum')  # Sum ALL trial phase scores
848      )
849      
850      # Check if target has any approved drug (any trial in phase 4)
851      approved_final = trials_global.groupby('target_gene_name', as_index=False).agg(
852         approved=('approved', lambda x: x.any())  # True if ANY trial is phase 4
853      )
854      
855      # Merge all components
856      res_final = merge(phase_scores_final, approved_final, on='target_gene_name')
857      res_final = merge(res_final, phase_counts, on='target_gene_name')
858      
859      # Step 9: Calculate final clinical scores with approval bonus
860      # Formula: clinical_score = phase_scores + (20 if approved else 0)
861      res_final['clinical_scores'] = res_final['phase_scores'] + 20 * res_final['approved'].astype(int)
862      
863      # Reorder columns: target_gene_name, phase_0-4, phase_scores, approved, clinical_scores
864      column_order = ['target_gene_name', 'phase_0', 'phase_1', 'phase_2', 'phase_3', 'phase_4',
865                      'phase_scores', 'approved', 'clinical_scores']
866      res_final = res_final[column_order]
867      
868      # Sort by clinical score descending (highest validation first)
869      res_final.sort_values('clinical_scores', ascending=False, inplace=True)
870      res_final.reset_index(drop=True, inplace=True)
871      
872      results[mesh_heading] = res_final
873      print(f"Completed {mesh_heading}: {len(res_final)} targets")
874   
875   return results

877def load(date = '2025-12-08'):
878   """
879   Load pre-computed clinical validation scores for seven major disease types from S3 storage.
880
881   This convenience function retrieves pre-calculated clinical target scores that have been
882   computed using :func:`compute` and stored in the Alethio Therapeutics public S3 bucket.
883   Each disease has its own CSV file containing target-level clinical validation scores.
884   
885   **Disease Coverage:**
886   
887   Returns scores for these seven diseases (in order):
888   
889   1. Breast cancer (Breast Neoplasms)
890   2. Lung cancer (Lung Neoplasms)
891   3. Prostate cancer (Prostatic Neoplasms)
892   4. Melanoma (Skin Neoplasms)
893   5. Bowel/colorectal cancer (Intestinal Neoplasms)
894   6. Type 2 diabetes (Diabetes Mellitus, Type 2)
895   7. Cardiovascular disease (Cardiovascular Diseases)
896
897   :param date: Date stamp identifying which pre-computed dataset to load. Format: ``YYYY-MM-DD``.
898                Different dates represent different ChEMBL versions, temporal filters, or
899                computation parameters. Defaults to ``2025-12-08``.
900   :type date: str, optional
901   :return: Tuple of 7 DataFrames containing clinical scores for each disease. Each DataFrame has:
902            
903   **Columns:**
904   
905   - ``Target Gene`` (str): HGNC gene symbol
906   - ``phase_0`` through ``phase_4`` (int): Count of trials in each phase
907   - ``Phase Score`` (int): Sum of trial phases (0-3)
908   - ``approved`` (bool): Whether target has any approved (phase 4) drug
909   - ``Clinical Score`` (int): Total validation score (Phase Score + 20 if approved)
910   
911   **Order:** 
912   
913   (breast, lung, prostate, melanoma, bowel, diabetes, cardiovascular)
914
915   :rtype: tuple[DataFrame, DataFrame, DataFrame, DataFrame, DataFrame, DataFrame, DataFrame]
916   :raises FileNotFoundError: If CSV files don't exist at the specified S3 paths for the given date
917   :raises ValueError: If CSV files are malformed or missing required columns
918   :raises botocore.exceptions.NoCredentialsError: If AWS credentials are not configured
919
920   .. note::
921      **S3 Storage Structure**: Files are located at:
922      ``s3://alethiotx-artemis/data/clinical_scores/{date}/{disease}.csv``
923      
924      Where ``{disease}`` is: breast, lung, prostate, melanoma, bowel, diabetes, or cardiovascular
925
926   .. note::
927      **AWS Configuration**: Requires AWS credentials with read access to the ``alethiotx-artemis``
928      S3 bucket. Configure via environment variables (``AWS_ACCESS_KEY_ID``, ``AWS_SECRET_ACCESS_KEY``)
929      or AWS CLI (``aws configure``).
930
931   .. warning::
932      Verify the date exists in S3 before calling. Available dates depend on when scores were
933      computed and uploaded. Contact Alethio Therapeutics for a list of available dates.
934
935   .. warning::
936      Each DataFrame can be large (thousands of rows). Consider filtering to specific targets
937      or approved-only drugs after loading to reduce memory usage.
938
939   **Examples**
940   
941   Load latest pre-computed scores::
942   
943      >>> from alethiotx.artemis.clinical import load
944      >>> breast, lung, prostate, melanoma, bowel, diabetes, cardio = load(date='2025-12-08')
945      >>> print(f"Breast cancer: {len(breast)} targets")
946      >>> print(breast.head())
947   
948   Load and filter to approved targets only::
949   
950      >>> breast, lung, prostate, melanoma, bowel, diabetes, cardio = load()
951      >>> breast_approved = breast[breast['approved'] == True]
952      >>> print(f"Approved breast cancer targets: {len(breast_approved)}")
953   
954   Compare scores across diseases for a specific gene::
955   
956      >>> breast, lung, prostate, melanoma, bowel, diabetes, cardio = load()
957      >>> diseases = [breast, lung, prostate, melanoma, bowel, diabetes, cardio]
958      >>> names = ['Breast', 'Lung', 'Prostate', 'Melanoma', 'Bowel', 'Diabetes', 'Cardio']
959      >>> target = 'EGFR'
960      >>> for name, df in zip(names, diseases):
961      ...     if target in df['Target Gene'].values:
962      ...         score = df[df['Target Gene'] == target]['Clinical Score'].iloc[0]
963      ...         print(f"{name}: {score}")
964   
965   Load older dataset version::
966   
967      >>> # Load scores computed on an earlier date
968      >>> breast, lung, prostate, melanoma, bowel, diabetes, cardio = load(date='2025-11-11')
969   
970   .. seealso::
971      - :func:`compute`: Function used to generate these pre-computed scores
972      - :func:`approved`: Filter loaded scores to approved targets only
973   """
974   breast = read_csv("s3://alethiotx-artemis/data/clinical_scores/" + date + "/breast.csv")
975   lung = read_csv("s3://alethiotx-artemis/data/clinical_scores/" + date + "/lung.csv")
976   prostate = read_csv("s3://alethiotx-artemis/data/clinical_scores/" + date + "/prostate.csv")
977   melanoma = read_csv("s3://alethiotx-artemis/data/clinical_scores/" + date + "/melanoma.csv")
978   bowel = read_csv("s3://alethiotx-artemis/data/clinical_scores/" + date + "/bowel.csv")
979   diabetes = read_csv("s3://alethiotx-artemis/data/clinical_scores/" + date + "/diabetes.csv")
980   cardiovascular = read_csv("s3://alethiotx-artemis/data/clinical_scores/" + date + "/cardiovascular.csv")
981
982   return(breast, lung, prostate, melanoma, bowel, diabetes, cardiovascular)

  8def lookup_drug_family_representation(chembl: DataFrame) -> DataFrame:
  9   """
 10   Create a comprehensive lookup table mapping drug-disease-family combinations to gene family representation.
 11   
 12   This function analyzes drug-target relationships from ChEMBL clinical trial data to identify potential
 13   bias from drugs that target many genes within the same family. For example, multi-kinase inhibitors
 14   targeting 10 out of 15 kinases in a family would represent 67% of that family, potentially inflating
 15   the importance of kinase targets in clinical scoring.
 16   
 17   The lookup table enables downstream filtering of over-represented families to produce unbiased
 18   target prioritization scores.
 19   
 20   **Processing Workflow:**
 21   
 22   1. **Download HGNC Data**: Retrieves complete gene family definitions from HGNC
 23   2. **Extract Relationships**: Identifies all unique drug-target-disease combinations from ChEMBL
 24   3. **Family Mapping**: Maps each drug target to its HGNC gene family (if assigned)
 25   4. **Calculate Representation**: For each drug-disease-family combination, computes:
 26      
 27      - Number of family genes targeted by the drug
 28      - Total genes in the family
 29      - Percentage representation (targets/total * 100)
 30   
 31   5. **Return Lookup**: Creates sorted lookup table for efficient querying
 32   
 33   **Why This Matters:**
 34   
 35   When a drug targets multiple members of the same gene family (e.g., multi-kinase inhibitors,
 36   pan-HDAC inhibitors), each family member contributes to clinical scores independently. This
 37   artificially inflates the importance of that family in target prioritization. By quantifying
 38   family representation, we can identify and filter these cases to prevent bias.
 39   
 40   :param chembl: ChEMBL molecule data containing drug-target-indication relationships.
 41                  Must include columns: ``chembl_id``, ``target_gene_name``, ``mesh_heading``
 42   :type chembl: DataFrame
 43   :return: Lookup table with one row per drug-disease-family combination. 
 44   
 45   **Columns:**
 46      
 47      - ``drug_chembl_id`` (str): ChEMBL identifier for the drug (e.g., 'CHEMBL123')
 48      - ``mesh_heading`` (str): MeSH disease term (e.g., 'Lung Neoplasms')
 49      - ``family_id`` (int): HGNC gene family identifier
 50      - ``family_abbreviation`` (str): Short family name (e.g., 'RTK' for receptor tyrosine kinases)
 51      - ``family_name`` (str): Full descriptive family name (e.g., 'Receptor tyrosine kinases')
 52      - ``total_genes_in_family`` (int): Total number of genes in this HGNC family
 53      - ``targets_in_family`` (int): Number of family genes targeted by this drug for this disease
 54      - ``target_genes_in_family`` (str): Comma-separated list of targeted gene symbols
 55      - ``pct_represented`` (float): Percentage of family genes targeted (0.0-100.0)
 56      
 57      Sorted by: drug_chembl_id, mesh_heading, then pct_represented (descending)
 58      
 59   :rtype: DataFrame
 60   
 61   .. note::
 62      Only includes targets with HGNC family assignments. Targets without family data are
 63      excluded from the lookup table (they pass through unfiltered in downstream processing).
 64   
 65   .. note::
 66      This function downloads fresh HGNC data on each call (~10-30 seconds). For production
 67      workflows processing multiple diseases, call this function once and reuse the lookup table.
 68   
 69   .. warning::
 70      Requires internet connection to download HGNC gene family data from public repositories.
 71      Ensure network access to HGNC/Google Cloud Storage endpoints.
 72   
 73   **Examples**
 74   
 75   Basic usage::
 76   
 77      >>> from alethiotx.artemis.chembl import molecules
 78      >>> chembl_data = molecules(version='36')
 79      >>> lookup = lookup_drug_family_representation(chembl_data)
 80      >>> print(f"Created lookup with {len(lookup)} drug-disease-family combinations")
 81   
 82   Find drugs over-representing kinase families::
 83   
 84      >>> kinase_overrep = lookup[
 85      ...     (lookup['family_name'].str.contains('kinase', case=False)) &
 86      ...     (lookup['pct_represented'] > 50)
 87      ... ]
 88      >>> print(kinase_overrep[['drug_chembl_id', 'family_name', 'pct_represented']].head())
 89   
 90   Identify drugs with many targets in small families::
 91   
 92      >>> small_families = lookup[
 93      ...     (lookup['total_genes_in_family'] < 10) &
 94      ...     (lookup['targets_in_family'] >= 4)
 95      ... ]
 96      >>> for _, row in small_families.iterrows():
 97      ...     print(f"{row['drug_chembl_id']}: {row['targets_in_family']}/{row['total_genes_in_family']} "
 98      ...           f"of {row['family_name']}")
 99   
100   .. seealso::
101      - :func:`filter_overrepresented_families`: Uses this lookup table to filter biased targets
102      - :func:`compute`: Integrates this lookup for unbiased clinical score calculation
103   """
104   # Download HGNC gene family data
105   gene_has_family, family, hgnc_data = hgnc_download()
106   family_summary, gene_family_mapping = hgnc_process(gene_has_family, family, hgnc_data)
107   
108   # Extract unique drug-target-disease combinations from ChEMBL
109   # Use vectorized operations for better performance
110   valid_rows = chembl[['chembl_id', 'target_gene_name', 'mesh_heading']].copy()
111   valid_rows = valid_rows[
112      notna(valid_rows['target_gene_name']) & notna(valid_rows['mesh_heading'])
113   ]
114   valid_rows.rename(columns={'chembl_id': 'drug_chembl_id'}, inplace=True)
115   targets_df = valid_rows.drop_duplicates()
116   
117   print(f"\nTotal drug-target-disease combinations: {len(targets_df)}")
118   print(f"Unique targets: {targets_df['target_gene_name'].nunique()}")
119   print(f"Unique drugs: {targets_df['drug_chembl_id'].nunique()}")
120   
121   # Map targets to their HGNC gene families
122   targets_with_families = targets_df.merge(
123      gene_family_mapping[['symbol', 'family_id', 'abbreviation', 'name_family']],
124      left_on='target_gene_name',
125      right_on='symbol',
126      how='left'
127   )
128   
129   # Report mapping success rate
130   with_family = targets_with_families['family_id'].notna().sum()
131   without_family = targets_with_families['family_id'].isna().sum()
132   print(f"  Targets with family: {with_family}, without family: {without_family}")
133   
134   # Keep only targets with family assignments for representation calculation
135   targets_with_families = targets_with_families[
136      targets_with_families['family_id'].notna()
137   ].copy()
138   
139   # Calculate family representation metrics
140   # Group by drug-disease-family to count targets per family
141   family_counts = targets_with_families.groupby(
142      ['drug_chembl_id', 'mesh_heading', 'family_id'],
143      as_index=False
144   ).agg(
145      targets_in_family=('target_gene_name', 'nunique'),
146      target_genes_in_family=('symbol', lambda x: ', '.join(sorted(x.unique())))
147   )
148   
149   # Add family metadata (name, abbreviation, total gene count)
150   family_counts = family_counts.merge(
151      family_summary[['id', 'abbreviation', 'name', 'total_genes_in_family']],
152      left_on='family_id',
153      right_on='id',
154      how='left'
155   )
156   
157   # Calculate percentage representation: (targets in family / total genes in family) * 100
158   family_counts['pct_represented'] = (
159      100 * family_counts['targets_in_family'] / family_counts['total_genes_in_family']
160   ).round(1)
161   
162   # Select and organize final columns
163   result = family_counts[[
164      'drug_chembl_id',
165      'mesh_heading',
166      'family_id',
167      'abbreviation',
168      'name',
169      'total_genes_in_family',
170      'targets_in_family',
171      'target_genes_in_family',
172      'pct_represented'
173   ]].copy()
174   
175   # Rename for clarity
176   result.rename(columns={
177      'abbreviation': 'family_abbreviation',
178      'name': 'family_name'
179   }, inplace=True)
180   
181   # Sort by drug, disease, then highest representation first
182   result.sort_values(
183      ['drug_chembl_id', 'mesh_heading', 'pct_represented'],
184      ascending=[True, True, False],
185      inplace=True
186   )
187   
188   return result

190def filter_overrepresented_families(
191   targets_df: DataFrame,
192   drug_chembl_id: str,
193   mesh_heading: str,
194   lookup_table: DataFrame,
195   pct_threshold: float = 40.0,
196   small_family_threshold: int = 5,
197   min_targets_threshold: int = 3
198) -> DataFrame:
199   """
200   Filter over-represented gene families from drug targets to eliminate family-based bias in clinical scoring.
201   
202   Multi-target drugs (e.g., multi-kinase inhibitors like sunitinib) often target multiple members of the
203   same gene family. When computing clinical validation scores, each family member contributes independently,
204   artificially inflating the perceived importance of that family. This function identifies over-represented
205   families and keeps only one representative target per family to produce unbiased scores.
206   
207   **Why Filtering is Needed:**
208   
209   Consider sunitinib, which targets 8 receptor tyrosine kinases (RTKs). Without filtering, all 8 RTKs
210   contribute to clinical scores, making RTKs appear highly validated. With filtering, we keep only 1
211   representative RTK (e.g., FLT3 alphabetically), preventing the family from dominating results.
212   
213   **Filtering Logic:**
214   
215   A gene family is considered "over-represented" and filtered if **ALL** of these conditions are met:
216   
217   1. **Sufficient targets**: Drug targets > ``min_targets_threshold`` genes in the family (default: >3)
218   2. **Size-adjusted criteria** (either condition triggers filtering):
219      
220      a. **Large family with high %**: Family size ≥ ``small_family_threshold`` (default: ≥5)
221         AND representation > ``pct_threshold``% (default: >40%)
222      b. **Absolute threshold**: Drug targets ≥ 4 genes in family (regardless of family size)
223   
224   **Exemptions from Filtering:**
225   
226   - **Small families** (<5 genes by default): Never filtered to preserve biological specificity
227     (e.g., the 4-gene HDAC class IIA family)
228   - **Low target count** (≤3 targets by default): Insufficient redundancy to justify filtering
229   - **Targets without family assignments**: Pass through unfiltered
230   
231   **Selection of Representative:**
232   
233   For each over-represented family, the first target alphabetically is kept as the representative.
234   All other family members are excluded from the targets DataFrame.
235   
236   :param targets_df: DataFrame of drug targets for a specific drug-disease combination.
237                      Must contain a ``target_gene_name`` column with gene symbols.
238   :type targets_df: DataFrame
239   :param drug_chembl_id: ChEMBL identifier for the drug being analyzed (e.g., 'CHEMBL1200')
240   :type drug_chembl_id: str
241   :param mesh_heading: MeSH disease heading for this drug-disease combination (e.g., 'Lung Neoplasms')
242   :type mesh_heading: str
243   :param lookup_table: Pre-computed family representation lookup from :func:`lookup_drug_family_representation`.
244                        Used to identify which families are over-represented for this drug-disease pair.
245   :type lookup_table: DataFrame
246   :param pct_threshold: Percentage threshold for filtering large families (0-100). Families with
247                         representation above this threshold are filtered (if other criteria met).
248                         Lower values = more aggressive filtering. Defaults to 40.0.
249   :type pct_threshold: float, optional
250   :param small_family_threshold: Minimum family size for percentage-based filtering. Families smaller
251                                  than this are exempt from filtering to preserve biological specificity.
252                                  Defaults to 5.
253   :type small_family_threshold: int, optional
254   :param min_targets_threshold: Minimum number of targets in a family to trigger filtering. Families with
255                                 fewer targets lack sufficient redundancy to justify filtering. Defaults to 3.
256   :type min_targets_threshold: int, optional
257   :return: Filtered DataFrame with same columns as input. Over-represented family members (except
258            the first alphabetical representative) are removed. Returns input unchanged if no filtering
259            is needed or if drug/disease not found in lookup table.
260   :rtype: DataFrame
261   
262   .. note::
263      **Alphabetical Selection**: Using the first alphabetical target as representative ensures
264      deterministic, reproducible filtering across runs.
265   
266   .. note::
267      **No Family Assignment**: Targets that don't map to HGNC gene families in the lookup table
268      pass through unfiltered. These are typically novel targets or non-gene targets.
269   
270   .. warning::
271      If ``drug_chembl_id`` and ``mesh_heading`` combination is not found in the lookup table,
272      the function returns the input unfiltered and prints a warning. This can occur if the
273      drug has no targets with family assignments.
274   
275   .. warning::
276      **Aggressive Filtering**: Lowering ``pct_threshold`` or ``min_targets_threshold`` makes
277      filtering more aggressive but may remove biologically important targets. Default values
278      balance bias reduction with information retention.
279   
280   **Examples**
281   
282   Basic usage with default thresholds::
283   
284      >>> lookup = lookup_drug_family_representation(chembl_data)
285      >>> filtered = filter_overrepresented_families(
286      ...     targets_df=drug_targets,
287      ...     drug_chembl_id='CHEMBL1200',  # Sunitinib
288      ...     mesh_heading='Renal Cell Carcinoma',
289      ...     lookup_table=lookup
290      ... )
291      >>> print(f"Filtered {len(drug_targets) - len(filtered)} targets")
292   
293   More aggressive filtering (lower thresholds)::
294   
295      >>> filtered = filter_overrepresented_families(
296      ...     targets_df,
297      ...     drug_chembl_id='CHEMBL123',
298      ...     mesh_heading='Lung Neoplasms',
299      ...     lookup_table=lookup,
300      ...     pct_threshold=25.0,  # Filter families >25% represented
301      ...     min_targets_threshold=2  # Filter families with >2 targets
302      ... )
303   
304   Preserve small families (e.g., HDACs) from filtering::
305   
306      >>> # Don't filter families with <8 genes
307      >>> filtered = filter_overrepresented_families(
308      ...     targets_df,
309      ...     drug_chembl_id='CHEMBL123',
310      ...     mesh_heading='Lymphoma',
311      ...     lookup_table=lookup,
312      ...     small_family_threshold=8
313      ... )
314   
315   Inspect which families were filtered::
316   
317      >>> original_targets = set(drug_targets['target_gene_name'])
318      >>> filtered_targets = set(filtered['target_gene_name'])
319      >>> removed = original_targets - filtered_targets
320      >>> print(f"Removed targets: {removed}")
321   
322   .. seealso::
323      - :func:`lookup_drug_family_representation`: Creates the required lookup table
324      - :func:`compute`: Integrates this filtering into the clinical score calculation pipeline
325   """
326   # Query lookup table for this drug's family representation
327   drug_families = lookup_table[
328      (lookup_table['drug_chembl_id'] == drug_chembl_id) &
329      (lookup_table['mesh_heading'] == mesh_heading)
330   ].copy()
331   
332   # Handle case where drug has no family data
333   if len(drug_families) == 0:
334      print(f"Warning: No family data found for drug {drug_chembl_id} and mesh {mesh_heading}, skipping filtering")
335      return targets_df
336   
337   # Identify over-represented families based on thresholds
338   # Family must be: 
339   # (1) large enough (≥ small_family_threshold) AND 
340   # (2) has enough targets from this drug (> min_targets_threshold) AND
341   # (3) EITHER highly represented (> pct_threshold%) OR has many absolute targets
342   #
343   # The last condition handles two cases:
344   # - High percentage: drug targets most of a small/medium family (e.g., 80% of 10 genes)
345   # - Many targets: drug targets a subfamily within a large family (e.g., 4 NDUFAF genes 
346   #   within a 72-gene family, which is only 5.6% but still represents redundancy)
347   #
348   # NOTE: When absolute threshold is met (targets >= 4), bypass small family check
349   # to handle cases like HDAC class IIA (4 genes total, all 4 targeted = 100%)
350   overrep_families = drug_families[
351      (drug_families['targets_in_family'] > min_targets_threshold) &
352      (
353         # Case 1: Large family with high percentage representation
354         (
355            (drug_families['total_genes_in_family'] >= small_family_threshold) &
356            (drug_families['pct_represented'] > pct_threshold)
357         ) |
358         # Case 2: Absolute threshold met (bypass small family check)
359         (drug_families['targets_in_family'] >= 4)
360      )
361   ]
362   
363   if len(overrep_families) > 0:
364      # Report filtering context
365      print(f"Found {len(overrep_families)} over-represented families for drug {drug_chembl_id} in {mesh_heading} (>{pct_threshold}%)")
366      
367      # Extract and deduplicate target lists from over-represented families
368      # Use list comprehension for better performance than iterrows()
369      all_family_targets = [
370         row['target_genes_in_family'].split(', ')
371         for _, row in overrep_families.iterrows()
372      ]
373      
374      # Keep first alphabetical target from each family as representative
375      targets_to_keep = [
376         sorted(targets)[0] for targets in all_family_targets
377      ]
378      
379      # Build exclusion list (all family targets except the kept representatives)
380      targets_to_exclude = [
381         target 
382         for targets in all_family_targets 
383         for target in targets 
384         if target not in targets_to_keep
385      ]
386      
387      print(f"  Excluding {len(targets_to_exclude)} targets, keeping {len(targets_to_keep)} (first alphabetical from each family)")
388      
389      # Apply filter and return
390      return targets_df[
391         ~targets_df['target_gene_name'].isin(targets_to_exclude)
392      ].reset_index(drop=True)
393   
394   # No filtering needed - return input unchanged
395   return targets_df

1052def unique(scores: list, overlap = 1, common_genes = []):
1053    """
1054    Remove cross-disease overlapping genes from clinical scores to create disease-specific target lists.
1055    
1056    When analyzing multiple diseases, some genes appear as targets across many indications (e.g., EGFR
1057    in lung, breast, and colorectal cancers). For disease-specific analyses or comparative studies,
1058    removing these "pan-disease" targets reveals indication-specific biology.
1059    
1060    This function identifies genes present in multiple disease DataFrames and removes them from all
1061    DataFrames, leaving only targets unique to fewer diseases.
1062    
1063    **Overlap Logic:**
1064    
1065    A gene is considered "overlapping" and removed if:
1066    
1067    - It appears in **more than** ``overlap`` diseases (not in all), OR
1068    - It's in the ``common_genes`` list (manual specification)
1069    
1070    **Important**: Genes appearing in **all** diseases are NOT removed by default (see Notes).
1071    
1072    **Use Cases:**
1073    
1074    - Identify disease-specific therapeutic targets
1075    - Remove broadly-acting targets (e.g., TP53, EGFR) for specialized analysis
1076    - Create non-overlapping target sets for comparative drug discovery
1077    - Study indication-specific versus pan-cancer biology
1078    
1079    :param scores: List of DataFrames containing clinical scores. Each must have a ``Target Gene``
1080                   column with HGNC gene symbols. Typically one DataFrame per disease.
1081    :type scores: list[DataFrame]
1082    :param overlap: Minimum number of diseases a gene must appear in to be removed (exclusive of all-disease genes).
1083                    - ``overlap=1``: Remove genes in 2+ diseases (but not all)
1084                    - ``overlap=2``: Remove genes in 3+ diseases (but not all)
1085                    - ``overlap=N``: Remove genes in (N+1)+ diseases (but not all)
1086                    Defaults to 1.
1087    :type overlap: int, optional
1088    :param common_genes: List of gene symbols to always remove regardless of overlap count.
1089                         Useful for manually excluding known pan-disease targets. Defaults to [].
1090    :type common_genes: list[str], optional
1091    :return: List of DataFrames (same length/order as input) with overlapping genes removed.
1092             Original DataFrames are **not** modified (returns filtered copies).
1093    :rtype: list[DataFrame]
1094
1095    .. note::
1096       **All-Disease Genes Preserved**: By default, genes appearing in ALL input DataFrames are
1097       NOT removed. This preserves truly universal targets. To change this behavior, modify the
1098       comparison in :func:`find_overlapping_genes` from ``<`` to ``<=``.
1099
1100    .. note::
1101       This function uses :func:`find_overlapping_genes` internally to identify overlap.
1102
1103    .. warning::
1104       This function creates copies of input DataFrames. Originals remain unchanged.
1105
1106    .. warning::
1107       Setting ``overlap`` too low (e.g., 0) removes genes present in only 1+ diseases,
1108       potentially eliminating most of your data. Use with caution.
1109
1110    **Examples**
1111    
1112    Remove genes appearing in 2+ diseases (but not all 7)::
1113    
1114       >>> from alethiotx.artemis.clinical import load, unique
1115       >>> breast, lung, prostate, melanoma, bowel, diabetes, cardio = load()
1116       >>> scores_list = [breast, lung, prostate, melanoma, bowel, diabetes, cardio]
1117       >>> unique_scores = unique(scores_list, overlap=1)
1118       >>> breast_unique, lung_unique, *rest = unique_scores
1119       >>> print(f"Breast: {len(breast)} -> {len(breast_unique)} unique targets")
1120    
1121    More aggressive filtering (remove genes in 3+ diseases)::
1122    
1123       >>> unique_scores = unique(scores_list, overlap=2)
1124    
1125    Manually exclude known pan-cancer targets::
1126    
1127       >>> pan_cancer_genes = ['TP53', 'EGFR', 'KRAS', 'PIK3CA']
1128       >>> unique_scores = unique(scores_list, overlap=1, common_genes=pan_cancer_genes)
1129    
1130    Compare disease-specific vs all targets::
1131    
1132       >>> all_scores = load()
1133       >>> unique_scores = unique(list(all_scores), overlap=1)
1134       >>> for i, disease in enumerate(['Breast', 'Lung', 'Prostate']):
1135       ...     original = len(all_scores[i])
1136       ...     unique_count = len(unique_scores[i])
1137       ...     removed = original - unique_count
1138       ...     print(f"{disease}: removed {removed} overlapping genes ({100*removed/original:.1f}%)")
1139    
1140    Identify what genes were removed::
1141    
1142       >>> original_breast_genes = set(breast['Target Gene'])
1143       >>> unique_breast_genes = set(unique_scores[0]['Target Gene'])
1144       >>> removed_genes = original_breast_genes - unique_breast_genes
1145       >>> print(f"Removed from breast: {sorted(removed_genes)}")
1146    
1147    .. seealso::
1148       - :func:`find_overlapping_genes`: Identifies which genes overlap across disease lists
1149       - :func:`load`: Load clinical scores for multiple diseases
1150    """
1151    genes = []
1152    for n, d in enumerate(scores):
1153        genes.append(d['Target Gene'].tolist())
1154
1155    overlapping_genes = find_overlapping_genes(genes, overlap = overlap, common_genes = common_genes)
1156    
1157    res = scores.copy()
1158
1159    for n, d in enumerate(res):
1160        res[n] = d[~d['Target Gene'].isin(overlapping_genes)]
1161    
1162    return(res)

 984def approved(scores: list):
 985   """
 986   Filter clinical score DataFrames to include only targets with approved drugs (Phase 4).
 987
 988   This utility function processes a list of clinical score DataFrames and retains only
 989   targets that have at least one approved drug (phase 4 trials). This is useful for
 990   focusing analysis on clinically validated targets with FDA/EMA-approved therapeutics.
 991   
 992   **Use Cases:**
 993   
 994   - Benchmark new targets against established approved targets
 995   - Study successful drug development patterns
 996   - Identify approved drugs for repurposing opportunities
 997   - Filter out experimental/developmental targets
 998
 999   :param scores: List of DataFrames containing clinical scores. Each DataFrame must have
1000                  an ``approved`` column (boolean) indicating whether the target has any
1001                  approved (phase 4) drugs. Typically obtained from :func:`load` or :func:`compute`.
1002   :type scores: list[DataFrame]
1003   :return: List of DataFrames (same length and order as input) with only approved targets.
1004            Original DataFrames are **not** modified (returns filtered copies).
1005   :rtype: list[DataFrame]
1006
1007   .. note::
1008      This function creates copies of the input DataFrames. Original DataFrames remain unchanged.
1009
1010   .. note::
1011      Empty DataFrames are returned for diseases with no approved targets (e.g., some rare cancers).
1012
1013   **Examples**
1014   
1015   Filter pre-computed scores to approved only::
1016   
1017      >>> from alethiotx.artemis.clinical import load, approved
1018      >>> breast, lung, prostate, melanoma, bowel, diabetes, cardio = load()
1019      >>> breast_app, lung_app, prostate_app, melanoma_app, bowel_app, diabetes_app, cardio_app = approved(
1020      ...     [breast, lung, prostate, melanoma, bowel, diabetes, cardio]
1021      ... )
1022      >>> print(f"Approved breast cancer targets: {len(breast_app)}")
1023   
1024   Compare approved vs all targets::
1025   
1026      >>> all_scores = load()
1027      >>> approved_scores = approved(list(all_scores))
1028      >>> for i, disease in enumerate(['Breast', 'Lung', 'Prostate', 'Melanoma', 'Bowel', 'Diabetes', 'Cardio']):
1029      ...     all_count = len(all_scores[i])
1030      ...     approved_count = len(approved_scores[i])
1031      ...     pct = 100 * approved_count / all_count if all_count > 0 else 0
1032      ...     print(f"{disease}: {approved_count}/{all_count} ({pct:.1f}%) approved")
1033   
1034   Identify top approved targets::
1035   
1036      >>> breast, lung, prostate, melanoma, bowel, diabetes, cardio = load()
1037      >>> breast_approved = approved([breast])[0]
1038      >>> top_5 = breast_approved.nlargest(5, 'Clinical Score')
1039      >>> print(top_5[['Target Gene', 'Clinical Score']])
1040   
1041   .. seealso::
1042      - :func:`load`: Load pre-computed scores (returns all targets)
1043      - :func:`compute`: Compute scores with ``approved`` column
1044   """
1045   res = scores.copy()
1046
1047   for n, d in enumerate(res):
1048      res[n] = d[d['approved'] == True]
1049
1050   return(res)