{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":118765,"databundleVersionId":15231210,"isSourceIdPinned":false,"sourceType":"competition"},{"sourceId":14467782,"sourceType":"datasetVersion","datasetId":9241051}],"dockerImageVersionId":30746,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install biopython","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:14:52.890315Z","iopub.execute_input":"2026-02-04T16:14:52.890808Z","iopub.status.idle":"2026-02-04T16:15:08.534288Z","shell.execute_reply.started":"2026-02-04T16:14:52.890751Z","shell.execute_reply":"2026-02-04T16:15:08.532935Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"In part 1, I did some basic stat t with R. Now I will try some better understanding with python.\nPart 1: https://www.kaggle.com/code/parsahriri/dataanalysiswithr","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.linalg import cholesky, LinAlgError\nimport time\nimport warnings\nimport os, sys\nfrom typing import Dict, Tuple, List\nimport re\nfrom Bio.PDB import MMCIFParser\nfrom Bio import SeqIO\nfrom mpl_toolkits.mplot3d import Axes3D\nimport glob\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nfrom sklearn.decomposition import PCA\nfrom sklearn.manifold import TSNE\nimport umap\nfrom tqdm import tqdm","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:15:08.536079Z","iopub.execute_input":"2026-02-04T16:15:08.536528Z","iopub.status.idle":"2026-02-04T16:15:10.676202Z","shell.execute_reply.started":"2026-02-04T16:15:08.536485Z","shell.execute_reply":"2026-02-04T16:15:10.674662Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# load the data","metadata":{}},{"cell_type":"code","source":"train_seqs = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/train_sequences.csv')\nvalid_seqs = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/validation_sequences.csv')\ntest_seqs = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/test_sequences.csv')\ntrain_labels = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/train_labels.csv')\nvalid_labels = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/validation_labels.csv')\n\nrna_metadata = pd.read_csv(\"/kaggle/input/stanford-rna-3d-folding-2/extra/rna_metadata.csv\")\n\n\ncif_folder = \"/kaggle/input/stanford-rna-3d-folding-2/PDB_RNA/\"\nmsa_folder = \"/kaggle/input/stanford-rna-3d-folding-2/MSA/\"\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:15:10.679688Z","iopub.execute_input":"2026-02-04T16:15:10.680317Z","iopub.status.idle":"2026-02-04T16:15:28.154865Z","shell.execute_reply.started":"2026-02-04T16:15:10.680288Z","shell.execute_reply":"2026-02-04T16:15:28.153278Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_seqs","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:15:28.156957Z","iopub.execute_input":"2026-02-04T16:15:28.157449Z","iopub.status.idle":"2026-02-04T16:15:28.198009Z","shell.execute_reply.started":"2026-02-04T16:15:28.157404Z","shell.execute_reply":"2026-02-04T16:15:28.196063Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_seqs.describe()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:15:28.199939Z","iopub.execute_input":"2026-02-04T16:15:28.200494Z","iopub.status.idle":"2026-02-04T16:15:28.28767Z","shell.execute_reply.started":"2026-02-04T16:15:28.200388Z","shell.execute_reply":"2026-02-04T16:15:28.285576Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"rna_metadata.describe()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:15:28.290188Z","iopub.execute_input":"2026-02-04T16:15:28.29117Z","iopub.status.idle":"2026-02-04T16:15:28.399865Z","shell.execute_reply.started":"2026-02-04T16:15:28.291121Z","shell.execute_reply":"2026-02-04T16:15:28.398293Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"rna_metadata.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:15:28.401452Z","iopub.execute_input":"2026-02-04T16:15:28.402171Z","iopub.status.idle":"2026-02-04T16:15:28.439655Z","shell.execute_reply.started":"2026-02-04T16:15:28.402141Z","shell.execute_reply":"2026-02-04T16:15:28.436906Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# plot one sample","metadata":{}},{"cell_type":"code","source":"cif_path = \"/kaggle/input/stanford-rna-3d-folding-2/PDB_RNA/100d.cif\"\n\nparser = MMCIFParser(QUIET=True)\nstructure = parser.get_structure(\"RNA\", cif_path)\n\n# Collect phosphorus (P) atoms for backbone\nxs, ys, zs = [], [], []\nfor model in structure:\n    for chain in model:\n        for residue in chain:\n            if \"P\" in residue:\n                atom = residue[\"P\"]\n                x, y, z = atom.get_coord()\n                xs.append(x)\n                ys.append(y)\n                zs.append(z)\n\n# 3D plot\nfig = plt.figure(figsize=(6,6))\nax = fig.add_subplot(111, projection=\"3d\")\nax.plot(xs, ys, zs, linewidth=1)\nax.scatter(xs, ys, zs, s=5)\nax.set_xlabel(\"X\")\nax.set_ylabel(\"Y\")\nax.set_zlabel(\"Z\")\nax.set_title(\"RNA 3D Structure (100d)\")\nplt.show()\n\n\nmsa_path = \"/kaggle/input/stanford-rna-3d-folding-2/MSA/157D.MSA.fasta\"\n\nmsa_sequences = list(SeqIO.parse(msa_path, \"fasta\"))\n\nprint(f\"Number of sequences in MSA: {len(msa_sequences)}\\n\")\n\n# Print first 5 sequences\nfor seq_record in msa_sequences[:5]:\n    print(f\">{seq_record.id}\")\n    print(seq_record.seq)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:15:28.442519Z","iopub.execute_input":"2026-02-04T16:15:28.443453Z","iopub.status.idle":"2026-02-04T16:15:28.87907Z","shell.execute_reply.started":"2026-02-04T16:15:28.443398Z","shell.execute_reply":"2026-02-04T16:15:28.87786Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# How to use parser\n### they have provided a function as you can see below for parsing fast files","metadata":{}},{"cell_type":"code","source":"def parse_fasta(fasta_content: str) -> Dict[str, Tuple[str, List[str]]]:\n    \"\"\"\n    Parse FASTA content into dictionary.\n\n    Args:\n        fasta_content: Multi-line FASTA string with format:\n        >1A1T_1|Chain A[auth B]|SL3 STEM-LOOP RNA|\n        or\n        >104D_1|Chains A[auth A], B[auth B]|DNA/RNA (...)|\n\n    Returns:\n        Dictionary mapping auth chain_id to (sequence, list_of_auth_chain_ids)\n        Example: {\"A\": (\"ACGT\", [\"A\", \"B\"]), \"C\": (\"UGCA\", [\"C\"])}\n        The key is the auth chain ID, and the list contains all auth chain IDs for this sequence\n    \"\"\"\n    result = {}\n    lines = fasta_content.strip().split(\"\\n\")\n\n    i = 0\n    while i < len(lines):\n        line = lines[i].strip()\n\n        if line.startswith(\">\"):\n            # Parse new format header: >104D_1|Chains A[auth A], B[auth B]|...| or >1A1T_1|Chain A[auth B]|...|\n            # Extract the chains part (between first | and second |)\n            parts = line.split(\"|\")\n            if len(parts) < 2:\n                print(\"Warning: Malformed FASTA header:\", line)\n                auth_chain_ids = []\n                chains_part = \"\"\n            else:\n                chains_part = parts[1].strip()\n\n                # Extract auth chain IDs from patterns like \"Chain A[auth B]\" or \"Chains A[auth A], B[auth B] or just \"Chain A\" or \"Chains A, B\"\n                auth_chain_ids = []\n                replaced_chains_part = re.sub(r\"^Chains? \", \"\", chains_part)\n                chains = replaced_chains_part.split(\",\")\n                for chain in chains:\n                    auth_match = re.search(r\"\\[auth ([^\\]]+)\\]\", chain)\n                    if auth_match:\n                        auth_chain_ids.append(auth_match.group(1).strip())\n                    else:\n                        c = chain.strip()\n                        if c:\n                            auth_chain_ids.append(c)\n\n            if not auth_chain_ids:\n                print(\"Warning: Empty chains part:\", chains_part)\n                primary_auth_chain = None\n            else:\n                # Use the first auth chain ID as the key\n                primary_auth_chain = auth_chain_ids[0]\n\n            # Read sequence (next lines until next header or end)\n            sequence = \"\"\n            while (i + 1) < len(lines) and lines[i + 1].startswith(\">\") is False:\n                sequence += lines[i + 1].strip()\n                i += 1\n            result[primary_auth_chain] = (sequence, auth_chain_ids)\n\n        i += 1\n\n    return result","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:15:28.882865Z","iopub.execute_input":"2026-02-04T16:15:28.883424Z","iopub.status.idle":"2026-02-04T16:15:28.896681Z","shell.execute_reply.started":"2026-02-04T16:15:28.883374Z","shell.execute_reply":"2026-02-04T16:15:28.895099Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"path = '/kaggle/input/stanford-rna-3d-folding-2/MSA/'\nfiles = os.listdir(path)\n\n\n\nfor i, file in enumerate(files[:2]):\n    with open(os.path.join(path, file), \"r\") as f:\n        content = f.read()\n        results = parse_fasta(content)\n        print(results)\n        first_id = list(results.keys())[0]\n        for key in results:\n            print(key)\n        try:\n            print(f\"File: {file} | Found ID: {first_id} | chain: {results['chain'][0]} | Length: {len(results[first_id][0])}\")\n        except:\n            print(f\"File: {file} | Found ID: {first_id}  | Length: {len(results[first_id][0])}\")\n            ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:16:38.551127Z","iopub.execute_input":"2026-02-04T16:16:38.552106Z","iopub.status.idle":"2026-02-04T16:16:38.619123Z","shell.execute_reply.started":"2026-02-04T16:16:38.552073Z","shell.execute_reply":"2026-02-04T16:16:38.618003Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# DR","metadata":{}},{"cell_type":"code","source":"# ============================================================\n# 1. Load all CIF files and extract backbone coordinates\n# ============================================================\n\nnp.random.seed(42)\ncif_files = glob.glob(os.path.join(cif_folder, \"*.cif\"))\ncif_files = np.random.choice(\n    cif_files,\n    size=max(1, int(0.03 * len(cif_files))),\n    replace=False\n)\nprint(f\"Found {len(cif_files)} CIF files\")\n\nparser = MMCIFParser(QUIET=True)\n\ndef extract_backbone_coords(cif_path):\n    \"\"\"Extract phosphorus (P) atom coordinates from RNA structure.\"\"\"\n    try:\n        structure = parser.get_structure(\"RNA\", cif_path)\n        coords = []\n        for model in structure:\n            for chain in model:\n                for residue in chain:\n                    if \"P\" in residue:\n                        atom = residue[\"P\"]\n                        coords.append(atom.get_coord())\n        return np.array(coords) if coords else None\n    except Exception as e:\n        print(f\"Error parsing {cif_path}: {e}\")\n        return None\n\n# Extract coordinates from all files\nall_structures = {}\nall_features = []\nstructure_names = []\n\nprint(\"Loading CIF files...\")\nfor cif_path in tqdm(cif_files):\n    name = os.path.basename(cif_path).replace(\".cif\", \"\")\n    coords = extract_backbone_coords(cif_path)\n    \n    if coords is not None and len(coords) > 0:\n        all_structures[name] = coords\n        \n        # Create feature vector: statistics of coordinates\n        features = [\n            len(coords),                    # Number of residues\n            coords[:, 0].mean(),            # Mean X\n            coords[:, 1].mean(),            # Mean Y\n            coords[:, 2].mean(),            # Mean Z\n            coords[:, 0].std(),             # Std X\n            coords[:, 1].std(),             # Std Y\n            coords[:, 2].std(),             # Std Z\n            np.linalg.norm(coords[-1] - coords[0]),  # End-to-end distance\n            np.mean(np.linalg.norm(np.diff(coords, axis=0), axis=1)),  # Mean bond length\n        ]\n        all_features.append(features)\n        structure_names.append(name)\n\nprint(f\"Successfully loaded {len(all_structures)} structures\")\n\n# Convert to numpy array\nX = np.array(all_features)\nprint(f\"Feature matrix shape: {X.shape}\")\n\n# ============================================================\n# 2. Apply Dimensionality Reduction\n# ============================================================\n\n# Standardize features\nfrom sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\nX_scaled = scaler.fit_transform(X)\n\n# --- PCA ---\nprint(\"\\nApplying PCA...\")\npca = PCA(n_components=2)\nX_pca = pca.fit_transform(X_scaled)\nprint(f\"PCA explained variance ratio: {pca.explained_variance_ratio_}\")\n\n# --- t-SNE ---\nprint(\"\\nApplying t-SNE...\")\ntsne = TSNE(n_components=2, random_state=42, perplexity=min(30, len(X)-1))\nX_tsne = tsne.fit_transform(X_scaled)\n\n# --- UMAP ---\nprint(\"\\nApplying UMAP...\")\nreducer = umap.UMAP(n_components=2, random_state=42)\nX_umap = reducer.fit_transform(X_scaled)\n\n# ============================================================\n# 3. Visualization\n# ============================================================\n\nfig, axes = plt.subplots(1, 3, figsize=(15, 5))\n\n# PCA plot\naxes[0].scatter(X_pca[:, 0], X_pca[:, 1], c=X[:, 0], cmap='viridis', alpha=0.7, s=20)\naxes[0].set_xlabel(\"PC1\")\naxes[0].set_ylabel(\"PC2\")\naxes[0].set_title(f\"PCA (Var: {sum(pca.explained_variance_ratio_):.2%})\")\n\n# t-SNE plot\naxes[1].scatter(X_tsne[:, 0], X_tsne[:, 1], c=X[:, 0], cmap='viridis', alpha=0.7, s=20)\naxes[1].set_xlabel(\"t-SNE 1\")\naxes[1].set_ylabel(\"t-SNE 2\")\naxes[1].set_title(\"t-SNE\")\n\n# UMAP plot\nscatter = axes[2].scatter(X_umap[:, 0], X_umap[:, 1], c=X[:, 0], cmap='viridis', alpha=0.7, s=20)\naxes[2].set_xlabel(\"UMAP 1\")\naxes[2].set_ylabel(\"UMAP 2\")\naxes[2].set_title(\"UMAP\")\n\nplt.colorbar(scatter, ax=axes[2], label='Num Residues')\nplt.tight_layout()\nplt.savefig(\"dimensionality_reduction_comparison.png\", dpi=150)\nplt.show()\n\n# ============================================================\n# 4. Load and Process MSA Files\n# ============================================================\n\nmsa_files = glob.glob(os.path.join(msa_folder, \"*.fasta\"))\nmsa_files = np.random.choice(\n    msa_files,\n    size=max(1, int(0.1 * len(msa_files))),\n    replace=False\n)\n\nprint(f\"\\nFound {len(msa_files)} MSA files\")\n\ndef msa_to_features(msa_path, max_seqs=100):\n    \"\"\"Convert MSA to feature vector using one-hot encoding statistics.\"\"\"\n    try:\n        sequences = list(SeqIO.parse(msa_path, \"fasta\"))[:max_seqs]\n        if not sequences:\n            return None\n        \n        # Get sequence length and count\n        n_seqs = len(sequences)\n        seq_len = len(sequences[0].seq)\n        \n        # Calculate conservation (frequency of most common nucleotide per position)\n        nucleotide_map = {'A': 0, 'C': 1, 'G': 2, 'U': 3, 'T': 3, '-': 4}\n        \n        conservation_scores = []\n        for pos in range(min(seq_len, 500)):  # Limit to first 500 positions\n            counts = np.zeros(5)\n            for seq in sequences:\n                if pos < len(seq.seq):\n                    nuc = str(seq.seq[pos]).upper()\n                    if nuc in nucleotide_map:\n                        counts[nucleotide_map[nuc]] += 1\n            conservation_scores.append(counts.max() / n_seqs)\n        \n        features = [\n            n_seqs,\n            seq_len,\n            np.mean(conservation_scores),\n            np.std(conservation_scores),\n            np.min(conservation_scores),\n            np.max(conservation_scores),\n        ]\n        return features\n    except Exception as e:\n        print(f\"Error processing {msa_path}: {e}\")\n        return None\n\n# Process MSA files\nmsa_features = []\nmsa_names = []\n\nprint(\"Processing MSA files...\")\nfor msa_path in tqdm(msa_files[:100]):  # Limit for speed\n    name = os.path.basename(msa_path).replace(\".MSA.fasta\", \"\")\n    features = msa_to_features(msa_path)\n    if features:\n        msa_features.append(features)\n        msa_names.append(name)\n\nif msa_features:\n    X_msa = np.array(msa_features)\n    X_msa_scaled = StandardScaler().fit_transform(X_msa)\n    \n    # Apply dimensionality reduction to MSA features\n    pca_msa = PCA(n_components=2)\n    X_msa_pca = pca_msa.fit_transform(X_msa_scaled)\n    \n    plt.figure(figsize=(8, 6))\n    plt.scatter(X_msa_pca[:, 0], X_msa_pca[:, 1], c=X_msa[:, 0], cmap='plasma', alpha=0.7)\n    plt.colorbar(label='Num Sequences')\n    plt.xlabel(\"PC1\")\n    plt.ylabel(\"PC2\")\n    plt.title(\"MSA Features - PCA\")\n    plt.savefig(\"msa_pca.png\", dpi=150)\n    plt.show()\n\nprint(\"\\nDone!\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:32:51.542483Z","iopub.execute_input":"2026-02-04T16:32:51.542959Z","iopub.status.idle":"2026-02-04T16:48:13.810789Z","shell.execute_reply.started":"2026-02-04T16:32:51.542929Z","shell.execute_reply":"2026-02-04T16:48:13.809482Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_seqs = train_seqs.sample(frac=0.1, random_state=42).reset_index(drop=True)\nprint(f\"Using 10% sample: {train_seqs.shape[0]} samples\")\nprint(f\"Columns: {train_seqs.columns.tolist()}\")\n\n# ============================================================\n# 2. Feature Engineering from train_seqs\n# ============================================================\n\ndef extract_ligand_features(df):\n    \"\"\"Extract features from ligand information.\"\"\"\n    df = df.copy()\n    \n    # Count number of ligands\n    df['n_ligands'] = df['ligand_ids'].fillna('').apply(\n        lambda x: len(x.split(';')) if x else 0\n    )\n    \n    # Check for specific common ligands\n    df['has_MG'] = df['ligand_ids'].fillna('').str.contains('MG').astype(int)\n    df['has_ZN'] = df['ligand_ids'].fillna('').str.contains('ZN').astype(int)\n    df['has_K'] = df['ligand_ids'].fillna('').str.contains('K').astype(int)\n    df['has_SPM'] = df['ligand_ids'].fillna('').str.contains('SPM').astype(int)\n    \n    return df\n\ndef extract_sequence_features(df):\n    \"\"\"Extract features from RNA sequences.\"\"\"\n    df = df.copy()\n    \n    df['seq_length'] = df['sequence'].fillna('').apply(len)\n    \n    # Nucleotide composition\n    for nuc in ['A', 'C', 'G', 'U']:\n        df[f'frac_{nuc}'] = df['sequence'].fillna('').apply(\n            lambda x: x.count(nuc) / len(x) if len(x) > 0 else 0\n        )\n    \n    # GC content\n    df['gc_content'] = df['frac_G'] + df['frac_C']\n    \n    # Purine/Pyrimidine ratio\n    df['purine_frac'] = df['frac_A'] + df['frac_G']\n    \n    return df\n\ndef extract_stoichiometry_features(df):\n    \"\"\"Extract features from stoichiometry.\"\"\"\n    df = df.copy()\n    \n    # Count number of chains\n    df['n_chains'] = df['stoichiometry'].fillna('').apply(\n        lambda x: len(x.split(';')) if x else 0\n    )\n    \n    # Total copy number\n    def get_total_copies(stoich):\n        if pd.isna(stoich) or stoich == '':\n            return 0\n        total = 0\n        parts = stoich.split(';')\n        for part in parts:\n            match = re.search(r':(\\d+)', part)\n            if match:\n                total += int(match.group(1))\n        return total\n    \n    df['total_copies'] = df['stoichiometry'].apply(get_total_copies)\n    \n    return df\n\ndef extract_temporal_features(df):\n    \"\"\"Extract year from temporal_cutoff.\"\"\"\n    df = df.copy()\n    df['year'] = pd.to_datetime(df['temporal_cutoff'], errors='coerce').dt.year\n    df['year'] = df['year'].fillna(df['year'].median())\n    return df\n\n# Apply all feature extractions\ntrain_seqs = extract_ligand_features(train_seqs)\ntrain_seqs = extract_sequence_features(train_seqs)\ntrain_seqs = extract_stoichiometry_features(train_seqs)\ntrain_seqs = extract_temporal_features(train_seqs)\n\nprint(\"\\nExtracted features:\")\nprint(train_seqs.head())\n\n# ============================================================\n# 3. Load CIF files and merge with metadata\n# ============================================================\n\ncif_folder = \"/kaggle/input/stanford-rna-3d-folding-2/PDB_RNA/\"\ncif_files = glob.glob(os.path.join(cif_folder, \"*.cif\"))\n\n# Filter CIF files to only those matching sampled target_ids\nsampled_ids = set(train_seqs['target_id'].str.upper().tolist())\ncif_files = [f for f in cif_files if os.path.basename(f).replace(\".cif\", \"\").upper() in sampled_ids]\nprint(f\"\\nUsing {len(cif_files)} CIF files (10% sample)\")\n\nparser = MMCIFParser(QUIET=True)\n\ndef extract_structure_features(cif_path):\n    \"\"\"Extract 3D structural features from CIF file.\"\"\"\n    try:\n        structure = parser.get_structure(\"RNA\", cif_path)\n        coords = []\n        for model in structure:\n            for chain in model:\n                for residue in chain:\n                    if \"P\" in residue:\n                        coords.append(residue[\"P\"].get_coord())\n        \n        if not coords:\n            return None\n        \n        coords = np.array(coords)\n        \n        # Calculate structural features\n        features = {\n            'n_residues': len(coords),\n            'mean_x': coords[:, 0].mean(),\n            'mean_y': coords[:, 1].mean(),\n            'mean_z': coords[:, 2].mean(),\n            'std_x': coords[:, 0].std(),\n            'std_y': coords[:, 1].std(),\n            'std_z': coords[:, 2].std(),\n            'end_to_end_dist': np.linalg.norm(coords[-1] - coords[0]),\n            'mean_bond_length': np.mean(np.linalg.norm(np.diff(coords, axis=0), axis=1)),\n            'radius_of_gyration': np.sqrt(np.mean(np.sum((coords - coords.mean(axis=0))**2, axis=1))),\n            'max_span': np.max(np.linalg.norm(coords[:, None] - coords[None, :], axis=2)),\n        }\n        return features\n    except Exception as e:\n        return None\n\n# Extract structural features\nprint(\"Extracting structural features from CIF files...\")\nstructure_data = []\nfor cif_path in tqdm(cif_files):\n    target_id = os.path.basename(cif_path).replace(\".cif\", \"\").upper()\n    features = extract_structure_features(cif_path)\n    if features:\n        features['target_id'] = target_id\n        structure_data.append(features)\n\nstructure_df = pd.DataFrame(structure_data)\nprint(f\"Extracted features for {len(structure_df)} structures\")\n\n# ============================================================\n# 4. Merge structural features with metadata\n# ============================================================\n\n# Normalize target_id for merging\ntrain_seqs['target_id_upper'] = train_seqs['target_id'].str.upper()\nstructure_df['target_id_upper'] = structure_df['target_id'].str.upper()\n\nmerged_df = pd.merge(\n    structure_df, \n    train_seqs, \n    left_on='target_id_upper', \n    right_on='target_id_upper', \n    how='inner'\n)\nprint(f\"\\nMerged dataset: {len(merged_df)} samples\")\n\n# ============================================================\n# 5. Create class labels for visualization\n# ============================================================\n\n# Create ligand-based classes\ndef get_primary_ligand(ligand_ids):\n    if pd.isna(ligand_ids) or ligand_ids == '':\n        return 'None'\n    ligands = ligand_ids.split(';')\n    return ligands[0] if ligands else 'None'\n\nmerged_df['primary_ligand'] = merged_df['ligand_ids'].apply(get_primary_ligand)\n\n# Create decade-based classes\nmerged_df['decade'] = (merged_df['year'] // 10 * 10).astype(int).astype(str) + 's'\n\n# Create sequence length classes\nmerged_df['length_class'] = pd.cut(\n    merged_df['seq_length'], \n    bins=[0, 50, 100, 200, 500, 10000],\n    labels=['<50', '50-100', '100-200', '200-500', '>500']\n)\n\n# Create GC content classes\nmerged_df['gc_class'] = pd.cut(\n    merged_df['gc_content'],\n    bins=[0, 0.4, 0.5, 0.6, 1.0],\n    labels=['Low (<40%)', 'Medium (40-50%)', 'High (50-60%)', 'Very High (>60%)']\n)\n\n# ============================================================\n# 6. Prepare features for dimensionality reduction\n# ============================================================\n\nfeature_cols = [\n    # Structural features\n    'n_residues', 'std_x', 'std_y', 'std_z', \n    'end_to_end_dist', 'mean_bond_length', 'radius_of_gyration', 'max_span',\n    # Sequence features\n    'seq_length', 'frac_A', 'frac_C', 'frac_G', 'frac_U', 'gc_content',\n    # Ligand features\n    'n_ligands', 'has_MG', 'has_ZN', 'has_K',\n    # Stoichiometry\n    'n_chains', 'total_copies',\n    # Temporal\n    'year'\n]\n\n# Filter to available columns\navailable_cols = [c for c in feature_cols if c in merged_df.columns]\nprint(f\"\\nUsing {len(available_cols)} features: {available_cols}\")\n\nX = merged_df[available_cols].fillna(0).values\n\n# Standardize\nscaler = StandardScaler()\nX_scaled = scaler.fit_transform(X)\n\n# ============================================================\n# 7. Apply Dimensionality Reduction\n# ============================================================\n\nprint(\"\\nApplying dimensionality reduction...\")\n\n# PCA\npca = PCA(n_components=2)\nX_pca = pca.fit_transform(X_scaled)\nmerged_df['pca_1'] = X_pca[:, 0]\nmerged_df['pca_2'] = X_pca[:, 1]\n\n# t-SNE\nperplexity = min(30, len(X) - 1)\ntsne = TSNE(n_components=2, random_state=42, perplexity=perplexity)\nX_tsne = tsne.fit_transform(X_scaled)\nmerged_df['tsne_1'] = X_tsne[:, 0]\nmerged_df['tsne_2'] = X_tsne[:, 1]\n\n# UMAP\nreducer = umap.UMAP(n_components=2, random_state=42, n_neighbors=min(15, len(X)-1))\nX_umap = reducer.fit_transform(X_scaled)\nmerged_df['umap_1'] = X_umap[:, 0]\nmerged_df['umap_2'] = X_umap[:, 1]\n\n# ============================================================\n# 8. Visualization Functions\n# ============================================================\n\ndef plot_by_class(df, x_col, y_col, class_col, title, ax, max_classes=10):\n    \"\"\"Plot dimensionality reduction colored by class.\"\"\"\n    classes = df[class_col].value_counts().head(max_classes).index.tolist()\n    df_plot = df[df[class_col].isin(classes)]\n    \n    colors = plt.cm.tab10(np.linspace(0, 1, len(classes)))\n    \n    for i, cls in enumerate(classes):\n        mask = df_plot[class_col] == cls\n        ax.scatter(\n            df_plot.loc[mask, x_col], \n            df_plot.loc[mask, y_col],\n            c=[colors[i]], label=str(cls), alpha=0.6, s=20\n        )\n    \n    ax.set_xlabel(x_col)\n    ax.set_ylabel(y_col)\n    ax.set_title(title)\n    ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=8)\n\n# ============================================================\n# 9. Plot by Primary Ligand\n# ============================================================\n\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\nplot_by_class(merged_df, 'pca_1', 'pca_2', 'primary_ligand', 'PCA by Ligand', axes[0])\nplot_by_class(merged_df, 'tsne_1', 'tsne_2', 'primary_ligand', 't-SNE by Ligand', axes[1])\nplot_by_class(merged_df, 'umap_1', 'umap_2', 'primary_ligand', 'UMAP by Ligand', axes[2])\nplt.suptitle('Dimensionality Reduction by Primary Ligand', fontsize=14)\nplt.tight_layout()\nplt.savefig('dim_reduction_by_ligand.png', dpi=150, bbox_inches='tight')\nplt.show()\n\n# ============================================================\n# 10. Plot by Decade\n# ============================================================\n\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\nplot_by_class(merged_df, 'pca_1', 'pca_2', 'decade', 'PCA by Decade', axes[0])\nplot_by_class(merged_df, 'tsne_1', 'tsne_2', 'decade', 't-SNE by Decade', axes[1])\nplot_by_class(merged_df, 'umap_1', 'umap_2', 'decade', 'UMAP by Decade', axes[2])\nplt.suptitle('Dimensionality Reduction by Publication Decade', fontsize=14)\nplt.tight_layout()\nplt.savefig('dim_reduction_by_decade.png', dpi=150, bbox_inches='tight')\nplt.show()\n\n# ============================================================\n# 11. Plot by Sequence Length Class\n# ============================================================\n\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\nplot_by_class(merged_df, 'pca_1', 'pca_2', 'length_class', 'PCA by Length', axes[0])\nplot_by_class(merged_df, 'tsne_1', 'tsne_2', 'length_class', 't-SNE by Length', axes[1])\nplot_by_class(merged_df, 'umap_1', 'umap_2', 'length_class', 'UMAP by Length', axes[2])\nplt.suptitle('Dimensionality Reduction by Sequence Length', fontsize=14)\nplt.tight_layout()\nplt.savefig('dim_reduction_by_length.png', dpi=150, bbox_inches='tight')\nplt.show()\n\n# ============================================================\n# 12. Plot by GC Content Class\n# ============================================================\n\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\nplot_by_class(merged_df, 'pca_1', 'pca_2', 'gc_class', 'PCA by GC Content', axes[0])\nplot_by_class(merged_df, 'tsne_1', 'tsne_2', 'gc_class', 't-SNE by GC Content', axes[1])\nplot_by_class(merged_df, 'umap_1', 'umap_2', 'gc_class', 'UMAP by GC Content', axes[2])\nplt.suptitle('Dimensionality Reduction by GC Content', fontsize=14)\nplt.tight_layout()\nplt.savefig('dim_reduction_by_gc.png', dpi=150, bbox_inches='tight')\nplt.show()\n\n# ============================================================\n# 13. Interactive Summary Plot with Multiple Class Views\n# ============================================================\n\nfig, axes = plt.subplots(2, 2, figsize=(14, 12))\n\n# UMAP colored by different features\nscatter1 = axes[0, 0].scatter(merged_df['umap_1'], merged_df['umap_2'], \n                               c=merged_df['seq_length'], cmap='viridis', alpha=0.6, s=15)\naxes[0, 0].set_title('UMAP - Sequence Length')\nplt.colorbar(scatter1, ax=axes[0, 0], label='Length')\n\nscatter2 = axes[0, 1].scatter(merged_df['umap_1'], merged_df['umap_2'],\n                               c=merged_df['gc_content'], cmap='coolwarm', alpha=0.6, s=15)\naxes[0, 1].set_title('UMAP - GC Content')\nplt.colorbar(scatter2, ax=axes[0, 1], label='GC%')\n\nscatter3 = axes[1, 0].scatter(merged_df['umap_1'], merged_df['umap_2'],\n                               c=merged_df['radius_of_gyration'], cmap='plasma', alpha=0.6, s=15)\naxes[1, 0].set_title('UMAP - Radius of Gyration')\nplt.colorbar(scatter3, ax=axes[1, 0], label='Rg')\n\nscatter4 = axes[1, 1].scatter(merged_df['umap_1'], merged_df['umap_2'],\n                               c=merged_df['n_ligands'], cmap='YlOrRd', alpha=0.6, s=15)\naxes[1, 1].set_title('UMAP - Number of Ligands')\nplt.colorbar(scatter4, ax=axes[1, 1], label='# Ligands')\n\nplt.suptitle('UMAP Embeddings with Different Feature Colorings', fontsize=14)\nplt.tight_layout()\nplt.savefig('umap_multi_feature.png', dpi=150)\nplt.show()\n\n# ============================================================\n# 14. PCA Feature Importance\n# ============================================================\n\npca_full = PCA(n_components=min(10, len(available_cols)))\npca_full.fit(X_scaled)\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\n# Explained variance\naxes[0].bar(range(1, len(pca_full.explained_variance_ratio_) + 1), \n            pca_full.explained_variance_ratio_ * 100)\naxes[0].set_xlabel('Principal Component')\naxes[0].set_ylabel('Explained Variance (%)')\naxes[0].set_title('PCA Explained Variance')\n\n# Feature loadings for PC1 and PC2\nloadings = pd.DataFrame(\n    pca_full.components_[:2].T,\n    columns=['PC1', 'PC2'],\n    index=available_cols\n)\nloadings_sorted = loadings.reindex(loadings['PC1'].abs().sort_values(ascending=True).index)\nloadings_sorted.plot(kind='barh', ax=axes[1])\naxes[1].set_title('Feature Loadings (PC1 & PC2)')\naxes[1].set_xlabel('Loading')\n\nplt.tight_layout()\nplt.savefig('pca_analysis.png', dpi=150)\nplt.show()\n\n# ============================================================\n# 15. Summary Statistics by Class\n# ============================================================\n\nprint(\"\\n\" + \"=\"*60)\nprint(\"SUMMARY STATISTICS BY CLASS\")\nprint(\"=\"*60)\n\n# By primary ligand\nprint(\"\\nBy Primary Ligand:\")\nprint(merged_df.groupby('primary_ligand')[['seq_length', 'gc_content', 'radius_of_gyration']].mean().round(2))\n\n# By length class\nprint(\"\\nBy Length Class:\")\nprint(merged_df.groupby('length_class')[['gc_content', 'n_ligands', 'radius_of_gyration']].mean().round(2))\n\n# By decade\nprint(\"\\nBy Decade:\")\nprint(merged_df.groupby('decade')[['seq_length', 'n_ligands']].mean().round(2))\n\nprint(\"\\nDone! All plots saved.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-04T16:54:05.328123Z","iopub.execute_input":"2026-02-04T16:54:05.328536Z","iopub.status.idle":"2026-02-04T17:33:47.230201Z","shell.execute_reply.started":"2026-02-04T16:54:05.328508Z","shell.execute_reply":"2026-02-04T17:33:47.228696Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}