基于单细胞参考图谱的分解：
构建参考图谱：首先，需要构建一个参考图谱，其中包含每种细胞类型的基因表达模式。参考图谱可以从已知的单细胞RNA-seq数据中获取。从单细胞数据中识别出不同的细胞类型，并计算每个细胞类型的平均表达水平，作为参考表达谱。
选择标记基因：只选择那些在至少一个细胞类型中有显著差异表达的基因作为标记基因，过滤掉方差为零、在bulk数据中未表达或者是线粒体基因的基因。
转换bulk数据：对于每个基因，根据单细胞数据和bulk数据之间的线性关系，学习一个转换函数，将bulk数据转换为与参考表达谱更一致的形式。如果有部分样本同时有单细胞和bulk数据，就用这些样本来拟合线性回归模型；如果没有，就用单细胞数据的均值和方差来估计转换函数。
估计细胞类型比例：对于每个样本，使用转换后的bulk数据和参考表达谱，求解一个带有非负约束和和为一约束的最小二乘问题，得到每个细胞类型的比例估计。其中，使用非负最小二乘回归方法来进行反卷积，这意味着在估计细胞类型比例时，结果必须保持非负且最小误差。非负最小二乘回归是一种优化算法，通过最小化实际数据与估计结果之间的误差来实现这一点。为了确保估计的细胞类型比例总和为一，还会添加一个约束条件。这样可以保证细胞类型比例之和等于100%，表示完整的细胞组成。
基于标记基因的分解：
获取标记：对于单细胞数据，需要每种细胞类型的一组marker基因，以及差异分析所得的倍数变化，以表示每个标记的特异性程度。差异倍数在后期可以作为一个权重，用来调整每个基因在PCA分析中的贡献。
选择标记基因：选择一组已知的marker基因，这些基因在某些细胞类型中有显著的差异表达。
估计细胞类型比例：使用加权PCA方法，使用已知的marker基因来估计相对细胞类型丰度。具体来说，对于每个细胞类型，选取其marker基因的子矩阵，这个子矩阵包含了所有样本在这些marker基因上的表达水平。然后对每个子矩阵进行加权PCA分析，第一主成分能解释最多变异性的方向，因此可以看作该细胞类型的代表性表达模式。所以，最终提取其第一主成分作为该细胞类型的比例估计。
1. 基于单细胞参考图谱的分解：
a. 准备bulk RNA-seq表达谱（例如：bulk_profile.txt）
b. 准备单细胞Seurat object的RData文件，包含count matrix信息和细胞类型注释信息（例如：sc_obj.RData）
c. 撰写代码并运行（以下代码写入一个R脚本，例如：Bisque_Reference.R）：
library(Biobase)library(BisqueRNA)library(Seurat)## 读取bulk表达谱，bulk_profile <- read.table("bulk_profile.txt", sep = "\t", quote = "\t", row.names = 1, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE)## 将bulk表达谱转换为ExpressionSet格式bulk.eset <- ExpressionSet(assayData = as.matrix(bulk_profile))## 读取单细胞Seurat object数据load("sc_obj.RData") # 此处RData里面存储的数据变量名称叫“sc_obj”，根据实际存储结果进行调用/重新赋值## 将单细胞数据转换为ExpressionSet格式# sc.eset <- SeuratToExpressionSet(sc_obj, delimiter="-", position=2, version="v3")# 如果Seurat object中以Idents的方式进行过细胞类型注释，并且格式标准，那么上述SeuratToExpressionSet的方式可以直接转换。如果转换后，有信息缺失，下游分析有问题，那么可以使用下述方法进行必要信息提取后，进行格式转换（注意：不同操作者/项目的单细胞Seurat object中存储的格式和信息可能有不同，提取信息时需要根据实际存储结构进行针对性提取）：sc.counts.matrix <- as.matrix(GetAssayData(sc_obj, assay = "RNA", slot = "counts"))sample.ids <- colnames(sc.counts.matrix)sc.pheno <- data.frame(check.names=F, check.rows=F,stringsAsFactors=F,row.names=sample.ids,SubjectName=sc_obj$orig.ident, # 每个细胞所属的样本（与samples.ids细胞顺序一致）cellType=Idents(sc_obj)) # 每个细胞所属的细胞类型（与samples.ids细胞顺序一致）sc.meta <- data.frame(labelDescription=c("SubjectName","cellType"),row.names=c("SubjectName","cellType"))sc.pdata <- new("AnnotatedDataFrame",data=sc.pheno,varMetadata=sc.meta)sc.eset <- Biobase::ExpressionSet(assayData=sc.counts.matrix,phenoData=sc.pdata)## 计算细胞类型比例res <- ReferenceBasedDecomposition(bulk.eset, sc.eset, use.overlap=FALSE)# 如果use.overlap设置为TRUE，那么函数将只使用在参考基因表达谱和目标基因表达谱中都存在的样本来进行分解。## 提取获取细胞类型比例，保存文件estimated.cell.proportions <- res$bulk.propswrite.table(estimated.cell.proportions, file = "estimated.cell.proportions.txt", sep = "\t", row.names = TRUE, col.names = NA, quote = FALSE)
2. 基于标记的分解：
a. 准备bulk RNA-seq表达谱（例如：bulk_profile.txt）
b. 准备marker基因表格（例如：markers.txt）
c. 撰写代码并运行（以下代码写入一个R脚本，例如：Bisque_Marker.R）
library(Biobase)library(BisqueRNA)bulk_profile <- read.table("bulk_profile.txt", sep = "\t", quote = "\t", row.names = 1, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE)markers <- read.table("markers.txt", sep = "\t", quote = "\t", header = TRUE, check.names = FALSE, stringsAsFactors = FALSE)res <- MarkerBasedDecomposition(bulk.eset, markers, weighted = FALSE)estimated.cell.proportions <- res$bulk.propswrite.table(estimated.cell.proportions, file = "estimated.cell.proportions.txt", sep = "\t", row.names = TRUE, col.names = NA, quote = FALSE)
Jew, B., Alvarez, M., Rahmani, E. et al. Accurate estimation of cell composition in bulk expression through robust integration of single-cell information. Nat Commun 11, 1971 (2020). https://doi.org/10.1038/s41467-020-15816-6
CytoTRACE：细胞分化潜能分析|生信开发实战
Visium HD数据分析之Bin2Cell |生信开发实战
ROC分析介绍|生信开发实战
如果使用find_circ来鉴定circRNA|生信开发实战
本文系联川生物公众号原创文章，未经授权禁止转载，侵权必究！ 扫描下方二维码 点分享
点点赞
点在看
