{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fname='Homo_sapiens.fa.gz'\n",
    "inDir='/cellar/users/btsui/Data/ensembl/snp_masked/'+fname\n",
    "import os\n",
    "withMicrobeDir=inDir.replace('.fa.gz','.microbe.fa.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r-- 1 btsui users 25 Dec 30 12:30 /cellar/users/btsui/Data/ensembl/snp_masked/Homo_sapiens.fa.gz\r\n"
     ]
    }
   ],
   "source": [
    "!ls -lah /cellar/users/btsui/Data/ensembl/snp_masked/Homo_sapiens.fa.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "specieName='Homo_sapiens'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/cellar/users/btsui/Data/ensembl/snp_masked/Homo_sapiens.microbe.fa.gz'"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "withMicrobeDir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### concat viral signatures from Methplan and 16s genes from mothur\n",
    "\n",
    "https://www.mothur.org/wiki/RDP_reference_files\n",
    "\n",
    "## bacteria and microbes reference mothur: \n",
    "1. trainset16_022016.rdp.tax : id to sequence tax mapping\n",
    "2. trainset16_022016.rdp.fasta.gz\n",
    "\n",
    "## viral \n",
    "1. ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz\n",
    "2. ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.2.1.genomic.fna.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!gunzip -c ../Microbiome/viral.1.1_2.1.genomic.fna.gz > ../Microbiome/viral.1.1_2.1.genomic.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!head ../Microbiome/viral.1.1_2.1.genomic.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/cellar/users/btsui/Project/METAMAP/notebook/RapMapTest/XGS_WGS/../Microbiome/viral.1.1_2.1.genomic.fa.fai\r\n"
     ]
    }
   ],
   "source": [
    "!echo $PWD/../Microbiome/viral.1.1_2.1.genomic.fa.fai"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "pd.read_csv('../Microbiome/viral.1.1_2.1.genomic.fa.fai',sep='\\t',header=None)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!samtools faidx ../Microbiome/viral.1.1_2.1.genomic.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!gunzip -c ../Microbiome/viral.1.1.genomic.fna.gz | head "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!echo $PWD/../Microbiome/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!gzip ../Microbiome/trainset16_022016.rdp.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!gunzip -c ../Microbiome/viral.1.1_2.1.genomic.fna.gz | head -n 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!cat ../Microbiome/trainset16_022016.rdp.fasta.gz ../Microbiome/viral.1.1_2.1.genomic.fna.gz > microbe.fa.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r-- 1 btsui users 86M Dec 30 12:51 /cellar/users/btsui/Data/ensembl/snp_masked/Homo_sapiens.fa.gz\r\n"
     ]
    }
   ],
   "source": [
    "!ls -lah /cellar/users/btsui/Data/ensembl/snp_masked/Homo_sapiens.fa.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.system('cat microbe.fa.gz '+inDir+' >'+withMicrobeDir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!ls -alh /cellar/users/btsui/Data/ensembl/snp_masked/Homo_sapiens.microbe.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r-- 1 btsui users 86M Dec 30 12:51 /cellar/users/btsui/Data/ensembl/snp_masked/Homo_sapiens.fa.gz\r\n"
     ]
    }
   ],
   "source": [
    "!ls -lah /cellar/users/btsui/Data/ensembl/snp_masked/Homo_sapiens.fa.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#withMicrobeDir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gunzipedFaDir=withMicrobeDir.replace('.fa.gz','.fa')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.system('rm '+gunzipedFaDir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.system('gunzip '+withMicrobeDir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.system('samtools faidx '+gunzipedFaDir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#import pandas as pd\n",
    "tmpDf=pd.read_csv(faDir+'.fai',sep='\\t',header=None)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"m1=tmpDf[0].str.contains('\\\\|')\\nm2=tmpDf[0].str.contains('_')\\ntmpDf.loc[m1,'state']='rdp'\\ntmpDf.loc[m2,'state']='ncbi_virus'\\ntmpDf.loc[~(m1|m2),'state']='human'\\ntmpDf\""
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"m1=tmpDf[0].str.contains('\\|')\n",
    "m2=tmpDf[0].str.contains('_')\n",
    "tmpDf.loc[m1,'state']='rdp'\n",
    "tmpDf.loc[m2,'state']='ncbi_virus'\n",
    "tmpDf.loc[~(m1|m2),'state']='human'\n",
    "tmpDf\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "myDir='/cellar/users/btsui/Data/BOWTIE_GENOME_SNP_INDEX/'+specieName+'/'\n",
    "faDir=gunzipedFaDir\n",
    "#tmpDir='/tmp/btsui/'+fname\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#os.system('cp '+faDir+' '+tmpDir)\n",
    "#os.system('gunzip '+tmpDir)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "cmd= 'bowtie2-build --threads 64 '+faDir+' '+myDir\n",
    "os.system('rm -r '+myDir)\n",
    "os.system('mkdir '+myDir)\n",
    "cmd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#os.system('mkdir '+myDir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tmpDf3=pd.read_csv('/cellar/users/btsui/per_fa_record_stat.txt',sep='\\t',header=None).sort_values(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>15804</th>\n",
       "      <td>NC_001357.1</td>\n",
       "      <td>7857</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 0     1  2  3\n",
       "15804  NC_001357.1  7857  0  0"
      ]
     },
     "execution_count": 171,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmpDf3[tmpDf3[0].str.contains('NC_001357')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>14162</th>\n",
       "      <td>NC_001493.2</td>\n",
       "      <td>134226</td>\n",
       "      <td>33</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14812</th>\n",
       "      <td>NC_024382.1</td>\n",
       "      <td>137090</td>\n",
       "      <td>33</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18228</th>\n",
       "      <td>NC_009823.1</td>\n",
       "      <td>9711</td>\n",
       "      <td>35</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17966</th>\n",
       "      <td>NC_004102.1</td>\n",
       "      <td>9646</td>\n",
       "      <td>39</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15950</th>\n",
       "      <td>NC_006641.1</td>\n",
       "      <td>15959</td>\n",
       "      <td>42</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17914</th>\n",
       "      <td>NC_030200.1</td>\n",
       "      <td>137448</td>\n",
       "      <td>43</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20658</th>\n",
       "      <td>NC_021858.1</td>\n",
       "      <td>1908524</td>\n",
       "      <td>59</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22235</th>\n",
       "      <td>NC_009127.1</td>\n",
       "      <td>295146</td>\n",
       "      <td>59</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13939</th>\n",
       "      <td>NC_020231.1</td>\n",
       "      <td>233501</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22138</th>\n",
       "      <td>NC_008912.1</td>\n",
       "      <td>3141</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14731</th>\n",
       "      <td>NC_020474.2</td>\n",
       "      <td>180421</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13843</th>\n",
       "      <td>NC_019491.1</td>\n",
       "      <td>291144</td>\n",
       "      <td>66</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20983</th>\n",
       "      <td>NC_005261.2</td>\n",
       "      <td>137821</td>\n",
       "      <td>68</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14853</th>\n",
       "      <td>NC_024709.1</td>\n",
       "      <td>33452</td>\n",
       "      <td>69</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15780</th>\n",
       "      <td>NC_001499.1</td>\n",
       "      <td>5894</td>\n",
       "      <td>77</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14004</th>\n",
       "      <td>NC_021312.1</td>\n",
       "      <td>459984</td>\n",
       "      <td>78</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15691</th>\n",
       "      <td>NC_002794.1</td>\n",
       "      <td>195859</td>\n",
       "      <td>102</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17235</th>\n",
       "      <td>NC_028834.1</td>\n",
       "      <td>48216</td>\n",
       "      <td>115</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20659</th>\n",
       "      <td>NC_022098.1</td>\n",
       "      <td>2473870</td>\n",
       "      <td>126</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14829</th>\n",
       "      <td>NC_024697.1</td>\n",
       "      <td>370920</td>\n",
       "      <td>358</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14211</th>\n",
       "      <td>NC_022518.1</td>\n",
       "      <td>9472</td>\n",
       "      <td>483</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15776</th>\n",
       "      <td>NC_001506.1</td>\n",
       "      <td>3811</td>\n",
       "      <td>506</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21974</th>\n",
       "      <td>NC_008168.1</td>\n",
       "      <td>104710</td>\n",
       "      <td>1535</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14240</th>\n",
       "      <td>NC_018464.1</td>\n",
       "      <td>927</td>\n",
       "      <td>3123</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18783</th>\n",
       "      <td>NC_032111.1</td>\n",
       "      <td>163005</td>\n",
       "      <td>8688</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22813</th>\n",
       "      <td>Y</td>\n",
       "      <td>57227415</td>\n",
       "      <td>16392</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22807</th>\n",
       "      <td>18</td>\n",
       "      <td>80373285</td>\n",
       "      <td>72541</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22810</th>\n",
       "      <td>21</td>\n",
       "      <td>46709983</td>\n",
       "      <td>79744</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22802</th>\n",
       "      <td>13</td>\n",
       "      <td>114364328</td>\n",
       "      <td>101329</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22811</th>\n",
       "      <td>22</td>\n",
       "      <td>50818468</td>\n",
       "      <td>152059</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 0          1       2  3\n",
       "14162  NC_001493.2     134226      33  0\n",
       "14812  NC_024382.1     137090      33  0\n",
       "18228  NC_009823.1       9711      35  0\n",
       "17966  NC_004102.1       9646      39  0\n",
       "15950  NC_006641.1      15959      42  0\n",
       "17914  NC_030200.1     137448      43  0\n",
       "20658  NC_021858.1    1908524      59  0\n",
       "22235  NC_009127.1     295146      59  0\n",
       "13939  NC_020231.1     233501      61  0\n",
       "22138  NC_008912.1       3141      61  0\n",
       "14731  NC_020474.2     180421      61  0\n",
       "13843  NC_019491.1     291144      66  0\n",
       "20983  NC_005261.2     137821      68  0\n",
       "14853  NC_024709.1      33452      69  0\n",
       "15780  NC_001499.1       5894      77  0\n",
       "14004  NC_021312.1     459984      78  0\n",
       "15691  NC_002794.1     195859     102  0\n",
       "17235  NC_028834.1      48216     115  0\n",
       "20659  NC_022098.1    2473870     126  0\n",
       "14829  NC_024697.1     370920     358  0\n",
       "14211  NC_022518.1       9472     483  0\n",
       "15776  NC_001506.1       3811     506  0\n",
       "21974  NC_008168.1     104710    1535  0\n",
       "14240  NC_018464.1        927    3123  0\n",
       "18783  NC_032111.1     163005    8688  0\n",
       "22813            Y   57227415   16392  0\n",
       "22807           18   80373285   72541  0\n",
       "22810           21   46709983   79744  0\n",
       "22802           13  114364328  101329  0\n",
       "22811           22   50818468  152059  0"
      ]
     },
     "execution_count": 177,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmpDf3[tmpDf3[2]>5].iloc[-50:-20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}