{
"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
}