單細胞交響樂21-實戰(zhàn)四批量處理并整合多個10X PBMC數(shù)據(jù)

劉小澤寫于2020.7.19
為何取名叫“交響樂”扩氢？因為單細胞分析就像一個大樂團，需要各個流程的協(xié)同配合
單細胞交響樂1-常用的數(shù)據(jù)結(jié)構(gòu)SingleCellExperiment
單細胞交響樂2-scRNAseq從實驗到下游簡介
 單細胞交響樂3-細胞質(zhì)控
 單細胞交響樂4-歸一化
 單細胞交響樂5-挑選高變化基因
 單細胞交響樂6-降維
 單細胞交響樂7-聚類分群
 單細胞交響樂8-marker基因檢測
 單細胞交響樂9-細胞類型注釋
 單細胞交響樂9-細胞類型注釋
 單細胞交響樂10-數(shù)據(jù)集整合后的批次矯正
 單細胞交響樂11-多樣本間差異分析
 單細胞交響樂12-檢測Doublet
單細胞交響樂13-細胞周期推斷
 單細胞交響樂14-細胞軌跡推斷
 單細胞交響樂15-scRNA與蛋白豐度信息結(jié)合
 單細胞交響樂16-處理大型數(shù)據(jù)
單細胞交響樂17-不同單細胞R包的數(shù)據(jù)格式相互轉(zhuǎn)換
 單細胞交響樂18-實戰(zhàn)一 Smart-seq2
單細胞交響樂19-實戰(zhàn)二 STRT-Seq
單細胞交響樂20-實戰(zhàn)三 10X 未過濾的PBMC數(shù)據(jù)

1 前言

前面的種種都是作為知識儲備捷绒，但是不實戰(zhàn)還是記不住前面的知識
這是第四個實戰(zhàn)練習(xí)

這次使用的數(shù)據(jù)是來自Zheng et al. 2017的三個PBMC數(shù)據(jù)，而且這個數(shù)據(jù)是經(jīng)過前期過濾的

準備數(shù)據(jù)

library(TENxPBMCData)
all.sce <- list(
    pbmc3k=TENxPBMCData('pbmc3k'),
    pbmc4k=TENxPBMCData('pbmc4k'),
    pbmc8k=TENxPBMCData('pbmc8k')
)

all.sce
# $pbmc3k
# class: SingleCellExperiment 
# dim: 32738 2700 
# metadata(0):
#   assays(1): counts
# rownames(32738): ENSG00000243485 ENSG00000237613 ...
# ENSG00000215616 ENSG00000215611
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
# colnames: NULL
# colData names(11): Sample Barcode ... Individual
# Date_published
# reducedDimNames(0):
#   altExpNames(0):
#   
#   $pbmc4k
# class: SingleCellExperiment 
# dim: 33694 4340 
# metadata(0):
#   assays(1): counts
# rownames(33694): ENSG00000243485 ENSG00000237613 ...
# ENSG00000277475 ENSG00000268674
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
# colnames: NULL
# colData names(11): Sample Barcode ... Individual
# Date_published
# reducedDimNames(0):
#   altExpNames(0):
#   
#   $pbmc8k
# class: SingleCellExperiment 
# dim: 33694 8381 
# metadata(0):
#   assays(1): counts
# rownames(33694): ENSG00000243485 ENSG00000237613 ...
# ENSG00000277475 ENSG00000268674
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
# colnames: NULL
# colData names(11): Sample Barcode ... Individual
# Date_published
# reducedDimNames(0):
#   altExpNames(0):

多個數(shù)據(jù)集的批量處理证鸥，重點就是列表list和for循環(huán)的熟練使用片部，還有相關(guān)的apply家族函數(shù)籽懦。而且每個結(jié)果數(shù)據(jù)也要對應(yīng)放在一個新列表中财骨，比如下面質(zhì)控使用的stats <- high.mito <- list() ，就是新建了兩個空列表具被，然后把結(jié)果放進去

2 批量質(zhì)控

依然是備份一下玻募，把unfiltered數(shù)據(jù)主要用在質(zhì)控的探索上

unfiltered <- all.sce

還是先通過線粒體含量計算質(zhì)控結(jié)果，然后根據(jù)這個結(jié)果進行過濾一姿，一個for循環(huán)搞定

library(scater)
stats <- high.mito <- list()

for (n in names(all.sce)) {
  current <- all.sce[[n]]
  is.mito <- grep("MT", rowData(current)$Symbol_TENx)
  stats[[n]] <- perCellQCMetrics(current, subsets=list(Mito=is.mito))
  high.mito[[n]] <- isOutlier(stats[[n]]$subsets_Mito_percent, type="higher")
  all.sce[[n]] <- current[,!high.mito[[n]]]
}

看一下根據(jù)線粒體過濾的結(jié)果

> lapply(high.mito, summary)
$pbmc3k
   Mode   FALSE    TRUE 
logical    2609      91 

$pbmc4k
   Mode   FALSE    TRUE 
logical    4182     158 

$pbmc8k
   Mode   FALSE    TRUE 
logical    8157     224

批量作圖（也是把作圖結(jié)果放進list，方便后期批量導(dǎo)出）

qcplots <- list()

for (n in names(all.sce)) {
  current <- unfiltered[[n]]
  colData(current) <- cbind(colData(current), stats[[n]])
  current$discard <- high.mito[[n]]
  qcplots[[n]] <- plotColData(current, x="sum", y="subsets_Mito_percent",
                              colour_by="discard") + scale_x_log10()
}
# do.call也是list處理中的常用函數(shù)
do.call(gridExtra::grid.arrange, c(qcplots, ncol=3))

3 批量歸一化

這里使用的是最簡單的方法logNormCounts()跃惫，就是用某個細胞中每個基因或spike-in轉(zhuǎn)錄本的表達量除以這個細胞計算的size factor叮叹，最后還進行一個log轉(zhuǎn)換，得到一個新的矩陣：logcounts 【不過這個名字并不是很貼切爆存，只是因為拼寫簡單蛉顽，真實應(yīng)該是：log-transformed normalized expression values。而不是單純字面意思取個log】

all.sce <- lapply(all.sce, logNormCounts)

看到這里先较，可能會想携冤，為什么沒計算size factor就直接進行了logNormCounts？

前面提到的操作闲勺，一般都是：

# 常規(guī)方法
lib.sf.zeisel <- librarySizeFactors(sce.zeisel)
lib.sf.zeisel <- logNormCounts(lib.sf.zeisel)
# 去卷積方法
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
sce.pbmc <- logNormCounts(sce.pbmc)

看一下logNormCounts的幫助文檔就能明白了曾棕，邏輯很清楚：

函數(shù)默認的參數(shù)是：size_factors=NULL，如果沒有計算size factor更新給函數(shù)菜循，那么函數(shù)會執(zhí)行normalizeCounts 的操作
再來看normalizeCounts會有什么操作：如果沒有提供size factor翘地，它會根據(jù)數(shù)據(jù)類型去自己計算size factor
- 對于count矩陣和SummarizedExperiment數(shù)據(jù)類型，會通過librarySizeFactors計算
- 對于SingleCellExperiment這種數(shù)據(jù)類型癌幕，它會首先在數(shù)據(jù)中尋找size factor是否存在衙耕，如果找不到，也是會使用librarySizeFactors

也就是說勺远，如果我們不提前計算橙喘，就會自動幫我們用最簡單的librarySizeFactors計算，并添加到我們的數(shù)據(jù)中胶逢。正是因為我們只需要最簡單的方法厅瞎，所以才可以不提供。如果要使用去卷積方法宪塔，還是要自己先計算好

最后看下結(jié)果

> lapply(all.sce, function(x) summary(sizeFactors(x)))
$pbmc3k
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2338  0.7478  0.9262  1.0000  1.1571  6.6042 

$pbmc4k
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3155  0.7109  0.8903  1.0000  1.1272 11.0267 

$pbmc8k
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2963  0.7043  0.8772  1.0000  1.1177  6.7942

4 批量找高變異基因

同樣也是使用了最簡單的計算方法

library(scran)
all.dec <- lapply(all.sce, modelGeneVar)
all.hvgs <- lapply(all.dec, getTopHVGs, prop=0.1)

作圖

par(mfrow=c(1,3))
for (n in names(all.dec)) {
    curdec <- all.dec[[n]]
    plot(curdec$mean, curdec$total, pch=16, cex=0.5, main=n,
        xlab="Mean of log-expression", ylab="Variance of log-expression")
    curfit <- metadata(curdec)
  # 可視化一條線（下圖的藍線）磁奖，這條線指所有的基因都會存在的一種偏差
    curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
}

每個點表示一個基因
圖中藍線指的是：技術(shù)因素導(dǎo)致的偏差
縱坐標表示總偏差：它等于技術(shù)偏差+生物因素偏差

因此，要衡量一個基因的生物因素偏差大小某筐，就看對應(yīng)的縱坐標減去對應(yīng)的藍線的值

5 批量降維

這里將每個PBMC數(shù)據(jù)單獨進行降維比搭，而并沒有把它們混合起來再分析

關(guān)于降維方法，這里選擇是與PCA近似的SVD算法（singular value decomposition，奇異值分解）身诺，scater或scran都可以直接通過函數(shù)計算SVD蜜托，利用參數(shù)BSPARAM=傳遞一個BiocSingularParam對象到runPCA中

關(guān)于SVD與PCA：https://www.cnblogs.com/bjwu/p/9280492.html
SVD是一種矩陣分解方法，相當于因式分解霉赡，目的純粹就是將一個矩陣拆分成多個矩陣相乘的形式
對于稀疏矩陣來說橄务，SVD算法更適用，這樣對于大數(shù)據(jù)來說節(jié)省了很大空間

library(BiocSingular)
set.seed(10000)
# 這里的runPCA是一個降維的函數(shù)名字穴亏，它可以用本身的PCA算法蜂挪，還可以用SVD算法
all.sce <- mapply(FUN=runPCA, x=all.sce, subset_row=all.hvgs, 
    MoreArgs=list(ncomponents=25, BSPARAM=RandomParam()), 
    SIMPLIFY=FALSE)

set.seed(100000)
all.sce <- lapply(all.sce, runTSNE, dimred="PCA")

set.seed(1000000)
all.sce <- lapply(all.sce, runUMAP, dimred="PCA")

這里使用SVD其實還是為了幫助更好地進行PCA，支持4種方式：

ExactParam: exact SVD with runExactSVD.
IrlbaParam: approximate SVD with irlba via runIrlbaSVD.
RandomParam: approximate SVD with rsvd via runRandomSVD.
FastAutoParam: fast approximate SVD, chosen based on the matrix representation.

6 批量聚類

使用基于圖形的聚類嗓化，最基礎(chǔ)的想法是：我們首先構(gòu)建一個圖棠涮，其中每個節(jié)點都是一個細胞，它與高維空間中最鄰近的細胞相連刺覆。連線基于細胞之間的相似性計算權(quán)重严肪，權(quán)重越高，表示細胞間關(guān)系更密切谦屑。如果一群細胞之間的權(quán)重高于另一群細胞驳糯，那么這一群細胞就被當做一個群體 “communiity”。

for (n in names(all.sce)) {
    g <- buildSNNGraph(all.sce[[n]], k=10, use.dimred='PCA')
    clust <- igraph::cluster_walktrap(g)$membership
    colLabels(all.sce[[n]])  <- factor(clust)
}
# 看看各自分了多少群
lapply(all.sce, function(x) table(colLabels(x)))
## $pbmc3k
## 
##   1   2   3   4   5   6   7   8   9  10 
## 487 154 603 514  31 150 179 333 147  11 
## 
## $pbmc4k
## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13 
##  497  185  569  786  373  232   44 1023   77  218   88   54   36 
## 
## $pbmc8k
## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
## 1004  759 1073 1543  367  150  201 2067   59  154  244   67   76  285   20   15 
##   17   18 
##   64    9

作圖

# 還是先新建一個空列表氢橙，方便保存數(shù)據(jù)
all.tsne <- list()

for (n in names(all.sce)) {
    all.tsne[[n]] <- plotTSNE(all.sce[[n]], colour_by="label") + ggtitle(n)
}
do.call(gridExtra::grid.arrange, c(all.tsne, list(ncol=3)))

7 數(shù)據(jù)整合

前面只是批量進行了各個數(shù)據(jù)集的分析酝枢，現(xiàn)在要把它們整合起來再分析一下

找共有基因

# 先看看各自多少基因
> lapply(all.sce, function(x) length(rownames(x)))
$pbmc3k
[1] 32738

$pbmc4k
[1] 33694

$pbmc8k
[1] 33694
         
# 找共同基因
universe <- Reduce(intersect, lapply(all.sce, rownames))
> length(universe)
[1] 31232

對每個數(shù)據(jù)批量取子集

# 這個操作可以記下來，lapply批量取子集
all.sce2 <- lapply(all.sce, "[", i=universe,)

> lapply(all.sce2, dim)
$pbmc3k
[1] 31232  2609

$pbmc4k
[1] 31232  4182

$pbmc8k
[1] 31232  8157

# 同樣的充蓝，對找高變異基因的結(jié)果取子集
all.dec2 <- lapply(all.dec, "[", i=universe,)

把三個數(shù)據(jù)當做一個數(shù)據(jù)的三個批次隧枫，重新進行歸一化

library(batchelor)
normed.sce <- do.call(multiBatchNorm, all.sce2)

根據(jù)重新歸一化的結(jié)果，再次找HVGs

這次是把3個批次放在一起再找的

combined.dec <- do.call(combineVar, all.dec2)
combined.hvg <- getTopHVGs(combined.dec, n=5000)

對一個大數(shù)據(jù)進行降維

結(jié)合：單細胞交響樂10-數(shù)據(jù)集整合后的批次矯正中的【第4部分 MNN矯正】

fastMNN先進行PCA降維谓苟，以加速下面的聚類環(huán)節(jié)

set.seed(1000101)
merged.pbmc <- do.call(fastMNN, c(normed.sce, 
    list(subset.row=combined.hvg, BSPARAM=RandomParam())))

merged.pbmc
# class: SingleCellExperiment 
# dim: 5000 14948 
# metadata(2): merge.info pca.info
# assays(1): reconstructed
# rownames(5000): ENSG00000090382 ENSG00000163220 ...
# ENSG00000122068 ENSG00000011132
# rowData names(1): rotation
# colnames: NULL
# colData names(2): batch label
# reducedDimNames(1): corrected
# altExpNames(0):

# 這時就出現(xiàn)了批次信息
> table(merged.pbmc$batch)
pbmc3k pbmc4k pbmc8k 
  2609   4182   8157

檢查一下結(jié)果官脓，使用lost.var ，值越大表示丟失的真實生物異質(zhì)性越多

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.

metadata(merged.pbmc)$merge.info$lost.var
##         pbmc3k    pbmc4k   pbmc8k
## [1,] 7.003e-03 3.126e-03 0.000000
## [2,] 7.137e-05 5.125e-05 0.003003

對一個大數(shù)據(jù)進行聚類

g <- buildSNNGraph(merged.pbmc, use.dimred="corrected")
colLabels(merged.pbmc) <- factor(igraph::cluster_louvain(g)$membership)
table(colLabels(merged.pbmc), merged.pbmc$batch)
##     
##      pbmc3k pbmc4k pbmc8k
##   1     113    387    825
##   2     507    395    806
##   3     175    344    581
##   4     295    539   1018
##   5     346    638   1210
##   6      11      3      9
##   7      17     27    111
##   8      33    113    185
##   9     423    754   1546
##   10      4     36     67
##   11    197    124    221
##   12    150    180    293
##   13    327    588   1125
##   14     11     54    160

可視化

set.seed(10101010)
merged.pbmc <- runTSNE(merged.pbmc, dimred="corrected")
# 查看3個數(shù)據(jù)混合后的聚類結(jié)果涝焙，以及有沒有批次效應(yīng)
gridExtra::grid.arrange(
    plotTSNE(merged.pbmc, colour_by="label", text_by="label", text_colour="red"),
    plotTSNE(merged.pbmc, colour_by="batch")
)

歡迎關(guān)注我們的公眾號～_～　　
我們是兩個農(nóng)轉(zhuǎn)生信的小碩卑笨，打造生信星球，想讓它成為一個不拽術(shù)語仑撞、通俗易懂的生信知識平臺赤兴。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!

單細胞交響樂21-實戰(zhàn)四 批量處理并整合多個10X PBMC數(shù)據(jù)