课程学习生信技能树单细胞转录组（基础）
现在福利来了，单细胞课程免费开放了：#完整版视频上线！全！面！免！费！#全网第一的单细胞转录组实战演练_ https://www.bilibili.com/video/av38741055?p=11
不同类型的细胞哪怕是在同一个时间点的细胞周期状态，它们的细胞周期相关基因表达也是不同的。更重要的是，还有很多非直接细胞周期相关基因也需要考虑。
细胞周期是利用基因来推断G、S、M期(https://en.wikipedia.org/wiki/Cell_cycle)
cell cycle
先复现一下视频里面的代码：
Scran包使用推断细胞周期
载入数据和
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = '../input.Rdata')
a[1:4,1:4]
head(df)
## 载入第0步准备好的表达矩阵，及细胞的一些属性（hclust分群，plate批次，检测到的基因数量）
# 注意 变量a是原始的counts矩阵，变量 dat是logCPM后的表达量矩阵。
group_list=df$g
plate=df$plate
table(plate)
a[1:4,1:4]
创建sce对象
library(scran)
# https://mp.weixin.qq.com/s/nFSa5hXuKHrGu_othopbWQ
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参数。
第二个pairs参数的意思可以看帮助文档
# 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()
探索：scran包中的cyclone函数细胞周期原理
主要利用了scran包中的cyclone函数
cyclone函数主要需要三个元素：一个是sce单细胞对象，一个是pairs参数，还有就是gene.names参数。第一个已准备好，第二个参数的意思可以看帮助文档，第三个参数要求是Ensembl ID
# 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"))
# 举个例子
head(mm.pairs$G1)
first second
1 ENSMUSG00000000001 ENSMUSG00000001785
2 ENSMUSG00000000001 ENSMUSG00000005470
3 ENSMUSG00000000001 ENSMUSG00000012443
4 ENSMUSG00000000001 ENSMUSG00000015120
5 ENSMUSG00000000001 ENSMUSG00000022033
6 ENSMUSG00000000001 ENSMUSG00000023015
注意：这里小鼠的训练数据集是利用胚胎干细胞数据得到的，但对于其他细胞类型也是准确的 (Scialdone et al. 2015)，可能是由于细胞周期相关转录的保守性 (Bertoli, Skotheim, and Bruin 2013; Conboy et al. 2007)。
具体的使用就很简单：
system.time(assignments <- cyclone(sce, mm.pairs,
gene.names=rowData(sce)$ENSEMBL))
## user system elapsed
## 21.740 0.376 26.856
save(sce,assignments,file = '416B_cell_cycle.Rdata')
这个结果并不是说某个细胞一定就处于哪个细胞周期，它也是根据大样本量背景进行概率估计，然后计算一个分值，分值越高，说明某个细胞更有可能属于哪个细胞周期，给我们一定的参考。
例如，看一下G1和G2/M的分值情况，其中每个点代表一个细胞：
plot(assignments$score$G1, assignments$score$G2M,
xlab="G1 score", ylab="G2/M score", pch=16)
计算完细胞分值，就该划分细胞类型了
具体规则是：如果一个细胞在G1中得分大于0.5，并且它高于G2/M的得分，那么这个细胞就被划分到G1期；如果细胞在G2/M中得分大于0.5，并且高于G1的得分，那么它就划为G2/M期；如果细胞的G1、G2/M得分都不大于0.5，那么它就划为S期。于是上面👆的结果就可以被划分成：
这里只是说明一下原理， 其实函数已经为我们划分好，存到了结果中：
sce$phases <- assignments$phases
table(sce$phases)
##
## G1 G2M S
## 98 62 23
作者的建议
- 为了去掉细胞周期带来的影响，我们一般只要某一群特定周期的细胞（一般是G1期）进行下游分析。当然如果其他时期的细胞数量不会造成明显的影响，我们也可以带着它们，到下游分析时直接将
assignments$phases作为一个批次效应因素考虑即可，这样既避免了其他周期细胞的干扰，又能避免丢失部分信息。 - 训练数据集虽然说对多种细胞类型都支持，如果本身的数据与训练数据相差太大（比如使用了不同的方法得到的数据），那么我们可以根据自己的数据去DIY一个分类器。使用
sandbag函数即可，同样如果研究其他的物种没有给定的分类器，我们就可以这样操作 - 在使用
cyclone之前不要过滤低丰度转录本。即使一个基因在任何一个细胞都不表达，但在周期推断环节这个基因名还是有作用的。因为这一步做的是一个基因对比较，所以它依然可以提供信息。
参考：
1.单细胞转录组学习笔记-11-生物学背景知识之细胞周期推断 - 简书 https://www.jianshu.com/p/46d597d21a16
2.关于细胞周期推断的知识更新 - 简书 https://www.jianshu.com/p/455047d7557c
3.在单细胞转录组表达矩阵里面去除细胞周期影响 - 简书 https://www.jianshu.com/p/aa867c3c12de
4.ccRemover说明书
5.清除单细胞数据细胞周期效应造成的影响
