單細胞轉錄組3大R包之monocle2

阿新 • • 發佈：2022-05-03

主要是針對單細胞轉錄組測序資料開發的，用來找不同細胞型別或者不同細胞狀態的差異表達基因。分析起始是表達矩陣，作者推薦用比較老舊的Tophat+Cufflinks流程，或者RSEM, eXpress,Sailfish,等等。需要的是基於轉錄本的表達矩陣，我一般用subjunc+featureCounts 來獲取表達矩陣。

2014年版本

由Cole Trapnell 於2014年在Nature Biotechnology 雜誌發表，是一個略微複雜的R包，並給出了一個測試資料，下載地址是：

Source code
HSMM expression data

安裝方法是：

install.packages(c("VGAM", "irlba", "matrixStats", "igraph", 
"combinat", "fastICA", "grid", "ggplot2", 
"reshape2", "plyr", "parallel", "methods"))
$ R CMD INSTALL HSMMSingleCell_0.99.0.tar.gz 
$ R CMD INSTALL monocle_0.99.0.tar.gz 
source("http://bioconductor.org/biocLite.R")
biocLite()
biocLite("monocle")
library(monocle)

這一版的教程有點過時了，還用的是tophat+cufflinks組合來計算表達量，就不過多介紹了。

2017年版本

在nature methods雜誌發表的文章，更新為monocle2版本並且更換了主頁，功能也不僅僅是差異分析那麼簡單。還包括pseudotime,clustering分析，而且還可以進行基於轉錄本的差異分析，其演算法是BEAM (used in branch analysis) and Census (the core of relative2abs)，也單獨發表了文章。

用了4個公共的資料來測試說明其軟體的用法和優點。

the HSMM data set, GSE52529 (ref. 1);

the lung data set, GSE52583 (ref. 8);
the Paul et al. data set ;
the Olsson data set9, synapse ID syn4975060.

也是有著非常詳細的使用教程 , 讀取表達矩陣和分組資訊，需要理解其定義好的一些S4物件。

還提出了好幾個演算法：

dpFeature: Selecting features from dense cell clusters
Reversed graph embedding
DRTree: Dimensionality Reduction via Learning a Tree

DDRTree: discriminative dimensionality reduction via learning a tree
Census: a normalization method to convert of single-cell mRNA transcript to relative transcript counts.
BEAM : to test for branch-dependent gene expression by formulating the problem as a contrast between two negative binomial GLMs.
Branch time point detection algorithm :

S4 物件

主要是基於 CellDataSet 物件來進行下游分析，繼承自ExpressionSet物件，也是常見的3個組成：

exprs, a numeric matrix of expression values, where rows are genes, and columns are cells
phenoData, an AnnotatedDataFrame object, where rows are cells, and columns are cell attributes (such as cell type, culture condition, day captured, etc.)
featureData, an AnnotatedDataFrame object, where rows are features (e.g. genes), and columns are gene attributes, such as biotype, gc content, etc.

可以從頭建立這樣的物件，程式碼如下：

#do not run
HSMM_expr_matrix <- read.table("fpkm_matrix.txt")
HSMM_sample_sheet <- read.delim("cell_sample_sheet.txt")
HSMM_gene_annotation <- read.delim("gene_annotations.txt")
pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)
HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix), phenoData = pd, featureData = fd)

建立物件的時候需要指定引入的表達矩陣的方法，monocle2推薦用基於轉錄本的counts矩陣，同時也是預設的引數 expressionFamily=negbinomial.size() ，如果是其它RPKM/TMP等等，需要找到對應的引數。

包的用法

monocle在bioconductor官網的主頁給出了比較詳盡的測試資料的示例程式碼：

PDF
R Script

基本上花上幾個小時執行該例子，一步步理解輸入輸出，就可以學會使用。當然，要看懂演算法就比較費勁了，需要仔細讀paper。

值得一提的是最新版的monocle（version 2.4.0）依賴於 R version 3.4.0 ，如果R沒有升級，即使強行安裝了最新版monocle也是無濟於事。

install.packages('https://www.bioconductor.org/packages/release/bioc/bin/macosx/el-capitan/contrib/3.4/monocle_2.4.0.tgz',
         repos=NULL, type="source")

載入表達矩陣並轉化為CellDataSet物件
對錶達矩陣進行基於基因和樣本的過濾並可視化
無監督的聚類
pseudotime分析
差異分析

下面是實戰演練：

初識monocle

monocle在bioconductor官網的主頁給出了比較詳盡的測試資料的示例程式碼：

PDF
R Script

基本上花上幾個小時執行該例子，一步步理解輸入輸出，就可以學會使用。當然，要看懂演算法就比較費勁了，需要仔細讀paper。

安裝並且載入包和測試資料

如果還沒安裝，就執行：

source("http://bioconductor.org/biocLite.R")
biocLite()
biocLite("monocle")
biocLite("HSMMSingleCell")

如果已經安裝，請直接載入

library(Biobase)
library(knitr)
library(reshape2)
library(ggplot2)
library(HSMMSingleCell)
library(monocle)
data(HSMM_expr_matrix) ## RPKM 矩陣,271個細胞，47192個基因
data(HSMM_gene_annotation)
data(HSMM_sample_sheet)
HSMM_expr_matrix[1:10,1:5]

##                     T0_CT_A01 T0_CT_A03 T0_CT_A05 T0_CT_A06 T0_CT_A07
## ENSG00000000003.10  21.984400  1.280040 43.461800   0.00000 39.807600
## ENSG00000000005.5    0.000000  0.000000  0.000000   0.00000  0.000000
## ENSG00000000419.8   40.059700 77.580800  6.496560   4.90934  1.156520
## ENSG00000000457.8    0.937081  0.729195  0.000000   0.00000  0.000000
## ENSG00000000460.12   0.740922 57.578500  3.935870   0.00000  0.000000
## ENSG00000000938.8    0.000000  0.000000  0.000000   0.00000  0.000000
## ENSG00000000971.11   3.002980 15.302400 50.804800   4.68513  0.000000
## ENSG00000001036.8  128.197000 16.086700 25.320900  10.66480 63.773500
## ENSG00000001084.6    7.619720  0.000000  0.000000   0.00000  0.000000
## ENSG00000001167.10  13.024900 24.777600  0.681409   1.36587  0.399352

head(HSMM_gene_annotation)

##                    gene_short_name        biotype num_cells_expressed
## ENSG00000000003.10          TSPAN6 protein_coding                 231
## ENSG00000000005.5             TNMD protein_coding                   0
## ENSG00000000419.8             DPM1 protein_coding                 275
## ENSG00000000457.8            SCYL3 protein_coding                  24
## ENSG00000000460.12        C1orf112 protein_coding                  78
## ENSG00000000938.8              FGR protein_coding                   0
##                    use_for_ordering
## ENSG00000000003.10            FALSE
## ENSG00000000005.5             FALSE
## ENSG00000000419.8             FALSE
## ENSG00000000457.8             FALSE
## ENSG00000000460.12             TRUE
## ENSG00000000938.8             FALSE

head(HSMM_sample_sheet)

##                Library Well Hours Media Mapped.Fragments Pseudotime State
## T0_CT_A01 SCC10013_A01  A01     0    GM          1958074  23.916673     1
## T0_CT_A03 SCC10013_A03  A03     0    GM          1930722   9.022265     1
## T0_CT_A05 SCC10013_A05  A05     0    GM          1452623   7.546608     1
## T0_CT_A06 SCC10013_A06  A06     0    GM          2566325  21.463948     1
## T0_CT_A07 SCC10013_A07  A07     0    GM          2383438  11.299806     1
## T0_CT_A08 SCC10013_A08  A08     0    GM          1472238  67.436042     2

構建S4物件，CellDataSet

主要是讀取表達矩陣和樣本描述資訊，這裡介紹兩種方式，一種是讀取基於 subjunc+featureCounts 分析後的reads counts矩陣，一種是讀取 tophat+cufflinks 得到的RPKM表達矩陣。

讀取上游分析的輸出檔案

library(monocle)
library(scater, quietly = TRUE)
library(knitr)
options(stringsAsFactors = FALSE)

# 這個檔案是表達矩陣，包括線粒體基因和 ERCC spike-ins 的表達量，可以用來做質控
molecules <- read.table("tung/molecules.txt", sep = "t")

## 這個檔案是表達矩陣涉及到的所有樣本的描述資訊，包括樣本來源於哪個細胞，以及哪個批次。
anno <- read.table("tung/annotation.txt", sep = "t", header = TRUE)
rownames(anno)=colnames(molecules)
library(org.Hs.eg.db)
eg2symbol=toTable(org.Hs.egSYMBOL)
eg2ensembl=toTable(org.Hs.egENSEMBL)
egid=eg2ensembl[ match(rownames(molecules),eg2ensembl$ensembl_id),'gene_id']
symbol=eg2symbol[match( egid ,eg2symbol$gene_id),'symbol']
gene_annotation = data.frame(ensembl=rownames(molecules),
                             gene_short_name=symbol,
                             egid=egid)
rownames(gene_annotation)=rownames(molecules)

pd <- new("AnnotatedDataFrame", data = anno)
fd <- new("AnnotatedDataFrame", data = gene_annotation)
#tung <- newCellDataSet(as.matrix(molecules), phenoData = pd, featureData = fd)
tung <- newCellDataSet(as(as.matrix(molecules), "sparseMatrix"),
                       phenoData = pd, 
                       featureData = fd,
                       lowerDetectionLimit=0.5,
                       expressionFamily=negbinomial.size())

tung

## CellDataSet (storageMode: environment)
## assayData: 19027 features, 864 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: NA19098.r1.A01 NA19098.r1.A02 ... NA19239.r3.H12
##     (864 total)
##   varLabels: individual replicate ... Size_Factor (6 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: ENSG00000237683 ENSG00000187634 ... ERCC-00171
##     (19027 total)
##   fvarLabels: ensembl gene_short_name egid
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:

可以看到物件已經構造成功，是一個包含了 19027 features, 864 samples 的表達矩陣，需要進行一系列的過濾之後，拿到高質量的單細胞轉錄組資料進行下游分析。

這些樣本來源於3個不同的人，每個人有3個批次的單細胞，每個批次單細胞都是96個。

或者使用內建資料個構建S4物件

pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)

# First create a CellDataSet from the relative expression levels

## 這裡僅僅是針對rpkm表達矩陣的讀取
HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),   
                       phenoData = pd, 
                       featureData = fd,
                       lowerDetectionLimit=0.1,
                       expressionFamily=tobit(Lower=0.1))

# Next, use it to estimate RNA counts
rpc_matrix <- relative2abs(HSMM)
rpc_matrix[1:10,1:5]

##                     T0_CT_A01  T0_CT_A03  T0_CT_A05  T0_CT_A06  T0_CT_A07
## ENSG00000000003.10 1.60309506 0.09929705 2.93679928 0.00000000 2.18692386
## ENSG00000000005.5  0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## ENSG00000000419.8  2.92113986 6.01820615 0.43898533 0.34343867 0.06353614
## ENSG00000000457.8  0.06833163 0.05656613 0.00000000 0.00000000 0.00000000
## ENSG00000000460.12 0.05402778 4.46655980 0.26595447 0.00000000 0.00000000
## ENSG00000000938.8  0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## ENSG00000000971.11 0.21897629 1.18705914 3.43298023 0.32775379 0.00000000
## ENSG00000001036.8  9.34808217 1.24789995 1.71098300 0.74606865 3.50354678
## ENSG00000001084.6  0.55562742 0.00000000 0.00000000 0.00000000 0.00000000
## ENSG00000001167.10 0.94977133 1.92208258 0.04604415 0.09555105 0.02193934

## rpkm格式的表達值需要轉換成reads counts之後才可以進行下游分析！

# Now, make a new CellDataSet using the RNA counts
HSMM <- newCellDataSet(as(as.matrix(rpc_matrix), "sparseMatrix"),
                       phenoData = pd, 
                       featureData = fd,
                       lowerDetectionLimit=0.5,
                       expressionFamily=negbinomial.size())

下面的分析，都基於內建資料構建的S4物件，HSMM

過濾低質量細胞和未檢測到的基因

基於基因的過濾

這裡只是把基因挑選出來，並沒有對S4物件進行過濾操作。這個 detectGenes 函式還計算了每個細胞裡面表達的基因數量。

HSMM <- estimateSizeFactors(HSMM)
HSMM <- estimateDispersions(HSMM)

## Warning: Deprecated, use tibble::rownames_to_column() instead.

## Removing 139 outliers

HSMM <- detectGenes(HSMM, min_expr = 0.1)
print(head(fData(HSMM)))

##                    gene_short_name        biotype num_cells_expressed
## ENSG00000000003.10          TSPAN6 protein_coding                 184
## ENSG00000000005.5             TNMD protein_coding                   0
## ENSG00000000419.8             DPM1 protein_coding                 211
## ENSG00000000457.8            SCYL3 protein_coding                  18
## ENSG00000000460.12        C1orf112 protein_coding                  47
## ENSG00000000938.8              FGR protein_coding                   0
##                    use_for_ordering
## ENSG00000000003.10            FALSE
## ENSG00000000005.5             FALSE
## ENSG00000000419.8             FALSE
## ENSG00000000457.8             FALSE
## ENSG00000000460.12             TRUE
## ENSG00000000938.8             FALSE

## 對每個基因都檢查一下在多少個細胞裡面是有表達量的。
## 只留下至少在10個細胞裡面有表達量的那些基因，做後續分析
expressed_genes <- row.names(subset(fData(HSMM), num_cells_expressed >= 10))
length(expressed_genes) ## 只剩下了14224個基因

## [1] 14224

print(head(pData(HSMM)))

##                Library Well Hours Media Mapped.Fragments Pseudotime State
## T0_CT_A01 SCC10013_A01  A01     0    GM          1958074  23.916673     1
## T0_CT_A03 SCC10013_A03  A03     0    GM          1930722   9.022265     1
## T0_CT_A05 SCC10013_A05  A05     0    GM          1452623   7.546608     1
## T0_CT_A06 SCC10013_A06  A06     0    GM          2566325  21.463948     1
## T0_CT_A07 SCC10013_A07  A07     0    GM          2383438  11.299806     1
## T0_CT_A08 SCC10013_A08  A08     0    GM          1472238  67.436042     2
##           Size_Factor num_genes_expressed
## T0_CT_A01    1.392811                6850
## T0_CT_A03    1.311607                6947
## T0_CT_A05    1.218922                7019
## T0_CT_A06    1.013981                5560
## T0_CT_A07    1.085580                5998
## T0_CT_A08    1.099878                6055

基於樣本表達量進行過濾

這裡選擇的是通過不同時間點取樣的細胞來進行分組檢視，把超過2個sd 的那些樣本的臨界值挑選出來，下一步過濾的時候使用。

pData(HSMM)$Total_mRNAs <- Matrix::colSums(exprs(HSMM))
HSMM <- HSMM[,pData(HSMM)$Total_mRNAs < 1e6]
upper_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) +
                     2*sd(log10(pData(HSMM)$Total_mRNAs)))
lower_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) -
                     2*sd(log10(pData(HSMM)$Total_mRNAs)))
table(pData(HSMM)$Hours)

## 
##  0 24 48 72 
## 69 74 79 49

qplot(Total_mRNAs, data = pData(HSMM), color = Hours, geom = "density") +
  geom_vline(xintercept = lower_bound) +
  geom_vline(xintercept = upper_bound)

執行過濾並可視化檢查一下

上面已經根據基因表達情況以及樣本的總測序資料選擇好了閾值，下面就可以視覺化並且對比檢驗一下執行過濾與否的區別。

HSMM <- HSMM[,pData(HSMM)$Total_mRNAs > lower_bound & 
               pData(HSMM)$Total_mRNAs < upper_bound]                                 
HSMM <- detectGenes(HSMM, min_expr = 0.1)

L <- log(exprs(HSMM[expressed_genes,]))

melted_dens_df <- melt(Matrix::t(scale(Matrix::t(L))))

qplot(value, geom="density", data=melted_dens_df) +  stat_function(fun = dnorm, size=0.5, color='red') + 
  xlab("Standardized log(FPKM)") +
  ylab("Density")

聚類

根據指定基因對單細胞轉錄組表達矩陣進行分類

下面這個程式碼只適用於這個測試資料，主要是生物學背景知識，用MYF5基因和ANPEP基因來對細胞進行分類，可以區分Myoblast和Fibroblast。如果是自己的資料，建議多讀讀paper看看如何選取合適的基因，或者乾脆跳過這個程式碼。

## 根據基因名字找到其在表達矩陣的ID，這裡是ENSEMBL資料庫的ID
MYF5_id <- row.names(subset(fData(HSMM), gene_short_name == "MYF5"))
ANPEP_id <- row.names(subset(fData(HSMM), gene_short_name == "ANPEP"))
## 這裡選取的基因取決於自己的單細胞實驗設計
cth <- newCellTypeHierarchy()

cth <- addCellType(cth, "Myoblast", classify_func = function(x) { x[MYF5_id,] >= 1 })
cth <- addCellType(cth, "Fibroblast", classify_func = function(x)
{ x[MYF5_id,] < 1 & x[ANPEP_id,] > 1 })

HSMM <- classifyCells(HSMM, cth, 0.1)

## Warning: Deprecated, use tibble::rownames_to_column() instead.

## Warning: Deprecated, use tibble::rownames_to_column() instead.

## 這個時候的HSMM已經被改變了，增加了屬性。

table(pData(HSMM)$CellType)

## 
## Fibroblast   Myoblast    Unknown 
##         60         87        124

pie <- ggplot(pData(HSMM), aes(x = factor(1), fill = factor(CellType))) +
  geom_bar(width = 1)
pie + coord_polar(theta = "y") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank())

可以看到還有很大一部分細胞僅僅是根據這兩個基因的表達量是無法成功的歸類的。這個是很正常的，因為單細胞轉錄組測序裡面的mRNA捕獲率不夠好。通過這個步驟成功的給HSMM這個S4物件增加了一個屬性，就是CellType，在下面的分析中會用得著。

無監督聚類

這裡需要安裝最新版R包才可以使用裡面的一些函式，因為上面的步驟基於指定基因的表達量進行細胞分組會漏掉很多資訊，所以需要更好的聚類方式。

disp_table <- dispersionTable(HSMM)
head(disp_table)

##              gene_id mean_expression dispersion_fit dispersion_empirical
## 1 ENSG00000000003.10      1.80534418       1.249323             1.215666
## 2  ENSG00000000419.8      2.17342979       1.099130             1.008759
## 3  ENSG00000000457.8      0.02518587      63.932303            23.177101
## 4 ENSG00000000460.12      0.15331486      10.805439            17.941440
## 5 ENSG00000000971.11      2.45231977       1.015354             1.287973
## 6  ENSG00000001036.8      1.04484075       1.894827             1.540376

## 只有滿足 條件的10198個基因才能進入聚類分析
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
HSMM <- setOrderingFilter(HSMM, unsup_clustering_genes$gene_id)
plot_ordering_genes(HSMM)

## 這裡看看基因的表達量和基因的變異度之間的關係
## 處在灰色陰影區域的基因會被拋棄掉，不進入聚類分析。
plot_pc_variance_explained(HSMM, return_all = F) # norm_method = 'log',

HSMM <- reduceDimension(HSMM, max_components=2, num_dim = 6, 
                        reduction_method = 'tSNE', verbose = T) 
HSMM <- clusterCells(HSMM, num_clusters=2)

## Distance cutoff calculated to 1.072748

## 這裡先用tSNE的聚類方法處理HSMM資料集，並可視化展示
plot_cell_clusters(HSMM, 1, 2, color="CellType", markers=c("MYF5", "ANPEP"))

## 可以看到並不能把細胞型別完全區分開，這個是完全有可能的，因為雖然是同一種細胞，但是有著不同的培養條件。
head(pData(HSMM))

##                Library Well Hours Media Mapped.Fragments Pseudotime State
## T0_CT_A01 SCC10013_A01  A01     0    GM          1958074  23.916673     1
## T0_CT_A03 SCC10013_A03  A03     0    GM          1930722   9.022265     1
## T0_CT_A05 SCC10013_A05  A05     0    GM          1452623   7.546608     1
## T0_CT_A06 SCC10013_A06  A06     0    GM          2566325  21.463948     1
## T0_CT_A07 SCC10013_A07  A07     0    GM          2383438  11.299806     1
## T0_CT_A08 SCC10013_A08  A08     0    GM          1472238  67.436042     2
##           Size_Factor num_genes_expressed Total_mRNAs CellType Cluster
## T0_CT_A01    1.392811                6850       39080 Myoblast       2
## T0_CT_A03    1.311607                6947       36720 Myoblast       1
## T0_CT_A05    1.218922                7019       34112 Myoblast       1
## T0_CT_A06    1.013981                5560       28384 Myoblast       2
## T0_CT_A07    1.085580                5998       30360 Myoblast       1
## T0_CT_A08    1.099878                6055       30808  Unknown       2
##           peaks  halo     delta      rho
## T0_CT_A01 FALSE FALSE 1.0694920 1.146961
## T0_CT_A03 FALSE FALSE 0.5544267 2.744092
## T0_CT_A05 FALSE FALSE 0.3270436 4.479191
## T0_CT_A06 FALSE FALSE 0.4767768 2.416054
## T0_CT_A07 FALSE FALSE 0.6011590 2.593689
## T0_CT_A08 FALSE FALSE 1.2702897 2.395104

head(fData(HSMM))

##                    gene_short_name        biotype num_cells_expressed
## ENSG00000000003.10          TSPAN6 protein_coding                 184
## ENSG00000000005.5             TNMD protein_coding                   0
## ENSG00000000419.8             DPM1 protein_coding                 211
## ENSG00000000457.8            SCYL3 protein_coding                  18
## ENSG00000000460.12        C1orf112 protein_coding                  47
## ENSG00000000938.8              FGR protein_coding                   0
##                    use_for_ordering
## ENSG00000000003.10             TRUE
## ENSG00000000005.5             FALSE
## ENSG00000000419.8              TRUE
## ENSG00000000457.8             FALSE
## ENSG00000000460.12             TRUE
## ENSG00000000938.8             FALSE

## 所以這裡也區分一下 培養基， a high-mitogen growth medium (GM) to a low-mitogen differentiation medium (DM). 
plot_cell_clusters(HSMM, 1, 2, color="Media")

## 因為我們假設就2種細胞型別，所以在做聚類的時候可以把這個引數新增進去，這樣可以去除無關變數的干擾。
HSMM <- reduceDimension(HSMM, max_components=2, num_dim = 2, reduction_method = 'tSNE', 
                        residualModelFormulaStr="~Media + num_genes_expressed", verbose = T) #
HSMM <- clusterCells(HSMM, num_clusters=2)

## Distance cutoff calculated to 1.284778

plot_cell_clusters(HSMM, 1, 2, color="CellType")

plot_cell_clusters(HSMM, 1, 2, color="Cluster") + facet_wrap(~CellType)

半監督聚類

## 這裡的差異分析非常耗時

marker_diff <- markerDiffTable(HSMM[expressed_genes,], 
                               cth, 
                               residualModelFormulaStr="~Media + num_genes_expressed",
                               cores=1)
head(marker_diff)

##                    status           family      pval      qval
## ENSG00000000003.10     OK negbinomial.size 0.8548230 1.0000000
## ENSG00000000419.8      OK negbinomial.size 0.9329316 1.0000000
## ENSG00000000457.8      OK negbinomial.size 0.7176166 0.9954975
## ENSG00000000460.12     OK negbinomial.size 0.2700496 0.8250088
## ENSG00000000971.11     OK negbinomial.size 0.4489895 0.9171190
## ENSG00000001036.8      OK negbinomial.size 0.5731998 0.9524046
##                    gene_short_name        biotype num_cells_expressed
## ENSG00000000003.10          TSPAN6 protein_coding                 184
## ENSG00000000419.8             DPM1 protein_coding                 211
## ENSG00000000457.8            SCYL3 protein_coding                  18
## ENSG00000000460.12        C1orf112 protein_coding                  47
## ENSG00000000971.11             CFH protein_coding                 198
## ENSG00000001036.8            FUCA2 protein_coding                 171
##                    use_for_ordering
## ENSG00000000003.10             TRUE
## ENSG00000000419.8              TRUE
## ENSG00000000457.8             FALSE
## ENSG00000000460.12             TRUE
## ENSG00000000971.11             TRUE
## ENSG00000001036.8              TRUE

## 就是對每個基因增加了pval和qval兩列資訊，挑選出那些在不同media培養條件下顯著差異表達的基因，310個，
candidate_clustering_genes <- row.names(subset(marker_diff, qval < 0.01))

## 計算這310個基因在不同的celltype的specificity值
marker_spec <- calculateMarkerSpecificity(HSMM[candidate_clustering_genes,], cth)
head(selectTopMarkers(marker_spec, 3))

##              gene_id   CellType specificity
## 1 ENSG00000019991.11 Fibroblast   0.9892130
## 2 ENSG00000128340.10 Fibroblast   0.9999602
## 3  ENSG00000163710.3 Fibroblast   0.9729971
## 4  ENSG00000111049.3   Myoblast   0.9743099
## 5  ENSG00000239922.1   Myoblast   0.9719681
## 6  ENSG00000270123.1   Myoblast   1.0000000

semisup_clustering_genes <- unique(selectTopMarkers(marker_spec, 500)$gene_id)
HSMM <- setOrderingFilter(HSMM, semisup_clustering_genes)
plot_ordering_genes(HSMM)

## 重新挑選基因，只用黑色高亮的基因來進行聚類。

plot_pc_variance_explained(HSMM, return_all = F) # norm_method = 'log',

HSMM <- reduceDimension(HSMM, max_components=2, num_dim = 2, reduction_method = 'tSNE', 
                        residualModelFormulaStr="~Media + num_genes_expressed", verbose = T) 
HSMM <- clusterCells(HSMM, num_clusters=2)

## Distance cutoff calculated to 1.02776

plot_cell_clusters(HSMM, 1, 2, color="CellType")

HSMM <- clusterCells(HSMM,
                     num_clusters=2, 
                     frequency_thresh=0.1,
                     cell_type_hierarchy=cth)

## Distance cutoff calculated to 1.02776

plot_cell_clusters(HSMM, 1, 2, color="CellType", markers = c("MYF5", "ANPEP"))

pie <- ggplot(pData(HSMM), aes(x = factor(1), fill = factor(CellType))) +
  geom_bar(width = 1)
pie + coord_polar(theta = "y") + 
  theme(axis.title.x=element_blank(), axis.title.y=element_blank())

Pseudotime分析

主要目的是：Constructing Single Cell Trajectories

發育過程中細胞狀態是不斷變化的，monocle包利用演算法學習所有基因的表達模式來把每個細胞安排到各各自的發展軌跡。在大多數生物學過程中，參與的細胞通常不是同步發展的，只有單細胞轉錄組技術才能把處於該過程中各個中間狀態的細胞分離開來，而monocle包裡面的pseudotime分析方法正是要探究這些。

choose genes that define a cell’s progress
reduce data dimensionality
order cells along the trajectory

其中第一個步驟挑選合適的基因有3種策略，分別是：

Ordering based on genes that differ between clusters
Selecting genes with high dispersion across cells
Ordering cells using known marker genes

無監督的Pseudotime分析

HSMM_myo <- HSMM[,pData(HSMM)$CellType == "Myoblast"]   
HSMM_myo <- estimateDispersions(HSMM_myo)

## Warning: Deprecated, use tibble::rownames_to_column() instead.

## Removing 143 outliers

## 策略1：  Ordering based on genes that differ between clusters
if(F){
  diff_test_res <- differentialGeneTest(HSMM_myo[expressed_genes,],
                                      fullModelFormulaStr="~Media")
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
}
## 策略2：Selecting genes with high dispersion across cells
disp_table <- dispersionTable(HSMM_myo)
ordering_genes <- subset(disp_table, 
                         mean_expression >= 0.5 & 
                           dispersion_empirical >= 1 * dispersion_fit)$gene_id

HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes)
plot_ordering_genes(HSMM_myo)

## Warning: Transformation introduced infinite values in continuous y-axis

## 挑選變異度大的基因，如圖所示

HSMM_myo <- reduceDimension(HSMM_myo, max_components=2)
HSMM_myo <- orderCells(HSMM_myo)

## 排序好的細胞可以直接按照發育順序視覺化
plot_cell_trajectory(HSMM_myo, color_by="State")

直接做差異分析

前面的聚類分析和Pseudotime分析都需要取基因子集，就已經利用過差異分析方法來挑選那些有著顯著表達差異的基因。如果對所有的基因來檢驗，非常耗時。

marker_genes <- row.names(subset(fData(HSMM_myo), 
                                 gene_short_name %in% c("MEF2C", "MEF2D", "MYF5", 
                                                        "ANPEP", "PDGFRA","MYOG", 
                                                        "TPM1",  "TPM2",  "MYH2", 
                                                        "MYH3",  "NCAM1", "TNNT1", 
                                                        "TNNT2", "TNNC1", "CDK1", 
                                                        "CDK2",  "CCNB1", "CCNB2", 
                                                        "CCND1", "CCNA1", "ID1")))

diff_test_res <- differentialGeneTest(HSMM_myo[marker_genes,], 
                                      fullModelFormulaStr="~Media")
# Select genes that are significant at an FDR < 10%
sig_genes <- subset(diff_test_res, qval < 0.1)
sig_genes[,c("gene_short_name", "pval", "qval")]

##                    gene_short_name         pval         qval
## ENSG00000081189.9            MEF2C 8.463396e-20 4.443283e-19
## ENSG00000105048.12           TNNT1 3.017738e-12 7.921562e-12
## ENSG00000109063.9             MYH3 4.105825e-33 4.311116e-32
## ENSG00000111049.3             MYF5 1.300906e-30 9.106344e-30
## ENSG00000114854.3            TNNC1 1.721612e-18 7.230769e-18
## ENSG00000118194.14           TNNT2 2.232213e-37 4.687647e-36
## ENSG00000122180.4             MYOG 2.532610e-12 7.597830e-12
## ENSG00000123374.6             CDK2 3.017043e-02 3.959868e-02
## ENSG00000125414.13            MYH2 6.221763e-06 1.005054e-05
## ENSG00000125968.7              ID1 1.734006e-05 2.601009e-05
## ENSG00000134057.10           CCNB1 4.502654e-11 1.050619e-10
## ENSG00000140416.15            TPM1 9.914869e-08 1.892839e-07
## ENSG00000149294.12           NCAM1 2.473279e-18 8.656478e-18
## ENSG00000157456.3            CCNB2 1.529020e-07 2.675785e-07
## ENSG00000170312.11            CDK1 5.316306e-08 1.116424e-07
## ENSG00000198467.8             TPM2 9.205156e-04 1.288722e-03

## 可以看到挑選的都是顯著差異表達的基因。

還可以挑選其中幾個基因來視覺化看看它們是如何在不同組差異表達的。這個畫圖函式自己都可以寫。

MYOG_ID1 <- HSMM_myo[row.names(subset(fData(HSMM_myo), 
                                      gene_short_name %in% c("MYOG", "CCNB2"))),]
plot_genes_jitter(MYOG_ID1, grouping="Media", ncol=2)

這樣就可以測試某些基因，是否能區分細胞群體的不同型別及狀態

to_be_tested <- row.names(subset(fData(HSMM), 
                                 gene_short_name %in% c("UBC", "NCAM1", "ANPEP"))) 
cds_subset <- HSMM[to_be_tested,]

diff_test_res <- differentialGeneTest(cds_subset, fullModelFormulaStr="~CellType")
diff_test_res[,c("gene_short_name", "pval", "qval")]

##                    gene_short_name         pval         qval
## ENSG00000149294.12           NCAM1 2.853848e-92 8.561545e-92
## ENSG00000150991.10             UBC 2.852264e-01 2.852264e-01
## ENSG00000166825.9            ANPEP 4.723193e-15 7.084790e-15

plot_genes_jitter(cds_subset, grouping="CellType", color_by="CellType", 
                  nrow=1, ncol=NULL, plot_trend=TRUE)

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

full_model_fits <- fitModel(cds_subset, modelFormulaStr="~CellType")
reduced_model_fits <- fitModel(cds_subset, modelFormulaStr="~1")
diff_test_res <- compareModels(full_model_fits, reduced_model_fits)
diff_test_res

##                    status           family         pval         qval
## ENSG00000149294.12     OK negbinomial.size 2.853848e-92 8.561545e-92
## ENSG00000150991.10     OK negbinomial.size 2.852264e-01 2.852264e-01
## ENSG00000166825.9      OK negbinomial.size 4.723193e-15 7.084790e-15

plot_genes_in_pseudotime(cds_subset, color_by="Hours")

演算法

dpFeature: Selecting features from dense cell clusters
Reversed graph embedding
DRTree: Dimensionality Reduction via Learning a Tree
DDRTree: discriminative dimensionality reduction via learning a tree
Census: a normalization method to convert of single-cell mRNA transcript to relative transcript counts.
BEAM : to test for branch-dependent gene expression by formulating the problem as a contrast between two negative binomial GLMs.
Branch time point detection algorithm :

演算法講起來，就複雜了，略過。