{"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":"# Converting Single-Cell Data\nThis notebook will teach you how to convert the files downloaded from Kaggle into AnnData, a common single-cell format that utilises sparse matrices to store the data. We similarly convert the pseudobulk data into AnnData, and multiome data is converted to the multiomics equivalent, called MuData.\n\nWhat you can expect to learn:\n- Memory-sensitive conversion of input files into efficient, compressed, indexed files that are interoperable with great packages\n- Annotation of genes with basic information about gene name, length, and gene type\n\nWhat you cannot expect to learn:\n- Predictive model building & submission\n\nTo learn more about AnnData and MuData, see the below links\n- AnnData: https://anndata.readthedocs.io/\n- MuData: https://mudata.readthedocs.io/\n\nFor packages that can use these files directly:\n- Scanpy: https://scanpy.readthedocs.io/\n- Muon: https://muon.readthedocs.io/\n\nFor example, muon can handle ATAC-seq in the here create MuData object with it's ATAC module: https://muon.readthedocs.io/en/latest/omics/atac.html\n\n## Versions\n- V1: Initial conversions\n- V2: bump dtypes to np.int16 to prevent overflow","metadata":{}},{"cell_type":"code","source":"!pip install anndata pybiomart mudata","metadata":{"execution":{"iopub.status.busy":"2023-10-12T16:47:03.905753Z","iopub.execute_input":"2023-10-12T16:47:03.906745Z","iopub.status.idle":"2023-10-12T16:47:17.492124Z","shell.execute_reply.started":"2023-10-12T16:47:03.906707Z","shell.execute_reply":"2023-10-12T16:47:17.490832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport anndata as ad\nimport mudata as md\nfrom gc import collect\nfrom pybiomart import Dataset\nfrom scipy.sparse import csr_matrix","metadata":{"execution":{"iopub.status.busy":"2023-10-12T16:47:17.494523Z","iopub.execute_input":"2023-10-12T16:47:17.494911Z","iopub.status.idle":"2023-10-12T16:47:22.909387Z","shell.execute_reply.started":"2023-10-12T16:47:17.494878Z","shell.execute_reply":"2023-10-12T16:47:22.908243Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Defining functions to download per-gene information from ENSEMBL","metadata":{}},{"cell_type":"code","source":"def get_biomart_query(symbols):\n    \"\"\"Query biomart for gene metadata given a list of gene symbols\"\"\"\n    dataset = Dataset(name=\"hsapiens_gene_ensembl\", host=\"http://www.ensembl.org\")\n    query = dataset.query(\n        attributes=[\n            \"chromosome_name\",\n            \"start_position\",\n            \"end_position\",\n            \"ensembl_gene_id\",\n            \"external_gene_name\",\n            \"gene_biotype\",\n        ]\n    )\n    query.columns = (\"chr\", \"start\", \"end\", \"gene_id\", \"gene_name\", \"biotype\")\n\n    # Remove duplicates and patches\n    query = (\n        query.loc[~query[\"chr\"].str.contains(\"_\")]\n        .sort_values([\"chr\", \"start\"])\n        .groupby(\"gene_name\")\n        .first()\n        .reset_index()\n        .sort_values([\"chr\", \"start\"])\n    )\n    idx = query[\"gene_name\"].isin(symbols)\n    matched = sum(idx)\n    if len(symbols) > matched:\n        print(f\"Warning: {len(symbols) - matched} symbols not found in biomart\")\n    return query[idx]\n\n\ndef gene_symbol_to_metadata(symbols):\n    \"\"\"Convert a list of gene symbols to a dataframe of gene metadata\"\"\"\n    query = get_biomart_query(symbols)\n    merge = pd.merge(\n        pd.Series(symbols, name=\"gene_name\"),\n        query,\n        how=\"left\",\n    )\n    if len(merge) != len(symbols):\n        raise ValueError(\"Potential duplicates in symbols or biomart query\")\n    return merge","metadata":{"execution":{"iopub.status.busy":"2023-10-12T16:47:22.910634Z","iopub.execute_input":"2023-10-12T16:47:22.911257Z","iopub.status.idle":"2023-10-12T16:47:22.920653Z","shell.execute_reply.started":"2023-10-12T16:47:22.911226Z","shell.execute_reply":"2023-10-12T16:47:22.919697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Converting Single-Cell RNA-seq Data","metadata":{}},{"cell_type":"code","source":"indir = \"/kaggle/input/open-problems-single-cell-perturbations/\"\noutdir = \"/kaggle/working/\"","metadata":{"execution":{"iopub.status.busy":"2023-10-12T16:47:22.923296Z","iopub.execute_input":"2023-10-12T16:47:22.923635Z","iopub.status.idle":"2023-10-12T16:47:22.937304Z","shell.execute_reply.started":"2023-10-12T16:47:22.923606Z","shell.execute_reply":"2023-10-12T16:47:22.935965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read in raw data\ndf = pd.read_parquet(f\"{indir}/adata_train.parquet\", \n                     columns=[\"obs_id\", \"gene\", \"count\"])\n\n# Observation metadata\nobs = pd.read_csv(f\"{indir}/adata_obs_meta.csv\")\n\n# Compress metadata\ndf[\"obs_id\"] = df[\"obs_id\"].astype(\"category\")\ndf[\"gene\"] = df[\"gene\"].astype(\"category\")\n\n# Build row and column metadata\nid_u = list(sorted(df[\"obs_id\"].unique()))\ngene_u = list(sorted(df[\"gene\"].unique()))\nvar = gene_symbol_to_metadata(gene_u).set_index(\"gene_name\")\n\n# Obtain compressed metadata\nrow = df[\"obs_id\"].cat.codes\ncol = df[\"gene\"].cat.codes\n\n# Build sparse, wide matrix\nmat = csr_matrix(\n    (df[\"count\"], (row, col)), shape=(len(id_u), len(gene_u)), dtype=np.int16\n)\n\n# Build final AnnData object\nadata = ad.AnnData(X=mat, obs=obs, var=var)\n\n# Export new data\nadata.write(f\"{outdir}/adata_train.h5ad\")\nprint(adata)\n\ndel adata, mat, row, col, var, obs, df\ncollect()\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T16:47:22.938506Z","iopub.execute_input":"2023-10-12T16:47:22.938811Z","iopub.status.idle":"2023-10-12T16:50:17.137214Z","shell.execute_reply.started":"2023-10-12T16:47:22.938786Z","shell.execute_reply":"2023-10-12T16:50:17.136101Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Converting Pseudobulk Data","metadata":{}},{"cell_type":"code","source":"df = pd.read_parquet(f\"{indir}/de_train.parquet\")\n\nmeta_cols = [\"cell_type\", \"sm_name\", \"sm_lincs_id\", \"SMILES\", \"control\"]\nobs = df[meta_cols]\nX = df.drop(columns=meta_cols).astype(np.float32)\nvar = gene_symbol_to_metadata(X.columns).set_index(\"gene_name\")\n\nadata = ad.AnnData(X=X, obs=obs, var=var)\nadata.write(f\"{outdir}/de_train.h5ad\")\nprint(adata)\ndel df, X, obs, var, adata\ncollect()","metadata":{"execution":{"iopub.status.busy":"2023-10-12T16:50:17.138710Z","iopub.execute_input":"2023-10-12T16:50:17.139008Z","iopub.status.idle":"2023-10-12T16:50:20.311829Z","shell.execute_reply.started":"2023-10-12T16:50:17.138984Z","shell.execute_reply":"2023-10-12T16:50:20.310864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Converting Multiome Data","metadata":{}},{"cell_type":"code","source":"\ndf = pd.read_parquet(\n    f\"{indir}/multiome_train.parquet\", \n    columns=[\"obs_id\", \"location\", \"count\"]\n)\nvar = pd.read_csv(f\"{indir}/multiome_var_meta.csv\")\nobs = pd.read_csv(f\"{indir}/multiome_obs_meta.csv\")\n\n\ndef df_to_adata(df, var_name=\"gene\"):\n    print(\"Converting columns\")\n    meta_cols = df.columns[:2]\n    for m in meta_cols:\n        df[m] = df[m].astype(\"category\")\n    row = df[meta_cols[0]].cat.codes\n    col = df[meta_cols[1]].cat.codes\n    print(\"Building sparse matrix\")\n    mat = csr_matrix(\n        (df[\"count\"], (row, col)),\n        shape=(len(row.unique()), len(col.unique())),\n        dtype=np.int16,\n    )\n    obs = df[meta_cols[0]].cat.categories\n    obs = pd.DataFrame({\"obs_id\": obs})\n    obs.index = obs[\"obs_id\"]\n    var = df[meta_cols[1]].cat.categories\n    var = pd.DataFrame({var_name: var})\n    var.index = var[var_name]\n    adata = ad.AnnData(X=mat, obs=obs, var=var)\n    return adata\n\n\natacidx = df[\"location\"].str.startswith(\"chr\")\nadata_rna = df_to_adata(df[~atacidx].copy())\nadata_atac = df_to_adata(df[atacidx].copy())\n\nadata_rna.var = var.set_index(\"location\").loc[adata_rna.var.index]\nadata_rna.obs = obs.set_index(\"obs_id\").loc[adata_rna.obs.index]\n\nadata_atac.var = var.set_index(\"location\").loc[adata_atac.var.index]\nadata_atac.obs = obs.set_index(\"obs_id\").loc[adata_atac.obs.index]\n\nmdata = md.MuData({\"rna\": adata_rna, \"atac\": adata_atac})\n\nmdata.write(f\"{outdir}/multiome_train.h5mu\")\nprint(mdata)\n# adata_rna is 25551 × 22847\n# adata_atac is 25551 × 135358\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T16:50:20.313603Z","iopub.execute_input":"2023-10-12T16:50:20.314121Z"},"trusted":true},"execution_count":null,"outputs":[]}]}