<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/creating-a-network-of-human-gene-homology-with-r-and-d3-without-html-widget-for-easy-loading/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://shiring.github.io/genome/2016/12/12/homologous_genes_post_noHTML"> Shirin's playgRound</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><em>The html widget I produced with networksD3 seems to create problems with loading the site for some people, so here is the same blog post as <a target="_blank" rel="nofollow" href="https://shiring.github.io/genome/2016/12/11/homologous_genes_post">yesterday</a> just without the rendered html widget in it.</em></p>
<hr>
<p></p>
<p>This week I want to explore how many of our human genes have homologs in other species. I will use this question to show how to visualize dendrograms and (D3-) networks.</p>
<p></p>
<h2>Accessing species information from Biomart</h2>
<p>I am using the <a target="_blank" rel="nofollow" href="https://bioconductor.org/packages/release/bioc/html/biomaRt.html">biomaRt</a> package to access all datasets available in the <a target="_blank" rel="nofollow" href="http://www.ensembl.org/biomart/martview/64258667e46f2ecac32be58d12091b87">Biomart</a> Ensembl database. The <em>listDatasets()</em> function shows that currently there is data for 69 species (including humans). This dataframe gives us the name of each dataset in Biomart Ensembl, the respective common species name and version number.</p>
<div>
<pre><code><span>library</span><span>(</span><span>biomaRt</span><span>)</span><span>
</span><span>ensembl</span><span> </span><span>=</span><span> </span><span>useMart</span><span>(</span><span>&quot;ensembl&quot;</span><span>)</span><span>
</span><span>datasets</span><span> </span><span>&lt;-</span><span> </span><span>listDatasets</span><span>(</span><span>ensembl</span><span>)</span><span>
</span><span>human</span><span> </span><span>=</span><span> </span><span>useMart</span><span>(</span><span>&quot;ensembl&quot;</span><span>,</span><span> </span><span>dataset</span><span> </span><span>=</span><span> </span><span>&quot;hsapiens_gene_ensembl&quot;</span><span>)</span><span>
</span></code></pre>
</div>
<p>I am then looping over each species to call <em>useMart()</em> and assign an object to each species’ Ensembl database.</p>
<div>
<pre><code><span>for</span><span> </span><span>(</span><span>i</span><span> </span><span>in</span><span> </span><span>1</span><span>:</span><span>nrow</span><span>(</span><span>datasets</span><span>))</span><span> </span><span>{</span><span>
</span><span>ensembl</span><span> </span><span>&lt;-</span><span> </span><span>datasets</span><span>[</span><span>i</span><span>,</span><span> </span><span>1</span><span>]</span><span>
</span><span>assign</span><span>(</span><span>paste0</span><span>(</span><span>ensembl</span><span>),</span><span> </span><span>useMart</span><span>(</span><span>&quot;ensembl&quot;</span><span>,</span><span> </span><span>dataset</span><span> </span><span>=</span><span> </span><span>paste0</span><span>(</span><span>ensembl</span><span>)))</span><span>
</span><span>}</span><span>
</span><span>specieslist</span><span> </span><span>&lt;-</span><span> </span><span>datasets</span><span>$</span><span>dataset</span><span>
</span></code></pre>
</div>
<p></p>
<h2>Identifying human gene homologs</h2>
<h3>Protein coding genes</h3>
<p>The majority of human genes are protein coding genes. This means that their DNA sequence will be translated into a protein with specific cellular functions. For an overview of all gene biotypes in the human genome, have a <a target="_blank" rel="nofollow" href="https://shiring.github.io/genome/2016/10/23/AnnotationDbi">look at this previous post</a>.</p>
<p>To identify all protein coding genes in the human genome, I am using <a target="_blank" rel="nofollow" href="https://bioconductor.org/packages/release/bioc/html/AnnotationDbi.html">AnnotationDbi</a> and the <a target="_blank" rel="nofollow" href="https://bioconductor.org/packages/release/data/annotation/html/EnsDb.Hsapiens.v79.html">EnsDb.Hsapiens.v79</a> database.</p>
<div>
<pre><code><span>library</span><span>(</span><span>AnnotationDbi</span><span>)</span><span>
</span><span>library</span><span>(</span><span>EnsDb.Hsapiens.v79</span><span>)</span><span>
</span><span># get all Ensembl gene IDs
</span><span>human_EnsDb</span><span> </span><span>&lt;-</span><span> </span><span>keys</span><span>(</span><span>EnsDb.Hsapiens.v79</span><span>,</span><span> </span><span>keytype</span><span> </span><span>=</span><span> </span><span>&quot;GENEID&quot;</span><span>)</span><span>
</span><span># get the biotype of each gene ID
</span><span>human_gene_EnsDb</span><span> </span><span>&lt;-</span><span> </span><span>ensembldb</span><span>::</span><span>select</span><span>(</span><span>EnsDb.Hsapiens.v79</span><span>,</span><span> </span><span>keys</span><span> </span><span>=</span><span> </span><span>human_EnsDb</span><span>,</span><span> </span><span>columns</span><span> </span><span>=</span><span> </span><span>&quot;GENEBIOTYPE&quot;</span><span>,</span><span> </span><span>keytype</span><span> </span><span>=</span><span> </span><span>&quot;GENEID&quot;</span><span>)</span><span>
</span><span># and keep only protein coding genes
</span><span>human_prot_coding_genes</span><span> </span><span>&lt;-</span><span> </span><span>human_gene_EnsDb</span><span>[</span><span>which</span><span>(</span><span>human_gene_EnsDb</span><span>$</span><span>GENEBIOTYPE</span><span> </span><span>==</span><span> </span><span>&quot;protein_coding&quot;</span><span>),</span><span> </span><span>]</span><span>
</span></code></pre>
</div>
<p>To get a dataframe of human genes and their homologs in other species, I am looping over all Biomart Ensembl databases using biomaRt’s <em>getLDS()</em> function. <em>getLDS()</em> links two datasets and is equivalent to homology mapping with Ensembl. I am also including the human database so that I can see how many genes from EnsDb.Hsapiens.v79 are found via Biomart.</p>
<div>
<pre><code><span>for</span><span> </span><span>(</span><span>species</span><span> </span><span>in</span><span> </span><span>specieslist</span><span>)</span><span> </span><span>{</span><span>
</span><span>print</span><span>(</span><span>species</span><span>)</span><span>
</span><span>assign</span><span>(</span><span>paste0</span><span>(</span><span>&quot;homologs_human_&quot;</span><span>,</span><span> </span><span>species</span><span>),</span><span> </span><span>getLDS</span><span>(</span><span>attributes</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;ensembl_gene_id&quot;</span><span>,</span><span> </span><span>&quot;chromosome_name&quot;</span><span>),</span><span> </span><span>filters</span><span> </span><span>=</span><span> </span><span>&quot;ensembl_gene_id&quot;</span><span>,</span><span> </span><span>values</span><span> </span><span>=</span><span> </span><span>human_prot_coding_genes</span><span>$</span><span>GENEID</span><span>,</span><span> </span><span>mart</span><span> </span><span>=</span><span> </span><span>human</span><span>,</span><span> </span><span>attributesL</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;ensembl_gene_id&quot;</span><span>,</span><span> </span><span>&quot;chromosome_name&quot;</span><span>),</span><span> </span><span>martL</span><span> </span><span>=</span><span> </span><span>get</span><span>(</span><span>species</span><span>)))</span><span>
</span><span>}</span><span>
</span></code></pre>
</div>
<p>Based on EnsDb.Hsapiens.v79 protein coding genes, I now have 69 individual datasets with homologous genes found in other species. I now want to combine these 69 datasets into one cooccurrence matrix. To do this, I am first joining the homologous gene subsets back to the original dataframe with all protein coding genes. Then, I convert existing homologous gene names to 1 and NAs to 0. These, I then merge into one table. This final table now contains the information whether each gene has a homolog in each of the other species (“1”) or not (“0”).</p>
<div>
<pre><code><span>library</span><span>(</span><span>dplyr</span><span>)</span><span>
</span><span>for</span><span> </span><span>(</span><span>i</span><span> </span><span>in</span><span> </span><span>1</span><span>:</span><span>length</span><span>(</span><span>specieslist</span><span>)){</span><span>
</span><span>species</span><span> </span><span>&lt;-</span><span> </span><span>specieslist</span><span>[</span><span>i</span><span>]</span><span>
</span><span>homologs_human</span><span> </span><span>&lt;-</span><span> </span><span>left_join</span><span>(</span><span>human_prot_coding_genes</span><span>,</span><span> </span><span>get</span><span>(</span><span>paste0</span><span>(</span><span>&quot;homologs_human_&quot;</span><span>,</span><span> </span><span>species</span><span>)),</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;GENEID&quot;</span><span> </span><span>=</span><span> </span><span>&quot;Ensembl.Gene.ID&quot;</span><span>))</span><span>
</span><span>homologs_human</span><span>[,</span><span> </span><span>paste0</span><span>(</span><span>species</span><span>)]</span><span> </span><span>&lt;-</span><span> </span><span>ifelse</span><span>(</span><span>is.na</span><span>(</span><span>homologs_human</span><span>$</span><span>Ensembl.Gene.ID.1</span><span>),</span><span> </span><span>0</span><span>,</span><span> </span><span>1</span><span>)</span><span>
</span><span>homologs_human</span><span> </span><span>&lt;-</span><span> </span><span>homologs_human</span><span>[,</span><span> </span><span>c</span><span>(</span><span>1</span><span>,</span><span> </span><span>6</span><span>)]</span><span>
</span><span>homologs_human</span><span> </span><span>&lt;-</span><span> </span><span>homologs_human</span><span>[</span><span>!</span><span>duplicated</span><span>(</span><span>homologs_human</span><span>$</span><span>GENEID</span><span>),</span><span> </span><span>]</span><span>
</span><span>if</span><span> </span><span>(</span><span>i</span><span> </span><span>==</span><span> </span><span>1</span><span>){</span><span>
</span><span>homologs_human_table</span><span> </span><span>&lt;-</span><span> </span><span>homologs_human</span><span>
</span><span>}</span><span> </span><span>else</span><span> </span><span>{</span><span>
</span><span>homologs_human_table</span><span> </span><span>&lt;-</span><span> </span><span>left_join</span><span>(</span><span>homologs_human_table</span><span>,</span><span> </span><span>homologs_human</span><span>,</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>&quot;GENEID&quot;</span><span>)</span><span>
</span><span>}</span><span>
</span><span>}</span><span>
</span></code></pre>
</div>
<p>I only want to keep the 21,925 human genes found via Biomart (77 genes were not found).</p>
<div>
<pre><code><span>homologs_human_table</span><span> </span><span>&lt;-</span><span> </span><span>homologs_human_table</span><span>[</span><span>which</span><span>(</span><span>homologs_human_table</span><span>[,</span><span> </span><span>grep</span><span>(</span><span>&quot;sapiens&quot;</span><span>,</span><span> </span><span>colnames</span><span>(</span><span>homologs_human_table</span><span>))]</span><span> </span><span>==</span><span> </span><span>1</span><span>),</span><span> </span><span>]</span><span>
</span></code></pre>
</div>
<p>The final steps to creating a cooccurrence matrix are</p>
<ul>
<li>removing the column with Ensembl gene IDs (because we have the Homo sapiens column from Biomart),</li>
<li>multiplying with the transposed matrix and</li>
<li>setting all rows and columns which don’t show cooccurrences with humans to 0 (because I only looked at homology from the human perspective).</li>
</ul>
<div>
<pre><code><span>gene_matrix</span><span> </span><span>&lt;-</span><span> </span><span>homologs_human_table</span><span>[,</span><span> </span><span>-1</span><span>]</span><span>
</span><span>co_occurrence</span><span> </span><span>&lt;-</span><span> </span><span>t</span><span>(</span><span>as.matrix</span><span>(</span><span>gene_matrix</span><span>))</span><span> </span><span>%*%</span><span> </span><span>as.matrix</span><span>(</span><span>gene_matrix</span><span>)</span><span>
</span><span>co_occurrence</span><span>[</span><span>-</span><span>grep</span><span>(</span><span>&quot;sapiens&quot;</span><span>,</span><span> </span><span>rownames</span><span>(</span><span>co_occurrence</span><span>)),</span><span> </span><span>-</span><span>grep</span><span>(</span><span>&quot;sapiens&quot;</span><span>,</span><span> </span><span>colnames</span><span>(</span><span>co_occurrence</span><span>))]</span><span> </span><span>&lt;-</span><span> </span><span>0</span><span>
</span></code></pre>
</div>
<p>The first information I want to extract from this data is how many and what proportion of human genes had a homolog in each of the other species.</p>
<div>
<pre><code><span>genes</span><span> </span><span>&lt;-</span><span> </span><span>data.frame</span><span>(</span><span>organism</span><span> </span><span>=</span><span> </span><span>colnames</span><span>(</span><span>gene_matrix</span><span>),</span><span>
</span><span>number_genes</span><span> </span><span>=</span><span> </span><span>colSums</span><span>(</span><span>gene_matrix</span><span>),</span><span>
</span><span>proportion_human_genes</span><span> </span><span>=</span><span> </span><span>colSums</span><span>(</span><span>gene_matrix</span><span>)</span><span>/</span><span>nrow</span><span>(</span><span>homologs_human_table</span><span>))</span><span>
</span></code></pre>
</div>
<p>I also want to know which species are similar regarding the gene homology they share with humans. This, I will visualize with a hierarchical clustering dendrogram using the <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/dendextend/vignettes/introduction.html">dendextend</a> and <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/circlize/index.html">circlize</a> packages.</p>
<p>Before I can produce plots however, I want to create annotation attributes for each species:</p>
<ul>
<li>their common name (because most people won’t know off the top of their heads what species T. nigroviridis is) and</li>
<li>a grouping attribute, e.g. mammal, fish, bird, etc.</li>
</ul>
<p>In order to have the correct order of species, I create the annotation dataframe based off the dendrogram.</p>
<p>Most of the common names I can extract from the list of Biomart Ensembl datasets but a few I have to change manually.</p>
<div>
<pre><code><span>library</span><span>(</span><span>dendextend</span><span>)</span><span>
</span><span>library</span><span>(</span><span>circlize</span><span>)</span><span>
</span><span># create a dendrogram
</span><span>h</span><span> </span><span>&lt;-</span><span> </span><span>hclust</span><span>(</span><span>dist</span><span>(</span><span>scale</span><span>(</span><span>as.matrix</span><span>(</span><span>t</span><span>(</span><span>gene_matrix</span><span>))),</span><span> </span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;manhattan&quot;</span><span>))</span><span>
</span><span>dend</span><span> </span><span>&lt;-</span><span> </span><span>as.dendrogram</span><span>(</span><span>h</span><span>)</span><span>
</span><span>library</span><span>(</span><span>dplyr</span><span>)</span><span>
</span><span>labels</span><span> </span><span>&lt;-</span><span> </span><span>as.data.frame</span><span>(</span><span>dend</span><span> </span><span>%>%</span><span> </span><span>labels</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>left_join</span><span>(</span><span>datasets</span><span>,</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;dend %>% labels&quot;</span><span> </span><span>=</span><span> </span><span>&quot;dataset&quot;</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>left_join</span><span>(</span><span>genes</span><span>,</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;dend %>% labels&quot;</span><span> </span><span>=</span><span> </span><span>&quot;organism&quot;</span><span>))</span><span>
</span><span>labels</span><span>[,</span><span> </span><span>2</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>gsub</span><span>(</span><span>&quot;(.*)( genes (.*))&quot;</span><span>,</span><span> </span><span>&quot;\\1&quot;</span><span>,</span><span> </span><span>labels</span><span>[,</span><span> </span><span>2</span><span>])</span><span>
</span><span>labels</span><span>$</span><span>group</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>rep</span><span>(</span><span>&quot;mammal&quot;</span><span>,</span><span> </span><span>2</span><span>),</span><span> </span><span>rep</span><span>(</span><span>&quot;fish&quot;</span><span>,</span><span> </span><span>8</span><span>),</span><span> </span><span>&quot;amphibia&quot;</span><span>,</span><span> </span><span>rep</span><span>(</span><span>&quot;fish&quot;</span><span>,</span><span> </span><span>4</span><span>),</span><span> </span><span>rep</span><span>(</span><span>&quot;bird&quot;</span><span>,</span><span> </span><span>5</span><span>),</span><span> </span><span>rep</span><span>(</span><span>&quot;reptile&quot;</span><span>,</span><span> </span><span>2</span><span>),</span><span> </span><span>rep</span><span>(</span><span>&quot;mammal&quot;</span><span>,</span><span> </span><span>41</span><span>),</span><span> </span><span>&quot;fungus&quot;</span><span>,</span><span> </span><span>&quot;lamprey&quot;</span><span>,</span><span> </span><span>rep</span><span>(</span><span>&quot;seasquirt&quot;</span><span>,</span><span> </span><span>2</span><span>),</span><span> </span><span>&quot;nematode&quot;</span><span>,</span><span> </span><span>&quot;insect&quot;</span><span>)</span><span>
</span><span>labels</span><span>$</span><span>description</span><span>[</span><span>grep</span><span>(</span><span>&quot;hedgehog&quot;</span><span>,</span><span> </span><span>labels</span><span>$</span><span>description</span><span>)]</span><span> </span><span>&lt;-</span><span> </span><span>&quot;Hedgehog Tenrec&quot;</span><span>
</span><span>labels</span><span>$</span><span>description</span><span>[</span><span>grep</span><span>(</span><span>&quot;Saccharomyces&quot;</span><span>,</span><span> </span><span>labels</span><span>$</span><span>description</span><span>)]</span><span> </span><span>&lt;-</span><span> </span><span>&quot;Yeast&quot;</span><span>
</span><span>labels</span><span>$</span><span>description</span><span>[</span><span>grep</span><span>(</span><span>&quot;savignyi&quot;</span><span>,</span><span> </span><span>labels</span><span>$</span><span>description</span><span>)]</span><span> </span><span>&lt;-</span><span> </span><span>&quot;C. savignyi&quot;</span><span>
</span><span>labels</span><span>$</span><span>description</span><span>[</span><span>grep</span><span>(</span><span>&quot;intestinalis&quot;</span><span>,</span><span> </span><span>labels</span><span>$</span><span>description</span><span>)]</span><span> </span><span>&lt;-</span><span> </span><span>&quot;C. intestinalis&quot;</span><span>
</span><span>labels</span><span>$</span><span>description</span><span>[</span><span>grep</span><span>(</span><span>&quot;elegans&quot;</span><span>,</span><span> </span><span>labels</span><span>$</span><span>description</span><span>)]</span><span> </span><span>&lt;-</span><span> </span><span>&quot;C. elegans&quot;</span><span>
</span><span>labels</span><span>$</span><span>description</span><span>[</span><span>grep</span><span>(</span><span>&quot;turtle&quot;</span><span>,</span><span> </span><span>labels</span><span>$</span><span>description</span><span>)]</span><span> </span><span>&lt;-</span><span> </span><span>&quot;Turtle&quot;</span><span>
</span><span>labels</span><span>$</span><span>description</span><span>[</span><span>grep</span><span>(</span><span>&quot;Vervet&quot;</span><span>,</span><span> </span><span>labels</span><span>$</span><span>description</span><span>)]</span><span> </span><span>&lt;-</span><span> </span><span>&quot;Vervet&quot;</span><span>
</span></code></pre>
</div>
<p>The first plot shows the proportion of human genes with a homolog in each of the other species.</p>
<div>
<pre><code><span>library</span><span>(</span><span>ggplot2</span><span>)</span><span>
</span><span>my_theme</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>base_size</span><span> </span><span>=</span><span> </span><span>12</span><span>,</span><span> </span><span>base_family</span><span> </span><span>=</span><span> </span><span>&quot;sans&quot;</span><span>){</span><span>
</span><span>theme_grey</span><span>(</span><span>base_size</span><span> </span><span>=</span><span> </span><span>base_size</span><span>,</span><span> </span><span>base_family</span><span> </span><span>=</span><span> </span><span>base_family</span><span>)</span><span> </span><span>+</span><span>
</span><span>theme</span><span>(</span><span>
</span><span>axis.text</span><span> </span><span>=</span><span> </span><span>element_text</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>12</span><span>),</span><span>
</span><span>axis.text.x</span><span> </span><span>=</span><span> </span><span>element_text</span><span>(</span><span>angle</span><span> </span><span>=</span><span> </span><span>90</span><span>,</span><span> </span><span>vjust</span><span> </span><span>=</span><span> </span><span>0.5</span><span>,</span><span> </span><span>hjust</span><span> </span><span>=</span><span> </span><span>1</span><span>),</span><span>
</span><span>axis.title</span><span> </span><span>=</span><span> </span><span>element_text</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>14</span><span>),</span><span>
</span><span>panel.grid.major</span><span> </span><span>=</span><span> </span><span>element_line</span><span>(</span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;grey&quot;</span><span>),</span><span>
</span><span>panel.grid.minor</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>(),</span><span>
</span><span>panel.background</span><span> </span><span>=</span><span> </span><span>element_rect</span><span>(</span><span>fill</span><span> </span><span>=</span><span> </span><span>&quot;aliceblue&quot;</span><span>),</span><span>
</span><span>strip.background</span><span> </span><span>=</span><span> </span><span>element_rect</span><span>(</span><span>fill</span><span> </span><span>=</span><span> </span><span>&quot;lightgrey&quot;</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;grey&quot;</span><span>,</span><span> </span><span>size</span><span> </span><span>=</span><span> </span><span>1</span><span>),</span><span>
</span><span>strip.text</span><span> </span><span>=</span><span> </span><span>element_text</span><span>(</span><span>face</span><span> </span><span>=</span><span> </span><span>&quot;bold&quot;</span><span>,</span><span> </span><span>size</span><span> </span><span>=</span><span> </span><span>12</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;navy&quot;</span><span>),</span><span>
</span><span>legend.position</span><span> </span><span>=</span><span> </span><span>&quot;bottom&quot;</span><span>,</span><span>
</span><span>legend.background</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>(),</span><span>
</span><span>panel.margin</span><span> </span><span>=</span><span> </span><span>unit</span><span>(</span><span>.5</span><span>,</span><span> </span><span>&quot;lines&quot;</span><span>),</span><span>
</span><span>panel.border</span><span> </span><span>=</span><span> </span><span>element_rect</span><span>(</span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;grey&quot;</span><span>,</span><span> </span><span>fill</span><span> </span><span>=</span><span> </span><span>NA</span><span>,</span><span> </span><span>size</span><span> </span><span>=</span><span> </span><span>0.5</span><span>)</span><span>
</span><span>)</span><span>
</span><span>}</span><span>
</span><span>labels_p</span><span> </span><span>&lt;-</span><span> </span><span>labels</span><span>[</span><span>-1</span><span>,</span><span> </span><span>]</span><span>
</span><span>f</span><span> </span><span>=</span><span> </span><span>as.factor</span><span>(</span><span>labels_p</span><span>[</span><span>order</span><span>(</span><span>labels_p</span><span>$</span><span>proportion_human_genes</span><span>,</span><span> </span><span>decreasing</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>),</span><span> </span><span>&quot;description&quot;</span><span>])</span><span>
</span><span>labels_p</span><span> </span><span>&lt;-</span><span> </span><span>within</span><span>(</span><span>labels_p</span><span>,</span><span> </span><span>description</span><span> </span><span>&lt;-</span><span> </span><span>factor</span><span>(</span><span>description</span><span>,</span><span> </span><span>levels</span><span> </span><span>=</span><span> </span><span>f</span><span>))</span><span>
</span><span>ggplot</span><span>(</span><span>labels_p</span><span>,</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>description</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>proportion_human_genes</span><span>,</span><span> </span><span>fill</span><span> </span><span>=</span><span> </span><span>group</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_bar</span><span>(</span><span>stat</span><span> </span><span>=</span><span> </span><span>&quot;identity&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>my_theme</span><span>()</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>
</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Proportion of homology&quot;</span><span>,</span><span>
</span><span>fill</span><span> </span><span>=</span><span> </span><span>&quot;&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Human gene homology&quot;</span><span>,</span><span>
</span><span>subtitle</span><span> </span><span>=</span><span> </span><span>&quot;Proportion of human protein coding genes with homologs in 68 other species&quot;</span><span>
</span><span>)</span><span>
</span></code></pre>
</div>
<p><img src="https://i0.wp.com/shiring.github.io/genome/2016/12/12/homologous_genes_post_files/figure-markdown_github/unnamed-chunk-12-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>As expected, the highest proportions of homologous genes are found between humans and most other mammals: between 81 and 70% of our genes have a homolog in most mammals. The lowest proportion of homologous genes are found in yeast, but even such a simply organism shares more than 20% of its genes with us. Think about that the next time you enjoy a delicious freshly baked bread or drink a nice cold glass of beer…</p>
<p>However, based on these Biomart Ensembl databases, more of our genes have a homolog in mice than in primates. I would have expected primates to come out on top, but what we see here might be confounded by the different amounts of information we have on species’ genomes and the accuracy of this information: there has been more research done on the mouse than e.g. on the alpaca. We therefore have more and more accurate information on the mouse genome and have identified more homologous genes than in other species.</p>
<p></p>
<p>For plotting the dendrograms, I am adding a column with color specifications for each group to the annotation table.</p>
<div>
<pre><code><span>labels</span><span>$</span><span>col</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>rep</span><span>(</span><span>&quot;aquamarine4&quot;</span><span>,</span><span> </span><span>2</span><span>),</span><span> </span><span>rep</span><span>(</span><span>&quot;darkorange3&quot;</span><span>,</span><span> </span><span>8</span><span>),</span><span> </span><span>&quot;blue4&quot;</span><span>,</span><span> </span><span>rep</span><span>(</span><span>&quot;darkorange3&quot;</span><span>,</span><span> </span><span>4</span><span>),</span><span> </span><span>rep</span><span>(</span><span>&quot;darkred&quot;</span><span>,</span><span> </span><span>5</span><span>),</span><span> </span><span>rep</span><span>(</span><span>&quot;darkmagenta&quot;</span><span>,</span><span> </span><span>2</span><span>),</span><span> </span><span>rep</span><span>(</span><span>&quot;aquamarine4&quot;</span><span>,</span><span> </span><span>41</span><span>),</span><span> </span><span>&quot;darkslateblue&quot;</span><span>,</span><span> </span><span>&quot;deepskyblue1&quot;</span><span>,</span><span> </span><span>rep</span><span>(</span><span>&quot;deeppink3&quot;</span><span>,</span><span> </span><span>2</span><span>),</span><span> </span><span>&quot;forestgreen&quot;</span><span>,</span><span> </span><span>&quot;brown3&quot;</span><span>)</span><span>
</span></code></pre>
</div>
<p>This dendrogram shows the hierarchical clustering of human gene homology with 68 other species. The species are grouped by color (see legend).</p>
<div>
<pre><code><span>labels</span><span>(</span><span>dend</span><span>)</span><span> </span><span>&lt;-</span><span> </span><span>labels</span><span>$</span><span>description</span><span>
</span><span>dend</span><span> </span><span>%>%</span><span> </span><span>set</span><span>(</span><span>&quot;labels_col&quot;</span><span>,</span><span> </span><span>labels</span><span>$</span><span>col</span><span>)</span><span> </span><span>%>%</span><span> </span><span>set</span><span>(</span><span>&quot;labels_cex&quot;</span><span>,</span><span> </span><span>2</span><span>)</span><span> </span><span>%>%</span><span> </span><span>circlize_dendrogram</span><span>(</span><span>labels_track_height</span><span> </span><span>=</span><span> </span><span>NA</span><span>,</span><span> </span><span>dend_track_height</span><span> </span><span>=</span><span> </span><span>0.3</span><span>)</span><span> </span><span>legend</span><span>(</span><span>&quot;topleft&quot;</span><span>,</span><span> </span><span>unique</span><span>(</span><span>labels</span><span>$</span><span>group</span><span>),</span><span> </span><span>pch</span><span> </span><span>=</span><span> </span><span>19</span><span>,</span><span>
</span><span>col</span><span> </span><span>=</span><span> </span><span>unique</span><span>(</span><span>labels</span><span>$</span><span>col</span><span>),</span><span> </span><span>pt.cex</span><span> </span><span>=</span><span> </span><span>3</span><span>,</span><span> </span><span>cex</span><span> </span><span>=</span><span> </span><span>3</span><span>,</span><span> </span><span>bty</span><span> </span><span>=</span><span> </span><span>&quot;n&quot;</span><span>,</span><span> </span><span>ncol</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span></code></pre>
</div>
<p><img src="https://i0.wp.com/shiring.github.io/genome/2016/12/12/homologous_genes_post_files/figure-markdown_github/unnamed-chunk-15-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>The most distant group from all other species include yeast, fruit fly, lamprey, C. elegans and the two seasquirt species, followed by tarsier, alpaca, sloth and shrew. The next group includes the fish, reptiles, amphibiens, birds and the platypus. Interestingly, the platypus the Xenopus are more similar to fish than to birds and reptiles.</p>
<p>And finally, we have the remaining mammals: here, the primates cluster nicely together.</p>
<p></p>
<p>In order to show the homology between human genes with the other species, we can also plot a network. The 2D-network is created from the cooccurrence matrix with node size and edge width representing the proportion of homology and colors depicting the species’ group.</p>
<div>
<pre><code><span># plot the network graph
</span><span>library</span><span>(</span><span>igraph</span><span>)</span><span>
</span><span>g</span><span> </span><span>&lt;-</span><span> </span><span>graph_from_adjacency_matrix</span><span>(</span><span>co_occurrence</span><span>,</span><span>
</span><span>weighted</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span>
</span><span>diag</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>,</span><span>
</span><span>mode</span><span> </span><span>=</span><span> </span><span>&quot;undirected&quot;</span><span>)</span><span>
</span><span>g</span><span> </span><span>&lt;-</span><span> </span><span>simplify</span><span>(</span><span>g</span><span>,</span><span> </span><span>remove.multiple</span><span> </span><span>=</span><span> </span><span>F</span><span>,</span><span> </span><span>remove.loops</span><span> </span><span>=</span><span> </span><span>T</span><span>,</span><span> </span><span>edge.attr.comb</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>weight</span><span> </span><span>=</span><span> </span><span>&quot;sum&quot;</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>&quot;ignore&quot;</span><span>))</span><span>
</span><span>labels_2</span><span> </span><span>&lt;-</span><span> </span><span>labels</span><span>[</span><span>match</span><span>(</span><span>rownames</span><span>(</span><span>co_occurrence</span><span>),</span><span> </span><span>labels</span><span>[,</span><span> </span><span>1</span><span>]),</span><span> </span><span>]</span><span>
</span><span>V</span><span>(</span><span>g</span><span>)</span><span>$</span><span>color</span><span> </span><span>&lt;-</span><span> </span><span>labels_2</span><span>$</span><span>col</span><span>
</span><span>V</span><span>(</span><span>g</span><span>)</span><span>$</span><span>label</span><span> </span><span>&lt;-</span><span> </span><span>labels_2</span><span>$</span><span>description</span><span>
</span><span>V</span><span>(</span><span>g</span><span>)</span><span>$</span><span>size</span><span> </span><span>&lt;-</span><span> </span><span>labels_2</span><span>$</span><span>proportion_human_genes</span><span>*</span><span>25</span><span>
</span><span>E</span><span>(</span><span>g</span><span>)</span><span>$</span><span>arrow.size</span><span> </span><span>&lt;-</span><span> </span><span>0.2</span><span>
</span><span>E</span><span>(</span><span>g</span><span>)</span><span>$</span><span>edge.color</span><span> </span><span>&lt;-</span><span> </span><span>&quot;gray80&quot;</span><span>
</span><span>E</span><span>(</span><span>g</span><span>)</span><span>$</span><span>width</span><span> </span><span>&lt;-</span><span> </span><span>labels_2</span><span>$</span><span>proportion_human_genes</span><span>*</span><span>10</span><span>
</span></code></pre>
</div>
<div>
<pre><code><span>plot</span><span>(</span><span>g</span><span>,</span><span>
</span><span>vertex.label.font</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span>
</span><span>vertex.shape</span><span> </span><span>=</span><span> </span><span>&quot;sphere&quot;</span><span>,</span><span>
</span><span>vertex.label.cex</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span>
</span><span>vertex.label.color</span><span> </span><span>=</span><span> </span><span>&quot;white&quot;</span><span>,</span><span>
</span><span>vertex.frame.color</span><span> </span><span>=</span><span> </span><span>NA</span><span>)</span><span>
</span><span>legend</span><span>(</span><span>&quot;topleft&quot;</span><span>,</span><span> </span><span>unique</span><span>(</span><span>labels_2</span><span>$</span><span>group</span><span>),</span><span> </span><span>pch</span><span> </span><span>=</span><span> </span><span>19</span><span>,</span><span>
</span><span>col</span><span> </span><span>=</span><span> </span><span>unique</span><span>(</span><span>labels_2</span><span>$</span><span>col</span><span>),</span><span> </span><span>pt.cex</span><span> </span><span>=</span><span> </span><span>2</span><span>,</span><span> </span><span>cex</span><span> </span><span>=</span><span> </span><span>1.5</span><span>,</span><span> </span><span>bty</span><span> </span><span>=</span><span> </span><span>&quot;n&quot;</span><span>,</span><span> </span><span>ncol</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span></code></pre>
</div>
<p><img src="https://i0.wp.com/shiring.github.io/genome/2016/12/12/homologous_genes_post_files/figure-markdown_github/unnamed-chunk-17-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Another nice way to visualize networks is to create a D3 JavaScript network graph with the <a target="_blank" rel="nofollow" href="https://christophergandrud.github.io/networkD3/">networkD3</a> package. It shows the same information as the 2D-network above but here you can interact with the nodes, which makes the species names easier to read.</p>
<div>
<pre><code><span>library</span><span>(</span><span>networkD3</span><span>)</span><span>
</span><span>net_d3</span><span> </span><span>&lt;-</span><span> </span><span>igraph_to_networkD3</span><span>(</span><span>g</span><span>,</span><span> </span><span>group</span><span> </span><span>=</span><span> </span><span>labels_2</span><span>$</span><span>group</span><span>)</span><span>
</span><span>net_d3</span><span>$</span><span>nodes</span><span> </span><span>&lt;-</span><span> </span><span>merge</span><span>(</span><span>net_d3</span><span>$</span><span>nodes</span><span>,</span><span> </span><span>genes</span><span>,</span><span> </span><span>by.x</span><span> </span><span>=</span><span> </span><span>&quot;name&quot;</span><span>,</span><span> </span><span>by.y</span><span> </span><span>=</span><span> </span><span>&quot;row.names&quot;</span><span>)</span><span>
</span><span>net_d3</span><span>$</span><span>nodes</span><span>$</span><span>proportion_human_genes</span><span> </span><span>&lt;-</span><span> </span><span>net_d3</span><span>$</span><span>nodes</span><span>$</span><span>proportion_human_genes</span><span>
</span><span>net_d3</span><span>$</span><span>nodes</span><span> </span><span>&lt;-</span><span> </span><span>merge</span><span>(</span><span>net_d3</span><span>$</span><span>nodes</span><span>,</span><span> </span><span>labels_2</span><span>[,</span><span> </span><span>1</span><span>:</span><span>2</span><span>],</span><span> </span><span>by.x</span><span> </span><span>=</span><span> </span><span>&quot;name&quot;</span><span>,</span><span> </span><span>by.y</span><span> </span><span>=</span><span> </span><span>&quot;dend %>% labels&quot;</span><span>)</span><span>
</span><span>net_d3</span><span>$</span><span>nodes</span><span>[,</span><span> </span><span>1</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>net_d3</span><span>$</span><span>nodes</span><span>$</span><span>description</span><span>
</span><span>net_d3</span><span>$</span><span>nodes</span><span> </span><span>&lt;-</span><span> </span><span>net_d3</span><span>$</span><span>nodes</span><span>[</span><span>match</span><span>(</span><span>labels_2</span><span>[,</span><span> </span><span>1</span><span>],</span><span> </span><span>net_d3</span><span>$</span><span>nodes</span><span>$</span><span>organism</span><span>),</span><span> </span><span>]</span><span>
</span><span># Create force directed network plot
</span><span>forceNetwork</span><span>(</span><span>Links</span><span> </span><span>=</span><span> </span><span>net_d3</span><span>$</span><span>links</span><span>,</span><span> </span><span>Nodes</span><span> </span><span>=</span><span> </span><span>net_d3</span><span>$</span><span>nodes</span><span>,</span><span>
</span><span>Nodesize</span><span> </span><span>=</span><span> </span><span>&quot;proportion_human_genes&quot;</span><span>,</span><span>
</span><span>radiusCalculation</span><span> </span><span>=</span><span> </span><span>JS</span><span>(</span><span>&quot; Math.sqrt(d.nodesize)*15&quot;</span><span>),</span><span>
</span><span>linkDistance</span><span> </span><span>=</span><span> </span><span>100</span><span>,</span><span>
</span><span>zoom</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>opacity</span><span> </span><span>=</span><span> </span><span>0.9</span><span>,</span><span>
</span><span>Source</span><span> </span><span>=</span><span> </span><span>'source'</span><span>,</span><span> </span><span>Target</span><span> </span><span>=</span><span> </span><span>'target'</span><span>,</span><span>
</span><span>linkWidth</span><span> </span><span>=</span><span> </span><span>networkD3</span><span>::</span><span>JS</span><span>(</span><span>&quot;function(d) { return d.value/5; }&quot;</span><span>),</span><span>
</span><span>bounded</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>,</span><span> </span><span>colourScale</span><span> </span><span>=</span><span> </span><span>JS</span><span>(</span><span>&quot;d3.scale.category10()&quot;</span><span>),</span><span>
</span><span>NodeID</span><span> </span><span>=</span><span> </span><span>'name'</span><span>,</span><span> </span><span>Group</span><span> </span><span>=</span><span> </span><span>'group'</span><span>,</span><span> </span><span>fontSize</span><span> </span><span>=</span><span> </span><span>10</span><span>,</span><span> </span><span>charge</span><span> </span><span>=</span><span> </span><span>-200</span><span>,</span><span>
</span><span>opacityNoHover</span><span> </span><span>=</span><span> </span><span>0.7</span><span>)</span><span>
</span></code></pre>
</div>
<p>If you want to see the rendered html widget, check out the <a target="_blank" rel="nofollow" href="https://shiring.github.io/genome/2016/12/11/homologous_genes_post">original post</a></p>
<hr>
<p></p>
<div>
<pre><code>## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Sierra 10.12
## 
## locale:
## [1] de_DE.UTF-8/de_DE.UTF-8/de_DE.UTF-8/C/de_DE.UTF-8/de_DE.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] networkD3_0.2.13         igraph_1.0.1            
##  [3] ggplot2_2.2.0            dplyr_0.5.0             
##  [5] circlize_0.3.9           dendextend_1.3.0        
##  [7] EnsDb.Hsapiens.v79_1.1.0 ensembldb_1.6.0         
##  [9] GenomicFeatures_1.26.0   GenomicRanges_1.26.0    
## [11] GenomeInfoDb_1.10.0      AnnotationDbi_1.36.0    
## [13] IRanges_2.8.0            S4Vectors_0.12.0        
## [15] Biobase_2.34.0           BiocGenerics_0.20.0     
## [17] biomaRt_2.30.0          
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.2.1                    jsonlite_1.1                 
##  [3] AnnotationHub_2.6.0           shiny_0.14.2                 
##  [5] assertthat_0.1                interactiveDisplayBase_1.12.0
##  [7] Rsamtools_1.26.0              yaml_2.1.14                  
##  [9] robustbase_0.92-6             RSQLite_1.0.0                
## [11] lattice_0.20-34               digest_0.6.10                
## [13] XVector_0.14.0                colorspace_1.3-1             
## [15] htmltools_0.3.5               httpuv_1.3.3                 
## [17] Matrix_1.2-7.1                plyr_1.8.4                   
## [19] XML_3.98-1.5                  zlibbioc_1.20.0              
## [21] xtable_1.8-2                  mvtnorm_1.0-5                
## [23] scales_0.4.1                  whisker_0.3-2                
## [25] BiocParallel_1.8.0            tibble_1.2                   
## [27] SummarizedExperiment_1.4.0    nnet_7.3-12                  
## [29] lazyeval_0.2.0                magrittr_1.5                 
## [31] mime_0.5                      mclust_5.2                   
## [33] evaluate_0.10                 MASS_7.3-45                  
## [35] class_7.3-14                  BiocInstaller_1.24.0         
## [37] tools_3.3.2                   GlobalOptions_0.0.10         
## [39] trimcluster_0.1-2             stringr_1.1.0                
## [41] kernlab_0.9-25                munsell_0.4.3                
## [43] cluster_2.0.5                 fpc_2.1-10                   
## [45] Biostrings_2.42.0             grid_3.3.2                   
## [47] RCurl_1.95-4.8                htmlwidgets_0.8              
## [49] bitops_1.0-6                  labeling_0.3                 
## [51] rmarkdown_1.1                 gtable_0.2.0                 
## [53] flexmix_2.3-13                DBI_0.5-1                    
## [55] R6_2.2.0                      GenomicAlignments_1.10.0     
## [57] knitr_1.15                    prabclus_2.2-6               
## [59] rtracklayer_1.34.0            shape_1.4.2                  
## [61] modeltools_0.2-21             stringi_1.1.2                
## [63] Rcpp_0.12.8                   DEoptimR_1.0-8               
## [65] diptest_0.75-7
</code></pre>
</div>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://shiring.github.io/genome/2016/12/12/homologous_genes_post_noHTML"> Shirin's playgRound</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/VFG1x8sj8Zw"></a></div>