課程學(xué)習(xí)生信技能樹單細(xì)胞轉(zhuǎn)錄組（基礎(chǔ)）
現(xiàn)在福利來了，單細(xì)胞課程免費(fèi)開放了：#完整版視頻上線！全卵洗！面！免弥咪！費(fèi)过蹂！#全網(wǎng)第一的單細(xì)胞轉(zhuǎn)錄組實(shí)戰(zhàn)演練_ https://www.bilibili.com/video/av38741055?p=11
不同類型的細(xì)胞哪怕是在同一個(gè)時(shí)間點(diǎn)的細(xì)胞周期狀態(tài)，它們的細(xì)胞周期相關(guān)基因表達(dá)也是不同的聚至。更重要的是酷勺，還有很多非直接細(xì)胞周期相關(guān)基因也需要考慮。
細(xì)胞周期是利用基因來推斷G扳躬、S脆诉、M期(https://en.wikipedia.org/wiki/Cell_cycle)
cell cycle
先復(fù)現(xiàn)一下視頻里面的代碼：
Scran包使用推斷細(xì)胞周期
載入數(shù)據(jù)和
rm(list = ls()) ## 魔幻操作甚亭，一鍵清空~
options(stringsAsFactors = F)
load(file = '../input.Rdata')
a[1:4,1:4]
head(df)
## 載入第0步準(zhǔn)備好的表達(dá)矩陣，及細(xì)胞的一些屬性（hclust分群库说，plate批次狂鞋，檢測(cè)到的基因數(shù)量）
# 注意 變量a是原始的counts矩陣，變量 dat是logCPM后的表達(dá)量矩陣潜的。
group_list=df$g
plate=df$plate
table(plate)
a[1:4,1:4]
創(chuàng)建sce對(duì)象
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函數(shù)
cyclone函數(shù)主要需要三個(gè)元素：
一個(gè)是sce單細(xì)胞對(duì)象骚揍，
一個(gè)是pairs參數(shù)，
還有就是gene.names參數(shù)啰挪。
第二個(gè)pairs參數(shù)的意思可以看幫助文檔
# scran包安裝好后信不，會(huì)在exdata文件夾中找到附件文件
library(org.Mm.eg.db)
# syste,.file會(huì)列出文件所在的路徑，下圖就是exdata文件夾下的文件亡呵，看到除了小鼠還有人的相關(guān)的RDS數(shù)據(jù)抽活。這個(gè)RDS其實(shí)和平常看到的Rdata差不多锰什，只不過Rdata是針對(duì)多個(gè)對(duì)象下硕，Rds是針對(duì)一個(gè)對(duì)象進(jìn)行存儲(chǔ)和讀取
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds",
package="scran"))
第三個(gè)參數(shù)：gene.names，cyclone函數(shù)需要使用ensembl基因名
# 將symbol轉(zhuǎn)為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"
三者齊全汁胆，可以進(jìn)行細(xì)胞周期計(jì)算：
system.time(assigned <- cyclone(sce, pairs=mm.pairs, gene.names=ensembl))
# 這一過程會(huì)比較慢梭姓，用system.time計(jì)算一下時(shí)間看看，大約一分半
# user system elapsed
# 96.229 0.767 104.666
save(assigned,file = 'cell_cycle_assigned.Rdata')
str(assigned) # 包含了phases嫩码、scores誉尖、normalized.scores三個(gè)元素
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 ...
下面就根據(jù)assigned進(jìn)行操作
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這兩個(gè)元素)
draw=cbind(assigned$score,assigned$phases)
attach(draw) #attach的目的就是現(xiàn)在加載，之后直接引用即可
library(scatterplot3d)
scatterplot3d(G1, S, G2M, angle=20,
color = rainbow(3)[as.numeric(as.factor(assigned$phases))],
grid=TRUE, box=FALSE)
detach(draw)
還能做個(gè)熱圖(就是在anno_col上不斷加內(nèi)容即可)
library(pheatmap)
# 取差異前100基因
cg=names(tail(sort(apply(dat,1,sd)),100))
# 矩陣歸一化
n=t(scale(t(dat[cg,])))
# 原來的樣本注釋信息 df中包含了 g铸题、plate 铡恕、n_g、all信息丢间，現(xiàn)在新增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函數(shù)細(xì)胞周期原理
主要利用了scran包中的cyclone函數(shù)
cyclone函數(shù)主要需要三個(gè)元素：一個(gè)是sce單細(xì)胞對(duì)象探熔，一個(gè)是pairs參數(shù)，還有就是gene.names參數(shù)千劈。第一個(gè)已準(zhǔn)備好祭刚，第二個(gè)參數(shù)的意思可以看幫助文檔，第三個(gè)參數(shù)要求是Ensembl ID
# scran包安裝好后墙牌，會(huì)在exdata文件夾中找到附件文件
library(org.Mm.eg.db)
# syste,.file會(huì)列出文件所在的路徑涡驮，下圖就是exdata文件夾下的文件，看到除了小鼠還有人的相關(guān)的RDS數(shù)據(jù)喜滨。這個(gè)RDS其實(shí)和平匙酵保看到的Rdata差不多，只不過Rdata是針對(duì)多個(gè)對(duì)象虽风，Rds是針對(duì)一個(gè)對(duì)象進(jìn)行存儲(chǔ)和讀取
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds",
package="scran"))
# 舉個(gè)例子
head(mm.pairs$G1)
first second
1 ENSMUSG00000000001 ENSMUSG00000001785
2 ENSMUSG00000000001 ENSMUSG00000005470
3 ENSMUSG00000000001 ENSMUSG00000012443
4 ENSMUSG00000000001 ENSMUSG00000015120
5 ENSMUSG00000000001 ENSMUSG00000022033
6 ENSMUSG00000000001 ENSMUSG00000023015
注意：這里小鼠的訓(xùn)練數(shù)據(jù)集是利用胚胎干細(xì)胞數(shù)據(jù)得到的棒口，但對(duì)于其他細(xì)胞類型也是準(zhǔn)確的 (Scialdone et al. 2015)寄月，可能是由于細(xì)胞周期相關(guān)轉(zhuǎn)錄的保守性 (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')
這個(gè)結(jié)果并不是說某個(gè)細(xì)胞一定就處于哪個(gè)細(xì)胞周期无牵，它也是根據(jù)大樣本量背景進(jìn)行概率估計(jì)漾肮，然后計(jì)算一個(gè)分值，分值越高茎毁，說明某個(gè)細(xì)胞更有可能屬于哪個(gè)細(xì)胞周期克懊，給我們一定的參考。
例如七蜘，看一下G1和G2/M的分值情況谭溉，其中每個(gè)點(diǎn)代表一個(gè)細(xì)胞：
plot(assignments$score$G1, assignments$score$G2M,
xlab="G1 score", ylab="G2/M score", pch=16)
計(jì)算完細(xì)胞分值，就該劃分細(xì)胞類型了
具體規(guī)則是：如果一個(gè)細(xì)胞在G1中得分大于0.5橡卤，并且它高于G2/M的得分扮念，那么這個(gè)細(xì)胞就被劃分到G1期；如果細(xì)胞在G2/M中得分大于0.5碧库，并且高于G1的得分柜与，那么它就劃為G2/M期；如果細(xì)胞的G1嵌灰、G2/M得分都不大于0.5旅挤，那么它就劃為S期。于是上面??的結(jié)果就可以被劃分成：
這里只是說明一下原理伞鲫， 其實(shí)函數(shù)已經(jīng)為我們劃分好，存到了結(jié)果中：
sce$phases <- assignments$phases
table(sce$phases)
##
## G1 G2M S
## 98 62 23
作者的建議
- 為了去掉細(xì)胞周期帶來的影響签舞，我們一般只要某一群特定周期的細(xì)胞（一般是G1期）進(jìn)行下游分析秕脓。當(dāng)然如果其他時(shí)期的細(xì)胞數(shù)量不會(huì)造成明顯的影響，我們也可以帶著它們儒搭，到下游分析時(shí)直接將
assignments$phases作為一個(gè)批次效應(yīng)因素考慮即可吠架，這樣既避免了其他周期細(xì)胞的干擾，又能避免丟失部分信息搂鲫。 - 訓(xùn)練數(shù)據(jù)集雖然說對(duì)多種細(xì)胞類型都支持傍药，如果本身的數(shù)據(jù)與訓(xùn)練數(shù)據(jù)相差太大（比如使用了不同的方法得到的數(shù)據(jù)），那么我們可以根據(jù)自己的數(shù)據(jù)去DIY一個(gè)分類器魂仍。使用
sandbag函數(shù)即可拐辽，同樣如果研究其他的物種沒有給定的分類器，我們就可以這樣操作 - 在使用
cyclone之前不要過濾低豐度轉(zhuǎn)錄本擦酌。即使一個(gè)基因在任何一個(gè)細(xì)胞都不表達(dá)俱诸，但在周期推斷環(huán)節(jié)這個(gè)基因名還是有作用的。因?yàn)檫@一步做的是一個(gè)基因?qū)Ρ容^赊舶，所以它依然可以提供信息睁搭。
參考：
1.單細(xì)胞轉(zhuǎn)錄組學(xué)習(xí)筆記-11-生物學(xué)背景知識(shí)之細(xì)胞周期推斷 - 簡書 http://www.reibang.com/p/46d597d21a16
2.關(guān)于細(xì)胞周期推斷的知識(shí)更新 - 簡書 http://www.reibang.com/p/455047d7557c
3.在單細(xì)胞轉(zhuǎn)錄組表達(dá)矩陣?yán)锩嫒コ?xì)胞周期影響 - 簡書 http://www.reibang.com/p/aa867c3c12de
4.ccRemover說明書
5.清除單細(xì)胞數(shù)據(jù)細(xì)胞周期效應(yīng)造成的影響
