生物学背景知识之细胞周期推断
- 2020 年 3 月 31 日
- 筆記
第二单元第九讲：生物学背景知识之细胞周期推断
上一次说到通过PAM50基因进行乳腺癌分型，利用的就是自己的表达矩阵和PAM50基因比较，看表达量变化进行分类。细胞周期分类和PAM50类似，也是利用基因来推断G、S、M期(https://en.wikipedia.org/wiki/Cell_cycle)
Scran包使用
依然第一步是加载矩阵
rm(list = ls()) options(stringsAsFactors = F) load(file = '../input.Rdata') a[1:4,1:4] head(df) # 放入分群、样本批次信息 group_list=df$g plate=df$plate table(plate)
然后创建sce对象
library(scran) sce <- SingleCellExperiment(list(counts=dat)) > sce class: SingleCellExperiment dim: 12198 768 metadata(0): assays(1): counts rownames(12198): 0610007P14Rik 0610009B22Rik ... ERCC-00170 ERCC-00171 rowData names(0): colnames(768): SS2_15_0048_A3 SS2_15_0048_A6 ... SS2_15_0049_P22 SS2_15_0049_P24 colData names(0): reducedDimNames(0): spikeNames(0):
主要使用cyclone函数
cyclone函数主要需要三个元素：一个是sce单细胞对象，一个是pairs参数，还有就是gene.names参数。第一个已准备好，第二个参数的意思可以看帮助文档
# scran包安装好后，会在exdata文件夹中找到附件文件 library(org.Mm.eg.db) # syste,.file会列出文件所在的路径，下图就是exdata文件夹下的文件，看到除了小鼠还有人的相关的RDS数据。这个RDS其实和平常看到的Rdata差不多，只不过Rdata是针对多个对象，Rds是针对一个对象进行存储和读取 mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))
然后是第三个参数：gene.names，cyclone函数需要使用ensembl基因名
# 将symbol转为ensembl基因 ensembl <- mapIds(org.Mm.eg.db, keys=rownames(sce), keytype="SYMBOL", column="ENSEMBL") > head(ensembl) 0610007P14Rik 0610009B22Rik 0610009L18Rik NA "ENSMUSG00000007777" "ENSMUSG00000043644" 0610009O20Rik 0610010F05Rik 0610010K14Rik NA "ENSMUSG00000042208" "ENSMUSG00000020831"
三者齐全，可以进行细胞周期计算：
system.time(assigned <- cyclone(sce, pairs=mm.pairs, gene.names=ensembl)) # 这一过程会比较慢，用system.time计算一下时间看看，大约一分半 # user system elapsed # 96.229 0.767 104.666 save(assigned,file = 'cell_cycle_assigned.Rdata') > str(assigned) # 包含了phases、scores、normalized.scores三个元素 List of 3 $ phases : chr [1:768] "G1" "G1" "G1" "G1" ... $ scores :'data.frame':768 obs. of 3 variables: ..$ G1 : num [1:768] 1 0.997 0.997 1 1 1 1 0.937 1 1 ... ..$ S : num [1:768] 0.119 0.002 0.039 0.011 0.395 0.009 0.011 0.008 0.04 0.013 ... ..$ G2M: num [1:768] 0.004 0.01 0.02 0.002 0 0 0.02 0.126 0 0.023 ... $ normalized.scores:'data.frame':768 obs. of 3 variables: ..$ G1 : num [1:768] 0.89 0.988 0.944 0.987 0.717 ... ..$ S : num [1:768] 0.10597 0.00198 0.03693 0.01086 0.28315 ... ..$ G2M: num [1:768] 0.00356 0.00991 0.01894 0.00197 0 ...
下面就根据assigned进行操作
> head(assigned$scores) G1 S G2M 1 1.000 0.119 0.004 2 0.997 0.002 0.010 3 0.997 0.039 0.020 4 1.000 0.011 0.002 5 1.000 0.395 0.000 6 1.000 0.009 0.000 > table(assigned$phases) G1 G2M S 723 34 11 # 作图(利用score和phases这两个元素) draw=cbind(assigned$score,assigned$phases) attach(draw) #attach的目的就是现在加载，之后直接引用即可 library(scatterplot3d) scatterplot3d(G1, S, G2M, angle=20, color = rainbow(3)[as.numeric(as.factor(assigned$phases))], grid=TRUE, box=FALSE) detach(draw)
还能做个热图(就是在anno_col上不断加内容即可)
library(pheatmap) # 取差异前100基因 cg=names(tail(sort(apply(dat,1,sd)),100)) # 矩阵归一化 n=t(scale(t(dat[cg,]))) # 原来的样本注释信息 df中包含了 g、plate 、n_g、all信息，现在新增phases信息 df$cellcycle=assigned$phases ac=df rownames(ac)=colnames(n) pheatmap(n,show_colnames =F,show_rownames = F, annotation_col=ac) dev.off()
