单细胞都有两个“对象”，你一个都没有

单细胞基础分析第二课：数据下载

大家好，我是风。欢迎来到风风的从零开始单细胞系列一。在上一讲我们已经介绍了如何在GEO数据库下载单细胞测序的数据集，有了数据，我们就可以正式进入分析环节。我们这个系列的内容只需要基于R语言，不涉及Linux和Python。一口吃不成个大胖子嘛，等到大家把系列一的内容消化了，有进一步的需要的话，我们再一起来学习Linux和Python的分析知识，当然我个人觉得如果你自己或者老板不考虑单细胞测序的话，学一学R语言的分析内容对整个流程有所了解也就可以了。

本次我们将从以下几个方面进行探讨：

使用工具：即单细胞分析过程中主要R包的介绍(当然这个过程不止这几个R包啦，还有一些常规的工具就没有介绍了)；
拟时序分析：对大家来说，聚类、降维等概念都常听到，但是拟时序分析可能比较陌生，所以给大家口语化地解释这个概念；
数据读取和构建分析对象：包括Smartseq2数据读取并构建scater对象和10x Genomics数据读取构建Seurat对象。注意：构建scater对象和构建Seurat对象是两个知识点，不管是Smartseq2或者10x Genomics，都可以构建这两种对象用于后续分析。

这里就涉及一个新的概念了——拟时序分析。

什么是拟时序分析呢？

拟时序分析，英文为pseudotime，前缀为pseudo，解释为假的、冒充的，那拟就是假的意思；时序分析好理解，时间序列分析嘛，集合到单细胞领域那就是描述细胞从一种状态转变成另一种状态的分析；ok，合起来拟时序分析就是：一种假的描述细胞从一种状态转变成另一种状态的分析。说假的不太好听，读书人嘛，哪来什么真的假的？那是“推断”，所以我们可以理解为：人为使用算法推断细胞从一种状态转变成另一种状态的分析，就叫“拟时序分析”。那什么是真时序分析呢？留给你们思考一下，留言区告诉我，等我们后面讲到dynverse的时候揭晓问题答案！

数据读取和构建分析对象

在下载完数据之后，我们需要使用R语言将数据读取，以便后续的分析。上一讲在数据下载的时候我们已经知道了单细胞分析的数据储存方式多种，如果我们按平台来区分，大概可以分为两种：Smartseq2和10x Genomics。同样地，两者对应的数据在读取时候也有一些小的区别，我们分别来看一下：

Smartseq2数据

为了操作的一致性，我把数据放在后台，大家回复关键词就可以获取。如果还有关于公共数据库单细胞数据下载的问题，请大家翻看上一讲的内容。接下来我们先安装对应的R包：

# BiocManager::install('SingleCellExperiment')# BiocManager::install('scater')# BiocManager::install('Seurat')

上面的代码去掉#号后运行就可以安装R包，这也是本节课我们会用到的3个R包，接下来读取数据：

rm(list=ls())# 读取数据dat <- read.delim("Smartseq2\\Kidney-counts.csv", sep=",", header=TRUE, row.names = 1)dat[1:3,1:3] #查看前3行3列## A14.MAA000545.3_8_M.1.1 E1.MAA000545.3_8_M.1.1## 0610005C13Rik 0 0## 0610007C21Rik 1 0## 0610007L01Rik 0 0## M4.MAA000545.3_8_M.1.1## 0610005C13Rik 0## 0610007C21Rik 0## 0610007L01Rik 0

观察数据，我们发现data的行名为基因名，由于这是一个Smartseq2数据集，它可能包含spike-in，所以我们需要对数据进行检查，并将列名分割，提取出我们要的元素：

# 查看数据的行列dim(dat) # 23433 865 表示有23433行，865列## [1] 23433 865# 检查spike-inspike_in <- rownames(dat)[grep("^ERCC-", rownames(dat))] length(spike_in) # 总共有92个## [1] 92# 提取元数据cellIDs <- colnames(dat)cell_info <- strsplit(cellIDs, "\\.") # 将列名按.进行分裂Well <- lapply(cell_info, function(x){x[1]}) # 保留列名第一个点前的元素Well <- unlist(Well) # 将list转为向量class(Well)## [1] "character"# 接下来把上面两步结合起来提取Plate和动物编号Plate <- unlist(lapply(cell_info, function(x){x[2]})) # 保留列名第一个点和第二个点之间的元素Mouse <- unlist(lapply(cell_info, function(x){x[3]})) # 保留列名第二个点和第三个点之间的元素summary(factor(Mouse)) # 检查数据的分类## 3_10_M 3_11_M 3_38_F 3_39_F 3_8_M 3_9_M ## 104 196 237 169 82 77table(Mouse, Plate)## Plate## Mouse B001717 B002775 MAA000545 MAA000752 MAA000801 MAA000922## 3_10_M 0 0 0 104 0 0## 3_11_M 0 0 0 0 196 0## 3_38_F 237 0 0 0 0 0## 3_39_F 0 169 0 0 0 0## 3_8_M 0 0 82 0 0 0## 3_9_M 0 0 0 0 0 77

结果显示数据主要来源于6只实验动物，并且给出了每只动物所占的列数。接下来我们就可以对矩阵进行注释了：

ann <- read.table("Smartseq2\\annotations_FACS.csv", sep=",", header=TRUE) # 读取注释信息head(ann) # 观察注释数据，发现注释信息第一列就是前面矩阵的列名## cell tissue cell_ontology_class## 1 A21.MAA000594.3_8_M.1.1 Aorta fibroblast## 2 F8.MAA000594.3_8_M.1.1 Aorta unknown## 3 H11.MAA000594.3_8_M.1.1 Aorta unknown## 4 A22.MAA000594.3_8_M.1.1 Aorta unknown## 5 H12.MAA000594.3_8_M.1.1 Aorta epicardial adipocyte## 6 L9.MAA000594.3_8_M.1.1 Aorta unknown## cell_ontology_term_iri cell_ontology_id## 1 http://purl.obolibrary.org/obo/CL_0000057 CL:0000057## 2 <NA> CL:.## 3 <NA> CL:.## 4 <NA> CL:.## 5 http://purl.obolibrary.org/obo/CL_1000309 CL:1000309## 6 <NA> CL:.ann <- ann[match(cellIDs, ann[,1]),] # 样本匹配,得到前面矩阵对应的注释信息celltype <- ann[,3] # 提取细胞类型

接下来我们就要构建一个scater对象了，也要用到了我们前面安装的R包：

library("SingleCellExperiment")## Loading required package: SummarizedExperiment## Loading required package: GenomicRanges## Loading required package: stats4## Loading required package: BiocGenerics## Loading required package: parallel## ## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':## ## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,## clusterExport, clusterMap, parApply, parCapply, parLapply,## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':## ## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':## ## anyDuplicated, append, as.data.frame, basename, cbind, colnames,## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,## union, unique, unsplit, which, which.max, which.min## Loading required package: S4Vectors## ## Attaching package: 'S4Vectors'## The following object is masked from 'package:base':## ## expand.grid## Loading required package: IRanges## ## Attaching package: 'IRanges'## The following object is masked from 'package:grDevices':## ## windows## Loading required package: GenomeInfoDb## Loading required package: Biobase## Welcome to Bioconductor## ## Vignettes contain introductory material; view with## 'browseVignettes()'. To cite Bioconductor, see## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Loading required package: DelayedArray## Loading required package: matrixStats## ## Attaching package: 'matrixStats'## The following objects are masked from 'package:Biobase':## ## anyMissing, rowMedians## Loading required package: BiocParallel## ## Attaching package: 'DelayedArray'## The following objects are masked from 'package:matrixStats':## ## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges## The following objects are masked from 'package:base':## ## aperm, apply, rowsumlibrary("scater")## Loading required package: ggplot2cell_anns <- data.frame(mouse = Mouse, well=Well, type=celltype)rownames(cell_anns) <- colnames(dat)head(cell_anns)## mouse well type## A14.MAA000545.3_8_M.1.1 3_8_M A14 kidney collecting duct cell## E1.MAA000545.3_8_M.1.1 3_8_M E1 kidney tubule cell## M4.MAA000545.3_8_M.1.1 3_8_M M4 kidney collecting duct cell## O21.MAA000545.3_8_M.1.1 3_8_M O21 kidney tubule cell## P4.MAA000545.3_8_M.1.1 3_8_M P4 kidney tubule cell## A16.MAA000545.3_8_M.1.1 3_8_M A16 kidney collecting duct cellsceset <- SingleCellExperiment(assays = list(counts = as.matrix(dat)), colData=cell_anns)head(sceset)## class: SingleCellExperiment ## dim: 6 865 ## metadata(0):## assays(1): counts## rownames(6): 0610005C13Rik 0610007C21Rik ... 0610007P08Rik## 0610007P14Rik## rowData names(0):## colnames(865): A14.MAA000545.3_8_M.1.1 E1.MAA000545.3_8_M.1.1 ...## N8.B002775.3_39_F.1.1 P7.B002775.3_39_F.1.1## colData names(3): mouse well type## reducedDimNames(0):## spikeNames(0):## altExpNames(0):# 还记得前面我们找出来的spike-ins吗？我们在scater对象中进行处理isSpike(sceset, "ERCC") <- grepl("ERCC-", rownames(sceset))## Warning: 'isSpike<-' is deprecated.## Use 'isSpike<-' instead.## See help("Deprecated")## Warning: 'spikeNames' is deprecated.## See help("Deprecated")## Warning: 'isSpike' is deprecated.## See help("Deprecated")

好了，到这里Smartseq2数据的读取就完成了，我知道你现在还有很多不懂，比如spike-ins？又比如为什么要对数据的列名进行处理？没关系，这些都不影响你对数据进行操作，等到学习了后面的内容你也自然就清楚了。我在后台文件中还提供了一份Liver的数据，所以今天作业之一就是用这份数据跑一遍上面的流程，看看有什么不一样？好好完成作业有惊喜哦！

10x Genomics数据

接下来到了10x数据的读取，10x Genomics数据作为当今单细胞测序的主流，当然也可以用我们上面说的方法来读取并且构建分析对象。但是，既然前面已经讲了一个知识点了，我也就不想再重复一次，这次我们不构建scater对象了，我们来构建更实用的Seurat对象，同样数据保存在后台，这份数据跟你从网站上下载的或者公司给到你手上的数据格式一致，也就是我们上一讲说到的第三种数据格式：

# 简单直接一行代码读取三个文件library(Seurat)## ## Attaching package: 'Seurat'## The following object is masked from 'package:SummarizedExperiment':## ## Assaysdat1 <- Read10X(data.dir = "10x\\seq") # 直接读取seq文件夹下三个文件

# 注意：seq下三个文件夹分别为matrix.mtx、genes.tsv和barcodes.tsv，名称和尾缀的数据格式都不能错，要一摸一样一摸一样一摸一样！重要事情说三次！

seuset <- CreateSeuratObject(counts = dat1, # 单细胞数据 min.cells = 3, # 保留在3个以上细胞中测到的基因 min.features = 200) # 保留至少测到了200个基因的细胞# 这样一个Seurat对象就构建完成，我们看一下数据大小dim(seuset) # 14627 1264 总共有14627行1264列## [1] 14627 1264head([email protected]) # 得到的metadata文件,也就是细胞信息## orig.ident nCount_RNA nFeature_RNA## AAACCTGCATGTTGAC-1 SeuratProject 1862 920## AAACGGGGTGCGGTAA-1 SeuratProject 8522 3212## AAACGGGTCATGCATG-1 SeuratProject 3161 1660## AAAGATGAGGACCACA-1 SeuratProject 4859 1771## AAAGATGTCAACACAC-1 SeuratProject 5293 1895## AAAGCAAAGCCACGTC-1 SeuratProject 2996 1237

这样10x Genomics的Seurat对象也构建完成了。注意哈，不管是什么平台测的数据类型都可以构建scater或者Seurat对象用于后续分析，我这里只是展示两个不同的知识点，当然如果公司给了你h5文件，那你直接用Read10X_h5就可以了，其他步骤不变。那么今天第二份作业就是从GEO数据库下载一份10x Genomics的数据，然后构建Seurat对象。

好啦，今天的内容就到这里，给大家留了两份作业，检验大家的学习结果。今天开始大家就不能只在上洗手间的时候学习单细胞了/(ㄒoㄒ)/~~

后台回复“feng41”获得数据和代码（记得写作业），我们下期见啦！

单细胞分析专栏传送门

大年初一，单细胞分析从nature communication开始

8+单细胞文章从解读到套路总结全部送给你

30分多的单细胞分析文章是怎么做的？

不用R，点点点就能get的单细胞公共数据下载，你“学废”了吗