单细胞测序R包PAGODA使用

pagoda是2016年在nature methods发表的一个分析单细胞测序的方法，这个包并不是很红火，只是最近在学习另一个基于RNA速率推断细胞演化的R包velocyto的时候发现，velocyto是基于pagoda的cluster和tsne做的，所以这里把pagoda学习一下并做记录。

高水平的技术噪声和对表达量的强烈依赖给主成分分析和其他降维方法带来了困难。因此，PCA的应用以及更灵活的方法如GP-LVM或tSNE通常仅限于高表达的基因。PAGODA全称pathway and gene set overdispersion analysis，在已知的重要信号通路基础上对细胞进行分类，以提高统计效力并揭示可能的功能性解释。例如，虽然单个神经元分化标记物如Neurod1的表达在细胞间差异可能过于嘈杂且不确定，但在相同细胞亚群中与神经元分化相关的许多基因的协调上调将提供区分的突出特征。

接下来，我对Multipotent Peripheral Glial Cells Generate Neuroendocrine Cells of the Adrenal Medulla这篇文章中的数据用pagoda进行分析。

下载pagoda

先在linux下安装依赖

sudo apt-get update
sudo apt-get -y install build-essential cmake gsl-bin libgsl0-dev libeigen3-dev libboost-all-dev libssl-dev libcurl4-openssl-dev libssl-dev libcairo2-dev libxt-dev libgtk2.0-dev libcairo2-dev xvfb xauth xfonts-base

接着在R中安装

# Install devtools
install.packages("devtools")
# Install Bioconductor dependencies
source("http://bioconductor.org/biocLite.R")
biocLite(c("GO.db", "org.Hs.eg.db","org.Mm.eg.db", "pcaMethods"), suppressUpdates=TRUE)
library(devtools)
install_github("igraph/rigraph") # Don't install with install.packages()
install_github("jkrijthe/Rtsne",ref="openmp")
install.packages(c("Cairo","urltools"))

# Install pagoda
install_github("hms-dbmi/pagoda2")
library('pagoda2')
# Pagoda2 is now ready to use

使用pagoda

1 载入数据

统计数据的整体表达情况，过滤掉表达基因数少于3000个的细胞以及在少于5个细胞中表达的基因,可以看到counts由23420行384列减少为16701行384列。

#Loading the libraries
library(Matrix)
library(pagoda2)
library(igraph)
#Loading the dataset
counts=read.table('E12.5_counts.txt',header = T,sep = '\t')
par(mfrow=c(1,2), mar = c(3.5,3.5,2.0,0.5), mgp = c(2,0.65,0), cex = 1.0)
hist(log10(colSums(cm)+1),main='molecules per cell',col='cornsilk',xlab='log10(molecules per cell)')
hist(log10(rowSums(cm)+1),main='molecules per gene',col='cornsilk',xlab='log10(molecules per gene])')
dim(counts)
[1] 23420   384
counts <- counts[rowSums(counts)>=5,colSums(counts)>=3000]
dim(counts)
[1] 16701   384

distribution.png

2 去掉重复的列名生成pagoda对象

因为我导入的counts文件rownames是基因名所以没有重复，如果是ensmble号的话是有重复的，去重之后才可以生成pagoda对象，之后所有的分析结果也都会储存在这个对象中。

rownames(counts) <- make.unique(rownames(counts))
#generate a pagoda object
r <- Pagoda2$new(as.matrix(counts),log.scale=TRUE, n.cores=2)
384 cells, 16701 genes; normalizing ... using plain model log scale ... done.

3 对表达量差异很大的基因对下游分析所占比重进行调整

r$adjustVariance(plot=T,gam.k=10)
calculating variance fit ... using gam 1640 overdispersed genes ... 1640 persisting ... done.

Rplot.png

4 降维聚类

pagoda中有许多可选的方法进行降维聚类，这里我采用默认的也是最简单的一种：首先用PCA降维的数据构建一个细胞间的k近邻稀疏矩阵，即将一个细胞与它距离上最近的k个细胞聚为一类，然后在此基础进行模块优化，找到高度连接的模块。最后通过层次聚类将位于同一区域内没有差异表达基因的cluster进一步融合，重复该过程直到没有clusters可以合并。可以看到，细胞被分为了8个cluster。

#First, the PCA reduction.
r$calculatePcaReduction(nPcs=50,n.odgenes=3e3)
running PCA using 3000 OD genes .... done
#generate a KNN graph
r$makeKnnGraph(k=35,type='PCA',center=T,distance='cosine')
creating space of type angular done
adding data ... done
building index ... done
querying ... done
#call clusters
r$getKnnClusters(method=infomap.community,type='PCA')
#tsne
r$getEmbedding(type='PCA',embeddingType='tSNE',perplexity=50,verbose=F)
calculating distance ... pearson ...running tSNE using 2 cores:
r$plotEmbedding(type='PCA',embeddingType='tSNE',show.legend=F,mark.clusters=T,min.group.size=20,shuffle.colors=F,mark.cluster.cex=1,alpha=0.8,main='clusters (tSNE)')

tsne.png

在tsne图的基础上观察某些特异基因的表达情况

par(mfrow=c(1,2))
gene <-"Chga"
r$plotEmbedding(type='PCA',embeddingType='tSNE',colors=r$counts[,gene],shuffle.colors=F,mark.cluster.cex=1,alpha=0.8,main=gene)
treating colors as a gradient with zlim: 0 1.61094 
gene <-"Serpine2"
r$plotEmbedding(type='PCA',embeddingType='tSNE',colors=r$counts[,gene],shuffle.colors=F,mark.cluster.cex=1,alpha=0.8,main=gene)
treating colors as a gradient with zlim: 0 1.61094

chga_serpine2.png

5 差异基因

计算每个cluster的差异基因并可视化，比如画出cluster2中高表达基因的热图

r$getDifferentialGenes(type='PCA',verbose=T,clusterType='community')
running differential expression with  8  clusters ... adjusting p-values ... done.
#visualise the top markers of a specific cluster
de=r$diffgenes$PCA$community$`2`
Warning messages:
1: In class(object) <- "environment" :
  把class(x)设成"environment"集空属性; 其结果不再是S4目标对象
2: In class(object) <- "environment" :
  把class(x)设成"environment"集空属性; 其结果不再是S4目标对象
r$plotGeneHeatmap(genes=rownames(de)[1:15],groups=r$clusters$PCA[[1]])

heatmap.png

好了，这个包的大概功能就是这样啦，搞懂了这个包就可以接着去搞明白velocyto了～～