{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59094,"databundleVersionId":7010844,"sourceType":"competition"},{"sourceId":7279447,"sourceType":"datasetVersion","datasetId":1608748},{"sourceId":7299702,"sourceType":"datasetVersion","datasetId":4234070}],"dockerImageVersionId":30626,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"transcription_factors = [\n    \"ATF1\", \"ATF2\", \"ATF4\", \"ATF6\", \"ATF7\", \"ATF7IP\", \"BTF3\", \"E2F4\", \"ERH\", \"HMGB1\", \n    \"ILF2\", \"IER2\", \"JUND\", \"TCEB2\"\n]\n\nrna_splicing_factors = [\n    \"BAT1\", \"HNRPD\", \"HNRPK\", \"PABPN1\", \"SRSF3\",\n    \"SNRPB\", \"SRSF1\", \"U2AF1\", \"PRPF8\", \"SF3B1\",\n    \"SNRPA\", \"SRSF7\", \"SNRPC\", \"U2AF2\", \"SRSF2\"\n]\n\ntranslation_factors = [\n    \"EIF1\", \"EIF1AD\", \"EIF1B\", \"EIF2A\", \"EIF2AK1\", \"EIF2AK3\", \"EIF2AK4\", \"EIF2B2\", \"EIF2B3\", \n    \"EIF2B4\", \"EIF2S2\", \"EIF3A\", \"EIF3B\", \"EIF3D\", \"EIF3G\", \"EIF3I\", \"EIF3H\", \"EIF3J\", \"EIF3K\", \n    \"EIF3L\", \"EIF3M\", \"EIF3S5\", \"EIF3S8\", \"EIF4A1\", \"EIF4A2\", \"EIF4A3\", \"EIF4E2\", \"EIF4G1\", \n    \"EIF4G2\", \"EIF4G3\", \"EIF4H\", \"EIF5\", \"EIF5A\", \"EIF5AL1\", \"EIF5B\", \"EIF6\", \"TUFM\"\n]\n\ntrna_synthesis_genes = [\n    \"AARS\", \"AARS2\", \"AARSD1\", \"CARS\", \"CARS2\", \"DARS\", \"DARS2\", \"EARS2\", \"FARS2\", \"FARSA\", \n    \"FARSB\", \"GARS\", \"HARS\", \"HARS2\", \"IARS\", \"IARS2\", \"KARS\", \"LARS2\", \"MARS\", \"MARS2\", \n    \"NARS\", \"NARS2\", \"QARS\", \"RARS\", \"RARS2\", \"SARS\", \"TARS\", \"VARS2\", \"WARS2\", \"YARS\", \"YARS2\"\n]\n\nribosomal_proteins = [\n    \"RPL5\", \"RPL8\", \"RPL9\", \"RPL10A\", \"RPL11\", \"RPL14\", \"RPL25\", \"RPL26L1\", \"RPL27\", \"RPL30\", \n    \"RPL32\", \"RPL34\", \"RPL35\", \"RPL35A\", \"RPL36AL\", \"RPS5\", \"RPS6\", \"RPS6KA3\", \"RPS6KB1\", \"RPS6KB2\", \n    \"RPS13\", \"RPS19BP1\", \"RPS20\", \"RPS23\", \"RPS24\", \"RPS27\", \"RPN1\"\n]\n\nmitochondrial_ribosomal_proteins = [\n    \"MRPL9\", \"MRPL1\", \"MRPL10\", \"MRPL11\", \"MRPL12\", \"MRPL13\", \"MRPL14\", \"MRPL15\", \"MRPL16\", \"MRPL17\", \n    \"MRPL18\", \"MRPL19\", \"MRPL2\", \"MRPL20\", \"MRPL21\", \"MRPL22\", \"MRPL23\", \"MRPL24\", \"MRPL27\", \"MRPL28\", \n    \"MRPL3\", \"MRPL30\", \"MRPL32\", \"MRPL33\", \"MRPL35\", \"MRPL36\", \"MRPL37\", \"MRPL38\", \"MRPL4\", \"MRPL40\", \n    \"MRPL41\", \"MRPL42\", \"MRPL43\", \"MRPL44\", \"MRPL45\", \"MRPL46\", \"MRPL47\", \"MRPL48\", \"MRPL49\", \"MRPL50\", \n    \"MRPL51\", \"MRPL52\", \"MRPL53\", \"MRPL54\", \"MRPL55\", \"MRPS10\", \"MRPS11\", \"MRPS12\", \"MRPS14\", \"MRPS15\", \n    \"MRPS16\", \"MRPS17\", \"MRPS18A\", \"MRPS18B\", \"MRPS18C\", \"MRPS2\", \"MRPS21\", \"MRPS22\", \"MRPS23\", \"MRPS24\", \n    \"MRPS25\", \"MRPS26\", \"MRPS27\", \"MRPS28\", \"MRPS30\", \"MRPS31\", \"MRPS33\", \"MRPS34\", \"MRPS35\", \"MRPS5\", \n    \"MRPS6\", \"MRPS7\", \"MRPS9\"\n]\n\nnadh_dehydrogenase_enzymes = [\n    \"NDUFA2\", \"NDUFA3\", \"NDUFA4\", \"NDUFA5\", \"NDUFA6\", \"NDUFA7\", \"NDUFA8\", \"NDUFA9\", \"NDUFA10\", \"NDUFA11\", \n    \"NDUFA12\", \"NDUFA13\", \"NDUFAF2\", \"NDUFAF3\", \"NDUFAF4\", \"NDUFB2\", \"NDUFB3\", \"NDUFB4\", \"NDUFB5\", \"NDUFB6\", \n    \"NDUFB7\", \"NDUFB10\", \"NDUFB11\", \"NDUFB8\", \"NDUFB9\", \"NDUFC1\", \"NDUFC2\", \"NDUFC2-KCTD14\", \"NDUFS2\", \n    \"NDUFS3\", \"NDUFS4\", \"NDUFS5\", \"NDUFS6\", \"NDUFS7\", \"NDUFS8\", \"NDUFV1\", \"NDUFV2\"\n]\n\ncytochrome_c_oxidase_enzymes = [\n    \"COX4I1\", \"COX5B\", \"COX6B1\", \"COX6C\", \"COX7A2\", \"COX7A2L\", \"COX7C\", \"COX8\", \"COX8A\", \"COX11\", \n    \"COX14\", \"COX15\", \"COX16\", \"COX19\", \"COX20\", \"CYC1\", \"UQCC\", \"UQCR10\", \"UQCR11\", \"UQCRB\", \"UQCRC1\", \n    \"UQCRC2\", \"UQCRHL\", \"UQCRQ\"\n]\n\ngene_groups = {\n    'transcription_factors': {'genes': transcription_factors, 't': 0.5},\n    'rna_splicing_factors': {'genes': rna_splicing_factors, 't': 0.5},\n    'translation_factors': {'genes': translation_factors, 't': 0.34},\n    'trna_synthesis_genes': {'genes': trna_synthesis_genes, 't': 0.51},\n    'ribosomal_proteins': {'genes': ribosomal_proteins, 't': 0.4},\n    'mitochondrial_ribosomal_proteins': {'genes': mitochondrial_ribosomal_proteins, 't': 0.4},\n    'nadh_dehydrogenase_enzymes': {'genes': nadh_dehydrogenase_enzymes, 't': 0.4},\n    'cytochrome_c_oxidase_enzymes': {'genes': cytochrome_c_oxidase_enzymes, 't': 0.4}\n}","metadata":{"execution":{"iopub.status.busy":"2023-12-30T19:39:43.248942Z","iopub.execute_input":"2023-12-30T19:39:43.249382Z","iopub.status.idle":"2023-12-30T19:39:43.266503Z","shell.execute_reply.started":"2023-12-30T19:39:43.249351Z","shell.execute_reply":"2023-12-30T19:39:43.265435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install anndata scanpy > /dev/null 2>&1","metadata":{"execution":{"iopub.status.busy":"2023-12-30T19:39:43.378525Z","iopub.execute_input":"2023-12-30T19:39:43.379418Z","iopub.status.idle":"2023-12-30T19:39:55.395710Z","shell.execute_reply.started":"2023-12-30T19:39:43.379387Z","shell.execute_reply":"2023-12-30T19:39:55.394418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import anndata\nimport scanpy as sc\nimport pandas as pd\nimport seaborn as sns\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.cluster import hierarchy\nfrom scipy.spatial import distance","metadata":{"execution":{"iopub.status.busy":"2023-12-30T19:39:55.398563Z","iopub.execute_input":"2023-12-30T19:39:55.399436Z","iopub.status.idle":"2023-12-30T19:39:55.404876Z","shell.execute_reply.started":"2023-12-30T19:39:55.399396Z","shell.execute_reply":"2023-12-30T19:39:55.403799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_de_train = pd.read_parquet('/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet')\n\ndf_hk_genes = pd.read_csv('/kaggle/input/genes-information/Human_housekeeping_genes_Eisenberg_Levanon_2013_TrendsGenetics.csv',index_col = 0)\nl = df_hk_genes.iloc[:,0]\ns = set(l) & set(df_de_train.columns)\nlist_not_hk = list( set(df_de_train.columns[5:]) - set(l) )\nlist_hk = list(s)\n\nadata = anndata.read_h5ad('/kaggle/input/open-problems-single-cell-perturbations-data/adata_OPSCP_FilteredNormalized_sparse_240074_21231.h5ad')","metadata":{"execution":{"iopub.status.busy":"2023-12-30T19:39:55.406183Z","iopub.execute_input":"2023-12-30T19:39:55.406473Z","iopub.status.idle":"2023-12-30T19:40:19.321368Z","shell.execute_reply.started":"2023-12-30T19:39:55.406449Z","shell.execute_reply":"2023-12-30T19:40:19.320562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_gene_correlation_clustermap(df, gene_group, t):\n    \n    pd.set_option('display.max_rows', None)\n    pd.set_option('display.max_columns', None)\n    pd.set_option('display.max_colwidth', None)\n    \n    available_genes = [gene for gene in gene_group if gene in df.columns]\n    df_filtered = df[available_genes]\n    correlation_matrix = df_filtered.corr()\n    clustergrid = sns.clustermap(correlation_matrix, cmap='coolwarm', center=0, robust=True, figsize=(5,5))\n    plt.show()\n\n    reordered_genes = [available_genes[i] for i in clustergrid.dendrogram_row.reordered_ind]\n    linkage_matrix = hierarchy.linkage(distance.pdist(correlation_matrix), 'average')\n    threshold = t * max(distance.pdist(correlation_matrix))\n    clusters = hierarchy.fcluster(linkage_matrix, threshold, criterion='distance')\n    \n    cluster_data = pd.DataFrame({'Gene': reordered_genes, 'Cluster': clusters})\n\n    grouped = cluster_data.groupby('Cluster')['Gene'].apply(list).reset_index()\n    return grouped\n\ndef plot_gene_expression(genes_list, group_name, cluster_number):\n    list_thresholds = [0.1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9, 1e-10, 1e-11, 1e-12, 1e-13, 1e-14, 1e-15, 1e-16]\n\n    def filter_existing_genes(genes):\n        return [gene for gene in genes if gene in df_de_train.columns]\n\n    filtered_genes = filter_existing_genes(genes_list)\n    filtered_not_hk = filter_existing_genes(list_not_hk)\n    filtered_hk = filter_existing_genes(list_hk)\n\n    d = pd.DataFrame()\n    IX = -1\n\n    group_label = f\"{group_name} Cluster {cluster_number}\"\n    for gene_group, gene_list in {group_label: filtered_genes, \"Not Housekeeping\": filtered_not_hk, \"Housekeeping\": filtered_hk}.items():\n        if not gene_list:\n            continue\n        for t in list_thresholds:\n            v = df_de_train[gene_list].values.ravel()\n            v = 10**(-np.abs(v))\n            m = v < t\n            IX += 1\n            d.loc[IX, 'Threshold'] = t\n            d.loc[IX, f'{gene_group} percent'] = 100 * m.sum() / len(v)\n\n    plt.figure(figsize=(15, 3))\n    for gene_group in [group_label, \"Not Housekeeping\", \"Housekeeping\"]:\n        plt.plot(d['Threshold'], d[f'{gene_group} percent'], label=gene_group, linewidth=1)\n\n    plt.grid()\n    plt.legend()\n    plt.xscale('log')\n    plt.xlabel('Threshold')\n    plt.ylabel('Percent DE stronger than threshold')\n    plt.title('Percent DE Expressed Stronger Than Threshold')\n    x_min, x_max = plt.xlim()\n    plt.xlim(x_max, x_min)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-30T19:40:19.322881Z","iopub.execute_input":"2023-12-30T19:40:19.323234Z","iopub.status.idle":"2023-12-30T19:40:19.341123Z","shell.execute_reply.started":"2023-12-30T19:40:19.323201Z","shell.execute_reply":"2023-12-30T19:40:19.340124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"gene_groups = {\n    'transcription_factors': {'genes': transcription_factors, 't': 0.5}}\"\"\"","metadata":{"execution":{"iopub.status.busy":"2023-12-30T19:40:19.343452Z","iopub.execute_input":"2023-12-30T19:40:19.344148Z","iopub.status.idle":"2023-12-30T19:40:19.355900Z","shell.execute_reply.started":"2023-12-30T19:40:19.344120Z","shell.execute_reply":"2023-12-30T19:40:19.355068Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from IPython.display import display\n    \nfor group_name, info in gene_groups.items():\n    genes = info['genes']\n    t_value = info['t']\n    \n    print()\n    print(f\"Group Name: {group_name}\")\n    \n    df = create_gene_correlation_clustermap(df_de_train, genes, t_value)\n    display(df)\n    \n    for cluster in df['Cluster'].unique():\n        cluster_genes = df[df['Cluster'] == cluster]['Gene'].iloc[0]\n        valid_genes = list(set(cluster_genes) & set(adata.var.index))\n\n        if valid_genes:\n            cluster_mean = np.asarray(adata[:, valid_genes].X.mean(axis=1)).ravel()\n            adata.obs[f'Cluster_{group_name}_{cluster}_mean'] = cluster_mean\n\n            plot_gene_expression(valid_genes, group_name, cluster)\n            \n        mean_columns = [col for col in adata.obs.columns if col.startswith(f'Cluster_{group_name}_')]\n\n    sc.pl.umap(adata, color=mean_columns)\n","metadata":{"execution":{"iopub.status.busy":"2023-12-30T19:40:19.357143Z","iopub.execute_input":"2023-12-30T19:40:19.357492Z"},"trusted":true},"execution_count":null,"outputs":[]}]}