Tabula MurisTabula 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./files/10038307
unzip 10038307
wget https://ndownloader./files/10038310
mv 10038310 FACS_metadata.csv
wget https://ndownloader./files/10039267
mv 10039267 FACS_annotations.csv
终端下载10X数据: wget https://ndownloader./files/10038325
unzip 10038325
wget https://ndownloader./files/10038328
mv 10038328 droplet_metadata.csv
wget https://ndownloader./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")
|