main_complete.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Written by Jeff Zhang"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scanpy as sc\n",
    "import scanpy.tools as tl \n",
    "import os\n",
    "\n",
    "sc.settings.verbosity = 3\n",
    "sc.settings.set_figure_params(dpi=80, facecolor='white')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def show_dup(data):\n",
    "    '''\n",
    "    Prints out number of duplicate genes and gene names in data file\n",
    "    '''\n",
    "    print(\"\\n----- Showing Duplicates -----\\n\")\n",
    "    dup_genes = data.var_names[data.var_names.duplicated()].tolist()\n",
    "    print(dup_genes)\n",
    "    print(len(dup_genes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_file(path,type=\"\"):\n",
    "    '''\n",
    "    Loads data files into Anndata format\n",
    "    For mtx data \n",
    "      - 3 files are required with the EXACT following filenames(\"matrix.mtx\",\"barcodes.tsv\",\"genes.tsv\") for the function to run. \n",
    "      - All 3 files must be in one folder which the \"path\" variable refers to.\n",
    "    For h5 data\n",
    "      - The \"path\" variable refers to the h5 file.\n",
    "    '''\n",
    "    if type == \"\":\n",
    "        filetype = input(\"What type of file are you loading? \\nThis program can only read h5 and mtx files \\nIf reading mtx files, Ensure the file_path provided is a folder that contains a 'barcodes.tsv', 'genes.tsv', and 'matrix.mtx' file\")\n",
    "    print(\"\\n----- Loading File -----\\n\")\n",
    "    if type == \"mtx\":\n",
    "        adata = sc.read_10x_mtx(path,var_names='gene_symbols')\n",
    "        show_dup(adata)\n",
    "        \n",
    "        return(adata)\n",
    "    elif type == \"h5\":\n",
    "        adata = sc.read_10x_h5(path)\n",
    "        show_dup(adata)\n",
    "        \n",
    "        return(adata)\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_filter(adata,gene_n=False,cell_n=False):\n",
    "    '''\n",
    "    Filters the anndata input and returns filtered data\n",
    "       - exclude cells with less than \"gene_n\" number of genes\n",
    "       - exclude genes which appear in less than \"cell_n\" number of cells\n",
    "       - calculate the percentage of mitocondrial genes per cell and add to the metadata\n",
    "    '''\n",
    "    print(\"\\n----- Filtering by minimum gene and cell count required -----\\n\")\n",
    "    if gene_n != False:\n",
    "        \n",
    "        sc.pp.filter_cells(adata, min_genes=gene_n)\n",
    "    elif cell_n != False:\n",
    "        \n",
    "        sc.pp.filter_genes(adata, min_cells=cell_n)\n",
    "    adata.var['mt'] = adata.var_names.str.startswith('MT-')  \n",
    "    sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)\n",
    "    return(adata)\n",
    "\n",
    "    \n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def max_filter(adata,gene_n=False,mito_n=False):\n",
    "    '''\n",
    "    Filters out cells above certain cut offs and returns filtered data\n",
    "        - filters out cells with more thant \"gene_n\" number of genes\n",
    "        - filters out cells with less than \"mito_n\"% mito gene umi count\n",
    "    '''\n",
    "    print(\"\\n----- Filtering by gene and mitochondria count cutoff -----\\n\")\n",
    "    if gene_n != False:\n",
    "        adata = adata[adata.obs.n_genes_by_counts < gene_n, :] \n",
    "    if mito_n != False:\n",
    "        adata = adata[adata.obs.pct_counts_mt < mito_n, :]\n",
    "    return adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hvg(adata, minM, maxM, minD):\n",
    "    '''\n",
    "    DEFUNCT\n",
    "    Differentiates and Isolates Highly Variable Genes from adata and returns filtered adata\n",
    "    \n",
    "    '''\n",
    "    print(\"\\n----- Differentiating Highly Variable Genes -----\\n\")\n",
    "    sc.pp.highly_variable_genes(adata, min_mean=minM, max_mean=maxM, min_disp = minD)\n",
    "    print(\"\\n----- Highly Variable Genes Count -----\\n\")\n",
    "    print(adata.var.highly_variable.value_counts())\n",
    "    print(\"\\n----- Printing chart of Highly/Non-Highly Variable Genes -----\\n\")\n",
    "    sc.pl.highly_variable_genes(adata)\n",
    "    print(\"\\n----- Isolating Highly Variable Genes -----\\n\")\n",
    "    adata.raw = adata\n",
    "    adata = adata[:, adata.var.highly_variable]\n",
    "    return adata\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_principal(adata,prin_comp):\n",
    "    '''\n",
    "    DEFUNCT\n",
    "    For every principal component in \"prin_comp\", this function plots it\n",
    "    '''\n",
    "    print(\"\\n----- Plotting Principal components -----\\n\")\n",
    "    for i in prin_comp:\n",
    "        try:\n",
    "            sc.pl.pca(adata,color=str(i))\n",
    "        except:\n",
    "            print(f\"an error occurred with {i}\")\n",
    "            pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def main(path,full_list,type,hvg_min_mean=0.0125,hvg_max_mean=3,hvg_min_D=0.5,top=20,gene_num=False,cell_num=False,max_gene=False,max_mito=False,variance_max=10,prin_gene=[],n_n=None,n_p=None):\n",
    "    '''\n",
    "    A reproducible function for the preprocessing of scRNA files followed by visualization for \n",
    "    genes of interest\n",
    "\n",
    "    -Explanation of variables tbc...\n",
    "    \n",
    "    '''\n",
    "    print(f\"one: {full_list}\")\n",
    "    adata = load_file(path,type)\n",
    "    adata.var_names_make_unique()\n",
    "    \n",
    "    #Plot \"top\" number of expressed genes\n",
    "    sc.pl.highest_expr_genes(adata,n_top=top)\n",
    "\n",
    "    adata = min_filter(adata,gene_num,cell_num)\n",
    "\n",
    "    #Plot violin and scatter plots for genes\n",
    "    sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'], jitter=0.4, multi_panel=True)\n",
    "    sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')\n",
    "    sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')\n",
    "\n",
    "    adata = max_filter(adata=adata,gene_n=max_gene,mito_n=max_mito)\n",
    "\n",
    "    #normalize and scale data\n",
    "    sc.pp.normalize_total(adata, target_sum=1e4)\n",
    "    sc.pp.log1p(adata)\n",
    "\n",
    "    #filter out highly variable genes and plot them\n",
    "    sc.pp.highly_variable_genes(adata, min_mean=hvg_min_mean, max_mean=hvg_max_mean, min_disp = hvg_min_D)\n",
    "    sc.pl.highly_variable_genes(adata)\n",
    "\n",
    "    #Include and filter out necessary data for plotting\n",
    "    adata.raw = adata\n",
    "    adata = adata[:, adata.var.highly_variable]\n",
    "\n",
    "    #Regress and scale data\n",
    "    sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])\n",
    "    sc.pp.scale(adata, max_value=10)\n",
    "\n",
    "    #Compute PCA Matrix\n",
    "    sc.tl.pca(adata, svd_solver='arpack')\n",
    "\n",
    "    #Plot PCA for principal genes\n",
    "    for i in prin_gene:\n",
    "        try:\n",
    "            sc.pl.pca(adata,color=str(i))\n",
    "        except:\n",
    "            print(f\"an error occurred with {i}\")\n",
    "            pass\n",
    "    \n",
    "    #Plot variance ratio\n",
    "    sc.pl.pca_variance_ratio(adata, log=True)\n",
    "\n",
    "    #Find Neighbours for data points\n",
    "    sc.pp.neighbors(adata, n_neighbors=n_n, n_pcs=n_p)\n",
    "\n",
    "    #Employ Leiden Algorithm on Data\n",
    "    sc.tl.leiden(adata)\n",
    "\n",
    "    # Mapping out the coarse-grained connectivity structures of complex manifolds using partition-based graph abstraction (PAGA)\n",
    "    tl.paga(adata) \n",
    "    sc.pl.paga(adata, plot=False)\n",
    "    tl.umap(adata, init_pos='paga')\n",
    "\n",
    "    #Plot UMAP of principal Gene\n",
    "    for i in prin_gene:\n",
    "        try:\n",
    "            sc.pl.umap(adata,color=str(i))\n",
    "        except:\n",
    "            print(f\"an error occurred with {i}\")\n",
    "            pass\n",
    "    \n",
    "    full_list_2 = full_list.copy()\n",
    "    full_list_2.append(\"leiden\")\n",
    "\n",
    "    #Plot UMAP and violin chart of Genes in question to visualize relative expression\n",
    "    sc.pl.umap(adata, color=full_list_2, wspace=0.4, color_map=\"Blues\")\n",
    "    sc.pl.violin(adata, full_list, groupby='leiden')\n",
    "    return adata\n",
    "    \n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3 Examples of data processing visualization from datasets in the following github repository:\n",
    "https://github.com/Jeff-ZYX/Scanpy_scRNA_Visualization.git\n",
    "\n",
    "The Github repository aims to replicate R processing of data in the following paper with Scanpy:\n",
    "https://link.springer.com/content/pdf/10.1007/s11684-020-0754-0.pdf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Lung data\n",
    "file_type =\"h5\"\n",
    "file_path = \".\\data\\Lung\\GSE122960_RAW\\GSM3489182_Donor_01_filtered_gene_bc_matrices_h5.h5\"\n",
    "gene_min = 200\n",
    "cell_min = 3\n",
    "gene_cutoff = 3000\n",
    "mito_cutoff = 10\n",
    "hvg_minMean=0.01,\n",
    "hvg_maxMean=8,\n",
    "hvg_minD=0.5\n",
    "principal_components = ['SFTPB']\n",
    "num_n = 20\n",
    "num_p = 30\n",
    "interest = ['SFTPB','SFTPC','ACE2']\n",
    "\n",
    "\n",
    "adata = main(path=file_path,type=file_type,gene_num=gene_min,cell_num=cell_min,max_gene=gene_cutoff,max_mito=mito_cutoff,prin_gene=principal_components,n_n=num_n,n_p=num_p,full_list=interest,hvg_min_mean=hvg_minMean,hvg_max_mean=hvg_maxMean,hvg_min_D=hvg_minD)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #Bladder Data\n",
    "file_type = \"mtx\"\n",
    "file_path = \".\\data\\Bladder\\collated_data\"\n",
    "gene_min = 200\n",
    "cell_min = 3\n",
    "gene_cutoff = 2500\n",
    "mito_cutoff = 5\n",
    "principal_components = ['ACE2']\n",
    "num_n = 10\n",
    "num_p = 40\n",
    "interest = ['SPINK1','CLDN4','ACE2']\n",
    "\n",
    "main(path=file_path,type=file_type,gene_num=gene_min,cell_num=cell_min,max_gene=gene_cutoff,max_mito=mito_cutoff,prin_gene=principal_components,n_n=num_n,n_p=num_p,full_list=interest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Ileum data\n",
    "file_type = \"mtx\"\n",
    "file_path = \".\\data\\Ileum\\collated_data\"\n",
    "gene_min = 200\n",
    "cell_min = 3\n",
    "gene_cutoff = 2500\n",
    "mito_cutoff = 5\n",
    "principal_components = ['CST3']\n",
    "num_n = 20\n",
    "num_p = 30\n",
    "interest = ['FABP6','ACE2','ANPEP']\n",
    "\n",
    "main(path=file_path,type=file_type,gene_num=gene_min,cell_num=cell_min,max_gene=gene_cutoff,max_mito=mito_cutoff,prin_gene=principal_components,n_n=num_n,n_p=num_p,full_list=interest)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.0 ('scan': venv)",
   "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.0"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "f7468b6969dd139b9bea76206cd6d2c132a2bef8d09a55ff75dacff243a86631"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}