差异分析及功能注释（上） / 开普饭

作者 | 单细胞天地小编刘小泽

课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
这次会介绍如何进行差异分析及功能注释，以及比较不同分组。对应视频第三单元9-11讲

前言

前面得到的6个发育时期和4个分群，而且还可视化了一些marker基因，那么现在就要对这4群细胞进行差异分析

将对应文章这张图：

1 使用monocle V2进行差异分析

1.1 准备表达矩阵和分群信息

1rm(list = ls()) 2options(warn=-1) 3options(stringsAsFactors = F) 4source("../analysis_functions.R") 5 6 # 差异分析一般都要求使用count矩阵 7load('../female_count.Rdata') 8 9# 6个发育时期获取 10head(colnames(female_count)) 11female_stages <- sapply(strsplit(colnames(female_count), "_"), `[`, 1) 12names(female_stages) <- colnames(female_count) 13table(female_stages) 14 15# 4个cluster获取 16cluster <- read.csv('female_clustering.csv') 17female_clustering=cluster[,2];names(female_clustering)=cluster[,1] 18table(female_clustering)

1.2 利用monocle V2构建对象

需要三样东西：表达矩阵、细胞信息、基因信息

1## 表达矩阵 2dim(female_count) 3## 细胞分群信息（包括6个stage和4个cluster） 4table(female_stages) 5table(female_clustering) 6## 基因信息 7gene_annotation <- as.data.frame(rownames(count_matrix))

开始构建

1library(monocle) 2# 直接使用作者包装的函数，代码更简洁 3DE_female <- prepare_for_DE ( 4 female_count, 5 female_clustering, 6 female_stages 7)

可以看到，包装好的函数只需要提供表达矩阵和细胞信息，其实它背后做的事情比较多，就像下面👇的代码

1 # 1.表达矩阵 2 expr_matrix <- as.matrix(count_matrix) 3 # 2.细胞信息 4 sample_ann <- data.frame(cells=names(count_matrix), 5 stages=stage, 6 cellType=cluster) 7 rownames(sample_ann)<- names(count_matrix) 8 # 3.基因信息 9 gene_ann <- as.data.frame(rownames(count_matrix)) 10 rownames(gene_ann)<- rownames(count_matrix) 11 colnames(gene_ann)<- "genes" 12 13 # 然后转换为AnnotatedDataFrame对象 14 pd <- new("AnnotatedDataFrame", 15 data=sample_ann) 16 fd <- new("AnnotatedDataFrame", 17 data=gene_ann) 18 19 # 最后构建CDS对象 20 sc_cds <- newCellDataSet( 21 expr_matrix, 22 phenoData = pd, 23 featureData =fd, 24 expressionFamily = negbinomial.size(), 25 lowerDetectionLimit=0.5) 26 27 sc_cds <- detectGenes(sc_cds, min_expr = 5) 28 sc_cds <- sc_cds[fData(sc_cds)$num_cells_expressed > 10, ] 29 sc_cds <- estimateSizeFactors(sc_cds) 30 sc_cds <- estimateDispersions(sc_cds)

1.3 进行差异分析

其实也是包装了differentialGeneTest函数

1start_time <- Sys.time() 2female_DE_genes <- findDEgenes( 3 DE_female, 4 qvalue=0.05 5) 6end_time <- Sys.time() 7end_time - start_time 8 9# 1] "4435 significantly DE genes (FDR<0.05)." 10# [1] "3268 significantly DE genes (FDR<0.01)." 11# Time difference of 1.796054 mins

我们最后保留的也就是4435 significantly DE genes (FDR<0.05).

实际上，这个包装的函数就是下面这样，我们只需要提供创建好的monocle对象和q值即可：

1function(HSMM=HSMM, qvalue=qvalue){ 2 diff_test_res <- differentialGeneTest( 3 HSMM, 4 fullModelFormulaStr="~cellType", 5 cores = 3 6 ) 7 8 sig_genes_0.05 <- subset(diff_test_res, qval < 0.05) 9 sig_genes_0.01 <- subset(diff_test_res, qval < 0.01) 10 # 帮我们统计了不同q值得到的差异基因数量 11 print(paste(nrow(sig_genes_0.05), " significantly DE genes (FDR<0.05).", sep="")) 12 print(paste(nrow(sig_genes_0.01), " significantly DE genes (FDR<0.01).", sep="")) 13 14 diff_test_res <- subset(diff_test_res, qval< qvalue) 15 16 return(diff_test_res) 17}

最后得到这么一个差异基因数据框：

上面得到的4435个差异基因，我们要知道这些基因在哪个cluster：

作者的做法是：先得到了每个差异基因在不同cluster的平均表达量，然后找平均表达量最大的那个cluster，就认为这个基因属于这个cluster

1# 有了差异基因以后，继续使用RPKM值 2load('../female_rpkm.Rdata') 3de_clusters <- get_up_reg_clusters( 4 females, 5 female_clustering, 6 female_DE_genes 7)

2 使用ROTS进行差异分析

2.1 准备表达矩阵和分群信息

1rm(list = ls()) 2options(warn=-1) 3options(stringsAsFactors = F) 4 5# 6个发育时期获取 6head(colnames(female_count)) 7female_stages <- sapply(strsplit(colnames(female_count), "_"), `[`, 1) 8names(female_stages) <- colnames(female_count) 9table(female_stages) 10# 4个cluster获取 11cluster <- read.csv('female_clustering.csv') 12female_clustering=cluster[,2];names(female_clustering)=cluster[,1] 13table(female_clustering)

2.2 使用ROTS

下面会使用ROTS包（Reproducibility-optimized test statistic），对每个群和其他几个群的共同体进行比较

ROTS可以使用RPKM值；它的速度可能会很慢
差异分析重点就在：表达矩阵和分组信息

配置最重要的分组信息

1library(ROTS) 2library(plyr) 3# 使用RPKM矩阵 4load('../female_rpkm.Rdata') 5 6# 当前分组是这样 7> table(female_clustering) 8female_clustering 9 C1 C2 C3 C4 10240 90 190 43 11 12# 我们只需要用数字来表示，于是可以用substring(x, start, stop)获取 13groups<-female_clustering 14groups<-as.numeric(substring(as.character(groups),2,2)) 15> table(groups) 16groups 17 1 2 3 4 18240 90 190 43 19 20# 如果要对第1群进行差异分析，那么就要把它和其他群比较(可将其他亚群定义为234) 21groups[groups!=1]<-234 22# 同理，对第2群进行差异分析 23groups[groups!=2]<-134

然后配置表达矩阵

1RPKM.full=females 2ROTS_input<-RPKM.full[rowMeans(RPKM.full)>=1,] 3ROTS_input<-as.matrix(log2(ROTS_input+1))

第1群

1start_time <- Sys.time() 2results_pop1 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234) 3summary_pop1<-data.frame(summary(results_pop1, fdr=1)) 4end_time <- Sys.time() 5(end_time - start_time) 6# Time difference of 18.40462 mins

然后第2群、第3群、第4群

1groups[groups!=2]<-134 2start_time <- Sys.time() 3results_pop2 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234) 4summary_pop2<-data.frame(summary(results_pop2, fdr=1)) 5end_time <- Sys.time() 6(end_time - start_time) 7 8# 第3群 9groups[groups!=3]<-124 10start_time <- Sys.time() 11results_pop3 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234) 12summary_pop3<-data.frame(summary(results_pop3, fdr=1)) 13end_time <- Sys.time() 14(end_time - start_time) 15 16# 第4群 17groups[groups!=4]<-123 18start_time <- Sys.time() 19results_pop4 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234) 20summary_pop3<-data.frame(summary(results_pop4, fdr=1)) 21end_time <- Sys.time() 22(end_time - start_time) 23 24# 都得到以后，共同保存 25save(summary_pop1,summary_pop2,summary_pop3,summary_pop4, 26 file = 'ROTS_summary.Rdata')

“

差异分析及功能注释（上）

前言