劉小澤寫于19.9.5-第三單元第十一講：使用monocle2分析文章數(shù)據(jù)
筆記目的：根據(jù)生信技能樹的單細胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
載入數(shù)據(jù)，創(chuàng)建對象
rm(list = ls())
Sys.setenv(R_MAX_NUM_DLLS=999)
## 首先載入文章的數(shù)據(jù)
load(file='../input.Rdata')
# 原始表達矩陣
counts=a
counts[1:4,1:4];dim(counts) # 2萬多個基因，768個細胞(需要下一步進行過濾)
library(stringr)
# 樣本信息
meta=df
> head(meta)
g plate n_g all
SS2_15_0048_A3 1 0048 2624 all
SS2_15_0048_A6 1 0048 2664 all
SS2_15_0048_A5 2 0048 3319 all
SS2_15_0048_A4 3 0048 4447 all
SS2_15_0048_A1 2 0048 4725 all
SS2_15_0048_A2 3 0048 5263 all
# 按行/基因檢查：每個基因在多少細胞中有表達量
fivenum(apply(counts,1,function(x) sum(x>0) ))
boxplot(apply(counts,1,function(x) sum(x>0) ))
# 按列/樣本檢查：每個細胞中存在多少表達的基因
fivenum(apply(counts,2,function(x) sum(x>0) ))
hist(apply(counts,2,function(x) sum(x>0) ))
按行+按列檢查結(jié)果
左圖顯示了：有50%的基因只在低于20個細胞中有表達量，還有許多沒有表達量的：
> table(apply(counts,1,function(x) sum(x>0) )==0)
FALSE TRUE
17429 7153
# 存在7000多個基因在任何一個細胞中都沒表達
右圖顯示了：大部分細胞都包含2000個以上的有表達的基因
開始使用newCellDataSet()構(gòu)建monocle對象：
# ---------首先構(gòu)建基因的注釋信息(feature_data)-----------
gene_ann <- data.frame(
gene_short_name = row.names(counts),
row.names = row.names(counts)
)
fd <- new("AnnotatedDataFrame",
data=gene_ann)
# ---------然后構(gòu)建樣本的注釋信息(sample_data)---------
sample_ann=meta
pd <- new("AnnotatedDataFrame",
data=sample_ann)
# ---------開始構(gòu)建對象---------
sc_cds <- newCellDataSet(
as.matrix(counts),
phenoData = pd,
featureData =fd,
expressionFamily = negbinomial.size(),
lowerDetectionLimit=1)
sc_cds
# CellDataSet (storageMode: environment)
# assayData: 24582 features, 768 samples
# element names: exprs
# protocolData: none
# phenoData
# sampleNames: SS2_15_0048_A3 SS2_15_0048_A6 ... SS2_15_0049_P24 (768 total)
# varLabels: g plate ... Size_Factor (5 total)
# varMetadata: labelDescription
# featureData
# featureNames: 0610005C13Rik 0610007P14Rik ... ERCC-00171 (24582 total)
# fvarLabels: gene_short_name
# fvarMetadata: labelDescription
# experimentData: use 'experimentData(object)'
# Annotation:
接下來質(zhì)控過濾
=> detectGenes()
cds=sc_cds
cds
## 起初是 24582 features, 768 samples
#---------首先是對基因的過濾-------------
cds <- detectGenes(cds, min_expr = 0.1)
print(head(cds@featureData@data))
expressed_genes <- row.names(subset(cds@featureData@data,
num_cells_expressed >= 5))
length(expressed_genes)
# 14442
# 這里需要去掉ERCC基因
# 去掉ERCC基因
is.ercc <- grepl("ERCC",expressed_genes)
length(expressed_genes[!is.ercc])
# 14362(看到去掉了80個ERCC)
cds <- cds[expressed_genes[!is.ercc],]
cds
# 過濾基因后是 14362 features, 768 samples
#---------然后是對細胞的過濾-------------
# 如果不支持使用pData()函數(shù)，可以使用cds@phenoData@data來獲得各種細胞注釋信息
cell_anno <- cds@phenoData@data
> head(cell_anno)
g plate n_g all Size_Factor num_genes_expressed
SS2_15_0048_A3 1 0048 2624 all 0.6693919 3069
SS2_15_0048_A6 1 0048 2664 all 0.6820532 3040
SS2_15_0048_A5 2 0048 3319 all 1.0759418 3743
SS2_15_0048_A4 3 0048 4447 all 0.7294812 5014
SS2_15_0048_A1 2 0048 4725 all 1.5658507 5128
SS2_15_0048_A2 3 0048 5263 all 1.6187569 5693
# 這里簡單過濾細胞
valid_cells <- row.names(cell_anno[cell_anno$num_genes_expressed>2000,] )
cds <- cds[,valid_cells]
cds
# 最后剩下：14362 features, 693 samples
然后進行歸一化
=> estimateSizeFactors()
library(dplyr)
colnames(phenoData(cds)@data)
## 必要的歸一化
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
然后降維聚類
step1：dispersionTable()
disp_table <- dispersionTable(cds)
unsup_clustering_genes <- subset(disp_table,
mean_expression >= 0.1)
cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id)
cds
plot_ordering_genes(cds)
# 圖中黑色的點就是被標(biāo)記出來一會要進行聚類的基因
step2：plot_pc_variance_explained()
plot_pc_variance_explained(cds, return_all = F) # norm_method='log'
step3:
關(guān)于聚類：monocle2一共有三種方法：densityPeak", "louvain", "DDRTree"
默認(rèn)使用density peak的方法，但是對于大型數(shù)據(jù)（例如有5萬細胞）推薦louvain方法
# ---------進行降維---------
cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
reduction_method = 'tSNE', verbose = T)
# ---------進行聚類---------
# (這里先設(shè)置num_clusters，不一定最后真就分成5群；我們這里設(shè)置5，最后只能得到4群；如果設(shè)成4，結(jié)果就得到3群)
cds <- clusterCells(cds, num_clusters = 5)
# Distance cutoff calculated to 1.818667
# color使用的這些數(shù)據(jù)就在：cds$Cluster
plot_cell_clusters(cds, 1, 2, color = "Cluster")
因為之前還做過層次聚類，所以還可以對比一下：
plot_cell_clusters(cds, 1, 2, color = "g")
看到monocle使用其他的聚類算法確實不如使用自己的結(jié)果得到的效果好
再次對比不同分群結(jié)果的基因數(shù)量差異：
boxplot(cds@phenoData@data$num_genes_expressed~cds@phenoData@data$Cluster)
boxplot(cds@phenoData@data$num_genes_expressed~cds@phenoData@data$g)
去除一些影響因素
因為這幾群的細胞中基因表達數(shù)量是有差別的，因此我們可以在聚類之前先去掉這部分影響，從而關(guān)注它們真正的生物學(xué)影響
## 去除檢測到基因數(shù)量效應(yīng)
cds <- reduceDimension(cds, max_components = 2, num_dim = 2,
reduction_method = 'tSNE',
residualModelFormulaStr = "~num_genes_expressed",
verbose = T)
cds <- clusterCells(cds, num_clusters = 5)
plot_cell_clusters(cds, 1, 2, color = "Cluster")
差異分析
=> differentialGeneTest()
start=Sys.time()
diff_test_res <- differentialGeneTest(cds,
fullModelFormulaStr = "~Cluster")
end=Sys.time()
end-start
# Time difference of 4.724045 mins
挑差異基因
# 選擇FDR-adjusted p-value(也就是q值) < 10%的基因作為差異基因
sig_genes <- subset(diff_test_res, qval < 0.1)
dim(sig_genes)
> head(sig_genes[,c("gene_short_name", "pval", "qval")] )
gene_short_name pval qval
0610007P14Rik 0610007P14Rik 3.377244e-03 1.277429e-02
0610010F05Rik 0610010F05Rik 1.243943e-02 3.924761e-02
0610011F06Rik 0610011F06Rik 2.338530e-03 9.285587e-03
1110004E09Rik 1110004E09Rik 2.487600e-02 7.003903e-02
1110020A21Rik 1110020A21Rik 2.318327e-02 6.618129e-02
1110059E24Rik 1110059E24Rik 5.193533e-09 4.856089e-08
作圖
cg=as.character(head(sig_genes$gene_short_name))
# 還能上色
plot_genes_jitter(cds[cg,],
grouping = "Cluster",
color_by = "Cluster",
nrow= 3,
ncol = NULL )
推斷發(fā)育軌跡
step1: 選合適基因 => setOrderingFilter()
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)
step2: 降維 => reduceDimension()
# 默認(rèn)使用DDRTree的方法
cds <- reduceDimension(cds, max_components = 2,
method = 'DDRTree')
step3: 細胞排序 => orderCells()
cds <- orderCells(cds)
最后可視化
plot_cell_trajectory(cds, color_by = "Cluster")
