快速可视化高表达基因的实现

生信技能树jimmy

发布于 2020-11-25 11:04:49

1.3K1

在单细胞测序下游分析中，当重点关注哪些基因在所有细胞平均表达显著时，可选取所选取的top基因进行可视化。

0、scater包方法

在scaterR包里，有个函数plotHighestExprs，可自动根据SingleCellExperiment对象数据进行绘制。但有个缺点就是太慢了，具体操作如下

library(SingleCellExperiment)
sce <- SingleCellExperiment(
  assays = list(counts = expr)
)
class(sce)
library(scater)
sce <- logNormCounts(sce) 
p1 <- plotHighestExprs(sce, exprs_values = "logcounts",n=20)
#因为出图比较慢，故保存为本地后再观察可提高效率
ggsave("./p1.pdf", plot = p1)

下面介绍一种手动绘制箱图的方法，不仅出图快，而且箱图反映的信息也更多一些

1、获取单细胞表达矩阵

具体可分为三种情况

直接获得原始表达矩阵；

dim(expr)
expr[1:4,1:4]

expr

从Seurat对象获取；

scRNA = CreateSeuratObject(counts=expr)
class(scRNA)
#[1] "Seurat"
#attr(,"package")
#[1] "Seurat"
scRNA@assays$RNA@counts[1:4,1:4]
counts <- scRNA@assays$RNA@counts[1:4,1:4]

scRNA

注意的是Seurat是以稀疏矩阵dgCMatrix格式储存的

从SingleCellExpriment对象获取；

sce <- SingleCellExperiment(
  assays = list(counts = expr)
)
class(sce)
#[1] "SingleCellExperiment"
#attr(,"package")
#[1] "SingleCellExperiment"
assay(sce, "counts")[1:4,1:4]
counts <- assay(sce, "counts")[1:4,1:4]

sce

2、原始表达矩阵标准化

最主要的原始是，不同细胞的文库大小可能因测序过程存在差异。通过标准化，可使基因在不同细胞的表达情况具有可比性。

hist(colSums(expr))

histogram

手动对原始矩阵标准化

cell_lib <- colSums(expr) #计算细胞文库
normalized <- t(t(expr)/cell_lib)*100

上述两步也可合并成一步，但可能跳的有点多，关键理解如下示例的第三步的除法：（1）R语言的计算是向量化的；例如可以进行向量间加减乘除运算，具体规则，自己尝试下理解更深刻；（2）应用到矩阵时，可以理解为一行代表向量的一个元素。

示例

# 最后乘100是将小数更友好得表示（百分比）
normalized[1:4,1:4]

normalized[1:4,1:4]

如果是Seurat对象、SingleCellExperiment对象，则可以运用相应的函数计算，再导出标准化矩阵。

scRNA <- NormalizeData(scRNA)
scRNA@assays$RNA@data[1:4,1:4]

如下图，结果与我们上面计算的有点不同；查看NormalizeData帮助文档可知，其默认方法是计算每个细胞中基因表达量与文库的比值，然后乘一个size.factor(一般是10000)，最后进行log转换（加1，避免0以及零点几的无意义结果log1p）

基于Seurat的NormalizeData

NormalizeData

library("scater")
sce <- logNormCounts(sce) 
assay(sce, "logcounts")[1:4,1:4]

logNormCounts

如上图结果，与之前两种方法结果还不一样。scater包的logNormCounts的主要思路是首先平衡每个细胞的文库大小，再计算相对于所在细胞文库的比例。具体流程可见下图的示例。

logNormCounts

综上不同的方法，可得到不同的标准化结果。具体采用哪种方法，可自己根据情况选择。

3、选取top基因，进行箱图可视化

这里我选择了第一种，自己简单进行的标准化的结果。

#选择平均在每个细胞表达最显著的基因
most_exp <- sort(rowSums(normalized),T)[20:1] / ncol(expr)
# 可以想一下为什么设置[20:1]，与下面绘图有关

boxplot

boxplot(t(normalized[names(most_exp),]),
        cex=.01, las=1, cex.axis = .5,
        #分别设置离群点大小，轴标签的大小及方向
        xlab="% total count per cell",
        col=scales::hue_pal()(20)[20:1], #配色
        horizontal=TRUE)

boxplot

ggplot2

library(ggplot2)
library(reshape2)
newdata <- melt(t(normalized[names(most_exp),]))
head(newdata)
ggplot(data=newdata, aes(x=Var2, y=value, fill=Var2)) +
  geom_boxplot(outlier.size = 0.5) +
  coord_flip() + 
  theme(legend.position="none") +
  labs(title="High expression gene", x="", 
       y="% total count per cell")