Hemberg-lab单细胞转录组数据分析（八）- Scater包输入导入和存储

生信宝典

发布于 2019-10-14 16:36:01

9030

发布于 2019-10-14 16:36:01

文章被收录于专栏：生信宝典生信宝典

表达矩阵质控

UMI表达定量 (UMI)

UMI表达定量简介

基因定量后会整理成一个行为基因（或转录本）列为细胞的表达矩阵。虽然前面做了原始数据质控和测序数据质控移除了一部分从reads数层面就不合格的细胞，还需要进一步根据表达矩阵移除其它类型低质量细胞。如果未能识别并移除低质量细胞会混淆下游分析中的有意义的生物信息。

鉴于scRNASeq还没有标准方法，后续用到的质控质量值的标准会因为实验方法不同而有很大变化。因此，执行质控时，我们是通过数据集内部比较找到异常细胞，而不是依赖于其它独立的质量标准。因此比较不同的建库方法获得的不同数据集时需要格外注意。

Tung数据集

我们使用芝加哥大学Yoav Gilad实验室的3个不同来源的诱导多能性干细胞 (iPSC)的数据集 (http://jdblischak.github.io/singleCellSeq/analysis/) in Yoav Gilad。细胞分选采用Fluidigm C1微流控台，同时使用UMIs和ERCC spike in进行质控为了保证可重复性，数据是2016年3月15生成的原始数据的拷贝，存储于tung文件夹下。

library(SingleCellExperiment)
library(scater)
options(stringsAsFactors = FALSE)

读入数据和注释：

molecules <- read.table("tung/molecules.txt", sep = "\t")
anno <- read.table("tung/annotation.txt", sep = "\t", header = TRUE)

查看读入的数据：

head(molecules[ , 1:3])

结果

##                 NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683              0              0              0
## ENSG00000187634              0              0              0
## ENSG00000188976              3              6              1
## ENSG00000187961              0              0              0
## ENSG00000187583              0              0              0
## ENSG00000187642              0              0              0

查看注释数据

head(anno)

结果

##   individual replicate well      batch      sample_id
## 1    NA19098        r1  A01 NA19098.r1 NA19098.r1.A01
## 2    NA19098        r1  A02 NA19098.r1 NA19098.r1.A02
## 3    NA19098        r1  A03 NA19098.r1 NA19098.r1.A03
## 4    NA19098        r1  A04 NA19098.r1 NA19098.r1.A04
## 5    NA19098        r1  A05 NA19098.r1 NA19098.r1.A05
## 6    NA19098        r1  A06 NA19098.r1 NA19098.r1.A06

数据包含 3个个体，3次重复和9个批次。

通过使用SingleCellExperiment (SCE) 和scater包标准化分析过程。首先构建SCE对象:

umi <- SingleCellExperiment(
    assays = list(counts = as.matrix(molecules)), 
    colData = anno
)

移除不在任何细胞表达的基因:

keep_feature <- rowSums(counts(umi) > 0) > 0
umi <- umi[keep_feature, ]

定义对照特征基因集：ERCC spike-ins和线粒体基因 (作者提供，见http://jdblischak.github.io/singleCellSeq/analysis/qc-filter-ipsc.html):

isSpike(umi, "ERCC") <- grepl("^ERCC-", rownames(umi))
isSpike(umi, "MT") <- rownames(umi) %in% 
    c("ENSG00000198899", "ENSG00000198727", "ENSG00000198888",
    "ENSG00000198886", "ENSG00000212907", "ENSG00000198786",
    "ENSG00000198695", "ENSG00000198712", "ENSG00000198804",
    "ENSG00000198763", "ENSG00000228253", "ENSG00000198938",
    "ENSG00000198840")

计算质量矩阵:

umi <- calculateQCMetrics(
    umi,
    feature_controls = list(
        ERCC = isSpike(umi, "ERCC"), 
        MT = isSpike(umi, "MT")
    )
)
str(umi)

SingleCellExperiment对象结构。str (structure)是一个很好的工具，可以用来查看数据的结构组成。(RStudio中的View对于较大的对象会给出更好的展现方式。)

如下面的展示SingleCellExperiment有10个数据槽 (slots)，使用@索引。比如想获得降维后的结果，使用umi@reduceDims，使用umi@colData可以获取质控信息，都有哪些质控变量，进一步的使用umi@colData@listData$total_features_by_counts可以获得每个细胞的基因数。umi@rowRanges@elementMetadata@listData基因的属性信息。

(在Rstudio中有个便利的地方，输入变量名再输@, $可以弹出对应的子属性，方便快速输入。)

# 获取原始读入的表达矩阵

head(umi@assays$data@listData$counts)[,1:3]
                NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
ENSG00000237683              0              0              0
ENSG00000187634              0              0              0
ENSG00000188976              3              6              1
ENSG00000187961              0              0              0
ENSG00000187583              0              0              0
ENSG00000187642              0              0              0

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2019-03-05，如有侵权请联系 cloudcommunity@tencent.com 删除

http