單細胞交響樂25-實戰(zhàn)八 Smart-seq2 胰腺細胞

劉小澤寫于2020.7.20
為何取名叫“交響樂”？因為單細胞分析就像一個大樂團撼港，需要各個流程的協(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ù)
單細胞交響樂21-實戰(zhàn)三批量處理并整合多個10X PBMC數(shù)據(jù)
單細胞交響樂22-實戰(zhàn)五 CEL-seq2
單細胞交響樂23-實戰(zhàn)六 CEL-seq
單細胞交響樂24-實戰(zhàn)七 SMARTer

1 前言

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

這是也是來自多個供體的人類胰腺細胞山上，使用Smart-seq2建庫技術(shù)套啤，數(shù)據(jù)來自Segerstolpe et al. (2016)

數(shù)據(jù)準(zhǔn)備

library(scRNAseq)
sce.seger <- SegerstolpePancreasData()
sce.seger
# class: SingleCellExperiment 
# dim: 26179 3514 
# metadata(0):
#   assays(1): counts
# rownames(26179): SGIP1 AZIN2 ... BIVM-ERCC5 eGFP
# rowData names(2): symbol refseq
# colnames(3514): HP1502401_N13 HP1502401_D14 ...
# HP1526901T2D_O11 HP1526901T2D_A8
# colData names(8): Source Name individual ... age
# body mass index
# reducedDimNames(0):
#   altExpNames(1): ERCC

看到3500多個細胞彤侍，包含ERCC肠缨，使用Symbol ID

看下樣本信息：

ID轉(zhuǎn)換

選擇的方式是：將沒有匹配的NA去掉，并且去掉重復(fù)的行

# 首先得到symbol ID和對應(yīng)的Ensembl ID（其中會存在無對應(yīng)的NA情況）
library(AnnotationHub)
edb <- AnnotationHub()[["AH73881"]]
symbols <- rowData(sce.seger)$symbol
ens.id <- mapIds(edb, keys=symbols, keytype="SYMBOL", column="GENEID")

# 之前見到的方法是：
# keep <- !is.na(gene.ids) & !duplicated(gene.ids)

# 這里使用了另一種方法（不是直接將NA去掉拥刻，而且替換成了symbol）
ens.id <- ifelse(is.na(ens.id), symbols, ens.id)
keep <- !duplicated(ens.id)

sce.seger <- sce.seger[keep,]
rownames(sce.seger) <- ens.id[keep]

小結(jié)一下：至此見到了三種ID轉(zhuǎn)換的方式怜瞒，根據(jù)最后保留的基因數(shù)量党窜，可以排個序：

保留基因最多（保留了NA和重復(fù)）：uniquifyFeatureNames
中等（保留了NA旅赢，去掉重復(fù)）：ifelse(is.na(ens.id), symbols, ens.id)
最少（去掉了NA以及重復(fù)）：!is.na(gene.ids) & !duplicated(gene.ids)

編輯樣本信息

之前有8列樣本的信息，有點冗余了喇闸。這里只保留3列關(guān)心的蒸眠，并重新命名

emtab.meta <- colData(sce.seger)[,c("cell type", 
    "individual", "single cell well quality")]
colnames(emtab.meta) <- c("CellType", "Donor", "Quality")
colData(sce.seger) <- emtab.meta

另外把細胞類型這一列中的“cell”字符去掉漾橙，并把首字母大寫

sce.seger$CellType <- gsub(" cell", "", sce.seger$CellType)
sce.seger$CellType <- paste0(
    toupper(substr(sce.seger$CellType, 1, 1)),
    substring(sce.seger$CellType, 2))

2 質(zhì)控

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

unfiltered <- sce.seger

之前作者在數(shù)據(jù)中已經(jīng)標(biāo)注了細胞質(zhì)量楞卡，可以看到有問題的細胞還是很多的：

table(sce.seger$Quality)
# 
# control, 2-cell well  control, empty well     low quality cell                   OK 
# 32                   96                 1177                 2209

因此就要注意了霜运，這里的數(shù)據(jù)會不會滿足“大部分細胞都是高質(zhì)量的”這個假設(shè)脾歇？

還是需要試一下，看看結(jié)果先

library(scater)
stats <- perCellQCMetrics(sce.seger)
qc1 <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
                     batch=sce.seger$Donor)


colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc1$discard

gridExtra::grid.arrange(
  plotColData(unfiltered, x="Donor", y="sum", colour_by="discard") +
    scale_y_log10() + ggtitle("Total count") +
    theme(axis.text.x = element_text(angle = 90)),
  plotColData(unfiltered, x="Donor", y="detected", colour_by="discard") +
    scale_y_log10() + ggtitle("Detected features") +
    theme(axis.text.x = element_text(angle = 90)),
  plotColData(unfiltered, x="Donor", y="altexps_ERCC_percent",
              colour_by="discard") + ggtitle("ERCC percent") +
    theme(axis.text.x = element_text(angle = 90)),
  ncol=3
)

看到HP1509101在過濾時存在過濾不完全的情況淘捡，HP1504901過濾的ERCC數(shù)量太多藕各，推測這兩個批次效果可能并不是很好，可能存在大量的低質(zhì)量細胞

因此焦除，再次指定subset 參數(shù)激况，重新畫圖

library(scater)
stats <- perCellQCMetrics(sce.seger)
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
    batch=sce.seger$Donor,
    subset=!sce.seger$Donor %in% c("HP1504901", "HP1509101"))

看看過濾掉多少

colSums(as.matrix(qc))
##              low_lib_size            low_n_features high_altexps_ERCC_percent 
##                       788                      1056                      1031 
##                   discard 
##                      1246

最后將qc過濾的與本來標(biāo)注低質(zhì)量的一同過濾

low.qual <- sce.seger$Quality == "low quality cell"
sce.seger <- sce.seger[,!(qc$discard | low.qual)]

# 過濾了大概1500個細胞
> dim(unfiltered);dim(sce.seger)
[1] 26179  3514
[1] 26179  2090

3 歸一化

此處會有一點小問題，值得注意膘魄！

本來有ERCC乌逐，操作應(yīng)該是：

library(scran)
sce.seger  = computeSpikeFactors(sce.seger, "ERCC")
sce.seger <- logNormCounts(sce.seger) 
# Error in .local(x, ...) : size factors should be positive

但由于存在幾個細胞中一個ERCC都沒有，所以會報錯
此時面臨兩個選擇：要么把這幾個細胞去掉创葡；要么就不借助ERCC浙踢，用另一種去卷積方法

> table(colSums(counts(altExp(sce.seger)))==0)

FALSE  TRUE 
 2087     3

如果要去掉這幾個細胞：

test=sce.seger[,!colSums(counts(altExp(sce.seger)))==0]
sce.test = computeSpikeFactors(test, "ERCC")
sce.test <- logNormCounts(test)

我們這里選擇保守的方法，不去掉細胞灿渴，使用另一種去卷積方法：

clusters <- quickCluster(sce.seger)
sce.seger <- computeSumFactors(sce.seger, clusters=clusters)
sce.seger <- logNormCounts(sce.seger) 

summary(sizeFactors(sce.seger))
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 0.0000  0.1832  0.4016  1.0000  1.0996 12.9607

4 找高變異基因

下面構(gòu)建模型想使用modelGeneVarWithSpikes洛波，于是首先應(yīng)該把那幾個沒有ERCC的細胞去掉；另外由于AZ這個批次相對其他批次的細胞數(shù)量過于少逻杖，因此在模型構(gòu)建中也把它去掉吧

for.hvg <- sce.seger[,librarySizeFactors(altExp(sce.seger)) > 0
    & sce.seger$Donor!="AZ"]
dec.seger <- modelGeneVarWithSpikes(for.hvg, "ERCC", block=for.hvg$Donor)
chosen.hvgs <- getTopHVGs(dec.seger, n=2000)

如果要批量作圖檢查的話

# 批次數(shù)量較多奋岁，因此設(shè)置多行多列顯示
par(mfrow=c(3,3))
blocked.stats <- dec.seger$per.block
for (i in colnames(blocked.stats)) {
    current <- blocked.stats[[i]]
    plot(current$mean, current$total, main=i, pch=16, cex=0.5,
        xlab="Mean of log-expression", ylab="Variance of log-expression")
    curfit <- metadata(current)
    points(curfit$mean, curfit$var, col="red", pch=16)
    curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
}

注意，這里在找完HVGs后荸百，沒有進行批次矯正闻伶，如果繼續(xù)向下做，會發(fā)現(xiàn)什么够话？

5 降維聚類

降維

library(BiocSingular)
set.seed(101011001)
sce.seger <- runPCA(sce.seger, subset_row=chosen.hvgs, ncomponents=25)
sce.seger <- runTSNE(sce.seger, dimred="PCA")

聚類

snn.gr <- buildSNNGraph(sce.seger, use.dimred="PCA")
colLabels(sce.seger) <- factor(igraph::cluster_walktrap(snn.gr)$membership)

檢查聚類分群與批次

tab <- table(Cluster=colLabels(sce.seger), Donor=sce.seger$Donor)
library(pheatmap)
pheatmap(log10(tab+10), color=viridis::viridis(100))

結(jié)果真的是：批次效應(yīng)影響了分群蓝翰，因此最好還是做一遍fastMNN操作

tSNE圖中也是顯示出了強烈的批次效應(yīng)

gridExtra::grid.arrange(
    plotTSNE(sce.seger, colour_by="label"),
    plotTSNE(sce.seger, colour_by="Donor"),
    ncol=2
)

6 補充矯正批次效應(yīng)

上圖看到很明顯的批次效應(yīng)，那么如果處理后女嘲，會有什么不同嗎畜份？

利用fastMNN矯正

library(batchelor)
set.seed(1001010)
merged.seger <- fastMNN(sce.seger, subset.row=chosen.hvgs, 
                         batch=sce.seger$Donor)
merged.seger
# class: SingleCellExperiment 
# dim: 2000 2090 
# metadata(2): merge.info pca.info
# assays(1): reconstructed
# rownames(2000): GCG TTR ... MAP6 LCP1
# rowData names(1): rotation
# colnames(2090): HP1502401_H13 HP1502401_J14 ...
# HP1526901T2D_N8 HP1526901T2D_A8
# colData names(2): batch label
# reducedDimNames(2): corrected TSNE
# altExpNames(0):

# metadata(merged.seger)$merge.info$lost.var
# lost.var ：值越大表示丟失的真實生物異質(zhì)性越多

因為fastMNN會包含PCA降維，所以下面繼續(xù)進行tSNE即可

降維聚類

library(BiocSingular)
set.seed(101011001)
merged.seger <- runTSNE(merged.seger, dimred="corrected")

snn.gr <- buildSNNGraph(merged.seger, use.dimred="corrected")
colLabels(merged.seger) <- factor(igraph::cluster_walktrap(snn.gr)$membership)

再次作圖欣尼，是不是明顯比之前好很多爆雹？

gridExtra::grid.arrange(
  plotTSNE(merged.seger, colour_by="label"),
  plotTSNE(merged.seger, colour_by="batch"),
  ncol=2
)

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

Welcome to our bioinfoplanet!