往期系列
收藏|北大生信平台"单细胞分析、染色质分析"视频和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")
练习
写一个函数读取任意组织的任意类型的数据。
