Cedric_gene_functions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[NbConvertApp] Converting notebook Cedric_gene_functions.ipynb to script\n",
      "[NbConvertApp] Writing 5955 bytes to Cedric_gene_functions.py\n"
     ]
    }
   ],
   "source": [
    "!jupyter nbconvert --to script Cedric_gene_functions.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## geneMNItoSchaefer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "def geneMNItoSchaefer (donor, numParcels, numNetworks, mm):\n",
    "\n",
    "    import itertools\n",
    "    import os\n",
    "    import subprocess\n",
    "    import numpy as np\n",
    "    import nibabel as nib\n",
    "    from nilearn import plotting\n",
    "    import pandas\n",
    "    import pickle\n",
    "    \n",
    "    print 'Start to map genes from donor %s to Schaefer at %s parcels and %s networks at %s mm resolution' \\\n",
    "        % (donor, numParcels, numNetworks, mm) + '\\n'\n",
    "    #files\n",
    "    wk_dir = '/Users/hxia/Desktop/BBL/' #change as appropriate\n",
    "    project_path = os.path.join(wk_dir,'data/joy/BBL/projects/prsConnectivity')\n",
    "    \n",
    "    schaefer_path = os.path.join(wk_dir,'data/joy/BBL/studies/pnc/template/Schaefer/MNI')\n",
    "    parcel_filename = os.path.join(schaefer_path, 'Schaefer2018_%sParcels_%sNetworks_order_FSLMNI152_%smm.nii.gz' \\\n",
    "                               % (numParcels, numNetworks, mm))\n",
    "    community_filename = os.path.join(schaefer_path,'Schaefer2018_%sParcels_%sNetworks_order.txt' \\\n",
    "                                     % (numParcels,numNetworks))\n",
    "    \n",
    "    # Genes\n",
    "    ABI_path = os.path.join(wk_dir,project_path,'ABI')\n",
    "    gene_filename = os.path.join(ABI_path,'normalized_microarray_donor%s' % donor,'SampleAnnot.csv')\n",
    "    \n",
    "    # Parcellation\n",
    "    img = nib.load(parcel_filename)\n",
    "    img_data = img.get_data()\n",
    "    plotting.plot_img(img,draw_cross=False,title='Schaefer %s Parcels' % numParcels,colorbar=True)\n",
    "\n",
    "    # Community Assignment\n",
    "    community = pandas.read_csv(community_filename,delim_whitespace=True,header=None)\n",
    "    \n",
    "    # Genes\n",
    "    gene = pandas.io.parsers.read_csv(gene_filename,delimiter=',')\n",
    "    gene_mni_coord = gene.iloc[:,-3:]\n",
    "    \n",
    "    # Parse Community file\n",
    "    community_assignment ={parcel[0]: parcel[1].split('_')[2] \n",
    "                       for idx, parcel in community.iterrows()}\n",
    "    community_assignment[0] = 'NotAssigned'\n",
    "    \n",
    "    # Initiate an empty dataframe of nrow = num gene, ncol = 3 for x y z coordinates\n",
    "    gene_vox_coord = pandas.DataFrame(np.zeros(gene_mni_coord.shape))\n",
    "    gene_vox_coord.columns = ['img_x','img_y','img_z']\n",
    "    \n",
    "    print 'Mapping genes to image coordinates' + '\\n'\n",
    "    # Loop through each gene coordinate in MNI space\n",
    "    # and then apply FSL std2imgcoord to transform the points\n",
    "    # to the image voxel coordiates\n",
    "\n",
    "    for ridx, row in gene_mni_coord.iterrows():\n",
    "        coord_cmd = \"echo %f %f %f|std2imgcoord -img %s -std %s -vox \" \\\n",
    "        % (row[0],row[1],row[2],parcel_filename,parcel_filename) # a shell command to apply FSL std2imgcoord\n",
    "        vox_coord_str = subprocess.check_output(coord_cmd, shell = True) # get the output from the shell\n",
    "        vox_coord_str_list = vox_coord_str.split()\n",
    "        for cidx, coord in enumerate(vox_coord_str_list): \n",
    "            vox_coord_float = float(coord)\n",
    "            gene_vox_coord.iloc[ridx][cidx] = vox_coord_float\n",
    "    gene_vox_coord_rd = np.round(gene_vox_coord).astype(int)\n",
    "    \n",
    "    gene_parcel_assignment ={idx: img_data[tuple(gene)].astype(int)\n",
    "                       for idx, gene in gene_vox_coord_rd.iterrows()}\n",
    "    \n",
    "    gene_community_assignment = {idx : community_assignment[gene_parcel]\n",
    "                            for idx, gene_parcel in gene_parcel_assignment.iteritems()}\n",
    "    \n",
    "    # Concat a master loopup table of gene assignment\n",
    "    gene_assignment = pandas.DataFrame([gene_parcel_assignment, gene_community_assignment]).T\n",
    "    gene_assignment.columns = ['d{}'.format(i) for i, col in enumerate(gene_assignment, 1)]\n",
    "    gene_assignment.columns = ['parcel','community']\n",
    "    gene_assignment = pandas.concat([gene_mni_coord,gene_vox_coord,gene_assignment],axis=1)\n",
    "    \n",
    "    print gene_assignment.head(10)\n",
    "    print '\\n'\n",
    "    #Save the variables\n",
    "    variables_filename = os.path.join(project_path,'ABI/gene_mapping','%sdonor_%sParcels_%sNetwork_%smm.pkl' \\\n",
    "                                      % (donor,numParcels,numNetworks,mm))\n",
    "    with open(variables_filename, 'w') as f:\n",
    "        pickle.dump([gene_mni_coord, gene_vox_coord, gene_vox_coord_rd,\\\n",
    "                     gene_parcel_assignment, gene_community_assignment,gene_assignment], f)\n",
    "    #write gene_assignment to a csv\n",
    "    csv_filename = os.path.join(project_path,'ABI/gene_mapping','%sdonor_%sParcels_%sNetwork_%smm.csv' \\\n",
    "                                % (donor,numParcels,numNetworks,mm))\n",
    "    gene_assignment.to_csv(path_or_buf=csv_filename)\n",
    "    \n",
    "    print 'Saved variables at ' + variables_filename + '\\n'\n",
    "    print 'Saved csv at ' + csv_filename + '\\n'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## selectProbe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def selectProbe (donor):\n",
    "    import os\n",
    "    import pandas as pd\n",
    "    \n",
    "    # Set up the working directory\n",
    "    # All paths after this are relative path to the root of chead\n",
    "    # this is where I have mounted chead to on my laptop. \n",
    "    # Change the above to an empty string if working directly on chead \n",
    "    wk_dir = '/Users/hxia/Desktop/BBL/'\n",
    "\n",
    "    # Set up the project directory\n",
    "    project_path = os.path.join(wk_dir,'data/joy/BBL/projects/prsConnectivity')\n",
    "    \n",
    "    # set paths\n",
    "    donor_path = os.path.join(project_path,'ABI','normalized_microarray_donor%s' % donor)\n",
    "    probes_info = pd.read_csv(os.path.join(donor_path,'Probes.csv'))\n",
    "    gene_expr = pd.read_csv(os.path.join(donor_path,'MicroarrayExpression.csv'),header=None)\n",
    "    sample_info = pd.read_csv(os.path.join(donor_path,'SampleAnnot.csv'))\n",
    "    \n",
    "    # get probes with available gene names\n",
    "    valid_probes = probes_info['entrez_id']>0\n",
    "    gene_expr = gene_expr[valid_probes]\n",
    "    probes_info = probes_info[valid_probes]\n",
    "    \n",
    "    # collapse rows\n",
    "    probe_names = probes_info['gene_symbol'][probes_info['probe_id'] == gene_expr[0]]\n",
    "    datCollapsed = gene_expr.groupby(probe_names).max()\n",
    "    \n",
    "    # create column names\n",
    "    gene_col_names = ['probe_id']\n",
    "    structure_ids = list(sample_info['structure_id'].astype(str))\n",
    "    gene_col_names.extend(structure_ids)\n",
    "    datCollapsed.columns = gene_col_names\n",
    "    \n",
    "    # output\n",
    "    return datCollapsed"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "myenv2",
   "language": "python",
   "name": "myenv2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}