生物学背景知识之细胞周期推断

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)

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):

# 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"))

# 将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 ...

> 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)

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()