专栏首页单细胞天地条条道路通罗马—单细胞分群分析

条条道路通罗马—单细胞分群分析

这一次的问题是:分析单细胞转录组一定要用R包吗?

之前在 差异分析及可视化 中使用monocle的plot_cell_clusters函数画出了PBMC的第4和第10群两种不同T细胞的差异。那么这个分析一定要用包装好的R包吗?不是的,即使不使用别人做的R包,自己也能利用作图原理画出来

问题的重点不在使用什么包,而是做出这个图的原理是什么

就上面这个图而言,为了分群,首先要有一个表达矩阵,还要设置分群信息(例如上面我们就明确知道要分成两群),然后还要进行PCA的线性降维和tSNE进一步非线性降维。有了这些认识,那么首先使用普通的函数来尝试一下

第一种方式:不使用R包,使用常规函数

先加载表达矩阵和分群信息
options(warn=-1) 
load('patient1.PBMC_RespD376_for_DEG.Rdata')
# 其中保存了表达矩阵(count_matrix)和分群的信息(cluster)
> count_matrix[1:4,1:4]
              AAACCTGCAACGATGG.3 AAACCTGGTCTCCATC.3 AAACGGGAGCTCCTTC.3 AAAGCAATCATATCGG.3
RP11-34P13.7                   0                  0                  0                  0
FO538757.2                     0                  0                  0                  0
AP006222.2                     0                  0                  0                  0
RP4-669L17.10                  0                  0                  0                  0
> table(cluster)
cluster
  4  10 
636 170 
根据分群因子,赋予不同颜色
# 根据分群因子,赋予不同颜色
color <- cluster
levels(color) <- rainbow(2)
> table(color)
color
#FF0000FF #00FFFFFF 
      636       170 
对表达矩阵过滤
# 首先是对基因过滤,根据标准差将那些在细胞中的表达量没有变化的基因去掉
choosed_count <- count_matrix
choosed_count <- choosed_count[apply(choosed_count, 1, sd)>0,]
# 看到过滤了5000多个基因
# 然后选择变化差异最明显的前1000个基因
choosed_count <- choosed_count[names(head(sort(apply(choosed_count, 1, sd),decreasing = T),1000)),]
然后进行PCA分析

对行进行操作,因此需要转置,另外进行标准化

pca_out <- prcomp(t(choosed_count),scale. = T)
library(ggfortify)
autoplot(pca_out, col=color) +theme_classic()+ggtitle('PCA plot')

PCA的全部结果可以通过str(pca_out)了解,其中坐标的信息在pca_out$x

>  pca_out$x[1:3,1:3]
                          PC1      PC2        PC3
AAACCTGCAACGATGG.3 -1.4940046 6.707691 -0.7655250
AAACCTGGTCTCCATC.3  0.1403445 3.239395  0.2672606
AAACGGGAGCTCCTTC.3 -1.9968613 3.924992 -0.4669826
再用tSNE继续降维
library(Rtsne)
# 利用PCA的前5个主成分
tsne_out <- Rtsne(pca_out$x[,1:6], perplexity = 10,
                    pca = F, max_iter = 2000,
                    verbose = T)
# tSNE的坐标结果保存在tsne_out$Y中
tsnes_cord <- tsne_out$Y
# ggplot作图
colnames(tsnes_cord) <- c('tSNE1','tSNE2')
ggplot(tsnes_cord, aes(x=tSNE1, y = tSNE2)) + geom_point(col=color) + theme_classic()+ggtitle('tSNE plot')

第二种方式:使用Seurat

Seurat V2
PBMC <- CreateSeuratObject(count_matrix, 
                           min.cells = 1, min.features = 0, 
                           project = '10x_PBMC')
PBMC <- AddMetaData(object = PBMC, metadata = apply(count_matrix, 2, sum), col.name = 'nUMI_raw')
PBMC <- AddMetaData(object = PBMC, metadata = cluster, col.name = 'cluster')
VlnPlot(PBMC, features.plot  = c("nUMI_raw", "nGene"), 
        group.by = 'cluster',nCol = 2)

image-20191018170939612

看到第10群细胞的基因数和UMI数比第4群多很多,也就是说它们本身的“起跑线”就有差异,那么最后的PCA结果差异难道是由于基因和UMI数量导致的?

那么如果我们接下来进行校正,会发生什么?

PBMC <- ScaleData(PBMC, vars.to.regress =  c("nUMI_raw", "nGene"),
                  model.use = 'linear',
                  use.umi = F)
# 这个FindVariableGenes其实和前面的找top1000 sd基因目的一样,就是为了增加分析的有效性,认为高变化的基因存储的信息更值得关注
PBMC <- FindVariableGenes(object = PBMC, mean.function = ExpMean, 
                          dispersion.function = LogVMR, 
                          x.low.cutoff = 0.0125, x.high.cutoff = 3, 
                          y.cutoff = 0.5)
length(PBMC@var.genes)
PBMC <- RunPCA(object = PBMC, pc.genes = PBMC@var.genes)
PBMC <- RunTSNE(object = PBMC, dims.use = 1:20)
TSNEPlot(PBMC, group.by='cluster')

这次结果发现,这两群就分不开了。从这里也反映出一些问题:本文的这个热图真的是由于生物学因素导致的吗?

猜想:可能这两群细胞本身表达的基因数量就不同,就是有一些基因在这群细胞表达,在那一群不表达

i

Seurat V3
PBMC_V3 <- CreateSeuratObject(counts = count_matrix,
                                 min.cells = 1, 
                                 min.features = 0, 
                                 project = "10x_PBMC_V3")
PBMC_V3 <- AddMetaData(object = PBMC_V3, metadata = apply(count_matrix, 2, sum), col.name = 'nUMI_raw')
PBMC_V3 <- AddMetaData(object = PBMC_V3, metadata = cluster, col.name = 'cluster')
VlnPlot(PBMC_V3, features =  c("nUMI_raw", "nFeature_RNA"), ncol = 2,group.by = 'cluster')
# 同样也进行校正一下
PBMC_V3 <- ScaleData(object = PBMC_V3, vars.to.regress =  c("nUMI_raw", "nFeature_RNA"), 
                     model.use = 'linear', use.umi = FALSE)

PBMC_V3 <- FindVariableFeatures(object = PBMC_V3, mean.function = ExpMean, 
                                dispersion.function = LogVMR, 
                                mean.cutoff = c(0.0125,3), 
                                dispersion.cutoff = c(0.5,Inf))
length(VariableFeatures(object = PBMC_V3))
PBMC_V3 <- RunPCA(PBMC_V3, features = VariableFeatures(object = PBMC_V3))
pbmc_tsne <- RunTSNE(PBMC_V3, dims = 1:20)
DimPlot(pbmc_tsne, reduction = "tsne",group.by='cluster')

第三种方式:使用monocle

library(monocle) 
# 1.表达矩阵
expr_matrix <- as.matrix(count_matrix)
# 2.细胞信息
sample_ann <- data.frame(cells=names(count_matrix),  
                         cellType=cluster)
rownames(sample_ann)<- names(count_matrix)
# 3.基因信息
gene_ann <- as.data.frame(rownames(count_matrix))
rownames(gene_ann)<- rownames(count_matrix)
colnames(gene_ann)<- "genes"
# 然后转换为AnnotatedDataFrame对象
pd <- new("AnnotatedDataFrame",
          data=sample_ann)
fd <- new("AnnotatedDataFrame",
          data=gene_ann)
# 最后构建CDS对象
sc_cds <- newCellDataSet(
  expr_matrix, 
  phenoData = pd,
  featureData =fd,
  expressionFamily = negbinomial.size(),
  lowerDetectionLimit=0.5)

cds=sc_cds
cds <- detectGenes(cds, min_expr = 1)
expressed_genes <- row.names(subset(cds@featureData@data,
                                    num_cells_expressed >= 1))
length(expressed_genes)
cds <- cds[expressed_genes,]

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
disp_table <- dispersionTable(cds) # 挑有差异的
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1) # 挑表达量不太低的
cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id)  # 准备聚类基因名单
# 进行降维
cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
                       reduction_method = 'tSNE', verbose = T)
# 进行聚类
cds <- clusterCells(cds, num_clusters = 4) 
# Distance cutoff calculated to 1.812595  
plot_cell_clusters(cds, 1, 2, color = "cellType")

第四种方式:使用scater

suppressMessages(library(scater))
## 创建 scater 要求的对象
# 在导入对象之前,最好是将表达量数据存为矩阵
sce <- SingleCellExperiment(
  assays = list(counts = as.matrix(count_matrix)), 
  colData = data.frame(cluster)
)
# 预处理
exprs(sce) <- log2(calculateCPM(sce ) + 1)
# 降维
sce <- runPCA(sce)
plotReducedDim(sce, use_dimred = "PCA", 
               colour_by= "cluster")
set.seed(1000)
sce <- runTSNE(sce, perplexity=10)
plotTSNE(sce, 
         colour_by= "cluster")

本文分享自微信公众号 - 单细胞天地(sc-ngs),作者:刘小泽

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

原始发表时间:2019-11-05

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 细胞亚群的生物学命名

    也就是把四个时间点映射到上面的tsne坐标中,并且理论上应该是:每群细胞都覆盖到四个时间点

    生信技能树jimmy
  • 差异分析及功能注释(下)

    因为有一些de_genes的SYMBOL没有对应的ENTREZID,因此看到少了100多个基因。

    生信技能树jimmy
  • 单细胞转录组技术梳理(一)

    谈起单细胞转录组测序就不得不提到北京大学汤富酬教授,2009年,汤富酬老师在博士后期间发表了世界上第一篇单细胞mRNA测序的文章“mRNA-Seq whole-...

    生信技能树jimmy
  • 直播系统源码需要不断更新迭代与升级

    直播行业发展至今仍处于高速成长期,从一开始的秀场直播到现在的游戏、教育、财经、音乐等多品类的直播类型,覆盖了各个行业和领域。随着科技不断进步,人们对直播的需求也...

    布谷安妮
  • Robin谈早期点石博客的优化策略

    有人很疑惑的问我,你们4个都是研究SEO的,为什么不在点石博客实施一些搜索引擎优化策略呢?如果大家只是把网页源码META标签部分的keywords和descri...

    蝙蝠侠IT
  • Android Studio使用ButterKnife和Zelezny的方法

    ButterKnife是一个专注于Android的View注入框架,可以减少大量的findViewById以及setOnClickListener代码,可视化一...

    砸漏
  • 【Java小工匠】JavaNIO-缓存区基础

      缓冲区(Buffer)就是在内存中预留指定大小的存储空间用来对输入/输出(I/O)的数据作临时存储,这部分预留的内存空间就叫做缓冲区。

    Java小工匠
  • Kubernetes Network Policy 101

    Network policy(下文简称为np)的本质是通过Kubernetes(下文简称k8s)的网络插件,创建一系列的网络规则,实现细粒度控制出入口流量,从而...

    nevermosby
  • 使用monocle做拟时序分析(单细胞谱系发育)

    因为是第一个课程,所以并没有提到单细胞转录组的部分新颖分析要点,比如构建细胞谱系发育,虽然我其实在课程里面也稍微提到过一点,不过怕大家印象不深刻,所以还是有必要...

    生信技能树jimmy
  • 一周播报|马斯克发帖称特斯拉破产,自己哭晕在Model 3车旁

    4月1日,愚人节当天,特斯拉首席执行官、“推特达人”埃隆•马斯克连发几条推特,戏称特斯拉已“完全破产”。

    养码场

扫码关注云+社区

领取腾讯云代金券