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

往期系列

表达矩阵质控

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()

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 202,406评论 5赞 475
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 84,976评论 2赞 379
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 149,302评论 0赞 335
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,366评论 1赞 273
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,372评论 5赞 363
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,457评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,872评论 3赞 395
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,521评论 0赞 256
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,717评论 1赞 295
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,523评论 2赞 319
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,590评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,299评论 4赞 318
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,859评论 3赞 306
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,883评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,127评论 1赞 259
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 42,760评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,290评论 2赞 342