||
收藏|北大生信平台"单细胞分析、染色质分析"视频和PPT分享
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
读入逗号分隔的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)
为了构建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
技术细胞通量高但测序覆盖度低,所以其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")
写一个函数读取任意组织的任意类型的数据。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-20 04:59
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社