劉小澤寫(xiě)于2020.7.21
為何取名叫“交響樂(lè)”为牍？因?yàn)閱渭?xì)胞分析就像一個(gè)大樂(lè)團(tuán)洞难，需要各個(gè)流程的協(xié)同配合
單細(xì)胞交響樂(lè)1-常用的數(shù)據(jù)結(jié)構(gòu)SingleCellExperiment
單細(xì)胞交響樂(lè)2-scRNAseq從實(shí)驗(yàn)到下游簡(jiǎn)介
單細(xì)胞交響樂(lè)3-細(xì)胞質(zhì)控
單細(xì)胞交響樂(lè)4-歸一化
單細(xì)胞交響樂(lè)5-挑選高變化基因
單細(xì)胞交響樂(lè)6-降維
單細(xì)胞交響樂(lè)7-聚類(lèi)分群
單細(xì)胞交響樂(lè)8-marker基因檢測(cè)
單細(xì)胞交響樂(lè)9-細(xì)胞類(lèi)型注釋
單細(xì)胞交響樂(lè)9-細(xì)胞類(lèi)型注釋
單細(xì)胞交響樂(lè)10-數(shù)據(jù)集整合后的批次矯正
單細(xì)胞交響樂(lè)11-多樣本間差異分析
單細(xì)胞交響樂(lè)12-檢測(cè)Doublet
單細(xì)胞交響樂(lè)13-細(xì)胞周期推斷
單細(xì)胞交響樂(lè)14-細(xì)胞軌跡推斷
單細(xì)胞交響樂(lè)15-scRNA與蛋白豐度信息結(jié)合
單細(xì)胞交響樂(lè)16-處理大型數(shù)據(jù)
單細(xì)胞交響樂(lè)17-不同單細(xì)胞R包的數(shù)據(jù)格式相互轉(zhuǎn)換
單細(xì)胞交響樂(lè)18-實(shí)戰(zhàn)一 Smart-seq2
單細(xì)胞交響樂(lè)19-實(shí)戰(zhàn)二 STRT-Seq
單細(xì)胞交響樂(lè)20-實(shí)戰(zhàn)三 10X 未過(guò)濾的PBMC數(shù)據(jù)
單細(xì)胞交響樂(lè)21-實(shí)戰(zhàn)三 批量處理并整合多個(gè)10X PBMC數(shù)據(jù)
單細(xì)胞交響樂(lè)22-實(shí)戰(zhàn)五 CEL-seq2
單細(xì)胞交響樂(lè)23-實(shí)戰(zhàn)六 CEL-seq
單細(xì)胞交響樂(lè)24-實(shí)戰(zhàn)七 SMARTer 胰腺細(xì)胞
單細(xì)胞交響樂(lè)25-實(shí)戰(zhàn)八 Smart-seq2 胰腺細(xì)胞
單細(xì)胞交響樂(lè)26-實(shí)戰(zhàn)九 胰腺細(xì)胞數(shù)據(jù)整合
單細(xì)胞交響樂(lè)27-實(shí)戰(zhàn)十 CEL-seq-小鼠造血干細(xì)胞
1 前言
前面的種種都是作為知識(shí)儲(chǔ)備，但是不實(shí)戰(zhàn)還是記不住前面的知識(shí)
這是第十一個(gè)實(shí)戰(zhàn)練習(xí)
數(shù)據(jù)來(lái)自 (Nestorowa et al. 2016) 的小鼠造血干細(xì)胞 haematopoietic stem cell (HSC) ，使用的技術(shù)是Smart-seq2
準(zhǔn)備數(shù)據(jù)
library(scRNAseq)
sce.nest <- NestorowaHSCData()
sce.nest
# class: SingleCellExperiment
# dim: 46078 1920
# metadata(0):
# assays(1): counts
# rownames(46078): ENSMUSG00000000001
# ENSMUSG00000000003 ... ENSMUSG00000107391
# ENSMUSG00000107392
# rowData names(0):
# colnames(1920): HSPC_007 HSPC_013 ... Prog_852
# Prog_810
# colData names(2): cell.type FACS
# reducedDimNames(1): diffusion
# altExpNames(1): ERCC
counts(sce.nest)[1:3,1:3]
# 3 x 3 sparse Matrix of class "dgCMatrix"
# HSPC_007 HSPC_013 HSPC_019
# ENSMUSG00000000001 . 7 1
# ENSMUSG00000000003 . . .
# ENSMUSG00000000028 4 1 2
看到使用了ERCC、Ensembl ID
ID轉(zhuǎn)換
library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
anno <- select(ens.mm.v97, keys=rownames(sce.nest),
keytype="GENEID", columns=c("SYMBOL", "SEQNAME"))
# 這里全部對(duì)應(yīng)
> sum(is.na(anno$SYMBOL))
[1] 0
> sum(is.na(anno$SEQNAME))
[1] 0
# 接下來(lái)只需要匹配順序即可
rowData(sce.nest) <- anno[match(rownames(sce.nest), anno$GENEID),]
sce.nest
# class: SingleCellExperiment
# dim: 46078 1920
# metadata(0):
# assays(1): counts
# rownames(46078): ENSMUSG00000000001
# ENSMUSG00000000003 ... ENSMUSG00000107391
# ENSMUSG00000107392
# rowData names(3): GENEID SYMBOL SEQNAME
# colnames(1920): HSPC_007 HSPC_013 ... Prog_852
# Prog_810
# colData names(2): cell.type FACS
# reducedDimNames(1): diffusion
# altExpNames(1): ERCC
2 質(zhì)控
依然是備份一下条霜，把unfiltered數(shù)據(jù)主要用在質(zhì)控的探索上
unfiltered <- sce.nest
這里沒(méi)有線粒體基因，因此只能用ERCC計(jì)算過(guò)濾條件
library(scater)
stats <- perCellQCMetrics(sce.nest)
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent")
sce.nest <- sce.nest[,!qc$discard]
# 看下過(guò)濾的細(xì)胞
colSums(as.matrix(qc))
# low_lib_size low_n_features
# 146 28
# high_altexps_ERCC_percent discard
# 241 264
做個(gè)圖
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc$discard
gridExtra::grid.arrange(
plotColData(unfiltered, y="sum", colour_by="discard") +
scale_y_log10() + ggtitle("Total count"),
plotColData(unfiltered, y="detected", colour_by="discard") +
scale_y_log10() + ggtitle("Detected features"),
plotColData(unfiltered, y="altexps_ERCC_percent",
colour_by="discard") + ggtitle("ERCC percent"),
ncol=2
)
最后對(duì)數(shù)據(jù)進(jìn)行過(guò)濾
sce.nest <- sce.nest[,!qc$discard]
# 過(guò)濾前后
> dim(unfiltered);dim(sce.nest)
[1] 46078 1920
[1] 46078 1656
3 歸一化
library(scran)
set.seed(101000110)
clusters <- quickCluster(sce.nest)
sce.nest <- computeSumFactors(sce.nest, clusters=clusters)
sce.nest <- logNormCounts(sce.nest)
4 找高變異基因
使用基于ERCC的構(gòu)建模型方法
set.seed(00010101)
dec.nest <- modelGeneVarWithSpikes(sce.nest, "ERCC")
top.nest <- getTopHVGs(dec.nest, prop=0.1)
class(dec.nest)
# [1] "DFrame"
# attr(,"package")
# [1] "S4Vectors"
# 其中ERCC的信息就存儲(chǔ)在dec.nest的metadata中
curfit <- metadata(dec.nest)
class(curfit)
# [1] "list"
names(curfit)
# [1] "mean" "var" "trend" "std.dev"
length(unique(names(curfit$mean))) # 一共92個(gè)ERCC spike-in
# [1] 92
# 其中的mean涵亏、var就定義了橫縱坐標(biāo)
head(curfit$mean)
# ERCC-00002 ERCC-00003 ERCC-00004 ERCC-00009 ERCC-00012 ERCC-00013
# 14.91183375 11.27060119 13.31197197 11.94866319 0.02211546 0.21249156
head(curfit$var)
# ERCC-00002 ERCC-00003 ERCC-00004 ERCC-00009 ERCC-00012 ERCC-00013
# 0.02375131 0.29308411 0.05376959 0.41814635 0.14928826 1.08599155
然后把基因（黑點(diǎn)）宰睡、ERCC（紅點(diǎn)）、根據(jù)ERCC擬合的線（藍(lán)線）畫(huà)出來(lái)
plot(dec.nest$mean, dec.nest$total, pch=16, cex=0.5,
xlab="Mean of log-expression", ylab="Variance of log-expression")
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
points(curfit$mean, curfit$var, col="red")
5 降維聚類(lèi)
降維
set.seed(101010011)
sce.nest <- denoisePCA(sce.nest, technical=dec.nest, subset.row=top.nest)
sce.nest <- runTSNE(sce.nest, dimred="PCA")
# 檢查PC的數(shù)量
ncol(reducedDim(sce.nest, "PCA"))
## [1] 9
聚類(lèi)
snn.gr <- buildSNNGraph(sce.nest, use.dimred="PCA")
colLabels(sce.nest) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
table(colLabels(sce.nest))
##
## 1 2 3 4 5 6 7 8 9
## 203 472 258 175 142 229 20 83 74
作圖看看
plotTSNE(sce.nest, colour_by="label")
6 marker基因檢測(cè)
markers <- findMarkers(sce.nest, colLabels(sce.nest),
test.type="wilcox", direction="up", lfc=0.5,
row.data=rowData(sce.nest)[,"SYMBOL",drop=FALSE])
比如檢測(cè)一下cluster8：
chosen <- markers[['8']]
best <- chosen[chosen$Top <= 10,]
length(best)
# [1] 13
# 將cluster8與其他clusters對(duì)比的AUC結(jié)果提取出來(lái)
aucs <- getMarkerEffects(best, prefix="AUC")
rownames(aucs) <- best$SYMBOL
library(pheatmap)
pheatmap(aucs, color=viridis::plasma(100))
看到其中血紅蛋白相關(guān)基因（Hba1气筋、Hba2拆内、Hbb）、Car2裆悄、Hebp1基因上調(diào)矛纹，說(shuō)明clsuter8可能包含紅細(xì)胞前體細(xì)胞
7 細(xì)胞類(lèi)型注釋
將會(huì)使用內(nèi)置的參考注釋數(shù)據(jù)臂聋，
SingleR中就包含了一些內(nèi)置數(shù)據(jù)集光稼，大部分是bulk RNA-Seq或芯片數(shù)據(jù)中經(jīng)過(guò)篩選的細(xì)胞類(lèi)型。
準(zhǔn)備參考數(shù)據(jù)
library(SingleR)
mm.ref <- MouseRNAseqData()
mm.ref
# class: SummarizedExperiment
# dim: 21214 358
# metadata(0):
# assays(1): logcounts
# rownames(21214): Xkr4 Rp1 ... LOC100039574
# LOC100039753
# rowData names(0):
# colnames(358): ERR525589Aligned
# ERR525592Aligned ... SRR1044043Aligned
# SRR1044044Aligned
# colData names(3): label.main label.fine
# label.ont
進(jìn)行轉(zhuǎn)換
renamed <- sce.nest
# 參考數(shù)據(jù)集中使用的是symbol name孩等，這里也轉(zhuǎn)換一下
rownames(renamed) <- uniquifyFeatureNames(rownames(renamed),
rowData(sce.nest)$SYMBOL)
# 然后把我們的細(xì)胞在參考數(shù)據(jù)集中找對(duì)應(yīng)的細(xì)胞類(lèi)型
# 返回的pred結(jié)果是一個(gè)數(shù)據(jù)框艾君，每行是我們自己數(shù)據(jù)的一個(gè)細(xì)胞
pred <- SingleR(test=renamed, ref=mm.ref, labels=mm.ref$label.fine)
table(pred$labels)
#
# B cells Endothelial cells Erythrocytes Granulocytes Macrophages Monocytes NK cells T cells
# 61 1 1005 1 2 500 1 85
這里也看到cluster8與紅細(xì)胞更相近
歡迎關(guān)注我們的公眾號(hào)～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩，打造生信星球肄方，想讓它成為一個(gè)不拽術(shù)語(yǔ)冰垄、通俗易懂的生信知識(shí)平臺(tái)。需要幫助或提出意見(jiàn)請(qǐng)后臺(tái)留言或發(fā)送郵件到jieandze1314@gmail.com
