Clustering_R.ipynb

{
 "cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["library(ggplot2)\n", "library(pheatmap)"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["source(\"SimData.R\")"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["set.seed(123)\n", "simDat = simulate2Group(n=40, p=2, effect=c(1, 0.75))\n", "xsim = simDat$x\n", "ysim = simDat$y"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["## -----------------------------------------------------------------\n", "## k-means\n", "## -----------------------------------------------------------------\n", "kmSim = kmeans(xsim, centers=2)"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["kmplot = function(xy) {\n", "    x = xy$x\n", "    y = xy$y\n", "    km = kmeans(x, centers=2)\n", "    ggdata = data.frame(x, cluster=factor(km$cluster), y=y)\n", "    ggobj = ggplot(data=ggdata, mapping=aes(\n", "        x = g1,\n", "        y = g2,\n", "        color = cluster,\n", "        shape = y\n", "    )) + theme_classic()\n", "    ggobj = ggobj + geom_point(size=3)\n", "    ggobj = ggobj + scale_shape_manual(values=c(6, 17))\n", "    print(ggobj)\n", "}"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["kmplot(simulate2Group(n=40, p=2, effect=c(10, 0)))"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["kmplot(simulate2Group(n=40, p=2, effect=c(1, 0.75)))"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["## -----------------------------------------------------------------\n", "## hierarchical clustering\n", "## -----------------------------------------------------------------\n", "simData2 = simulate2Group(n=40, p=20, effect=c(2, 1, 1))\n", "xsim2 = simData2$x\n", "ysim2 = simData2$y"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["## cluster pseudosamples\n", "ihcSim = hclust(xsim2)  ## generates error -- hclust wants distance matrix\n", "                        ## not raw data!\n", "xdist = dist(xsim2, method=\"euclidean\")\n", "ihcSim = hclust(xdist, method=\"average\")\n", "plot(ihcSim)"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["## cluster pseudogenes\n", "ghcSim = hclust(dist(t(xsim2), method=\"euclidean\"), method=\"average\")\n", "plot(ghcSim)"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["## -----------------------------------------------------------------\n", "## clustered heatmap\n", "## -----------------------------------------------------------------\n", "heatY = data.frame(row.names=rownames(xsim2), group=ysim2)\n", "pheatmap(t(xsim2), annotation=heatY,\n", "         annotation_colors=list(group=c(A='black', B=rgb(1, 0, 0.4))))"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["## -----------------------------------------------------------------\n", "## on real data...\n", "## -----------------------------------------------------------------\n", "load(\"prepared_datasets.RData\")"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["hcShen = hclust(dist(xnorms$shen), method=\"complete\")\n", "## pdf(\"ShenHClust.pdf\", h=6, w=7.25)\n", "plot(hcShen)\n", "## garbage = dev.off()"]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": ["load(\"shenGeneAnnot.RData\")\n", "shenHighVar = colnames(xnorms$shen)[apply(xnorms$shen, 2, sd) > 2]\n", "heatX = t(xnorms$shen[ , shenHighVar])\n", "rownames(heatX) = ifelse(\n", "    rownames(heatX) %in% as.character(shenGeneAnnot$gene),\n", "    shenGeneSyms[rownames(heatX)],\n", "    rownames(heatX)\n", ")\n", "## remove overall gene-means from data for more useful plot\n", "heatX = data.frame(sweep(heatX, 1, rowMeans(heatX)), check.names=FALSE)\n", "## pay attention to changes around mean, not far from it\n", "maxLogFoldChange = 2.5\n", "heatX[heatX > maxLogFoldChange] = maxLogFoldChange\n", "heatX[heatX < -maxLogFoldChange] = -maxLogFoldChange\n", "heatY = data.frame(row.names=colnames(heatX), System=annots$shen$System)\n", "## pdf(\"ShenHighVarHeatmap.pdf\", h=8, w=8*1.3, onefile=FALSE)\n", "pheatmap(\n", "    heatX,\n", "    annotation_col = heatY,\n", "    annotation_colors = list(System = c(\n", "        'circulatory' = 'firebrick',\n", "        'digestive/excretory' = 'goldenrod',\n", "        'lymphatic' = 'lightseagreen',\n", "        'nervous' = 'darkorchid',\n", "        'other' = 'darkslategray',\n", "        'respiratory' = 'dodgerblue'\n", "    )),\n", "    show_rownames = FALSE\n", ")\n", "## garbage = dev.off()"]
}
 ],
 "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.2.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}