{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scanpy==1.6.0 anndata==0.7.4 umap==0.4.6 numpy==1.18.5 scipy==1.4.1 pandas==1.1.3 scikit-learn==0.23.2 statsmodels==0.12.0 python-igraph==0.7.1 leidenalg==0.8.2\n"
     ]
    }
   ],
   "source": [
    "# import the transcriptome data\n",
    "import os\n",
    "import dandelion as ddl\n",
    "import pandas as pd\n",
    "import scanpy as sc\n",
    "sc.logging.print_header()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dandelion==0.0.27.post2 pandas==1.1.3 numpy==1.18.5 matplotlib==3.3.2 networkx==2.5 scipy==1.4.1 skbio==0.5.6\n"
     ]
    }
   ],
   "source": [
    "ddl.logging.print_versions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir('/lustre/scratch117/cellgen/team297/kt16/newcastle_covid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in full vdj object\n",
    "vdj = ddl.read_h5('dandelion_output/covid_jan_2021_bcells_vdj.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dandelion class object with n_obs = 47489 and n_contigs = 96500\n",
       "    data: 'sequence_id', 'sequence', 'rev_comp', 'productive', 'v_call', 'd_call', 'j_call', 'sequence_alignment', 'germline_alignment', 'junction', 'junction_aa', 'v_cigar', 'd_cigar', 'j_cigar', 'stop_codon', 'vj_in_frame', 'locus', 'junction_length', 'np1_length', 'np2_length', 'v_sequence_start', 'v_sequence_end', 'v_germline_start', 'v_germline_end', 'd_sequence_start', 'd_sequence_end', 'd_germline_start', 'd_germline_end', 'j_sequence_start', 'j_sequence_end', 'j_germline_start', 'j_germline_end', 'v_score', 'v_identity', 'v_support', 'd_score', 'd_identity', 'd_support', 'j_score', 'j_identity', 'j_support', 'fwr1', 'fwr2', 'fwr3', 'fwr4', 'cdr1', 'cdr2', 'cdr3', 'cell_id', 'c_call', 'consensus_count', 'umi_count', 'v_call_10x', 'd_call_10x', 'j_call_10x', 'junction_10x', 'junction_10x_aa', 'v_call_genotyped', 'germline_alignment_d_mask', 'sample_id', 'c_sequence_alignment', 'c_germline_alignment', 'c_sequence_start', 'c_sequence_end', 'c_score', 'c_identity', 'c_support', 'c_call_10x', 'junction_aa_length', 'fwr1_aa', 'fwr2_aa', 'fwr3_aa', 'fwr4_aa', 'cdr1_aa', 'cdr2_aa', 'cdr3_aa', 'sequence_alignment_aa', 'v_sequence_alignment_aa', 'd_sequence_alignment_aa', 'j_sequence_alignment_aa', 'mu_freq', 'duplicate_count'\n",
       "    metadata: 'sample_id', 'isotype', 'lightchain', 'status', 'vdj_status', 'productive', 'umi_counts_heavy', 'umi_counts_light', 'v_call_heavy', 'v_call_light', 'j_call_heavy', 'j_call_light', 'c_call_heavy', 'c_call_light'\n",
       "    distance: None\n",
       "    edges: None\n",
       "    layout: None\n",
       "    graph: None"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vdj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AnnData object with n_obs × n_vars = 74437 × 1051\n",
       "    obs: 'Site', 'doublet', 'patient_id', 'sample_id', 'batch', 'Resample', 'Collection_Day', 'Sex', 'Age', 'Swab_result', 'Status', 'Smoker', 'Status_on_day_collection', 'Status_on_day_collection_summary', 'Days_from_onset', 'time_after_LPS', 'Worst_Clinical_Status', 'Outcome', 'leiden', 'consensus', 'initial_clustering', 'study_id', 'initial_clustering_B', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'leiden_B', 'celltype_B', 'celltype_B_v2', 'Status_on_day_collection_summary_v2', 'filter_rna', 'has_bcr', 'filter_bcr_quality', 'filter_bcr_heavy', 'filter_bcr_light', 'bcr_QC_pass', 'filter_bcr'\n",
       "    var: 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'\n",
       "    uns: 'Site_colors', 'Status_on_day_collection_summary_colors', 'Status_on_day_collection_summary_v2_colors', 'bcr_QC_pass_colors', 'celltype_B_colors', 'celltype_B_v2_colors', 'consensus_colors', \"dendrogram_['leiden_B']\", 'filter_bcr_colors', 'hvg', 'initial_clustering_colors', 'leiden', 'leiden_B_colors', 'leiden_colors', 'neighbors', 'pca', 'rank_genes_groups', 'umap'\n",
       "    obsm: 'X_pca', 'X_pca_harmony', 'X_umap'\n",
       "    varm: 'PCs'\n",
       "    obsp: 'connectivities', 'distances'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# read in B cell object\n",
    "bdata = sc.read_h5ad('h5ad/covid_jan_2021_bcells.h5ad')\n",
    "bdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nfs/users/nfs_k/kt16/miniconda2/envs/dandelion/lib/python3.7/site-packages/anndata/_core/anndata.py:1094: FutureWarning: is_categorical is deprecated and will be removed in a future version.  Use is_categorical_dtype instead\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "AnnData object with n_obs × n_vars = 68609 × 1051\n",
       "    obs: 'Site', 'doublet', 'patient_id', 'sample_id', 'batch', 'Resample', 'Collection_Day', 'Sex', 'Age', 'Swab_result', 'Status', 'Smoker', 'Status_on_day_collection', 'Status_on_day_collection_summary', 'Days_from_onset', 'time_after_LPS', 'Worst_Clinical_Status', 'Outcome', 'leiden', 'consensus', 'initial_clustering', 'study_id', 'initial_clustering_B', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'leiden_B', 'celltype_B', 'celltype_B_v2', 'Status_on_day_collection_summary_v2', 'filter_rna', 'has_bcr', 'filter_bcr_quality', 'filter_bcr_heavy', 'filter_bcr_light', 'bcr_QC_pass', 'filter_bcr'\n",
       "    var: 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'\n",
       "    uns: 'Site_colors', 'Status_on_day_collection_summary_colors', 'Status_on_day_collection_summary_v2_colors', 'celltype_B_colors', 'celltype_B_v2_colors', 'consensus_colors', \"dendrogram_['leiden_B']\", 'filter_bcr_colors', 'hvg', 'initial_clustering_colors', 'leiden', 'leiden_B_colors', 'leiden_colors', 'neighbors', 'pca', 'rank_genes_groups', 'umap'\n",
       "    obsm: 'X_pca', 'X_pca_harmony', 'X_umap'\n",
       "    varm: 'PCs'\n",
       "    obsp: 'connectivities', 'distances'"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# subset to only D0 samples\n",
    "bdata2 = bdata[bdata.obs['Collection_Day'] == 'D0'].copy()\n",
    "bdata2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AnnData object with n_obs × n_vars = 43917 × 1051\n",
       "    obs: 'Site', 'doublet', 'patient_id', 'sample_id', 'batch', 'Resample', 'Collection_Day', 'Sex', 'Age', 'Swab_result', 'Status', 'Smoker', 'Status_on_day_collection', 'Status_on_day_collection_summary', 'Days_from_onset', 'time_after_LPS', 'Worst_Clinical_Status', 'Outcome', 'leiden', 'consensus', 'initial_clustering', 'study_id', 'initial_clustering_B', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'leiden_B', 'celltype_B', 'celltype_B_v2', 'Status_on_day_collection_summary_v2', 'filter_rna', 'has_bcr', 'filter_bcr_quality', 'filter_bcr_heavy', 'filter_bcr_light', 'bcr_QC_pass', 'filter_bcr'\n",
       "    var: 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'\n",
       "    uns: 'Site_colors', 'Status_on_day_collection_summary_colors', 'Status_on_day_collection_summary_v2_colors', 'celltype_B_colors', 'celltype_B_v2_colors', 'consensus_colors', \"dendrogram_['leiden_B']\", 'filter_bcr_colors', 'hvg', 'initial_clustering_colors', 'leiden', 'leiden_B_colors', 'leiden_colors', 'neighbors', 'pca', 'rank_genes_groups', 'umap'\n",
       "    obsm: 'X_pca', 'X_pca_harmony', 'X_umap'\n",
       "    varm: 'PCs'\n",
       "    obsp: 'connectivities', 'distances'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bdata2 = bdata2[bdata2.obs['bcr_QC_pass'] == 'True'].copy()\n",
    "bdata2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dandelion class object with n_obs = 43917 and n_contigs = 89311\n",
       "    data: 'sequence_id', 'sequence', 'rev_comp', 'productive', 'v_call', 'd_call', 'j_call', 'sequence_alignment', 'germline_alignment', 'junction', 'junction_aa', 'v_cigar', 'd_cigar', 'j_cigar', 'stop_codon', 'vj_in_frame', 'locus', 'junction_length', 'np1_length', 'np2_length', 'v_sequence_start', 'v_sequence_end', 'v_germline_start', 'v_germline_end', 'd_sequence_start', 'd_sequence_end', 'd_germline_start', 'd_germline_end', 'j_sequence_start', 'j_sequence_end', 'j_germline_start', 'j_germline_end', 'v_score', 'v_identity', 'v_support', 'd_score', 'd_identity', 'd_support', 'j_score', 'j_identity', 'j_support', 'fwr1', 'fwr2', 'fwr3', 'fwr4', 'cdr1', 'cdr2', 'cdr3', 'cell_id', 'c_call', 'consensus_count', 'umi_count', 'v_call_10x', 'd_call_10x', 'j_call_10x', 'junction_10x', 'junction_10x_aa', 'v_call_genotyped', 'germline_alignment_d_mask', 'sample_id', 'c_sequence_alignment', 'c_germline_alignment', 'c_sequence_start', 'c_sequence_end', 'c_score', 'c_identity', 'c_support', 'c_call_10x', 'junction_aa_length', 'fwr1_aa', 'fwr2_aa', 'fwr3_aa', 'fwr4_aa', 'cdr1_aa', 'cdr2_aa', 'cdr3_aa', 'sequence_alignment_aa', 'v_sequence_alignment_aa', 'd_sequence_alignment_aa', 'j_sequence_alignment_aa', 'mu_freq', 'duplicate_count'\n",
       "    metadata: 'sample_id', 'isotype', 'lightchain', 'status', 'vdj_status', 'productive', 'umi_counts_heavy', 'umi_counts_light', 'v_call_heavy', 'v_call_light', 'j_call_heavy', 'j_call_light', 'c_call_heavy', 'c_call_light'\n",
       "    distance: None\n",
       "    edges: None\n",
       "    layout: None\n",
       "    graph: None"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# subset vdj data to only cells found in the B cell transcriptome data\n",
    "vdj2 = ddl.Dandelion(vdj.data[vdj.data['cell_id'].isin(list(bdata2.obs_names))].copy())\n",
    "vdj2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Finding clones based on heavy chains : 100%|██████████| 472/472 [00:05<00:00, 79.58it/s] \n",
      "Refining clone assignment based on light chain pairing : 100%|██████████| 40259/40259 [02:24<00:00, 279.34it/s]\n"
     ]
    }
   ],
   "source": [
    "ddl.tl.find_clones(vdj2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "vdj2.write_h5('dandelion_output/covid_jan_2021_bcells_vdj.h5', complib = 'blosc:lz4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['LPS',\n",
       " 'Non_covid',\n",
       " 'Critical',\n",
       " 'Healthy',\n",
       " 'Asymptomatic',\n",
       " 'Moderate',\n",
       " 'Mild',\n",
       " 'Severe']"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# partition into Status_on_day_collection_summary\n",
    "list(set(bdata2.obs['Status_on_day_collection_summary_v2']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Moderate        11707\n",
       "Mild            10776\n",
       "Severe           7131\n",
       "Healthy          5372\n",
       "Critical         4047\n",
       "LPS              2148\n",
       "Asymptomatic     1762\n",
       "Non_covid         974\n",
       "Name: Status_on_day_collection_summary_v2, dtype: int64"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bdata2.obs['Status_on_day_collection_summary_v2'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AnnData object with n_obs × n_vars = 43917 × 1051\n",
       "    obs: 'Site', 'doublet', 'patient_id', 'sample_id', 'batch', 'Resample', 'Collection_Day', 'Sex', 'Age', 'Swab_result', 'Status', 'Smoker', 'Status_on_day_collection', 'Status_on_day_collection_summary', 'Days_from_onset', 'time_after_LPS', 'Worst_Clinical_Status', 'Outcome', 'leiden', 'consensus', 'initial_clustering', 'study_id', 'initial_clustering_B', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'leiden_B', 'celltype_B', 'celltype_B_v2', 'Status_on_day_collection_summary_v2', 'filter_rna', 'has_bcr', 'filter_bcr_quality', 'filter_bcr_heavy', 'filter_bcr_light', 'bcr_QC_pass', 'filter_bcr'\n",
       "    var: 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'\n",
       "    uns: 'Site_colors', 'Status_on_day_collection_summary_colors', 'Status_on_day_collection_summary_v2_colors', 'celltype_B_colors', 'celltype_B_v2_colors', 'consensus_colors', \"dendrogram_['leiden_B']\", 'filter_bcr_colors', 'hvg', 'initial_clustering_colors', 'leiden', 'leiden_B_colors', 'leiden_colors', 'neighbors', 'pca', 'rank_genes_groups', 'umap'\n",
       "    obsm: 'X_pca', 'X_pca_harmony', 'X_umap'\n",
       "    varm: 'PCs'\n",
       "    obsp: 'connectivities', 'distances'"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bdata2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'MH9179824': 'Moderate',\n",
       " 'newcastle65': 'Healthy',\n",
       " 'MH9143327': 'Moderate',\n",
       " 'MH9143326': 'Moderate',\n",
       " 'MH9143325': 'Severe',\n",
       " 'MH9143320': 'Severe',\n",
       " 'MH9143276': 'Mild',\n",
       " 'MH9143274': 'Severe',\n",
       " 'MH8919226': 'Healthy',\n",
       " 'MH9179821': 'Moderate',\n",
       " 'MH8919230': 'LPS',\n",
       " 'MH9179823': 'Moderate',\n",
       " 'MH8919327': 'Severe',\n",
       " 'newcastle20': 'Mild',\n",
       " 'MH9143370': 'Non_covid',\n",
       " 'MH8919331': 'Moderate',\n",
       " 'MH9143420': 'Mild',\n",
       " 'MH9143373': 'Non_covid',\n",
       " 'MH9143372': 'Non_covid',\n",
       " 'MH9143371': 'Moderate',\n",
       " 'MH9143323': 'Moderate',\n",
       " 'MH9143277': 'Critical',\n",
       " 'MH8919277': 'LPS',\n",
       " 'MH9143270': 'Moderate',\n",
       " 'MH8919333': 'Healthy',\n",
       " 'MH8919332': 'Healthy',\n",
       " 'MH9179828': 'Non_covid',\n",
       " 'newcastle21': 'Critical',\n",
       " 'MH8919281': 'LPS',\n",
       " 'MH8919280': 'LPS',\n",
       " 'MH8919278': 'LPS',\n",
       " 'MH8919276': 'LPS',\n",
       " 'MH8919233': 'LPS',\n",
       " 'MH8919227': 'Healthy',\n",
       " 'MH8919326': 'Mild',\n",
       " 'newcastle59': 'Mild',\n",
       " 'MH9143422': 'Mild',\n",
       " 'MH9143271': 'Mild',\n",
       " 'newcastle49': 'Severe',\n",
       " 'MH9143324': 'Moderate',\n",
       " 'MH9143322': 'Moderate',\n",
       " 'MH8919330': 'Moderate',\n",
       " 'MH8919329': 'Critical',\n",
       " 'MH9179822': 'Severe',\n",
       " 'MH9143426': 'Mild',\n",
       " 'MH8919283': 'Healthy',\n",
       " 'MH9179827': 'Non_covid',\n",
       " 'MH8919178': 'Healthy',\n",
       " 'MH8919177': 'Healthy',\n",
       " 'MH8919232': 'LPS',\n",
       " 'MH8919229': 'LPS',\n",
       " 'MH8919176': 'Healthy',\n",
       " 'MH9143424': 'Moderate',\n",
       " 'MH8919228': 'LPS',\n",
       " 'MH8919179': 'Healthy',\n",
       " 'MH9143423': 'Moderate',\n",
       " 'MH9143421': 'Mild',\n",
       " 'MH8919231': 'LPS',\n",
       " 'MH9179826': 'Severe',\n",
       " 'MH9143321': 'Mild',\n",
       " 'MH9143427': 'Moderate',\n",
       " 'MH9143425': 'Asymptomatic',\n",
       " 'MH9143275': 'Moderate',\n",
       " 'MH8919279': 'LPS',\n",
       " 'MH9179825': 'Moderate',\n",
       " 'newcastle74': 'Healthy',\n",
       " 'MH9143272': 'Mild',\n",
       " 'MH8919282': 'Healthy',\n",
       " 'AP3': 'Moderate',\n",
       " 'AP8': 'Severe',\n",
       " 'AP12': 'Critical',\n",
       " 'AP10': 'Mild',\n",
       " 'AP1': 'Severe',\n",
       " 'AP2': 'Critical',\n",
       " 'AP4': 'Mild',\n",
       " 'AP9': 'Moderate',\n",
       " 'AP5': 'Severe',\n",
       " 'AP11': 'Mild',\n",
       " 'CV0904': 'Healthy',\n",
       " 'CV0902': 'Healthy',\n",
       " 'CV0025': 'Moderate',\n",
       " 'CV0144': 'Severe',\n",
       " 'CV0068': 'Asymptomatic',\n",
       " 'CV0059': 'Asymptomatic',\n",
       " 'CV0052': 'Moderate',\n",
       " 'CV0074': 'Mild',\n",
       " 'CV0200': 'Critical',\n",
       " 'CV0176': 'Severe',\n",
       " 'CV0084': 'Moderate',\n",
       " 'CV0100': 'Mild',\n",
       " 'CV0262': 'Asymptomatic',\n",
       " 'CV0164': 'Asymptomatic',\n",
       " 'CV0911': 'Healthy',\n",
       " 'CV0120': 'Moderate',\n",
       " 'CV0128': 'Mild',\n",
       " 'CV0929': 'Healthy',\n",
       " 'CV0284': 'Critical',\n",
       " 'CV0178': 'Severe',\n",
       " 'CV0037': 'Moderate',\n",
       " 'CV0234': 'Moderate',\n",
       " 'CV0137': 'Moderate',\n",
       " 'CV0104': 'Mild',\n",
       " 'CV0094': 'Moderate',\n",
       " 'CV0134': 'Moderate',\n",
       " 'CV0073': 'Mild',\n",
       " 'CV0155': 'Asymptomatic',\n",
       " 'CV0915': 'Healthy',\n",
       " 'CV0201': 'Critical',\n",
       " 'CV0917': 'Healthy',\n",
       " 'CV0062': 'Mild',\n",
       " 'CV0171': 'Asymptomatic',\n",
       " 'CV0279': 'Critical',\n",
       " 'CV0160': 'Mild',\n",
       " 'CV0231': 'Asymptomatic',\n",
       " 'CV0939': 'Healthy',\n",
       " 'CV0257': 'Moderate',\n",
       " 'CV0058': 'Asymptomatic',\n",
       " 'CV0050': 'Mild',\n",
       " 'CV0180': 'Critical',\n",
       " 'CV0926': 'Healthy',\n",
       " 'CV0326': 'Mild',\n",
       " 'CV0934': 'Healthy',\n",
       " 'CV0940': 'Healthy',\n",
       " 'CV0944': 'Healthy'}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "severity_dict = dict(zip([x.split('_')[1] if '_' in x else x for x in bdata2.obs['sample_id']], bdata2.obs['Status_on_day_collection_summary_v2']))\n",
    "severity_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# separate to the classification groups\n",
    "from collections import defaultdict\n",
    "sample_dict = defaultdict(list)\n",
    "statuses = ['Healthy', 'Asymptomatic', 'Mild', 'Moderate', 'Severe', 'Critical', 'Non_covid', 'LPS']\n",
    "for s in statuses:\n",
    "    for k in severity_dict.keys():\n",
    "        if severity_dict[k] == s:\n",
    "            sample_dict[s].append(k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dandelion class object with n_obs = 43917 and n_contigs = 89311\n",
       "    data: 'sequence_id', 'sequence', 'rev_comp', 'productive', 'v_call', 'd_call', 'j_call', 'sequence_alignment', 'germline_alignment', 'junction', 'junction_aa', 'v_cigar', 'd_cigar', 'j_cigar', 'stop_codon', 'vj_in_frame', 'locus', 'junction_length', 'np1_length', 'np2_length', 'v_sequence_start', 'v_sequence_end', 'v_germline_start', 'v_germline_end', 'd_sequence_start', 'd_sequence_end', 'd_germline_start', 'd_germline_end', 'j_sequence_start', 'j_sequence_end', 'j_germline_start', 'j_germline_end', 'v_score', 'v_identity', 'v_support', 'd_score', 'd_identity', 'd_support', 'j_score', 'j_identity', 'j_support', 'fwr1', 'fwr2', 'fwr3', 'fwr4', 'cdr1', 'cdr2', 'cdr3', 'cell_id', 'c_call', 'consensus_count', 'umi_count', 'v_call_10x', 'd_call_10x', 'j_call_10x', 'junction_10x', 'junction_10x_aa', 'v_call_genotyped', 'germline_alignment_d_mask', 'sample_id', 'c_sequence_alignment', 'c_germline_alignment', 'c_sequence_start', 'c_sequence_end', 'c_score', 'c_identity', 'c_support', 'c_call_10x', 'junction_aa_length', 'fwr1_aa', 'fwr2_aa', 'fwr3_aa', 'fwr4_aa', 'cdr1_aa', 'cdr2_aa', 'cdr3_aa', 'sequence_alignment_aa', 'v_sequence_alignment_aa', 'd_sequence_alignment_aa', 'j_sequence_alignment_aa', 'mu_freq', 'duplicate_count', 'clone_id'\n",
       "    metadata: 'sample_id', 'clone_id', 'clone_id_by_size', 'isotype', 'lightchain', 'status', 'vdj_status', 'productive', 'umi_counts_heavy', 'umi_counts_light', 'v_call_heavy', 'v_call_light', 'j_call_heavy', 'j_call_light', 'c_call_heavy', 'c_call_light'\n",
       "    distance: None\n",
       "    edges: None\n",
       "    layout: None\n",
       "    graph: None"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vdj = ddl.read_h5('dandelion_output/covid_jan_2021_bcells_vdj.h5')\n",
    "vdj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "vdjs = defaultdict(dict)\n",
    "\n",
    "for s in statuses:\n",
    "    vdjs[s] = ddl.Dandelion(vdj.data[vdj.data['sample_id'].isin(sample_dict[s])].copy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating distances... : 100%|██████████| 6/6 [02:50<00:00, 28.49s/it]\n",
      "Generating edge list : 100%|██████████| 57/57 [00:00<00:00, 1435.99it/s]\n",
      "Linking edges : 100%|██████████| 5292/5292 [00:22<00:00, 233.28it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "generating network layout\n",
      "CPU times: user 5min, sys: 2min 19s, total: 7min 19s\n",
      "Wall time: 5min 1s\n"
     ]
    }
   ],
   "source": [
    "%time ddl.tl.generate_network(vdjs['Healthy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating distances... : 100%|██████████| 5/5 [00:16<00:00,  3.37s/it]\n",
      "Generating edge list : 100%|██████████| 72/72 [00:00<00:00, 1370.41it/s]\n",
      "Linking edges : 100%|██████████| 1685/1685 [00:01<00:00, 1042.00it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "generating network layout\n",
      "CPU times: user 31.3 s, sys: 514 ms, total: 31.8 s\n",
      "Wall time: 33.3 s\n"
     ]
    }
   ],
   "source": [
    "%time ddl.tl.generate_network(vdjs['Asymptomatic'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating distances... : 100%|██████████| 5/5 [10:38<00:00, 127.80s/it]\n",
      "Generating edge list : 100%|██████████| 160/160 [00:00<00:00, 1286.98it/s]\n",
      "Linking edges : 100%|██████████| 10396/10396 [01:53<00:00, 91.77it/s] \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "generating network layout\n",
      "CPU times: user 17min 51s, sys: 3min 29s, total: 21min 20s\n",
      "Wall time: 19min 4s\n"
     ]
    }
   ],
   "source": [
    "%time ddl.tl.generate_network(vdjs['Mild'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating distances... : 100%|██████████| 7/7 [14:44<00:00, 126.29s/it]\n",
      "Generating edge list : 100%|██████████| 400/400 [00:00<00:00, 1153.47it/s]\n",
      "Linking edges : 100%|██████████| 10302/10302 [02:17<00:00, 74.83it/s]  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "generating network layout\n",
      "CPU times: user 23min 16s, sys: 3min 52s, total: 27min 9s\n",
      "Wall time: 24min 53s\n"
     ]
    }
   ],
   "source": [
    "%time ddl.tl.generate_network(vdjs['Moderate'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating distances... : 100%|██████████| 5/5 [04:38<00:00, 55.61s/it] \n",
      "Generating edge list : 100%|██████████| 219/219 [00:00<00:00, 1222.04it/s]\n",
      "Linking edges : 100%|██████████| 6468/6468 [00:35<00:00, 182.34it/s]  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "generating network layout\n",
      "CPU times: user 7min 53s, sys: 2min 34s, total: 10min 27s\n",
      "Wall time: 8min 11s\n"
     ]
    }
   ],
   "source": [
    "%time ddl.tl.generate_network(vdjs['Severe'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating distances... : 100%|██████████| 4/4 [01:20<00:00, 20.01s/it]\n",
      "Generating edge list : 100%|██████████| 76/76 [00:00<00:00, 1232.57it/s]\n",
      "Linking edges : 100%|██████████| 3864/3864 [00:10<00:00, 357.23it/s] \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "generating network layout\n",
      "CPU times: user 2min 19s, sys: 12.4 s, total: 2min 32s\n",
      "Wall time: 2min 33s\n"
     ]
    }
   ],
   "source": [
    "%time ddl.tl.generate_network(vdjs['Critical'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating distances... : 100%|██████████| 3/3 [00:04<00:00,  1.34s/it]\n",
      "Generating edge list : 100%|██████████| 13/13 [00:00<00:00, 899.53it/s]\n",
      "Linking edges : 100%|██████████| 670/670 [00:00<00:00, 920.51it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "generating network layout\n",
      "CPU times: user 9.94 s, sys: 115 ms, total: 10.1 s\n",
      "Wall time: 11.5 s\n"
     ]
    }
   ],
   "source": [
    "%time ddl.tl.generate_network(vdjs['Non_covid'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating distances... : 100%|██████████| 5/5 [00:24<00:00,  4.96s/it]\n",
      "Generating edge list : 100%|██████████| 31/31 [00:00<00:00, 1235.32it/s]\n",
      "Linking edges : 100%|██████████| 2101/2101 [00:02<00:00, 731.69it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "generating network layout\n",
      "CPU times: user 44 s, sys: 3.64 s, total: 47.6 s\n",
      "Wall time: 49.1 s\n"
     ]
    }
   ],
   "source": [
    "%time ddl.tl.generate_network(vdjs['LPS'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "for v in vdjs:\n",
    "    ddl.tl.clone_centrality(vdjs[v])\n",
    "    ddl.tl.clone_degree(vdjs[v])\n",
    "    ddl.tl.clone_size(vdjs[v])\n",
    "    ddl.tl.clone_size(vdjs[v], max_size=3)\n",
    "    ddl.update_metadata(vdjs[v], retrieve='mu_freq')\n",
    "    vdjs[v].metadata['Age'] = pd.Series(bdata2.obs['Age'])\n",
    "    vdjs[v].metadata['Sex'] = pd.Series(bdata2.obs['Sex'])\n",
    "    vdjs[v].metadata['Site'] = pd.Series(bdata2.obs['Site'])\n",
    "    vdjs[v].metadata['celltype_B'] = pd.Series(bdata2.obs['celltype_B'])\n",
    "    vdjs[v].metadata['celltype_B_v2'] = pd.Series(bdata2.obs['celltype_B_v2'])\n",
    "    vdjs[v].metadata['Collection_Day'] = pd.Series(bdata2.obs['Collection_Day'])\n",
    "    vdjs[v].metadata['patient_id'] = pd.Series(bdata2.obs['patient_id'])\n",
    "    vdjs[v].metadata['Status_on_day_collection_summary'] = pd.Series(bdata2.obs['Status_on_day_collection_summary'])\n",
    "    vdjs[v].metadata['Status_on_day_collection_summary_v2'] = pd.Series(bdata2.obs['Status_on_day_collection_summary_v2'])\n",
    "    vdjs[v].metadata['Worst_Clinical_Status'] = pd.Series(bdata2.obs['Worst_Clinical_Status'])\n",
    "    vdjs[v].metadata['Outcome'] = pd.Series(bdata2.obs['Outcome'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 34 s, sys: 650 ms, total: 34.7 s\n",
      "Wall time: 34.7 s\n"
     ]
    }
   ],
   "source": [
    "%time vdjs['Healthy'].write_h5('dandelion_output/covid_jan_2021_bcells_vdj_Healthy.h5', compression = 'blosc:lz4')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.25 s, sys: 60.6 ms, total: 2.31 s\n",
      "Wall time: 2.37 s\n"
     ]
    }
   ],
   "source": [
    "%time vdjs['Asymptomatic'].write_h5('dandelion_output/covid_jan_2021_bcells_vdj_Asymptomatic.h5', compression = 'blosc:lz4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4min 13s, sys: 2.27 s, total: 4min 16s\n",
      "Wall time: 4min 16s\n"
     ]
    }
   ],
   "source": [
    "%time vdjs['Mild'].write_h5('dandelion_output/covid_jan_2021_bcells_vdj_Mild.h5', compression = 'blosc:lz4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 5min 25s, sys: 2.96 s, total: 5min 28s\n",
      "Wall time: 5min 28s\n"
     ]
    }
   ],
   "source": [
    "%time vdjs['Moderate'].write_h5('dandelion_output/covid_jan_2021_bcells_vdj_Moderate.h5', compression = 'blosc:lz4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 16s, sys: 984 ms, total: 1min 17s\n",
      "Wall time: 1min 18s\n"
     ]
    }
   ],
   "source": [
    "%time vdjs['Severe'].write_h5('dandelion_output/covid_jan_2021_bcells_vdj_Severe.h5', compression = 'blosc:lz4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 15.5 s, sys: 348 ms, total: 15.8 s\n",
      "Wall time: 15.8 s\n"
     ]
    }
   ],
   "source": [
    "%time vdjs['Critical'].write_h5('dandelion_output/covid_jan_2021_bcells_vdj_Critical.h5', compression = 'blosc:lz4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.15 s, sys: 42.5 ms, total: 1.2 s\n",
      "Wall time: 1.2 s\n"
     ]
    }
   ],
   "source": [
    "%time vdjs['Non_covid'].write_h5('dandelion_output/covid_jan_2021_bcells_vdj_Non_covid.h5', compression = 'blosc:lz4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.4 s, sys: 67.3 ms, total: 3.46 s\n",
      "Wall time: 3.46 s\n"
     ]
    }
   ],
   "source": [
    "%time vdjs['LPS'].write_h5('dandelion_output/covid_jan_2021_bcells_vdj_LPS.h5', compression = 'blosc:lz4')"
   ]
  },
  {
   "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": 4
}
