testData/BD_results/CombinedAnalysis_scanpy_v1.2.1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a23ff812-9ecf-4d88-99f3-b335e9d18bf2",
   "metadata": {},
   "source": [
    "# Analyze the BD results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fecabdf4-3f41-4e67-80e9-0105e492e8ef",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import scvelo as scv\n",
    "import scanpy\n",
    "import igraph\n",
    "import glob, os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import re\n",
    "from collections import Counter\n",
    "import anndata\n",
    "\n",
    "import shutil\n",
    "    \n",
    "import h5py\n",
    "from shutil import copyfile\n",
    "\n",
    "\n",
    "def copyFiles(files, to):\n",
    "    for f in files:\n",
    "        name = os.path.basename( f )\n",
    "        print( f\"copy {f} to {to}\" )\n",
    "        copyfile( f, os.path.join(to, name ) )\n",
    "    print( \"all copied\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af033cc9-5df3-4380-998c-557991ffccdf",
   "metadata": {},
   "source": [
    "## This now (2023.05.20)\n",
    "\n",
    "is based on the Rhapsody 1.2.1 (first commit).\n",
    "The data was created in the TestData folder using the commands\n",
    "\n",
    "```\n",
    "../target/release/quantify_rhapsody -r cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz -f cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz -o BD_results/Rustody_S1 -s mouse  -e 2276_20220531_chang_to_rpl36a_amplicons.fasta -a MyAbSeqP\n",
    "anel.fasta -m 200 -v v2.96 --gene-kmers 16\n",
    "\n",
    "../target/release/quantify_rhapsody -r cells.1.Rhapsody_SV_index2_S2_R1_001.fastq.gz -f cells.1.Rhapsody_SV_index2_S2_R2_001.fastq.gz -o BD_results/Rustody_S2 -s mouse  -e 2276_20220531_chang_to_rpl36a_amplicons.fasta -a MyAbSeqPanel.fasta -m 200 -v v2.96 --gene-kmers 16\n",
    "```\n",
    "\n",
    "The --gene_kmers makes a huge difference - still.\n",
    "This is the setting where I hat more antibody tags out.\n",
    "\n",
    "I am shure this mapper is more stringent than my first implementation as this here actually checks two consecutive DNA fragments 8bp +16bp. Hence the matching are is bigger than in version 0.x. \n",
    "\n",
    "## Overall result\n",
    "\n",
    "Looks acceptable. I identify less reads as containing sequence, but the result still looks more or less the same.\n",
    "I need to look into that in more detail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c7f984e2-12a2-40c4-9142-cf4e0bc462b6",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scvelo==0.2.4\n",
      "scanpy==1.8.2\n",
      "igraph==0.9.9\n",
      "pandas==1.5.0.dev0+268.gbe8d1ec88\n",
      "numpy==1.21.6\n",
      "re==2.2.1\n",
      "anndata==0.8.0\n",
      "h5py==3.7.0\n"
     ]
    }
   ],
   "source": [
    "print('\\n'.join(f'{m.__name__}=={m.__version__}' for m in globals().values() if getattr(m, '__version__', None)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0afb3d34-3683-4621-b44f-e8c26e3222eb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "files = [ 'Combined_S1Subset_DBEC_MolsPerCell.csv', 'Combined_S2Subset_DBEC_MolsPerCell.csv' ]\n",
    "sinfoF = [ 'S1Subset_Sample_Tag_Calls.csv', 'S2Subset_Sample_Tag_Calls.csv']\n",
    "sname = ['Index1', 'Index2' ]\n",
    "\n",
    "ofiile = \"BD_analyzed.h5ad\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a6119570-5b3e-4dac-aca9-ede6be6813c6",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['CD117:2B8|Kit|AMM2023|pAbO', 'CD19|Cd19|AMM2007|pAbO',\n",
       "       'CD25:PC61|Il2ra|AMM2012|pAbO', 'CD5|Cd5|AMM2043|pAbO',\n",
       "       'IgM|Ighm|AMM2031|pAbO', '2810417H13Rik', 'Ada', 'Adgre1', 'Adgrg3',\n",
       "       'Aicda', 'Alas2', 'Anxa5', 'Apoe', 'Aqp9', 'Arg1', 'Arg2', 'Arid3a',\n",
       "       'Arl4c', 'Atf6b', 'Atg5'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(files[0], skiprows=7, index_col=0)\n",
    "df.columns[0:20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "888aee57-b0dc-4f5a-94ea-257e1c1a96e1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "m = np.array([ s.__contains__ (\"|\") for s in df.columns ]).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bf07d0df-3b6a-467b-86d8-d1bc2d00a4a1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata = anndata.AnnData(X = df.iloc[0:,m:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f724c794-5c1d-48fe-898c-918dc3d4f969",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata.var_names\n",
    "adata.obs['sample'] = sname[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "90f3cc42-36bb-431e-a242-712978acd5f4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(files[1], skiprows=7, index_col=0)\n",
    "tmp = anndata.AnnData(X = df.iloc[0:,m:])\n",
    "tmp.obs['sample'] = sname[1]\n",
    "adata = adata.concatenate( tmp )\n",
    "BD_Data = adata\n",
    "del adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1fe3e876-0295-4eb3-b72a-dedc04041b92",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'posix'"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bc0516ff-3188-4e68-8951-363d9967f50c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[1m\u001b[33mwarning\u001b[0m\u001b[0m\u001b[1m: function `check_file_existence` is never used\u001b[0m\n",
      "\u001b[0m  \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m--> \u001b[0m\u001b[0msrc/bin/quantify_rhapsody_multi.rs:82:4\u001b[0m\n",
      "\u001b[0m   \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\n",
      "\u001b[0m\u001b[1m\u001b[38;5;12m82\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0mfn check_file_existence(file_path: &str, option: &str, errors: &mut Vec<String>) {\u001b[0m\n",
      "\u001b[0m   \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m| \u001b[0m\u001b[0m   \u001b[0m\u001b[0m\u001b[1m\u001b[33m^^^^^^^^^^^^^^^^^^^^\u001b[0m\n",
      "\u001b[0m   \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\n",
      "\u001b[0m   \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m= \u001b[0m\u001b[0m\u001b[1mnote\u001b[0m\u001b[0m: `#[warn(dead_code)]` on by default\u001b[0m\n",
      "\n",
      "\u001b[0m\u001b[0m\u001b[1m\u001b[33mwarning\u001b[0m\u001b[1m:\u001b[0m `rustody` (bin \"quantify_rhapsody_multi\") generated 1 warning\n",
      "\u001b[0m\u001b[0m\u001b[1m\u001b[32m    Finished\u001b[0m release [optimized] target(s) in 0.04s\n"
     ]
    }
   ],
   "source": [
    "! source ~/.cargo/env &&  cargo build -r  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0331fb56-8b18-4e58-8223-00bef94a11aa",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "../../target/release/quantify_rhapsody_multi -r ../cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz -f ../cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz -o Rustody_S1 -s mouse  -e ../2276_20220531_chang_to_rpl36a_amplicons.fasta -a ../MyAbSeqPanel.fasta -m 200 -v v2.96\n",
      "New output directory created successfully!\n",
      "Analysis will stop after having processed 18446744073709551615 fastq entries containing a cell info\n",
      "\n",
      "init models\n",
      "the log file: Mapping_log.txt\n",
      "Changing the expression start gene id to 462\n",
      "After indexing all fastq files we have the following indices:\n",
      "the mRNA index:\n",
      "I have 41039 kmers for 462 genes with 0.16375615% duplicate entries\n",
      "gene names like 'Rpl36a'\n",
      "gene_ids range from Some(0) to Some(461)\n",
      "\n",
      "the sample id index:\n",
      "I have 334 kmers for 12 genes with 0.09677419% duplicate entries\n",
      "gene names like 'Sample12'\n",
      "gene_ids range from Some(0) to Some(11)\n",
      "\n",
      "and the antibodies index:\n",
      "I have 143 kmers for 5 genes with 0% duplicate entries\n",
      "gene names like 'IgM'\n",
      "gene_ids range from Some(0) to Some(4)\n",
      "\n",
      "Writing index version 6\n",
      "with kmer_len 32\n",
      "And a total of 41039 data entries\n",
      "\n",
      "Parsing file pair 1\n",
      "\n",
      "I am using 12 cpus\n",
      "\u001b[2K\u001b[1m\u001b[2m\u001b[0m   0.50 mio reads (95.09% with cell_id, 93.94% with gene_id 0.00% multimapper)  \n",
      "\n",
      "Writing outfiles ...\n",
      "filtering cells\n",
      "Dropping cell with too little counts (n=59)\n",
      "54 cells have passed the cutoff of 200 umi counts per cell.\n",
      "\n",
      "\n",
      "writing gene expression\n",
      "sparse Matrix: 54 cell(s), 264 gene(s) and 3885 entries written to path Ok(\"Rustody_S1/BD_Rhapsody_expression\"); \n",
      "Writing Antibody counts\n",
      "sparse Matrix: 54 cell(s), 3 gene(s) and 127 entries written to path Ok(\"Rustody_S1/BD_Rhapsody_antibodies\"); \n",
      "Writing samples table\n",
      "dense matrix: 54 cell written\n",
      "\n",
      "Summary:\n",
      "cellular   reads  : 475429 reads (95.09% of total)\n",
      "no cell ID reads  : 21542 reads (4.31% of total)\n",
      "no gene ID reads  : 0 reads (0.00% of total)\n",
      "filtered   reads  : 3029 reads (0.61% of total)\n",
      " ->  multimapper  : 0 reads (0.00% of total)\n",
      " -> bad qualiity  : 3028 reads (0.61% of total)\n",
      " ->    too short  : 0 reads (0.00% of total)\n",
      " ->          N's  : 1 reads (0.00% of total)\n",
      "\n",
      "total      reads  : 500000 reads\n",
      "\n",
      "collected read counts:\n",
      "expression reads  : 465099 reads (97.83% of cellular)\n",
      "antibody reads    : 3625 reads (0.76% of cellular)\n",
      "sample reads      : 960 reads (0.20% of cellular)\n",
      "\n",
      "reported UMI counts:\n",
      "expression reads  : 190938 UMIs (40.16% of cellular)\n",
      "antibody reads    : 2950 UMIs (0.62% of cellular)\n",
      "sample reads      : 820 UMIs (0.17% of cellular)\n",
      "\n",
      "PCR duplicates or bad cells: 280721 reads (59.05% of cellular)\n",
      "\n",
      "timings:\n",
      "   overall run time 0 h 0 min 8 sec 979 millisec\n",
      "   file-io run time 0 h 0 min 0 sec 773 millisec\n",
      "single-cpu run time 0 h 0 min 0 sec 50 millisec\n",
      " multi-cpu run time 0 h 0 min 7 sec 895 millisec\n",
      "\n",
      "\n",
      "Cell->Sample table written to \"Rustody_S1/SampleCounts.tsv\"\n",
      "\n",
      "quantify_rhapsody finished in 0h 0min 8 sec 980milli sec\n",
      "\n"
     ]
    }
   ],
   "source": [
    "if os.path.exists(\"Rustody_S1\"):\n",
    "    shutil.rmtree('Rustody_S1')\n",
    "if not os.path.exists(\"Rustody_S1\"):\n",
    "    f1 =  os.path.join( \"..\", \"cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz\")\n",
    "    f2 =  os.path.join( \"..\", \"cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz\")\n",
    "    ab =  os.path.join( \"..\", \"MyAbSeqPanel.fasta\")\n",
    "    ex =  os.path.join( \"..\", \"2276_20220531_chang_to_rpl36a_amplicons.fasta\")\n",
    "    exe = os.path.join( \"..\", '..', 'target', 'release', 'quantify_rhapsody_multi' )\n",
    "    if os.name == 'nt': # windows...\n",
    "        exe = exe + '.exe'\n",
    "    print ( f\"{exe} -r {f1} -f {f2} -o Rustody_S1 -s mouse  -e {ex} -a {ab} -m 200 -v v2.96\")\n",
    "    ! {exe} -r {f1} -f {f2} -o Rustody_S1 -s mouse  -e {ex} -a {ab} -m 200 -v v2.96"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "9dec48d3-6690-4cc8-be29-77986360bc39",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'../../target/release/quantify_rhapsody_multi'"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4cd3cf0f-d7f3-4e31-9fc6-289f921a7aba",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "../../target/release/quantify_rhapsody_multi -r ../cells.1.Rhapsody_SV_index2_S2_R1_001.fastq.gz -f ../cells.1.Rhapsody_SV_index2_S2_R2_001.fastq.gz -o Rustody_S2 -s mouse  -e ../2276_20220531_chang_to_rpl36a_amplicons.fasta -a ../MyAbSeqPanel.fasta -m 200 -v v2.96\n",
      "New output directory created successfully!\n",
      "Analysis will stop after having processed 18446744073709551615 fastq entries containing a cell info\n",
      "\n",
      "init models\n",
      "the log file: Mapping_log.txt\n",
      "Changing the expression start gene id to 462\n",
      "After indexing all fastq files we have the following indices:\n",
      "the mRNA index:\n",
      "I have 41039 kmers for 462 genes with 0.16375615% duplicate entries\n",
      "gene names like 'Rpl36a'\n",
      "gene_ids range from Some(0) to Some(461)\n",
      "\n",
      "the sample id index:\n",
      "I have 334 kmers for 12 genes with 0.09677419% duplicate entries\n",
      "gene names like 'Sample12'\n",
      "gene_ids range from Some(0) to Some(11)\n",
      "\n",
      "and the antibodies index:\n",
      "I have 143 kmers for 5 genes with 0% duplicate entries\n",
      "gene names like 'IgM'\n",
      "gene_ids range from Some(0) to Some(4)\n",
      "\n",
      "Writing index version 6\n",
      "with kmer_len 32\n",
      "And a total of 41039 data entries\n",
      "\n",
      "Parsing file pair 1\n",
      "\n",
      "I am using 12 cpus\n",
      "\u001b[2K\u001b[1m\u001b[2m\u001b[0m   0.50 mio reads (94.53% with cell_id, 93.77% with gene_id 0.00% multimapper)  \n",
      "\n",
      "Writing outfiles ...\n",
      "filtering cells\n",
      "Dropping cell with too little counts (n=86)\n",
      "33 cells have passed the cutoff of 200 umi counts per cell.\n",
      "\n",
      "\n",
      "writing gene expression\n",
      "sparse Matrix: 33 cell(s), 266 gene(s) and 2530 entries written to path Ok(\"Rustody_S2/BD_Rhapsody_expression\"); \n",
      "Writing Antibody counts\n",
      "sparse Matrix: 33 cell(s), 3 gene(s) and 88 entries written to path Ok(\"Rustody_S2/BD_Rhapsody_antibodies\"); \n",
      "Writing samples table\n",
      "dense matrix: 33 cell written\n",
      "\n",
      "Summary:\n",
      "cellular   reads  : 472656 reads (94.53% of total)\n",
      "no cell ID reads  : 24114 reads (4.82% of total)\n",
      "no gene ID reads  : 0 reads (0.00% of total)\n",
      "filtered   reads  : 3230 reads (0.65% of total)\n",
      " ->  multimapper  : 0 reads (0.00% of total)\n",
      " -> bad qualiity  : 3228 reads (0.65% of total)\n",
      " ->    too short  : 0 reads (0.00% of total)\n",
      " ->          N's  : 2 reads (0.00% of total)\n",
      "\n",
      "total      reads  : 500000 reads\n",
      "\n",
      "collected read counts:\n",
      "expression reads  : 466435 reads (98.68% of cellular)\n",
      "antibody reads    : 1785 reads (0.38% of cellular)\n",
      "sample reads      : 643 reads (0.14% of cellular)\n",
      "\n",
      "reported UMI counts:\n",
      "expression reads  : 118453 UMIs (25.06% of cellular)\n",
      "antibody reads    : 1430 UMIs (0.30% of cellular)\n",
      "sample reads      : 553 UMIs (0.12% of cellular)\n",
      "\n",
      "PCR duplicates or bad cells: 352220 reads (74.52% of cellular)\n",
      "\n",
      "timings:\n",
      "   overall run time 0 h 0 min 8 sec 366 millisec\n",
      "   file-io run time 0 h 0 min 0 sec 791 millisec\n",
      "single-cpu run time 0 h 0 min 0 sec 39 millisec\n",
      " multi-cpu run time 0 h 0 min 7 sec 269 millisec\n",
      "\n",
      "\n",
      "Cell->Sample table written to \"Rustody_S2/SampleCounts.tsv\"\n",
      "\n",
      "quantify_rhapsody finished in 0h 0min 8 sec 366milli sec\n",
      "\n"
     ]
    }
   ],
   "source": [
    "if os.path.exists(\"Rustody_S2\"):\n",
    "    shutil.rmtree('Rustody_S2')\n",
    "if not os.path.exists(\"Rustody_S2\"):\n",
    "    f1 =  os.path.join( \"..\", \"cells.1.Rhapsody_SV_index2_S2_R1_001.fastq.gz\")\n",
    "    f2 =  os.path.join( \"..\", \"cells.1.Rhapsody_SV_index2_S2_R2_001.fastq.gz\")\n",
    "    ab =  os.path.join( \"..\", \"MyAbSeqPanel.fasta\")\n",
    "    ex =  os.path.join( \"..\", \"2276_20220531_chang_to_rpl36a_amplicons.fasta\")\n",
    "    print ( f\"{exe} -r {f1} -f {f2} -o Rustody_S2 -s mouse  -e {ex} -a {ab} -m 200 -v v2.96\")\n",
    "    ! {exe} -r {f1} -f {f2} -o Rustody_S2 -s mouse  -e {ex} -a {ab} -m 200 -v v2.96"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b674bcaf-cc4a-4eb4-ac9a-5bf821417fb9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def readRustodyExpression(path, name):\n",
    "    print(f\"reading Rustody expression from path {path}/BD_Rhapsody_expression/\")\n",
    "    this = scanpy.read_10x_mtx( path+'/BD_Rhapsody_expression/' )\n",
    "    this.obs['sample'] = name\n",
    "    obs1 = pd.read_csv( path+'/SampleCounts.tsv', sep=\"\\t\")\n",
    "    this.obs = this.obs.merge( obs1, left_index= True, right_on = 'CellID' )\n",
    "    this.obs_names = this.obs['CellID'] + \"_\" +  this.obs['sample']\n",
    "    this = this[this.obs['AssignedSampleName'] != \"na\"]\n",
    "    this.obs['AssignedSampleName'] = this.obs['AssignedSampleName'] + \"_\" + this.obs['sample']\n",
    "    return(this)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a9975528-4aac-415b-a518-946f1d735b94",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading Rustody expression from path Rustody_S1/BD_Rhapsody_expression/\n",
      "reading Rustody expression from path Rustody_S2/BD_Rhapsody_expression/\n"
     ]
    }
   ],
   "source": [
    "Rustody = readRustodyExpression( 'Rustody_S1', 'Index1')\n",
    "Rustody = Rustody.concatenate(readRustodyExpression( 'Rustody_S2', 'Index2'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "47a69629-d5d1-4181-afff-31cd756eb5c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "Rustody.obs['cellID'] =  [ re.findall(r'\\d+', n)[0] for n in Rustody.obs_names]\n",
    "BD_Data.obs['cellID'] =  [ re.findall(r'\\d+', n)[0] for n in BD_Data.obs_names]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "0ae0a645-7904-44a4-bf1c-6e69c62c2526",
   "metadata": {},
   "outputs": [],
   "source": [
    "Rustody.obs['cellname'] = [Rustody.obs['cellID'][id] + \"_\" + Rustody.obs['sample'][id] for id in range( len(Rustody.obs_names))]\n",
    "BD_Data.obs['cellname'] = [BD_Data.obs['cellID'][id] + \"_\" + BD_Data.obs['sample'][id] for id in range( len(BD_Data.obs_names))]\n",
    "Rustody.obs['source'] = \"Rustody\"\n",
    "BD_Data.obs['source'] = \"BD\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "83b65f9e-c2fc-437e-afbc-4c73475b5fbb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Cell_Index\n",
       "7243072-0      7243072_Index1\n",
       "4753157-0      4753157_Index1\n",
       "8267932-0      8267932_Index1\n",
       "7405530-0      7405530_Index1\n",
       "7983433-0      7983433_Index1\n",
       "                   ...       \n",
       "1509915-1      1509915_Index2\n",
       "10329694-1    10329694_Index2\n",
       "11079622-1    11079622_Index2\n",
       "2522963-1      2522963_Index2\n",
       "1799438-1      1799438_Index2\n",
       "Name: cellname, Length: 130, dtype: object"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BD_Data.obs['cellname']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "742bc8bd-dad8-49f4-be5d-ea162906e15a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'4753157_Index1' in BD_Data.obs['cellname']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "360caf58-6324-41be-9294-c8163947faa0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Cell1164_Index1-0            1164_Index1\n",
       "Cell30353_Index1-0          30353_Index1\n",
       "Cell305363_Index1-0        305363_Index1\n",
       "Cell314929_Index1-0        314929_Index1\n",
       "Cell1509915_Index1-0      1509915_Index1\n",
       "                              ...       \n",
       "Cell7972283_Index2-1      7972283_Index2\n",
       "Cell8116618_Index2-1      8116618_Index2\n",
       "Cell8143512_Index2-1      8143512_Index2\n",
       "Cell9732142_Index2-1      9732142_Index2\n",
       "Cell10329694_Index2-1    10329694_Index2\n",
       "Name: cellname, Length: 76, dtype: object"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Rustody.obs['cellname']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "cea8c3f5-ff53-4c74-ac7b-c2c8b8ceaac6",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " False,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True,\n",
       " True]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " [name  in BD_Data.obs['cellname'].values for name in Rustody.obs['cellname'].values ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "d742c735-1288-4fb2-be05-fb52b25be4e6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "View of AnnData object with n_obs × n_vars = 75 × 244\n",
       "    obs: 'sample', 'CellID', 'Sample1', 'Sample10', 'Sample11', 'Sample12', 'Sample2', 'Sample3', 'Sample4', 'Sample5', 'Sample6', 'Sample7', 'Sample8', 'Sample9', 'AssignedSampleName', 'FractionTotal', 'n', 'dist to nr.2 [%max]', 'batch', 'cellID', 'cellname', 'source'\n",
       "    var: 'gene_ids', 'feature_types'"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Rustody_4test = Rustody[  [name  in BD_Data.obs['cellname'].values for name in Rustody.obs['cellname'].values ] ]\n",
    "Rustody_4test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "5e089fbb-b325-424f-9068-fa784ffc8a7e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "View of AnnData object with n_obs × n_vars = 75 × 462\n",
       "    obs: 'sample', 'batch', 'cellID', 'cellname', 'source'"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BD_Data_4test = BD_Data[  [name  in Rustody.obs['cellname'].values for name in BD_Data.obs['cellname'].values ] ]\n",
    "BD_Data_4test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "314bca8f-13f8-4a9b-a0e9-ea90ae03ca5b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Ccnd2', 'Ctsd', 'Il15ra', 'Itgam', 'Cdkn1a', 'Mcl1', 'Ctla4', 'Anxa5',\n",
       "       'Ccr7', 'Il18',\n",
       "       ...\n",
       "       'Ybx3', 'Thbd', 'Cd28', 'Tigit', 'Fyn', 'Dock8', 'Cd22', 'Fbl', 'Ebf1',\n",
       "       'H2-DMa'],\n",
       "      dtype='object', length=244)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Rustody_4test.var_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "25e540a8-3376-424f-ab61-32c7582d846c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.04845015831115095, 0.6797635360781131)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene =  'Il18'\n",
    "pearsonr( np.array(BD_Data_4test[:,gene].X).flatten(), \n",
    "               np.array( Rustody_4test[:,gene].X.todense()).flatten() ) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "2d079d8e-5abc-44ad-951a-d5d5acf41026",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['10175248-0', '10329694-1', '10330078-0', '10481686-0', '1164-0',\n",
       "       '14029829-0', '1509915-0', '1509915-1', '157841-1', '1787564-1',\n",
       "       '1797122-0', '1799438-0', '1799438-1', '1949217-0', '2079807-0',\n",
       "       '2090208-0', '2099766-0', '25040-1', '2510987-0', '2522515-0',\n",
       "       '2522963-1', '2538660-0', '2543289-0', '30353-0', '305363-0',\n",
       "       '3103568-1', '314929-0', '3569368-1', '4283522-1', '4305474-1',\n",
       "       '4753157-0', '4885340-1', '5018926-1', '5160998-0', '5166047-1',\n",
       "       '5328779-0', '5331545-1', '5343016-0', '5463585-0', '5463585-1',\n",
       "       '5753122-1', '5784982-0', '5926310-0', '6069190-0', '6193184-1',\n",
       "       '6215132-0', '6228103-1', '6348767-1', '6367528-0', '6658620-0',\n",
       "       '6785717-0', '6800349-0', '7107103-1', '7236933-1', '7243072-0',\n",
       "       '7243433-1', '7405530-0', '7546029-1', '7693075-0', '7695391-0',\n",
       "       '7972283-1', '7983433-0', '8116618-0', '8116618-1', '8143512-1',\n",
       "       '8267932-0', '8405806-0', '8416901-0', '8713787-0', '9152684-0',\n",
       "       '9446463-0', '9473754-0', '9613135-0', '9732142-1', '9892316-0'],\n",
       "      dtype='object', name='Cell_Index')"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BD_Data_4test.obs_names.sort_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "745108ae-a757-4ac4-8600-5239466db774",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['10175248-0', '10329694-1', '10330078-0', '10481686-0', '1164-0',\n",
       "       '14029829-0', '1509915-0', '1509915-1', '157841-1', '1787564-1',\n",
       "       '1797122-0', '1799438-0', '1799438-1', '1949217-0', '2079807-0',\n",
       "       '2090208-0', '2099766-0', '25040-1', '2510987-0', '2522515-0',\n",
       "       '2522963-1', '2538660-0', '2543289-0', '30353-0', '305363-0',\n",
       "       '3103568-1', '314929-0', '3569368-1', '4283522-1', '4305474-1',\n",
       "       '4753157-0', '4885340-1', '5018926-1', '5160998-0', '5166047-1',\n",
       "       '5328779-0', '5331545-1', '5343016-0', '5463585-0', '5463585-1',\n",
       "       '5753122-1', '5784982-0', '5926310-0', '6069190-0', '6193184-1',\n",
       "       '6215132-0', '6228103-1', '6348767-1', '6367528-0', '6658620-0',\n",
       "       '6785717-0', '6800349-0', '7107103-1', '7236933-1', '7243072-0',\n",
       "       '7243433-1', '7405530-0', '7546029-1', '7693075-0', '7695391-0',\n",
       "       '7972283-1', '7983433-0', '8116618-0', '8116618-1', '8143512-1',\n",
       "       '8267932-0', '8405806-0', '8416901-0', '8713787-0', '9152684-0',\n",
       "       '9446463-0', '9473754-0', '9613135-0', '9732142-1', '9892316-0'],\n",
       "      dtype='object', name='Cell_Index')"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BD_Data_4test = BD_Data_4test[ BD_Data_4test.obs_names.sort_values() ]\n",
    "BD_Data_4test.obs_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "f8aa1449-0693-4370-ab48-b61eb64ef200",
   "metadata": {},
   "outputs": [],
   "source": [
    "Rustody_4test = Rustody_4test[Rustody_4test.obs_names.sort_values()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "e975d77c-38b2-409c-a0ed-f543454d8074",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0.], dtype=float32)"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(BD_Data_4test[:,gene].X).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "f22ee3b9-159c-41d4-9e0e-6e5a39e6bd09",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 0., 1., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0.], dtype=float32)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(Rustody_4test[:,gene].X.todense()).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "4f204295-2604-494a-bdbf-464ce16a7e11",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>gname</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>Ccnd2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.973065</td>\n",
       "      <td>2.879407e-48</td>\n",
       "      <td>Ctsd</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.929403</td>\n",
       "      <td>2.502492e-33</td>\n",
       "      <td>Il15ra</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>Itgam</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>Cdkn1a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>239</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>Dock8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>240</th>\n",
       "      <td>0.985371</td>\n",
       "      <td>7.531320e-58</td>\n",
       "      <td>Cd22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>241</th>\n",
       "      <td>0.983262</td>\n",
       "      <td>9.896981e-56</td>\n",
       "      <td>Fbl</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>242</th>\n",
       "      <td>0.739830</td>\n",
       "      <td>3.391389e-14</td>\n",
       "      <td>Ebf1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>243</th>\n",
       "      <td>0.996584</td>\n",
       "      <td>8.011100e-81</td>\n",
       "      <td>H2-DMa</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>244 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            0             1   gname\n",
       "0    1.000000  0.000000e+00   Ccnd2\n",
       "1    0.973065  2.879407e-48    Ctsd\n",
       "2    0.929403  2.502492e-33  Il15ra\n",
       "3    1.000000  0.000000e+00   Itgam\n",
       "4    1.000000  0.000000e+00  Cdkn1a\n",
       "..        ...           ...     ...\n",
       "239  1.000000  0.000000e+00   Dock8\n",
       "240  0.985371  7.531320e-58    Cd22\n",
       "241  0.983262  9.896981e-56     Fbl\n",
       "242  0.739830  3.391389e-14    Ebf1\n",
       "243  0.996584  8.011100e-81  H2-DMa\n",
       "\n",
       "[244 rows x 3 columns]"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats.stats import pearsonr   \n",
    "cmp = pd.DataFrame(\n",
    "    [ pearsonr( np.array(BD_Data_4test[:,gene].X).flatten(), \n",
    "               np.array( Rustody_4test[:,gene].X.todense()).flatten() ) for gene in Rustody_4test.var_names if gene in BD_Data_4test.var_names ]  )\n",
    "cmp['gname'] = [gene  for gene in Rustody_4test.var_names if gene in BD_Data_4test.var_names ]\n",
    "cmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "6a57c025-ff59-4ea0-8f08-f5b741d79bc3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[59]"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vals = np.array(BD_Data_4test[:,'Ebf1'].X).flatten() == 2\n",
    "[ id for id in range( len(vals) ) if vals[id] ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "28f2cd5f-74e7-4a77-a39f-9cb41439ac58",
   "metadata": {},
   "outputs": [],
   "source": [
    "sorted_var_names = sorted(Rustody_4test.var_names)\n",
    "BD_Data_4test = BD_Data_4test[:, sorted_var_names]\n",
    "Rustody_4test = Rustody_4test[:, sorted_var_names]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "933ef027-4a64-4dca-93c8-a8ade33d31d6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.9604759666203988, 2.7698182615922647e-42)"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pearsonr( np.array(Rustody_4test [:, 59].X.toarray().flatten()),\n",
    "       np.array(BD_Data_4test [:, 59].X.toarray().flatten()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "6407fd57-bb02-42c0-b21d-582121838d52",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 1., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 2., 0., 0.,\n",
       "       0., 0., 0., 1., 0., 0., 1., 2., 0., 1., 0., 2., 0., 0., 0., 1., 0.,\n",
       "       0., 0., 0., 0., 1., 2., 0., 1., 0., 0., 2., 0., 0., 1., 0., 0., 1.,\n",
       "       0., 0., 0., 2., 0., 2., 0., 1., 0., 1., 1., 0., 2., 0., 1., 0., 0.,\n",
       "       0., 0., 0., 0., 2., 0., 0.], dtype=float32)"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(Rustody_4test [:, 59].X.toarray().flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "17bfa13c-08be-4f1e-a5d8-3220ec68fcc8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 0., 0., 2., 0., 0.,\n",
       "       0., 0., 0., 1., 0., 0., 1., 2., 0., 1., 0., 3., 0., 0., 0., 1., 0.,\n",
       "       0., 0., 0., 0., 1., 2., 0., 1., 0., 0., 2., 0., 0., 1., 0., 0., 1.,\n",
       "       0., 0., 0., 1., 0., 2., 0., 1., 0., 1., 1., 0., 2., 0., 1., 0., 0.,\n",
       "       0., 0., 0., 0., 2., 0., 0.], dtype=float32)"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(BD_Data_4test [:, 59].X.toarray().flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "026eedfb-5c73-4246-b610-fab9ca87cff3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>gname</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.994028</td>\n",
       "      <td>5.487680e-72</td>\n",
       "      <td>2810417H13Rik</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.990703</td>\n",
       "      <td>5.384115e-65</td>\n",
       "      <td>Ada</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.955216</td>\n",
       "      <td>2.413487e-40</td>\n",
       "      <td>Anxa5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.997755</td>\n",
       "      <td>1.805195e-87</td>\n",
       "      <td>Apoe</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.980786</td>\n",
       "      <td>1.458155e-53</td>\n",
       "      <td>Arid3a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>239</th>\n",
       "      <td>0.996302</td>\n",
       "      <td>1.447880e-79</td>\n",
       "      <td>Vpreb3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>240</th>\n",
       "      <td>0.979166</td>\n",
       "      <td>2.720688e-52</td>\n",
       "      <td>Vps28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>241</th>\n",
       "      <td>0.998216</td>\n",
       "      <td>4.153936e-91</td>\n",
       "      <td>Xbp1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>242</th>\n",
       "      <td>0.979472</td>\n",
       "      <td>1.594057e-52</td>\n",
       "      <td>Ybx3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>243</th>\n",
       "      <td>0.919087</td>\n",
       "      <td>3.020247e-31</td>\n",
       "      <td>tdTomato</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>244 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            0             1          gname\n",
       "0    0.994028  5.487680e-72  2810417H13Rik\n",
       "1    0.990703  5.384115e-65            Ada\n",
       "2    0.955216  2.413487e-40          Anxa5\n",
       "3    0.997755  1.805195e-87           Apoe\n",
       "4    0.980786  1.458155e-53         Arid3a\n",
       "..        ...           ...            ...\n",
       "239  0.996302  1.447880e-79         Vpreb3\n",
       "240  0.979166  2.720688e-52          Vps28\n",
       "241  0.998216  4.153936e-91           Xbp1\n",
       "242  0.979472  1.594057e-52           Ybx3\n",
       "243  0.919087  3.020247e-31       tdTomato\n",
       "\n",
       "[244 rows x 3 columns]"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats.stats import pearsonr   \n",
    "cmp = pd.DataFrame(\n",
    "    [ pearsonr( np.array(BD_Data_4test[:,gene].X).flatten(), \n",
    "               np.array( Rustody_4test[:,gene].X.todense()).flatten() ) for gene in Rustody_4test.var_names if gene in BD_Data_4test.var_names ]  )\n",
    "cmp['gname'] = [gene  for gene in Rustody_4test.var_names if gene in BD_Data_4test.var_names ]\n",
    "cmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "e02b0742-c497-4ef1-b3aa-9cbfcd50b6ed",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>gname</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>0.336416</td>\n",
       "      <td>3.167011e-03</td>\n",
       "      <td>Bcl2l11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>0.216720</td>\n",
       "      <td>6.181983e-02</td>\n",
       "      <td>Cd36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>0.127551</td>\n",
       "      <td>2.754801e-01</td>\n",
       "      <td>Cd37</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>0.008318</td>\n",
       "      <td>9.435359e-01</td>\n",
       "      <td>Cd74</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>60</th>\n",
       "      <td>0.866443</td>\n",
       "      <td>1.018606e-23</td>\n",
       "      <td>Cnot2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>65</th>\n",
       "      <td>0.218218</td>\n",
       "      <td>5.999690e-02</td>\n",
       "      <td>Cxcr4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70</th>\n",
       "      <td>0.755025</td>\n",
       "      <td>5.079868e-15</td>\n",
       "      <td>Dpp4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>73</th>\n",
       "      <td>0.739830</td>\n",
       "      <td>3.391389e-14</td>\n",
       "      <td>Ebf1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>83</th>\n",
       "      <td>0.696733</td>\n",
       "      <td>3.866365e-12</td>\n",
       "      <td>Fbxl12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>125</th>\n",
       "      <td>0.859990</td>\n",
       "      <td>5.060482e-23</td>\n",
       "      <td>Il12a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>127</th>\n",
       "      <td>0.859990</td>\n",
       "      <td>5.060482e-23</td>\n",
       "      <td>Il18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>157</th>\n",
       "      <td>0.489760</td>\n",
       "      <td>8.214055e-06</td>\n",
       "      <td>Ltb</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>184</th>\n",
       "      <td>0.888106</td>\n",
       "      <td>2.372411e-26</td>\n",
       "      <td>Pmaip1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>208</th>\n",
       "      <td>0.893764</td>\n",
       "      <td>3.957065e-27</td>\n",
       "      <td>Syk</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>217</th>\n",
       "      <td>0.772904</td>\n",
       "      <td>4.523571e-16</td>\n",
       "      <td>Tlr2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>218</th>\n",
       "      <td>0.862062</td>\n",
       "      <td>3.051655e-23</td>\n",
       "      <td>Tlr4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>221</th>\n",
       "      <td>0.874022</td>\n",
       "      <td>1.388275e-24</td>\n",
       "      <td>Tmem97</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>229</th>\n",
       "      <td>0.810885</td>\n",
       "      <td>1.175957e-18</td>\n",
       "      <td>Trac</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>235</th>\n",
       "      <td>0.847148</td>\n",
       "      <td>9.812321e-22</td>\n",
       "      <td>Tyms</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            0             1    gname\n",
       "14   0.336416  3.167011e-03  Bcl2l11\n",
       "41   0.216720  6.181983e-02     Cd36\n",
       "42   0.127551  2.754801e-01     Cd37\n",
       "50   0.008318  9.435359e-01     Cd74\n",
       "60   0.866443  1.018606e-23    Cnot2\n",
       "65   0.218218  5.999690e-02    Cxcr4\n",
       "70   0.755025  5.079868e-15     Dpp4\n",
       "73   0.739830  3.391389e-14     Ebf1\n",
       "83   0.696733  3.866365e-12   Fbxl12\n",
       "125  0.859990  5.060482e-23    Il12a\n",
       "127  0.859990  5.060482e-23     Il18\n",
       "157  0.489760  8.214055e-06      Ltb\n",
       "184  0.888106  2.372411e-26   Pmaip1\n",
       "208  0.893764  3.957065e-27      Syk\n",
       "217  0.772904  4.523571e-16     Tlr2\n",
       "218  0.862062  3.051655e-23     Tlr4\n",
       "221  0.874022  1.388275e-24   Tmem97\n",
       "229  0.810885  1.175957e-18     Trac\n",
       "235  0.847148  9.812321e-22     Tyms"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "not_good = [id for id in range( len( cmp[0])) if cmp[0][id] < 0.9]\n",
    "cmp.iloc[not_good, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "9a798dce-6496-4821-abda-5fa762e94357",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BD\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0.]\n",
      "Rustody\n",
      "[0. 0. 0. 0. 1. 0. 0. 1. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 4. 0. 0. 1. 0. 2. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 2. 0. 0. 0. 0. 0. 1. 0. 0. 0. 1. 0.\n",
      " 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "gene  = 'Bcl2l11'\n",
    "print(\"BD\")\n",
    "print (np.array(BD_Data_4test[:,gene].X).flatten() )\n",
    "print(\"Rustody\")\n",
    "print(np.array(Rustody_4test[:,gene].X.todense()).flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "0223688b-ced2-4176-b4b1-a761a167514a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({'Index1': 47, 'Index2': 28})"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Counter( Rustody_4test.obs['sample'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "51eded4a-b800-4241-b446-10ae72a7ab18",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'adata' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Input \u001b[0;32mIn [80]\u001b[0m, in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mplotly\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mexpress\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpx\u001b[39;00m\n\u001b[1;32m      2\u001b[0m df \u001b[38;5;241m=\u001b[39m pd\u001b[38;5;241m.\u001b[39mDataFrame( \n\u001b[1;32m      3\u001b[0m     { \n\u001b[0;32m----> 4\u001b[0m         \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mBD_I1\u001b[39m\u001b[38;5;124m'\u001b[39m :\u001b[43madata\u001b[49m[adata\u001b[38;5;241m.\u001b[39mobs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msample\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mS1_BD\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mobs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlog1p_total_counts\u001b[39m\u001b[38;5;124m'\u001b[39m],\n\u001b[1;32m      5\u001b[0m         \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mBD_I2\u001b[39m\u001b[38;5;124m'\u001b[39m :adata[adata\u001b[38;5;241m.\u001b[39mobs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msample\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mS2_BD\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mobs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlog1p_total_counts\u001b[39m\u001b[38;5;124m'\u001b[39m],\n\u001b[1;32m      6\u001b[0m         \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mRustody_I1\u001b[39m\u001b[38;5;124m'\u001b[39m :adata[adata\u001b[38;5;241m.\u001b[39mobs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msample\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mS1_Rustody\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mobs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlog1p_total_counts\u001b[39m\u001b[38;5;124m'\u001b[39m],\n\u001b[1;32m      7\u001b[0m         \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mRustody_I2\u001b[39m\u001b[38;5;124m'\u001b[39m :adata[adata\u001b[38;5;241m.\u001b[39mobs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msample\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mS2_Rustody\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mobs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlog1p_total_counts\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[1;32m      8\u001b[0m     })\n\u001b[1;32m      9\u001b[0m fig \u001b[38;5;241m=\u001b[39m px\u001b[38;5;241m.\u001b[39mviolin(df, y\u001b[38;5;241m=\u001b[39m df\u001b[38;5;241m.\u001b[39mcolumns )\n\u001b[1;32m     10\u001b[0m fig\u001b[38;5;241m.\u001b[39mshow()\n",
      "\u001b[0;31mNameError\u001b[0m: name 'adata' is not defined"
     ]
    }
   ],
   "source": [
    "import plotly.express as px\n",
    "df = pd.DataFrame( \n",
    "    { \n",
    "        'BD_I1' :adata[adata.obs['sample'] == 'S1_BD'].obs['log1p_total_counts'],\n",
    "        'BD_I2' :adata[adata.obs['sample'] == 'S2_BD'].obs['log1p_total_counts'],\n",
    "        'Rustody_I1' :adata[adata.obs['sample'] == 'S1_Rustody'].obs['log1p_total_counts'],\n",
    "        'Rustody_I2' :adata[adata.obs['sample'] == 'S2_Rustody'].obs['log1p_total_counts']\n",
    "    })\n",
    "fig = px.violin(df, y= df.columns )\n",
    "fig.show() \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72f5c108-9e72-4054-8736-c93617a08946",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata[:,\"Cnot2\"][adata.obs['sample'] == 'S1_Rustody'].X.todense().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00dc46e1-5d13-4a05-b38e-b59464dd61c0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "scanpy.pp.filter_genes(adata, min_counts=1 )\n",
    "scanpy.pp.filter_cells(adata, min_genes=20 )\n",
    "Counter( adata.obs['sample'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da9e0338-ef65-4c6f-a28c-e676fc63fe68",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata.obs['cellID'] =  [ re.findall(r'\\d+', n)[0] for n in adata.obs_names] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66afb40d-67a8-4617-a2e3-dfcd8debdb96",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs['cellID']  = adata.obs['cellID'].astype('str')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f611ccf-bb4d-41be-b50d-3814f17c4e68",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "Counter([ name in adata[adata.obs['sample'] == 'S1_BD'].obs['cellID'].values for name in adata[adata.obs['sample'] == 'S1_Rustody'].obs['cellID'] ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edf25c31-7e4d-472e-a85f-5fcfc730932e",
   "metadata": {},
   "outputs": [],
   "source": [
    "Counter([ name in adata[adata.obs['sample'] == 'S1_BD'].obs['cellID'].values \n",
    "         for name in adata[adata.obs['sample'] == 'S1_Rustody'].obs['cellID'] ])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c28d7db6-9598-49bf-bbdb-7f17ba28aa51",
   "metadata": {},
   "outputs": [],
   "source": [
    "Counter([ name in adata[adata.obs['sample'] == 'S1_Rustody'].obs['cellID'].values \n",
    "         for name in adata[adata.obs['sample'] == 'S1_BD'].obs['cellID'] ])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7db9f5b1-33c9-441a-b8ef-6570e847f851",
   "metadata": {},
   "outputs": [],
   "source": [
    "[[ name, name in adata[adata.obs['sample'] == 'S1_Rustody'].obs['cellID'].values ]\n",
    "         for name in adata[adata.obs['sample'] == 'S1_BD'].obs['cellID'] ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bddb2a8d-157c-4dad-aafb-4f34d2a7c006",
   "metadata": {},
   "outputs": [],
   "source": [
    "! bd_cell_id_2_seq --id 12842538 --version v2.96"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a18000d-af90-4ed6-af27-6a79cf2578f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata[adata.obs['sample'] == 'S1_Rustody'].obs['cellID'] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6c33230-54e2-4c1b-99f9-2909af337adf",
   "metadata": {},
   "outputs": [],
   "source": [
    ".astype('str').astype('str').astype('str')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c332c81-6379-435a-bd00-13cf8d29f24e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"1164\" in adata[adata.obs['sample'] == 'S1_BD'].obs['cellID'].values "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1f07051-1523-4fc7-92dd-a6935c2e2001",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dat = adata[ adata.obs['cellID'] == '4753157' ].X.todense()\n",
    "ids = [id  for id in range(dat.shape[1]) if not dat[0,id] == dat[1,id] ]\n",
    "bools = [ not dat[0,id] == dat[1,id]  for id in range(dat.shape[1]) ]\n",
    "adata.var.loc[bools,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2a79f3e-1736-41b9-888e-af806f845770",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f1f0cb6-4360-4388-8fd2-098575db4ebc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "diff_val_BD = [ int(n) for n in dat[0,bools].transpose()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2efd7e6-d4ef-45eb-acb7-7734c1482539",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "diff_val_Rhapsody = [ int(n) for n in dat[1,bools].transpose()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ac465ec-279c-412b-99c5-eebbd9fee348",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import plotly.express as px\n",
    "df = pd.DataFrame( { 'BD' :np.log1p(diff_val_BD), 'Rhapsody' : np.log1p(diff_val_Rhapsody)})\n",
    "fig = px.violin(df, y= df.columns )\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe7438f6-49cc-4014-b08c-c7a7aef8d212",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata[ adata.obs['cellID'] == '4753157' ].obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f54db747-adeb-488e-a20b-941662a72dc5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dat.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56ef3f02-af9d-49b1-a30c-80d3a880e426",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "scanpy.pp.calculate_qc_metrics( adata, inplace=True, log1p=True, percent_top= [10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69221251-f1bd-4be8-9c04-cd650be0582b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "scanpy.pp.filter_cells(adata, min_counts=400 )\n",
    "scanpy.pp.downsample_counts(adata, counts_per_cell= 400 )\n",
    "scv.pp.log1p(adata)\n",
    "Counter( adata.obs['sample'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dd341b8-f1f2-4879-aeaf-51c923878d32",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "scanpy.pp.neighbors(adata)\n",
    "dimensions = 2\n",
    "scanpy.tl.umap(adata,n_components= dimensions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06c206b5-4710-48a3-84d6-c7f3be3ee84e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "scv.pl.scatter(adata, color='sample', figsize =(7,5), legend_loc='right margin')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f158ec37-87e1-482c-8be7-e6d660f5b4ff",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "scanpy.tl.louvain(adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0312703-8da0-4758-a17b-8eb2ef50b1c3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "scv.pl.scatter(adata, color='louvain', figsize =(7,5), legend_loc='right margin')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da18417a-caf6-4eff-8150-66f1d97f0928",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62f72d31-4be7-4c67-8452-79e249230440",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata.obs.pivot_table(values = \"n_counts\", index = \"louvain\", columns=\"sample\", aggfunc='count')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea692deb-f54d-4a8f-b903-f702acfd2dcf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata.obs['n_genes' ] = adata.obs['n_genes' ].astype('int32')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0876b73e-0656-4e31-a6dd-2bb12940f67e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa7182fc-33e4-4c37-a798-fa61222aad18",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "scanpy.pl.violin(adata, [ 'n_genes', 'log1p_total_counts' ], 'sample')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd04f124-2b4e-4556-8363-8ed2949f0098",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "key_added = \"louvain\"\n",
    "scanpy.tl.rank_genes_groups(\n",
    "    adata, \n",
    "    groupby   = 'louvain',\n",
    "    key_added = key_added,\n",
    "    method    = 'wilcoxon',\n",
    ")\n",
    "\n",
    "scanpy.pl.rank_genes_groups(adata, key = key_added )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fbad0e3-5df0-40d6-b74f-5b107c281a34",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "scv.pl.scatter(adata, color=['Igkc', 'Igha', 'Ighm'], figsize =(7,5), legend_loc='right margin')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12a08c30-f4ad-4079-b4e3-3eb976ca3aba",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "scv.pl.scatter(adata, color=[ 'Ighg1', 'Ighg2b', 'Ighg3' ], figsize =(7,5), legend_loc='right margin')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7398ab7-8745-417a-ae7f-a3f059256dee",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "names = np.array(adata.obs['sample'].unique()).astype('str')\n",
    "names.sort()\n",
    "names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f7263c6-5b9c-4c91-9385-f214caf88692",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for n in names:\n",
    "    print (n)\n",
    "    scv.pl.scatter(adata[adata.obs['sample'] == n], color=[ 'Ighg1', 'Ighg2b', 'Ighg3' ], figsize =(7,5), legend_loc='right margin')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8400da5-f8a8-4242-9d5b-c6f3bccfa162",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for n in names:\n",
    "    print (n)\n",
    "    scv.pl.scatter(adata[adata.obs['sample'] == n], color=[ 'Top2a', 'Tmem97', 'Jchain' ], figsize =(7,5), legend_loc='right margin')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "021a7881-9681-48d2-9454-57971416f4f7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e13bf18a-e405-4a2e-b656-99363a0eefea",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "415bd241-abe1-40ef-92e4-718de0717890",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}