{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"This is a short introduction of [Mapper Method](https://research.math.osu.edu/tgda/mapperPBG.pdf). \n\n**What is Mapper?**\n\nMapper is a method that can be use to reduce high dimensional datasets into multiple sets of compoused of points which that capture topological and geometric information at a specified resolution. This method outputs a coordinatiozation that uses a discrete and combinatorial object to which the dataset maps and which can represent the dataset in a useful way. Futhurmore, mapper can represent higher dimensional objects, such as spheres, tori, etc.","metadata":{}},{"cell_type":"markdown","source":"The basic idea of Mapper can be referred to as partial clustering. If $U$ and $V$ are subsets of the dataset and $U\\cap V\\neq\\phi$, the clusters obtained from $U$ and $V$ have non-empty intersections, and these intersections can be used in constructing a simplicial complex. ","metadata":{}},{"cell_type":"markdown","source":"#### Example\n\nLet $X$ debe the unit circle $\\{(x,y)\\,:\\,x^2+y^2=1\\}$. \n\nTo apply Mapper on $X$, we first transforms $X$ to a simpler set. For the shake of simplicity, let $Z=[-1,1]$ and define $f:X\\to Z$ such that $f(x,y)=y$. Additionally, $Z$ needs to be equipped with a convering $\\mathcal{U}=\\{U_{\\alpha}\\}_{\\alpha\\in A}$ where $A$ is a finite indexing set. Mathematically, we need $Z\\subseteq\\mathcal{U}$ and $\\cap_{\\alpha\\in A}U_{\\alpha}\\neq\\phi$. A simple convering for $Z$ can be $\\{[-1,-2/3),[-1/2, 1/2],(2/3,1]\\}$. For clustering, we can simply use the inverse function $f^{-1}$. We have $f^{-1}([-1,-2/3))=[0, \\sqrt{5}/3)\\times[-1,-2/3),f^{-1}([-1/2,1/2])=[-\\sqrt{3}/2, -1]\\times[-1/2,1/2]\\cup[\\sqrt{3}/2, 1]\\times[-1/2,1/2]$ and $f^{-1}((2/3,1])=[0, \\sqrt{5}/3)\\times[-1,-2/3)$\n\nThe simplicial complex can be viewed as a dimond where the top node is $f^{-1}((2/3,1])$, two nodes on the sides are the two disjoint sets of $f^{-1}([-1/2,1/2])$, and the bottom node is $f^{-1}([-1,-2/3))$.","metadata":{}},{"cell_type":"code","source":"import numpy as np, pandas as pd, seaborn as sns, matplotlib.pyplot as plt\nimport warnings, time, gc, os\nfrom sklearn.feature_extraction.text import TfidfVectorizer, CountVectorizer\n\nfrom sklearn.decomposition import TruncatedSVD, LatentDirichletAllocation\nfrom sklearn.manifold import TSNE, Isomap\nfrom sklearn.decomposition import PCA\n\nimport kmapper as km\nfrom kmapper import Cover, jupyter\n\nfrom sklearn.cluster import AgglomerativeClustering\nfrom sklearn.feature_extraction.text import TfidfVectorizer\nfrom sklearn.decomposition import TruncatedSVD, LatentDirichletAllocation\nfrom sklearn import cluster\nfrom sklearn.preprocessing import MinMaxScaler, StandardScaler\nimport sklearn.covariance, scipy.stats, math\n\nimport igraph\n\nimport bokeh.plotting as bp\nfrom bokeh.models import HoverTool, BoxSelectTool\nfrom bokeh.models import ColumnDataSource\nfrom bokeh.plotting import figure, show, output_notebook, reset_output\nfrom bokeh.palettes import d3\nimport bokeh.models as bmo\nfrom bokeh.io import save, output_file\n\nfrom bokeh.io import output_notebook\nfrom bokeh.resources import INLINE\noutput_notebook(INLINE)\n\nnp.random.seed(1024)","metadata":{"execution":{"iopub.status.busy":"2023-11-11T04:27:32.422086Z","iopub.execute_input":"2023-11-11T04:27:32.422509Z","iopub.status.idle":"2023-11-11T04:27:36.532568Z","shell.execute_reply.started":"2023-11-11T04:27:32.422475Z","shell.execute_reply":"2023-11-11T04:27:36.530648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_DIR = \"/kaggle/input/stanford-ribonanza-rna-folding\"\n\n# DEBUG to only open a fraction\nDEBUG = True\n\n# Define the paths to the batches of train and test data respectively\nTRAIN_CSV_PATH = os.path.join(DATA_DIR, \"train_data.csv\")\n\ntrain_df = pd.read_csv(TRAIN_CSV_PATH, nrows=10000 if DEBUG else None)\ntrain_df = train_df[['sequence_id', 'sequence', 'experiment_type', 'dataset_name', 'reads','signal_to_noise', 'SN_filter']]\ntrain_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-11-11T04:27:36.536257Z","iopub.execute_input":"2023-11-11T04:27:36.537118Z","iopub.status.idle":"2023-11-11T04:27:37.398745Z","shell.execute_reply.started":"2023-11-11T04:27:36.537045Z","shell.execute_reply":"2023-11-11T04:27:37.397411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df[\"seq_len\"] = train_df[\"sequence\"].apply(lambda x:len(x))","metadata":{"execution":{"iopub.status.busy":"2023-11-11T04:30:31.677859Z","iopub.execute_input":"2023-11-11T04:30:31.678387Z","iopub.status.idle":"2023-11-11T04:30:32.243585Z","shell.execute_reply.started":"2023-11-11T04:30:31.678349Z","shell.execute_reply":"2023-11-11T04:30:32.242181Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub =  train_df.sample(n = 5000, random_state = 2048)","metadata":{"execution":{"iopub.status.busy":"2023-11-11T04:27:37.420643Z","iopub.execute_input":"2023-11-11T04:27:37.421242Z","iopub.status.idle":"2023-11-11T04:27:37.432864Z","shell.execute_reply.started":"2023-11-11T04:27:37.421182Z","shell.execute_reply":"2023-11-11T04:27:37.430877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del train_df; gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-11-11T04:27:37.436807Z","iopub.execute_input":"2023-11-11T04:27:37.437842Z","iopub.status.idle":"2023-11-11T04:27:37.563757Z","shell.execute_reply.started":"2023-11-11T04:27:37.437805Z","shell.execute_reply":"2023-11-11T04:27:37.562250Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mapper = km.KeplerMapper(verbose = 1)\n\nprojected_X = mapper.fit_transform(np.array(sub.sequence),\n                                   projection = [CountVectorizer(analyzer = \"char\",\n                                                                 ngram_range = (3,6),\n                                                                 min_df = 5,\n                                                                 max_features = 700),\n                                               TruncatedSVD(n_components = 50,\n                                                            random_state = 2048),\n                                               TSNE(n_components = 2,\n                                                      n_jobs = -1)],\n                                   scaler = [None, None, StandardScaler()])\n\n\ngraph = mapper.map(projected_X,\n                   X = None,\n                   clusterer = cluster.AgglomerativeClustering(n_clusters = 5,\n                                                             linkage = \"complete\",\n                                                             metric = \"cosine\"),\n                   cover = Cover(n_cubes = 10, perc_overlap = 0.33))","metadata":{"execution":{"iopub.status.busy":"2023-11-11T04:35:01.489471Z","iopub.execute_input":"2023-11-11T04:35:01.490142Z","iopub.status.idle":"2023-11-11T04:35:33.056851Z","shell.execute_reply.started":"2023-11-11T04:35:01.490096Z","shell.execute_reply":"2023-11-11T04:35:33.055587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_ = mapper.visualize(graph,\n                     path_html = \"./RNA_cluster.html\",\n                     lens_names = [\"TSNE_X\", \"TSNE_Y\"],\n                     title = \"RNA Sequence Analysis\")\n\njupyter.display(\"./RNA_cluster.html\")","metadata":{"execution":{"iopub.status.busy":"2023-11-11T04:39:28.366033Z","iopub.execute_input":"2023-11-11T04:39:28.366727Z","iopub.status.idle":"2023-11-11T04:39:29.698329Z","shell.execute_reply.started":"2023-11-11T04:39:28.366681Z","shell.execute_reply":"2023-11-11T04:39:29.697113Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub.iloc[graph[\"nodes\"][\"cube2_cluster0\"]]","metadata":{"execution":{"iopub.status.busy":"2023-11-11T04:54:41.151150Z","iopub.execute_input":"2023-11-11T04:54:41.151599Z","iopub.status.idle":"2023-11-11T04:54:41.179266Z","shell.execute_reply.started":"2023-11-11T04:54:41.151565Z","shell.execute_reply":"2023-11-11T04:54:41.178286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub.iloc[graph[\"nodes\"][\"cube1_cluster2\"]]","metadata":{"execution":{"iopub.status.busy":"2023-11-11T04:54:20.485103Z","iopub.execute_input":"2023-11-11T04:54:20.485552Z","iopub.status.idle":"2023-11-11T04:54:20.511304Z","shell.execute_reply.started":"2023-11-11T04:54:20.485518Z","shell.execute_reply":"2023-11-11T04:54:20.509709Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nodes = list(graph[\"nodes\"])\nnode_size = [0]*len(nodes)\n\nfor i, node in enumerate(nodes):\n    node_size[i] = len(graph[\"nodes\"][node])\n    \nplt.figure()\nsns.histplot(node_size, kde = True)\nplt.xlabel(\"Node Size\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-11T04:50:48.244804Z","iopub.execute_input":"2023-11-11T04:50:48.245291Z","iopub.status.idle":"2023-11-11T04:50:48.605838Z","shell.execute_reply.started":"2023-11-11T04:50:48.245257Z","shell.execute_reply":"2023-11-11T04:50:48.604632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}