万水千山粽是情，点开看看行不行 / 四六文摘

为了让您亲身体验从头到尾分析单细胞RNASeq数据集，我们将使用Tabula Muris初始版本的数据作为示例。Tabula Muris是一项国际合作项目，旨在使用标准化方法对鼠中的每种细胞类型进行分析。它们将高通量但低覆盖率的10X数据与低通量但高覆盖率的FACS分选细胞+ Smartseq2相结合。

最初发布的数据（2017年12月20日）包含20个不同组织/器官的近100,000个细胞。您可以使用这些组织中的一个作为例子来处理这个过程。您可以选择一个组织去作为一个例子来进行本节课的工作。

6.2

下载数据

与大多数单细胞RNASeq数据不同，Tabula Muris通过figshare平台发布数据，而不是将其上传到GEO或ArrayExpress。：您可以通过使用在他们的论文中的DOI找到数据，用于FACS/Smartseq2的10.6084 / m9.figshare.5715040 和10.6084 / m9.figshare.5715025的10X数据。可以通过使用doi链接或使用以下命令行命令来手动下载数据：

基于终端的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”以及每个文件夹会有一个注释和元数据文件。要检查这些文件，您可以使用head来查看文本文件的前几行（按“q”退出）：

head -n 10 droplet_metadata.csv## channel,mouse.id,tissue,subtissue,mouse.sex
## 10X_P4_0,3-M-8,Tongue,NA,M
## 10X_P4_1,3-M-9,Tongue,NA,M
## 10X_P4_2,3-M-8/9,Liver,hepatocytes,M
## 10X_P4_3,3-M-8,Bladder,NA,M
## 10X_P4_4,3-M-9,Bladder,NA,M
## 10X_P4_5,3-M-8,Kidney,NA,M
## 10X_P4_6,3-M-9,Kidney,NA,M
## 10X_P4_7,3-M-8,Spleen,NA,M
## 10X_P7_0,3-F-56,Liver,NA,F

您还可以使用以下方法检查每个文件中的行数：

wc -l droplet_annotation.csv
## 54838 droplet_annotation.csv

练习我们有多少个细胞来自FACS吗？10X呢？

答 FACS：54,838；10X：42,193

6.3

读取数据（Smartseq2）

我们现在可以从逗号分隔文件中读取相关的计数矩阵。然后检查作为结果的数据框：

dat = read.delim("FACS/Kidney-counts.csv", sep=",", header=TRUE)
dat[1:5,1:5]## X A14.MAA000545.3_8_M.1.1 E1.MAA000545.3_8_M.1.1
## 1 0610005C13Rik 0 0
## 2 0610007C21Rik 1 0
## 3 0610007L01Rik 0 0
## 4 0610007N19Rik 0 0
## 5 0610007P08Rik 0 0
## M4.MAA000545.3_8_M.1.1 O21.MAA000545.3_8_M.1.1
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0

我们可以看到数据框中的第一列是基因名称，所以首先我们将它们移动到rownames，所以我们就有了一个数字矩阵：

dim(dat)## [1] 23433 866rownames(dat) <- dat[,1]
dat <- dat[,-1]

由于这是一个Smartseq2数据集，它可能包含spike-ins，所以我们可以检查：

rownames(dat)[grep("^ERCC-", rownames(dat))]## [1] "ERCC-00002" "ERCC-00003" "ERCC-00004" "ERCC-00009" "ERCC-00012"
## [6] "ERCC-00013" "ERCC-00014" "ERCC-00016" "ERCC-00017" "ERCC-00019"
## [11] "ERCC-00022" "ERCC-00024" "ERCC-00025" "ERCC-00028" "ERCC-00031"
## [16] "ERCC-00033" "ERCC-00034" "ERCC-00035" "ERCC-00039" "ERCC-00040"
## [21] "ERCC-00041" "ERCC-00042" "ERCC-00043" "ERCC-00044" "ERCC-00046"
## [26] "ERCC-00048" "ERCC-00051" "ERCC-00053" "ERCC-00054" "ERCC-00057"
## [31] "ERCC-00058" "ERCC-00059" "ERCC-00060" "ERCC-00061" "ERCC-00062"
## [36] "ERCC-00067" "ERCC-00069" "ERCC-00071" "ERCC-00073" "ERCC-00074"
## [41] "ERCC-00075" "ERCC-00076" "ERCC-00077" "ERCC-00078" "ERCC-00079"
## [46] "ERCC-00081" "ERCC-00083" "ERCC-00084" "ERCC-00085" "ERCC-00086"
## [51] "ERCC-00092" "ERCC-00095" "ERCC-00096" "ERCC-00097" "ERCC-00098"
## [56] "ERCC-00099" "ERCC-00104" "ERCC-00108" "ERCC-00109" "ERCC-00111"
## [61] "ERCC-00112" "ERCC-00113" "ERCC-00116" "ERCC-00117" "ERCC-00120"
## [66] "ERCC-00123" "ERCC-00126" "ERCC-00130" "ERCC-00131" "ERCC-00134"
## [71] "ERCC-00136" "ERCC-00137" "ERCC-00138" "ERCC-00142" "ERCC-00143"
## [76] "ERCC-00144" "ERCC-00145" "ERCC-00147" "ERCC-00148" "ERCC-00150"
## [81] "ERCC-00154" "ERCC-00156" "ERCC-00157" "ERCC-00158" "ERCC-00160"
## [86] "ERCC-00162" "ERCC-00163" "ERCC-00164" "ERCC-00165" "ERCC-00168"
## [91] "ERCC-00170" "ERCC-00171"

现在我们可以从列名中提取这些数据的大部分元数据：

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]}))

我们可以检查每个元数据分类的分布：

summary(factor(Mouse))## 3_10_M 3_11_M 3_38_F 3_39_F 3_8_M 3_9_M
## 104 196 237 169 82 77

我们还可以检查是否有任何技术因素混淆：

table(Mouse, Plate)## Plate
## Mouse B001717 B002775 MAA000545 MAA000752 MAA000801 MAA000922
## 3_10_M 0 0 0 104 0 0
## 3_11_M 0 0 0 0 196 0
## 3_38_F 237 0 0 0 0 0
## 3_39_F 0 169 0 0 0 0
## 3_8_M 0 0 82 0 0 0
## 3_9_M 0 0 0 0 0 77

最后，我们将读取用于推断的细胞类型注释，并将它们与表达矩阵中的细胞进行匹配：

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

6.4

构建一个scater对象

要创建一个SingleCellExperiment对象，我们必须将所有细胞注释放在一个数据框中，因为实验批处理（pcr板）与供体小鼠完全混淆，所以我们只保留其中一个。

require("SingleCellExperiment")
## Loading required package: SingleCellExperiment
## Loading required package: SummarizedExperiment
## Loading required package: methods
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, cbind, colMeans,
## colnames, colSums, do.call, duplicated, eval, evalq, Filter,
## Find, get, grep, grepl, intersect, is.unsorted, lapply,
## lengths, Map, mapply, match, mget, order, paste, pmax,
## pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce,
## rowMeans, rownames, rowSums, sapply, setdiff, sort, table,
## tapply, union, unique, unsplit, which, which.max, which.min## Loading required package: S4Vectors##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:base':
##
## expand.grid## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Loading required package: DelayedArray
## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians
##
## Attaching package: 'DelayedArray'
## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
## The following object is masked from 'package:base':
##
## apply
require("scater")
## Loading required package: scater
## Loading required package: ggplot2
##
## Attaching package: 'scater'
## The following object is masked from 'package:S4Vectors':
##
## rename
## The following object is masked from 'package:stats':
##
## filter
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)

最后，如果数据集包含spike-ins，我们会在SingleCellExperiment对象中隐藏一个隐藏变量来跟踪它们：

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

6.5

读取数据（10X）

由于10X数据的大尺寸和稀疏性（表达矩阵可能高达90％也可能是0），它通常存储为稀疏矩阵。CellRanger的默认输出格式是.mtx文件，CellRanger的默认输出格式是.mtx文件，该文件会将该稀疏矩阵存储其行坐标列，列坐标列和表达式值>0的列。注意，如果查看.mtx文件，您将看到两个标题行后面跟着一行，用来详细说明完整矩阵的行数，列数和计数（counts）。由于只有坐标存储在.mtx文件中，因此每行和列的名称必须分别存储在“genes.tsv”和“barcodes.tsv”文件中。

我们将使用“Matrix”包在R中以稀疏矩阵的格式存储矩阵。

require("Matrix")
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following object is masked from 'package:S4Vectors':
##
## expand
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的细胞barcode，您将看到它们只是与每个细胞关联的barcode序列。这是一个问题，因为每批10X数据使用相同的barcode池，因此如果我们需要组合来自多个10X批次的数据，则细胞的barcode将不相同。因此，我们把每一批次的ID附加到每个细胞的barcode：

head(cellbarcodes)
## V1
## 1 AAACCTGAGATGCCAG-1
## 2 AAACCTGAGTGTCCAT-1
## 3 AAACCTGCAAGGCTCC-1
## 4 AAACCTGTCCTTGCCA-1
## 5 AAACGGGAGCTGAACG-1
## 6 AAACGGGCAGGACCCT-1
rownames(molecules) <- genenames[,1]
colnames(molecules) <- paste("10X_P4_5", cellbarcodes[,1], sep="_")

现在让我们获取这些数据的元数据和计算注释：

meta <- read.delim("droplet_metadata.csv", sep=",", header=TRUE)
head(meta)
## channel mouse.id tissue subtissue mouse.sex
## 1 10X_P4_0 3-M-8 Tongue <NA> M
## 2 10X_P4_1 3-M-9 Tongue <NA> M
## 3 10X_P4_2 3-M-8/9 Liver hepatocytes M
## 4 10X_P4_3 3-M-8 Bladder <NA> M
## 5 10X_P4_4 3-M-9 Bladder <NA> M
## 6 10X_P4_5 3-M-8 Kidney <NA> M

在这里我们可以看到我们需要使用“10X_P4_5”来查找该批次的元数据，还要注意小鼠ID的格式在此元数据表中是不同的，带有连字符而不是下划线，并且性别位于ID的中间。通过检查随附论文的方法部分，我们知道相同的8只小鼠都用了droplet和plate的技术。因此，我们需要修复小鼠ID，以与FACS实验中使用的ID一致。

meta[meta$channel == "10X_P4_5",]
## channel mouse.id tissue subtissue mouse.sex
## 6 10X_P4_5 3-M-8 Kidney <NA> M
mouseID <- "3_8_M"

注意：您有可能从混合样本中获得10X数据：例如，小鼠id = 3-M-5/6。您仍应将这些重新格式化为一致，但它们不会匹配可能影响下游分析的FACS数据中的小鼠ID。如果小鼠不是来自近交系，那么可以使用exonic-SNP将单个细胞分配给特定小鼠，但这超出了本课程的范围。

ann <- read.delim("droplet_annotation.csv", sep=",", header=TRUE)
head(ann)
## cell tissue cell_ontology_class
## 1 10X_P4_3_AAAGTAGAGATGCCAG Bladder mesenchymal cell
## 2 10X_P4_3_AACCGCGTCCAACCAA Bladder mesenchymal cell
## 3 10X_P4_3_AACTCCCGTCGGGTCT Bladder mesenchymal cell
## 4 10X_P4_3_AACTCTTAGTTGCAGG Bladder bladder cell
## 5 10X_P4_3_AACTCTTTCATAACCG Bladder mesenchymal cell
## 6 10X_P4_3_AAGACCTAGATCCGAG Bladder endothelial cell
## cell_ontology_term_iri cell_ontology_id
## 1 http://purl.obolibrary.org/obo/CL_0008019 CL:0008019
## 2 http://purl.obolibrary.org/obo/CL_0008019 CL:0008019
## 3 http://purl.obolibrary.org/obo/CL_0008019 CL:0008019
## 4 http://purl.obolibrary.org/obo/CL_1001319 CL:1001319
## 5 http://purl.obolibrary.org/obo/CL_0008019 CL:0008019
## 6 http://purl.obolibrary.org/obo/CL_0000115 CL:0000115

再一次，您会发现注释中的cellID与我们在匹配之前必须更正的细胞barcode之间存在轻微的格式差异。

ann[,1] <- paste(ann[,1], "-1", sep="")
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);

练习对你选择的组织的其他10X批次重复上述操作。

6.6

构建scater对象

现在我们已经在多个批次中读取了10X数据，我们需要将它们组合成一个SingleCellExperiment对象。首先，我们将检查所有批次中基因名称是否相同且顺序是否相同：

identical(rownames(molecules1), rownames(molecules2))
## [1] TRUE
identical(rownames(molecules1), rownames(molecules3))
## [1] TRUE

现在我们将检查没有任何重复的cellID：

sum(colnames(molecules1) %in% colnames(molecules2))
## [1] 0
sum(colnames(molecules1) %in% colnames(molecules3))
## [1] 0
sum(colnames(molecules2) %in% colnames(molecules3))
## [1] 0

一切都很好，所以我们可以继续并将它们结合起来：

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

练习整个数据集中有多少个细胞？

现在构建SingleCellExperiment对象。SingleCellExperiment类的一个优点是它能够以正常矩阵或稀疏矩阵格式存储数据，以及HDF5格式，它允许以有效的方式在磁盘上存储和访问大型非稀疏矩阵，而不是加载所有到RAM上。

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

由于这是10X数据，因此不会包含spike-ins，因此我们只保存数据：

saveRDS(sceset, "kidney_droplet.rds")

6.7

高级练习

编写一个R函数/脚本，它将为任何组织的每种数据类型都可以完全自动化此过程。

万水千山粽是情，点开看看行不行

相关推荐