劉小澤寫于2020.7.20
為何取名叫“交響樂”颤专？因?yàn)閱渭?xì)胞分析就像一個(gè)大樂團(tuán)，需要各個(gè)流程的協(xié)同配合
單細(xì)胞交響樂1-常用的數(shù)據(jù)結(jié)構(gòu)SingleCellExperiment
單細(xì)胞交響樂2-scRNAseq從實(shí)驗(yàn)到下游簡介
單細(xì)胞交響樂3-細(xì)胞質(zhì)控
單細(xì)胞交響樂4-歸一化
單細(xì)胞交響樂5-挑選高變化基因
單細(xì)胞交響樂6-降維
單細(xì)胞交響樂7-聚類分群
單細(xì)胞交響樂8-marker基因檢測
單細(xì)胞交響樂9-細(xì)胞類型注釋
單細(xì)胞交響樂9-細(xì)胞類型注釋
單細(xì)胞交響樂10-數(shù)據(jù)集整合后的批次矯正
單細(xì)胞交響樂11-多樣本間差異分析
單細(xì)胞交響樂12-檢測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
單細(xì)胞交響樂19-實(shí)戰(zhàn)二 STRT-Seq
單細(xì)胞交響樂20-實(shí)戰(zhàn)三 10X 未過濾的PBMC數(shù)據(jù)
單細(xì)胞交響樂21-實(shí)戰(zhàn)三 批量處理并整合多個(gè)10X PBMC數(shù)據(jù)
1 前言
前面的種種都是作為知識儲備围俘，但是不實(shí)戰(zhàn)還是記不住前面的知識
這是第五個(gè)實(shí)戰(zhàn)練習(xí)
第一代CEL-Seq由以色列理工學(xué)院研究團(tuán)隊(duì)于2012發(fā)表在Cell Reports
CEL-Seq的全稱是：Cell expression by linear amplification and sequencing胧奔，采用線性擴(kuò)增的測序方法哆姻，其主要優(yōu)勢在于錯(cuò)誤率比較低，但是和PCR一樣都存在序列偏向性诲祸。另外還具有以下優(yōu)點(diǎn)：
- 使用barcode允許多種類型的細(xì)胞同時(shí)分析；
- 每個(gè)細(xì)胞在一個(gè)管中進(jìn)行處理憨愉，避免了樣本之間的污染
后來在2016年烦绳，CEL-seq2誕生，與早期版本相比具有低成本配紫、高靈敏度径密，并縮短了實(shí)驗(yàn)操作實(shí)驗(yàn)時(shí)間
數(shù)據(jù)準(zhǔn)備
這里使用的數(shù)據(jù)是： Grun et al. (2016) 中不同人類供體的胰腺細(xì)胞
library(scRNAseq)
sce.grun <- GrunPancreasData()
sce.grun
# class: SingleCellExperiment
# dim: 20064 1728
# metadata(0):
# assays(1): counts
# rownames(20064): A1BG-AS1__chr19 A1BG__chr19 ...
# ZZEF1__chr17 ZZZ3__chr1
# rowData names(2): symbol chr
# colnames(1728): D2ex_1 D2ex_2 ... D17TGFB_95
# D17TGFB_96
# colData names(2): donor sample
# reducedDimNames(0):
# altExpNames(1): ERCC
rowData(sce.grun)
# DataFrame with 20064 rows and 2 columns
# symbol chr
# <character> <character>
# A1BG-AS1__chr19 A1BG-AS1 chr19
# A1BG__chr19 A1BG chr19
# A1CF__chr10 A1CF chr10
# A2M-AS1__chr12 A2M-AS1 chr12
# A2ML1__chr12 A2ML1 chr12
# ... ... ...
# ZYG11A__chr1 ZYG11A chr1
# ZYG11B__chr1 ZYG11B chr1
# ZYX__chr7 ZYX chr7
# ZZEF1__chr17 ZZEF1 chr17
# ZZZ3__chr1 ZZZ3 chr1
ID轉(zhuǎn)換
先將symbol轉(zhuǎn)為Ensembl ID
library(org.Hs.eg.db)
gene.ids <- mapIds(org.Hs.eg.db, keys=rowData(sce.grun)$symbol,
keytype="SYMBOL", column="ENSEMBL")
接下來有兩種處理方式：
第一種：將Ensembl ID加到原來數(shù)據(jù)中
library(scater)
rowData(sce.grun)$ensembl <- uniquifyFeatureNames(
rowData(sce.grun)$symbol, gene.ids)
rowData(sce.grun)
# DataFrame with 20064 rows and 3 columns
# symbol chr ensembl
# <character> <character> <character>
# A1BG-AS1__chr19 A1BG-AS1 chr19 ENSG00000268895
# A1BG__chr19 A1BG chr19 ENSG00000121410
# A1CF__chr10 A1CF chr10 ENSG00000148584
# A2M-AS1__chr12 A2M-AS1 chr12 ENSG00000245105
# A2ML1__chr12 A2ML1 chr12 ENSG00000166535
# ... ... ... ...
# ZYG11A__chr1 ZYG11A chr1 ENSG00000203995
# ZYG11B__chr1 ZYG11B chr1 ENSG00000162378
# ZYX__chr7 ZYX chr7 ENSG00000159840
# ZZEF1__chr17 ZZEF1 chr17 ENSG00000074755
# ZZZ3__chr1 ZZZ3 chr1 ENSG00000036549
第二種：將沒有匹配的NA去掉，并且去掉重復(fù)的行（因?yàn)榭赡艽嬖谝粋€(gè)symbol對應(yīng)多個(gè)Ensembl ID的情況）【這里就試一下這種方法】
keep <- !is.na(gene.ids) & !duplicated(gene.ids)
> sum(is.na(gene.ids));sum(duplicated(gene.ids))
[1] 2487
[1] 2515
sce.grun <- sce.grun[keep,]
rownames(sce.grun) <- gene.ids[keep]
> dim(sce.grun) # 過濾掉2000多基因
[1] 17548 1728
2 質(zhì)控
前幾次實(shí)戰(zhàn)的質(zhì)控都很順利躺孝，但并不意味著實(shí)際操作中這一步就可以跑跑代碼解決
這一次享扔，就遇到了一個(gè)常規(guī)代碼解決不了的問題，需要思考+調(diào)整
依然是備份一下植袍，把unfiltered數(shù)據(jù)主要用在質(zhì)控的探索上
unfiltered <- sce.grun
看到這里的數(shù)據(jù)中沒有線粒體基因
table(rowData(sce.grun)$chr)
#
# chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18
# 1766 672 1071 947 328 551 532 740 1050 256
# chr19 chr2 chr20 chr21 chr22 chr3 chr4 chr5 chr6 chr7
# 1245 1142 484 188 400 1014 693 791 907 825
# chr8 chr9 chrX chrY
# 604 664 658 20
倒是有幾個(gè)donor的批次信息
table(colData(sce.grun)$donor)
#
# D10 D17 D2 D3 D7
# 288 480 96 480 384
我們之前進(jìn)行質(zhì)控的時(shí)候惧眠，一般采用的先使用perCellQCMetrics 指定線粒體計(jì)算，然后用 quickPerCellQC 指定ERCC+線粒體進(jìn)行過濾的策略∮诟觯現(xiàn)在既然沒有線粒體的信息氛魁，那么在第一步的perCellQCMetrics就不需要加subset參數(shù)了
stats <- perCellQCMetrics(sce.grun)
# 參數(shù)subsets的含義: used to identify interesting feature subsets such as ERCC spike-in transcripts or mitochondrial genes.
# > colnames(stats)
# [1] "sum" "detected"
# [3] "percent_top_50" "percent_top_100"
# [5] "percent_top_200" "percent_top_500"
# [7] "altexps_ERCC_sum" "altexps_ERCC_detected"
# [9] "altexps_ERCC_percent" "total"
然后進(jìn)行過濾，并讓函數(shù)注意批次信息
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
batch=sce.grun$donor)
colSums(as.matrix(qc), na.rm=TRUE)
# low_lib_size low_n_features high_altexps_ERCC_percent
# 85 93 88
# discard
# 129
作圖
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc$discard
gridExtra::grid.arrange(
plotColData(unfiltered, x="donor", y="sum", colour_by="discard") +
scale_y_log10() + ggtitle("Total count"),
plotColData(unfiltered, x="donor", y="detected", colour_by="discard") +
scale_y_log10() + ggtitle("Detected features"),
plotColData(unfiltered, x="donor", y="altexps_ERCC_percent",
colour_by="discard") + ggtitle("ERCC percent"),
ncol=3
)
【重點(diǎn)】發(fā)現(xiàn)了問題
質(zhì)控并非是一帆風(fēng)順的厅篓，遇到問題應(yīng)該知道如何探索解決
看到ERCC那個(gè)圖秀存，很明顯D10和D3沒有被正確過濾，它們中高ERCC的細(xì)胞還是沒有被去掉
質(zhì)控過濾有一個(gè)前提條件：大部分的細(xì)胞都是高質(zhì)量的羽氮，但顯然這里的D10或链、D3不符合這個(gè)要求，因此我們需要再次只根據(jù)D17档押、D2澳盐、D7這三個(gè)“好一點(diǎn)的樣本”進(jìn)行過濾
# 重新計(jì)算過濾條件
qc2 <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
batch=sce.grun$donor,
subset=sce.grun$donor %in% c("D17", "D7", "D2"))
colSums(as.matrix(qc2), na.rm=TRUE)
# low_lib_size low_n_features high_altexps_ERCC_percent discard
# 450 511 606 665
這次再作圖看看
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc2$discard
gridExtra::grid.arrange(
plotColData(unfiltered, x="donor", y="sum", colour_by="discard") +
scale_y_log10() + ggtitle("Total count"),
plotColData(unfiltered, x="donor", y="detected", colour_by="discard") +
scale_y_log10() + ggtitle("Detected features"),
plotColData(unfiltered, x="donor", y="altexps_ERCC_percent",
colour_by="discard") + ggtitle("ERCC percent"),
ncol=3
)
看到這次對D10、D3的過濾力度大了很多令宿，基本滿意
最后把過濾條件應(yīng)用在原數(shù)據(jù)
sce.grun <- sce.grun[,!qc$discard]
3 歸一化
也是使用預(yù)分群+去卷積計(jì)算size factor的方法
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.grun)
sce.grun <- computeSumFactors(sce.grun, cluster=clusters)
sce.grun <- logNormCounts(sce.grun)
summary(sizeFactors(sce.grun))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.09581 0.50459 0.79505 1.00000 1.22212 12.16491
看看兩種歸一化方法的差異
plot(librarySizeFactors(sce.grun), sizeFactors(sce.grun), pch=16,
col=as.integer(factor(sce.grun$donor)),
xlab="Library size factors", ylab="Deconvolution factors", log="xy")
abline(a=0, b=1, col="red")
也是再一次證明了去卷積方法比常規(guī)的方法更加精細(xì)叼耙，對批次層面也會納入考量
4 找高變異基因
既然有批次的信息，那就加到構(gòu)建模型中去掀淘，而且有ERCC旬蟋，就選用modelGeneVarWithSpikes這個(gè)方法
block <- paste0(sce.grun$sample, "_", sce.grun$donor)
dec.grun <- modelGeneVarWithSpikes(sce.grun, spikes="ERCC", block=block)
top.grun <- getTopHVGs(dec.grun, prop=0.1)
5 矯正批次效應(yīng)
使用FastMNN方法：Correct for batch effects in single-cell expression data using a fast version of the mutual nearest neighbors (MNN) method.
library(batchelor)
set.seed(1001010)
merged.grun <- fastMNN(sce.grun, subset.row=top.grun, batch=sce.grun$donor)
merged.grun
# class: SingleCellExperiment
# dim: 1467 1063
# metadata(2): merge.info pca.info
# assays(1): reconstructed
# rownames(1467): ENSG00000172023 ENSG00000172016 ...
# ENSG00000144867 ENSG00000179820
# rowData names(1): rotation
# colnames(1063): D2ex_1 D2ex_2 ... D17TGFB_94
# D17TGFB_95
# colData names(1): batch
# reducedDimNames(2): corrected
# altExpNames(0):
檢查一下結(jié)果，使用lost.var 革娄，值越大表示丟失的真實(shí)生物異質(zhì)性越多
metadata(merged.grun)$merge.info$lost.var
# D10 D17 D2 D3 D7
# [1,] 0.030843209 0.030956515 0.000000000 0.00000000 0.00000000
# [2,] 0.007129994 0.011275730 0.036836491 0.00000000 0.00000000
# [3,] 0.003528589 0.005027722 0.006846244 0.05020422 0.00000000
# [4,] 0.012070814 0.014673302 0.013972521 0.01267586 0.05281354
- It contains a matrix of the variance lost in each batch (column) at each merge step (row).
- Large proportions of lost variance (>10%) suggest that correction is removing genuine biological heterogeneity.
6 降維 + 聚類
降維
set.seed(100111)
merged.grun <- runTSNE(merged.grun, dimred="corrected")
聚類
snn.gr <- buildSNNGraph(merged.grun, use.dimred="corrected")
colLabels(merged.grun) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
# 聚類的分群與批次之間的對比
table(Cluster=colLabels(merged.grun), Donor=merged.grun$batch)
## Donor
## Cluster D10 D17 D2 D3 D7
## 1 32 71 33 80 29
## 2 5 7 3 4 4
## 3 3 44 0 0 13
## 4 11 119 0 0 55
## 5 12 73 30 3 78
## 6 14 29 3 2 64
## 7 1 9 0 0 7
## 8 2 5 2 3 4
## 9 5 13 0 0 10
## 10 15 33 11 10 37
## 11 5 18 0 1 33
## 12 4 13 0 0 1
最后作圖看看批次效應(yīng)矯正如何
gridExtra::grid.arrange(
plotTSNE(merged.grun, colour_by="label"),
plotTSNE(merged.grun, colour_by="batch"),
ncol=2
)
歡迎關(guān)注我們的公眾號～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩倾贰，打造生信星球冕碟，想讓它成為一個(gè)不拽術(shù)語、通俗易懂的生信知識平臺匆浙。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com
