GMM_TI.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(monocle)\n",
    "library(mvtnorm)\n",
    "library(igraph)\n",
    "library(reshape2)\n",
    "#导入函数\n",
    "setwd(\"C:/Users/jliu25/Desktop/Unknown/functions\")\n",
    "source('functions.r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#构建cellData类型数据\n",
    "library(HSMMSingleCell)\n",
    "data(HSMM_expr_matrix) ## RPKM 矩阵,271个细胞，47192个基因\n",
    "data(HSMM_gene_annotation)\n",
    "data(HSMM_sample_sheet)\n",
    "pd <- new(\"AnnotatedDataFrame\", data = HSMM_sample_sheet)\n",
    "fd <- new(\"AnnotatedDataFrame\", data = HSMM_gene_annotation)\n",
    "HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),   \n",
    "                       phenoData = pd, \n",
    "                       featureData = fd,\n",
    "                       lowerDetectionLimit=0.1,\n",
    "                       expressionFamily=tobit(Lower=0.1))\n",
    "rpc_matrix <- relative2abs(HSMM)\n",
    "HSMM <- newCellDataSet(as(as.matrix(rpc_matrix), \"sparseMatrix\"),\n",
    "                       phenoData = pd, \n",
    "                       featureData = fd,\n",
    "                       lowerDetectionLimit=0.5,\n",
    "                       expressionFamily=negbinomial.size())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#monocle2复现\n",
    "#数据预处理\n",
    "HSMM <- estimateSizeFactors(HSMM)\n",
    "HSMM <- estimateDispersions(HSMM)\n",
    "HSMM <- detectGenes(HSMM, min_expr = 0.1)\n",
    "## 对每个基因都检查一下在多少个细胞里面是有表达量的。\n",
    "## 只留下至少在10个细胞里面有表达量的那些基因，做后续分析\n",
    "expressed_genes <- row.names(subset(fData(HSMM), num_cells_expressed >= 10))\n",
    "\n",
    "pData(HSMM)$Total_mRNAs <- Matrix::colSums(exprs(HSMM))\n",
    "HSMM <- HSMM[,pData(HSMM)$Total_mRNAs < 1e6]\n",
    "upper_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) +\n",
    "                     2*sd(log10(pData(HSMM)$Total_mRNAs)))\n",
    "lower_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) -\n",
    "                     2*sd(log10(pData(HSMM)$Total_mRNAs)))\n",
    "table(pData(HSMM)$Hours)\n",
    "qplot(Total_mRNAs, data = pData(HSMM), color = Hours, geom = \"density\") +\n",
    "  geom_vline(xintercept = lower_bound) +\n",
    "  geom_vline(xintercept = upper_bound)\n",
    "HSMM <- HSMM[,pData(HSMM)$Total_mRNAs > lower_bound & \n",
    "               pData(HSMM)$Total_mRNAs < upper_bound]                                 \n",
    "HSMM <- detectGenes(HSMM, min_expr = 0.1)\n",
    "L <- log(exprs(HSMM[expressed_genes,]))\n",
    "melted_dens_df <- melt(Matrix::t(scale(Matrix::t(L))))\n",
    "qplot(value, geom=\"density\", data=melted_dens_df) +  stat_function(fun = dnorm, size=0.5, color='red') + \n",
    "  xlab(\"Standardized log(FPKM)\") +\n",
    "  ylab(\"Density\")\n",
    "\n",
    "#细胞分类，挑出Myoblast\n",
    "## 根据基因名字找到其在表达矩阵的ID，这里是ENSEMBL数据库的ID\n",
    "MYF5_id <- row.names(subset(fData(HSMM), gene_short_name == \"MYF5\"))\n",
    "ANPEP_id <- row.names(subset(fData(HSMM), gene_short_name == \"ANPEP\"))\n",
    "## 这里选取的基因取决于自己的单细胞实验设计\n",
    "cth <- newCellTypeHierarchy()\n",
    "\n",
    "cth <- addCellType(cth, \"Myoblast\", classify_func = function(x) { x[MYF5_id,] >= 1 })\n",
    "cth <- addCellType(cth, \"Fibroblast\", classify_func = function(x)\n",
    "{ x[MYF5_id,] < 1 & x[ANPEP_id,] > 1 })\n",
    "\n",
    "HSMM <- classifyCells(HSMM, cth, 0.1)\n",
    "## 这个时候的HSMM已经被改变了，增加了属性。\n",
    "\n",
    "table(pData(HSMM)$CellType)\n",
    "\n",
    "pie <- ggplot(pData(HSMM), aes(x = factor(1), fill = factor(CellType))) +\n",
    "  geom_bar(width = 1)\n",
    "pie + coord_polar(theta = \"y\") +\n",
    "  theme(axis.title.x = element_blank(), axis.title.y = element_blank())\n",
    "\n",
    "marker_diff <- markerDiffTable(HSMM[expressed_genes,],\n",
    "                               cth,\n",
    "                               residualModelFormulaStr = \"~Media + num_genes_expressed\",\n",
    "                               cores = 1)\n",
    "candidate_clustering_genes <- row.names(subset(marker_diff, qval < 0.01))\n",
    "marker_spec <- calculateMarkerSpecificity(HSMM[candidate_clustering_genes,], cth)\n",
    "head(selectTopMarkers(marker_spec, 3))\n",
    "semisup_clustering_genes <- unique(selectTopMarkers(marker_spec, 500)$gene_id)\n",
    "HSMM <- setOrderingFilter(HSMM, semisup_clustering_genes)\n",
    "plot_ordering_genes(HSMM)\n",
    "plot_pc_variance_explained(HSMM, return_all = F)\n",
    "HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 3,\n",
    "                        norm_method = 'log',\n",
    "                        reduction_method = 'tSNE',\n",
    "                        residualModelFormulaStr = \"~Media + num_genes_expressed\",\n",
    "                        verbose = T)\n",
    "HSMM <- clusterCells(HSMM, num_clusters = 2)\n",
    "plot_cell_clusters(HSMM, 1, 2, color = \"CellType\")\n",
    "\n",
    "HSMM_myo <- HSMM\n",
    "diff_test_res <- differentialGeneTest(HSMM_myo[expressed_genes,],\n",
    "                                      fullModelFormulaStr = \"~Media\")\n",
    "ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))\n",
    "HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes)\n",
    "plot_ordering_genes(HSMM_myo)\n",
    "\n",
    "HSMM_myo <- reduceDimension(HSMM_myo, max_components = 2,\n",
    "                            method = 'DDRTree')\n",
    "HSMM_myo <- orderCells(HSMM_myo)\n",
    "plot_cell_trajectory(HSMM_myo, color_by = \"Hours\")\n",
    "plot_cell_trajectory(HSMM_myo, color_by = \"State\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HSMM_myo <- reduceDimension(HSMM, max_components=2)\n",
    "k=4\n",
    "pi = c(1/4,1/4,1/4,1/4)\n",
    "data<-t(reducedDimS(HSMM_myo))\n",
    "initialpara<-initialGMMs(data,k)\n",
    "mu<-initialpara[[\"mu\"]]\n",
    "sigma<-initialpara[[\"sig\"]]\n",
    "\n",
    "cds<-HSMM_myo\n",
    "count=1\n",
    "while (TRUE){\n",
    "    cds <- project2MST(cds, project_point_to_line_segment)\n",
    "    newdata<-cds@auxOrderingData[[\"DDRTree\"]]$pr_graph_cell_proj_dist\n",
    "    updatepara<-mulGaussCluster(t(newdata),k,pi,mu,sigma)\n",
    "    mu<-updatepara[[\"mu\"]]\n",
    "    sigma<-updatepara[[\"sigma\"]]\n",
    "    pi<-updatepara[[\"pi\"]]\n",
    "    count<-count+1\n",
    "    if(count == 100){\n",
    "      break\n",
    "    }\n",
    "}\n",
    "\n",
    "r<-updatepara[[\"r\"]]\n",
    "cluster <- which(r == apply(r,1,max),arr.ind = TRUE)\n",
    "cluster <- cluster[order(cluster[,1]),]\n",
    "#pData(cds)$State<-cluster[,2]\n",
    "pData(cds)$cluster<-as.character(cluster[,2])\n",
    "plot_cell_trajectory(cds, color_by=\"cluster\")\n",
    "\n",
    "#true state\n",
    "plot_cell_trajectory(cds, color_by=\"State\")\n",
    "plot_cell_trajectory(HSMM_myo, color_by = \"Hours\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data 2 test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#x<-read.table(\"C:/Users/jliu25/Desktop/Unknown/data/GSE60781_RPKM.txt\", header = TRUE, sep = \"\\t\", row.names=1)\n",
    "load('C:/Users/jliu25/Desktop/Unknown/data/gene_annotation.Rdata')\n",
    "load('C:/Users/jliu25/Desktop/Unknown/data/cell_annotation.Rdata')\n",
    "load('C:/Users/jliu25/Desktop/Unknown/data/expression.Rdata')\n",
    "\n",
    "#构建cellData类型数据\n",
    "pd <- new(\"AnnotatedDataFrame\", data = data.frame(c))\n",
    "fd <- new(\"AnnotatedDataFrame\", data = data.frame(g))\n",
    "cds <- newCellDataSet(as.matrix(x), phenoData = pd, featureData = fd,expressionFamily = tobit())\n",
    "\n",
    "cds <- estimateSizeFactors(cds)\n",
    "cds<- detectGenes(cds, min_expr = 0.1)\n",
    "expressed_genes <- row.names(subset(fData(cds), num_cells_expressed >= 10))\n",
    "diff_test_res <- differentialGeneTest(cds[expressed_genes,],\n",
    "+                    fullModelFormulaStr = \"~celltype\")\n",
    "ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))\n",
    "cds <- setOrderingFilter(cds, ordering_genes)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ----------以下是functions--------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mulGaussCluster <- function(data,k,pi,mu,sigma){\n",
    "  #data:数据集,向量或数据框\n",
    "  #k:聚类的个数，或高斯分布的个数\n",
    "  #pi:高斯分布的先验概率,选择各个高斯分布的概率,向量\n",
    "  #mu:高斯分布的初始均值\n",
    "  #sigma:多元高斯分布的初始协方差阵\n",
    "    N = nrow(data)\n",
    "    u0 <- mu\n",
    "    p = c()\n",
    "    for(j in 1:k){\n",
    "      prob = apply(data,1,dmvnorm,mean=mu[j,],sigma = sigma[[j]])\n",
    "      p = cbind(p,prob)\n",
    "    }\n",
    "    amatrix <- matrix(rep(pi,N),nrow=N,byrow = TRUE)\n",
    "    r = (p * amatrix) / rowSums(p * amatrix)  \n",
    "    #updata mu\n",
    "    mu <- t(r) %*% data / colSums(r) #求和可以转化为向量相乘的形式，简化计算\n",
    "    #updata sigma\n",
    "    for(j in 1:k){\n",
    "      upsigma = matrix(0,ncol(data),ncol(data))\n",
    "      for(i in 1:N){\n",
    "        upsigma = upsigma + r[i,j] * (data[i,]-u0[j,]) %*% t(data[i,]-u0[j,])  #R中向量默认为列向量\n",
    "      }\n",
    "      sigma[[j]]= upsigma/sum(r[,j]) \n",
    "    }\n",
    "    #updata pi\n",
    "    pi <- colSums(r)/N\n",
    "    return(list(mu = mu,sigma = sigma,pi = pi,r=r))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use kmeans to get initial parameters of GMMs\n",
    "initialGMMs<-function(data,k){\n",
    "  km<-kmeans(data,k)\n",
    "  cluster<-km[[\"cluster\"]]\n",
    "  mu<-matrix(0,k,ncol(data))\n",
    "  sigma<-list()\n",
    "  for (i in 1:k){\n",
    "    subdata<-data[which(cluster==i),]\n",
    "    m<-apply(subdata,2,mean)\n",
    "    s<-cov(subdata)\n",
    "    mu[i,]<-m\n",
    "    sigma[[i]]<-s\n",
    "  }\n",
    "  initialpara<-list(mu=mu,sig=sigma)\n",
    "  return (initialpara)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ----------以下未用到----------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "####################\n",
    "library(mclust) #diabetes data\n",
    "data(diabetes)\n",
    "data <- diabetes[,-1]\n",
    "k<-3\n",
    "initialpara<-initialGMMs(data,3)\n",
    "pi = c(1/3,1/3,1/3)\n",
    "mu<-initialpara[[\"mu\"]]\n",
    "sigma<-initialpara[[\"sig\"]]\n",
    "\n",
    "paras<-mulGaussCluster(data,k,pi,mu,sigma)\n",
    "centers<-ncenter(data,cluster)\n",
    "ew<-Edge_weight(centers,k)\n",
    "g <- make_graph(t(ew[,1:2]),directed = FALSE)\n",
    "m<-mst(g)\n",
    "#####################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#计算聚类中心\n",
    "ncenter<-function(data,cluster){\n",
    "    centers<-matrix(0,k,ncol(data))\n",
    "    for (i in 1:k){\n",
    "    m<-which(cluster[,2]==i)\n",
    "    centers[i，]<-apply(data[m,],2,mean)\n",
    "    }\n",
    "    return(centers)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#计算边权矩阵\n",
    "Edge_weight<-function(centers,k){\n",
    "    row<-nrow(centers)*(nrow(centers)-1)/2\n",
    "    ew<-matrix(0,row,3)\n",
    "    count<-0\n",
    "    dis<-as.matrix(dist(centers))\n",
    "    for (i in 1:(k-1)){\n",
    "        for (j in (i+1):k){\n",
    "            count<-count+1\n",
    "            ew[count,1]=i\n",
    "            ew[count,2]=j\n",
    "            ew[count,3]=dis[i,j]\n",
    "        }\n",
    "    }\n",
    "    return(ew)\n",
    "}"
   ]
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}