PermutationTest_Expression_kstest.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatial Expression of One Gene\n",
    "Tests if the spatial expression of a gene is significantly different from the spatial distribution of random genes. The p-vale and effect size are calculated using a Kolmogorov–Smirnov test.<br>\n",
    "7/29/19"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#packages\n",
    "import argparse \n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import csv\n",
    "import scipy.stats\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Block creates parser to interpret command line info and make arguments into variables.\n",
    "parser = argparse.ArgumentParser(description = \"handle inputs from SOMETHING script to run permutation tests\" +\n",
    "                                 \"on puck data\")\n",
    "parser.add_argument(\"-ern\", type=int,\n",
    "                   help = \"Enforced Read Number: enter 1 or 0. if 1, enforces that all random samples have the \"+\n",
    "                    \"same number of positive as the test sample at the cost of some computational time.\")\n",
    "parser.add_argument(\"-fg\", type=int,\n",
    "                   help = \"Filter Genes:input values to determine how genes are filtered: 0, in which case all \"+\n",
    "                    \"genes are analyzed (not recommended, due to false positives); 1, in which case genes are \"+\n",
    "                    \"filtered by within-dropseq-cluster expression; 2, in which case it's filtered by within-dropseq-\"+\n",
    "                    \"cluster variance; or 3, in which case genes either match the expression cutoff or variance \"+\n",
    "                    \"cutoff, and are labeled according to which they pass (or both)\")\n",
    "parser.add_argument(\"-pg\", type=int,\n",
    "                   help = \"Plot Genes: enter 1 or 0. If 1, will output a pdf with the significant genes at the \"+\n",
    "                    \"0.005 level plotted\")\n",
    "parser.add_argument(\"--ns\", type=int,\n",
    "                   help = \"Number of Samples: Number of samples for the null model. 1000 by default. Note: Runtime\"+\n",
    "                    \"scales linearly with numsamples\")\n",
    "parser.add_argument(\"--bc\", type=int,\n",
    "                   help = \"Bead Cutoff: minimum number of beads needed to assess significance of a gene. 15 by \"+\n",
    "                    \"default\") \n",
    "parser.add_argument(\"-bmf\", type=str,\n",
    "                   help = \"Bead Mapping File: enter data path for puck gene expression/cluster data file. If no \"+ \n",
    "                   \"extension, assumes csv\")\n",
    "parser.add_argument(\"-pn\", type=str,\n",
    "                   help = \"Puck Number: enter puck identifier here\")\n",
    "parser.add_argument('--clust', nargs='*', type=int,\n",
    "                   help=\"enter the cluster numbers you wish to analyze. Multiple arguments allowed\")\n",
    "\n",
    "print(parser)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In the block below, you may change the following for your data:\n",
    " -ern    : to force sample distribution to have the same number of beads as the data <br>\n",
    " -fg     : how genes are filtered <br>\n",
    " -pg     : if you want to plot significant genes<br>\n",
    " --ns    : the number of samples for the null model<br>\n",
    " --bc    : minimum number of beads expressing the gene to assess it<br>\n",
    " -bmf    : data path for puck data<br>\n",
    " -pn     : puck number<br>\n",
    " --clust : cluster(s) to analyze"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Block utilizes parser\n",
    "args = parser.parse_args('-ern 0 -fg 0 -pg 0 -bmf /broad/thechenlab/breanna/permutation_test_data -pn Puck_181206_3'.split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Block formalizes variables from parser info\n",
    "EnforcedReadNumbers = args.ern\n",
    "\n",
    "FilterGenes = args.fg\n",
    "\n",
    "PlotGenes = args.pg\n",
    "\n",
    "if args.ns is None:\n",
    "    NumSamples = 1000\n",
    "else: \n",
    "    NumSamples = args.ns\n",
    "\n",
    "if args.bc is None:\n",
    "    BeadCutoff = 15\n",
    "else: \n",
    "    BeadCutoff = args.bc\n",
    "    \n",
    "BeadMappingFile = args.bmf\n",
    "\n",
    "PuckNumber = args.pn\n",
    "\n",
    "if args.clust is None:\n",
    "    ClustertoAnalyze=[]\n",
    "else:\n",
    "    ClustertoAnalyze=args.clust\n",
    "\n",
    "DataPath = \"{}/{}.csv\".format(BeadMappingFile,PuckNumber) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Read in and save data\n",
    "AllMappedBeads=pd.read_csv(DataPath, header = 0, index_col = 0)\n",
    "if ClustertoAnalyze:\n",
    "    UniqueMappedBeads=AllMappedBeads[AllMappedBeads[\"cluster\"].isin(ClustertoAnalyze)]\n",
    "else: UniqueMappedBeads=AllMappedBeads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Count number of reads per bead\n",
    "genes_only=UniqueMappedBeads.iloc[:,0:-3]#exclude cluster number and coordinates\n",
    "NumReadsPerBead = genes_only.sum(axis=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Calculate pair-wise distances between each bead\n",
    "x=UniqueMappedBeads.as_matrix(columns=['xcoord'])\n",
    "BeadXCoordMatrix=x*np.ones((1,UniqueMappedBeads.shape[0]))\n",
    "y=UniqueMappedBeads.as_matrix(columns=['ycoord'])\n",
    "BeadYCoordMatrix=y*np.ones((1,UniqueMappedBeads.shape[0]))\n",
    "\n",
    "BeadPairwiseXValDifferences=BeadXCoordMatrix-np.transpose(BeadXCoordMatrix)\n",
    "BeadPairwiseYValDifferences=BeadYCoordMatrix-np.transpose(BeadYCoordMatrix)\n",
    "\n",
    "dist = lambda x,y: ((x**2)+(y**2))**(1/2)\n",
    "BeadPairwiseDistanceMat = dist(BeadPairwiseXValDifferences,BeadPairwiseYValDifferences)\n",
    "\n",
    "#set up bins for histograms later\n",
    "num_bins=100\n",
    "Triu = np.triu_indices(BeadPairwiseDistanceMat.shape[0],1)\n",
    "hist,bins=np.histogram(BeadPairwiseDistanceMat[Triu],num_bins)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Determines the probability of picking each bead for the null distribution based on the number of reads per bead\n",
    "NumReadsPerBead = UniqueMappedBeads.iloc[:,0:-3].sum(axis=1)\n",
    "ProbabilityPerBead=NumReadsPerBead/NumReadsPerBead.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Note: Filtering by variance has not yet been implemented"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#filtering genes according to \"FilterGenes\" input\n",
    "\n",
    "ExpressionGenes=[]\n",
    "VarianceGenes=[]\n",
    "\n",
    "for cluster in ClustertoAnalyze:\n",
    "    tmp=genes_only\n",
    "    #Filter by within-cluster expression\n",
    "    if FilterGenes==1 or FilterGenes==3:\n",
    "        GoodGenes=list(np.array(tmp.columns)[np.array(tmp.sum(axis=0)>0)])\n",
    "        ExpressionGenes=list(np.unique(ExpressionGenes+GoodGenes))\n",
    "    #Filter by within-cluster variance\n",
    "    #if FilterGenes==2 or FilterGenes==3:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Selecting your genes\n",
    "Here, the GeneNames variable is a list of genes that will be analyzed. Replace the selected kidney genes currently in the GeneNames list with your desired genes, or uncomment the lines below to analyze all genes in your sample.\n",
    "\n",
    "#### Note:\n",
    "If the genes you select are not in the sample, the method will fail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Perform permutation test and save data to file path provided to the arg parser.\n",
    "\n",
    "#Select kidney genes\n",
    "GeneNames=['Slc12a1', 'Umod', 'Ctgf', 'Nphs1','Nphs2','Wt1', 'Synpo', 'Itga8','Ptn',  'Plvap', 'Ehd3', 'Ren1','Slc27a2', 'Aqp2', 'Napsa', 'Aqp3', 'Pck1', 'Miox', 'Acsm2', 'Acsm3', 'Rarres2', 'Col4a2', 'C1qc', 'C1qa']\n",
    "\n",
    "#Uncomment the following two lines to analyze every gene in the sample\n",
    "#GeneNames=list(UniqueMappedBeads.columns)[:-3]\n",
    "#GeneNames=list(genes_only.columns)\n",
    "\n",
    "pvals=np.zeros(len(GeneNames))\n",
    "effectsize=np.zeros(len(GeneNames))\n",
    "\n",
    "for geneval in GeneNames:\n",
    "    if ClustertoAnalyze:\n",
    "        #filter genes based on filter genes\n",
    "        if FilterGenes==3:\n",
    "            PassingVariance=0\n",
    "            PassingExpression=0\n",
    "            if geneval in ExpressionGenes:\n",
    "                PassingExpression=1\n",
    "            if geneval in VarianceGenes:\n",
    "                PassingVariance=1\n",
    "            if not ExpressionGenes and not VarianceGenes:\n",
    "                pvals[GeneNames.index(geneval)]=-1\n",
    "                continue          \n",
    "        elif FilterGenes==1 and not geneval in ExpressionGenes:\n",
    "            pvals[GeneNames.index(geneval)]=-1\n",
    "            continue\n",
    "        elif FilterGenes==2 and not geneval in VarianceGenes:\n",
    "            pvals[GeneNames.index(geneval)]=-1\n",
    "            continue \n",
    "        \n",
    "    #filter out genes expressed in too few beads    \n",
    "    NumBeads=sum(UniqueMappedBeads[geneval]>0)\n",
    "    print('NumBeads')\n",
    "    print(NumBeads)\n",
    "    if NumBeads<BeadCutoff:\n",
    "        pvals[GeneNames.index(geneval)]=-1\n",
    "        print('here')\n",
    "        continue \n",
    "        \n",
    "    #This will give the true distribution\n",
    "    #Find which beads have gene expressed\n",
    "    A=np.matrix(UniqueMappedBeads[geneval]>0)\n",
    "    NonzeroBeads=BeadPairwiseDistanceMat[np.matmul(np.transpose(A),A)]\n",
    "    #create boolean matrix with locations of beads with genes expressed\n",
    "    NonzeroMat=np.reshape(NonzeroBeads,(int(NonzeroBeads.size**(1/2)),int(NonzeroBeads.size**(1/2))))\n",
    "    Triu = np.triu_indices(NonzeroMat.shape[0],1)\n",
    "    #Matrix with distances between beads with gene expressed \n",
    "    #Triu to take only one half of the symmetric matrix (other half is identical, avoid double-counting)\n",
    "    PairWiseDistances=NonzeroMat[Triu]\n",
    "    \n",
    "    #plot true distribution\n",
    "    n=plt.hist(PairWiseDistances, bins=num_bins, color='royalblue') \n",
    "    plt.rcParams.update({'font.size': 20})\n",
    "    # arguments are passed to np.histogram\n",
    "    plt.xlabel('Distance',fontsize=24)\n",
    "    plt.ylabel('Count',fontsize=24)\n",
    "    plt.title(\"TAL True Distance {}\".format(geneval),fontsize=29)\n",
    "    plt.show()\n",
    "    \n",
    "    DistanceDist=n[0]/n[0].sum()\n",
    "    \n",
    "    #Generate permuted distribution. There is a ton of duplication here, because this calculation is the \n",
    "    #same regardless of geneval. It only depends on the NUMBER of beads in which geneval appears.\n",
    "    AverageDistribution=np.zeros(num_bins)\n",
    "    RandomDists=[]\n",
    "    pvals_tmp=np.zeros(NumSamples)\n",
    "    effectsize_tmp=np.zeros(NumSamples)\n",
    "    for p in range(NumSamples):\n",
    "        if EnforcedReadNumbers:\n",
    "            NonzeroBeadsRandomTmp=np.random.choice(np.arange(len(NumReadsPerBead)),NumBeads,replace=False,p=ProbabilityPerBead)\n",
    "            NonzeroBeadsRandom=np.zeros(len(NumReadsPerBead))\n",
    "            for elem in NonzeroBeadsRandomTmp:\n",
    "                NonzeroBeadsRandom[elem]=1\n",
    "        else: NonzeroBeadsRandom=np.random.uniform(size=len(NumReadsPerBead))/NumBeads<ProbabilityPerBead\n",
    "        #Find which beads have gene expressed\n",
    "        A=np.matrix(NonzeroBeadsRandom)\n",
    "        NonzeroBeadsRandomTmp=BeadPairwiseDistanceMat[np.matmul(np.transpose(A),A)]\n",
    "        #create boolean matrix with locations of beads with genes expressed\n",
    "        NonzeroBeadsRandomMat=np.reshape(NonzeroBeadsRandomTmp,(int(NonzeroBeadsRandomTmp.size**(1/2)),int(NonzeroBeadsRandomTmp.size**(1/2))))\n",
    "        Triu = np.triu_indices(NonzeroBeadsRandomMat.shape[0],1)\n",
    "        #Matrix with distances between beads with gene expressed \n",
    "        #Triu to take only one half of the symmetric matrix (other half is identical, avoid double-counting)\n",
    "        RandomDistTmp=NonzeroBeadsRandomMat[Triu]\n",
    "        #save only the counts in each bin\n",
    "        RandomDists.append(np.histogram(RandomDistTmp,num_bins, range=(bins[0],bins[num_bins]))[0])\n",
    "        AverageDistribution+=np.array(RandomDists[p])\n",
    "        RandomDists[p]=np.divide(RandomDists[p],RandomDists[p].sum())\n",
    "        effectsize_tmp[p][0]=scipy.stats.ks_2samp(PairWiseDistances,RandomDistTmp)[0]\n",
    "        pvals_tmp[p][1]=scipy.stats.ks_2samp(PairWiseDistances,RandomDistTmp)[p][1]\n",
    "        \n",
    "           \n",
    "    pvals[GeneNames.index(geneval)]=pvals_tmp.max(axis=0)[1]\n",
    "    effectsize[GeneNames.index(geneval)]=effectsize_tmp.max(axis=0)[0]\n",
    "    \n",
    "    with open('kidney_expression_data_fixed.csv', 'a') as csvfile:\n",
    "        my_writer = csv.writer(csvfile, delimiter=',')\n",
    "        my_writer.writerow((geneval, pvals[GeneNames.index(geneval)], effectsize[GeneNames.index(geneval)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Selecting your genes\n",
    "Here, the GeneNames variable is a list of genes that will be analyzed. The lines that are currently commented out will run the method for every gene detected in the sample that meets the filter gene critera and is expressed in more beads than the bead cutoff number. \n",
    "\n",
    "Otherwise, replace the selected kidney genes currently in the GeneNames list with your desired genes. \n",
    "\n",
    "#### Note:\n",
    "If the genes you select are not in the sample, the method will fail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#p-values of genes analyzed\n",
    "pvals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#effect size of genes analyzed\n",
    "effectsize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Find genes that yield significant results\n",
    "SignificantGenes=list(np.array(GeneNames)[(pvals>0) & (pvals<.25)])\n",
    "SignificantGeneNames=[]\n",
    "for gene in SignificantGenes:\n",
    "    SignificantGeneNames.append(GeneNames[GeneNames.index(gene)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#genes that yielded significant results\n",
    "SignificantGeneNames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#function to plot significant genes (currently supports 19. For more, increase colors in color list)\n",
    "def plot_significant_genes(user_genes):\n",
    "    color_list=['b','darkorange','green','r','darkviolet','saddlebrown','magenta','orchid','k','olive','teal','cornflowerblue','gold','lawngreen','lightsalmon','plum','peru','y','turquoise']\n",
    "    c_num=0\n",
    "    empty_beads=pcounts_and_coords[pcounts_and_coords[user_gene].sum(axis=1)==0]\n",
    "    ax=empty_beads.plot(kind='scatter',x='xcoord',y='ycoord',color='lightgrey',alpha=.2,figsize=(10,10))\n",
    "    for gene in user_genes:\n",
    "        tmp=pcounts_and_coords[pcounts_and_coords[gene]>0]\n",
    "        ax = tmp.plot(kind=\"scatter\", x=\"xcoord\",y=\"ycoord\", color=color_list[c_num], label=gene,alpha=.5,ax=ax)\n",
    "        c_num+=1\n",
    "    plt.rcParams.update({'font.size': 20})\n",
    "    plt.legend(bbox_to_anchor=(1.04,1), loc=\"upper left\",fontsize='x-large')\n",
    "    ax.set_xlabel(\"x\")\n",
    "    ax.set_ylabel(\"y\")\n",
    "    ax.set_title(\"Location of Gene Expression\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plot_significant_genes(SignificantGenes)"
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}