劉小澤寫(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ì)胞
單細(xì)胞交響樂(lè)28-實(shí)戰(zhàn)十一 Smart-seq2-小鼠造血干細(xì)胞
單細(xì)胞交響樂(lè)29-實(shí)戰(zhàn)十二 10X 小鼠嵌合體胚胎
1 前言
前面的種種都是作為知識(shí)儲(chǔ)備序宦，但是不實(shí)戰(zhàn)還是記不住前面的知識(shí)
這是第十三個(gè)實(shí)戰(zhàn)練習(xí)
數(shù)據(jù)來(lái)自Bach et al. (2017)，使用的是10X的妊娠期小鼠乳腺上皮細(xì)胞
數(shù)據(jù)準(zhǔn)備
library(scRNAseq)
sce.mam <- BachMammaryData(samples="G_1")
sce.mam
# class: SingleCellExperiment
# dim: 27998 2915
# metadata(0):
# assays(1): counts
# rownames: NULL
# rowData names(2): Ensembl Symbol
# colnames: NULL
# colData names(3): Barcode Sample Condition
# reducedDimNames(0):
# altExpNames(0):
數(shù)據(jù)初探
# 樣本信息
sapply(names(colData(sce.mam)), function(x) head(colData(sce.mam)[,x]))
# Barcode Sample Condition
# [1,] "AAACCTGAGGATGCGT-1" "G_1" "Gestation"
# [2,] "AAACCTGGTAGTAGTA-1" "G_1" "Gestation"
# [3,] "AAACCTGTCAGCATGT-1" "G_1" "Gestation"
# [4,] "AAACCTGTCGTCCGTT-1" "G_1" "Gestation"
# [5,] "AAACGGGCACGAAATA-1" "G_1" "Gestation"
# [6,] "AAACGGGCAGACGCTC-1" "G_1" "Gestation"
ID轉(zhuǎn)換
依然是整合行名 + 添加染色體信息
library(scater)
rownames(sce.mam) <- uniquifyFeatureNames(
rowData(sce.mam)$Ensembl, rowData(sce.mam)$Symbol)
library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
rowData(sce.mam)$SEQNAME <- mapIds(ens.mm.v97, keys=rowData(sce.mam)$Ensembl,
keytype="GENEID", column="SEQNAME")
# 總共有13個(gè)線(xiàn)粒體基因
sum(grepl("MT",rowData(sce.mam)$SEQNAME))
# [1] 13
2 質(zhì)控
依然是備份一下，把unfiltered數(shù)據(jù)主要用在質(zhì)控的探索上
unfiltered <- sce.mam
使用線(xiàn)粒體信息進(jìn)行過(guò)濾
is.mito <- rowData(sce.mam)$SEQNAME == "MT"
stats <- perCellQCMetrics(sce.mam, subsets=list(Mito=which(is.mito)))
qc <- quickPerCellQC(stats, percent_subsets="subsets_Mito_percent")
colSums(as.matrix(qc))
## low_lib_size low_n_features high_subsets_Mito_percent
## 0 0 143
## discard
## 143
作圖
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="subsets_Mito_percent",
colour_by="discard") + ggtitle("Mito percent"),
ncol=3
)
再看看線(xiàn)粒體含量與文庫(kù)大小的關(guān)系
plotColData(unfiltered, x="sum", y="subsets_Mito_percent",
colour_by="discard") + scale_x_log10()
最后過(guò)濾
dim(unfiltered);dim(sce.mam)
# [1] 27998 2915
# [1] 27998 2772
3 歸一化
使用去卷積的方法
library(scran)
set.seed(101000110)
clusters <- quickCluster(sce.mam)
sce.mam <- computeSumFactors(sce.mam, clusters=clusters)
sce.mam <- logNormCounts(sce.mam)
4 找高變異基因
這里由于是10X的數(shù)據(jù)派歌，所以會(huì)有UMI信息，因此可以用基于泊松分布的模型構(gòu)建方法
set.seed(00010101)
dec.mam <- modelGeneVarByPoisson(sce.mam)
top.mam <- getTopHVGs(dec.mam, prop=0.1)
最后做個(gè)圖
plot(dec.mam$mean, dec.mam$total, pch=16, cex=0.5,
xlab="Mean of log-expression", ylab="Variance of log-expression")
curfit <- metadata(dec.mam)
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
5 降維聚類(lèi)
降維
library(BiocSingular)
set.seed(101010011)
sce.mam <- denoisePCA(sce.mam, technical=dec.mam, subset.row=top.mam)
sce.mam <- runTSNE(sce.mam, dimred="PCA")
# 檢查PC的數(shù)量
ncol(reducedDim(sce.mam, "PCA"))
## [1] 15
聚類(lèi)
有一個(gè)很重要的參數(shù)是
k，含義是：the number of nearest neighbors used to construct the graph胶果。如果k設(shè)置越大匾嘱，得到的圖之間聯(lián)通程度越高，cluster也越大早抠。因此這個(gè)參數(shù)也是可以不斷嘗試的
我們這里由于細(xì)胞數(shù)量比較多霎烙，所以設(shè)置的k就比較大，得到的cluster就少而大
snn.gr <- buildSNNGraph(sce.mam, use.dimred="PCA", k=25)
colLabels(sce.mam) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
table(colLabels(sce.mam))
##
## 1 2 3 4 5 6 7 8 9 10
## 550 799 716 452 24 84 52 39 32 24
最后作圖
plotTSNE(sce.mam, colour_by="label")
歡迎關(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
