劉小澤寫于19.7.9-第二單元第九講：生物學(xué)背景知識之細胞周期推斷
筆記目的：根據(jù)生信技能樹的單細胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
前言
上一次說到通過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)
然后創(chuàng)建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函數(shù)
cyclone函數(shù)主要需要三個元素：一個是sce單細胞對象，一個是pairs參數(shù)钢属，還有就是gene.names參數(shù)徘熔。第一個已準備好，第二個參數(shù)的意思可以看幫助文檔
# scran包安裝好后淆党，會在exdata文件夾中找到附件文件
library(org.Mm.eg.db)
# syste,.file會列出文件所在的路徑酷师，下圖就是exdata文件夾下的文件，看到除了小鼠還有人的相關(guān)的RDS數(shù)據(jù)染乌。這個RDS其實和平成娇祝看到的Rdata差不多，只不過Rdata是針對多個對象荷憋，Rds是針對一個對象進行存儲和讀取
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds",
package="scran"))
然后是第三個參數(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"
三者齊全，可以進行細胞周期計算：
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 ...
下面就根據(jù)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的目的就是現(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)
還能做個熱圖(就是在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()
