前言
前面得到的6个发育时期和4个分群，而且还可视化了一些marker基因，那么现在就要对这4群细胞进行差异分析
将对应文章这张图：
1 使用monocle V2进行差异分析
1.1 准备表达矩阵和分群信息
代码语言：javascript复制rm(list = ls())
options(warn=-1)
options(stringsAsFactors = F)
source("../analysis_functions.R")
# 差异分析一般都要求使用count矩阵
load('../female_count.Rdata')
# 6个发育时期获取
head(colnames(female_count))
female_stages <- sapply(strsplit(colnames(female_count), "_"), `[`, 1)
names(female_stages) <- colnames(female_count)
table(female_stages)
# 4个cluster获取
cluster <- read.csv('female_clustering.csv')
female_clustering=cluster[,2];names(female_clustering)=cluster[,1]
table(female_clustering)
1.2 利用monocle V2构建对象
需要三样东西：表达矩阵、细胞信息、基因信息
代码语言：javascript复制## 表达矩阵
dim(female_count)
## 细胞分群信息（包括6个stage和4个cluster）
table(female_stages)
table(female_clustering)
## 基因信息
gene_annotation <- as.data.frame(rownames(count_matrix))
开始构建
代码语言：javascript复制library(monocle)
# 直接使用作者包装的函数，代码更简洁
DE_female <- prepare_for_DE (
female_count,
female_clustering,
female_stages
)
可以看到，包装好的函数只需要提供表达矩阵和细胞信息，其实它背后做的事情比较多，就像下面?的代码
代码语言：javascript复制 # 1.表达矩阵
expr_matrix <- as.matrix(count_matrix)
# 2.细胞信息
sample_ann <- data.frame(cells=names(count_matrix),
stages=stage,
cellType=cluster)
rownames(sample_ann)<- names(count_matrix)
# 3.基因信息
gene_ann <- as.data.frame(rownames(count_matrix))
rownames(gene_ann)<- rownames(count_matrix)
colnames(gene_ann)<- "genes"
# 然后转换为AnnotatedDataFrame对象
pd <- new("AnnotatedDataFrame",
data=sample_ann)
fd <- new("AnnotatedDataFrame",
data=gene_ann)
# 最后构建CDS对象
sc_cds <- newCellDataSet(
expr_matrix,
phenoData = pd,
featureData =fd,
expressionFamily = negbinomial.size(),
lowerDetectionLimit=0.5)
sc_cds <- detectGenes(sc_cds, min_expr = 5)
sc_cds <- sc_cds[fData(sc_cds)$num_cells_expressed > 10, ]
sc_cds <- estimateSizeFactors(sc_cds)
sc_cds <- estimateDispersions(sc_cds)
1.3 进行差异分析
其实也是包装了differentialGeneTest函数
代码语言：javascript复制start_time <- Sys.time()
female_DE_genes <- findDEgenes(
DE_female,
qvalue=0.05
)
end_time <- Sys.time()
end_time - start_time
# 1] "4435 significantly DE genes (FDR<0.05)."
# [1] "3268 significantly DE genes (FDR<0.01)."
# Time difference of 1.796054 mins
我们最后保留的也就是4435 significantly DE genes (FDR<0.05).
实际上，这个包装的函数就是下面这样，我们只需要提供创建好的monocle对象和q值即可：
代码语言：javascript复制function(HSMM=HSMM, qvalue=qvalue){
diff_test_res <- differentialGeneTest(
HSMM,
fullModelFormulaStr="~cellType",
cores = 3
)
sig_genes_0.05 <- subset(diff_test_res, qval < 0.05)
sig_genes_0.01 <- subset(diff_test_res, qval < 0.01)
# 帮我们统计了不同q值得到的差异基因数量
print(paste(nrow(sig_genes_0.05), " significantly DE genes (FDR<0.05).", sep=""))
print(paste(nrow(sig_genes_0.01), " significantly DE genes (FDR<0.01).", sep=""))
diff_test_res <- subset(diff_test_res, qval< qvalue)
return(diff_test_res)
}
最后得到这么一个差异基因数据框：
上面得到的4435个差异基因，我们要知道这些基因在哪个cluster：
作者的做法是：先得到了每个差异基因在不同cluster的平均表达量，然后找平均表达量最大的那个cluster，就认为这个基因属于这个cluster
代码语言：javascript复制# 有了差异基因以后，继续使用RPKM值
load('../female_rpkm.Rdata')
de_clusters <- get_up_reg_clusters(
females,
female_clustering,
female_DE_genes
)
2 使用ROTS进行差异分析
2.1 准备表达矩阵和分群信息
代码语言：javascript复制rm(list = ls())
options(warn=-1)
options(stringsAsFactors = F)
# 6个发育时期获取
head(colnames(female_count))
female_stages <- sapply(strsplit(colnames(female_count), "_"), `[`, 1)
names(female_stages) <- colnames(female_count)
table(female_stages)
# 4个cluster获取
cluster <- read.csv('female_clustering.csv')
female_clustering=cluster[,2];names(female_clustering)=cluster[,1]
table(female_clustering)
2.2 使用ROTS
下面会使用ROTS包（Reproducibility-optimized test statistic），对每个群和其他几个群的共同体进行比较
- ROTS可以使用RPKM值；它的速度可能会很慢
- 差异分析重点就在：表达矩阵和分组信息
配置最重要的分组信息
代码语言：javascript复制library(ROTS)
library(plyr)
# 使用RPKM矩阵
load('../female_rpkm.Rdata')
# 当前分组是这样
> table(female_clustering)
female_clustering
C1 C2 C3 C4
240 90 190 43
# 我们只需要用数字来表示，于是可以用substring(x, start, stop)获取
groups<-female_clustering
groups<-as.numeric(substring(as.character(groups),2,2))
> table(groups)
groups
1 2 3 4
240 90 190 43
# 如果要对第1群进行差异分析，那么就要把它和其他群比较(可将其他亚群定义为234)
groups[groups!=1]<-234
# 同理，对第2群进行差异分析
groups[groups!=2]<-134
然后配置表达矩阵
代码语言：javascript复制RPKM.full=females
ROTS_input<-RPKM.full[rowMeans(RPKM.full)>=1,]
ROTS_input<-as.matrix(log2(ROTS_input 1))
第1群
代码语言：javascript复制start_time <- Sys.time()
results_pop1 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234)
summary_pop1<-data.frame(summary(results_pop1, fdr=1))
end_time <- Sys.time()
(end_time - start_time)
# Time difference of 18.40462 mins
然后第2群、第3群、第4群
代码语言：javascript复制groups[groups!=2]<-134
start_time <- Sys.time()
results_pop2 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234)
summary_pop2<-data.frame(summary(results_pop2, fdr=1))
end_time <- Sys.time()
(end_time - start_time)
# 第3群
groups[groups!=3]<-124
start_time <- Sys.time()
results_pop3 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234)
summary_pop3<-data.frame(summary(results_pop3, fdr=1))
end_time <- Sys.time()
(end_time - start_time)
# 第4群
groups[groups!=4]<-123
start_time <- Sys.time()
results_pop4 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234)
summary_pop3<-data.frame(summary(results_pop4, fdr=1))
end_time <- Sys.time()
(end_time - start_time)
# 都得到以后，共同保存
save(summary_pop1,summary_pop2,summary_pop3,summary_pop4,
file = 'ROTS_summary.Rdata')
