02. Calculate triad category.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Relative expression levels of the A, B, and D genome homoeologs across triads\n",
    "This analysis focused exclusively on the gene triads which had a 1:1:1 correspondence across the three homoeologous genomes, including 17,400 syntenic and 1,074 non-syntenic triads (total of 18,474 triads or 55,422 genes). Again, starting from the subset of genes considered expressed using the initial 850 filter criterion, we defined a triad as expressed when the sum of the A, B, and D genome homoeologs was > 0.5 TPM. This allowed us to include triads in which, for example, only a single homoeolog was expressed, and which could later be classified as a dominant triad. Using this criteria, we defined a total of 53,259 genes (17,753 triads) which were considered expressed (Table S3).\n",
    "\n",
    "To standardize the relative expression of each homoeolog across the triad, we normalized the absolute TPM for each gene within the triad as follows:$$expression(A) = \\frac{TPM(A)}{TPM(A)+TPM(B)+TPM(D)}$$\n",
    "             $$expression(B) = \\frac{TPM(A)}{TPM(A)+TPM(B)+TPM(D)}$$\n",
    "             $$expression(D) = \\frac{TPM(A)}{TPM(A)+TPM(B)+TPM(D)}$$\n",
    "\n",
    "   Where A, B, and D represent the gene corresponding to the A, B, and D homoeologs in the triad. The normalized expression was calculated for each one of the intermediate tissues and for the average across all expressed tissues (“combined analysis” as described previously). Fig. S6 shows an example of these calculations for the roots and the combined analysis across three triads. The values of the relative contributions of each genome per triad were used to plot the ternary diagrams using the R package ggtern (57)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: gsubfn\n",
      "Loading required package: proto\n",
      "Loading required package: RSQLite\n",
      "Loading required package: spam\n",
      "Loading required package: grid\n",
      "Spam version 1.4-0 (2016-08-29) is loaded.\n",
      "Type 'help( Spam)' or 'demo( spam)' for a short introduction \n",
      "and overview of this package.\n",
      "Help for individual functions is also obtained by adding the\n",
      "suffix '.spam' to the function name, e.g. 'help( chol.spam)'.\n",
      "\n",
      "Attaching package: ‘spam’\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    backsolve, forwardsolve\n",
      "\n",
      "Loading required package: maps\n",
      "--\n",
      "Consider donating at: http://ggtern.com\n",
      "Even small amounts (say $10-50) are very much appreciated!\n",
      "Remember to cite, run citation(package = 'ggtern') for further info.\n",
      "--\n",
      "\n",
      "Attaching package: ‘ggtern’\n",
      "\n",
      "The following objects are masked from ‘package:gridExtra’:\n",
      "\n",
      "    arrangeGrob, grid.arrange\n",
      "\n",
      "The following objects are masked from ‘package:ggplot2’:\n",
      "\n",
      "    %+%, aes, annotate, calc_element, ggplot, ggplot_build,\n",
      "    ggplot_gtable, ggplotGrob, ggsave, layer_data, theme, theme_bw,\n",
      "    theme_classic, theme_dark, theme_gray, theme_light, theme_linedraw,\n",
      "    theme_minimal, theme_void\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(sqldf)\n",
    "library(ggplot2)\n",
    "library(reshape2)\n",
    "library(fields)\n",
    "library(\"gridExtra\")\n",
    "library(ggtern)\n",
    "library(clue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "850"
      ],
      "text/latex": [
       "850"
      ],
      "text/markdown": [
       "850"
      ],
      "text/plain": [
       "[1] 850"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "850"
      ],
      "text/latex": [
       "850"
      ],
      "text/markdown": [
       "850"
      ],
      "text/plain": [
       "[1] 850"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "metadata<-read.csv(\"./TablesForExploration/metadata.csv\")\n",
    "nrow(metadata)\n",
    "loadValuesFromExperiment<-function(metadata, folder, unit=\"tpm\", values=c(\"Development\")){\n",
    "    metadata$Sample.IDs <- gsub(\"-\",\".\",metadata$Sample.IDs)\n",
    "    \n",
    "    v<-values[1]\n",
    "    v<-gsub(\" \",\"_\",v)\n",
    "    v<-gsub(\",\",\".\",v)\n",
    "    path<-paste0(folder,\"/\",v,\"_\",unit,\".tsv\")\n",
    "    ret<-read.table(path, row.names = 1, header= TRUE)\n",
    "    if(length(values) > 1){\n",
    "      for(i in 2:length(values)){\n",
    "        v<-values[i]\n",
    "        v<-gsub(\" \",\"_\",v)\n",
    "        v<-gsub(\",\",\".\",v)\n",
    "        path<-paste0(folder,\"/\",v,\"_\",unit,\".tsv\")\n",
    "        tmp<-read.table(path, row.names = 1, header= TRUE)\n",
    "        ret<-cbind(ret,tmp)\n",
    "      }  \n",
    "    }\n",
    "    \n",
    "    md<-metadata[metadata$Sample.IDs%in%colnames(ret),]\n",
    "    ret<-ret[,as.character(md$Sample.IDs),]\n",
    "    list(ret,md)\n",
    "}\n",
    "folder<-\"./expressionValuesPerGene\"\n",
    "tpms  <-loadValuesFromExperiment(metadata, folder, unit=\"tpm\",  values=unique(metadata$study_title))\n",
    "\n",
    "\n",
    "metadata_used<-tpms[[2]]\n",
    "tpms<-tpms[[1]]\n",
    "\n",
    "nrow(metadata_used)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loading the homologies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>group_id</th><th scope=col>family_name</th><th scope=col>origin</th><th scope=col>chrs</th><th scope=col>cardinality_formal</th><th scope=col>cardinality_abs</th><th scope=col>A</th><th scope=col>B</th><th scope=col>D</th><th scope=col>synteny</th><th scope=col>TE</th><th scope=col>HC.LC</th><th scope=col>source</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>227                                                 </td><td>OG0021509                                           </td><td>ABD                                                 </td><td>7                                                   </td><td>1:1:1                                               </td><td>1:1:1                                               </td><td>TraesCS7A01G243100                                  </td><td>TraesCS7B01G148400                                  </td><td>TraesCS7D01G241900                                  </td><td>segmental homeologs                                 </td><td>NA                                                  </td><td>HC-only                                             </td><td>OrthoFinder SO,BLAST RBH,DAGchainer,i-ADHoRe,MCScanX</td></tr>\n",
       "\t<tr><td>253                                                 </td><td>OG0021212                                           </td><td>ABD                                                 </td><td>7                                                   </td><td>1:1:1                                               </td><td>1:1:1                                               </td><td>TraesCS7A01G360600                                  </td><td>TraesCS7B01G267100                                  </td><td>TraesCS7D01G362400                                  </td><td>segmental homeologs                                 </td><td>NA                                                  </td><td>HC-only                                             </td><td>OrthoFinder SO,BLAST RBH,i-ADHoRe,MCScanX           </td></tr>\n",
       "\t<tr><td>255                                                 </td><td>OG0021210                                           </td><td>ABD                                                 </td><td>7                                                   </td><td>1:1:1                                               </td><td>1:1:1                                               </td><td>TraesCS7A01G267100                                  </td><td>TraesCS7B01G165300                                  </td><td>TraesCS7D01G267800                                  </td><td>segmental homeologs                                 </td><td>NA                                                  </td><td>HC-only                                             </td><td>OrthoFinder SO,BLAST RBH,i-ADHoRe                   </td></tr>\n",
       "\t<tr><td>259                                                 </td><td>OG0021208                                           </td><td>ABD                                                 </td><td>7                                                   </td><td>1:1:1                                               </td><td>1:1:1                                               </td><td>TraesCS7A01G235400                                  </td><td>TraesCS7B01G133600                                  </td><td>TraesCS7D01G235500                                  </td><td>segmental homeologs                                 </td><td>NA                                                  </td><td>HC-only                                             </td><td>OrthoFinder SO,BLAST RBH,i-ADHoRe,MCScanX           </td></tr>\n",
       "\t<tr><td>263                                                 </td><td>OG0021206                                           </td><td>ABD                                                 </td><td>7                                                   </td><td>1:1:1                                               </td><td>1:1:1                                               </td><td>TraesCS7A01G185200                                  </td><td>TraesCS7B01G090200                                  </td><td>TraesCS7D01G186600                                  </td><td>segmental homeologs                                 </td><td>NA                                                  </td><td>HC-only                                             </td><td>OrthoFinder SO,BLAST RBH,DAGchainer                 </td></tr>\n",
       "\t<tr><td>269                                                 </td><td>OG0021202                                           </td><td>ABD                                                 </td><td>7                                                   </td><td>1:1:1                                               </td><td>1:1:1                                               </td><td>TraesCS7A01G155800                                  </td><td>TraesCS7B01G059900                                  </td><td>TraesCS7D01G156200                                  </td><td>segmental homeologs                                 </td><td>NA                                                  </td><td>HC-only                                             </td><td>OrthoFinder SO,BLAST RBH,DAGchainer,MCScanX         </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllllllllllll}\n",
       " group\\_id & family\\_name & origin & chrs & cardinality\\_formal & cardinality\\_abs & A & B & D & synteny & TE & HC.LC & source\\\\\n",
       "\\hline\n",
       "\t 227                                                  & OG0021509                                            & ABD                                                  & 7                                                    & 1:1:1                                                & 1:1:1                                                & TraesCS7A01G243100                                   & TraesCS7B01G148400                                   & TraesCS7D01G241900                                   & segmental homeologs                                  & NA                                                   & HC-only                                              & OrthoFinder SO,BLAST RBH,DAGchainer,i-ADHoRe,MCScanX\\\\\n",
       "\t 253                                                  & OG0021212                                            & ABD                                                  & 7                                                    & 1:1:1                                                & 1:1:1                                                & TraesCS7A01G360600                                   & TraesCS7B01G267100                                   & TraesCS7D01G362400                                   & segmental homeologs                                  & NA                                                   & HC-only                                              & OrthoFinder SO,BLAST RBH,i-ADHoRe,MCScanX           \\\\\n",
       "\t 255                                                  & OG0021210                                            & ABD                                                  & 7                                                    & 1:1:1                                                & 1:1:1                                                & TraesCS7A01G267100                                   & TraesCS7B01G165300                                   & TraesCS7D01G267800                                   & segmental homeologs                                  & NA                                                   & HC-only                                              & OrthoFinder SO,BLAST RBH,i-ADHoRe                   \\\\\n",
       "\t 259                                                  & OG0021208                                            & ABD                                                  & 7                                                    & 1:1:1                                                & 1:1:1                                                & TraesCS7A01G235400                                   & TraesCS7B01G133600                                   & TraesCS7D01G235500                                   & segmental homeologs                                  & NA                                                   & HC-only                                              & OrthoFinder SO,BLAST RBH,i-ADHoRe,MCScanX           \\\\\n",
       "\t 263                                                  & OG0021206                                            & ABD                                                  & 7                                                    & 1:1:1                                                & 1:1:1                                                & TraesCS7A01G185200                                   & TraesCS7B01G090200                                   & TraesCS7D01G186600                                   & segmental homeologs                                  & NA                                                   & HC-only                                              & OrthoFinder SO,BLAST RBH,DAGchainer                 \\\\\n",
       "\t 269                                                  & OG0021202                                            & ABD                                                  & 7                                                    & 1:1:1                                                & 1:1:1                                                & TraesCS7A01G155800                                   & TraesCS7B01G059900                                   & TraesCS7D01G156200                                   & segmental homeologs                                  & NA                                                   & HC-only                                              & OrthoFinder SO,BLAST RBH,DAGchainer,MCScanX         \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "group_id | family_name | origin | chrs | cardinality_formal | cardinality_abs | A | B | D | synteny | TE | HC.LC | source | \n",
       "|---|---|---|---|---|---|\n",
       "| 227                                                  | OG0021509                                            | ABD                                                  | 7                                                    | 1:1:1                                                | 1:1:1                                                | TraesCS7A01G243100                                   | TraesCS7B01G148400                                   | TraesCS7D01G241900                                   | segmental homeologs                                  | NA                                                   | HC-only                                              | OrthoFinder SO,BLAST RBH,DAGchainer,i-ADHoRe,MCScanX | \n",
       "| 253                                                  | OG0021212                                            | ABD                                                  | 7                                                    | 1:1:1                                                | 1:1:1                                                | TraesCS7A01G360600                                   | TraesCS7B01G267100                                   | TraesCS7D01G362400                                   | segmental homeologs                                  | NA                                                   | HC-only                                              | OrthoFinder SO,BLAST RBH,i-ADHoRe,MCScanX            | \n",
       "| 255                                                  | OG0021210                                            | ABD                                                  | 7                                                    | 1:1:1                                                | 1:1:1                                                | TraesCS7A01G267100                                   | TraesCS7B01G165300                                   | TraesCS7D01G267800                                   | segmental homeologs                                  | NA                                                   | HC-only                                              | OrthoFinder SO,BLAST RBH,i-ADHoRe                    | \n",
       "| 259                                                  | OG0021208                                            | ABD                                                  | 7                                                    | 1:1:1                                                | 1:1:1                                                | TraesCS7A01G235400                                   | TraesCS7B01G133600                                   | TraesCS7D01G235500                                   | segmental homeologs                                  | NA                                                   | HC-only                                              | OrthoFinder SO,BLAST RBH,i-ADHoRe,MCScanX            | \n",
       "| 263                                                  | OG0021206                                            | ABD                                                  | 7                                                    | 1:1:1                                                | 1:1:1                                                | TraesCS7A01G185200                                   | TraesCS7B01G090200                                   | TraesCS7D01G186600                                   | segmental homeologs                                  | NA                                                   | HC-only                                              | OrthoFinder SO,BLAST RBH,DAGchainer                  | \n",
       "| 269                                                  | OG0021202                                            | ABD                                                  | 7                                                    | 1:1:1                                                | 1:1:1                                                | TraesCS7A01G155800                                   | TraesCS7B01G059900                                   | TraesCS7D01G156200                                   | segmental homeologs                                  | NA                                                   | HC-only                                              | OrthoFinder SO,BLAST RBH,DAGchainer,MCScanX          | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  group_id family_name origin chrs cardinality_formal cardinality_abs\n",
       "1 227      OG0021509   ABD    7    1:1:1              1:1:1          \n",
       "2 253      OG0021212   ABD    7    1:1:1              1:1:1          \n",
       "3 255      OG0021210   ABD    7    1:1:1              1:1:1          \n",
       "4 259      OG0021208   ABD    7    1:1:1              1:1:1          \n",
       "5 263      OG0021206   ABD    7    1:1:1              1:1:1          \n",
       "6 269      OG0021202   ABD    7    1:1:1              1:1:1          \n",
       "  A                  B                  D                  synteny            \n",
       "1 TraesCS7A01G243100 TraesCS7B01G148400 TraesCS7D01G241900 segmental homeologs\n",
       "2 TraesCS7A01G360600 TraesCS7B01G267100 TraesCS7D01G362400 segmental homeologs\n",
       "3 TraesCS7A01G267100 TraesCS7B01G165300 TraesCS7D01G267800 segmental homeologs\n",
       "4 TraesCS7A01G235400 TraesCS7B01G133600 TraesCS7D01G235500 segmental homeologs\n",
       "5 TraesCS7A01G185200 TraesCS7B01G090200 TraesCS7D01G186600 segmental homeologs\n",
       "6 TraesCS7A01G155800 TraesCS7B01G059900 TraesCS7D01G156200 segmental homeologs\n",
       "  TE HC.LC   source                                              \n",
       "1 NA HC-only OrthoFinder SO,BLAST RBH,DAGchainer,i-ADHoRe,MCScanX\n",
       "2 NA HC-only OrthoFinder SO,BLAST RBH,i-ADHoRe,MCScanX           \n",
       "3 NA HC-only OrthoFinder SO,BLAST RBH,i-ADHoRe                   \n",
       "4 NA HC-only OrthoFinder SO,BLAST RBH,i-ADHoRe,MCScanX           \n",
       "5 NA HC-only OrthoFinder SO,BLAST RBH,DAGchainer                 \n",
       "6 NA HC-only OrthoFinder SO,BLAST RBH,DAGchainer,MCScanX         "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "17036"
      ],
      "text/latex": [
       "17036"
      ],
      "text/markdown": [
       "17036"
      ],
      "text/plain": [
       "[1] 17036"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "homologies<-read.csv(\"./TablesForExploration/HCTriads.csv\", sep=\",\")\n",
    "head(homologies)\n",
    "nrow(homologies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>gene</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>TraesCS1A01G000100  </td></tr>\n",
       "\t<tr><td>TraesCS1A01G000100LC</td></tr>\n",
       "\t<tr><td>TraesCS1A01G000200  </td></tr>\n",
       "\t<tr><td>TraesCS1A01G000200LC</td></tr>\n",
       "\t<tr><td>TraesCS1A01G000300  </td></tr>\n",
       "\t<tr><td>TraesCS1A01G000300LC</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|l}\n",
       " gene\\\\\n",
       "\\hline\n",
       "\t TraesCS1A01G000100  \\\\\n",
       "\t TraesCS1A01G000100LC\\\\\n",
       "\t TraesCS1A01G000200  \\\\\n",
       "\t TraesCS1A01G000200LC\\\\\n",
       "\t TraesCS1A01G000300  \\\\\n",
       "\t TraesCS1A01G000300LC\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "gene | \n",
       "|---|---|---|---|---|---|\n",
       "| TraesCS1A01G000100   | \n",
       "| TraesCS1A01G000100LC | \n",
       "| TraesCS1A01G000200   | \n",
       "| TraesCS1A01G000200LC | \n",
       "| TraesCS1A01G000300   | \n",
       "| TraesCS1A01G000300LC | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  gene                \n",
       "1 TraesCS1A01G000100  \n",
       "2 TraesCS1A01G000100LC\n",
       "3 TraesCS1A01G000200  \n",
       "4 TraesCS1A01G000200LC\n",
       "5 TraesCS1A01G000300  \n",
       "6 TraesCS1A01G000300LC"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "expressed_genes<-read.csv(\"./TablesForExploration/expressed_genes_tpmsOver0.5AtLeast8Samples.csv\")\n",
    "colnames(expressed_genes) <- c(\"gene\",\"tpm\", \"count\")\n",
    "genes_to_use <- data.frame(gene=expressed_genes$gene)\n",
    "head(genes_to_use)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "172034"
      ],
      "text/latex": [
       "172034"
      ],
      "text/markdown": [
       "172034"
      ],
      "text/plain": [
       "[1] 172034"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>gene</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>TraesCS1A01G000100</td></tr>\n",
       "\t<tr><td>TraesCS1A01G000200</td></tr>\n",
       "\t<tr><td>TraesCS1A01G000300</td></tr>\n",
       "\t<tr><td>TraesCS1A01G000400</td></tr>\n",
       "\t<tr><td>TraesCS1A01G000500</td></tr>\n",
       "\t<tr><td>TraesCS1A01G000600</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|l}\n",
       " gene\\\\\n",
       "\\hline\n",
       "\t TraesCS1A01G000100\\\\\n",
       "\t TraesCS1A01G000200\\\\\n",
       "\t TraesCS1A01G000300\\\\\n",
       "\t TraesCS1A01G000400\\\\\n",
       "\t TraesCS1A01G000500\\\\\n",
       "\t TraesCS1A01G000600\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "gene | \n",
       "|---|---|---|---|---|---|\n",
       "| TraesCS1A01G000100 | \n",
       "| TraesCS1A01G000200 | \n",
       "| TraesCS1A01G000300 | \n",
       "| TraesCS1A01G000400 | \n",
       "| TraesCS1A01G000500 | \n",
       "| TraesCS1A01G000600 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  gene              \n",
       "1 TraesCS1A01G000100\n",
       "2 TraesCS1A01G000200\n",
       "3 TraesCS1A01G000300\n",
       "4 TraesCS1A01G000400\n",
       "5 TraesCS1A01G000500\n",
       "6 TraesCS1A01G000600"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "94114"
      ],
      "text/latex": [
       "94114"
      ],
      "text/markdown": [
       "94114"
      ],
      "text/plain": [
       "[1] 94114"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nrow(genes_to_use)\n",
    "hc_genes_to_use<-data.frame(gene=genes_to_use[!grepl(\"LC$\",genes_to_use$gene),])\n",
    "head(hc_genes_to_use)\n",
    "nrow(hc_genes_to_use)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "getSamplesForFactor<- function(metadata, type=\"High.level.tissue\",factor=\"roots\"){\n",
    "    ret<-\"\"\n",
    "    if(type != \"all\"){\n",
    "        ret<-as.character(metadata[metadata[,type] == factor,]$Sample.IDs)\n",
    "    }else{\n",
    "        ret<-as.character(unique(metadata$Sample.IDs))\n",
    "    }\n",
    "    ret\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "getMeansPerFactor<- function(values, metadata,  type=\"High.level.tissue\",factor=\"roots\"){\n",
    "    samples <- getSamplesForFactor(metadata, type, factor)\n",
    "    vals <- values[,samples]\n",
    "    mean<-0\n",
    "    if(length(samples) == 1){\n",
    "        print(\"This factor only has one sample!\")\n",
    "        print(factor)\n",
    "        mean<-vals\n",
    "    }else{\n",
    "        mean<-rowMeans(vals)\n",
    "    }\n",
    "    \n",
    "    \n",
    "    mean<-sort(mean,decreasing=T)\n",
    "    cumulative <- cumsum(mean)\n",
    "    \n",
    "    cumulative<-data.frame(cumulative)\n",
    "    mean<-data.frame(mean)\n",
    "    \n",
    "    mean$gene <- rownames(mean)\n",
    "    cumulative$gene <- rownames(cumulative)\n",
    "    mean$total_samples <- length(samples)\n",
    "    \n",
    "    n <-merge(mean,cumulative, by='gene', all=T)\n",
    "    n <- n[order(n$cumulative,decreasing = F),]\n",
    "    n$seq <- seq(from = 1, to = nrow(n))\n",
    "    n$factor = factor\n",
    "    n\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "getMeansForAllFactors<-function(values, metadata,  type=\"High.level.tissue\"){\n",
    "    factors<-unique(metadata[,type])\n",
    "    f<-factors[1]\n",
    "    meansDFs <- getMeansPerFactor(tpms,metadata,type=type, factor=f)\n",
    "    for (i in 2:length(factors)){\n",
    "        f<-factors[i]\n",
    "        localDF<-getMeansPerFactor(tpms,metadata,type=type, factor=f)\n",
    "        meansDFs <- rbind(meansDFs,localDF)\n",
    "    }\n",
    "    meansDFs\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "isExpressedPerFactor <- function(values, metadata,  type=\"High.level.tissue\",factor=\"roots\", minTPM=0.5){\n",
    "    samples <- getSamplesForFactor(metadata, type, factor)\n",
    "    vals <- values[,samples]\n",
    "    means <- rowMeans(vals)\n",
    "    expr <- means > minTPM\n",
    "    m2 <- data.frame( expressed = expr)\n",
    "    m2$factor<-factor\n",
    "    m2$transcript<-rownames(m2)\n",
    "    m2$total_samples <- length(samples)\n",
    "    m2\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "getExclusiveExpression<-function(values, metadata, minTPM=0.5, type=\"High.level.tissue\"){\n",
    "    means <- getMeansForAllFactors(values, metadata,type=type)\n",
    "    means$expressed<-means$mean > minTPM\n",
    "    exclusiveExpresison<-sqldf(\"SELECT gene, factor, mean, total_samples \n",
    "        FROM means \n",
    "        WHERE expressed \n",
    "        GROUP BY gene HAVING count(factor) = 1 \")\n",
    "   list(means, exclusiveExpresison )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "reshape_triad_groups<-function(triads){\n",
    "    values<-data.frame(gene=triads$A, chr_group=\"A\", group_id=triads$group_id)\n",
    "    values<-rbind(values,data.frame(gene=triads$B, chr_group=\"B\", group_id=triads$group_id))\n",
    "    values<-rbind(values,data.frame(gene=triads$D, chr_group=\"D\", group_id=triads$group_id))\n",
    "    values\n",
    "}\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "get_triad_matrix<-function(normalised_triad, factor=\"all\"){\n",
    "    tmp<-normalised_triad[factor==normalised_triad$factor,]\n",
    "    tmp<-acast(tmp, group_id~chr_group, value.var=\"normalised_triad\") \n",
    "    tmp<-tmp[complete.cases(tmp),]\n",
    "    tmp\n",
    "}\n",
    "plot_normalized_triads<-function(triads){\n",
    "    \n",
    "    p <- ggplot(triads, aes(chr_group, normalised_triad))\n",
    "    p <- p + geom_boxplot(outlier.alpha = 0.05) \n",
    "    p <- p + ylab(\"Contribution\") + xlab(\"Chromosome group\")\n",
    "    p\n",
    "}\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "get_centroids<-function(clust_df){\n",
    "    \n",
    "    clusters<-sort(unique(clust_df$clust))\n",
    "    n_clust<-length(clusters)\n",
    "    dat <- data.frame(\n",
    "        A=numeric(n_clust),B=numeric(n_clust), D=numeric(n_clust), size=numeric(n_clust), \n",
    "        clust=numeric(n_clust), description=character(n_clust),stringsAsFactors=FALSE  )\n",
    "    \n",
    "    for(c in clusters){\n",
    "        tmp_df_clust<-clust_df[clust_df$clust==c,]\n",
    "             \n",
    "        dat[c,1] <- mean(tmp_df_clust$A)\n",
    "        dat[c,2] <- mean(tmp_df_clust$B)\n",
    "        dat[c,3] <- mean(tmp_df_clust$D)\n",
    "        dat[c,4] <- nrow(tmp_df_clust)\n",
    "        dat[c,5] <- c\n",
    "        dat[c,6] <- paste0(\"other.\",c)\n",
    "    }\n",
    "    rownames(dat)<-clusters\n",
    "    dat\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_clust_dist<-function(clust_df, title=\"All\"){\n",
    "    tmp_df<-clust_df[,c(\"A\",\"B\",\"D\",\"group_id\",\"clust\",\"description\")]\n",
    "    tmp_df<-melt(tmp_df,id=c(\"group_id\",\"clust\",\"description\"))\n",
    "    colnames(tmp_df)<-c(\"group_id\",\"clust\",\"description\",\"chr_group\",\"normalised_triad\")\n",
    "    \n",
    "    clusters<-sort(unique(tmp_df$description))\n",
    "    tern <- ggtern(clust_df,aes(A,B,D,color=description, alpha=0.01)) + \n",
    "       geom_point() + theme_legend_position(x = \"topleft\") + scale_alpha(guide = 'none')\n",
    "    gs<-list(tern)\n",
    "    dat <- data.frame(\n",
    "        A=numeric(0),B=numeric(0), D=numeric(0), size=numeric(0),stringsAsFactors=FALSE ) \n",
    "    \n",
    "    \n",
    "    rownames(dat)<-rownames(clusters)\n",
    "    for(c in clusters){\n",
    "        tmp_df_clust<-tmp_df[tmp_df$description==c,]\n",
    "        p<-plot_normalized_triads(tmp_df_clust)\n",
    "        p<- p + ylim(0,1)\n",
    "        p<- p + ylab(\"\") + xlab(\"\")\n",
    "        p<- p + ggtitle(c)\n",
    "        dat[c,1] <- round(100*mean(tmp_df_clust[tmp_df_clust$chr_group==\"A\",\"normalised_triad\"]),digits=2)\n",
    "        dat[c,2] <- round(100*mean(tmp_df_clust[tmp_df_clust$chr_group==\"B\",\"normalised_triad\"]),digits=2)\n",
    "        dat[c,3] <- round(100*mean(tmp_df_clust[tmp_df_clust$chr_group==\"D\",\"normalised_triad\"]),digits=2)\n",
    "        dat[c,4] <- nrow(tmp_df_clust)\n",
    "        gs[[length(gs)+1]] <- p\n",
    "    }\n",
    "    \n",
    "    total_size<-sum(dat$size)\n",
    "    dat$percentage<-round(100*dat$size/total_size,digits=2)\n",
    "    \n",
    "    gs[[length(gs)+1]]<-tableGrob(dat)\n",
    "    lay <- rbind(c( 1, 1, 1, 2, 4, 7),\n",
    "                 c( 1, 1, 1, 3, 5, 8),\n",
    "                 c( 9, 9, 9, 6,NA,NA)\n",
    "                 )\n",
    "\n",
    "    g2 <- arrangeGrob(grobs = gs, layout_matrix = lay, top = title)\n",
    "    g2\n",
    "}\n",
    "\n",
    "\n",
    "#head(mat_centroid_df)\n",
    "#head(mat_centroid_df[,c(\"A\",\"B\",\"D\")])\n",
    "#grid.draw(plot_clust_dist(mat_centroid_df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "get_means_df<-function(metadata, tpms, type=\"High.level.tissue\", min_mean_tpm=0.5){\n",
    "    samples<-getSamplesForFactor(metadata, type=\"all\",factor=\"all\")\n",
    "    values<-data.frame(value=numeric(nrow(tpms)),stringsAsFactors=FALSE)\n",
    "    if(length(samples) > 1){\n",
    "        print(\"All samples:\")\n",
    "        #print(samples)\n",
    "        #print(colnames(tpms))\n",
    "        values$value<-rowMeans(tpms[,samples])\n",
    "    }else{\n",
    "        values$value<-tpms[,samples]\n",
    "    }\n",
    "     \n",
    "    values$factor<-\"all\"\n",
    "    values$gene<-rownames(tpms)\n",
    "    values$samples<-length(samples)\n",
    "    \n",
    "    print(unique(metadata[,type]))\n",
    "    \n",
    "    for(f in as.character(unique(metadata[,type]))){\n",
    "        print(f)\n",
    "        samples<-getSamplesForFactor(metadata, type=type,factor=f)\n",
    "        \n",
    "        tmp<-data.frame(value=numeric(nrow(tpms)),stringsAsFactors=FALSE)\n",
    "        \n",
    "        if(length(samples) > 1){\n",
    "            tmp$value<-rowMeans(tpms[,samples])\n",
    "        }else{\n",
    "            tmp$value<-tpms[,samples]\n",
    "        }\n",
    "        \n",
    "        tmp$factor<-f\n",
    "        tmp$gene<-rownames(tpms)\n",
    "        tmp$samples<-length(samples)\n",
    "        print(colnames(values))\n",
    "        print(colnames(tmp))\n",
    "        values<-rbind(values,tmp)\n",
    "    }\n",
    "    \n",
    "    casted<-dcast(values, gene~factor, value.var=\"value\")\n",
    "    casted$all<-NULL\n",
    "   \n",
    "    rownames(casted)<-casted$gene\n",
    "    casted$gene<-NULL\n",
    "    casted<-as.matrix(casted)\n",
    "    print(head(casted))\n",
    "    tmp<-data.frame(value=numeric(nrow(casted)),stringsAsFactors=FALSE)\n",
    "    \n",
    "    tmp$value<-rowMeans(casted)\n",
    "    tmp$factor<-\"all_means\"\n",
    "    tmp$gene<-rownames(casted)\n",
    "    tmp$samples<-ncol(casted)\n",
    "    print(\"...\")\n",
    "    print(colnames(values))\n",
    "    print(colnames(tmp))\n",
    "    values<-rbind(values,tmp)\n",
    "    \n",
    "    print(head(casted))\n",
    "    print(head(min_mean_tpm))\n",
    "    casted<-ifelse(casted < min_mean_tpm, NA, casted)\n",
    "    print(head(casted))\n",
    "    \n",
    "    tmp<-data.frame(value=numeric(nrow(casted)),stringsAsFactors=FALSE)\n",
    "    tmp$value<-rowMeans(casted, na.rm = TRUE)\n",
    "   \n",
    "    tmp$factor<-\"all_mean_filter\"\n",
    "    tmp$gene<-rownames(casted)\n",
    "    tmp$samples<-rowSums(!is.na(casted))\n",
    "    print(\"...\")\n",
    "    print(colnames(values))\n",
    "    print(colnames(tmp))\n",
    "    \n",
    "    values<-rbind(values,tmp)\n",
    "    values\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "get_means_df<-function(metadata, tpms, triads, type=\"High.level.tissue\"){\n",
    "    samples<-getSamplesForFactor(metadata, type=\"all\",factor=\"all\")\n",
    "    values<-data.frame(value=numeric(nrow(tpms)),stringsAsFactors=FALSE)\n",
    "    if(length(samples) > 1){\n",
    "        values$value<-rowMeans(tpms[,samples])\n",
    "    }else{\n",
    "        values$value<-tpms[,samples]\n",
    "    }\n",
    "     \n",
    "    values$factor<-\"all\"\n",
    "    values$gene<-rownames(tpms)\n",
    "    values$samples<-length(samples)\n",
    "    \n",
    "    #print(unique(metadata[,type]))\n",
    "    \n",
    "    for(f in unique(metadata[,type])){\n",
    "        #print(f)\n",
    "        samples<-getSamplesForFactor(metadata, type=type,factor=f)\n",
    "        \n",
    "        tmp<-data.frame(value=numeric(nrow(tpms)),stringsAsFactors=FALSE)\n",
    "        \n",
    "        if(length(samples) > 1){\n",
    "            tmp$value<-rowMeans(tpms[,samples])\n",
    "        }else{\n",
    "            tmp$value<-tpms[,samples]\n",
    "        }\n",
    "        \n",
    "        tmp$factor<-f\n",
    "        tmp$gene<-rownames(tpms)\n",
    "        tmp$samples<-length(samples)\n",
    "        values<-rbind(values,tmp)\n",
    "    }\n",
    "    \n",
    "    casted<-dcast(values, gene~factor, value.var=\"value\")\n",
    "    casted$all<-NULL\n",
    "   \n",
    "    rownames(casted)<-casted$gene\n",
    "    casted$gene<-NULL\n",
    "    casted<-as.matrix(casted)\n",
    " \n",
    "    tmp<-data.frame(value=numeric(nrow(casted)),stringsAsFactors=FALSE)\n",
    "    tmp$value<-rowMeans(casted)\n",
    "    tmp$factor<-\"all_means\"\n",
    "    tmp$gene<-rownames(casted)\n",
    "    tmp$samples<-ncol(casted)\n",
    "    values<-rbind(values,tmp)\n",
    "    \n",
    "    triads_flat<-reshape_triad_groups(triads)\n",
    "    #print(nrow(values))\n",
    "    vals<-sqldf(\"SELECT `values`.*, chr_group, group_id FROM `values` \n",
    "LEFT JOIN triads_flat on `values`.gene = triads_flat.gene\")\n",
    "    vals\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "123"
      ],
      "text/latex": [
       "123"
      ],
      "text/markdown": [
       "123"
      ],
      "text/plain": [
       "[1] 123"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 51108\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1653 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 6178 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1653 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1653 rows containing non-finite values (stat_boxplot).”"
     ]
    }
   ],
   "source": [
    "\n",
    "get_triad_details<-function(normalised_triad, factor=\"all\",min_triad_sum=1){\n",
    "    tmp<-normalised_triad[factor==normalised_triad$factor,]\n",
    "    tmp<-tmp[tmp$triad_sum>min_triad_sum,]\n",
    "    norm_triad_mat<-get_triad_matrix(tmp, factor=factor)\n",
    "    centroid<-t(as.matrix(colMeans(norm_triad_mat,)))\n",
    "    \n",
    "    dists<-rdist(norm_triad_mat,centroid)\n",
    "    rownames(dists)<-rownames(norm_triad_mat)\n",
    "    colnames(dists)<-c(\"Distance\")\n",
    "    \n",
    "    ranked_dists<-as.matrix(rank(dists)/length(dists))\n",
    "    rownames(ranked_dists)<-rownames(dists)\n",
    "    colnames(ranked_dists)<-c(\"P rank\")\n",
    "    rank_df<-data.frame(dists,ranked_dists)  \n",
    "    rank_df$group_id<-rownames(rank_df)\n",
    "    ranked_triads<-sqldf(\"SELECT tmp.*, Distance, `P.rank` FROM tmp NATURAL JOIN rank_df \")\n",
    "    \n",
    "    t_centroid <- tableGrob(round(centroid*100,digits=2))\n",
    "    p_dist<-ggplot(data=rank_df, aes(Distance)) + geom_histogram(bins=100,aes(y =..density..)) +geom_density(col=2) \n",
    "    p_rank<-ggplot(rank_df,aes(Distance,P.rank))+  geom_bin2d(bins = 75)+theme(legend.position=\"bottom\")\n",
    "    p_genome_dist<-plot_normalized_triads(tmp)\n",
    "    \n",
    "    s<-sample(1:nrow(rank_df), 1500,  replace=FALSE)\n",
    "    rank_s <- rank_df[s,]\n",
    "    title<-paste0(\"Genome contribution in triads: \", factor)\n",
    "    p<-arrangeGrob(p_genome_dist, p_dist, t_centroid, p_rank\n",
    "             , ncol=2, top = title)\n",
    "    \n",
    "    list(distances=rank_df,triad=tmp ,centroid=centroid, p_dist=p_dist, dist_sample=rank_s, \n",
    "         p_rank=p_rank, p_genome_dist=p_genome_dist, plot=p, ranked_triads=ranked_triads, matrix=norm_triad_mat) \n",
    "}\n",
    "\n",
    "\n",
    "get_clusters_by_distance_and_plots<-function(normalized_triads, \n",
    "                                             factor=\"all\", \n",
    "                                             output_prefix=\"\",\n",
    "                                             min_triad_sum=1,\n",
    "                                             title=\"Triads\"){\n",
    "    f<-factor\n",
    "    triad_test<-get_triad_details(normalized_triads,factor=factor,min_triad_sum=min_triad_sum)\n",
    "    triad<-triad_test$ranked_triads\n",
    "    \n",
    "    test_mat<-triad_test$matrix\n",
    "    mat_df<-data.frame(test_mat)\n",
    "    mat_df$group_id<-rownames(test_mat)\n",
    "\n",
    "    title<-paste0(title, \": \" , f)\n",
    "    \n",
    "    centers<-t(matrix(c(0.33,0.33,0.33,1,0,0,0,1,0,0,0,1,0,0.5,0.5,0.5,0,0.5,0.5,0.5,0), nrow=3))\n",
    "    colnames(centers)<-c(\"A\",\"B\",\"D\")\n",
    "    rownames(centers)<-c(\"Central\",\"A.dominant\",\"B.dominant\",\"D.dominant\",\"A.suppressed\",\"B.suppressed\",\"D.suppressed\")\n",
    "    \n",
    "    expectation_distance<-rdist(test_mat,centers)\n",
    "    colnames(expectation_distance)<-c(\"Central\",\n",
    "                                      \"A.dominant\",  \"B.dominant\",  \"D.dominant\",\n",
    "                                      \"A.suppressed\",\"B.suppressed\",\"D.suppressed\")\n",
    "    rownames(expectation_distance)<-rownames(test_mat)\n",
    "    triad_test$expectation_distance<-expectation_distance\n",
    "    mins<-apply( expectation_distance, 1, which.min)\n",
    "    clust_desc<-colnames(expectation_distance)\n",
    "    name_mins<-clust_desc[mins]\n",
    "    \n",
    "    \n",
    "    general_desc<-c(\"Central\",\"Dominant\",  \"Dominant\",  \"Dominant\",\n",
    "                    \"Suppressed\",\"Suppressed\",\"Suppressed\")\n",
    "    \n",
    "    general_name_mins<-general_desc[mins]\n",
    "    \n",
    "    mat_df$clust<-mins\n",
    "    mat_df$description<-name_mins\n",
    "    mat_df$general_description<-general_name_mins\n",
    "    \n",
    "    \n",
    "    mat_df<-cbind(mat_df,expectation_distance)\n",
    "    triad_test$triad<-mat_df\n",
    "    centroids<-get_centroids(mat_df)\n",
    "    \n",
    "    centroids$description<-c(\"Central\",\"A.dominant\",\"B.dominant\",\"D.dominant\",\n",
    "                             \"A.suppressed\",\"B.suppressed\",\"D.suppressed\" \n",
    "                             )\n",
    "    \n",
    "    total_genes<-sum(centroids$size)\n",
    "    centroids[8,\"A\"]<-triad_test$centroid[,\"A\"]\n",
    "    centroids[8,\"B\"]<-triad_test$centroid[,\"B\"]\n",
    "    centroids[8,\"D\"]<-triad_test$centroid[,\"D\"]\n",
    "    centroids[8,\"description\"] <- \"Global\"\n",
    "    centroids[8,\"size\"]<-total_genes\n",
    "    centroids[8,\"clust\"]<-8\n",
    "    centroids$factor<-factor\n",
    "    centroids$percentage<-centroids$size/total_genes\n",
    "    \n",
    "    p<-plot_clust_dist(mat_df, title=title)\n",
    "    \n",
    "    table_save<-sqldf(\"SELECT * FROM mat_df\n",
    "    NATURAL LEFT JOIN triad   \n",
    "    ORDER BY triad.group_id, triad.chr_group\")\n",
    "    table_save<-table_save[ , !(names(table_save) %in% c(\"A\",\"B\",\"D\"))]\n",
    "    table_save$min_triad_sum<-min_triad_sum\n",
    "    triad_test$triad<-table_save\n",
    "    triad_test$centroids<-centroids\n",
    "    triad_test$plot_cluster<-p\n",
    "    centroids$samples<-mean(triad_test$triad$samples)\n",
    "    centroids$min_triad_sum<-min_triad_sum\n",
    "    if(length(output_prefix) > 0){\n",
    "        f2<-gsub(\"\\n\",\"-\",f)\n",
    "        f2<-gsub(\"/\",\"_\",f)\n",
    "        p_filename<-paste0(output_prefix,\"_\",f2,\"_triad_dist.pdf\")\n",
    "        ggsave(p_filename,plot=triad_test$plot,  width = 30, height = 25, units = c(\"cm\"))\n",
    "                            \n",
    "        p_filename<-paste0(output_prefix,\"_\",f2,\"_triad_cluster.pdf\")\n",
    "        ggsave(p_filename,plot=p,  width = 30, height = 25, units = c(\"cm\"))\n",
    "                            \n",
    "        t_filename<-paste0(output_prefix,\"_\",f2,\"_triad_cluster.txt\")\n",
    "        write.table(table_save, file=t_filename,\n",
    "            sep=\"\\t\",quote=TRUE,row.names=FALSE, na=\"\" )\n",
    "        \n",
    "        t_filename<-paste0(output_prefix,\"_\",f2,\"_triad_centroids.txt\")\n",
    "        write.table(centroids, file=t_filename,\n",
    "            sep=\"\\t\",quote=TRUE,row.names=FALSE, na=\"\")\n",
    "    }\n",
    "    triad_test\n",
    "}\n",
    "\n",
    "\n",
    "ks.central<-function(meanTPMS, factor=\"all\"){\n",
    "    m<-meanTPMS[meanTPMS$factor==factor,c(\"description\",\"chr_group\",\"value\")]\n",
    "    \n",
    "    test_df<-NULL\n",
    "    for(gr in c(\"A\",\"B\",\"D\")){\n",
    "        central<-m[m$description==\"Central\" & m$chr_group==gr, \"value\"]\n",
    "        for(d in unique(meanTPMS$description)){\n",
    "            compare<-m[m$description==d & m$chr_group==gr, \"value\"]\n",
    "            for(alt in c(\"less\",\"greater\")){\n",
    "                test<-ks.test(central, compare, alternative=alt)\n",
    "                test$chr_grp<-gr\n",
    "                test$compare<-d\n",
    "                test$alternative<-alt\n",
    "                tmp<-data.frame(t(unlist(test)),stringsAsFactors = FALSE)\n",
    "                if(is.null(test_df)){\n",
    "                    test_df<-tmp\n",
    "                }else{\n",
    "                    colnames(tmp)<-colnames(test_df)\n",
    "                    test_df<-rbind(test_df,tmp)\n",
    "                }\n",
    "            }\n",
    "            \n",
    "        }\n",
    "    }\n",
    "    #print(test_df$p.value)\n",
    "    test_df$fdr_by<-p.adjust(as.numeric(test_df$p.value),method=\"BY\")\n",
    "    test_df$bonferroni<-p.adjust(as.numeric(test_df$p.value),method=\"bonferroni\")\n",
    "    test_df$factor<-factor\n",
    "    test_df\n",
    "}\n",
    "\n",
    "plotMeansTPM<-function(meansTPM, title=\"Test\", filename=\"test\"){\n",
    "    m1<-meansTPM[meansTPM$factor==\"all_mean_filter\",]\n",
    "    m2<-meansTPM[meansTPM$factor!=\"all\",]\n",
    "    m2<-m2[m2$factor!=\"all_means\",]\n",
    "    m2<-m2[m2$factor!=\"all_mean_filter\",]\n",
    "    \n",
    "    p <- ggplot(m1, aes(description, value, fill=chr_group))\n",
    "    p <- p + geom_boxplot(outlier.alpha = 0.05) \n",
    "    p <- p + ylim(c(0, 50))\n",
    "    p <- p + theme(axis.text.x = element_text(angle = 90, hjust = 1),\n",
    "                  legend.position=\"bottom\")\n",
    "    p1 <- p + ylab(\"TPM\") + xlab(\"Category\") \n",
    "\n",
    "    \n",
    "    p <- ggplot(m2, aes(description, value))\n",
    "    p <- p + geom_boxplot(outlier.alpha = 0.05) \n",
    "    p <- p + ylim(c(0, 50)) + theme(axis.text.x = element_text(angle = 90, hjust = 1))\n",
    "    p <- p + facet_wrap(~ factor,  drop = TRUE)\n",
    "    \n",
    "    p2 <- p + ylab(\"TPM\") + xlab(\"Category\") \n",
    "    \n",
    "    p <- ggplot(m1, aes(description, value))\n",
    "    p <- p + geom_boxplot(outlier.alpha = 0.05) \n",
    "    p <- p + ylim(c(0, 50)) \n",
    "    p  <-p + theme(axis.text.x = element_text(angle = 90, hjust = 1),\n",
    "                  legend.position=\"bottom\")\n",
    "    p3 <- p + ylab(\"TPM\") + xlab(\"\") \n",
    "\n",
    "    \n",
    "    p <- ggplot(m1, aes(general_description, value, fill=chr_group))\n",
    "    p <- p + geom_boxplot(outlier.alpha = 0.05) \n",
    "    p <- p + ylim(c(0, 50))\n",
    "    p <- p + theme(axis.text.x = element_text(angle = 90, hjust = 1),\n",
    "                  legend.position=\"bottom\")\n",
    "    p4 <- p + ylab(\"TPM\") + xlab(\"Category\") \n",
    "    \n",
    "    lay <- rbind(c( 3,2,2,4),\n",
    "                 c( 1,2,2,NA))\n",
    "    \n",
    "    g2 <- arrangeGrob(p1,p2, p3,p4,layout_matrix = lay, top = title)\n",
    "    \n",
    "    \n",
    "    p_filename<-paste0(filename,\"_TPM_by_triad_disribution.pdf\")\n",
    "    ggsave(p_filename,plot=g2,  width = 30, height = 25, units = c(\"cm\"))\n",
    "}\n",
    "\n",
    "\n",
    "get_normalized_triad_from_clusters<-function(meansTPM, triads){\n",
    "    \n",
    "    m2<-meansTPM[meansTPM$factor!=\"all\",]\n",
    "    m2<-m2[m2$factor!=\"all_means\",]\n",
    "    tpms<-dcast(m2, gene~factor, value.var=\"value\",drop=FALSE, fill=NA)\n",
    "    rownames(tpms)<-tpms$gene\n",
    "    tpms$gene<-NULL\n",
    "    \n",
    "    values<-data.frame(value=numeric(nrow(tpms)),stringsAsFactors=FALSE)\n",
    "   \n",
    "    values$value<-rowMeans(tpms, na.rm = TRUE)\n",
    "   \n",
    "    values$factor<-\"all_mean_filter\"\n",
    "    values$gene<-rownames(tpms)\n",
    "    values$samples<-rowSums(!is.na(tpms))\n",
    "    \n",
    "    triads_flat<-reshape_triad_groups(triads)\n",
    "    \n",
    "    means_group<-sqldf(\"SELECT `values`.*, chr_group, group_id FROM `values` \n",
    "LEFT JOIN triads_flat on `values`.gene = triads_flat.gene\")\n",
    "    \n",
    "    sums_per_group<-sqldf(\"SELECT group_id , factor, sum(value) as triad_sum \n",
    "FROM means_group GROUP BY group_id , factor\")\n",
    "    \n",
    "    normalized_triads<-sqldf(\"SELECT means_group.*, triad_sum, value/triad_sum as normalised_triad \n",
    "FROM means_group LEFT JOIN sums_per_group \n",
    "ON  sums_per_group.group_id = means_group.group_id \n",
    "AND sums_per_group.factor = means_group.factor \")\n",
    "    \n",
    "    normalized_triads\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "get_triads_and_plots<-function(metadata,tpms, homologies, genes_to_use, \n",
    "                               type=\"High.level.tissue\", \n",
    "                               dataset=\"All\",\n",
    "                               folder=\"./Figures/GenomeDominance/latest\", \n",
    "                               min_triad_sums=c(10), \n",
    "                               min_mean_triad_sums=c(0.5) ){\n",
    "    triadas_with_genes<-sqldf(\"SELECT * from homologies \n",
    "WHERE synteny='segmental homeologs' AND cardinality_abs = '1:1:1'\n",
    "AND \n",
    "(A in genes_to_use \n",
    "OR B in genes_to_use \n",
    "OR D in genes_to_use )\n",
    "\")\n",
    "\n",
    "    #print(head(homologies))\n",
    "    #print(nrow(triadas_with_genes))\n",
    "    #print(head(genes_to_use))\n",
    "    tpms_for_triads<-tpms[c(as.character(triadas_with_genes$A), \n",
    "         as.character(triadas_with_genes$B),\n",
    "         as.character(triadas_with_genes$D)),]\n",
    "    print(nrow(tpms_for_triads))\n",
    "    #print(head(tpms_for_triads))\n",
    "    means_group<-get_means_df(metadata, tpms_for_triads,triadas_with_genes, type=type)\n",
    "    #means_group<-get_means_df(metadata, tpms_for_triads,type=type, min_mean_tpm=0.5)\n",
    "    \n",
    "    sums_per_group<-sqldf(\"SELECT group_id , factor, sum(value) as triad_sum \n",
    "FROM means_group GROUP BY group_id , factor\")\n",
    "    \n",
    "    normalized_triads<-sqldf(\"SELECT means_group.*, triad_sum, value/triad_sum as normalised_triad \n",
    "FROM means_group LEFT JOIN sums_per_group \n",
    "ON  sums_per_group.group_id = means_group.group_id \n",
    "AND sums_per_group.factor = means_group.factor \")\n",
    "    \n",
    "    \n",
    "    path<-paste0(folder,\"/\",dataset,\"/\",type,\"/\")\n",
    "    \n",
    "    for(min_triad_sum in min_triad_sums){\n",
    "        local_triads<-NULL\n",
    "        ks_ret<-NULL\n",
    "        centroids<-NULL\n",
    "        triads<-NULL\n",
    "        path<-paste0(folder,\"/\",dataset,\"/\",type,\"/min_tpm_sum_\",min_triad_sum,\"/\")\n",
    "        path<-gsub(\" \",\"_\",path)\n",
    "        dir.create(path, showWarnings = TRUE, recursive = TRUE, mode = \"0777\")\n",
    "        path<-paste0(folder,\"/\",dataset,\"/\",type,\"/min_tpm_sum_\",min_triad_sum,\"/min_tpm_sum_\",min_triad_sum)\n",
    "        \n",
    "        for(f in unique(normalized_triads$factor)){\n",
    "            title<-paste0(dataset,\" triads:\" , f, \".\\nMinimum triad TPM sum: \", min_triad_sum)\n",
    "            clusters<-get_clusters_by_distance_and_plots(normalized_triads, factor=f, title=title, \n",
    "                                                         min_triad_sum=min_triad_sum,\n",
    "                           output_prefix=path)\n",
    "            suppressWarnings(\n",
    "                tmp<-ks.central(clusters$triad, factor=f)\n",
    "            )\n",
    "            if(is.null(centroids)){\n",
    "                centroids<-clusters$centroids\n",
    "                triads<-clusters$triad\n",
    "               \n",
    "                ks_ret<-tmp\n",
    "            }else{\n",
    "                centroids<-rbind(centroids,clusters$centroids)\n",
    "                triads<-rbind(triads,clusters$triad) \n",
    "                 \n",
    "                ks_ret<-rbind(ks_ret, tmp)\n",
    "            }\n",
    "           \n",
    "        }\n",
    "        \n",
    "        \n",
    "        normalized_triads_filter<-get_normalized_triad_from_clusters(triads,triadas_with_genes )\n",
    "        title<-paste0(dataset,\" triads:\" , f, \".\\nMinimum triad TPM sum: \", min_triad_sum)\n",
    "        \n",
    "        \n",
    "        clusters<-get_clusters_by_distance_and_plots(normalized_triads_filter, factor=\"all_mean_filter\", \n",
    "                                                     title=title, \n",
    "                                                     min_triad_sum=min_triad_sum,\n",
    "                                                     output_prefix=path)\n",
    "        \n",
    "        centroids<-rbind(centroids,clusters$centroids)\n",
    "        triads<-rbind(triads,clusters$triad) \n",
    "        \n",
    "        suppressWarnings(\n",
    "            tmp<-ks.central(clusters$triad, factor=\"all_mean_filter\")\n",
    "        )\n",
    "        ks_ret<-rbind(ks_ret, tmp)\n",
    "        \n",
    "        t_filename<-paste0(path,\"min_tpm_sum\",min_triad_sum)\n",
    "        title<-paste0(\"Category distribution \", dataset, \"\\n\", type, \"\\n\", \"Min TPM average sum: \", min_triad_sum)\n",
    "        plotMeansTPM(triads, filename=t_filename, title=title)\n",
    "        \n",
    "        t_filename<-paste0(path, \"min_tpm_sum_\", min_triad_sum, \"_ks.txt\")\n",
    "        write.table(ks_ret, file=t_filename, sep=\"\\t\",quote=TRUE,row.names=FALSE, na=\"\")\n",
    "        t_filename<-paste0(path,\"summary_triad_centroids.txt\")\n",
    "        write.table(centroids, file=t_filename, sep=\"\\t\",quote=TRUE,row.names=FALSE, na=\"\")\n",
    "        t_filename<-paste0(path,\"summary_triad_cluster.txt\")\n",
    "        write.table(triads, file=t_filename, sep=\"\\t\",quote=TRUE,row.names=FALSE, na=\"\")\n",
    "    }\n",
    "}\n",
    "\n",
    "metadata<-read.csv(\"./TablesForExploration/metadata.csv\", row.names = 1, sep=\"\\t\")\n",
    "metadata <- metadata_used[metadata_used[,\"Variety\"] == \"Chinese Spring\",]\n",
    "metadata <- metadata[metadata[,\"High.level.stress.disease\"] == \"none\",]\n",
    "      \n",
    "nrow(metadata)\n",
    "get_triads_and_plots(metadata,tpms,homologies,hc_genes_to_use, type=\"High.level.tissue\",\n",
    "                               dataset=\"HC_CS_no_stress\" , min_triad_sum=c(1))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>Sample.IDs</th><th scope=col>subset</th><th scope=col>type_to_use</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>Sample_10   </td><td>Development </td><td>Intermediate</td></tr>\n",
       "\t<tr><td>Sample_18   </td><td>Development </td><td>Intermediate</td></tr>\n",
       "\t<tr><td>Sample_3A   </td><td>Development </td><td>Intermediate</td></tr>\n",
       "\t<tr><td>Sample_12   </td><td>Development </td><td>Intermediate</td></tr>\n",
       "\t<tr><td>Sample_26   </td><td>Development </td><td>Intermediate</td></tr>\n",
       "\t<tr><td>Sample_6A   </td><td>Development </td><td>Intermediate</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       " Sample.IDs & subset & type\\_to\\_use\\\\\n",
       "\\hline\n",
       "\t Sample\\_10  & Development  & Intermediate\\\\\n",
       "\t Sample\\_18  & Development  & Intermediate\\\\\n",
       "\t Sample\\_3A  & Development  & Intermediate\\\\\n",
       "\t Sample\\_12  & Development  & Intermediate\\\\\n",
       "\t Sample\\_26  & Development  & Intermediate\\\\\n",
       "\t Sample\\_6A  & Development  & Intermediate\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "Sample.IDs | subset | type_to_use | \n",
       "|---|---|---|---|---|---|\n",
       "| Sample_10    | Development  | Intermediate | \n",
       "| Sample_18    | Development  | Intermediate | \n",
       "| Sample_3A    | Development  | Intermediate | \n",
       "| Sample_12    | Development  | Intermediate | \n",
       "| Sample_26    | Development  | Intermediate | \n",
       "| Sample_6A    | Development  | Intermediate | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  Sample.IDs subset      type_to_use \n",
       "1 Sample_10  Development Intermediate\n",
       "2 Sample_18  Development Intermediate\n",
       "3 Sample_3A  Development Intermediate\n",
       "4 Sample_12  Development Intermediate\n",
       "5 Sample_26  Development Intermediate\n",
       "6 Sample_6A  Development Intermediate"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "datasets<-read.csv(\"./TablesForExploration/samples_to_use.csv\")\n",
    "head(datasets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "subsets<-unique(datasets$subset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Development\"\n",
      "[1] Intermediate\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 209\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1873 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 34562 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1873 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1873 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"850_samples\"\n",
      "[1] Intermediate\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 850\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1553 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 44924 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1553 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1553 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"CS_no_stress\"\n",
      "[1] Intermediate\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in dir.create(path):\n",
      "“'./Figures/GenomeDominance/latest/HC_CS_no_stress' already exists”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 123\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1549 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 20930 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1549 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1549 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"CS_NB_inc_stress\"\n",
      "[1] Intermediate\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 144\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1539 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 21208 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1539 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1539 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"abiotic\"\n",
      "[1] Intermediate_Stress\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 50\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1411 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 13782 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1411 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1411 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"disease\"\n",
      "[1] Intermediate_Stress\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 163\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1777 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 24223 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1777 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1777 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"grain\"\n",
      "[1] Intermediate\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 119\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 875 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 7701 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 875 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 875 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"leaf\"\n",
      "[1] Intermediate\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 245\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1649 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 21178 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1649 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1649 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"root\"\n",
      "[1] Intermediate\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 45\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1850 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 3753 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1850 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1850 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"spike\"\n",
      "[1] Intermediate\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 128\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1806 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 13419 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1806 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1806 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"abiotic_merged_control\"\n",
      "[1] Intermediate_Stress_merged_control\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 50\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1446 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 9863 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1446 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1446 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"disease_merged_control\"\n",
      "[1] Intermediate_Stress_merged_control\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 163\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1832 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 16025 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1832 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1832 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"stress_control\"\n",
      "[1] Intermediate_Stress\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 77\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1562 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 15126 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1562 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1562 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"abiotic_stress_control\"\n",
      "[1] Intermediate_Stress\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 13\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1226 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 3655 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1226 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1226 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"disease_stress_control\"\n",
      "[1] Intermediate_Stress\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 64\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1692 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 11471 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1692 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1692 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"abiotic_stress\"\n",
      "[1] Intermediate_Stress\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 34\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1450 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 8472 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1450 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1450 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"disease_stress\"\n",
      "[1] Intermediate_Stress\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 106\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1808 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 14032 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1808 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1808 rows containing non-finite values (stat_boxplot).”"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"stress\"\n",
      "[1] Intermediate_Stress\n",
      "3 Levels: Intermediate ... Intermediate_Stress_merged_control\n",
      "[1] 140\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 1659 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 22504 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1659 rows containing non-finite values (stat_boxplot).”Warning message:\n",
      "“Removed 1659 rows containing non-finite values (stat_boxplot).”"
     ]
    }
   ],
   "source": [
    "\n",
    "get_triads_and_plots<-function(metadata,tpms, homologies, genes_to_use, \n",
    "                               type=\"High.level.tissue\", \n",
    "                               dataset=\"All\",\n",
    "                               folder=\"./Figures/GenomeDominance/latest\", \n",
    "                               min_triad_sums=c(10), \n",
    "                               min_mean_triad_sums=c(0.5) ){\n",
    "    triadas_with_genes<-sqldf(\"SELECT * from homologies \n",
    "WHERE synteny='segmental homeologs' AND cardinality_abs = '1:1:1'\n",
    "AND \n",
    "(A in genes_to_use \n",
    "OR B in genes_to_use \n",
    "OR D in genes_to_use )\n",
    "\")\n",
    "    type<-as.character(type)\n",
    "   \n",
    "    tpms_for_triads<-tpms[c(as.character(triadas_with_genes$A), \n",
    "         as.character(triadas_with_genes$B),\n",
    "         as.character(triadas_with_genes$D)),]\n",
    "  \n",
    "    means_group<-get_means_df(metadata, tpms_for_triads,triadas_with_genes, type=type)\n",
    "    \n",
    "    sums_per_group<-sqldf(\"SELECT group_id , factor, sum(value) as triad_sum \n",
    "FROM means_group GROUP BY group_id , factor\")\n",
    "    \n",
    "    normalized_triads<-sqldf(\"SELECT means_group.*, triad_sum, value/triad_sum as normalised_triad \n",
    "FROM means_group LEFT JOIN sums_per_group \n",
    "ON  sums_per_group.group_id = means_group.group_id \n",
    "AND sums_per_group.factor = means_group.factor \")\n",
    "    \n",
    "    \n",
    "    path<-paste0(folder,\"/\",dataset,\"/\",type,\"/\")\n",
    "    \n",
    "    for(min_triad_sum in min_triad_sums){\n",
    "        local_triads<-NULL\n",
    "        ks_ret<-NULL\n",
    "        centroids<-NULL\n",
    "        triads<-NULL\n",
    "        path<-paste0(folder,\"/\",dataset,\"/\",type,\"/min_tpm_sum_\",min_triad_sum,\"/\")\n",
    "        path<-gsub(\" \",\"_\",path)\n",
    "        dir.create(path, showWarnings = TRUE, recursive = TRUE, mode = \"0777\")\n",
    "        path<-paste0(folder,\"/\",dataset,\"/\",type,\"/min_tpm_sum_\",min_triad_sum,\"/min_tpm_sum_\",min_triad_sum)\n",
    "        \n",
    "        for(f in unique(normalized_triads$factor)){\n",
    "            title<-paste0(dataset,\" triads:\" , f, \".\\nMinimum triad TPM sum: \", min_triad_sum)\n",
    "            clusters<-get_clusters_by_distance_and_plots(normalized_triads, factor=f, title=title, \n",
    "                                                         min_triad_sum=min_triad_sum,\n",
    "                           output_prefix=path)\n",
    "            suppressWarnings(\n",
    "                tmp<-ks.central(clusters$triad, factor=f)\n",
    "            )\n",
    "            if(is.null(centroids)){\n",
    "                centroids<-clusters$centroids\n",
    "                triads<-clusters$triad\n",
    "               \n",
    "                ks_ret<-tmp\n",
    "            }else{\n",
    "                centroids<-rbind(centroids,clusters$centroids)\n",
    "                triads<-rbind(triads,clusters$triad) \n",
    "                 \n",
    "                ks_ret<-rbind(ks_ret, tmp)\n",
    "            }\n",
    "           \n",
    "        }\n",
    "        \n",
    "        \n",
    "        normalized_triads_filter<-get_normalized_triad_from_clusters(triads,triadas_with_genes )\n",
    "        title<-paste0(dataset,\" triads:\" , f, \".\\nMinimum triad TPM sum: \", min_triad_sum)\n",
    "        \n",
    "        \n",
    "        clusters<-get_clusters_by_distance_and_plots(normalized_triads_filter, factor=\"all_mean_filter\", \n",
    "                                                     title=title, \n",
    "                                                     min_triad_sum=min_triad_sum,\n",
    "                                                     output_prefix=path)\n",
    "        \n",
    "        centroids<-rbind(centroids,clusters$centroids)\n",
    "        triads<-rbind(triads,clusters$triad) \n",
    "        \n",
    "        suppressWarnings(\n",
    "            tmp<-ks.central(clusters$triad, factor=\"all_mean_filter\")\n",
    "        )\n",
    "        ks_ret<-rbind(ks_ret, tmp)\n",
    "        \n",
    "        t_filename<-paste0(path,\"min_tpm_sum\",min_triad_sum)\n",
    "        title<-paste0(\"Category distribution \", dataset, \"\\n\", type, \"\\n\", \"Min TPM average sum: \", min_triad_sum)\n",
    "        plotMeansTPM(triads, filename=t_filename, title=title)\n",
    "        \n",
    "        t_filename<-paste0(path, \"min_tpm_sum_\", min_triad_sum, \"_ks.txt\")\n",
    "        write.table(ks_ret, file=t_filename, sep=\"\\t\",quote=TRUE,row.names=FALSE, na=\"\")\n",
    "        t_filename<-paste0(path,\"summary_triad_centroids.txt\")\n",
    "        write.table(centroids, file=t_filename, sep=\"\\t\",quote=TRUE,row.names=FALSE, na=\"\")\n",
    "        t_filename<-paste0(path,\"summary_triad_cluster.txt\")\n",
    "        write.table(triads, file=t_filename, sep=\"\\t\",quote=TRUE,row.names=FALSE, na=\"\")\n",
    "    }\n",
    "}\n",
    "\n",
    "for(s in subsets){\n",
    "    print(s)\n",
    "    samples<-datasets[datasets$subset == s, ]\n",
    "    type_to_use <- unique(samples$type_to_use)\n",
    "    print(type_to_use)\n",
    "    metadata <- metadata_used[ metadata_used$Sample.IDs %in% samples$Sample.IDs, ]\n",
    "    path <- paste0(\"./Figures/GenomeDominance/latest/HC_\", s)\n",
    "    dir.create(path)\n",
    "    path <- paste0(path, \"/\",type_to_use,\".txt\")\n",
    "    write.table(metadata, file=path,\n",
    "            sep=\"\\t\",quote=TRUE,row.names=FALSE, na=\"\")\n",
    "    print(nrow(metadata))\n",
    "    ds<-paste0(\"HC_\",s)\n",
    "    get_triads_and_plots(metadata,tpms,homologies,hc_genes_to_use, \n",
    "                         type=type_to_use, \n",
    "                               dataset=ds, min_triad_sum=c(0.5))\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "read_cluster_triads<-function(dir=\"./Figures/GenomeDominance/latest\", min_tpm_sum=1, dataset=\"HC_CS_no_stress\", type=\"High.level.tissue\"){\n",
    "    path<-paste0(dir,\"/\",dataset,\"/\",type,\"/min_tpm_sum_\", min_tpm_sum,\n",
    "               \"/min_tpm_sum_\",min_tpm_sum,\n",
    "               \"summary_triad_cluster.txt\")\n",
    "    table<-read.csv(path, sep=\"\\t\")  \n",
    "    table<-table[table$factor!=\"all\",]\n",
    "    table<-table[table$factor!=\"all_means\",]\n",
    "    table\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Definition of homoeolog expression bias categories\n",
    "\n",
    "The ideal normalized expression bias for the seven categories was defined as:\n",
    "\n",
    "  Category      | A | B | D\n",
    "  ------------- | -------------\n",
    "  Balanced      | 0.33 | 0.33 | 0.33\n",
    "  A suppressed  | 0 | 0.5 | 0.5\n",
    "  B suppressed  | 0.5 | 0 | 0.5\n",
    "  D suppressed  | 0.5 | 0.5 | 0\n",
    "  A dominant    | 1 | 0 | 0\n",
    "  B dominant    | 0 | 1 | 0\n",
    "  D dominant    | 0 | 0 | 1  \n",
    "  \n",
    "  \n",
    "We calculated the Euclidean distance (rdist function from R 3.3.2) from the observed normalized expression of each triad to each of the seven ideal categories listed above. We assigned the homoeolog expression bias category for each triad by selecting the shortest distance. This was done for each of the intermediate tissue as well as for the average across all expressed tissues (combined analysis)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "getMovementForTriad<-function(triad_mat, group_id){\n",
    "    local_triad<-triad_mat[triad_mat$group_id==group_id,]\n",
    "    local_triad<-local_triad[local_triad$factor==\"all_mean_filter\",]\n",
    "    category<-local_triad$description\n",
    "    local_centroid<-as.matrix(local_triad[local_triad$factor==\"all_mean_filter\",c(\"A\",\"B\",\"D\")])\n",
    "    \n",
    "    local_triad<-triad_mat[triad_mat$group_id==group_id,]\n",
    "    local_triad<-local_triad[local_triad$factor!=\"all_mean_filter\" ,]\n",
    "    local_mat<-as.matrix(local_triad[,c(\"A\",\"B\",\"D\")])\n",
    "    \n",
    "    rownames(local_mat)<-local_triad$factor\n",
    "    \n",
    "    \n",
    "    dists<-rdist(local_mat, local_centroid)\n",
    "    f_dists<-rdist(local_mat)\n",
    "   # print(local_mat)\n",
    "    #hull<-convhulln(local_mat,options=\"Qc\")\n",
    "    #print(hull)\n",
    "    \n",
    "    ret<-list(group_id=group_id, \n",
    "              factor_count=nrow(local_triad),\n",
    "              central_total_distance=sum(dists),  \n",
    "              central_mean_distance=mean(dists),\n",
    "              central_max_distance=max(dists),\n",
    "              central_sd_distance=sd(dists),\n",
    "              central_max_over_mean=max(dists)/mean(dists),\n",
    "              \n",
    "              factor_total_distance=sum(f_dists),  \n",
    "              factor_mean_distance=mean(f_dists),\n",
    "              factor_max_distance=max(f_dists),\n",
    "              factor_sd_distance=sd(f_dists),\n",
    "              factor_max_over_mean=max(f_dists)/mean(f_dists),\n",
    "              category=as.character(category),\n",
    "              sum_mean_tpm=sum(local_triad$triad_sum),\n",
    "              total_categories=length(unique(local_triad$description)),\n",
    "              categories=    paste(sort(unique(local_triad$description)), collapse=\", \"),\n",
    "              Central       = sum(local_triad$description == \"Central\"),\n",
    "              A.dominant    = sum(local_triad$description == \"A.dominant\"),\n",
    "              A.suppressed  = sum(local_triad$description == \"A.suppressed\"),\n",
    "              B.dominant    = sum(local_triad$description == \"B.dominant\"),\n",
    "              B.suppressed  = sum(local_triad$description == \"B.suppressed\"),\n",
    "              D.dominant    = sum(local_triad$description == \"D.dominant\"),\n",
    "              D.suppressed  = sum(local_triad$description == \"D.suppressed\")\n",
    "              \n",
    "             )\n",
    "    ret\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Development\"\n",
      "[1] \"850_samples\"\n",
      "[1] \"CS_no_stress\"\n",
      "[1] \"CS_NB_inc_stress\"\n",
      "[1] \"abiotic\"\n",
      "[1] \"disease\"\n",
      "[1] \"grain\"\n",
      "[1] \"leaf\"\n",
      "[1] \"root\"\n",
      "[1] \"spike\"\n",
      "[1] \"abiotic_merged_control\"\n",
      "[1] \"disease_merged_control\"\n",
      "[1] \"stress_control\"\n",
      "[1] \"abiotic_stress_control\"\n",
      "[1] \"disease_stress_control\"\n",
      "[1] \"abiotic_stress\"\n",
      "[1] \"disease_stress\"\n",
      "[1] \"stress\"\n"
     ]
    }
   ],
   "source": [
    "res<-NULL\n",
    "triads<-NULL\n",
    "base_path<-\"./Figures/GenomeDominance/latest/\"\n",
    "for(s in subsets){\n",
    "    print(s)\n",
    "    \n",
    "    samples<-datasets[datasets$subset == s, ]\n",
    "    type <- unique(samples$type_to_use)\n",
    "    min <-0.5\n",
    "    dataset<-paste0(\"HC_\",s)\n",
    "    triads_t<-read_cluster_triads(min=min, type=type,dataset=dataset)\n",
    "    triads_t$dataset<-dataset\n",
    "    triads_df<-dcast(triads_t, group_id+factor+description+triad_sum~chr_group, value.var=\"normalised_triad\")\n",
    "    res_t<-sapply(unique(triads_df$group_id), \n",
    "                  function(x){\n",
    "                      getMovementForTriad(triads_df,x)\n",
    "                  })\n",
    "    res_t<-data.frame(t(res_t))\n",
    "    res_t$dataset<-dataset\n",
    "    res_t$dataset<-as.factor(res_t$dataset)\n",
    "    for(c in colnames(res_t)){\n",
    "        res_t[,c]<-unlist(res_t[,c])\n",
    "    }\n",
    "    \n",
    "    head(res_t)    \n",
    "    #path<-paste0(base_path,dataset,\"/\",type,\"/min_tpm_sum_\",min,\"/triadMovement.csv\")\n",
    "    #write.csv(res_t, path, row.names=FALSE)\n",
    "    if(is.null(res)){\n",
    "        res<-res_t\n",
    "        triads<-triads_t\n",
    "    }else{\n",
    "        res<-rbind(res,res_t)\n",
    "        triads<-rbind(triads,triads_t)\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "saveRDS(triads, \"./TablesForExploration/Triads.rds\")\n",
    "saveRDS(res,\"./TablesForExploration/TriadMovement.rds\" )\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}