pagoda是2016年在nature methods發(fā)表的一個分析單細胞測序的方法交汤，這個包并不是很紅火叶沛，只是最近在學(xué)習(xí)另一個基于RNA速率推斷細胞演化的R包velocyto的時候發(fā)現(xiàn)图贸，velocyto是基于pagoda的cluster和tsne做的新症，所以這里把pagoda學(xué)習(xí)一下并做記錄罢低。
高水平的技術(shù)噪聲和對表達量的強烈依賴給主成分分析和其他降維方法帶來了困難说墨。因此，PCA的應(yīng)用以及更靈活的方法如GP-LVM或tSNE通常僅限于高表達的基因与斤。PAGODA全稱pathway and gene set overdispersion analysis肪康，在已知的重要信號通路基礎(chǔ)上對細胞進行分類，以提高統(tǒng)計效力并揭示可能的功能性解釋撩穿。例如磷支，雖然單個神經(jīng)元分化標(biāo)記物如Neurod1的表達在細胞間差異可能過于嘈雜且不確定，但在相同細胞亞群中與神經(jīng)元分化相關(guān)的許多基因的協(xié)調(diào)上調(diào)將提供區(qū)分的突出特征食寡。
接下來齐唆，我對Multipotent Peripheral Glial Cells Generate Neuroendocrine Cells of the Adrenal Medulla這篇文章中的數(shù)據(jù)用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 載入數(shù)據(jù)
統(tǒng)計數(shù)據(jù)的整體表達情況冻河，過濾掉表達基因數(shù)少于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
2 去掉重復(fù)的列名生成pagoda對象
因為我導(dǎo)入的counts文件rownames是基因名所以沒有重復(fù)，如果是ensmble號的話是有重復(fù)的叨叙，去重之后才可以生成pagoda對象锭弊，之后所有的分析結(jié)果也都會儲存在這個對象中。
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 對表達量差異很大的基因?qū)ο掠畏治鏊急戎剡M行調(diào)整
r$adjustVariance(plot=T,gam.k=10)
calculating variance fit ... using gam 1640 overdispersed genes ... 1640 persisting ... done.
4 降維聚類
pagoda中有許多可選的方法進行降維聚類擂错，這里我采用默認的也是最簡單的一種：首先用PCA降維的數(shù)據(jù)構(gòu)建一個細胞間的k近鄰稀疏矩陣味滞，即將一個細胞與它距離上最近的k個細胞聚為一類，然后在此基礎(chǔ)進行模塊優(yōu)化钮呀，找到高度連接的模塊剑鞍。最后通過層次聚類將位于同一區(qū)域內(nèi)沒有差異表達基因的cluster進一步融合，重復(fù)該過程直到?jīng)]有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圖的基礎(chǔ)上觀察某些特異基因的表達情況
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
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)設(shè)成"environment"集空屬性; 其結(jié)果不再是S4目標(biāo)對象
2: In class(object) <- "environment" :
把class(x)設(shè)成"environment"集空屬性; 其結(jié)果不再是S4目標(biāo)對象
r$plotGeneHeatmap(genes=rownames(de)[1:15],groups=r$clusters$PCA[[1]])
好了，這個包的大概功能就是這樣啦遂赠，搞懂了這個包就可以接著去搞明白velocyto了～～
