{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import warnings\n", "warnings.simplefilter(action='ignore', category=FutureWarning)\n", "\n", "import scanpy as sc\n", "import anndata as ad\n", "import numpy as np\n", "import pandas as pd\n", "import anndata2ri\n", "import scipy\n", "import scipy.io\n", "from rpy2.robjects import r\n", "\n", "anndata2ri.activate()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%load_ext rpy2.ipython" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "%%R\n", "suppressPackageStartupMessages({\n", " library(SingleCellExperiment)\n", " library(Seurat)\n", " library(SeuratDisk)\n", "})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Load data" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 10849 × 15792\n", " obs: 'n_genes', 'percent_mito', 'n_counts', 'batch', 'labels'\n", " var: 'encode', 'n_cells-0', 'highly_variable-0', 'hvg_encode-0', 'n_cells-1', 'highly_variable-1', 'hvg_encode-1'\n", " obsm: 'protein_expression'" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata = sc.read('../data/10x-cite/10x-cite.h5ad')\n", "adata" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Seurat categories have to be strings\n", "adata.obs.batch = adata.obs.batch.astype('category')\n", "adata.obs.batch = adata.obs.batch.cat.rename_categories([\"0\", \"1\"])" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "View of AnnData object with n_obs × n_vars = 3994 × 15792\n", " obs: 'n_genes', 'percent_mito', 'n_counts', 'batch', 'labels'\n", " var: 'encode', 'n_cells-0', 'highly_variable-0', 'hvg_encode-0', 'n_cells-1', 'highly_variable-1', 'hvg_encode-1'\n", " obsm: 'protein_expression'" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "query = adata[adata.obs.batch == '1']\n", "query" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "View of AnnData object with n_obs × n_vars = 6855 × 15792\n", " obs: 'n_genes', 'percent_mito', 'n_counts', 'batch', 'labels'\n", " var: 'encode', 'n_cells-0', 'highly_variable-0', 'hvg_encode-0', 'n_cells-1', 'highly_variable-1', 'hvg_encode-1'\n", " obsm: 'protein_expression'" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata = adata[adata.obs.batch == '0']\n", "adata" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "adata_ = ad.AnnData(adata.X)\n", "adata_.obs_names = adata.obs_names\n", "adata_.var_names = adata.var_names\n", "adata_.obs['cell_type'] = adata.obs['labels']" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "An object of class Seurat \n", "15792 features across 6855 samples within 1 assay \n", "Active assay: RNA (15792 features, 0 variable features)\n" ] } ], "source": [ "%%R -i adata_\n", "rna = as.Seurat(adata_, counts='X', data=NULL)\n", "rna" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "adata_ = ad.AnnData(adata.obsm['protein_expression'])\n", "adata_.obs['cell_type'] = adata.obs['labels']" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Warning:\n", "R[write to console]: Feature names cannot have underscores ('_'), replacing with dashes ('-')\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "An object of class Seurat \n", "14 features across 6855 samples within 1 assay \n", "Active assay: RNA (14 features, 0 variable features)\n" ] } ], "source": [ "%%R -i adata_\n", "cite = as.Seurat(adata_, counts='X', data=NULL)\n", "cite" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] \"Preprocessing RNA...\"\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Warning:\n", "R[write to console]: Not all features provided are in this Assay object, removing the following feature(s): AL034417.3, FHAD1, AL121992.2, SLC30A2, AL512408.1, SFN, CFAP57, AL139220.2, RAB3B, TMEM56, AL603832.1, AL390066.1, SPAG17, PDZK1, NUDT4B, AL450992.1, AL139011.2, AL356441.1, GPR52, AL596220.1, C4BPB, TGFB2, CAPN8, MIR3681HG, TCF7L1, NT5DC4, CKAP2L, LINC01412, SATB2, AC007879.3, MAP2, ABCB6, CHL1, ARPC4-TTLL3, AC011816.2, AC104184.1, AC134772.1, FAM3D, C3ORF67, SCAANT1, ZBTB20-AS4, SLC12A8, LINC02029, AC078802.1, AC008040.1, AC131235.3, AC092546.1, AC114811.2, AC097460.1, AC109347.2, GYPE, WWC2-AS2, AC084871.2, CCDC110, AC092354.1, MAST4-AS1, AC010273.1, ADGRV1, AC008522.1, GPRIN1, AL138831.2, HIST1H2BI, AL358933.1, TRIM31, KCNK5, LRFN2, GUCA1B, AL136304.1, RSPH9, AL133406.3, ENPP3, TAB2-AS1, AL078581.1, ULBP1, SOD2-OT1, AC004948.1, RAPGEF5, ADCY1, CLDN12, AC002456.1, AC006329.1, AC087071.1, AC008264.2, TMEM139, GHET1, AC093458.1, INE2, FAM156A, TRO, OTUD6A, CENPI, ALG13-AS1, KLHL13, AFF2, L1CAM, AC2 [... truncated]\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[1] \"Preprocessing ADT...\"\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Warning in irlba(A = t(x = object), nv = npcs, ...) :\n", "R[write to console]: \n", " \n", "R[write to console]: You're computing too large a percentage of total singular values, use a standard svd instead.\n", "\n" ] } ], "source": [ "%%R\n", "bm <- rna\n", "bm[[\"ADT\"]] <- CreateAssayObject(counts = cite@assays$RNA@counts)\n", "\n", "rm(rna)\n", "rm(cite)\n", "\n", "print('Preprocessing RNA...')\n", "DefaultAssay(bm) <- \"RNA\"\n", "bm <- SCTransform(bm, verbose = FALSE)\n", "VariableFeatures(bm) <- rownames(bm[[\"RNA\"]])\n", "bm <- ScaleData(bm, verbose=FALSE)\n", "bm <- RunPCA(bm, verbose=FALSE)\n", "\n", "print('Preprocessing ADT...')\n", "DefaultAssay(bm) <- \"ADT\"\n", "VariableFeatures(bm) <- rownames(bm[[\"ADT\"]])\n", "bm <- NormalizeData(bm, normalization.method = 'CLR', margin = 2, verbose=FALSE)\n", "bm <- ScaleData(bm, verbose=FALSE)\n", "bm <- RunPCA(bm, reduction.name = \"apca\", verbose=FALSE)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Calculating cell-specific modality weights\n", "\n", "R[write to console]: Finding 20 nearest neighbors for each modality.\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " | | 0 % ~calculating |+++++++++++++++++++++++++ | 50% ~04s |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=07s \n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Calculating kernel bandwidths\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " | | 0 % ~calculating |+++++++++++++++++++++++++ | 50% ~00s |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s \n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Finding multimodal neighbors\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " | | 0 % ~calculating |+++++++++++++++++++++++++ | 50% ~11s |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=23s \n", " | | 0 % ~calculating |+++++++++++++++++++++++++ | 50% ~01s |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s \n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Constructing multimodal KNN graph\n", "\n", "R[write to console]: Constructing multimodal SNN graph\n", "\n" ] } ], "source": [ "%%R\n", "bm <- FindMultiModalNeighbors(\n", " bm, reduction.list = list(\"pca\", \"apca\"), \n", " dims.list = list(1:30, 1:13), modality.weight.name = \"RNA.weight\"\n", " )" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "adata_ = ad.AnnData(query.X)\n", "adata_.obs_names = query.obs_names\n", "adata_.var_names = query.var_names\n", "adata_.obs['cell_type'] = query.obs['labels']" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "An object of class Seurat \n", "15792 features across 3994 samples within 1 assay \n", "Active assay: RNA (15792 features, 0 variable features)\n" ] } ], "source": [ "%%R -i adata_\n", "query = as.Seurat(adata_, counts='X', data=NULL)\n", "query" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "%%R\n", "query <- SCTransform(query, verbose = FALSE)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "An object of class Seurat \n", "31380 features across 6855 samples within 3 assays \n", "Active assay: RNA (15792 features, 15792 variable features)\n", " 2 other assays present: ADT, SCT\n", " 2 dimensional reductions calculated: pca, apca\n" ] } ], "source": [ "%%R\n", "DefaultAssay(bm) <- \"RNA\"\n", "VariableFeatures(bm) <- rownames(bm[[\"RNA\"]])\n", "bm" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Computing sPCA transformation\n", "\n" ] } ], "source": [ "%%R\n", "bm <- RunSPCA(bm, assay = 'SCT', graph = 'wsnn')" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "An object of class Seurat \n", "31380 features across 6855 samples within 3 assays \n", "Active assay: RNA (15792 features, 15792 variable features)\n", " 2 other assays present: ADT, SCT\n", " 3 dimensional reductions calculated: pca, apca, spca\n" ] } ], "source": [ "%%R\n", "bm" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SCTAssay data with 15574 features for 6855 cells, and 1 SCTModel(s) \n", "Top 10 variable features:\n", " AL627309.1, AL669831.5, LINC00115, FAM41C, NOC2L, KLHL17, PLEKHN1,\n", "AL645608.8, HES4, ISG15 \n" ] } ], "source": [ "%%R\n", "DefaultAssay(bm) <- \"RNA\"\n", "bm[['SCT']]" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Normalizing query using reference SCT model\n", "\n", "R[write to console]: Projecting cell embeddings\n", "\n", "R[write to console]: Finding neighborhoods\n", "\n", "R[write to console]: Finding anchors\n", "\n", "R[write to console]: \tFound 2203 anchors\n", "\n" ] } ], "source": [ "%%R\n", "anchors <- FindTransferAnchors(\n", " reference = bm,\n", " reference.assay = 'SCT',\n", " query = query,\n", " normalization.method = \"SCT\",\n", " reference.reduction = \"spca\"\n", ")" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: Finding integration vectors\n", "\n", "R[write to console]: Finding integration vector weights\n", "\n", "R[write to console]: 0% 10 20 30 40 50 60 70 80 90 100%\n", "\n", "R[write to console]: [----|----|----|----|----|----|----|----|----|----|\n", "\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: *\n", "R[write to console]: |\n", "\n", "R[write to console]: Transfering 14 features onto reference data\n", "\n", "R[write to console]: Warning:\n", "R[write to console]: Keys should be one or more alphanumeric characters followed by an underscore, setting key from predicted_ADT_ to predictedADT_\n", "\n" ] } ], "source": [ "%%R\n", "query <- TransferData(\n", " anchorset = anchors, \n", " reference = bm,\n", " query = query,\n", " refdata = list(\n", " predicted_ADT = \"ADT\")\n", ")" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "An object of class Seurat \n", "31329 features across 3994 samples within 3 assays \n", "Active assay: predicted_ADT (14 features, 0 variable features)\n", " 2 other assays present: RNA, SCT\n" ] } ], "source": [ "%%R\n", "DefaultAssay(query) <- \"predicted_ADT\"\n", "query" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "%%R -o query_sc\n", "query_sc <- as.SingleCellExperiment(query)\n", "query_sc <- as(query_sc, 'SingleCellExperiment')" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "query_sc.obsm['predicted_protein_expression'] = query_sc.X.A" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['CD3: Corr=0.937',\n", " 'CD4: Corr=0.9',\n", " 'CD8a: Corr=0.737',\n", " 'CD14: Corr=0.956',\n", " 'CD15: Corr=0.216',\n", " 'CD16: Corr=0.813',\n", " 'CD56: Corr=0.766',\n", " 'CD19: Corr=0.888',\n", " 'CD25: Corr=0.555',\n", " 'CD45RA: Corr=0.836',\n", " 'CD45RO: Corr=0.792',\n", " 'PD-1: Corr=0.455',\n", " 'TIGIT: Corr=0.571',\n", " 'CD127: Corr=0.849']" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy.stats import pearsonr\n", "\n", "protein_corrs = []\n", "for i, protein in enumerate(query_sc.var_names):\n", " protein = protein[:-10]\n", " protein_corrs.append(protein + ': Corr=' + str(np.round(pearsonr(np.log1p(query_sc.obsm['predicted_protein_expression'][:, i]), np.log1p(query.obsm['protein_expression'].iloc[:, i].values))[0], 3)))\n", " \n", "protein_corrs" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" } }, "nbformat": 4, "nbformat_minor": 5 }