专栏首页单细胞天地Seurat包的findmarkers函数只能根据划分好的亚群进行差异分析吗

Seurat包的findmarkers函数只能根据划分好的亚群进行差异分析吗

虽然是免费的, 但是学员关于课程的提问,我还是会尽我所能帮忙的。当然,受限于时间和精力,只能是挑选重点和普适性的有价值的问题,比如其中一个问题是:

我发现我的目的基因主要分布在B细胞群,想问一下能不能通过是否表达目的基因将B细胞分成两群,再寻找差异表达基因呀?我看seurat包中,findmarkers的函数只要能找不同cluster 间的差异基因?

这个问题有两个解决方案,第一个把已经划分为B细胞群的那些细胞的表达矩阵,重新走seurat流程,看看这个时候它们是否是否根据有没有表达目的基因来进行分群,如果有,就可以使用 findmarkers 函数。

另外一个解决方案,就是人为划分,然后强行走seurat标准流程。

首先看看 seurat标准流程

首先看看我们的seurat标准流程,自己去GEO数据库下载GSE125616的GSE125616_READ_COUNT_table_lvbo_TE_gencode.txt.gz文件,然后走下面的代码即可:

rm(list=ls())
options(stringsAsFactors = F)
library(data.table)
a=fread('GSE125616_READ_COUNT_table_lvbo_TE_gencode.txt.gz',
        header = T,sep = '\t',data.table = F)
hg=a$Geneid
dat=a[,7:ncol(a)]
rownames(dat)=hg
hg[grepl('^MT-',hg)]
head(hg)
library(stringr)
colnames(dat)
colnames(dat)=str_split(colnames(dat),'_',simplify = T)[,1]
meta=as.data.frame(str_split(colnames(dat),'_',simplify = T)[,1])

colnames(meta)=c('cell name') 
rownames(meta)=colnames(dat)
library(Seurat)
pbmc  <- CreateSeuratObject(counts = dat, 
                            meta.data = meta ,
                            min.cells = 5, min.features = 1000, 
                            project = "test")
pbmc
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) 
# 步骤 ScaleData 的耗时取决于电脑系统配置(保守估计大于一分钟)
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)
DimPlot(sce,reduction = "tsne",label=T)
phe=data.frame(cell=rownames(sce@meta.data),
               cluster =sce@meta.data$seurat_clusters)
head(phe)
table(phe$cluster)

可以看到, 走默认代码,划分的是5个亚群:

这个时候,学员的问题可以简化为,不满意上图里面的0那个亚群,觉得它可以继续划分,比如根据感兴趣的目的基因,强行把0那个红色亚群继续划分为高低表达(或者是否表达)目的基因的亚群。

先熟悉FindMarkers函数

通常,我们使用FindMarkers函数针对感兴趣的细胞亚群,去寻找它与其它所有的亚群,表达有差异的基因,代码如下:

markers_df <- FindMarkers(object = sce, ident.1 = 0, min.pct = 0.25)
print(x = head(markers_df))
markers_genes =  rownames(head(x = markers_df, n = 5))
VlnPlot(object = sce, features =markers_genes,log =T )
FeaturePlot(object = sce, features=markers_genes )

markers_df <- FindMarkers(object = sce, ident.1 = 1, min.pct = 0.25)
print(x = head(markers_df))
markers_genes =  rownames(head(x = markers_df, n = 5))
VlnPlot(object = sce, features =markers_genes,log =T )
FeaturePlot(object = sce, features=markers_genes )

如果你仔细看这个FindMarkers函数的说明书,它其实是可以给定两个亚群,然后单独比较这两个亚群的。

## Default S3 method:
FindMarkers(
  object,
  slot = "data",
  counts = numeric(),
  cells.1 = NULL,
  cells.2 = NULL,
  features = NULL,
  reduction = NULL,
  logfc.threshold = 0.25,
  test.use = "wilcox",
  min.pct = 0.1,
  min.diff.pct = -Inf,
  verbose = TRUE,
  only.pos = FALSE,
  max.cells.per.ident = Inf,
  random.seed = 1,
  latent.vars = NULL,
  min.cells.feature = 3,
  min.cells.group = 3,
  pseudocount.use = 1,
  ...
)

所以,我们的目标就是把两个亚群信息,给到这个FindMarkers函数即可。当然了,这个差异分析结果表格也是需要理解的:

avg_logFC: log fold-chage of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group
pct.1: The percentage of cells where the gene is detected in the first group
pct.2: The percentage of cells where the gene is detected in the second group
p_val_adj: Adjusted p-value, based on bonferroni correction using all genes in the dataset

主要是借助统计可视化来。

根据高低表达(或者是否表达)目的基因划分亚群

其实在这个FindMarkers函数的说明书里面,就有一个现成的例子:

# Take all cells in cluster 2, and find markers that separate cells in the 'g1' group (metadata
# variable 'group')
markers <- FindMarkers(pbmc_small, ident.1 = "g1", group.by = 'groups', subset.ident = "2")
head(x = markers)

所以我们也是需要给单细胞数据对象,制作一个新的分组即可;

# 下面的 TP53 > 10 仅仅是我举例方便,你有自己感兴趣的基因,自己感兴趣的基因表达量阈值
highCells=colnames(subset(x = sce, subset = TP53 > 10, slot = 'counts'))
highORlow=ifelse(colnames(sce) %in% highCells,'high','low')
table(highORlow)
sce@meta.data$highORlow=highORlow
markers <- FindMarkers(sce, ident.1 = "high", 
                       group.by = 'highORlow', 
                       subset.ident = "0")
head(x = markers)

是不是很简单啊?

为什么我会挑选这个问题来解答呢,其实就是督促大家读说明书。

本文分享自微信公众号 - 单细胞天地(sc-ngs),作者:生信技能树

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-04-27

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 用Seurat包分析文章数据(二)

    下面这个图第一张和第三张就明显是负相关:也说明了检测的ERCC占比越多,检测到的基因数目就越少(因为第二张图中总共的reads数是一定的嘛,大约就0.5M)

    生信技能树jimmy
  • 用Scater包分析文章数据

    Scater需要利用SingleCellExperiment这个对象(需要注意的是,官方友情提示,在导入对象之前,最好是将表达量数据存为矩阵)

    生信技能树jimmy
  • 使用scran包的MNN算法来去除多个单细胞转录组数据批次效应

    这里我们使用一下scran包的 mutual nearest neighbors (MNNs)方法吧,主要就是读文档而已:https://bioconducto...

    生信技能树jimmy
  • 用Seurat包分析文章数据(二)

    下面这个图第一张和第三张就明显是负相关:也说明了检测的ERCC占比越多,检测到的基因数目就越少(因为第二张图中总共的reads数是一定的嘛,大约就0.5M)

    生信技能树jimmy
  • AngularJS 使用$sce控制代码安全检查

    由于浏览器都有同源加载策略,不能加载不同域下的文件、也不能使用不合要求的协议比如file进行访问。 在angularJs中为了避免安全漏洞,一些ng-src...

    用户1154259
  • 使用scran包的MNN算法来去除多个单细胞转录组数据批次效应

    这里我们使用一下scran包的 mutual nearest neighbors (MNNs)方法吧,主要就是读文档而已:https://bioconducto...

    生信技能树jimmy
  • 用Scater包分析文章数据

    Scater需要利用SingleCellExperiment这个对象(需要注意的是,官方友情提示,在导入对象之前,最好是将表达量数据存为矩阵)

    生信技能树jimmy
  • 如果传统bulk转录组数据队列足够大也可以使用单细胞流程

    还给出了一些简单代码,就是看看样本聚类情况,然后留成作业给另外一个学徒,看单细胞R包Seurat的FindAllMarkers函数对7个亚型找到的marker基...

    生信技能树
  • 使用seurat3里面计算线粒体基因含量的2个方法

    展示的非常清楚啦,因为每个教程想说明的情况不一样,所以需要重新把计算线粒体基因含量讲解一下。 为了维持教程的统一性,我这里一直使用 sce 来代表构建好的seu...

    生信技能树jimmy
  • seurat3的merge功能和cellranger的aggr整合多个10X单细胞转录组对比

    我们得比较一下,作者的ellranger的aggr整合多个10X单细胞转录组得到的表达矩阵,跟我们使用seurat3的merge功能整合8个10X单细胞转录组样...

    生信技能树jimmy

扫码关注云+社区

领取腾讯云代金券