12_CFScrossdataset.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nre-writing CFS cross dataset for re-do\\nmodified from: allbyallCFS.py and debugCFSallbyall_twentytwo.ipynb\\nmajor changes:\\n\\nShaina Lu\\nZador + Gillis Labs\\nMay 2020\\n'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "re-writing CFS cross dataset for re-do\n",
    "modified from: allenadultmouseISH/allbyallCFS.py, allenadultmouseISH/debugCFSallbyall_twentytwo.ipynb, and spatial/10_crossdataset.py\n",
    "major changes: MWU vectorized, getting and testing feat sets using pd.apply\n",
    "\n",
    "Shaina Lu\n",
    "Zador + Gillis Labs\n",
    "May 2020\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy import stats\n",
    "from sklearn import metrics\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import train_test_split #stratify train/test split\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import logging\n",
    "import timeit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in files and pre-processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#file paths\n",
    "ALLEN_FILT_PATH = \"/home/slu/spatial/data/ABAISH_filt_v6_avgdup.h5\"\n",
    "ONTOLOGY_PATH = \"/data/slu/allen_adult_mouse_ISH/ontologyABA.csv\"\n",
    "ST_CANTIN_FILT_PATH = \"/home/slu/spatial/data/cantin_ST_filt_v2.h5\"\n",
    "\n",
    "#outfiles\n",
    "FEATSETS_OUT = \"STtoABA_featsetsSTtrain_f1_CFS_052120.csv\"\n",
    "MOD_TEST_OUT = \"STtoABA_STtest_f1_CFS_052120.csv\"\n",
    "MOD_TRAIN_OUT = \"STtoABA_STtrain_f1_CFS_052120.csv\"\n",
    "CROSS_ALL_OUT = \"STtoABA_ABAall_f1_CFS_052120.csv\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read functions copied from 10_crossdataset.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def read_ABAdata():\n",
    "    \"\"\"read in all ABA datasets needed using pandas\"\"\"\n",
    "    metabrain = pd.read_hdf(ALLEN_FILT_PATH, key='metabrain', mode='r')\n",
    "    voxbrain = pd.read_hdf(ALLEN_FILT_PATH, key='avgvoxbrain', mode='r')\n",
    "    propontvox = pd.read_hdf(ALLEN_FILT_PATH, key='propontology', mode='r')\n",
    "    #geneIDName = pd.read_hdf(ALLEN_FILT_PATH, key='geneIDName', mode='r')\t\n",
    "\n",
    "    return metabrain, voxbrain, propontvox\n",
    "\n",
    "def read_STdata():\n",
    "    \"\"\"read in all ST datasets needed using pandas\"\"\"\n",
    "    STspotsmeta = pd.read_hdf(ST_CANTIN_FILT_PATH, key='STspotsmeta', mode='r')\n",
    "    STspots = pd.read_hdf(ST_CANTIN_FILT_PATH, key='STspots', mode='r')\n",
    "    STpropont = pd.read_hdf(ST_CANTIN_FILT_PATH, key='STpropont', mode='r')\n",
    "    \n",
    "    return STspotsmeta, STspots, STpropont\n",
    "\n",
    "def read_ontology():\n",
    "    ontology = pd.read_csv(ONTOLOGY_PATH)\n",
    "    ontology = ontology.drop([ontology.columns[5], ontology.columns[6]], axis=1)\n",
    "    ontology = ontology.fillna(-1)  #make root's parent -1\n",
    "\n",
    "    return ontology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pre-processing functions copied from 10_crossdataset.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def filterproponto(sampleonto):\n",
    "    \"\"\"pre-processing for propogated ontology\"\"\"\n",
    "    #remove brain areas that don't have any samples\n",
    "    sampleonto_sums = sampleonto.apply(lambda col: col.sum(), axis=0)\n",
    "    sampleonto = sampleonto.loc[:,sampleonto_sums > 5] #greater than 5 becuase less is not enough for train/test split to have non-zero areas\n",
    "    \n",
    "    return sampleonto\n",
    "\n",
    "def getleaves(propontvox, ontology):\n",
    "    \"\"\"helper function to get only leaf brain areas\"\"\"\n",
    "    #leaves are brain areas in the ontology that never show up in the parent column\n",
    "    allareas = list(propontvox)\n",
    "    parents = list(ontology.parent)\n",
    "    for i in range(len(parents)): #convert parents from float to int, ids are ints\n",
    "        parents[i] = int(parents[i])\n",
    "    \n",
    "    #remove parents from all areas\n",
    "    leaves = []\n",
    "    for area in allareas:\n",
    "        if int(area) not in parents:\n",
    "            leaves.append(area)\n",
    "    \n",
    "    print(\"number of leaf areas: %d\" %len(leaves))\n",
    "    return leaves\n",
    "\n",
    "def findoverlapareas(STonto, propontvox, ontology):\n",
    "    \"\"\"find leaf brain areas overlapping between the two datasets\"\"\"\n",
    "    leafST = getleaves(STonto, ontology)\n",
    "    leafABA = getleaves(propontvox, ontology)\n",
    "\n",
    "    leafboth = [] \n",
    "    for i in range(len(leafABA)):\n",
    "        if leafABA[i] in leafST:\n",
    "            leafboth.append(leafABA[i])\n",
    "    \n",
    "    STonto = STonto.loc[:,leafboth]\n",
    "    propontvox = propontvox.loc[:,leafboth]\n",
    "    \n",
    "    return STonto, propontvox    \n",
    "\n",
    "def zscore(voxbrain):\n",
    "    \"\"\"zscore voxbrain or subsets of voxbrain (rows: voxels, cols: genes)\"\"\"\n",
    "    #z-score on whole data set before splitting into test and train\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    scaler.fit(voxbrain)\n",
    "    z_voxbrain = scaler.transform(voxbrain)\n",
    "    \n",
    "    #store z-scored voxbrain as pandas dataframe\n",
    "    z_voxbrain = pd.DataFrame(z_voxbrain)\n",
    "    z_voxbrain.columns = voxbrain.columns\n",
    "    z_voxbrain.index = voxbrain.index\n",
    "    \n",
    "    return z_voxbrain\n",
    "\n",
    "def analytical_auroc(featurevector, binarylabels):\n",
    "    \"\"\"analytical calculation of auroc\n",
    "       inputs: feature (mean rank of expression level), binary label (ctxnotctx)\n",
    "       returns: auroc\n",
    "    \"\"\"\n",
    "    #sort ctxnotctx binary labels by mean rank, aescending\n",
    "    s = sorted(zip(featurevector, binarylabels))\n",
    "    feature_sort, binarylabels_sort = map(list, zip(*s))\n",
    "\n",
    "    #get the sum of the ranks in feature vector corresponding to 1's in binary vector\n",
    "    sumranks = 0\n",
    "    for i in range(len(binarylabels_sort)):\n",
    "        if binarylabels_sort[i] == 1:\n",
    "            sumranks = sumranks + feature_sort[i]\n",
    "    \n",
    "    poslabels = binarylabels.sum()\n",
    "    neglabels = (len(binarylabels) - poslabels)\n",
    "    \n",
    "    auroc = ((sumranks/(neglabels*poslabels)) - ((poslabels+1)/(2*neglabels)))\n",
    "    \n",
    "    return auroc\n",
    "\n",
    "def getoverlapgenes(STspots, ABAvox):\n",
    "    ABAgenes = list(ABAvox)\n",
    "    STgenes = list(STspots)\n",
    "    \n",
    "    #get overlapping genes\n",
    "    overlap = []\n",
    "    for i in range(len(ABAgenes)):\n",
    "        if ABAgenes[i] in STgenes:\n",
    "            overlap.append(ABAgenes[i])\n",
    "    \n",
    "    print(\"number of overlapping genes: %d\" %len(overlap))\n",
    "    \n",
    "    #index datasets to keep only genes that are overlapping\n",
    "    STspots = STspots.loc[:,overlap]\n",
    "    ABAvox = ABAvox.loc[:,overlap]\n",
    "    \n",
    "    return STspots, ABAvox\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### testing vectorized MWU against scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.09 s ± 6.99 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "#def getDEgenes2(Xtrain, ytrain):  \n",
    "#modified vectorized MWU from Ben and scipy.stats.mannwhitneyu source code\n",
    "Xtrain_ranked = Xtrain.apply(lambda col: sp.stats.mstats.rankdata(col), axis=0)\n",
    "n1 = ytrain.sum() #instances of brain area marked as 1\n",
    "n2 = len(ytrain) - n1\n",
    "U = Xtrain_ranked.loc[ytrain==1, :].sum() - ((n1*(n1+1))/2)\n",
    "\n",
    "T = Xtrain_ranked.apply(lambda col: sp.stats.tiecorrect(col), axis=0)\n",
    "sd = np.sqrt(T * n1 * n2 * (n1+n2+1) / 12.0)\n",
    "meanrank = n1*n2/2.0 + 0.5\n",
    "z = (U - meanrank) / sd\n",
    "#p = distributions.norm.sf(z)\n",
    "\n",
    "#Z = np.abs(U - ((n1*n2)/2)) / np.sqrt(n1*n2*(n1+n2+1)/12)\n",
    "pvals = pd.Series(stats.norm.sf(z), index=list(Xtrain))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11.5 s ± 2.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "cols = list(Xtrain)\n",
    "    \n",
    "#one-sided Mann-Whitney, Ha: areaofinterest < not areaofinterest\n",
    "u2 = []\n",
    "pvals2 = []\n",
    "genes2 = []\n",
    "errors2 = []\n",
    "for i in range(len(cols)):\n",
    "    try:\n",
    "        curr_u, curr_pval = sp.stats.mannwhitneyu(Xtrain.loc[ytrain ==1,\\\n",
    "                            cols[i]],Xtrain.loc[ytrain == 0,cols[i]],alternative='greater')\n",
    "        u2.append(curr_u)\n",
    "        pvals2.append(curr_pval)\n",
    "        genes2.append(cols[i])\n",
    "    except:   #some genes raise the error that \"all numbers are identical in mwu\"\n",
    "        u2.append(1)\n",
    "        pvals2.append(1)\n",
    "        errors2.append(cols[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11050\n",
      "(14299,)\n"
     ]
    }
   ],
   "source": [
    "#U stats are not equal because scipy implementation has a correction when all are tied\n",
    "print((U == u2).sum())\n",
    "print(U.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14255\n",
      "(14299,)\n"
     ]
    }
   ],
   "source": [
    "#p-vals are actually all equal despite the difference here due to differences in floats\n",
    "print((pvals == pvals2).sum())\n",
    "print(pvals.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.288289  |  original: 0.288289\n",
      "vector: 0.009699  |  original: 0.009699\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.009914  |  original: 0.009914\n",
      "vector: 0.419548  |  original: 0.419548\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.791791  |  original: 0.791791\n",
      "vector: 0.791791  |  original: 0.791791\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.745346  |  original: 0.745346\n",
      "vector: 0.008798  |  original: 0.008798\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.008798  |  original: 0.008798\n",
      "vector: 0.419548  |  original: 0.419548\n",
      "vector: 0.791791  |  original: 0.791791\n",
      "vector: 0.070494  |  original: 0.070494\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.953471  |  original: 0.953471\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.867956  |  original: 0.867956\n",
      "vector: 0.920523  |  original: 0.920523\n",
      "vector: 0.380344  |  original: 0.380344\n",
      "vector: 0.275606  |  original: 0.275606\n",
      "vector: 0.775916  |  original: 0.775916\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.180430  |  original: 0.180430\n",
      "vector: 0.071765  |  original: 0.071765\n",
      "vector: 0.419548  |  original: 0.419548\n"
     ]
    }
   ],
   "source": [
    "#few pvals that are different, are just due to slight differences in floats\n",
    "for i in range(len(np.where(pvals !=pvals2)[0])):\n",
    "    index = np.where(pvals!=pvals2)[0][i]\n",
    "    print(\"vector: %f  |  original: %f\" %(pvals.iloc[index], pvals2[index]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "show that flipping n1 and n2 and setting ytrain == 0 is truely alternative = less"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#my 'vectorized' version\n",
    "Xtrain_ranked = Xtrain.apply(lambda col: sp.stats.mstats.rankdata(col), axis=0)\n",
    "n2 = ytrain.sum() #instances of brain area marked as 1\n",
    "n1 = len(ytrain) - n1\n",
    "U = Xtrain_ranked.loc[ytrain==0, :].sum() - ((n1*(n1+1))/2)\n",
    "\n",
    "T = Xtrain_ranked.apply(lambda col: sp.stats.tiecorrect(col), axis=0)\n",
    "sd = np.sqrt(T * n1 * n2 * (n1+n2+1) / 12.0)\n",
    "meanrank = n1*n2/2.0 + 0.5\n",
    "z = (U - meanrank) / sd\n",
    "\n",
    "pvals = pd.Series(stats.norm.sf(z), index=list(Xtrain))\n",
    "\n",
    "#scipy versoin\n",
    "cols = list(Xtrain)\n",
    "    \n",
    "#one-sided Mann-Whitney, Ha: areaofinterest < not areaofinterest\n",
    "u2 = []\n",
    "pvals2 = []\n",
    "genes2 = []\n",
    "errors2 = []\n",
    "for i in range(len(cols)):\n",
    "    try:\n",
    "        curr_u, curr_pval = sp.stats.mannwhitneyu(Xtrain.loc[ytrain ==1,\\\n",
    "                            cols[i]],Xtrain.loc[ytrain == 0,cols[i]],alternative='less')\n",
    "        u2.append(curr_u)\n",
    "        pvals2.append(curr_pval)\n",
    "        genes2.append(cols[i])\n",
    "    except:   #some genes raise the error that \"all numbers are identical in mwu\"\n",
    "        u2.append(1)\n",
    "        pvals2.append(1)\n",
    "        errors2.append(cols[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14283\n",
      "(14299,)\n",
      "vector: 0.342334  |  original: 0.342334\n",
      "vector: 0.745346  |  original: 0.745346\n",
      "vector: 0.193018  |  original: 0.193018\n",
      "vector: 0.342334  |  original: 0.342334\n",
      "vector: 0.041490  |  original: 0.041490\n",
      "vector: 0.288289  |  original: 0.288289\n",
      "vector: 0.106370  |  original: 0.106370\n",
      "vector: 0.293989  |  original: 0.293989\n",
      "vector: 0.694708  |  original: 0.694708\n",
      "vector: 0.058077  |  original: 0.058077\n",
      "vector: 0.154992  |  original: 0.154992\n",
      "vector: 0.096759  |  original: 0.096759\n",
      "vector: 0.657666  |  original: 0.657666\n",
      "vector: 0.106370  |  original: 0.106370\n",
      "vector: 0.305292  |  original: 0.305292\n",
      "vector: 0.845008  |  original: 0.845008\n"
     ]
    }
   ],
   "source": [
    "print((pvals == pvals2).sum())\n",
    "print(pvals.shape)\n",
    "\n",
    "for i in range(len(np.where(pvals !=pvals2)[0])):\n",
    "    index = np.where(pvals!=pvals2)[0][i]\n",
    "    print(\"vector: %f  |  original: %f\" %(pvals.iloc[index], pvals2[index]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In conclusion, the vectorized version of MWU from deconstructing the function takes less than half the time of the non-vectorized version. Some steps are still a bit linear (tiecorrect, ranking, but pd. apply helps with that). \\\n",
    "\\\n",
    "The p-values between the my version and the scipy version are the same. My vectorized version implements alternative = \"greater\" for x > y. Just flip order of inputs to do alternative = \"less.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CFS functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modified from previous for speed-ups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcDE(Xtrain, ytrain):\n",
    "    #Ha: areaofinterest > not areaofinterest; i.e. alternative = greater\n",
    "    Xtrain_ranked = Xtrain.apply(lambda col: sp.stats.mstats.rankdata(col), axis=0)\n",
    "    n1 = ytrain.sum() #instances of brain area marked as 1\n",
    "    n2 = len(ytrain) - n1\n",
    "    U = Xtrain_ranked.loc[ytrain==1, :].sum() - ((n1*(n1+1))/2)\n",
    "\n",
    "    T = Xtrain_ranked.apply(lambda col: sp.stats.tiecorrect(col), axis=0)\n",
    "    sd = np.sqrt(T * n1 * n2 * (n1+n2+1) / 12.0)\n",
    "    meanrank = n1*n2/2.0 + 0.5\n",
    "    z = (U - meanrank) / sd\n",
    "\n",
    "    pvals_greater = pd.Series(stats.norm.sf(z), index=list(Xtrain), name='pvals_greater')\n",
    "    \n",
    "    #Ha: areaofinterest < notareaofinterest; i.e. alternative = less\n",
    "    Xtrain_ranked = Xtrain.apply(lambda col: sp.stats.mstats.rankdata(col), axis=0)\n",
    "    n2 = ytrain.sum() #instances of brain area marked as 1\n",
    "    n1 = len(ytrain) - n1\n",
    "    U = Xtrain_ranked.loc[ytrain==0, :].sum() - ((n1*(n1+1))/2)\n",
    "\n",
    "    T = Xtrain_ranked.apply(lambda col: sp.stats.tiecorrect(col), axis=0)\n",
    "    sd = np.sqrt(T * n1 * n2 * (n1+n2+1) / 12.0)\n",
    "    meanrank = n1*n2/2.0 + 0.5\n",
    "    z = (U - meanrank) / sd\n",
    "\n",
    "    pvals_less = pd.Series(stats.norm.sf(z), index=list(Xtrain), name='pvals_less')\n",
    "    \n",
    "    allpvals = pd.concat([pvals_greater, pvals_less], axis=1)\n",
    "    return allpvals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getDEgenes(allpvals, numtotal):\n",
    "    #melt\n",
    "    allpvals['gene'] = allpvals.index\n",
    "    allpvals_melt = allpvals.melt(id_vars='gene')\n",
    "    #sort by p-value\n",
    "    allpvals_melt = allpvals_melt.sort_values(by='value', ascending=True)\n",
    "    #get top X number of DE genes\n",
    "    topDEgenes = allpvals_melt.iloc[0:numtotal, :]\n",
    "    \n",
    "    return topDEgenes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_featset(featurecorrs, ranksdf, ylabels, seedgene):\n",
    "    \"\"\"picking feature set based on fwd sellection, random seed, and lowest possible corr,\n",
    "       stop when average auroc prediction is no longer improving\n",
    "       inputs: featurecorrs - correlation matrix of features being considered; seedgene - gene to start CFS\n",
    "       returns: feature set\"\"\"\n",
    "    #start with passed in randomly picked gene\n",
    "    featset = [seedgene]\n",
    "    #get start performance\n",
    "    curr_auroc = analytical_auroc(sp.stats.mstats.rankdata(ranksdf.loc[:,featset].mean(axis=1)), ylabels)\n",
    "    improving = True\n",
    "    while improving:\n",
    "        #look at all other possible features and take lowest correlated to seed, others in feat set\n",
    "        means = featurecorrs.loc[:,featset].mean(axis=1)  #get average corr across choosen features\n",
    "        featset.append(means.idxmin())                  #gets row name of min mean corrs, picks first of ties\n",
    "        #check featset performance\n",
    "        new_auroc = analytical_auroc(sp.stats.mstats.rankdata(ranksdf.loc[:,featset].mean(axis=1)),ylabels)\n",
    "        if new_auroc <= curr_auroc:  #if not improved, stop\n",
    "            featset.pop(len(featset)-1)\n",
    "            final_auroc = curr_auroc\n",
    "            improving = False\n",
    "        else:\n",
    "            curr_auroc = new_auroc\n",
    "            \n",
    "    return featset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [],
   "source": [
    "def applyCFS(zXtrain, zXtest, zXcross, ytrain, ytest, ycross):\n",
    "    #calculate DE for all genes across the two brain areas\n",
    "    allpvals = calcDE(zXtrain, ytrain)\n",
    "    #get top X DE genes\n",
    "    topDEgenes = getDEgenes(allpvals, 500)\n",
    "\n",
    "    #ranks DE genes\n",
    "    #train\n",
    "    rankedXtrain = zXtrain.loc[:, topDEgenes.gene]\n",
    "    rankedXtrain.loc[:,(topDEgenes.variable=='pvals_less').values] = \\\n",
    "                     -1 * rankedXtrain.loc[:,(topDEgenes.variable=='pvals_less').values]\n",
    "    rankedXtrain = rankedXtrain.apply(lambda col: sp.stats.mstats.rankdata(col), axis=0)\n",
    "    #test\n",
    "    rankedXtest = zXtest.loc[:, topDEgenes.gene]\n",
    "    rankedXtest.loc[:,(topDEgenes.variable=='pvals_less').values] = \\\n",
    "                    -1 * rankedXtest.loc[:,(topDEgenes.variable=='pvals_less').values]\n",
    "    rankedXtest = rankedXtest.apply(lambda col: sp.stats.mstats.rankdata(col), axis=0)\n",
    "    #cross\n",
    "    rankedXcross = zXcross.loc[:, topDEgenes.gene]\n",
    "    rankedXcross.loc[:,(topDEgenes.variable=='pvals_less').values] = \\\n",
    "                    -1 * rankedXcross.loc[:,(topDEgenes.variable=='pvals_less').values]\n",
    "    rankedXcross = rankedXcross.apply(lambda col: sp.stats.mstats.rankdata(col), axis=0)\n",
    "\n",
    "    #correlation matrix (spearman, b/c already ranked)\n",
    "    traincorrs = np.corrcoef(rankedXtrain.values.T)\n",
    "    traincorrs = pd.DataFrame(traincorrs, index=topDEgenes.gene.values, columns=topDEgenes.gene.values)\n",
    "\n",
    "    #get 100 feature sets using CFS\n",
    "    random.seed(0)\n",
    "    random.seed(42)\n",
    "    startingpts = pd.Series(random.sample(list(traincorrs),100))\n",
    "    featsets = startingpts.apply(lambda x: get_featset(traincorrs, rankedXtrain, ytrain, x))\n",
    "    trainaurocs = featsets.apply(lambda x: analytical_auroc(sp.stats.mstats.rankdata(rankedXtrain.loc[:,x].mean(axis=1)), ytrain))\n",
    "    testaurocs = featsets.apply(lambda x: analytical_auroc(sp.stats.mstats.rankdata(rankedXtest.loc[:,x].mean(axis=1)), ytest))\n",
    "    crossaurocs = featsets.apply(lambda x: analytical_auroc(sp.stats.mstats.rankdata(rankedXcross.loc[:,x].mean(axis=1)), ycross))\n",
    "\n",
    "    #return all 100 feature sets and aurocs\n",
    "    return featsets, trainaurocs, testaurocs, crossaurocs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mod_propont = STpropont\n",
    "#mod_data = STspots\n",
    "#cross_propont = ABApropont\n",
    "#cross_data = ABAvox\n",
    "def getallbyall(mod_data, mod_propont, cross_data, cross_propont):\n",
    "    #initialize zeros dataframe to store entries\n",
    "    allbyall_featsets = pd.DataFrame(index=list(mod_propont), columns=list(mod_propont))\n",
    "    allbyall_test = pd.DataFrame(index=list(mod_propont), columns=list(mod_propont))\n",
    "    allbyall_train = pd.DataFrame(index=list(mod_propont), columns=list(mod_propont))\n",
    "    allbyall_cross = pd.DataFrame(index=list(mod_propont), columns=list(mod_propont))\n",
    "\n",
    "    areas = list(mod_propont)\n",
    "    #for each column, brain area\n",
    "    for i in range(mod_propont.shape[1]):\n",
    "        logging.debug(\"starting col %d\" %i)\n",
    "        #for each row in each column\n",
    "        for j in range(i+1,mod_propont.shape[1]): #upper triangular!\n",
    "            #print(\"brain area j: %d\" %j)\n",
    "            area1 = areas[i]\n",
    "            area2 = areas[j]\n",
    "            #get binary label vectors\n",
    "            ylabels = mod_propont.loc[mod_propont[area1]+mod_propont[area2] != 0, area1]\n",
    "            ycross = cross_propont.loc[cross_propont[area1]+cross_propont[area2] != 0, area1]\n",
    "            #ylabels = pd.Series(np.random.permutation(ylabels1),index=ylabels1.index) #try permuting\n",
    "            #subset train and test sets for only samples in the two areas\n",
    "            Xcurr = mod_data.loc[mod_propont[area1]+mod_propont[area2] != 0, :]\n",
    "            Xcrosscurr = cross_data.loc[cross_propont[area1]+cross_propont[area2] != 0, :]\n",
    "            #split train test for X data and y labels\n",
    "            #split data function is seeded so all will split the same way\n",
    "            Xtrain, Xtest, ytrain, ytest = train_test_split(Xcurr, ylabels, test_size=0.5,\\\n",
    "                                                            random_state=42, shuffle=True,\\\n",
    "                                                            stratify=ylabels)\n",
    "            #z-score train and test folds\n",
    "            zXtrain = zscore(Xtrain)\n",
    "            zXtest = zscore(Xtest)\n",
    "            zXcross = zscore(Xcrosscurr)\n",
    "\n",
    "            featsets, currauroc_train, currauroc_test, currauroc_cross = applyCFS(zXtrain, zXtest, zXcross, ytrain, ytest, ycross)\n",
    "            allbyall_featsets.iloc[i,j] = featsets.values\n",
    "            allbyall_train.iloc[i,j] = currauroc_train.values\n",
    "            allbyall_test.iloc[i,j] = currauroc_test.values\n",
    "            allbyall_cross.iloc[i,j] = currauroc_cross.values\n",
    "\n",
    "            break\n",
    "\n",
    "        #periodically save\n",
    "        #if i%10 == 0:\n",
    "        #    logging.debug(\"saving\")\n",
    "        #    allbyall_featsets.to_csv(FEATSETS_OUT, sep=',', header=True, index=False)\n",
    "        #    allbyall_train.to_csv(ST_TRAIN_OUT, sep=',', header=True, index=False)\n",
    "        #    allbyall_test.to_csv(ST_TEST_OUT, sep=',', header=True, index=False)\n",
    "        #    allbyall_cross.to_csv(ABA_ALL_OUT, sep=',', header=True, index=False)\n",
    "\n",
    "        break\n",
    "\n",
    "\n",
    "    return allbyall_featsets, allbyall_train, allbyall_test, allbyall_cross"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Main"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read in data\n",
    "ontology = read_ontology()\n",
    "ABAmeta, ABAvox, ABApropont = read_ABAdata()\n",
    "STmeta, STspots, STpropont = read_STdata()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of leaf areas: 461\n",
      "number of leaf areas: 560\n"
     ]
    }
   ],
   "source": [
    "#filter brain areas for those that have at least x samples\n",
    "STpropont = filterproponto(STpropont)\n",
    "ABApropont = filterproponto(ABApropont)\n",
    "#filter brain areas for overlapping leaf areas\n",
    "STpropont, ABApropont = findoverlapareas(STpropont, ABApropont, ontology)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of overlapping genes: 14299\n"
     ]
    }
   ],
   "source": [
    "#keep only genes that are overlapping between the two datasets\n",
    "STspots, ABAvox = getoverlapgenes(STspots, ABAvox)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {},
   "outputs": [],
   "source": [
    "#predictability matrix using CFS\n",
    "allbyall_featsets, allbyall_train, allbyall_test, allbyall_cross = getallbyall(STspots, STpropont, ABAvox, ABApropont)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#write AUROC matrices to outfiles\n",
    "#allbyall_featsets.to_csv(FEATSETS_OUT, sep=',', header=True, index=False)\n",
    "#allbyall_train.to_csv(MOD_TRAIN_OUT, sep=',', header=True, index=False)\n",
    "#allbyall_test.to_csv(MOD_TEST_OUT, sep=',', header=True, index=False)\n",
    "#allbyall_cross.to_csv(CROSS_ALL_OUT, sep=',', header=True, index=False)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}