一、seurat执行差异表达分析
根据单细胞数据预处理的方式不同（lognormalize和sctransform），执行差异表达分析代码有所不同，这个需要注意。我当时请教过生信分析同行，也跟人交流过，再此记录一下。
1.1 lognormalize预处理
seurat官网给出的标准分析流程，通过lognormalize预处理，归一化的数据储存在seurat_obj[['RNA']]@data，我们用归一化后的data数据进行差异基因分析。对于多样本整合分析，同样是在RNA assay上做差异表达分析。
# Normalizing the data
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
# Identification of highly variable features
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Scaling the data
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
pbmc <- RunUMAP(pbmc, dims = 1:10)
# find markers for every cluster compared to all remaining cells, report only the positive
DefaultAssay(pbmc) <- "RNA"
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>%
group_by(cluster) %>%
slice_max(n = 2, order_by = avg_log2FC)
1.2 sctransform预处理
在前面sctransform的文章中，通过sctransform预处理，标准化的结果保存在一个新的assay中（默认命名为 SCT），包含：
- seurat_obj[['SCT']]@counts：矫正后的counts；
- seurat_obj[['SCT']]@data：矫正后的counts通过 log1p(counts)处理后的data;
- seurat_obj[['SCT']]@scale.data：pearson残差；
问题：通过sctransform预处理，我们究竟怎么选择assay进行差异基因分析？
如果依旧应用seurat_obj[['RNA']]@data进行分析，显然行不通，在seurat官网sctransform教程中没有明确的指导，这也是让我当时感到很困惑的一点。
在sctransform的教材中，作者提到
You can use the corrected log-normalized counts for differential expression and integration. However, in principle, it would be most optimal to perform these calculations directly on the residuals (stored in the scale.data slot) themselves. This is not currently supported in Seurat v3, but will be soon.
我们可以这样理解，校正的log-normalized归一化计数（seurat_obj[['SCT']]@data）可以用于差异表达分析。但是，原则上，最好直接对残差（seurat_obj[['SCT']]@scale.data）本身执行这些计算。目前，Seurat v3 不支持此功能，但很快就会支持。
但是在Seurat v4，我们依旧没有等到seurat官网对此分析流程的官宣。
但这也说明了seurat_obj[['SCT']]@data是可行的。
DefaultAssay(seurat_obj) <- "SCT"
all.markers.SCT <- FindAllMarkers(seurat_obj, only.pos = T, min.pct = 0.5, logfc.threshold = 0.5)
作者在discussions中也给出了一段说明：
We had anticipated extending Seurat to actively support DE using the pearson residuals of sctransform, but have decided not to do so. In some cases, Pearson residuals may not be directly comparable across different datasets, particularly if there are batch effects that are unrelated to sequencing depth. While it is possible to correct these differences using the SCTransform-based integration workflow for the purposes of visualization/clustering/etc., we do not recommend running differential expression directly on Pearson residuals. Instead, we recommend running DE on the standard RNA assay.
我们曾期望使用sctransform的皮尔逊残差扩展 Seurat 以积极支持DE，但我们决定不这样做。 在某些情况下，Pearson 残差可能无法在不同数据集之间直接进行比较，尤其是当存在与测序深度无关的批次效应时。 虽然可以使用基于SCTransform 的分析工作流来矫正这些实验差异以实现可视化/聚类等目的，但我们不建议直接在Pearson残差上运行差异表达分享。 相反，我们建议在标准RNA assay中运行DE分析。
我们可以这样理解，Satija团队本来期望用sctransform的皮尔逊残差进行DE分析，但是他们发现，Pearson残差可能无法在不同数据集之间直接进行比较，决定不这么做。基于SCTransform 的分析工作流可以剔除实验差异以实现可基因可视化和聚类等功能。但是他们建议在标准RNA assay中运行DE分析。
我们用一个具体的案例来实践一下。
step1：下载一个10X单细胞数据，进行基本的数据质控；
library(Seurat)
library(ggplot2)
library(SingleR)
library(dplyr)
library(celldex)
library(RColorBrewer)
library(SingleCellExperiment)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# part1: Basic quality control and filtering
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
adj.matrix <- Read10X("data/update/soupX_pbmc10k_filt")
srat <- CreateSeuratObject(adj.matrix,project = "pbmc10k")
head(meta)
# orig.ident nCount_RNA nFeature_RNA
# AAACCCACATAACTCG-1 pbmc10k 22196 4734
# AAACCCACATGTAACC-1 pbmc10k 7630 2191
# AAACCCAGTGAGTCAG-1 pbmc10k 21358 4246
# AAACCCAGTGCTTATG-1 pbmc10k 857 342
# AAACGAACAGTCAGTT-1 pbmc10k 15007 4075
# AAACGAACATTCGGGC-1 pbmc10k 9855 2285
srat[["percent.mt"]] <- PercentageFeatureSet(srat, pattern = "^MT-")
srat[["percent.rb"]] <- PercentageFeatureSet(srat, pattern = "^RP[SL]")
# 去除双细胞
doublets <- read.table("data/update/scrublet_calls.tsv",header = F,row.names = 1)
colnames(doublets) <- c("Doublet_score","Is_doublet")
srat <- AddMetaData(srat,doublets)
head(srat[[]])
srat[['QC']] <- ifelse(srat@meta.data$Is_doublet == 'True','Doublet','Pass')
srat[['QC']] <- ifelse(srat@meta.data$nFeature_RNA < 500 & srat@meta.data$QC == 'Pass','Low_nFeature',srat@meta.data$QC)
srat[['QC']] <- ifelse(srat@meta.data$nFeature_RNA < 500 & srat@meta.data$QC != 'Pass' & srat@meta.data$QC != 'Low_nFeature',paste('Low_nFeature',srat@meta.data$QC,sep = ','),srat@meta.data$QC)
srat[['QC']] <- ifelse(srat@meta.data$percent.mt > 15 & srat@meta.data$QC == 'Pass','High_MT',srat@meta.data$QC)
srat[['QC']] <- ifelse(srat@meta.data$nFeature_RNA < 500 & srat@meta.data$QC != 'Pass' & srat@meta.data$QC != 'High_MT',paste('High_MT',srat@meta.data$QC,sep = ','),srat@meta.data$QC)
table(srat[['QC']])
# Doublet High_MT High_MT,Low_nFeature High_MT,Low_nFeature,Doublet Pass
# 546 548 275 1 8824
step2：数据进行SCTransform预处理
srat <- subset(srat, subset = QC == 'Pass')
dim(srat)
# [1] 36601 8824
srat <- SCTransform(srat, method = "glmGamPoi", ncells = 8824,
vars.to.regress = c("percent.mt"), verbose = F)
srat
接下来，我们比较下，RNA assay和SCT assay矩阵数据的一致性。
# RNA/SCT counts
gene_exp_matrix_RNA_counts <- srat@assays$RNA@counts
gene_exp_matrix_RNA_counts[1:20,1:20]
gene_exp_matrix_SCT_counts <- srat@assays$SCT@counts
gene_exp_matrix_SCT_counts[1:20,1:20]
testthat::expect_equal(gene_exp_matrix_RNA_counts,gene_exp_matrix_SCT_counts)
# Error: `gene_exp_matrix_RNA_counts` not equal to `gene_exp_matrix_SCT_counts`
# RNA/SCT data
gene_exp_matrix_RNA_data <- srat@assays$RNA@data
gene_exp_matrix_RNA_data[1:20,1:20]
gene_exp_matrix_SCT_data <- srat@assays$SCT@data
gene_exp_matrix_SCT_data[1:20,1:20]
testthat::expect_equal(gene_exp_matrix_RNA_data,gene_exp_matrix_SCT_data)
# Error: `gene_exp_matrix_RNA_data` not equal to `gene_exp_matrix_SCT_data`.
testthat::expect_equal(gene_exp_matrix_RNA_counts,gene_exp_matrix_RNA_data)
RNA assay中的data和SCT assay中的data矩阵数据不相同，但是RNA assay中的data等于RNA assay中的counts，也就是说通过SCTransform预处理，seurat_obj[['RNA']]@data为原始矩阵数据，是不能进行下游分析的。
step3：完成后面的降维聚类
# clustering
srat <- RunPCA(srat, verbose = F)
srat <- RunUMAP(srat, dims = 1:30, verbose = F)
srat <- FindNeighbors(srat, dims = 1:30, verbose = F)
srat <- FindClusters(srat, verbose = F)
table(srat[[]]$seurat_clusters)
# 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
# 1720 1503 1214 1090 955 326 324 315 271 250 246 239 135 104 97 35
step4：进行差异基因分析
推荐使用在RNA assay上做差异表达分析，而不是 SCTransform。
切换到RNA assay，对原始的基因counts矩阵进行NormalizeData和ScaleData，获取归一化的seurat_obj[['RNA']]@data后，进行差异基因分析。
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# part3: Differential expression and marker selection
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DefaultAssay(srat) <- "RNA"
srat <- NormalizeData(srat)
srat <- FindVariableFeatures(srat, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(srat)
srat <- ScaleData(srat, features = all.genes)
all.markers.RNA <- FindAllMarkers(srat, only.pos = T, min.pct = 0.5, logfc.threshold = 0.5)
dim(all.markers.RNA)
# [1] 4026 7
table(all.markers.RNA$cluster)
# 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
# 541 136 129 75 531 448 65 182 261 155 139 142 250 524 121 327
另外，我们使用SCT assay进行差异表达分析，与上面的结果做对比。
DefaultAssay(srat) <- "SCT"
all.markers.SCT <- FindAllMarkers(srat, only.pos = T, min.pct = 0.5, logfc.threshold = 0.5)
dim(all.markers.SCT)
# [1] 3458 7
table(all.markers.SCT$cluster)
# 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
# 464 118 127 65 453 398 50 151 228 134 122 138 209 459 101 241
接下来，我们可视化每个cluster的marker基因差异。
gene_df <- data.frame()
for (i in levels(all.markers.RNA$cluster)){
share_num <- length(intersect(x=all.markers.RNA[all.markers.RNA$cluster==i,]$gene, y = all.markers.SCT[all.markers.SCT$cluster==i,]$gene))
rna_num <- length(setdiff(x=all.markers.RNA[all.markers.RNA$cluster==i,]$gene, y = all.markers.SCT[all.markers.SCT$cluster==i,]$gene))
sct_num <- length(setdiff(x=all.markers.SCT[all.markers.SCT$cluster==i,]$gene, y = all.markers.RNA[all.markers.RNA$cluster==i,]$gene))
gene_df <- rbind(gene_df, data.frame("cluster"=i, "share_Gene"=share_num, "RNA_Gene" = rna_num, "SCT_Gene"=sct_num))
}
gene_data <- reshape2::melt(gene_df,id.vars=c("cluster"), variable.name="Type",value.name="gene_num")
head(gene_data)
# cluster Type gene_num
# 1 0 share_Gene 464
# 2 1 share_Gene 118
# 3 2 share_Gene 125
# 4 3 share_Gene 65
# 5 4 share_Gene 451
# 6 5 share_Gene 393
gene_data$cluster <- factor(gene_data$cluster,levels = order(as.numeric(unique(gene_data$cluster))))
ggplot(data=gene_data, aes(x=cluster, y=gene_num, fill=Type))+ geom_bar(stat="identity") + theme_bw()
用SCT assay查找的差异基因记录为all.markers.SCT，用RNA assay查找的差异基因记录为all.markers.RNA，对每个cluster的差异基因进行分类："share_Gene"表示markers.SCT和markers.RNA共有的基因；"RNA_Gene"表示markers.RNA独有的基因；"SCT_Gene"表示markers.SCT独有的基因。从下图可以看出，cluster0-15，红色条块的"share_Gene"占主要部分，"RNA_Gene"和"SCT_Gene"占小部分。这个例子说明，用SCT assay查找的差异基因和用RNA assay查找的差异基因大致是相同的。
1.3 小结
问：sctransform预处理后，如何进行差异表达分析？
答：虽然普遍认为，SCTransform和RNA assay计算出的差异基因大致是一致的。
但是作者的答复是：推荐使用在RNA assay上做差异表达分析，而不是SCTransform。切记一定要进行NormalizeData和ScaleData标准化处理。
也就是说，SCTransform预处理适合降维聚类，基因可视化；做差异表达分析，lognormalize预处理更适宜一些。
参考代码链接：https://pan.baidu.com/s/1wJFhPQoKSKTjqt3Qhd4a7w
提取码：1234
