scRNA-跟着Bioconductor一步一步学习scRNA（三）

Sep 25, 2019 6 min read cnposts

刘小泽开始写于19.9.25 内容来自Bioconductor的scRNA教程第三讲: Analyzing single-cell RNA-seq data containing UMI counts 上一次我们用了 少量的永生型小鼠骨髓祖细胞系+Smart-seq2，这次会用不同的细胞类型、数量以及不同的技术和分析方法

官方更新于2019-05-03

1 前言

依旧是我的传统：不🙅‍♀️翻译 并结合自己的知识尝试加入自己的理解

这次使用的是小鼠大脑细胞(Zeisel et al. 2015)，这个数据集中包含大约3000个细胞，异质性很强（包含的类型比如：少突细胞oligodendrocytes, 小神经胶质细胞microglia 以及神经元 neurons）。细胞是通过Fluidigm C1微流控系统分离的(Pollen et al. 2014)，另外文库制备过程中加入了UMI。mRNA、线粒体基因以及spike-in的定量结果保存在http://linnarssonlab.org/cortex/

下面的代码会自己在当前目录创建一个文件夹，然后我们只要给定路径（注意：file.path的链接是可以多段的，会自动paste，方便以后修改某一段链接），然后bfcrpath会根据链接下载，并将结果文件的路径保存下来

# 使用最方便的方法下载、保存数据
library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
base.url <- file.path("https://storage.googleapis.com",
    "linnarsson-lab-www-blobs/blobs/cortex")
mRNA.path <- bfcrpath(bfc, file.path(base.url, 
    "expression_mRNA_17-Aug-2014.txt"))
mito.path <- bfcrpath(bfc, file.path(base.url, 
    "expression_mito_17-Aug-2014.txt"))
spike.path <- bfcrpath(bfc, file.path(base.url, 
    "expression_spikes_17-Aug-2014.txt"))

2 数据的设置与调整

2.1 数据拆分

下载后的文件大小是：

下载后的文件大小

挑一个小一点的线粒体数据，看看长什么样子
tmp <- read.csv('raw_data/bc8538f89f6b_expression_mito_17-Aug-2014.txt',sep = '\t')
dim(tmp)
## 44 3007
# 这个线粒体有44行，3007列，因此肯定不能直接看，否则电脑就会卡死(一般列数超过200列就要小心了)
可以用命令行来看一下

先看一下数据结构

可以清楚看到，这个文件中同时包含了注释信息（1-10行）和表达矩阵（11行以后）；另外其他两个文件也是如此

所以明确了任务：第一步拆开注释和表达矩阵
# 因为要对三个文件都做一样的处理，所以可以用一个函数来解决(注意：函数不是通用的，一定是根据数据结构来创造的)
readFormat <- function(infile) { 
   ## 首先对注释信息处理 
   # 提取前10行；列名变表头；第一列是空的，所以去掉第一列
    metadata <- read.delim(infile, stringsAsFactors=FALSE, header=FALSE, nrow=10)[,-1] 
   # 将tissue、group等信息放到行名 
    rownames(metadata) <- metadata[,1]
    metadata <- metadata[,-1]
   # 转置一下，让细胞变行，group、tissue等统计信息变列；保存为数据框
    metadata <- as.data.frame(t(metadata))
    
    ## 然后对表达矩阵处理 
    # 跳过前11行；将基因名变到行名；保存为矩阵
    counts <- read.delim(infile, stringsAsFactors=FALSE, 
        header=FALSE, row.names=1, skip=11)[,-1] 
    counts <- as.matrix(counts)
    
    # 最后将注释和表达矩阵存到一个list中
    return(list(metadata=metadata, counts=counts))
}
endo.data <- readFormat(mRNA.path)
spike.data <- readFormat(spike.path)
mito.data <- readFormat(mito.path)

2.2 数据检查

想一下，数据框三要素：行、列、值。这三个数据如果要放在一起分析，它们的行名都是基因或转录本名，这个肯定不一样；它们的值也肯定不一样；那么我们唯一要关心的就是它们都包含的这3000多个细胞的ID号顺序是不是一样。只有顺序一样的细胞ID，才能在后续分析中进行整合

identical(spike.data$metadata$cell_id,mito.data$metadata$cell_id)
## FALSE
identical(spike.data$metadata$cell_id,endo.data$metadata$cell_id)
## TRUE
发现线粒体的细胞ID顺序和另外两个不同，需要按照另外两个顺序重新排序

于是match函数再一次派上用场。match函数：想按谁进行匹配，谁就放在第一个位置， 例如：match(A,B) 就实现了B按照A排列

m <- match(endo.data$metadata$cell_id, mito.data$metadata$cell_id)
mito.data$metadata <- mito.data$metadata[m,]
mito.data$counts <- mito.data$counts[,m]

另外这个内源基因数据集中的基因名也需要做点调整【这个就需要仔细检查名字了】：

rownames(endo.data$counts)
# 看到其中某些基因存在很多行，用loc1、loc2表示：说明某个基因存在于很多位置
## BC002163_loc2
## BC002163_loc1

对于这种情况，这里作者的建议是：将它们合并汇总，表达量相加（只是为了方便数据解释）

raw.names <- sub("_loc[0-9]+$", "", rownames(endo.data$counts))
new.counts <- rowsum(endo.data$counts, group=raw.names, reorder=FALSE)
endo.data$counts <- new.counts
接着将三个表达矩阵整合起来，来构建sce对象
library(SingleCellExperiment)
all.counts <- rbind(endo.data$counts, mito.data$counts, spike.data$counts)
sce <- SingleCellExperiment(list(counts=all.counts), colData=endo.data$metadata)
dim(sce)
## [1] 19896  3005

有了细胞的注释，还需要添加基因注释

nrows <- c(nrow(endo.data$counts), nrow(mito.data$counts), nrow(spike.data$counts))
is.spike <- rep(c(FALSE, FALSE, TRUE), nrows)
is.mito <- rep(c(FALSE, TRUE, FALSE), nrows)
isSpike(sce, "Spike") <- is.spike

添加Ensembl ID

library(org.Mm.eg.db)
ensembl <- mapIds(org.Mm.eg.db, keys=rownames(sce), keytype="SYMBOL", column="ENSEMBL")
rowData(sce)$ENSEMBL <- ensembl
到此，sce对象构建完成
sce
## class: SingleCellExperiment 
## dim: 19896 3005 
## metadata(0):
## assays(1): counts
## rownames(19896): Tspan12 Tshz1 ... ERCC-00170 ERCC-00171
## rowData names(1): ENSEMBL
## colnames(3005): V3 V4 ... V3006 V3007
## colData names(10): tissue group # ... level1class level2class
## reducedDimNames(0):
## spikeNames(1): Spike

3 细胞质控

3.1 先作图检查

我们这里使用的数据集是作者已经过滤好的、已发表的。不过依然可以来检查一下细胞质量

library(scater)
sce <- calculateQCMetrics(sce, feature_controls=list(Mt=is.mito)) 
# 然后作图:主要看文库大小、基因数量、ERCC占比、mt占比
par(mfrow=c(2,2), mar=c(5.1, 4.1, 0.1, 0.1)) # mar指定c(bottom, left, top, right) 
hist(sce$total_counts/1e3, xlab="Library sizes (thousands)", main="", 
    breaks=20, col="grey80", ylab="Number of cells")
hist(sce$total_features_by_counts, xlab="Number of expressed genes", main="", 
    breaks=20, col="grey80", ylab="Number of cells")
hist(sce$pct_counts_Spike, xlab="ERCC proportion (%)",
    ylab="Number of cells", breaks=20, main="", col="grey80")
hist(sce$pct_counts_Mt, xlab="Mitochondrial proportion (%)", 
    ylab="Number of cells", breaks=20, main="", col="grey80")

细胞质控

  • **文库大小：**相比于上一次分析的416B数据集，文库大小低了一个数量级，上一次的结果是：

    这是因为这里统计的是UMI数量，而不是read count（UMI就是为了去除PCR扩增偏差的。一般一个基因对应多个UMI时，出现多个reads含有同一个UMI时，这里只计数一次）英文解释就是：Each transcript molecule can only produce one UMI count but can yield many reads after fragmentation

    416B数据集文库大小

  • spike-in：它的波动比416B 数据集有点大，很有可能是因为这次的数据集细胞类型较多，每种细胞的内源RNA数量差异导致的

3.2 过滤细胞

依然是去除小的细胞文库和少的基因数量，去除大的spike-in占比。另外有spike-in存在就可以不用线粒体占比

libsize.drop <- isOutlier(sce$total_counts, nmads=3, type="lower", log=TRUE)
feature.drop <- isOutlier(sce$total_features_by_counts, nmads=3, type="lower", log=TRUE)
spike.drop <- isOutlier(sce$pct_counts_Spike, nmads=3, type="higher")

然后进行过滤，看到基本保留了全部的细胞，说明这个数据原来的质控确实不错

sce <- sce[,!(libsize.drop | feature.drop | spike.drop)]
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop), 
    BySpike=sum(spike.drop), Remaining=ncol(sce))
##   ByLibSize ByFeature BySpike Remaining
## 1         8         3       8      2989

其实isOutlier可以设置不同的因子型批次信息，来细化过滤条件。但是这里简化一下流程，因为前面质控结果还不错。

4 细胞周期推断

需要注意的一点就是：使用Ensembl ID，也就是前面用mapIds和rowData(sce)$ENSEMBL得到的结果

library(scran)
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))
assignments <- cyclone(sce, mm.pairs, gene.names=rowData(sce)$ENSEMBL)
table(assignments$phase)
## 
##   G1  G2M    S 
## 2980    8    1
plot(assignments$score$G1, assignments$score$G2M, xlab="G1 score", ylab="G2/M score", pch=16)

细胞周期推断

可以看到，基本就在G1期。**但是，cyclone不是万能的。**毕竟本身的数据集和cyclone提供的训练数据集存在不同，因此这个分期结果需要注意一下。

训练数据集使用的是C1 SMARTer方法得到的数据，它考虑的是那种方法包含的一些特性；而这里的数据集利用了UMI counts（比如10X数据就是用了UMI），这种数据的特性就是：3’-end coverage only, no length bias, no amplification noise。另外，许多神经细胞类型应该停留在 G0期暂时停止分裂，这个时期和其他细胞时期应该区别很明显 (Coller, Sang, and Roberts 2006)。cyclone可能只是将属于G0期的那部分细胞分配给了临近的G1期

5 基因层面检查

5.1 检查高表达基因

fontsize <- theme(axis.text=element_text(size=12), axis.title=element_text(size=16))
plotHighestExprs(sce, n=50) + fontsize

top50高表达基因

5.2 平均表达量=》Gene abundance

前面提到，UMI的表达量要比read count低

ave.counts <- calcAverage(sce, use_size_factors=FALSE)
hist(log10(ave.counts), breaks=100, main="", col="grey",
    xlab=expression(Log[10]~"average count"))

平均表达量

这里设置的基因过滤条件宽松一些，平均表达量不为0即可（因为这种数据特点是：细胞数量多但基因数量不多）

rowData(sce)$ave.count <- ave.counts
to.keep <- ave.counts > 0
sce <- sce[to.keep,]
summary(to.keep)
##    Mode   FALSE    TRUE 
## logical       2   19894

6 对细胞进行Normalization

因为这个数据中细胞异质性比较强，因此在用内源基因计算文库size factors 之前，先提前聚类一下(Lun, Bach, and Marioni 2016)。这个操作的结果就是： Cells in each cluster are normalized separately, and the size factors are rescaled to be comparable across clusters.

提前聚类的作用就是：不需要像上一个 416B数据（第7.1部分）一样，假设大部分基因在总体是无差异的，只需要假设

scRNA
Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related