往期系列

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_annscellbarcodes <- 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_annscellbarcodes <- 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")

练习

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

Hemberg-lab单细胞转录组数据分析（七）-导入10X和SmartSeq2数据Tabula Muris相关推荐

  1. Hemberg-lab单细胞转录组数据分析（八） - Scater包输入导入和存储

    往期系列 亨贝格实验室单细胞转录组数据分析(一) 亨贝格实验室单细胞转录组数据分析(二) 亨贝格实验室单细胞转录组数据分析(三) 亨贝格实验室单细胞转录组数据分析(四) 亨贝格实验室单细胞转录组数据分 ...

  2. Hemberg-lab单细胞转录组数据分析（九）- Scater包单细胞过滤

    往期系列 Hemberg-lab单细胞转录组数据分析(一) Hemberg-lab单细胞转录组数据分析(二) Hemberg-lab单细胞转录组数据分析(三) Hemberg-lab单细胞转录组数据分 ...

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

    Hemberg-lab单细胞转录组数据分析(一) Hemberg-lab单细胞转录组数据分析(二) Hemberg-lab单细胞转录组数据分析(三) Hemberg-lab单细胞转录组数据分析(四) ...

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

    Hemberg-lab单细胞转录组数据分析(一) Hemberg-lab单细胞转录组数据分析(二) Hemberg-lab单细胞转录组数据分析(三) Hemberg-lab单细胞转录组数据分析(四) ...

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

    Hemberg-lab单细胞转录组数据分析(一) Hemberg-lab单细胞转录组数据分析(二) Hemberg-lab单细胞转录组数据分析(三) 收藏|北大生信平台"单细胞分析.染色质分 ...

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

    Hemberg-lab单细胞转录组数据分析(一) Hemberg-lab单细胞转录组数据分析(二) 收藏|北大生信平台"单细胞分析.染色质分析"视频和PPT分享 生信老司机以中心法 ...

  7. seurat提取表达矩阵_Hemberg-lab单细胞转录组数据分析

    单细胞RNA-seq简介 混合RNA-seq2000年末的重大技术突破,取代微阵列表达芯片被广泛使用 通过混合大量细胞获取足够RNA用于建库测序,来定量每个基因的平均表达水平 用于比较转录组,例如比较 ...

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

    在别人的电子书,你的电子书,都在bookdown一文中推荐过这一篇教程(https://hemberg-lab.github.io/scRNA.seq.course),从2016年一直更新到2018年 ...

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

    单细胞测序实验方法 开发scRNA-seq的新实验方法和操作手册是目前很火的研究领域,而且最近这些年已经发表了一些改进方法.上图可以看出检测的细胞数目以指数形式增加,以下是一份不完全的列表: STRT ...

最新文章

  1. 一块只要4美元，超廉价版树莓派诞生，还用上了自研芯片
  2. Silverlight RIA Services基础专题
  3. c语言数组用户注册登入管理系统_学生成绩管理系统案例
  4. OpenJ_Bailian - 1088 滑雪（记忆化搜索）
  5. windows 2003 上Lotus Notes 客户端无法运行的解决办法
  6. HOG特征，LBP特征，Haar特征（图像特征提取）
  7. python实现雪花飘落效果_python实现雪花飘落效果实例讲解
  8. web页面-JS/DOM/BOM/窗口滚动/修改内容/上传文件
  9. 汇率换算自然语言理解功能IOS DEMO
  10. 微信小程序之获取百度语音合成
  11. IE不能下载MSG文件的解决方案
  12. 我发表的论文，怎么跑到百度文库中了
  13. windows下tomcat8启动脚本代码剖析--catalina.bat
  14. vue下载和预览word
  15. “家贼”倒卖“征途”源代码 13万元卖给识货人
  16. 扒一扒坑人的“微信支付”SDK开发文档
  17. 设计模式：备忘录模式
  18. STM32串口通信（串口中断、FIFO机制）之安富莱代码学习笔记
  19. XP安装SQLSERVER企业版
  20. DJI OSDK开发笔记（N3飞控）（1）——开发工作流程

热门文章

  1. 【操作系统】Semaphore处理哲学家就餐问题
  2. 软件项目管理第三课—如何应对投标书的软件功能报价
  3. BGP 路由表即将突破 768k
  4. const、extern、static的使用不再神秘
  5. saltstack常用模块
  6. ifconfig route 手动设置网卡route路由 示例
  7. 老罗的情怀“被”与廉价划等号 锤子“被”选错了对手
  8. oracle 之 内存—鞭辟近里（一）
  9. XSL样式，分页方法
  10. Python基础学习笔记三