此文章是通过学习瑞典国家生物信息学基础设施(NBIS) 所开放的单细胞分析教程加上网上所查找的资料，自身的理解所形成的，可能会有不足之处。该部分是对下机处理完成后的数据进行Seurat分析的质控。
参考来源：https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/seurat/seurat_01_qc.html
感兴趣的话可以阅读原文。
安装相关包并下载参考数据
参考数据下载
原文中提供了演示数据的下载方式，是通过R下载的，但我觉得下载方式无所谓，你熟悉什么方式就用什么方式下载就好，这里贴上网址，一共是3个covid-19患者和3个健康对照的一共6个PBMC的10x数据集：
Normal_PBMC_13:"https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/new_dataset/labs/data/covid_data_GSE149689/sub/Normal_PBMC_13.h5"
Normal_PBMC_14:"https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/new_dataset/labs/data/covid_data_GSE149689/sub/Normal_PBMC_14.h5"
Normal_PBMC_5:"https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/new_dataset/labs/data/covid_data_GSE149689/sub/Normal_PBMC_5.h5"
nCoV_PBMC_15:"https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/new_dataset/labs/data/covid_data_GSE149689/sub/nCoV_PBMC_15.h5"
nCoV_PBMC_17:"https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/new_dataset/labs/data/covid_data_GSE149689/sub/nCoV_PBMC_17.h5"
nCoV_PBMC_1:"https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/new_dataset/labs/data/covid_data_GSE149689/sub/nCoV_PBMC_1.h5"
原文中的下载方式,下载之前最好先设置工作目录：
setwd("/Volumes/Flower0501/研究学习目录/单细胞转录组/分析流程/Seurat下游分析/")
webpath <- "https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/new_dataset/labs/data/covid_data_GSE149689/sub/"
dir.create("./data/raw", recursive = T)
file_list <- c("Normal_PBMC_13.h5", "Normal_PBMC_14.h5", "Normal_PBMC_5.h5", "nCoV_PBMC_15.h5",
"nCoV_PBMC_17.h5", "nCoV_PBMC_1.h5")
for (i in file_list) {
download.file(url = paste0(webpath, i), destfile = paste0("./data/raw/", i))
}
相关包下载
# install.packages('Seurat')
# install.packages('Matrix')
# remotes::install_github("chris-mcginnis-ucsf/DoubletFinder", upgrade = F)
suppressMessages(require(Seurat))
suppressMessages(require(Matrix))
suppressMessages(require(DoubletFinder))
读取数据
先使用Read10X_h5读取文件，在使用CreateSeuratObject创建Seurat对象,使用循环，一步到位，当然也可以一个一个读取，然后在合并。
挨个读取
# 使用Read10X_h5函数读取h5文件
cov.15 <- Seurat::Read10X_h5(filename = "data/raw/nCoV_PBMC_15.h5", use.names = T)
cov.1 <- Seurat::Read10X_h5(filename = "data/raw/nCoV_PBMC_1.h5", use.names = T)
cov.17 <- Seurat::Read10X_h5(filename = "data/raw/nCoV_PBMC_17.h5", use.names = T)
ctrl.5 <- Seurat::Read10X_h5(filename = "data/raw/Normal_PBMC_5.h5", use.names = T)
ctrl.13 <- Seurat::Read10X_h5(filename = "data/raw/Normal_PBMC_13.h5", use.names = T)
ctrl.14 <- Seurat::Read10X_h5(filename = "data/raw/Normal_PBMC_14.h5", use.names = T)
# 使用CreateSeuratObject函数构建seurat对象
sdata.cov15 <- CreateSeuratObject(cov.15, project = "covid_15")
sdata.cov1 <- CreateSeuratObject(cov.1, project = "covid_1")
sdata.cov17 <- CreateSeuratObject(cov.17, project = "covid_17")
sdata.ctrl5 <- CreateSeuratObject(ctrl.5, project = "ctrl_5")
sdata.ctrl13 <- CreateSeuratObject(ctrl.13, project = "ctrl_13")
sdata.ctrl14 <- CreateSeuratObject(ctrl.14, project = "ctrl_14")
# 添加样本分组信息
sdata.cov1$type = "Covid"
sdata.cov15$type = "Covid"
sdata.cov17$type = "Covid"
sdata.ctrl5$type = "Ctrl"
sdata.ctrl13$type = "Ctrl"
sdata.ctrl14$type = "Ctrl"
循环读取
data_list <- list.files(".")
data <- list()
cell_ids <- c()
for(i in data_list){
if (grepl("CoV",i)){
Nam <- paste("covid",str_match(i,".*_(\\d+).*")[2],sep = "_")
type <- "Covid"
} else {
Nam <- paste("ctrl",str_match(i,".*_(\\d+).*")[2],sep = "_")
type <- "Ctrl"
}
data[[i]] <- CreateSeuratObject(Seurat::Read10X_h5(filename = i, use.names = T),project = Nam)
data[[i]]$type <- type
cell_ids <- c(cell_ids,Nam)
}
PS:使用循环读取，需要在读取前对文件名进行重命名。让文件名变得有规律。
合并数据
# 合并所有样本，然后删除没有必要的数据，清楚内存
# 挨个读取
alldata <- merge(x = sdata.cov15,
y = c(sdata.cov1, sdata.cov17, sdata.ctrl5,
sdata.ctrl13, sdata.ctrl14),
add.cell.ids = c("covid_15", "covid_1", "covid_17",
"ctrl_5", "ctrl_13", "ctrl_14"))
rm(cov.15, cov.1, cov.17,
ctrl.5, ctrl.13, ctrl.14,
sdata.cov15 , sdata.cov1, sdata.cov17,
sdata.ctrl5, sdata.ctrl13, sdata.ctrl14)
# 循环读取
alldata <- merge(x = data[[1]],y=c(data[[2]],data[[3]],data[[4]],data[[5]],data[[6]]),add.cell.ids=cell_ids)
rm(data)
PS：合并数据也不一定要用Seurat进行合并，可以看到合并两三个样本还好，多个样本就会导致操作繁琐，推荐在读取进入Seurat之前，就使用
Cellranger aggr把样本进行合并。
质控指标的计算
质控的目的是为了去掉低质量的数据，并且检查细胞的分布是否符合预期，基本检查的指标有基因总数，细胞总数，线粒体相关基因含量，核糖体相关基因含量，血红蛋白相关基因含量等，还可以有其他的检查指标，可以根据自己课题需要进行设置。
在上述的基本检查指标中，基因总数，细胞总数在创建Seurat对象的时候就已经计算好了，剩下的线粒体，核糖体和血红蛋白数量可以通过Seurat自带函数PercentageFeatureSet进行计算，也可以手动计算。
| 检测指标 | 检测目的 | 参考QC值 | 基因前缀 |
|---|---|---|---|
| 线粒体相关基因含量 | 异常：过高的线粒体含量可能是因为细胞质RNA穿过细胞膜导致，或者在取样的过程中死亡细胞过多 正常：心肌等耗能较多的地方，线粒体基因含量也可能会比较高 |
人：MT 小鼠：mt |
|
| 核糖体相关基因含量 | 人：RP[SL] 小鼠：rp[sl] |
||
| 血红蛋白相关基因含量 | 检测是否存在血红蛋白污染 | 人：HB[ ^§] 小鼠：hb] ^§] |
PS：一般情况下，小鼠的相关基因的名称是小写，人相关基因的名称是大写的
- 线粒体相关基因含量计算
# 线粒体相关基因含量计算，计算线粒体相关基因百分比
# 方法1:使用Seurat自带函数计算
# 使用PercentageFeatureSet函数计算线粒体基因含量
alldata <- PercentageFeatureSet(alldata, "^MT-", col.name = "percent_mito")
# 方法2:手动计算线粒体含量百分比
total_counts_per_cell <- colSums(alldata@assays$RNA@counts) ## 首先计算每个细胞中所有基因的总量
# 获取线粒体相关基因
mito_genes <- rownames(alldata)[grep("^MT-", rownames(alldata))]
# head(mito_genes)
# [1] "MT-ND1" "MT-ND2" "MT-CO1" "MT-CO2" "MT-ATP8" "MT-ATP6"
# 计算线粒体基因含量
alldata$percent_mito <- colSums(alldata@assays$RNA@counts[mito_genes, ])/total_counts_per_cell*100
- 核糖体相关基因的计算
# 核糖体相关基因计算，也和上面线粒体计算一样，可以直接使用Seurat中的自带函数进行计算，也可以自行计算，结果是一致的
# 自带函数PercentageFeatureSet
alldata <- PercentageFeatureSet(alldata, "^RP[SL]", col.name = "percent_ribo")
# 自行计算
total_counts_per_cell <- colSums(alldata@assays$RNA@counts) ## 首先计算每个细胞中所有基因的总量
ribo_genes <- rownames(alldata)[grep("^RP[SL]", rownames(alldata))]
head(ribo_genes, 10)
# [1] "RPL22" "RPL11" "RPS6KA1" "RPS8" "RPL5" "RPS27" "RPS6KC1" "RPS7"
# [9] "RPS27A" "RPL31"
alldata$percent_ribo <- colSums(alldata@assays$RNA@counts[ribo_genes, ])/total_counts_per_cell*100
- 血红蛋白相关基因的计算
# 血红细胞相关基因计算
# 自带函数PercentageFeatureSet
alldata <- PercentageFeatureSet(alldata, "^HB[^(P)]", col.name = "percent_hb")
# 自行计算
total_counts_per_cell <- colSums(alldata@assays$RNA@counts) ## 首先计算每个细胞中所有基因的总量
hb_genes <- rownames(alldata)[grep("^HB[^(P)]", rownames(alldata))]
head(hb_genes)
# [1] "HBEGF" "HBS1L" "HBB" "HBD" "HBG1" "HBG2"
alldata$percent_hb <- colSums(alldata@assays$RNA@counts[hb_genes, ])/total_counts_per_cell*100
- 质控结果可视化
# 质控结果的可视化
feats <- c("nFeature_RNA", "nCount_RNA", "percent_mito", "percent_ribo", "percent_hb")
# 一般用小提琴图对结果进行质控
VlnPlot(alldata, group.by = "orig.ident", features = feats, pt.size = 0.1, ncol = 3) +
NoLegend()
同时，也需要查看基因数量随细胞数量的变化关系，与线粒体基因数量岁细胞数量的变化关系等，这样可以看出是否存在实验上的错误。
# 基因数量随细胞数量的变化关系
FeatureScatter(alldata, "nCount_RNA", "nFeature_RNA", group.by = "orig.ident", pt.size = 0.5)
# 线粒体基因数量随细胞数量的变化关系
FeatureScatter(alldata, "nCount_RNA", "percent_mito", group.by = "orig.ident", pt.size = 0.5)
质控过滤
其实过滤的指标并不固定，需要根据自身课题的需要进行分析。
基于检测的过滤指标
教程中采用的过滤条件是，至少检出200个基因的细胞，与至少要在三个细胞中存在的基因。
selected_c <- WhichCells(alldata, expression = nFeature_RNA > 200)
selected_f <- rownames(alldata)[Matrix::rowSums(alldata) > 3]
data.filt <- subset(alldata, features = selected_f, cells = selected_c)
dim(data.filt)
## [1] 18147 7973
初筛之后，还可以看看过滤之后，哪些基因的占比比较高，从而辅助判断是否过滤成功。注意，这一步的时间随着你数据体量的增大而增大。
par(mar = c(4, 8, 2, 1))
C <- data.filt@assays$RNA@counts
C <- Matrix::t(Matrix::t(C)/Matrix::colSums(C)) * 100
most_expressed <- order(apply(C, 1, median), decreasing = T)[20:1]
boxplot(as.matrix(t(C[most_expressed, ])), cex = 0.1, las = 1, xlab = "% total count per cell",
col = (scales::hue_pal())(20)[20:1], horizontal = TRUE)
可以看出MALAT1基因的数量占了基因总数的30%，这是一个非常不正常的值，可能是由于实验或测序过程中的技术错误导致的，所以在后续的分析中需要讲这个基因从基因列表中删除，在质控中，这一步很重要，是对质控结果的一次检验，也就是说，质控的结果中，不应该存在一个基因占基因总数较大的百分比。
基于核糖体与线粒体的过滤指标
可以从上图看出，在过滤结果中占比较多的基因依旧是与核糖体和线粒体相关的基因，对这些基因进行下一步的过滤的时候，有三种选择可以考虑：
- 第一种：如果在第一步基本过滤后，还存在较多的细胞，那么可以把线粒体基因过滤的阈值降到最低，直接过滤到含有线粒体基因的细胞。
- 第二种：删除所有与线粒体相关的基因，
- 第三种：对每个细胞的线粒体基因含量占比进行计算，删除那些线粒体相关基因含量占比较高的细胞。
一般采用较低的阈值先过滤异常的细胞，然后在除去线粒体与核糖体的基因，这篇教程里就是这么做的。
selected_mito <- WhichCells(data.filt, expression = percent_mito < 0.2)
selected_ribo <- WhichCells(data.filt, expression = percent_ribo > 0.05)
# and subset the object to only keep those cells
data.filt <- subset(data.filt, cells = selected_mito)
data.filt <- subset(data.filt, cells = selected_ribo)
dim(data.filt)
## [1] 18147 5762
table(data.filt$orig.ident)
## covid_1 covid_15 covid_17 ctrl_13 ctrl_14 ctrl_5
## 878 585 1042 1154 1063 1040
# 在对其进行质控绘图
feats <- c("nFeature_RNA", "nCount_RNA", "percent_mito", "percent_ribo", "percent_hb")
VlnPlot(data.filt, group.by = "orig.ident", features = feats, pt.size = 0.1, ncol = 3) +
NoLegend()
可以看到具有较高占比线粒体相关基因的细胞已经被删除了。
过滤基因
如前所述，除了线粒体与血红蛋白相关基因外，MALAT1也是一个异常基因，在这里也要剔除。
dim(data.filt)
## [1] 18147 5762
# Filter MALAT1
data.filt <- data.filt[!grepl("MALAT1", rownames(data.filt)), ]
# Filter Mitocondrial
data.filt <- data.filt[!grepl("^MT-", rownames(data.filt)), ]
# Filter Ribossomal gene (optional if that is a problem on your data) data.filt
# 是否剔除核糖体相关基因需要看你的课题要求来
# <- data.filt[ ! grepl('^RP[SL]', rownames(data.filt)), ]
# Filter Hemoglobin gene (optional if that is a problem on your data)
data.filt <- data.filt[!grepl("^HB[^(P)]", rownames(data.filt)), ]
dim(data.filt)
## [1] 18121 5762
样本性别控制
这也是要考虑的一个因素吧，因为在处理人类和动物样本时，样本的性别会对实验的恶结果产生较大的影响，所以应该在实验设计中尽量避免包含性别的差异，只保留单一性别，但是现实是残酷的。通过对于Y染色体（男性）和XIST（X非活性特异性转录本）（主要是女性）表达的分析，可以确定样本的性别，并且与metadata的信息进行相互验证。
一般要获取某个染色体上的所有基因时，通常会去查询对应的gtf文件，因为其具有相同注释版本的基因名，如果没有，就使用biomart来获取相关信息。在这篇教程中，他们给了相应的注释文件，下载地址：https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/master/labs/misc/genes.table.csv
genes.file = "genes.table.csv"
if (!file.exists(genes.file)) {
suppressMessages(require(biomaRt))
# initialize connection to mart, may take some time if the sites are
# unresponsive.
mart <- useMart("ENSEMBL_MART_ENSEMBL", dataset = "hsapiens_gene_ensembl")
# fetch chromosome info plus some other annotations
genes.table <- try(biomaRt::getBM(attributes = c("ensembl_gene_id", "external_gene_name",
"description", "gene_biotype", "chromosome_name", "start_position"), mart = mart,
useCache = F))
if (!dir.exists("data/results")) {
dir.create("data/results")
}
if (is.data.frame(genes.table)) {
write.csv(genes.table, file = genes.file)
}
if (!file.exists(genes.file)) {
download.file("https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/master/labs/misc/genes.table.csv",
destfile = "data/results/genes.table.csv")
genes.table = read.csv(genes.file)
}
} else {
genes.table = read.csv(genes.file)
}
genes.table <- genes.table[genes.table$external_gene_name %in% rownames(data.filt),]
# 提取相应的基因进行绘图
chrY.gene = genes.table$external_gene_name[genes.table$chromosome_name == "Y"]
data.filt$pct_chrY = colSums(data.filt@assays$RNA@counts[chrY.gene, ])/colSums(data.filt@assays$RNA@counts)
- XIST与Y的相关性散点图
FeatureScatter(data.filt, feature1 = "XIST", feature2 = "pct_chrY")
可以看出，点基本位于图像的左下角，证明基本上没有在样本中检测出性别染色体上的基因。
- XIST与Y染色体相关基因的样本分布图
VlnPlot(data.filt, features = c("XIST", "pct_chrY"))
可以看出有4个女性，两个男性，而且两个男性都是新冠患者，其实，我感觉在做分析的时候可以考虑去做XIST与Y染色体相关基因的相关性散点图，证明性别对实验没有影响。
细胞周期分数计算
通过比较实验得到的基因与不同细胞周期的相关基因的平均表达量的差异，从而对给定的基因列表进行打分，由此评估对应样本细胞所在的细胞周期。具体原理，很感兴趣的可以自行百度。
# 首先进行归一化，在使用自带CellCycleScoring函数进行细胞周期评分的计算
data.filt = NormalizeData(data.filt)
data.filt <- CellCycleScoring(object = data.filt, g2m.features = cc.genes$g2m.genes,
s.features = cc.genes$s.genes)
# 绘制小提琴图进行分析
VlnPlot(data.filt, features = c("S.Score", "G2M.Score"), group.by = "orig.ident",
ncol = 4, pt.size = 0.1)
Doublet去除
双细胞指的是在测量的过程中，进行测序的液滴中包含了两个细胞导致所测得的基因数量异常，所以需要去除这些双细胞的测序结果，尽量不对后面的聚类分析造成影响。
正如 Chromium 用户指南中所指出的，10X genomics 双细胞率大约如下：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-YGrYhjNV-1661081172343)(/Volumes/Flower0501/研究学习目录/单细胞转录组/分析流程/Seurat下游分析/10x_doublet_rate.png)]
大多数Doublet检测器通过合并细胞计数来模拟Doublet，并将Doublet预测为与模拟Doublet具有相似含量的细胞。大多数此类包需要对数据集中预期Doublet的数量/比例进行假设。比如，若你使用的数据是二次抽样的，但原始数据集每个样本包含大约 5 000 个细胞，因此我们可以假设它们加载了大约 9 000 个细胞，并且应该有大约 4% 的Doublet率。 理想情况下，应分别对每个样本进行Doublet预测，尤其是当你的不同样本具有不同比例的细胞类型时。在这种情况下，数据是二次抽样的，因此每个样本的细胞非常少，且所有样本的基因数均是经过排序的，因此可以一起运行它们。
# 在这里，我们将使用 DoubletFinder 来预测双重细胞。但在进行双峰检测之前，我们需要进行基因含量标准化，可变基因选择和 pca，以及UMAP
data.filt <- FindVariableFeatures(data.filt, verbose = F)
data.filt <- ScaleData(data.filt, vars.to.regress = c("nFeature_RNA", "percent_mito"),
verbose = F)
data.filt <- RunPCA(data.filt, verbose = F, npcs = 20)
data.filt <- RunUMAP(data.filt, dims = 1:10, verbose = F)
# DoubletFinder
# 可以使用 paramSweep 运行参数优化，但暂时跳过。
# sweep.res <- paramSweep_v3(data.filt)
# sweep.stats <-summarizeSweep(sweep.res, GT = FALSE) bcmvn <- find.pK(sweep.stats)
# barplot(bcmvn$BCmetric, names.arg = bcmvn$pK, las=2)
# 定义双重细胞细胞的预期数量。
nExp <- round(ncol(data.filt) * 0.04) # expect 4% doublets
data.filt <- doubletFinder_v3(data.filt, pN = 0.25, pK = 0.09, nExp = nExp, PCs = 1:10)
# DoubletFinder预测的名称可以更改，因此提取正确的列名称。
DF.name <- colnames(data.filt@meta.data)[grepl("DF.classification", colnames(data.filt@meta.data))]
# 绘图查看双细胞所在的区域
cowplot::plot_grid(ncol = 2, DimPlot(data.filt, group.by = "orig.ident") + NoAxes(),
DimPlot(data.filt, group.by = DF.name) + NoAxes())
按照一般逻辑，含有双细胞的数据应该比正常的数据含有更多的基因的数量，此时可以从过绘制基因数量的小提琴图进行检验。
VlnPlot(data.filt, features = "nFeature_RNA", group.by = DF.name, pt.size = 0.1)
# 最终出去所有的Doublet，得到最后的质控数据
data.filt = data.filt[, data.filt@meta.data[, DF.name] == "Singlet"]
dim(data.filt)
# [1] 18133 6374
保存数据
dir.create("../results", showWarnings = F)
saveRDS(data.filt, "../results/seurat_covid_qc.rds")
本文含有隐藏内容，请 开通VIP 后查看