往期系列

Hemberg-lab单细胞转录组数据分析（一）

Hemberg-lab单细胞转录组数据分析（二）

Hemberg-lab单细胞转录组数据分析（三）

Hemberg-lab单细胞转录组数据分析（四）

Hemberg-lab单细胞转录组数据分析（五）

Hemberg-lab单细胞转录组数据分析（六）

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

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

生物信息之程序学习

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

Tabula Muris

Tabula Muris是测序小鼠20个器官和组织的单细胞转录组图谱的国际合作项目 (Transcriptomic characterization of 20 organs and tissues from mouse at single cell resolution creates a Tabula Muris)。

简介

我们使用 Tabula Muris最开始释放的数据做为测试数据来完成完整的单细胞数据分析。The Tabula Muris是一个国际合作组织，目的是采用标准方法生成小鼠每个细胞的图谱。建库测序方法包括通量高覆盖率低的10X数据和通量低覆盖率高的FACS筛选+Smartseq2建库技术。

起始数据于2017年12月20日释放，包含20个组织/器官的100,000细胞的转录组图谱。

下载数据

与其它sc-RNASeq数据上传到GEO或ArrayExpress不同，Tabula Muris通过figshare平台释放数据。可以通过搜索文章的doi号10.6084/m9.figshare.5715040下载FACS/Smartseq2数据和10.6084/m9.figshare.5715025下载10X数据。数据可以直接点击链接下载，或使用下面的wget命令：

终端下载FACS数据:

wget https://ndownloader.figshare.com/files/10038307
unzip 10038307
wget https://ndownloader.figshare.com/files/10038310
mv 10038310 FACS_metadata.csv
wget https://ndownloader.figshare.com/files/10039267
mv 10039267 FACS_annotations.csv

终端下载10X数据:

wget https://ndownloader.figshare.com/files/10038325
unzip 10038325
wget https://ndownloader.figshare.com/files/10038328
mv 10038328 droplet_metadata.csv
wget https://ndownloader.figshare.com/files/10039264
mv 10039264 droplet_annotation.csv

如果数据是手动下载的，也需要像上面一样解压和重命名。

现在应该有两个文件夹: FACS和droplet，每个对应一个annotation和metadata文件。使用head命令查看前10行：

head -n 10 droplet_metadata.csv

使用wc -l查看文件的行数:

wc -l droplet_annotation.csv

练习：FACS和10X数据中各有多少细胞有注释信息？

答案: FACS : 54,838 cells; Droplet : 42,193 cells

读入数据 (Smartseq2)

读入逗号分隔的count matrix，存储为数据框:

dat = read.delim("FACS/Kidney-counts.csv", sep=",", header=TRUE)
dat[1:5,1:5]

数据库第一列是基因名字，把它移除作为列名字：

dim(dat)
rownames(dat) <- dat[,1]
dat <- dat[,-1]

这是Smartseq2数据集，可能含有spike-ins:

rownames(dat)[grep("^ERCC-", rownames(dat))]

从列名字中提取metadata信息:

cellIDs <- colnames(dat)
cell_info <- strsplit(cellIDs, "\\.")
Well <- lapply(cell_info, function(x){x[1]})
Well <- unlist(Well)
Plate <- unlist(lapply(cell_info, function(x){x[2]}))
Mouse <- unlist(lapply(cell_info, function(x){x[3]}))

检测每种metadata类型的数据分布：

summary(factor(Mouse))

查看有没有技术因子是cofounded，实验批次与供体小鼠批次一致:

table(Mouse, Plate)

最后读入计算预测的细胞类型注释，并与表达矩阵中的细胞注释做比较:

ann <- read.table("FACS_annotations.csv", sep=",", header=TRUE)
ann <- ann[match(cellIDs, ann[,1]),]
celltype <- ann[,3]
table(celltype)

构建scater对象

为了构建SingleCellExperiment对象，先把所有的细胞注释放到一个数据框中。因为实验批次(pcr plate)和供体小鼠完全重合，所以只保留一个信息。

suppressMessages(require("SingleCellExperiment"))
suppressMessages(require("scater"))
cell_anns <- data.frame(mouse = Mouse, well=Well, type=celltype)
rownames(cell_anns) <- colnames(dat)
sceset <- SingleCellExperiment(assays = list(counts = as.matrix(dat)), colData=cell_anns)

str(sceset)

如果数据集包含spike-ins，我们在SingleCellExperiment对象中定义一个变量记录它们。

isSpike(sceset, "ERCC") <- grepl("ERCC-", rownames(sceset))
str(sceset)

读入10X的数据

因为10X技术细胞通量高但测序覆盖度低，所以其count matrix是一个大的稀疏矩阵(矩阵中高达90%的数据的数值为0)。CellRanger默认的输出格式是.mtx文件用于存储这个稀疏矩阵，第一列是基因的坐标（0-based），第二列是细胞的坐标（0-based），第三列是大于0的表达值 (长表格形式)。 打开.mtx文件会看到两行标题行后面是包含总行数 (基因数)、列数 (样本数)和稀疏矩阵总行数 (生信宝典注：所有细胞中表达不为0的基因的总和)的一行数据。

%%MatrixMarket matrix coordinate integer general
%
23433 610 1392643 5 1 1
28 1 1
40 1 2

鉴于.mtx文件中只存储了基因和样品名字的坐标，而实际的基因和样品的名字必须单独存储到文件genes.tsv和barcodes.tsv。

下面使用Matrix包读入稀疏矩阵：

suppressMessages(require("Matrix"))
cellbarcodes <- read.table("droplet/Kidney-10X_P4_5/barcodes.tsv")
genenames <- read.table("droplet/Kidney-10X_P4_5/genes.tsv")
molecules <- Matrix::readMM("droplet/Kidney-10X_P4_5/matrix.mtx")

下一步增加合适的行或列的名字。首先查看read的cellbarcode信息会发现这个文件只有barcode序列。考虑到10X数据每一批的cellbarcode是有重叠的，所以在合并数据前，需要把批次信息与barcode信息合并一起。

head(cellbarcodes)

rownames(molecules) <- genenames[,1]
colnames(molecules) <- paste("10X_P4_5", cellbarcodes[,1], sep="_")

读入计算注释的细胞类型信息:

meta <- read.delim("droplet_metadata.csv", sep=",", header=TRUE)
head(meta)

我们需要用10X_P4_5获得这批数据对应的metadata信息。这时需要注意metadata表格中mouse ID与前面plate-based (FACS SmartSeq2)数据集的mouse ID不同，这里用-而非_作为分隔符，并且性别在中间。通过查阅文献中的描述得知droplet (10X)和plate-based (FACS SmartSeq2)的技术用了同样的8只老鼠。所以对数据做下修正，使得10X与FACS的数据一致。

meta[meta$channel == "10X_P4_5",]
mouseID <- "3_8_M"

注意：有些组织的10X数据可能来源于多个小鼠的样品，如mouse id = 3-M-5/6。也需要格式化这些信息，但可能这些与FACS数据的mouse id会不一致，进而影响下游分析。如果小鼠不是纯系，可能需要通过exonic-SNP把细胞和对应的小鼠联系起来 (本课程不会涉及)。

ann <- read.delim("droplet_annotation.csv", sep=",", header=TRUE)
head(ann)

注释中的cellID和cellbarcodes也存在细微差别，少了最后的-1，在匹配前需要做下校正。(生信宝典注：这种数据不一致是经常要处理的问题，每一步检查结果。如果与预期不符，考虑有没有未考虑到的数据不一致的地方。)

ann[,1] <- paste(ann[,1], "-1", sep="")
ann_subset <- ann[match(colnames(molecules), ann[,1]),]
celltype <- ann_subset[,3]

构建cell-metadata数据框:

cell_anns <- data.frame(mouse = rep(mouseID, times=ncol(molecules)), type=celltype)
rownames(cell_anns) <- colnames(molecules)
head(cell_anns)

练习 用组织的其它批次重复上面的处理。

答案

molecules1 <- molecules
cell_anns1 <- cell_anns

cellbarcodes <- read.table("droplet/Kidney-10X_P4_6/barcodes.tsv")
genenames <- read.table("droplet/Kidney-10X_P4_6/genes.tsv")
molecules <- Matrix::readMM("droplet/Kidney-10X_P4_6/matrix.mtx")
rownames(molecules) <- genenames[,1]
colnames(molecules) <- paste("10X_P4_6", cellbarcodes[,1], sep="_")
mouseID <- "3_9_M"
ann_subset <- ann[match(colnames(molecules), ann[,1]),]
celltype <- ann_subset[,3]
cell_anns <- data.frame(mouse = rep(mouseID, times=ncol(molecules)), type=celltype)
rownames(cell_anns) <- colnames(molecules)

molecules2 <- molecules
cell_anns2 <- cell_anns

cellbarcodes <- read.table("droplet/Kidney-10X_P7_5/barcodes.tsv")
genenames <- read.table("droplet/Kidney-10X_P7_5/genes.tsv")
molecules <- Matrix::readMM("droplet/Kidney-10X_P7_5/matrix.mtx")
rownames(molecules) <- genenames[,1]
colnames(molecules) <- paste("10X_P7_5", cellbarcodes[,1], sep="_")
mouseID <- "3_57_F"
ann_subset <- ann[match(colnames(molecules), ann[,1]),]
celltype <- ann_subset[,3]
cell_anns <- data.frame(mouse = rep(mouseID, times=ncol(molecules)), type=celltype)
rownames(cell_anns) <- colnames(molecules)

molecules3 <- molecules
cell_anns3 <- cell_anns

创建scater对象

现在读入了多个批次的10X数据，把它们组合成一个SingleCellExperiment object对象。首先检查不同批次数据的基因名字是否一致:

identical(rownames(molecules1), rownames(molecules2))
identical(rownames(molecules1), rownames(molecules3))

确认没有重复的细胞ID:

sum(colnames(molecules1) %in% colnames(molecules2))
sum(colnames(molecules1) %in% colnames(molecules3))
sum(colnames(molecules2) %in% colnames(molecules3))

检查无误，把它们组合起来:

# 获得大的表达矩阵
all_molecules <- cbind(molecules1, molecules2, molecules3)
# 获得大的数据矩阵
all_cell_anns <- as.data.frame(rbind(cell_anns1, cell_anns2, cell_anns3))
# 增加批次信息
all_cell_anns$batch <- rep(c("10X_P4_5", "10X_P4_6","10X_P7_5"), times = c(nrow(cell_anns1), nrow(cell_anns2), nrow(cell_anns3)))

练习: 全部数据集有多少细胞？

答案：

dim(all_molecules)[2]

现在创建SingleCellExperiment对象。SingleCellExperiment对象的优势是可以正常矩阵、稀疏矩阵格式存储数据，还可以以HDF5格式在磁盘存储和访问大的非稀疏矩阵而不用全部加载到内存中。

require("SingleCellExperiment")
require("scater")
all_molecules <- as.matrix(all_molecules)
sceset <- SingleCellExperiment(assays = list(counts = all_molecules), colData=all_cell_anns)

这是10X的数据，不包含spike-ins直接存储数据:

saveRDS(sceset, "kidney_droplet.rds")

练习

写一个函数读取任意组织的任意类型的数据。