快速可视化高表达基因的实现
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]
从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]
注意的是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]
2、原始表达矩阵标准化
最主要的原始是，不同细胞的文库大小可能因测序过程存在差异。通过标准化，可使基因在不同细胞的表达情况具有可比性。
hist(colSums(expr))
手动对原始矩阵标准化
cell_lib <- colSums(expr) #计算细胞文库
normalized <- t(t(expr)/cell_lib)*100
上述两步也可合并成一步，但可能跳的有点多，关键理解如下示例的第三步的除法：（1）R语言的计算是向量化的；例如可以进行向量间加减乘除运算，具体规则，自己尝试下理解更深刻；（2）应用到矩阵时，可以理解为一行代表向量的一个元素。
# 最后乘100是将小数更友好得表示（百分比）
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）
library("scater")
sce <- logNormCounts(sce)
assay(sce, "logcounts")[1:4,1:4]
如上图结果，与之前两种方法结果还不一样。scater包的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)
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")
最后总结下：主要根据单细胞表达矩阵，绘制箱图来可视化高表达基因。中间介绍了下如何让同一基因在不同细胞表达具有可比性的几种方法。
References
[1] https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/seurat/seurat_01_qc.html