劉小澤寫于2020.7.19
為何取名叫“交響樂”埠啃？因?yàn)閱渭?xì)胞分析就像一個(gè)大樂團(tuán)，需要各個(gè)流程的協(xié)同配合
單細(xì)胞交響樂1-常用的數(shù)據(jù)結(jié)構(gòu)SingleCellExperiment
單細(xì)胞交響樂2-scRNAseq從實(shí)驗(yàn)到下游簡(jiǎn)介
單細(xì)胞交響樂3-細(xì)胞質(zhì)控
單細(xì)胞交響樂4-歸一化
單細(xì)胞交響樂5-挑選高變化基因
單細(xì)胞交響樂6-降維
單細(xì)胞交響樂7-聚類分群
單細(xì)胞交響樂8-marker基因檢測(cè)
單細(xì)胞交響樂9-細(xì)胞類型注釋
單細(xì)胞交響樂9-細(xì)胞類型注釋
單細(xì)胞交響樂10-數(shù)據(jù)集整合后的批次矯正
單細(xì)胞交響樂11-多樣本間差異分析
單細(xì)胞交響樂12-檢測(cè)Doublet
單細(xì)胞交響樂13-細(xì)胞周期推斷
單細(xì)胞交響樂14-細(xì)胞軌跡推斷
單細(xì)胞交響樂15-scRNA與蛋白豐度信息結(jié)合
單細(xì)胞交響樂16-處理大型數(shù)據(jù)
單細(xì)胞交響樂17-不同單細(xì)胞R包的數(shù)據(jù)格式相互轉(zhuǎn)換
單細(xì)胞交響樂18-實(shí)戰(zhàn)一 Smart-seq2
1 前言
前面的種種都是作為知識(shí)儲(chǔ)備，但是不實(shí)戰(zhàn)還是記不住前面的知識(shí)
這是第二個(gè)實(shí)戰(zhàn)練習(xí)
使用了一個(gè)存在異質(zhì)性的數(shù)據(jù)集，是研究小鼠大腦的 (Zeisel et al. 2015)
其中大約包含3000個(gè)細(xì)胞，包括少突膠質(zhì)細(xì)胞琉兜，小膠質(zhì)細(xì)胞和神經(jīng)元等，使用的細(xì)胞分離平臺(tái)是Fluidigm C1微流控系統(tǒng)毙玻，屬于比較早期的系統(tǒng)【單細(xì)胞測(cè)序的知識(shí)】
文庫(kù)制備時(shí)加入了UMI
UMI簡(jiǎn)單解釋：
UMI就是為了去除PCR擴(kuò)增偏差的豌蟋。一般一個(gè)基因?qū)?yīng)多個(gè)UMI時(shí)，出現(xiàn)多個(gè)reads含有同一個(gè)UMI時(shí)桑滩，這里只計(jì)數(shù)一次梧疲。UMI英文解釋：
Each transcript molecule can only produce one UMI count but can yield many reads after fragmentationUMI詳細(xì)解釋：
不管是bulk RNA還是scRNA，都需要進(jìn)行PCR擴(kuò)增运准，但是不可避免有一些轉(zhuǎn)錄本會(huì)被擴(kuò)增太多次幌氮，超過了真實(shí)表達(dá)量。當(dāng)起始文庫(kù)很小時(shí)（比如單細(xì)胞數(shù)據(jù)）胁澳，就需要更多次的PCR過程该互，這個(gè)次數(shù)越多，引入的誤差就越大韭畸。UMI就是Unique Molecular Identifier宇智，由4-10個(gè)隨機(jī)核苷酸組成，在mRNA反轉(zhuǎn)錄后胰丁，進(jìn)入到文庫(kù)中普筹，每一個(gè)mRNA隨機(jī)連上一個(gè)UMI，根據(jù)PCR結(jié)果可以計(jì)數(shù)不同的UMI隘马，最終統(tǒng)計(jì)mRNA的數(shù)量。UMI圖片解釋：
UMI有幾個(gè)要求：
- 不能是均聚物 妻顶，如AAAAAAAAAA
- 不能有N堿基
- 不能包含堿基質(zhì)量低于10的堿基
2 數(shù)據(jù)準(zhǔn)備
# 自己下載
library(scRNAseq)
sce.zeisel <- ZeiselBrainData()
# 或者使用之前分享的RData
load('sce.zeisel.RData')
sce.zeisel
# class: SingleCellExperiment
# dim: 20006 3005
# metadata(0):
# assays(1): counts
# rownames(20006): Tspan12 Tshz1 ... mt-Rnr1
# mt-Nd4l
# rowData names(1): featureType
# colnames(3005): 1772071015_C02 1772071017_G12
# ... 1772066098_A12 1772058148_F03
# colData names(10): tissue group # ...
# level1class level2class
# reducedDimNames(0):
# altExpNames(2): ERCC repeat
看到這么幾個(gè)信息：2萬(wàn)多基因酸员，3005個(gè)樣本蜒车；只有原始count矩陣；使用了symbol ID幔嗦；加入了ERCC
# 有57個(gè)ERCC
> dim(altExp(sce.zeisel,'ERCC'))
[1] 57 3005
# 而且已經(jīng)標(biāo)注了線粒體基因
> table(rowData(sce.zeisel))
endogenous mito
19972 34
一個(gè)重要操作：aggregateAcrossFeatures
英文解釋是：Sum together expression values (by default, counts) for each feature set in each cell.
但是只看說(shuō)明還是不好理解酿愧，舉個(gè)例子：
可以看到下面會(huì)有很多基因具有多個(gè)loc
比如OTTMUSG00000016609_loc4、OTTMUSG00000016609_loc3 其實(shí)可以算作一個(gè)基因
head(rownames(sce.zeisel)[grep("_loc[0-9]+$",rownames(sce.zeisel))])
# [1] "Syne1_loc2" "Hist1h2ap_loc1"
# [3] "Inadl_loc1" "OTTMUSG00000016609_loc4"
# [5] "OTTMUSG00000016609_loc3" "Gm5643_loc2"
# 這樣的有300多個(gè)
> length(grep("_loc[0-9]+$",rownames(sce.zeisel)))
[1] 330
如果拿一個(gè)基因來(lái)看
Syne1有Syne1_loc1和Syne1_loc2
> length(grep("Syne1",rownames(sce.zeisel)))
[1] 2
counts(sce.zeisel)[grep("Syne1",rownames(sce.zeisel)),][1:2,1:3]
# 1772071015_C02 1772071017_G12 1772071017_A05
# Syne1_loc2 11 2 4
# Syne1_loc1 0 0 4
如果使用這個(gè)函數(shù)邀泉，會(huì)有怎樣效果
test <- aggregateAcrossFeatures(sce.zeisel,
id=sub("_loc[0-9]+$", "", rownames(sce.zeisel)))
# 只剩一個(gè)了嬉挡，也就是合二為一
> length(grep("Syne1",rownames(test)))
[1] 1
# 看表達(dá)量，也是合二為一
> counts(test)[grep("Syne1",rownames(test)),][1:3]
1772071015_C02 1772071017_G12 1772071017_A05
11 2 8
因此汇恤，明白了庞钢，這個(gè)函數(shù)就是處理相同行：把幾個(gè)相同的行的值加在一起變?yōu)橐恍?/strong>
也就明白了，下面??為什么要進(jìn)行sub操作因谎，其實(shí)就是為了把loc去掉基括，暴露出相同的基因名，才能執(zhí)行aggregateAcrossFeatures函數(shù)
sce.zeisel <- aggregateAcrossFeatures(sce.zeisel,
id=sub("_loc[0-9]+$", "", rownames(sce.zeisel)))
> dim(sce.zeisel)
[1] 19839 3005
再添加Ensembl ID
library(org.Mm.eg.db)
rowData(sce.zeisel)$Ensembl <- mapIds(org.Mm.eg.db,
keys=rownames(sce.zeisel), keytype="SYMBOL", column="ENSEMBL")
3 質(zhì)控
還是備份一下财岔，把unfiltered數(shù)據(jù)主要用在質(zhì)控的探索上
unfiltered <- sce.zeisel
這個(gè)公共數(shù)據(jù)的作者在發(fā)表文章時(shí)將數(shù)據(jù)的低質(zhì)量細(xì)胞去掉了风皿，但并不妨礙我們做個(gè)質(zhì)控，也可以看看它去除的怎樣
stats <- perCellQCMetrics(sce.zeisel, subsets=list(
Mt=rowData(sce.zeisel)$featureType=="mito"))
qc <- quickPerCellQC(stats, percent_subsets=c("altexps_ERCC_percent",
"subsets_Mt_percent"))
sce.zeisel <- sce.zeisel[,!qc$discard]
> sum(qc$discard)
[1] 189
> dim(sce.zeisel)
[1] 19839 2816
根據(jù)原來(lái)的數(shù)據(jù)匠璧，加上質(zhì)控標(biāo)準(zhǔn)作圖
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc$discard
# 做個(gè)圖
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"),
plotColData(unfiltered, y="subsets_Mt_percent",
colour_by="discard") + ggtitle("Mito percent"),
ncol=2
)
再看下文庫(kù)大小和ERCC分別和線粒體含量的關(guān)系
gridExtra::grid.arrange(
plotColData(unfiltered, x="sum", y="subsets_Mt_percent",
colour_by="discard") + scale_x_log10(),
plotColData(unfiltered, x="altexps_ERCC_percent", y="subsets_Mt_percent",
colour_by="discard"),
ncol=2
)
然后檢查一下被過濾的原因
## low_lib_size low_n_features high_altexps_ERCC_percent
## 0 3 65
## high_subsets_Mt_percent discard
## 128 189
4 歸一化
這里細(xì)胞數(shù)量較多桐款，因此需要預(yù)先分群+去卷積計(jì)算size factor
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.zeisel)
sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clusters)
sce.zeisel <- logNormCounts(sce.zeisel)
summary(sizeFactors(sce.zeisel))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.119 0.486 0.831 1.000 1.321 4.509
看看兩種歸一化方法的差異
# 常規(guī)：最簡(jiǎn)單的只考慮文庫(kù)大小
summary(librarySizeFactors(sce.zeisel))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.1757 0.5680 0.8680 1.0000 1.2783 4.0839
plot(librarySizeFactors(sce.zeisel), sizeFactors(sce.zeisel), pch=16,
xlab="Library size factors", ylab="Deconvolution factors", log="xy")
5 找高變異基因
理論上，應(yīng)該對(duì)每個(gè)細(xì)胞都標(biāo)記批次信息夷恍，添加block信息魔眨。但是這里由于技術(shù)不同，每個(gè)板子上只有20-40個(gè)細(xì)胞裁厅，并且細(xì)胞群體具有高度異質(zhì)性冰沙，不能假設(shè)每個(gè)板上的細(xì)胞類型的分布是相同的，因此這里不使用block將批次信息“鎖住”是合適的
既然有ERCC执虹，就可以用第三種方法【在之前單細(xì)胞交響樂5-挑選高變化基因的2.3 考慮技術(shù)噪音】：
dec.zeisel <- modelGeneVarWithSpikes(sce.zeisel, "ERCC")
top.hvgs <- getTopHVGs(dec.zeisel, prop=0.1)
> length(top.hvgs)
[1] 1816
同樣使用了spike-in港庄，對(duì)比一下，看到這里不管總體方差還是技術(shù)因素方差都要比之前smart-seq2要小棵譬。
smart-seq2是按read計(jì)數(shù)凳宙，這里由于添加了UMI，是按molecule計(jì)數(shù)茬故，也就是說(shuō)盖灸，UMI的加入，確實(shí)減少了PCR擴(kuò)增的偏差影響
另外圖中看到磺芭，這里STRT-seq的spike-in方差一直要比內(nèi)源基因的方差小赁炎，也就是說(shuō)內(nèi)源基因的變化幅度一直保持高位，體現(xiàn)了數(shù)據(jù)中包含多種細(xì)胞類型而導(dǎo)致的異質(zhì)性钾腺，異質(zhì)性導(dǎo)致了基因表達(dá)極度不均衡
6 降維
library(BiocSingular)
set.seed(101011001)
sce.zeisel <- denoisePCA(sce.zeisel, technical=dec.zeisel, subset.row=top.hvgs)
sce.zeisel <- runTSNE(sce.zeisel, dimred="PCA")
看一下得到的PC數(shù)
ncol(reducedDim(sce.zeisel, "PCA"))
## [1] 50
7 聚類
snn.gr <- buildSNNGraph(sce.zeisel, use.dimred="PCA")
colLabels(sce.zeisel) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
看一下結(jié)果
table(colLabels(sce.zeisel))
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 283 451 114 143 599 167 191 128 350 70 199 58 39 24
畫圖
plotTSNE(sce.zeisel, colour_by="label")
8 找marker基因并解釋結(jié)果
主要還是關(guān)注上調(diào)基因徙垫，可以幫我們快速判斷出異質(zhì)性群體中各個(gè)細(xì)胞類型的差異
比如還是針對(duì)cluster1看看
markers <- findMarkers(sce.zeisel, direction="up")
marker.set <- markers[["1"]]
> ncol(marker.set)
[1] 17
> colnames(marker.set)
[1] "Top" "p.value" "FDR" "summary.logFC" "logFC.2" "logFC.3" "logFC.4"
[8] "logFC.5" "logFC.6" "logFC.7" "logFC.8" "logFC.9" "logFC.10" "logFC.11"
[15] "logFC.12" "logFC.13" "logFC.14"
head(marker.set[,1:8], 10)
使用cluster1的Top10基因（但不一定只是10個(gè)）畫熱圖
top.markers <- rownames(marker.set)[marker.set$Top <= 10]
> length(top.markers)
[1] 58
plotHeatmap(sce.zeisel, features=top.markers, order_columns_by="label")
接下來(lái)就是根據(jù)背景知識(shí)了讥裤，比如看到Gad1、Slc6a1表達(dá)量都很高姻报，可能表明cluster1屬于中間神經(jīng)元
另一種方法：基于logFC
比如可以挑出cluster1的計(jì)算結(jié)果marker.set中前50個(gè)基因（這里就是50個(gè)己英，而不是Top50），然后根據(jù)cluster1與其他clusters的logFC吴旋，對(duì)每個(gè)基因表達(dá)量做熱圖
library(pheatmap)
logFCs <- getMarkerEffects(marker.set[1:50,])
pheatmap(logFCs, breaks=seq(-5, 5, length.out=101))
那么這個(gè)函數(shù)到底做了什么呢损肛？看圖就知道：
也就是把每個(gè)基因在其他clusters的logFC結(jié)果挑出來(lái)，匯集成了一個(gè)新矩陣荣瑟，我們自己手動(dòng)也是可以做到的
歡迎關(guān)注我們的公眾號(hào)～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩治拿，打造生信星球，想讓它成為一個(gè)不拽術(shù)語(yǔ)褂傀、通俗易懂的生信知識(shí)平臺(tái)忍啤。需要幫助或提出意見請(qǐng)后臺(tái)留言或發(fā)送郵件到jieandze1314@gmail.com
