单细胞DEGs、markers、celltypes | 专辑完结篇 by 生信菜鸟团

这周转录组周更憋个大招，细心的小伙伴可能也发现了，初学者暑期搞定单细胞这个专辑也是我写的，暑期磨磨蹭蹭总归也还是出来了点东西，但是开学后就没有时间搞单细胞了，和目前做的方向也不相干，一些我想做的，比如细胞通讯、拟时序分析，都还没来得及，不知道有没有朋友真的暑期实现了我立的flag“搞定”了单细胞。所以这期周更算是这个单细胞专辑的完结篇，我将系统地总结：单细胞水平样本samples间的差异分析、clusters间的差异分析鉴定marker、一些检索细胞亚型markers的数据库以及如何使用一些打分算法结合markers辅助确定细胞亚型，这四个问题。

内容很长，干货很多，小张让我分成两期推送，我想着反正最后一篇单细胞了，那就毕其功于一役吧。上期推文使用snakemake管理转录组流程，有一些小伙伴在评论区给我提出问题，这是我非常欣喜看到的，因为我本身是个生信菜狗，很多时候曾老师不一定有时间看我的推送帮我检查，我更多的时候是抱着学习的态度分享笔记，而不是教程。

单细胞水平样本间的差异分析

主要参考：

生物学重复的不同条件下的scrna-seq数据的差异基因表达分析（https://www.10xgenomics.com/resources/analysis-guides/differential-gene-expression-analysis-in-scrna-seq-data-between-conditions-with-biological-replicates）

需要注意的是，这部分差异表达分析不应与下一部分，计算来自两个簇的细胞之间的差异表达基因以用于标记基因鉴定或细胞类型注释混淆，它们是两种不同目的的不同分析

background

Samples as biological replicates for DGE between conditions

对于任何给定的细胞类型和条件，来自一个样本的细胞可能比来自不同样本的细胞更相似

Multi-condition DGE analysis approaches in scRNAseq data

为了获得有效的假设检验结果，已经开发了几种DGE算法来解释差异表达分析中来自同一样本的细胞的分组性质。这些算法包括拟合混合效应模型、对样本内细胞计数求和（pseudobulk）以及测试分布差异。

General analysis workflow

需要关注批次效应问题，就这个问题我们前面已经有所讨论harmony、不harmony，这是个问题

个人感觉单细胞是否需要去除批次效应比bulk更难判断，因为单细胞我们已经进入到了细胞亚型的层次，任何校正影响都很大，bulk这个时候对于case、control来说更“粗枝大叶”一些

一般来说，原始计数应用于DGE分析，而不是标准化或批量校正的数据

但是在细胞聚类鉴定亚型时往往需要将原始单细胞RNA测序数据标准化，以消除技术差异

参考 Batch Effect Correction Analysis Guide：

https://www.10xgenomics.com/resources/analysis-guides/introduction-batch-effect-correction

Common multi-condition DGE analysis approaches

Mixed-effects model混合效应模型：处理相关数据的经典方法是使用混合效应模型，该模型使用随机效应来解释观测值之间的相关性。感兴趣的条件/组被建模为固定效应，并且样本特定的随机截距对来自相同样本（和细胞类型）的细胞之间的相关性进行建模。

Pseudobulk：使用标准统计软件在单细胞数据（数千次观测）的规模上拟合混合效应模型的计算量非常大。一种更简单的方法是总结每个样本的细胞类型内的UMI计数，并使用为bulk RNAseq开发的方法进行多条件DE。

Differential distribution test 差异分布测试：单细胞数据使我们能够估计和了解每个样本中相同细胞类型的细胞中基因表达的分布，而不是测试每个基因在不同条件下的平均表达差异。正如前两种方法所解决的那样，这些分布可能具有不同的手段，但差异分布测试也可以检测到更细微的差异，例如方差或模式数量的变化。我们可以测试这些分布在条件/组之间是否不同。

关于不同类型方法的介绍，参考：

各种方法代表性工具描述性总结表格

https://www.10xgenomics.com/resources/analysis-guides/differential-gene-expression-analysis-in-scrna-seq-data-between-conditions-with-biological-replicates

下面我就将挑选每种方法的代表性工具进行展示

在正式分析前，我还想提醒一个装包的问题，有的bioC的包在线安装会出现网络问题，我们往往会去bioconductor官网下载源码安装，有的时候bioC3.16 对应的R包源码无法访问，比如下面即将使用的muscat

我这里提供两种解决方法：

bioC3.17对于的R包并不是非4.3的R不能使用，只是可能会出现不能用的情况，主要还是去看看有哪些依赖
如果确实在线、bioconductor源码都无法访问，最后就是去github看看。但是github上当前托管的仓库代码应该是development的，并不是我们希望使用的stable版本，这个时候可以参考具体要使用的版本号，去仓库历史里返回之前commit的版本下载安装即可，参考：如何安装Github上特定版本号（每次提交的唯一标识）的R包？

muscat

参考：

单细胞差异基因分析：pseudobulk（muscat包）
Differential state analysis with muscat [http://www.bioconductor.org/packages/release/bioc/vignettes/muscat/inst/doc/analysis.html]

来自8名狼疮患者的10x基于液滴的scRNA-seq PBCM数据

加载数据

数据来源：

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96583

使用官方教程提供的方法无法连接访问

去GEO查看提供的数据：

发现好像每个样本没有对应的genes信息，汇总的genes文件也不对，无法用 read10X读取：

于是回到我的本地小电脑开启魔法，使用教程提供的方法获取数据：

# install.packages("interactiveDisplayBase_1.38.0.tar.gz",repos = NULL)
# install.packages("AnnotationHub_3.8.0.tar.gz",repos = NULL)
# install.packages("ExperimentHub_2.8.1.tar.gz",repos = NULL)

library(ExperimentHub)
eh <- ExperimentHub()
query(eh, "Kang")

options()$repos
# CRAN 
# "https://cran.rstudio.com/" 
# attr(,"RStudio")
# [1] TRUE

# install.packages("SingleCellExperiment_1.22.0.tar.gz",repos = NULL)
# install.packages("muscData_1.14.0.tar.gz",repos = NULL)

sce <- eh[["EH2259"]]
save(sce,file = "sce.Rdata")

保存好sce传给服务器

我这里把R的下载镜像从清华退回到了CRAN，因为我发现使用镜像会导致魔法失效

包含来自8名狼疮患者的10x基于液滴的scRNA-seq PBCM数据，这些数据在用IFN-β治疗6小时之前和之后获得

library(ExperimentHub)
eh <- ExperimentHub()
query(eh, "Kang")

ExperimentHub with 3 records
# snapshotDate(): 2022-10-31
# $dataprovider: NCI_GDC, GEO
# $species: Homo sapiens
# $rdataclass: character, SingleCellExperiment, BSseq
# additional mcols(): taxonomyid, genome, description, coordinate_1_based, maintainer,
#   rdatadateadded, preparerclass, tags, rdatapath, sourceurl, sourcetype 
# retrieve records with, e.g., 'object[["EH1661"]]' 

           title                                         
  EH1661 | Whole Genome Bisulfit Sequencing Data for 47 samples
  EH1662 | Whole Genome Bisulfit Sequencing Data for 47 samples
  EH2259 | Kang18_8vs8

数据集包含>35，000个基因和~29，000个细胞

EH2259 | Kang18_8vs8

sce <- eh[["EH2259"]]

load("sce.Rdata")
sce

class: SingleCellExperiment 
dim: 35635 29065 
metadata(0):
assays(1): counts
rownames(35635): MIR1302-10 FAM138A ... MT-ND6 MT-CYB
rowData names(2): ENSEMBL SYMBOL
colnames(29065): AAACATACAATGCC-1 AAACATACATTTCC-1 ... TTTGCATGGTTTGG-1 TTTGCATGTCTTAC-1
colData names(5): ind stim cluster cell multiplets
reducedDimNames(1): TSNE
mainExpName: NULL
altExpNames(0):

预处理

去除未检测到的基因
去除检测到很少或许多基因的细胞
去除非常低表达的基因
计算标准化表达式值以进行可视化

# remove undetected genes
sce <- sce[rowSums(counts(sce) > 0) > 0, ]
dim(sce)

使用perCellQCMetrics来计算各种每个细胞的质量控制指标，并如上所述继续过滤细胞和基因：

# calculate per-cell quality control (QC) metrics
library(scater)
qc <- perCellQCMetrics(sce)

# remove cells with few or many detected genes
ol <- isOutlier(metric = qc$detected, nmads = 2, log = TRUE)
sce <- sce[, !ol]
dim(sce)

# remove lowly expressed genes
sce <- sce[rowSums(counts(sce) > 1) >= 10, ]
dim(sce)

最后使用logNormCounts，通过将每个计数除以其大小因子、添加伪计数 1 和对数转换，来计算log2转换后的标准表达量

# compute sum-factors & normalize
sce <- computeLibraryFactors(sce)
sce <- logNormCounts(sce)

数据准备

muscat期望输入 SCE 的特定格式。具体而言，必须提供以下单元格元数据（colData）列：

"sample_id"：唯一的样本标识符（例如，PeterPan_ref1、Nautilus_trt3等）
"cluster_id"：亚群（集群）分配（例如，T细胞，单核细胞等）
"group_id"：实验组/条件（例如，对照/治疗，健康/患病等）

sce$id <- paste0(sce$stim, sce$ind)
(sce <- prepSCE(sce, 
    kid = "cell", # subpopulation assignments
    gid = "stim",  # group IDs (ctrl/stim)
    sid = "id",   # sample IDs (ctrl/stim.1234)
    drop = TRUE))  # drop all other colData columns

class: SingleCellExperiment 
dim: 7118 26820 
metadata(1): experiment_info
assays(2): counts logcounts
rownames(7118): NOC2L HES4 ... S100B PRMT2
rowData names(2): ENSEMBL SYMBOL
colnames(26820): AAACATACAATGCC-1 AAACATACATTTCC-1 ... TTTGCATGGTTTGG-1 TTTGCATGTCTTAC-1
colData names(3): cluster_id sample_id group_id
reducedDimNames(1): TSNE
mainExpName: NULL
altExpNames(0):

将聚类和样本 ID 以及聚类和样本的数量存储到以下简单变量中：

nk <- length(kids <- levels(sce$cluster_id))
ns <- length(sids <- levels(sce$sample_id))
names(kids) <- kids; names(sids) <- sids

nk;ns
kids;sids

[1] 8
[1] 16
            B cells     CD14+ Monocytes         CD4 T cells         CD8 T cells     Dendritic cells   FCGR3A+ Monocytes 
          "B cells"   "CD14+ Monocytes"       "CD4 T cells"       "CD8 T cells"   "Dendritic cells" "FCGR3A+ Monocytes" 
     Megakaryocytes            NK cells 
   "Megakaryocytes"          "NK cells" 
   ctrl101   ctrl1015   ctrl1016   ctrl1039    ctrl107   ctrl1244   ctrl1256   ctrl1488    stim101   stim1015   stim1016   stim1039 
 "ctrl101" "ctrl1015" "ctrl1016" "ctrl1039"  "ctrl107" "ctrl1244" "ctrl1256" "ctrl1488"  "stim101" "stim1015" "stim1016" "stim1039" 
   stim107   stim1244   stim1256   stim1488 
 "stim107" "stim1244" "stim1256" "stim1488"

聚类样本大小

# nb. of cells per cluster-sample
t(table(sce$cluster_id, sce$sample_id))

           B cells CD14+ Monocytes CD4 T cells CD8 T cells Dendritic cells FCGR3A+ Monocytes Megakaryocytes NK cells
  ctrl101      113             186         336          95               5                81             12       84
  ctrl1015     476             783         919         226              11               232             25      208
  ctrl1016     144             419         526         671              10               126             15      151
  ctrl1039      30             116         202          30               5                28              5       20
  ctrl107       51             222         197          31               2                29              5       49
  ctrl1244     134             429        1215          82              28                53             19      131
  ctrl1256     240             383        1136         156              12                50             15      275
  ctrl1488     234             317        1343          78              17                99             21      120
  stim101      144             222         437         121              20               126              7      120
  stim1015     357             683         814         153              17               222             24      224
  stim1016     129             361         426         600              10               124             12      239
  stim1039      39             154         318          40               7                36             13       32
  stim107       56             185         217          22               7                36              4       51
  stim1244      94             318         980          46              19                35             14      136
  stim1256     211             369        1047         133              16                73             21      257
  stim1488     283             370        1658          73              34               139             35      187

降维

可以发现直接提供的数据已经包含t-SNE降维坐标信息

# compute UMAP using 1st 20 PCs
sce <- runUMAP(sce, pca = 20)

使用 scat 的 plotReducedDim 函数，我们可以分别绘制由clusters和groups ID 着色的 t-SNE 和 UMAP 表示。我们还创建了一个小的包装函数 .plot_dr()以提高颜色图例的可读性并简化绘图主题：

# wrapper to prettify reduced dimension plots
.plot_dr <- function(sce, dr, col)
  plotReducedDim(sce, dimred = dr, colour_by = col) +
    guides(fill = guide_legend(override.aes = list(alpha = 1, size = 3))) +
    theme_minimal() + theme(aspect.ratio = 1)

对于我们的数据集，由cluster_ids着色的t-SNE和UMAP表明细胞群彼此分离良好。IFN-β刺激表现为当着色group_ids时细胞低维投影的剧烈变化，表明广泛的基因组规模发生转录变化。

# downsample to max. 100 cells per cluster
cs_by_k <- split(colnames(sce), sce$cluster_id)
cs100 <- unlist(sapply(cs_by_k, function(u) 
  sample(u, min(length(u), 100))))

# plot t-SNE & UMAP colored by cluster & group ID
for (dr in c("TSNE", "UMAP"))
  for (col in c("cluster_id", "group_id"))
    .plot_dr(sce[, cs100], dr, col)

差异状态（Differential State，DS）分析

为了测试不同条件下的状态变化，我们将考虑两种类型的方法：

i）直接作用于细胞水平测量的混合模型;
ii）作用于伪批量数据的基于聚合的方法。

1.单细胞数据到pseudo（伪） bulk数据的聚合

为了利用现有的强大的bulk RNA-seq DE框架，如edgeR，DESeq2和limma，我们首先汇总每个样本（在每个分组中）的测量值以获得伪批量数据。

通常，aggregateData() 将通过参数 by 指定的 colData 变量聚合数据，并返回包含伪批量数据的 SingleCellExperiment。

对于 DS 分析，必须在聚类样本级别聚合测量值（默认值by = c("cluster_id", "sample_id"）。

在这种情况下，返回的 SingleCellExperiment 将包含每个簇的一个测定，其中行 = 基因，列 = 样本。

可以发现和bulk RNAseq的输入表达矩阵一致

参数assay和fun分别指定用于聚合的输入数据和汇总统计量。

虽然原则上可以应用输入数据（原始/（对数）标准化计数，CPM等）和汇总统计数据（总和，平均值，中位数）的各种组合，但我们在这里默认为原始计数的总和：

pb <- aggregateData(sce,
    assay = "counts", fun = "sum",
    by = c("cluster_id", "sample_id"))
# one sheet per subpopulation
assayNames(pb)  # cluster

[1] "B cells"           "CD14+ Monocytes"   "CD4 T cells"  
[4] "CD8 T cells"       "Dendritic cells"   "FCGR3A+ Monocytes"
[7] "Megakaryocytes"    "NK cells"

# pseudobulks for 1st subpopulation
t(head(assay(pb)))  # sample

         NOC2L HES4 ISG15 TNFRSF18 TNFRSF4 SDF4
ctrl101     12    0    13        8       1   13
ctrl1015    37    4   136       55      16   42
ctrl1016    13    3    42       10       0   15
ctrl1039     2    1    16        3       1    1
ctrl107      7    0     6        3       4    5
ctrl1244    24    7    25       30      11    7
ctrl1256    22    1    65       30      11   17
ctrl1488    19    1    34       27       6   23
stim101     12    8  1362        3       1    9
stim1015    34   27  4022       21       2   19
stim1016     7    8  1271        2       4    6
stim1039     3    3   393        1       1    2
stim107      4    2   556        1       0    3
stim1244    13   10   830        8       2    5
stim1256    14   10  2235       10       2    6
stim1488    18    5  2927       13       1   19

Pseudobulk-level MDS plot

multi-dimensional scaling (MDS) plot

在进行任何正式测试之前，我们可以计算聚合信号的多维缩放（MDS）图，以探索整体样本的相似性。

pbMDS 将 aggregateData 返回的任何包含 PB 数据的 SCE 作为输入，并使用 edgeR 计算 MDS 维度。理想情况下，数据的这种表示形式应将聚类和组彼此分开。反之亦然，来自同一聚类或组的样本应聚类在一起。

在我们的关于伪散装计数的 MDS 图中可以观察到，第一维（MDS1）清楚地分离了细胞群（簇），而第二维（MDS2）分离了对照和受刺激的样本（组）。此外，两个T细胞簇彼此靠近。

(pb_mds <- pbMDS(pb))

Sample-level analysis: Pseudobulk methods

一旦我们组装了伪批量数据，我们就可以使用 pbDS 测试 DS。默认情况下，a∼group_id模型拟合，并且测试线性模型的最后一个系数等于零。

# run DS analysis
res <- pbDS(pb, verbose = FALSE)
res$table %>% str()
res$table %>% View()

# access results table for 1st comparison
tbl <- res$table[[1]]
# one data.frame per cluster
names(tbl)

[1] "B cells"           "CD14+ Monocytes"   "CD4 T cells"       "CD8 T cells"       "Dendritic cells"   "FCGR3A+ Monocytes"
[7] "Megakaryocytes"    "NK cells"

# view results for 1st cluster
k1 <- tbl[[1]]  # B cells
head(format(k1[, -ncol(k1)], digits = 2))
B_deg_pseudobulk <- format(k1[, -ncol(k1)], digits = 2)

选取Bcells进行展示

pb$group_id

 [1] ctrl ctrl ctrl ctrl ctrl ctrl ctrl ctrl stim stim stim stim stim stim stim stim
Levels: ctrl stim

2.Cell-level analysis: Mixed models

作为上述样本水平方法的替代方案，我们将（针对每个基因）混合模型（MM）拟合到细胞水平的测量数据中。muscat提供了使用 3 种主要方法的 MM 实现：

对数归一化数据与观测权重的线性混合模型（LMM）拟合；
在方差稳定数据上拟合LMM;
直接在计数上拟合广义线性混合模型（GLMM）

在每种情况下，a∼1+group_id+(1|sample_id)模型适合每个基因，优化对数似然（即REML = FALSE，P 值使用参数 df（默认"Satterthwaite") 指定的自由度估计值计算。拟合、测试和适度按亚群维度被应用。对于差异分析，mmDS将仅考虑：

在至少 n_samples 个样本中至少有 n_cells 个细胞（默认为 10 个）的亚群
计数为 >= min_count（默认值 1）的基因至少在 min_cells（默认值 20）中

耗时占用资源过多，未跑完结果

FindMarkers

参考：

任意细胞亚群的差异分析
单细胞差异基因分析：pseudobulk（muscat包）

sce_seurat <- as.Seurat(sce)
sce_seurat
View(sce_seurat)

An object of class Seurat 
7118 features across 26820 samples within 1 assay 
Active assay: originalexp (7118 features, 0 variable features)
 2 dimensional reductions calculated: TSNE, UMAP

Bcells <- subset(sce_seurat, cluster_id=="B cells")
Bcells

An object of class Seurat 
7118 features across 2735 samples within 1 assay 
Active assay: originalexp (7118 features, 0 variable features)
 2 dimensional reductions calculated: TSNE, UMAP

# Idents(sce_seurat) %>% head()
# sce_seurat@meta.data %>% colnames()
# sce_seurat <- SetIdent(sce_seurat,value = "cluster_id")
# Idents(sce_seurat) %>% head()

B_markers_found <- FindMarkers(Bcells,
                             min.pct = 0,logfc.threshold = 0,
                             slot = "counts",
                             group.by = "group_id",ident.1 = "ctrl",ident.2 = "stim")

比较Bcells结果FindMarkers vs pseudoBulk（muscat）

B_deg_pseudobulk$abs_logfc <- B_deg_pseudobulk$logFC %>% as.numeric() %>% abs()
B_markers_found$abs_logfc <- B_markers_found$avg_log2FC %>% as.numeric() %>% abs()

B_deg_pseudobulk$reg <- ifelse(B_deg_pseudobulk$abs_logfc>1&as.numeric(B_deg_pseudobulk$p_val)<0.05,"deg","normal")
table(B_deg_pseudobulk$reg)

B_markers_found$reg <- ifelse(B_markers_found$abs_logfc>1&as.numeric(B_markers_found$p_val)<0.05,"deg","normal")
table(B_markers_found$reg)

   deg normal 
   238   1729
   deg normal 
    17   7101

Pseudobulk鉴定了更多的DEGs

看看FindMakers的结果在不在Pseudobulk的结果中：

table(rownames(B_markers_found)[B_markers_found$reg=="deg"] %in% B_deg_pseudobulk$gene[B_deg_pseudobulk$reg=="deg"])

FALSE  TRUE 
    1    16

MAST

单细胞分析工具--MAST差异基因分析

Using MAST with RNASeq: MAIT Analysis. （http://bioconductor.org/packages/release/bioc/vignettes/MAST/inst/doc/MAITAnalysis.html）

1、示例数据

经细胞因子刺激前后的MAITs细胞的单细胞表达数据，使用MAST包分析其差异表达基因。

data(maits, package='MAST')
dim(maits$expressionmat) #共96个细胞
# [1]    96 16302
maits$expressionmat[1:4,1:4]
head(maits$cdat) #细胞注释数据
head(maits$fdat) #基因注释数据

##整理表达矩阵，使用symbol作为基因名
idx = -which(duplicated(maits$fdat$symbolid))
maits$fdat = maits$fdat[idx,]
maits$fdat = maits$fdat[,-1]
rownames(maits$fdat) = maits$fdat$symbolid
maits$expressionmat = maits$expressionmat[,idx]
colnames(maits$expressionmat) = rownames(maits$fdat)

maits$expressionmat[1:4,1:4]

head(maits$fdat)

2、构建SingleCellAssay对象

library(MAST)
sca = FromMatrix(t(maits$expressionmat), maits$cdat, maits$fdat)
## 方便演示，随机筛选出2000个基因
sca = sca[unique(sample(rownames(sca),2000)),]
dim(sca)
# [1] 2000    96
sca

## (1) 目标分组
colData(sca)
table(colData(sca)$condition)
# Stim Unstim 
# 47     49 
## 将Unstim设置为参考组
cond<-factor(colData(sca)$condition)
cond<-relevel(cond,"Unstim")
cond  # Levels: Unstim Stim
colData(sca)$condition<-cond

## (2) 考虑一个协变量因素
summary(colData(sca)$cngeneson)
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# -0.17159 -0.02578  0.01408  0.00000  0.05943  0.15415

## (3) 模拟一个分类变量分组因素
set.seed(123)
colData(sca)$region = sample(letters[1:3], ncol(sca), replace=T)
table(colData(sca)$region)
# a  b  c 
#32 31 33

3、回归分析

MAST做差异分析的主要优势在于可以校正多种协变量因素，如下演示几种常用的用法。

## (1) 校正cngeneson协变量:默认参数
zlmCond <- zlm(~condition + cngeneson, sca,
    method="bayesglm", ebayes=TRUE)
summaryCond <- summary(zlmCond, doLRT='conditionStim')
summaryDt <- summaryCond$datatable
levels(summaryDt$contrast)
# [1] "conditionStim" "(Intercept)"   "cngeneson"

## (2) 校正cngeneson与region协变量:默认参数
zlmCond <- zlm(~condition + cngeneson + region, sca, 
    method = "bayesglm", ebayes=TRUE)
summaryCond <- summary(zlmCond, doLRT='conditionStim')
summaryDt <- summaryCond$datatable
levels(summaryDt$contrast)
# [1] "conditionStim" "(Intercept)"   "cngeneson"     "regionb"       "regionc"

## (3) 校正cngeneson协变量，将region认为是随机因素
zlmCond <- zlm(~condition + cngeneson + (1|region), sca, 
    method = "glmer", ebayes=FALSE)
summaryCond <- summary(zlmCond, doLRT='conditionStim')
summaryDt <- summaryCond$datatable
levels(summaryDt$contrast)
# [1] "conditionStim" "(Intercept)"   "cngeneson"

汇总差异基因结果

df_pval = summaryDt %>% 
    filter(contrast=='conditionStim') %>% 
    filter(component=='H') %>% 
    select(primerid, `Pr(>Chisq)`)

df_logfc = summaryDt %>% 
    filter(contrast=='conditionStim') %>% 
    filter(component=='logFC') %>% 
    select(primerid, coef, ci.hi, ci.lo)

df_stat = inner_join(df_logfc, df_pval) %>% 
    rename("primerid"="symbol") %>% 
    rename("Pr(>Chisq)"="pval","coef"="logFC") %>% 
    mutate("fdr" = p.adjust(pval)) %>% 
    arrange(fdr)

dim(df_stat)
dim(na.omit(df_stat))
df_stat <- na.omit(df_stat) %>% arrange(fdr) %>% mutate("abs_logfc"=abs(logFC))
head(df_stat)

   symbol     logFC     ci.hi     ci.lo         pval          fdr abs_logfc
1:     U1 -3.980627 -3.248247 -4.713007 1.857517e-16 3.694601e-13  3.980627
2:  PSAT1  7.090388  8.553495  5.627281 8.466767e-16 1.683193e-12  7.090388
3:  DUSP1 -6.513010 -5.260700 -7.765321 8.668938e-15 1.722518e-11  6.513010
4: L26953  4.839456  6.001081  3.677831 1.201518e-11 2.386214e-08  4.839456
5:  PSMA5  4.989426  6.218109  3.760743 3.346791e-11 6.643381e-08  4.989426
6:    JUN -4.171474 -3.213458 -5.129490 4.173966e-11 8.281148e-08  4.171474

FindMarkers

maits$expressionmat[1:4,1:4]
counts <- maits$expressionmat %>% t()
counts[1:4,1:4]
mysce <- CreateSeuratObject(counts = counts[rownames(counts) %in% rownames(sca),])

dim(mysce)
colData(sca)$condition %>% table()
mysce$group_id <- colData(sca)$condition
mysce$group_id %>% table()

[1] 2000   96
.
Unstim   Stim 
    49     47 
.
Unstim   Stim 
    49     47

markers_found <- FindMarkers(mysce,min.pct = 0,logfc.threshold = 0,
                             slot = "counts",group.by = "group_id",
                             ident.1 = "Unstim", ident.2 = "Stim")

|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s

markers_found <- markers_found %>% arrange(p_val) %>% mutate("abs_logfc"=abs(avg_log2FC))
head(markers_found)

             p_val avg_log2FC pct.1 pct.2    p_val_adj abs_logfc
DUSP1 2.478956e-14  2.2429608 0.980 0.766 4.957912e-11 2.2429608
JUND  4.125845e-14  1.1074815 0.959 0.936 8.251689e-11 1.1074815
PSMA5 1.810076e-11 -1.5517307 0.714 1.000 3.620153e-08 1.5517307
RPL34 3.458046e-11  0.2375686 1.000 1.000 6.916093e-08 0.2375686
SMS   6.266779e-11 -1.7052944 0.143 0.787 1.253356e-07 1.7052944
SRP19 1.410429e-10 -1.5377399 0.510 0.957 2.820858e-07 1.5377399

比较结果Findmarkers vs MAST

注意我们只随机取了两千个基因（MAST，MM模型，耗时）

df_stat$reg <- ifelse(df_stat$abs_logfc>1&as.numeric(df_stat$pval)<0.05,"deg","normal")
table(df_stat$reg)

markers_found$reg <- ifelse(markers_found$abs_logfc>1&as.numeric(markers_found$p_val)<0.05,"deg","normal")
table(markers_found$reg)

   deg normal 
   381   1378 

   deg normal 
   149   1851

可以发现MAST鉴定了更多的差异表达基因

table(rownames(markers_found)[markers_found$reg=="deg"] %in% df_stat$symbol[df_stat$reg=="deg"])

FALSE  TRUE 
   23   126

FindMarkers的大部分结果还是在MAST中的

distinct

参考：

distinct: a method for differential analyses via hierarchical permutation tests(https://github.com/SimoneTiberi/distinct) github

https://bioconductor.org/packages/release/bioc/vignettes/distinct/inst/doc/distinct.html

distinct是一种在两组或多组分布之间执行差异分析的统计方法；差异分析是通过对每个样本的累积分布函数（cdfs）进行非参数排列检验来进行的。distinct是一种通用且灵活的工具：由于其完全非参数性质，对数据的生成方式没有任何假设，因此它可以应用于各种数据集。特别适合对单细胞数据进行差异状态分析（即，细胞亚群内的差异分析），例如单细胞RNA测序（scRNA-seq）和高维流式或质谱仪（HDCyto）数据。该方法还考虑了干扰协变量（如批处理效应）。

rm(list=ls())
# BiocManager::install("distinct")
library(Seurat)
library(distinct)

1.加载示例数据

由 6 个样本的子集（在 2 个条件下观察到的 3 个个体）和从 muscData 包的 Kang18_8vs8() 对象中选择的 100 个基因组成

library(SingleCellExperiment)
data("Kang_subset", package = "distinct")
Kang_subset

class: SingleCellExperiment 
dim: 100 9517 
metadata(1): experiment_info
assays(2): logcounts cpm
rownames(100): ISG15 SYF2 ... MX2 PDXK
rowData names(0):
colnames: NULL
colData names(3): stim cell sample_id
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

colData的ind列和stim表示每个细胞的个体ID和实验条件（对照或处理），而sample_id列显示差异分析所需的样品ID。对每个cell聚类分别执行条件之间的差异测试。请注意，如果细胞聚类标签未知，我们需要通过某种聚类算法将细胞聚类到组中。

colData(Kang_subset)

DataFrame with 9517 rows and 3 columns
         stim            cell sample_id
     <factor>        <factor>  <factor>
1        ctrl CD4 T cells     ctrl_107 
2        ctrl CD14+ Monocytes ctrl_1015
3        ctrl NK cells        ctrl_1015
4        ctrl CD4 T cells     ctrl_107 
5        ctrl CD14+ Monocytes ctrl_1015
...       ...             ...       ...
9513     stim CD14+ Monocytes stim_1015
9514     stim CD4 T cells     stim_101 
9515     stim CD14+ Monocytes stim_101 
9516     stim CD14+ Monocytes stim_107 
9517     stim CD4 T cells     stim_107

实验设计比较两组（stim与ctrl），每组有3个生物学重复

Kang_subset@metadata$experiment_info

  sample_id stim
1  ctrl_107 ctrl
2 ctrl_1015 ctrl
3  ctrl_101 ctrl
4  stim_101 stim
5 stim_1015 stim
6  stim_107 stim

2.设计实验

samples = Kang_subset@metadata$experiment_info$sample_id
group = Kang_subset@metadata$experiment_info$stim
design = model.matrix(~group)
# rownames of the design must indicate sample ids:
rownames(design) = samples
design

          (Intercept) groupstim
ctrl_107            1         0
ctrl_1015           1         0
ctrl_101            1         0
stim_101            1         1
stim_1015           1         1
stim_107            1         1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

在条件之间执行差异分析。参数 name_assays_expression 指定assays(x)）中的输入数据（对数counts），而name_cluster和name_sample定义包含细胞聚类（细胞）和单个样本 ID （sample_id）的 colData(x) 的列名。我们要测试的组位于设计的第二列中，因此我们将指定：column_to_test = 2。

请注意，colData(x)$name_sample 中的样本名称必须与 rownames(design) 中的样本名称相同（尽管顺序不一定相同）。

rownames(design)

unique(colData(Kang_subset)$sample_id)

[1] "ctrl_107"  "ctrl_1015" "ctrl_101"  "stim_101"  "stim_1015" "stim_107" 
[1] ctrl_107  ctrl_1015 ctrl_101  stim_101  stim_1015 stim_107 
Levels: ctrl_101 ctrl_1015 ctrl_107 stim_101 stim_1015 stim_107

3.差异分析

为了获得最重要基因的更精细排名，如果计算资源允许，我们鼓励用户增加P_4（即，原始p值<0.001时的排列数量）并设置P_4 = 20，000（默认情况下P_4 = 10,000).

强烈建议使用标准化数据，例如每百万次counts （CPM）或 log2-CPM（例如，通过 scater::logNormCounts 创建的logcounts）).

前面的方法用的都是raw counts

set.seed(61217)

res = distinct_test(x = Kang_subset, 
                    name_assays_expression = "logcounts",
                    name_cluster = "cell",
                    name_sample = "sample_id",
                    design = design,
                    column_to_test = 2,
                    min_non_zero_cells = 20,
                    n_cores = 2)

2 groups of samples provided
Data loaded, starting differential testing
Differential testing completed, returning results

head(res)
table(res$cluster_id)

     gene cluster_id      p_val  p_adj.loc    p_adj.glb filtered  mean_ctrl  mean_stim FC_ctrl/stim log2FC_ctrl/stim
1   ISG15    B cells 0.00009999 0.00139986 0.0007650398    FALSE  198.05692 8577.04339   0.02309151       -5.4364934
2    SYF2    B cells 0.41584158 0.60824590 0.6162697351    FALSE  333.73137  285.50405   1.16891993        0.2251761
3  LAPTM5    B cells 0.01349325 0.10171837 0.0551463399    FALSE 1495.39376 1185.58215   1.26131602        0.3349298
4   EIF3I    B cells 0.91089109 0.94965241 0.9528876576    FALSE  259.08484  233.50382   1.10955291        0.1499785
5 PRPF38A    B cells 0.04590818 0.17303854 0.1379341773    FALSE   61.89912   78.05793   0.79298948       -0.3346264
6   GTF2B    B cells 0.38613861 0.59127475 0.5881463146    FALSE   98.64613  109.72591   0.89902308       -0.1535699

          B cells   CD14+ Monocytes       CD4 T cells       CD8 T cells   Dendritic cells FCGR3A+ Monocytes 
              100               100               100               100               100               100 
   Megakaryocytes          NK cells 
              100               100

结果以data.frame的形式报告，其中gene和cluster_id列包含基因和细胞簇名称，而p_val，p_adj.loc通过p_adj.glbBenjamini和Hochberg（BH）校正报告原始p值，局部和全局调整的p值。在局部调整的 p 值（p_adj.loc）中，BH 校正分别应用于每个聚类，而在全局调整的 p 值（p_adj.glb）中，对所有聚类的结果执行 BH 校正。

我们可以进一步计算组之间的倍数变化（FC）和log2-FC。要计算FCs，请使用标准化数据，例如cpm;不要使用对数转换数据（例如，logcounts）。

res = log2_FC(res = res,
              x = Kang_subset, 
              name_assays_expression = "cpm",
              name_group = "stim",
              name_cluster = "cell")

FC and log2_FC computed, returning results

取出Bcells的deg结果：

bcells_distinct <- res[res$cluster_id=="B cells",]
head(bcells_distinct)

     gene cluster_id      p_val  p_adj.loc    p_adj.glb filtered  mean_ctrl  mean_stim FC_ctrl/stim log2FC_ctrl/stim abs_logfc    reg
1   ISG15    B cells 0.00009999 0.00139986 0.0007650398    FALSE  198.05692 8577.04339   0.02309151       -5.4364934 5.4364934    deg
2    SYF2    B cells 0.41584158 0.60824590 0.6162697351    FALSE  333.73137  285.50405   1.16891993        0.2251761 0.2251761 normal
3  LAPTM5    B cells 0.01349325 0.10171837 0.0551463399    FALSE 1495.39376 1185.58215   1.26131602        0.3349298 0.3349298 normal
4   EIF3I    B cells 0.91089109 0.94965241 0.9528876576    FALSE  259.08484  233.50382   1.10955291        0.1499785 0.1499785 normal
5 PRPF38A    B cells 0.04590818 0.17303854 0.1379341773    FALSE   61.89912   78.05793   0.79298948       -0.3346264 0.3346264 normal
6   GTF2B    B cells 0.38613861 0.59127475 0.5881463146    FALSE   98.64613  109.72591   0.89902308       -0.1535699 0.1535699 normal

此外，distinct还可以处理协变量和批处理效应

FindMarkers

distinct用的logcounts这里我们也用logcounts

Kang_subset
Kang_subset@assays@data$logcounts %>% colnames() %>% head()
counts <- Kang_subset@assays@data$logcounts %>% SparseM::as.matrix()
counts[1:4,1:4]
dim(counts)
table(Kang_subset$stim)
colnames(counts) <- Kang_subset$sample_id
kang_seurat <- CreateSeuratObject(counts = counts)

感觉SingleCellExperiment对象 colData函数取出来的像是Seurat对象的metadata

colData(Kang_subset)
kang_seurat$group_id <- Kang_subset$stim
kang_seurat$celltype <- Kang_subset$cell

取出Bcells

kang_seurat_bcells <- subset(kang_seurat,celltype=="B cells")
table(kang_seurat_bcells$celltype)

          B cells   CD14+ Monocytes       CD4 T cells       CD8 T cells   Dendritic cells FCGR3A+ Monocytes 
             1246                 0                 0                 0                 0                 0 
   Megakaryocytes          NK cells 
                0                 0

kangB_markers_found <- FindMarkers(kang_seurat_bcells,min.pct = 0,logfc.threshold = 0,
                             slot = "counts",group.by = "group_id",
                             ident.1 = "ctrl", ident.2 = "stim")

head(kangB_markers_found)

              p_val avg_log2FC pct.1 pct.2     p_val_adj abs_logfc    reg
ISG15 1.768192e-216 -1.9237559 0.168 0.991 1.768192e-214 1.9237559    deg
MX2    2.722777e-62 -0.6442930 0.038 0.431  2.722777e-60 0.6442930 normal
PSME2  3.444846e-42 -0.5201504 0.497 0.812  3.444846e-40 0.5201504 normal
RPL7   7.601723e-31  0.1705121 0.994 0.969  7.601723e-29 0.1705121 normal
SPI1   2.603640e-16  0.2425778 0.165 0.026  2.603640e-14 0.2425778 normal
PRDX5  1.645329e-14  0.3127005 0.333 0.153  1.645329e-12 0.3127005 normal

结果比较FindMarkers vs distinct

bcells_distinct$abs_logfc <- abs(bcells_distinct$`log2FC_ctrl/stim`)
kangB_markers_found$abs_logfc <- abs(kangB_markers_found$avg_log2FC)

bcells_distinct$reg <- ifelse(bcells_distinct$abs_logfc>1&as.numeric(bcells_distinct$p_val)<0.05,"deg","normal")
table(bcells_distinct$reg)

kangB_markers_found$reg <- ifelse(kangB_markers_found$abs_logfc>1&as.numeric(kangB_markers_found$p_val)<0.05,"deg","normal")
table(kangB_markers_found$reg)

   deg normal 
    16     84 

   deg normal 
     1     99

看一下是不是全部基因数量size的影响

kang_markers_found <- FindMarkers(kang_seurat,min.pct = 0,logfc.threshold = 0,
                             slot = "counts",group.by = "group_id",
                             ident.1 = "ctrl", ident.2 = "stim")
kang_markers_found$abs_logfc <- abs(kang_markers_found$avg_log2FC)
kang_markers_found$reg <- ifelse(kang_markers_found$abs_logfc>1&as.numeric(kang_markers_found$p_val)<0.05,"deg","normal")
table(kang_markers_found$reg)

   deg normal 
     1     99

> kang_markers_found[kang_markers_found$reg=="deg",]
      p_val avg_log2FC pct.1 pct.2 p_val_adj abs_logfc reg
ISG15     0  -2.008996 0.258 0.989         0  2.008996 deg
> kangB_markers_found[kangB_markers_found$reg=="deg",]
              p_val avg_log2FC pct.1 pct.2     p_val_adj abs_logfc reg
ISG15 1.768192e-216  -1.923756 0.168 0.991 1.768192e-214  1.923756 deg

可以发现参与差异分析的基因数量size对结果有影响但不大

"ISG15" %in% bcells_distinct$gene
# TRUE

distinct找到了更多的差异表达基因，结果的格式类似前面的muscat

就在我写这篇推送的时候发现最近我们自己和一些其他的公众号也写了推文介绍一系列做单细胞样本间差异表达分析的方法，可惜我检索整理资料基本上是在8月中旬完成的，没有覆盖到这部分

这里贴出，作为参考：

muscat：专注于多样本多分组的单细胞差异分析

muscat代码实操（一）：模拟复杂的 scRNA-seq 数据以及评估不同的差异表达分析方法的性能

Libra：单细胞差异分析算法的全家桶

开发这个包的灵感来源于之前Seurat和Muscat包（muscat：专注于多样本多分组的单细胞差异分析），但是这两个包都是从不同的单一角度着手，虽各有所长，但适用场景有限。而Libra着手于解决一站式解决scRNA-seq差异分析，其总共集成了22种不同的差异分析算法，这些方法都可以从一个函数中访问，其中包括传统的单细胞方法以及包括pseudobulk和mixed model方法在内的所有主流算法。

clusters间的差异分析鉴定marker

区别于上一部分，我们这里第二部分是“计算来自两个簇的细胞之间的差异表达基因以用于标记基因鉴定或细胞类型注释”，所以可以用官方提供的pbmc单样本数据集

rm(list=ls())
# devtools::install_github('satijalab/seurat-data')
library(SeuratData) #加载seurat数据集  
getOption('timeout')
options(timeout=10000)
InstallData("pbmc3k")
data("pbmc3k")  
sce <- pbmc3k.final   
library(Seurat)
DimPlot(sce,label = T,repel = T)

这里直接获取降维等处理好的数据

> sce
An object of class Seurat 
13714 features across 2638 samples within 1 assay 
Active assay: RNA (13714 features, 2000 variable features)
 2 dimensional reductions calculated: pca, umap

需要理解具体走了哪些处理可以参考初探单细胞下游，同样也是官方的pbmc示例数据

FindAllMarkers

参考：

整合素β1基因敲除前后小鼠肺腺癌上皮细胞水平变化
速度上吊打FindAllMarkers的单细胞亚群特异性高表达基因查询算法

这里开了4个线程，主要是针对 FindAllMarkers 函数

library(future)
# check the current active plan
plan()
plan("multisession", workers = 4)
plan()
start = Sys.time()
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, 
                              min.pct = 0.25, 
                              thresh.use = 0.25)
end = Sys.time() 
dur = end-start 
dur

这里开了4个线程，主要是针对 FindAllMarkers 函数

每个单细胞亚群都有几百个左右的基因被挑选出来，9个单细胞亚群就是2885个基因，当然了，每个基因的特异性并不一样

COSG

参考：

速度上吊打FindAllMarkers的单细胞亚群特异性高表达基因查询算法
单细胞分析工具||COSG鉴定marker基因
【单细胞】上下游，合体！

#remotes::install_github(repo = 'genecell/COSGR')
start = Sys.time()
library(COSG)
marker_cosg <- cosg(
  sce,
  groups='all',
  assay='RNA',
  slot='data',
  mu=1,
  n_genes_user=100)
end = Sys.time() 
dur = end-start 
dur

参数说明如下：

sce是一个SingleCellExperiment对象，其中包含了RNA-seq数据。
groups='all'用于定义要分析的细胞组。在这里，我们使用了数据集中的所有细胞组。
assay='RNA'指定要分析的assay类型为RNA。
slot='data'表示从SingleCellExperiment对象的"data" slot中获取标准化后的数据。
mu=1是cosg函数中的一个参数，用于调整标记基因识别的准确性和灵敏度。默认值为1。
n_genes_user=100是cosg函数中的一个参数，指定在计算标记基因时要考虑的top基因数量。

每个单细胞亚群默认就找到了 n_genes_user参数指定数量单细胞亚群特异性高表达基因：

可以看到，这个结果跟前面的Seurat包的FindAllMarkers函数结果的形式就不一样的，前面的是长数据，现在这个 cosg 函数的结果是宽数据，有点类似于前面我们提到的：基于非负矩阵分解的单细胞降维聚类分群，找各个亚群的特征基因

简单的比较两个函数的结果：

findG = split(sce.markers$gene,sce.markers$cluster)
str(findG)
names(findG)
tmp = do.call(rbind ,
        lapply(names(findG), function(x){
          table(marker_cosg$names[,x] %in% findG[[x]])
        }))
rownames(tmp) = names(findG)
tmp

             FALSE TRUE
Naive CD4 T     14   86
Memory CD4 T    27   73
CD14+ Mono      15   85
B               51   49
CD8 T           46   54
FCGR3A+ Mono    25   75
NK              34   66
DC              54   46
Platelet        37   63

虽然没曾老师推文中那么多，我们的结果中也是，DC和CD8 T 的表现比较差

试试看仅仅是看top10的基因，它们交集如何：

             FALSE TRUE
Naive CD4 T     10   0
Memory CD4 T     5    5
CD14+ Mono       3    7
B                2    8
CD8 T            2    8
FCGR3A+ Mono     3    7
NK               1    9
DC               4    6
Platelet         5    5

这时DC也较差，Platelet的表现也不行，最主要是Naive CD4 T top10竟然一个都没有overlap

和曾老师的推文也有部分出入

翻到评论区：

查看 FindAllmarkers和 cosg 函数的参数：

我一开始以为 FindAllmarkers 函数 min.pct 参数表达该基因的细胞百分比，这个细胞是所有的细胞，经曾老师提醒，这里的细胞是当前差异对照分组细胞类型对应的细胞

添加COSG过滤参数：

#remotes::install_github(repo = 'genecell/COSGR')
start = Sys.time()
library(COSG)
marker_cosg_fit <- cosg(
  sce,
  groups='all',
  assay='RNA',
  slot='data',
  mu=1,
  remove_lowly_expressed = TRUE,
  expressed_pct = 0.25,
  n_genes_user=100)
end = Sys.time() 
dur = end-start 
dur

看看 cosg 函数的结果是否都在前面的 Seurat包的FindAllMarkers函数结果里面：

findG = split(sce.markers$gene,sce.markers$cluster)
str(findG)
names(findG)
tmp = do.call(rbind ,
        lapply(names(findG), function(x){
          table(marker_cosg_fit$names[,x] %in% findG[[x]])
        }))
rownames(tmp) = names(findG)
tmp[,2]

 Naive CD4 T Memory CD4 T   CD14+ Mono            B        CD8 T FCGR3A+ Mono           NK           DC 
         100          100          100          100           99          100          100          100 
    Platelet 
         100

可以发现将过滤阈值调整为一致后，基本上所有的cosg 函数的结果是否都在前面的 Seurat包的FindAllMarkers函数结果中

一些检索细胞亚型markers的数据库

上一部分是自己鉴定不同clusters之间的DEGs检查markers

这里我们也提供一些可以检索的数据资源

https://cellxgene.cziscience.com/datasets

CELL×GENE | Datasets 可能是目前的最大的统一的单细胞数据资源中心，横跨多种单细胞技术，全部的细胞都注释清楚了亚群，以层次结构组织展示

http://117.50.127.228/CellMarker/CellMarker_help.html

细胞类型注释，用它就够：CellMarker 2.0

下面第四部分，我将介绍如何使用一些打分算法结合markers辅助确定细胞亚型，这些方法大部分都用R封装好了，并且有自己使用的细胞亚型的markers
初步定义细胞亚群往往用的是通用规则，我们过去的推文也有很多介绍，如听说你还缺单细胞亚群标记基因

由于篇幅问题，第四部分“一些辅助划分细胞亚型的工具”将挪到本期推送的次版里

ixxmu / mp_duty

单细胞DEGs、markers、celltypes | 专辑完结篇 #3871

单细胞DEGs、markers、celltypes | 专辑完结篇 by 生信菜鸟团

单细胞水平样本间的差异分析

background

muscat

加载数据

预处理

数据准备

聚类样本大小

降维

差异状态（Differential State，DS）分析

1.单细胞数据到pseudo（伪） bulk数据的聚合

Pseudobulk-level MDS plot

Sample-level analysis: Pseudobulk methods

2.Cell-level analysis: Mixed models

FindMarkers

比较Bcells结果FindMarkers vs pseudoBulk（muscat）

MAST

FindMarkers

比较结果Findmarkers vs MAST

distinct

FindMarkers

结果比较FindMarkers vs distinct

clusters间的差异分析鉴定marker

FindAllMarkers

COSG

简单的比较两个函数的结果：

一些检索细胞亚型markers的数据库