R语言学习笔记：从数据清洗到高级统计学和生物信息学 - 寻找不同样本类型间同一细胞类型内的差异基因

[giscus[bot]]

R语言学习笔记：从数据清洗到高级统计学和生物信息学 - 寻找不同样本类型间同一细胞类型内的差异基因

https://djhcod.github.io/r-notes/single_cell/seurat/de_vignette.html

[ShengXincainiao]

我想请教您，可不可以将数据提取出来重新构建DEseq对象，进行差异分析呢？如果这样做跟Findmarkers有什么区别吗？感谢！

[djhcod]

个人认为应该是可以的。Findmarkers就是一个包装了很多差异分析方法的对单细胞分析场景进行了优化的函数，其中的test.use参数即为指定差异分析方法的参数。默认的是“Wilcoxon Rank Sum test”，但是你也可以指定为"DESeq2"，即通过调用DESeq2包基于“negative binomial distribution”来进行差异分析，具体可以详见Findmarkers的帮助文档。同时，从Findmarkers的源码中也可以看到Findmarkers如何调用DESeq2包进行差异分析的具体代码（如下所示）。

Findmarkers部分源码摘录：

DESeq2DETest <- function(
data.use,
cells.1,
cells.2,
verbose = TRUE,
...
) {
if (!PackageCheck('DESeq2', error = FALSE)) {
stop("Please install DESeq2 - learn more at https://bioconductor.org/packages/release/bioc/html/DESeq2.html")
}
CheckDots(..., fxns = 'DESeq2::results')
group.info <- data.frame(row.names = c(cells.1, cells.2))
group.info[cells.1, "group"] <- "Group1"
group.info[cells.2, "group"] <- "Group2"
group.info[, "group"] <- factor(x = group.info[, "group"])
group.info$wellKey <- rownames(x = group.info)
dds1 <- DESeq2::DESeqDataSetFromMatrix(
countData = data.use,
colData = group.info,
design = ~ group
)
dds1 <- DESeq2::estimateSizeFactors(object = dds1)
dds1 <- DESeq2::estimateDispersions(object = dds1, fitType = "local")
dds1 <- DESeq2::nbinomWaldTest(object = dds1)
res <- DESeq2::results(
object = dds1,
contrast = c("group", "Group1", "Group2"),
alpha = 0.05,
...
)
to.return <- data.frame(p_val = res$pvalue, row.names = rownames(res))
return(to.return)
}

[ShengXincainiao]

非常详细的指导，感谢您！赞赞赞

[alice683]

您好，遇见一个问题想请教您，我有4个样本，CK1, CK2, ST1, ST2, 我进行了数据整合（细胞分类分了7类），整合后进行差异分析，差异分析报错了，我的数据分析流程大致如下（sc_inted，整合后的数据）：
head(sc_inted, 3)
orig.ident nCount_RNA nFeature_RNA celltype rep log10GenesPerUMI sample
CKrep1_AAACCCAAGCATGCAG-1 CK_rep1 1185 1025 Cortex rep1 0.9795054 CK
CKrep1_AAACCCAAGCTGAAAT-1 CK_rep1 867 753 Cortex rep1 0.9791614 CK
CKrep1_AAACCCACAACAGCCC-1 CK_rep1 1350 1076 Cortex rep1 0.9685268 CK

table(sc_inted$rep, sc_inted$sample)

CK ST

#rep1 12656 21215
#rep2 11413 18285

pseudo_sc_inted <- AggregateExpression(sc_inted, assays = "RNA", return.seurat = T,
group.by = c("sample", "rep", "celltype"))
pseudo_sc_inted
#An object of class Seurat
#30036 features across 28 samples within 1 assay
#Active assay: RNA (30036 features, 0 variable features)

3 layers present: counts, data, scale.data

pseudo_sc_inted$sample.celltype <- paste(pseudo_sc_inted$sample,
pseudo_sc_inted$celltype,
sep = "_")

head(pseudo_sc_inted)
orig.ident sample rep celltyper sample.celltype
CK_rep1_Cortex CK_rep1_Cortex CK rep1 Cortex CK_Cortex
CK_rep1_Initials CK_rep1_Initials CK rep1 Initials CK_Initials
CK_rep1_Root-cap CK_rep1_Root-cap CK rep1 Root-cap CK_Root-cap
CK_rep1_Phloem CK_rep1_Phloem CK rep1 Phloem CK_Phloem
CK_rep1_Stele CK_rep1_Stele CK rep1 Stele CK_Stele
CK_rep1_Epidermis CK_rep1_Epidermis CK rep1 Epidermis CK_Epidermis
CK_rep1_Endodermis CK_rep1_Endodermis CK rep1 Endodermis CK_Endodermis
CK_rep2_Cortex CK_rep2_Cortex CK rep2 Cortex CK_Cortex
CK_rep2_Initials CK_rep2_Initials CK rep2 Initials CK_Initials
CK_rep2_Root-cap CK_rep2_Root-cap CK rep2 Root-cap CK_Root-cap

Idents(pseudo_sc_inted) <- "sample.celltype"
library(DESeq2)
bulk.Cortex.de <- FindMarkers(pseudo_sc_inted,
ident.1 = "ST_Cortex",
ident.2 = "CK_Cortex",
test.use = "DESeq2")
#Error in ValidateCellGroups(object = object, cells.1 = cells.1, cells.2 = cells.2, :

Cell group 1 has fewer than 3 cells

请问您遇见过这种问题吗？如何解决呢？先谢谢啦~~~

[djhcod]

之前没有遇到过这个问题。可以检查一下sample.celltype各个组别里面的细胞数

[alice683]

好的谢谢，检查了sample.celltype各个组别的细胞数，最少的在138个，我再找找原因，谢谢呀~~

[alice683]

我找到问题所在了，我的数据有2次重复，执行pseudobulk analysis时，需要添加一个参数 min.cells.group = 2 ，这样可以运行成功。代码如下：
bulk.Cortex.de <- FindMarkers(pseudo_sc_inted,
ident.1 = "ST_Cortex",
ident.2 = "CK_Cortex",
test.use = "DESeq2",
min.cells.group = 2)

[djhcod]

学习了👍

[alice683]

你好，请问近期可以不可以出一个 scSTAR分析讲解的帖子呀？

[vaelgc7]

老师，您好，想请教一下我下面的代码是正确的吗？data_merged是将16个样本合并为一个seurat对象，然后使用SCTransform归一化，使用harmony去除批次效应得到的seurat对象。在进行decontx去除RNA污染时，默认的assay为SCT assay，counts数据使用的RNA assay的counts，跑了5天都没跑出结果。
另外，还想请教一下，我参考您的文章，也是想只使用contamination，不使用矫正后的矩阵，那么我应该使用过滤后的counts数据重新建一个对象还是使用过滤掉污染细胞的SCT assay呢？我思考的是还是使用SCT assay，因为之前做过去除doublet和harmony的操作，如果使用counts那么这些操作都相当于没有进行，您觉得呢？

sel.clust = "SCT_snn_res.0.3"
table(data_merged[[sel.clust]][,1])
raw_counts <- GetAssayData(data_merged, assay = "RNA", slot = "counts")
#36601 x 96696 sparse Matrix of class "dgCMatrix"

cluster_labels <- data_merged[[sel.clust]][,1]

filtered_counts <- raw_counts[rowSums(raw_counts) > 0, ]
#31018 x 96696 sparse Matrix of class "dgCMatrix"

decont_res <- decontX(counts = filtered_counts, z = cluster_labels, verbose = TRUE, seed = 123)

添加 contamination 分数

data_merged$Contamination <- decont_res$contamination

Filter 高污染细胞 (Contamination < 0.2)

low_con_seurat <- subset(data_merged, subset = Contamination < 0.2)
paste0("通过decontX过滤掉了", ncol(data_merged) - ncol(low_con_seurat ), "个细胞")

data_clean <- CreateSeuratObject(
counts = low_con_seurat [["RNA"]]$counts,
meta.data = low_con_seurat @meta.data)
saveRDS(data_clean, "../output/4-decontX/data_Harmony_decontX_clean.rds")

[djhcod]

尝试①在harmony整合之前执行decontX（首选）；②按照decontX官方vignettes的做法将Seurat对象转换成SingleCellExperiment对象后执行decontX，然后将其重新添加到 Seurat 对象中。
无论是doublet还是环境污染的评估及去除实际上都还是属于细胞的质控，所以这个步骤理论上应该是在整个单细胞处理管线中比较靠前完成的，所以我个人感觉应该使用过滤后的counts数据重新建质控后的Seurat对象，然后重新进行后续的归一化、整合、降维、分群等操作。