Extract.probe_id.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from google.cloud import bigquery as bq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chrom</th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>sig-DMD</th>\n",
       "      <th>DMV</th>\n",
       "      <th>pairs</th>\n",
       "      <th>numberOfPairs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>862600</td>\n",
       "      <td>865000</td>\n",
       "      <td>sig-DMD_1</td>\n",
       "      <td>DMV_4</td>\n",
       "      <td>nontrip_trip,normal_trip</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr1</td>\n",
       "      <td>930000</td>\n",
       "      <td>934000</td>\n",
       "      <td>sig-DMD_2</td>\n",
       "      <td>DMV_6</td>\n",
       "      <td>nontrip_trip,normal_nontrip,normal_trip</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr1</td>\n",
       "      <td>1209000</td>\n",
       "      <td>1211000</td>\n",
       "      <td>sig-DMD_3</td>\n",
       "      <td>DMV_7</td>\n",
       "      <td>normal_nontrip,normal_trip</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>chr1</td>\n",
       "      <td>1286600</td>\n",
       "      <td>1296000</td>\n",
       "      <td>sig-DMD_4</td>\n",
       "      <td>DMV_9</td>\n",
       "      <td>ers_err,normal_nontrip,normal_trip</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chr1</td>\n",
       "      <td>1365000</td>\n",
       "      <td>1369600</td>\n",
       "      <td>sig-DMD_5</td>\n",
       "      <td>DMV_11</td>\n",
       "      <td>ers_err,normal_nontrip,normal_trip</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  chrom    start      end    sig-DMD     DMV  \\\n",
       "0  chr1   862600   865000  sig-DMD_1   DMV_4   \n",
       "1  chr1   930000   934000  sig-DMD_2   DMV_6   \n",
       "2  chr1  1209000  1211000  sig-DMD_3   DMV_7   \n",
       "3  chr1  1286600  1296000  sig-DMD_4   DMV_9   \n",
       "4  chr1  1365000  1369600  sig-DMD_5  DMV_11   \n",
       "\n",
       "                                     pairs  numberOfPairs  \n",
       "0                 nontrip_trip,normal_trip              2  \n",
       "1  nontrip_trip,normal_nontrip,normal_trip              3  \n",
       "2               normal_nontrip,normal_trip              2  \n",
       "3       ers_err,normal_nontrip,normal_trip              3  \n",
       "4       ers_err,normal_nontrip,normal_trip              3  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigDMD = pd.read_table(\"./preliminary/sig-DMD.union.includeERR.reference\", header=None)\n",
    "sigDMD.columns = ['chrom', 'start', 'end', 'sig-DMD', 'DMV', 'pairs', 'numberOfPairs']\n",
    "sigDMD.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4375, 7)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigDMD.shape #4375 signatures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "client = bq.Client()\n",
    "hm450_anno_ref = 'isb-cgc.platform_reference.methylation_annotation'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\\\n",
    "SELECT \n",
    "    CONCAT(\"chr\",CHR) AS chrom, \n",
    "    (MAPINFO-1) AS start, (MAPINFO+1) AS `end`, \n",
    "    IlmnID AS probe_id, UCSC\n",
    "FROM\n",
    "    `{}`\n",
    "WHERE\n",
    "    CHR IS NOT NULL\n",
    "    AND MAPINFO IS NOT NULL\n",
    "ORDER BY \n",
    "    CHR, MAPINFO\n",
    "ASC\n",
    "\"\"\".format(hm450_anno_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(485512, 5)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chrom</th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>probe_id</th>\n",
       "      <th>UCSC</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>15864</td>\n",
       "      <td>15866</td>\n",
       "      <td>cg13869341</td>\n",
       "      <td>[{'RefGene_Group': 'Body', 'RefGene_Accession'...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr1</td>\n",
       "      <td>18826</td>\n",
       "      <td>18828</td>\n",
       "      <td>cg14008030</td>\n",
       "      <td>[{'RefGene_Group': 'Body', 'RefGene_Accession'...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr1</td>\n",
       "      <td>29406</td>\n",
       "      <td>29408</td>\n",
       "      <td>cg12045430</td>\n",
       "      <td>[{'RefGene_Group': 'TSS200', 'RefGene_Accessio...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>chr1</td>\n",
       "      <td>29424</td>\n",
       "      <td>29426</td>\n",
       "      <td>cg20826792</td>\n",
       "      <td>[{'RefGene_Group': 'TSS200', 'RefGene_Accessio...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chr1</td>\n",
       "      <td>29434</td>\n",
       "      <td>29436</td>\n",
       "      <td>cg00381604</td>\n",
       "      <td>[{'RefGene_Group': 'TSS200', 'RefGene_Accessio...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  chrom  start    end    probe_id  \\\n",
       "0  chr1  15864  15866  cg13869341   \n",
       "1  chr1  18826  18828  cg14008030   \n",
       "2  chr1  29406  29408  cg12045430   \n",
       "3  chr1  29424  29426  cg20826792   \n",
       "4  chr1  29434  29436  cg00381604   \n",
       "\n",
       "                                                UCSC  \n",
       "0  [{'RefGene_Group': 'Body', 'RefGene_Accession'...  \n",
       "1  [{'RefGene_Group': 'Body', 'RefGene_Accession'...  \n",
       "2  [{'RefGene_Group': 'TSS200', 'RefGene_Accessio...  \n",
       "3  [{'RefGene_Group': 'TSS200', 'RefGene_Accessio...  \n",
       "4  [{'RefGene_Group': 'TSS200', 'RefGene_Accessio...  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hm450_anno = client.query(query).to_dataframe()\n",
    "\n",
    "print(hm450_anno.shape)\n",
    "\n",
    "hm450_anno.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cg13869341\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "str"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(hm450_anno.iloc[0,3])\n",
    "type(hm450_anno.iloc[0,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "str"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(hm450_anno.iloc[0,3][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "RefGene_Group=[] \n",
    "RefGene_Accession=[] \n",
    "RefGene_Name=[]\n",
    "\n",
    "def transform(l):\n",
    "    d=l[0]\n",
    "    global RefGene_Group, RefGene_Accession, RefGene_Name\n",
    "    RefGene_Group.append(d['RefGene_Group'])\n",
    "    RefGene_Accession.append(d['RefGene_Accession'])\n",
    "    RefGene_Name.append(d['RefGene_Name'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "485512\n",
      "485512\n",
      "485512\n"
     ]
    }
   ],
   "source": [
    "hm450_anno['UCSC'].apply(transform)\n",
    "\n",
    "print(len(RefGene_Group))\n",
    "print(len(RefGene_Accession))\n",
    "print(len(RefGene_Name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chrom</th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>probe_id</th>\n",
       "      <th>UCSC.RefGene_Group</th>\n",
       "      <th>UCSC.RefGene_Accession</th>\n",
       "      <th>UCSC.RefGene_Name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>15864</td>\n",
       "      <td>15866</td>\n",
       "      <td>cg13869341</td>\n",
       "      <td>Body</td>\n",
       "      <td>NR_024540</td>\n",
       "      <td>WASH5P</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr1</td>\n",
       "      <td>18826</td>\n",
       "      <td>18828</td>\n",
       "      <td>cg14008030</td>\n",
       "      <td>Body</td>\n",
       "      <td>NR_024540</td>\n",
       "      <td>WASH5P</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr1</td>\n",
       "      <td>29406</td>\n",
       "      <td>29408</td>\n",
       "      <td>cg12045430</td>\n",
       "      <td>TSS200</td>\n",
       "      <td>NR_024540</td>\n",
       "      <td>WASH5P</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>chr1</td>\n",
       "      <td>29424</td>\n",
       "      <td>29426</td>\n",
       "      <td>cg20826792</td>\n",
       "      <td>TSS200</td>\n",
       "      <td>NR_024540</td>\n",
       "      <td>WASH5P</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chr1</td>\n",
       "      <td>29434</td>\n",
       "      <td>29436</td>\n",
       "      <td>cg00381604</td>\n",
       "      <td>TSS200</td>\n",
       "      <td>NR_024540</td>\n",
       "      <td>WASH5P</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  chrom  start    end    probe_id UCSC.RefGene_Group UCSC.RefGene_Accession  \\\n",
       "0  chr1  15864  15866  cg13869341               Body              NR_024540   \n",
       "1  chr1  18826  18828  cg14008030               Body              NR_024540   \n",
       "2  chr1  29406  29408  cg12045430             TSS200              NR_024540   \n",
       "3  chr1  29424  29426  cg20826792             TSS200              NR_024540   \n",
       "4  chr1  29434  29436  cg00381604             TSS200              NR_024540   \n",
       "\n",
       "  UCSC.RefGene_Name  \n",
       "0            WASH5P  \n",
       "1            WASH5P  \n",
       "2            WASH5P  \n",
       "3            WASH5P  \n",
       "4            WASH5P  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hm450_anno.drop(columns=['UCSC'], inplace=True)\n",
    "hm450_anno['UCSC.RefGene_Group'] = RefGene_Group\n",
    "hm450_anno['UCSC.RefGene_Accession'] = RefGene_Accession\n",
    "hm450_anno['UCSC.RefGene_Name'] = RefGene_Name\n",
    "\n",
    "hm450_anno.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm450_anno.to_csv(\"./metadata/hm450_annotation\", sep=\"\\t\", header=None, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash \n",
    "sort -k1,1 -k2,2n datasets/hm450_annotation > /tmp/temp.bed\n",
    "\n",
    "mv /tmp/temp.bed datasets/hm450_annotation\n",
    "\n",
    "bedtools intersect -wa -wb \\\n",
    "-a preliminary/sig-DMD.union.includeERR.reference \\\n",
    "-b metadata/hm450_annotation | \\\n",
    "awk '{OFS=\"\\t\"} {print $8,$9,$10,$11,$12,$13,$14,$4,$5,$6,$7}' > \\\n",
    "metadata/sig-DMD.union.includeERR.hm450probes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25329 datasets/sig-DMD.union.includeERR.hm450probes\r\n"
     ]
    }
   ],
   "source": [
    "! wc -l metadata/sig-DMD.union.includeERR.hm450probes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr1\t864878\t864880\tcg02896266\tBody\tNM_152486\tSAMD11\tsig-DMD_1\tDMV_4\tnontrip_trip,normal_trip\t2\r\n",
      "chr1\t931326\t931328\tcg03648020\t\t\t\tsig-DMD_2\tDMV_6\tnontrip_trip,normal_nontrip,normal_trip\t3\r\n",
      "chr1\t933305\t933307\tcg01729262\t\t\t\tsig-DMD_2\tDMV_6\tnontrip_trip,normal_nontrip,normal_trip\t3\r\n",
      "chr1\t933387\t933389\tcg15882305\t\t\t\tsig-DMD_2\tDMV_6\tnontrip_trip,normal_nontrip,normal_trip\t3\r\n",
      "chr1\t933684\t933686\tcg15713103\t\t\t\tsig-DMD_2\tDMV_6\tnontrip_trip,normal_nontrip,normal_trip\t3\r\n",
      "chr1\t933986\t933988\tcg26003967\t\t\t\tsig-DMD_2\tDMV_6\tnontrip_trip,normal_nontrip,normal_trip\t3\r\n",
      "chr1\t1209195\t1209197\tcg02581421\t5'UTR\tNM_194315\tUBE2J2\tsig-DMD_3\tDMV_7\tnormal_nontrip,normal_trip\t2\r\n",
      "chr1\t1209672\t1209674\tcg03322129\tTSS1500\tNM_194458\tUBE2J2\tsig-DMD_3\tDMV_7\tnormal_nontrip,normal_trip\t2\r\n",
      "chr1\t1209695\t1209697\tcg01794932\tTSS1500\tNM_194458\tUBE2J2\tsig-DMD_3\tDMV_7\tnormal_nontrip,normal_trip\t2\r\n",
      "chr1\t1209719\t1209721\tcg26163375\tTSS1500\tNM_194458\tUBE2J2\tsig-DMD_3\tDMV_7\tnormal_nontrip,normal_trip\t2\r\n"
     ]
    }
   ],
   "source": [
    "! head metadata/sig-DMD.union.includeERR.hm450probes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Among 4375 signatures, only 3054 regions have HM450k probes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3054\r\n"
     ]
    }
   ],
   "source": [
    "! cut -f8 metadata/sig-DMD.union.includeERR.hm450probes | uniq | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 144/149 ERR signatures have probes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "144\r\n"
     ]
    }
   ],
   "source": [
    "! grep \"ers_err\" metadata/sig-DMD.union.includeERR.hm450probes | cut -f8 | uniq | wc -l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "echo -e Class\"\\t\"No. signatures\"\\t\" No. signatures have HM450k probes\"\\t\"Percentage >\\\n",
    "metadata/sig-DMD.union.includeERR.hm450probes.count\n",
    "\n",
    "class=\"err trip nontrip normal\"\n",
    "\n",
    "for c in $class; do \\\n",
    "    if [ \"$c\" = \"err\" ]; then\n",
    "        exp=\"$c\"\n",
    "    elif [ \"$c\" = \"trip\" ]; then\n",
    "        exp=\"_trip[[:print:]]*_trip\"\n",
    "    else\n",
    "        exp=\"$c\"\"[[:print:]]*\"\"$c\"\n",
    "    fi;\n",
    "    countp=$(grep \"$exp\" metadata/sig-DMD.union.includeERR.hm450probes | \\\n",
    "    cut -f8 | uniq | wc -l | awk '{print $1}');\n",
    "    count=$(grep -c \"$exp\" preliminary/sig-DMD.union.includeERR.reference);\n",
    "    percent=$(echo $countp / $count | bc -l);\n",
    "    echo -e $c\"\\t\"$count\"\\t\"$countp\"\\t\"$percent >> metadata/sig-DMD.union.includeERR.hm450probes.count;\n",
    "done\n",
    "totalp=$(cut -f8 metadata/sig-DMD.union.includeERR.hm450probes | uniq | wc -l)\n",
    "total=$(wc -l preliminary/sig-DMD.union.includeERR.reference | awk '{print $1}')\n",
    "percent=$(echo $totalp / $total | bc -l)\n",
    "echo -e total\"\\t\"$total\"\\t\"$totalp\"\\t\"$percent >> metadata/sig-DMD.union.includeERR.hm450probes.count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Class\tNo. signatures\t No. signatures have HM450k probes\tPercentage\r\n",
      "err\t149\t144\t.96644295302013422818\r\n",
      "trip\t950\t891\t.93789473684210526315\r\n",
      "nontrip\t926\t897\t.96868250539956803455\r\n",
      "normal\t3286\t2043\t.62172854534388314059\r\n",
      "total\t4375\t3054\t.69805714285714285714\r\n"
     ]
    }
   ],
   "source": [
    "! cat metadata/sig-DMD.union.includeERR.hm450probes.count"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}