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

往期系列

亨贝格实验室单细胞转录组数据分析（一）

亨贝格实验室单细胞转录组数据分析（二）

亨贝格实验室单细胞转录组数据分析（三）

亨贝格实验室单细胞转录组数据分析（四）

亨贝格实验室单细胞转录组数据分析（五）

亨贝格实验室单细胞转录组数据分析（六）

Hemberg-lab单细胞转录组数据分析（七） - 导入10X和SmartSeq2数据Tabula Muris

收藏|北大生信平台 “单细胞分析，染色质分析” 视频和PPT分享

该如何自学入门生物信息学

生物信息之程序学习

收藏|你想要的生信学习系列教程 - 宝典在手，生信无忧

表达矩阵质控

UMI表达定量（UMI）

UMI表达定量简介

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

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

董数据集

Yoav Gilad在Yoav Gilad 我们使用芝加哥大学实验室的3个不同来源的诱导多能性干细胞（iPSC）的数据集（http://jdblischak.github.io/singleCellSeq/analysis/）。细胞分选采用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（结构）是一个很好的工具，可以用来查看数据的结构组成。（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


str(umi)


## Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 10 slots

##   ..@ int_elementMetadata:Formal class 'DataFrame' [package "S4Vectors"] with 6 slots
##   .. .. ..@ rownames       : NULL
##   .. .. ..@ nrows          : int 18726
##   .. .. ..@ listData       :List of 3
##   .. .. .. ..$ is_spike_ERCC: logi [1:18726] FALSE FALSE FALSE FALSE FALSE FALSE ...
##   .. .. .. ..$ is_spike     : logi [1:18726] FALSE FALSE FALSE FALSE FALSE FALSE ...
##   .. .. .. ..$ is_spike_MT  : logi [1:18726] FALSE FALSE FALSE FALSE FALSE FALSE ...
##   .. .. ..@ elementType    : chr "ANY"
##   .. .. ..@ elementMetadata: NULL
##   .. .. ..@ metadata       : list()
##   ..@ int_colData        :Formal class 'DataFrame' [package "S4Vectors"] with 6 slots
##   .. .. ..@ rownames       : NULL
##   .. .. ..@ nrows          : int 864
##   .. .. ..@ listData       : Named list()
##   .. .. ..@ elementType    : chr "ANY"
##   .. .. ..@ elementMetadata: NULL
##   .. .. ..@ metadata       : list()
##   ..@ int_metadata       :List of 3
##   .. ..$ version          :Classes 'package_version', 'numeric_version'  hidden list of 1
##   .. .. ..$ : int [1:3] 1 4 1
##   .. ..$ spike_names      : chr [1:2] "ERCC" "MT"
##   .. ..$ size_factor_names: chr(0) 
##   ..@ reducedDims        :Formal class 'SimpleList' [package "S4Vectors"] with 4 slots
##   .. .. ..@ listData       : Named list()
##   .. .. ..@ elementType    : chr "ANY"
##   .. .. ..@ elementMetadata: NULL
##   .. .. ..@ metadata       : list()
##   ..@ rowRanges          :Formal class 'CompressedGRangesList' [package "GenomicRanges"] with 5 slots
##   .. .. ..@ unlistData     :Formal class 'GRanges' [package "GenomicRanges"] with 7 slots
##   .. .. .. .. ..@ seqnames       :Formal class 'Rle' [package "S4Vectors"] with 4 slots
##   .. .. .. .. .. .. ..@ values         : Factor w/ 0 levels: 
##   .. .. .. .. .. .. ..@ lengths        : int(0) 
##   .. .. .. .. .. .. ..@ elementMetadata: NULL
##   .. .. .. .. .. .. ..@ metadata       : list()
##   .. .. .. .. ..@ ranges         :Formal class 'IRanges' [package "IRanges"] with 6 slots
##   .. .. .. .. .. .. ..@ start          : int(0) 
##   .. .. .. .. .. .. ..@ width          : int(0) 
##   .. .. .. .. .. .. ..@ NAMES          : NULL
##   .. .. .. .. .. .. ..@ elementType    : chr "ANY"
##   .. .. .. .. .. .. ..@ elementMetadata: NULL
##   .. .. .. .. .. .. ..@ metadata       : list()
##   .. .. .. .. ..@ strand         :Formal class 'Rle' [package "S4Vectors"] with 4 slots
##   .. .. .. .. .. .. ..@ values         : Factor w/ 3 levels "+","-","*": 
##   .. .. .. .. .. .. ..@ lengths        : int(0) 
##   .. .. .. .. .. .. ..@ elementMetadata: NULL
##   .. .. .. .. .. .. ..@ metadata       : list()
##   .. .. .. .. ..@ seqinfo        :Formal class 'Seqinfo' [package "GenomeInfoDb"] with 4 slots
##   .. .. .. .. .. .. ..@ seqnames   : chr(0) 
##   .. .. .. .. .. .. ..@ seqlengths : int(0) 
##   .. .. .. .. .. .. ..@ is_circular: logi(0) 
##   .. .. .. .. .. .. ..@ genome     : chr(0) 
##   .. .. .. .. ..@ elementMetadata:Formal class 'DataFrame' [package "S4Vectors"] with 6 slots
##   .. .. .. .. .. .. ..@ rownames       : NULL
##   .. .. .. .. .. .. ..@ nrows          : int 0
##   .. .. .. .. .. .. ..@ listData       : Named list()
##   .. .. .. .. .. .. ..@ elementType    : chr "ANY"
##   .. .. .. .. .. .. ..@ elementMetadata: NULL
##   .. .. .. .. .. .. ..@ metadata       : list()
##   .. .. .. .. ..@ elementType    : chr "ANY"
##   .. .. .. .. ..@ metadata       : list()
##   .. .. ..@ elementMetadata:Formal class 'DataFrame' [package "S4Vectors"] with 6 slots
##   .. .. .. .. ..@ rownames       : NULL
##   .. .. .. .. ..@ nrows          : int 18726
##   .. .. .. .. ..@ listData       :List of 9
##   .. .. .. .. .. ..$ is_feature_control     : logi [1:18726] FALSE FALSE FALSE FALSE FALSE FALSE ...
##   .. .. .. .. .. ..$ is_feature_control_ERCC: logi [1:18726] FALSE FALSE FALSE FALSE FALSE FALSE ...
##   .. .. .. .. .. ..$ is_feature_control_MT  : logi [1:18726] FALSE FALSE FALSE FALSE FALSE FALSE ...
##   .. .. .. .. .. ..$ mean_counts            : num [1:18726] 0.2824 0.0301 2.6389 0.2384 0.0116 ...
##   .. .. .. .. .. ..$ log10_mean_counts      : num [1:18726] 0.108 0.0129 0.561 0.0929 0.005 ...
##   .. .. .. .. .. ..$ n_cells_by_counts      : int [1:18726] 201 24 728 178 10 11 20 632 798 19 ...
##   .. .. .. .. .. ..$ pct_dropout_by_counts  : num [1:18726] 76.7 97.2 15.7 79.4 98.8 ...
##   .. .. .. .. .. ..$ total_counts           : int [1:18726] 244 26 2280 206 10 11 21 1638 2873 19 ...
##   .. .. .. .. .. ..$ log10_total_counts     : num [1:18726] 2.39 1.43 3.36 2.32 1.04 ...
##   .. .. .. .. ..@ elementType    : chr "ANY"
##   .. .. .. .. ..@ elementMetadata: NULL
##   .. .. .. .. ..@ metadata       : list()
##   .. .. ..@ elementType    : chr "GRanges"
##   .. .. ..@ metadata       : list()
##   .. .. ..@ partitioning   :Formal class 'PartitioningByEnd' [package "IRanges"] with 5 slots
##   .. .. .. .. ..@ end            : int [1:18726] 0 0 0 0 0 0 0 0 0 0 ...
##   .. .. .. .. ..@ NAMES          : chr [1:18726] "ENSG00000237683" "ENSG00000187634" "ENSG00000188976" "ENSG00000187961" ...
##   .. .. .. .. ..@ elementType    : chr "ANY"
##   .. .. .. .. ..@ elementMetadata: NULL
##   .. .. .. .. ..@ metadata       : list()
##   ..@ colData            :Formal class 'DataFrame' [package "S4Vectors"] with 6 slots
##   .. .. ..@ rownames       : chr [1:864] "NA19098.r1.A01" "NA19098.r1.A02" "NA19098.r1.A03" "NA19098.r1.A04" ...
##   .. .. ..@ nrows          : int 864
##   .. .. ..@ listData       :List of 50
##   .. .. .. ..$ individual                                    : chr [1:864] "NA19098" "NA19098" "NA19098" "NA19098" ...
##   .. .. .. ..$ replicate                                     : chr [1:864] "r1" "r1" "r1" "r1" ...
##   .. .. .. ..$ well                                          : chr [1:864] "A01" "A02" "A03" "A04" ...
##   .. .. .. ..$ batch                                         : chr [1:864] "NA19098.r1" "NA19098.r1" "NA19098.r1" "NA19098.r1" ...
##   .. .. .. ..$ sample_id                                     : chr [1:864] "NA19098.r1.A01" "NA19098.r1.A02" "NA19098.r1.A03" "NA19098.r1.A04" ...
##   .. .. .. ..$ is_cell_control                               : logi [1:864] FALSE FALSE FALSE FALSE FALSE FALSE ...
##   .. .. .. ..$ total_features_by_counts                      : int [1:864] 8368 8234 7289 7985 8619 8659 8054 9429 8243 9407 ...
##   .. .. .. ..$ log10_total_features_by_counts                : num [1:864] 3.92 3.92 3.86 3.9 3.94 ...
##   .. .. .. ..$ total_counts                                  : int [1:864] 63322 63976 43630 53922 70887 68051 58069 89082 60053 87489 ...
##   .. .. .. ..$ log10_total_counts                            : num [1:864] 4.8 4.81 4.64 4.73 4.85 ...
##   .. .. .. ..$ pct_counts_in_top_50_features                 : num [1:864] 17.6 16.5 18.2 16.8 16.2 ...
##   .. .. .. ..$ pct_counts_in_top_100_features                : num [1:864] 25.2 24.4 26.1 24.4 24.1 ...
##   .. .. .. ..$ pct_counts_in_top_200_features                : num [1:864] 34.9 34.5 36 34.1 34 ...
##   .. .. .. ..$ pct_counts_in_top_500_features                : num [1:864] 50.2 50.5 51.5 49.8 50 ...
##   .. .. .. ..$ total_features_by_counts_endogenous           : int [1:864] 8324 8190 7248 7942 8573 8616 8009 9384 8201 9361 ...
##   .. .. .. ..$ log10_total_features_by_counts_endogenous     : num [1:864] 3.92 3.91 3.86 3.9 3.93 ...
##   .. .. .. ..$ total_counts_endogenous                       : int [1:864] 57252 58967 39420 49076 65244 63508 53456 83920 56280 82287 ...
##   .. .. .. ..$ log10_total_counts_endogenous                 : num [1:864] 4.76 4.77 4.6 4.69 4.81 ...
##   .. .. .. ..$ pct_counts_endogenous                         : num [1:864] 90.4 92.2 90.4 91 92 ...
##   .. .. .. ..$ pct_counts_in_top_50_features_endogenous      : num [1:864] 12.5 13 13.1 12.3 12.5 ...
##   .. .. .. ..$ pct_counts_in_top_100_features_endogenous     : num [1:864] 20 20.7 21 19.8 20.2 ...
##   .. .. .. ..$ pct_counts_in_top_200_features_endogenous     : num [1:864] 29.8 30.8 31.1 29.5 30.1 ...
##   .. .. .. ..$ pct_counts_in_top_500_features_endogenous     : num [1:864] 45.9 47.3 47.4 45.9 46.7 ...
##   .. .. .. ..$ total_features_by_counts_feature_control      : int [1:864] 44 44 41 43 46 43 45 45 42 46 ...
##   .. .. .. ..$ log10_total_features_by_counts_feature_control: num [1:864] 1.65 1.65 1.62 1.64 1.67 ...
##   .. .. .. ..$ total_counts_feature_control                  : int [1:864] 6070 5009 4210 4846 5643 4543 4613 5162 3773 5202 ...
##   .. .. .. ..$ log10_total_counts_feature_control            : num [1:864] 3.78 3.7 3.62 3.69 3.75 ...
##   .. .. .. ..$ pct_counts_feature_control                    : num [1:864] 9.59 7.83 9.65 8.99 7.96 ...
##   .. .. .. ..$ pct_counts_in_top_50_features_feature_control : num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ pct_counts_in_top_100_features_feature_control: num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ pct_counts_in_top_200_features_feature_control: num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ pct_counts_in_top_500_features_feature_control: num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ total_features_by_counts_ERCC                 : int [1:864] 31 31 28 30 33 30 32 32 29 33 ...
##   .. .. .. ..$ log10_total_features_by_counts_ERCC           : num [1:864] 1.51 1.51 1.46 1.49 1.53 ...
##   .. .. .. ..$ total_counts_ERCC                             : int [1:864] 1187 1277 1121 1240 1262 1308 1203 1217 1217 1339 ...
##   .. .. .. ..$ log10_total_counts_ERCC                       : num [1:864] 3.07 3.11 3.05 3.09 3.1 ...
##   .. .. .. ..$ pct_counts_ERCC                               : num [1:864] 1.87 2 2.57 2.3 1.78 ...
##   .. .. .. ..$ pct_counts_in_top_50_features_ERCC            : num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ pct_counts_in_top_100_features_ERCC           : num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ pct_counts_in_top_200_features_ERCC           : num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ pct_counts_in_top_500_features_ERCC           : num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ total_features_by_counts_MT                   : int [1:864] 13 13 13 13 13 13 13 13 13 13 ...
##   .. .. .. ..$ log10_total_features_by_counts_MT             : num [1:864] 1.15 1.15 1.15 1.15 1.15 ...
##   .. .. .. ..$ total_counts_MT                               : int [1:864] 4883 3732 3089 3606 4381 3235 3410 3945 2556 3863 ...
##   .. .. .. ..$ log10_total_counts_MT                         : num [1:864] 3.69 3.57 3.49 3.56 3.64 ...
##   .. .. .. ..$ pct_counts_MT                                 : num [1:864] 7.71 5.83 7.08 6.69 6.18 ...
##   .. .. .. ..$ pct_counts_in_top_50_features_MT              : num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ pct_counts_in_top_100_features_MT             : num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ pct_counts_in_top_200_features_MT             : num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. .. ..$ pct_counts_in_top_500_features_MT             : num [1:864] 100 100 100 100 100 100 100 100 100 100 ...
##   .. .. ..@ elementType    : chr "ANY"
##   .. .. ..@ elementMetadata: NULL
##   .. .. ..@ metadata       : list()
##   ..@ assays             :Reference class 'ShallowSimpleListAssays' [package "SummarizedExperiment"] with 1 field
##   .. ..$ data: NULL
##   .. ..and 14 methods.
##   ..@ NAMES              : NULL
##   ..@ elementMetadata    :Formal class 'DataFrame' [package "S4Vectors"] with 6 slots
##   .. .. ..@ rownames       : NULL
##   .. .. ..@ nrows          : int 18726
##   .. .. ..@ listData       : Named list()
##   .. .. ..@ elementType    : chr "ANY"
##   .. .. ..@ elementMetadata: NULL
##   .. .. ..@ metadata       : list()
##   ..@ metadata           : list()