FindAllMarkers, FindMarkers 以及 FindConservedMarkers 的区别

在seurat中,如果运行了RunUMAP或者RunTSNE后自动分群后,FindAllMarkers和FindMarkers基本就是一样的;如果没有进行RunUMAP或者RunTSNE分群,那么需要先运行BuildClusterTree(object)函数,利用树聚类先分群

FindAllMarkers and FindMarkers

示例:
首先加载数据

pbmc.data <- read.csv("GSE117988_raw.expMatrix_PBMC.csv",header=TRUE,row.names = 1)
pbmc.data <- log2(1 + sweep(pbmc.data, 2, median(colSums(pbmc.data))/colSums(pbmc.data), '*'))
pbmc<- CreateSeuratObject(counts = pbmc.data, project = "10X pbmc", min.cells = 1, min.features = 0)

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

sce=pbmc
sce <- NormalizeData(sce, normalization.method =  "LogNormalize", scale.factor = 10000)
GetAssay(sce,assay = "RNA")
sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000) 

sce <- ScaleData(sce) 
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce)) 
sce <- FindNeighbors(sce, dims = 1:15)
sce <- FindClusters(sce, resolution = 0.2)
table(sce@meta.data$RNA_snn_res.0.2) 

set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)
## 运行 FindAllMarkers
all.markers <- FindAllMarkers(object = sce)

先看FindAllMarkers:

object = sce
assay = NULL
features = NULL
logfc.threshold = 0.25
test.use = 'wilcox'
slot = 'data'
min.pct = 0.1
min.diff.pct = -Inf
node = NULL
verbose = TRUE
only.pos = FALSE
max.cells.per.ident = Inf
random.seed = 1
latent.vars = NULL
min.cells.feature = 3
min.cells.group = 3
mean.fxn = NULL
fc.name = NULL
base = 2
return.thresh = 1e-2
densify = FALSE

if (is.null(x = node)) {
  ## 将分群的结果进行排序
  idents.all <- sort(x = unique(x = Idents(object = object)))
  ## 一共 7 个类群
  ##[1] 0 1 2 3 4 5 6
  ##Levels: 0 1 2 3 4 5 6
} 

genes.de <- list()
messages <- list()

## 分别遍历每一个细胞类群( 0 1 2 3 4 5 6 )
for (i in 1:length(x = idents.all)) {
  if (verbose) {
    message("Calculating cluster ", idents.all[i])
  }
  genes.de[[i]] <- tryCatch(
    expr = {
      FindMarkers(
        object = object,
        assay = assay,
        ident.1 = if (is.null(x = node)) {
          idents.all[i]
        } else {
          tree
        },
        ident.2 = if (is.null(x = node)) {
          NULL
        } else {
          idents.all[i]
        },
        features = features,
        logfc.threshold = logfc.threshold,
        test.use = test.use,
        slot = slot,
        min.pct = min.pct,
        min.diff.pct = min.diff.pct,
        verbose = verbose,
        only.pos = only.pos,
        max.cells.per.ident = max.cells.per.ident,
        random.seed = random.seed,
        latent.vars = latent.vars,
        min.cells.feature = min.cells.feature,
        min.cells.group = min.cells.group,
        mean.fxn = mean.fxn,
        fc.name = fc.name,
        base = base,
        densify = densify
      )
    },
    error = function(cond) {
      return(cond$message)
    }
  )
}
gde.all <- data.frame()

# 对每一类群找到的 marker 基因进行阈值过滤
for (i in 1:length(x = idents.all)) {
  gde <- genes.de[[i]]
  gde <- gde[order(gde$p_val, -gde[, 2]), ]
  gde <- subset(x = gde, subset = p_val < return.thresh)
  gde$cluster <- idents.all[i]
  gde$gene <- rownames(x = gde)
  gde.all <- rbind(gde.all, gde)

# 重命名  
rownames(x = gde.all) <- make.unique(names = as.character(x = gde.all$gene))

我们可以知道,当运行完RunUMAP或者RunTSNE后自动分群后,FindAllMarkers 内部其实就是调用了 FindMarkers,并且:

# FindMarkers 的 ident.1 为RunUMAP或者RunTSNE后自动分群后每一个类群的细胞
ident.1 = if (is.null(x = node)) 
{
    idents.all[i]
} 

# FindMarkers 的 ident.2 默认为空 
ident.2 = if (is.null(x = node))
{
     NULL
} 

而 FindMarkers 在计算差异基因的时候需要先运行 IdentsToCells() 函数

IdentsToCells <- function(
  object,
  ident.1,
  ident.2,
  cellnames.use
) {
  #
  if (is.null(x = ident.1)) {
    stop("Please provide ident.1")
  } else if ((length(x = ident.1) == 1 && ident.1[1] == 'clustertree') || is(object = ident.1, class2 = 'phylo')) {
    if (is.null(x = ident.2)) {
      stop("Please pass a node to 'ident.2' to run FindMarkers on a tree")
    }
    tree <- if (is(object = ident.1, class2 = 'phylo')) {
      ident.1
    } else {
      Tool(object = object, slot = 'BuildClusterTree')
    }
    if (is.null(x = tree)) {
      stop("Please run 'BuildClusterTree' or pass an object of class 'phylo' as 'ident.1'")
    }
    ident.1 <- tree$tip.label[GetLeftDescendants(tree = tree, node = ident.2)]
    ident.2 <- tree$tip.label[GetRightDescendants(tree = tree, node = ident.2)]
  }
  if (length(x = as.vector(x = ident.1)) > 1 &&
      any(as.character(x = ident.1) %in% cellnames.use)) {
    bad.cells <- cellnames.use[which(x = !as.character(x = ident.1) %in% cellnames.use)]
    if (length(x = bad.cells) > 0) {
      stop(paste0("The following cell names provided to ident.1 are not present in the object: ", paste(bad.cells, collapse = ", ")))
    }
  } else {
    ident.1 <- WhichCells(object = object, idents = ident.1)
  }
  # if NULL for ident.2, use all other cells
  if (length(x = as.vector(x = ident.2)) > 1 &&
      any(as.character(x = ident.2) %in% cellnames.use)) {
    bad.cells <- cellnames.use[which(!as.character(x = ident.2) %in% cellnames.use)]
    if (length(x = bad.cells) > 0) {
      stop(paste0("The following cell names provided to ident.2 are not present in the object: ", paste(bad.cells, collapse = ", ")))
    }
  } else {
    if (is.null(x = ident.2)) {
      ident.2 <- setdiff(x = cellnames.use, y = ident.1)
    } else {
      ident.2 <- WhichCells(object = object, idents = ident.2)
    }
  }
  return(list(cells.1 = ident.1, cells.2 = ident.2))
}

对于 FindMarkers 若 ident.2 = NULL ,那么 IdentsToCells() 函数将会选择除 ident.1 类群以外的剩下所有细胞

# select which data to use
assay <- assay %||% DefaultAssay(object = sce)
data.use <- object[[assay]]
# cellnames.use 代表所有的细胞
cellnames.use <-  colnames(x = data.use)

 # if NULL for ident.2, use all other cells
 if (is.null(x = ident.2)) {
      ## 选择除 ident.1 类群以外的其他所有细胞,setdiff 为取差集
      ident.2 <- setdiff(x = cellnames.use, y = ident.1)
 } else {
      ident.2 <- WhichCells(object = object, idents = ident.2)
   }
}

因此,对于 FindMarkers 若 ident.2 = NULL ,计算的是 ident.1 细胞类群与除 ident.1 类群以外的剩下所有细胞相比的差异基因;如果 ident.2 = 某一细胞类群 那么计算的是 ident.1 细胞类群与ident.2 类群相比的差异基因

FindConservedMarkers

# 安装必要的包
BiocManager::install('multtest')
install.packages('metap')
install.packages('sn')
library(Seurat)

# 加载示例数据
data("pbmc_small")
# 假设所有细胞一共分为 'g1' 和 'g2' 两种细胞类群
pbmc_small[['groups']] <- sample(x = c('g1', 'g2'), size = ncol(x = pbmc_small), replace = TRUE)
FindConservedMarkers(pbmc_small, ident.1 = 0, ident.2 = 1, grouping.var = "groups")

pbmc_small[['groups']]
> pbmc_small[['groups']]
               groups
ATGCCAGAACGACT     g1
CATGGCCTGTGCAT     g2
GAACCTGATGAACC     g1
TGACTGGATTCTCA     g1
AGTCAGACTGCACA     g2
TCTGATACACGTGT     g1
TGGTATCTAAACAG     g1
GCAGCTCTGTTTCT     g1
GATATAACACGCAT     g2
AATGTTGACAGTCA     g2


### FindConservedMarkers 源代码分析
object = pbmc_small
ident.1 = 0
ident.2 = 1
grouping.var = "groups"
assay = 'RNA'
slot = 'data'
min.cells.group = 3
meta.method = metap::minimump
verbose = TRUE

# 假设有两个分组, 分别为 'g1' 和 'g2' 
# ident.1 和 ident.2 为每一组('g1' 和 'g2') 的细胞分群, ident.1 = 0, ident.2 = 1
# 提取每个细胞对应分组的信息,即每个细胞分为 'g1' 还是 'g2'
object.var <- FetchData(object = object, vars = grouping.var)
# 提取 'g1' 和 'g2' 的细胞
object <- SetIdent(
  object = object,
  cells = colnames(x = object),
  value = paste(Idents(object = object), object.var[, 1], sep = "_")
)
levels.split <- names(x = sort(x = table(object.var[, 1])))
num.groups <- length(levels.split)
cells <- list()

# 分别将 'g1' 和 'g2' 的细胞单独取出来,  cells[[1]] 为 'g1' 的细胞, cells[[2]] 为 'g2' 的细胞
for (i in 1:num.groups) {
  cells[[i]] <- rownames(
    x = object.var[object.var[, 1] == levels.split[i], , drop = FALSE]
  )
}
marker.test <- list()
# do marker tests
ident.2.save <- ident.2
# 分别遍历 'g1' 和 'g2' 
for (i in 1:num.groups) {
  level.use <- levels.split[i]
  ident.use.1 <- paste(ident.1, level.use, sep = "_")
  ident.use.1.exists <- ident.use.1 %in% Idents(object = object)

  ident.2 <- ident.2.save
  cells.1 <- WhichCells(object = object, idents = ident.use.1)
  ident.use.2 <- paste(ident.2, level.use, sep = "_")
  ident.use.2.exists <- ident.use.2 %in% Idents(object = object)
  
  ## ident.use.1 = "0_g1" or "0_g2"; ident.use.2 = "1_g1" or "1_g2"
  ## 分别对 'g1' 计算类群 0 与类群 1 的差异基因; 对 'g2' 计算类群 0 与类群 1 的差异基因
  marker.test[[i]] <- FindMarkers(
    object = object,
    assay = assay,
    slot = slot,
    ident.1 = ident.use.1,
    ident.2 = ident.use.2,
    verbose = verbose
  )
  names(x = marker.test)[i] <- levels.split[i]
}
## 过滤 NULL 值
marker.test <- Filter(f = Negate(f = is.null), x = marker.test)
## 对 'g1' 中类群 0 与类群 1 差异基因item和 'g2' 中类群 0 与类群 1 差异基因item取交集
genes.conserved <- Reduce(
  f = intersect,
  x = lapply(
    X = marker.test,
    FUN = function(x) {
      return(rownames(x = x))
    }
  )
)

## 将剩下的 logFC, p_val 整理进去
markers.conserved <- list()
for (i in 1:length(x = marker.test)) {
  markers.conserved[[i]] <- marker.test[[i]][genes.conserved, ]
  colnames(x = markers.conserved[[i]]) <- paste(
    names(x = marker.test)[i],
    colnames(x = markers.conserved[[i]]),
    sep = "_"
  )
}

## 合并对 'g1' 中类群 0 与类群 1 差异基因item和 'g2' 中类群 0 与类群 1 差异基因item取交集后的表格
## 包括 logFC, p_val
markers.combined <- Reduce(cbind, markers.conserved)
pval.codes <- colnames(x = markers.combined)[grepl(pattern = "*_p_val$", x = colnames(x = markers.combined))]
if (length(x = pval.codes) > 1) {
  markers.combined$max_pval <- apply(
    X = markers.combined[, pval.codes, drop = FALSE],
    MARGIN = 1,
    FUN = max
  )
  combined.pval <- data.frame(cp = apply(
    X = markers.combined[, pval.codes, drop = FALSE],
    MARGIN = 1,
    FUN = function(x) {
      return(meta.method(x)$p)
    }
  ))
  meta.method.name <- as.character(x = formals()$meta.method)
  colnames(x = combined.pval) <- paste0(meta.method.name, "_p_val")
  markers.combined <- cbind(markers.combined, combined.pval)
  markers.combined <- markers.combined[order(markers.combined[, paste0(meta.method.name, "_p_val")]), ]
}

所谓的 FindConservedMarkers 其实是在单细胞分析中,有两个处理组 'g1' 和 'g2' ,而 'g1' 和 'g2' 组中又恰好有相同的细胞亚群 0 和 1,求在 'g1' 中细胞亚群 0 和细胞亚群 1 的差异基因,并且与在 'g2' 中细胞亚群 0 和细胞亚群 1 的差异基因取交集;这部分基因定义为 在 'g1' 和 'g2' 组中,细胞亚群 0 和细胞亚群 1 差异保守的基因

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,547评论 6 477
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,399评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,428评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,599评论 1 274
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,612评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,577评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,941评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,603评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,852评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,605评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,693评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,375评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,955评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,936评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,172评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 43,970评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,414评论 2 342

推荐阅读更多精彩内容